alg111.txt

Numerical Anaysis 8th Edition Burden and Faires (Maple Source)

> 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.
> alg111 := proc() local P, Q, R, OK, A, B, ALPHA, BETA, N, FLAG, NAME, OUP, H, U1, U2, V1, V2, I, X, T, K11, K12, K21, K22, K31, K32, K41, K42, U, V, W1, Z, W2;
> printf(`This is the Linear Shooting Method.\n`);
> printf(`Input the functions P(X), Q(X) and R(X) in terms of x, separated by spaces.\n`);
> printf(`For example: -2/x 2/(x^2) sin(log(x))/(x^2)\n`);
> P := scanf(`%a`)[1];
> Q := scanf(`%a`)[1];
> R := scanf(`%a`)[1];
> P := unapply(P,x);
> Q := unapply(Q,x);
> R := unapply(R,x);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input left and right endpoints separated by blank.\n`);
> A := scanf(`%f`)[1];
> B := scanf(`%f`)[1];
> if A >= B then
> printf(`Left endpoint must be less than right endpoint.\n`);
> else
> OK := TRUE;
> fi;
> od;
> printf(`Input Y(  %.10e).\n`, A);
> ALPHA := scanf(`%a`)[1];
> printf(`Input Y(  %.10e).\n`, B);
> BETA := scanf(`%a`)[1];
> OK := FALSE;
> while OK = FALSE do
> printf(`Input a positive integer for the number of subintervals.\n`);
> N := scanf(`%d`)[1];
> if N <= 0 then
> printf(`Number must be a positive integer.\n`);
> else
> OK := TRUE;
> fi;
> od;
> if OK = TRUE then
> printf(`Choice of output method:\n`);
> printf(`1. Output to screen\n`);
> printf(`2. Output to text File\n`);
> printf(`Please enter 1 or 2.\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 2 then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example  A:\\OUTPUT.DTA\n`);
> NAME := scanf(`%s`)[1];
> OUP := fopen(NAME,WRITE,TEXT);
> else
> OUP := default;
> fi;
> fprintf(OUP, `LINEAR SHOOTING METHOD\n\n`);
> fprintf(OUP, `  I  X(I)        W(1,I)      W(2,I)\n`);
> # Step 1
> H := (B-A)/N;
> U1 := ALPHA;
> U2 := 0;
> V1 := 0;
> V2 := 1;
> # Step 2
> for I from 1 to N do
> # Step 3
> X := A+(I-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,I-1] := U1;
> U[1,I-1] := U2;
> V[0,I-1] := V1;
> V[1,I-1] := V2;
> od;
> # Step 5
> W1 := ALPHA;
> Z := (BETA-U[0,N-1])/V[0,N-1];
> X := A;
> I := 0;
> fprintf(OUP, `%3d %11.8f %11.8f %11.8f\n`, I, X, W1, Z);
> # Step 6
> for I from 1 to N do
> X := A+I*H;
> W1 := U[0,I-1]+Z*V[0,I-1];
> W2 := U[1,I-1]+Z*V[1,I-1];
> fprintf(OUP, `%3d %11.8f %11.8f %11.8f\n`, I, X, W1, W2);
> od;
> # Step 7
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg111();