restart:
# LINEAR SHOOTING ALGORITHM 11.1
#
# To approximate the solution of the boundary-value problem
#
# -Y'' + P(X)Y' + Q(X)Y + R(X) = 0, A<=X<=B, Y(A)=ALPHA, Y(B)=BETA:
#
#
# INPUT: Endpoints A,B: boundary conditions ALPHA, BETA: number of
# subintervals N.
#
# OUTPUT: Approximations W(1,I) to Y(X(I)): W(2,I) to Y'(X(I))
# for each I=0,1,...,N.
print(`This is the Linear Shooting Method.`):
print(`Input the functions P(X), Q(X) and R(X) in terms of x, separated by spaces.`):
print(`For example: -2/x 2/(x^2) sin(log(x))/(x^2)`):
P := scanf(`%a`)[1]: print(`P(x) = `):print(P):
Q := scanf(`%a`)[1]: print(`Q(x) = `):print(P):
R := scanf(`%a`)[1]: print(`R(x) = `):print(R):
P := unapply(P,x):
Q := unapply(Q,x):
R := unapply(R,x):
OK := FALSE:
while OK = FALSE do
print(`Input left and right endpoints separated by blank.`):
A := scanf(`%f`)[1]: print(`Left endpoint = `):print(A):
B := scanf(`%f`)[1]: print(`Right endpoint = `):print(B):
if A >= B then
print(`Left endpoint must be less than right endpoint.`):
else
OK := TRUE:
fi:
od:
print(`Input Y(A).`):
ALPHA := scanf(`%a`)[1]: print(`alpha = `):print(ALPHA):
print(`Input Y(B).`):
BETA := scanf(`%a`)[1]: print(`beta = `):print(BETA):
OK := FALSE:
while OK = FALSE do
print(`Input a positive integer for the number of subintervals.`):
N := scanf(`%d`)[1]: print(`Input = `): print(N):
if N <= 0 then
print(`Number must be a positive integer.`):
else
OK := TRUE:
fi:
od:
if OK = TRUE then
print(`Choice of output method:`):
print(`1. Output to screen`):
print(`2. Output to text File`):
print(`Please enter 1 or 2.`):
FLAG := scanf(`%d`)[1]: print(`Input = `): print(FLAG):
if FLAG = 2 then
print(`Input the file name in the form - drive:\134\134name.ext`):
print(`for example A:\134\134OUTPUT.DTA`):
NAME := scanf(`%s`)[1]: print(`Output file is `): print(NAME):
OUP := fopen(NAME,WRITE,TEXT):
else
OUP := default:
fi:
fprintf(OUP, `LINEAR SHOOTING METHOD\134n\134n`):
fprintf(OUP, ` I X(I) W(1,I) W(2,I)\134n`):
# Step 1
H := (B-A)/N:
U1 := ALPHA:
U2 := 0:
V1 := 0:
V2 := 1:
# Step 2
for I2 from 1 to N do
# Step 3
X := A+(I2-1)*H:
T := X+0.5*H:
# Step 4
K11 := H*U2:
K12 := H*(P(X)*U2+Q(X)*U1+R(X)):
K21 := H*(U2+0.5*K12):
K22 := H*(P(T)*(U2+0.5*K12)+Q(T)*(U1+0.5*K11)+R(T)):
K31 := H*(U2+0.5*K22):
K32 := H*(P(T)*(U2+0.5*K22)+Q(T)*(U1+0.5*K21)+R(T)):
T := X+H:
K41 := H*(U2+K32):
K42 := H*(P(T)*(U2+K32)+Q(T)*(U1+K31)+R(T)):
U1 := U1+(K11+2*(K21+K31)+K41)/6:
U2 := U2+(K12+2*(K22+K32)+K42)/6:
K11 := H*V2:
K12 := H*(P(X)*V2+Q(X)*V1):
T := X+0.5*H:
K21 := H*(V2+0.5*K12):
K22 := H*(P(T)*(V2+0.5*K12)+Q(T)*(V1+0.5*K11)):
K31 := H*(V2+0.5*K22):
K32 := H*(P(T)*(V2+0.5*K22)+Q(T)*(V1+0.5*K21)):
T := X+H:
K41 := H*(V2+K32):
K42 := H*(P(T)*(V2+K32)+Q(T)*(V1+K31)):
V1 := V1+(K11+2*(K21+K31)+K41)/6:
V2 := V2+(K12+2*(K22+K32)+K42)/6:
U[0,I2-1] := U1:
U[1,I2-1] := U2:
V[0,I2-1] := V1:
V[1,I2-1] := V2:
od:
# Step 5
W1 := ALPHA:
Z := (BETA-U[0,N-1])/V[0,N-1]:
X := A:
I2 := 0:
fprintf(OUP, `%3d %11.8f %11.8f %11.8f\134n`, I2, X, W1, Z):
# Step 6
for I2 from 1 to N do
X := A+I2*H:
W1 := U[0,I2-1]+Z*V[0,I2-1]:
W2 := U[1,I2-1]+Z*V[1,I2-1]:
fprintf(OUP, `%3d %11.8f %11.8f %11.8f\134n`, I2, X, W1, W2):
od:
# Step 7
if OUP <> default then
fclose(OUP):
print(`Output file `,NAME,` created successfully`):
fi:
fi:
JSFH