alg112.txt

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

> restart;
> # NONLINEAR SHOOTING ALGORITHM 11.2
> #
> # To approximate the solution of the nonlinear boundary-value problem
> #
> #          Y'' = F(X,Y,Y'), A<=X<=B, Y(A) = ALPHA, Y(B) = BETA:
> #
> #
> # INPUT:   Endpoints A,B; boundary conditions ALPHA, BETA; number of
> #          subintervals N; tolerance TOL; maximum number of iterations M.
> #
> # OUTPUT:  Approximations W(1,I) TO Y(X(I)); W(2,I) TO Y'(X(I))
> #          for each I=0,1,...,N or a message that the maximum
> #          number of iterations was exceeded.
> alg112 := proc() local F, FY, FYP, OK, A, B, ALPHA, BETA, TK, AA, N, TOL, NN, FLAG, NAME, OUP, H, K, W1, W2, U1, U2, I, X, T, K11, K12, K21, K22, K31, K32, K41, K42, J;
> printf(`This is the Nonlinear Shooting Method.\n`);
> printf(`Input the function F(X,Y,Z) in terms of x, y, z.\n`);
> printf(`For example:   (32+2*x^3-y*z)/8\n`);
> F := scanf(`%a`)[1];
> FY := diff(F,y);
> FYP := diff(F,z);
> F := unapply(F,x,y,z);
> FY := unapply(FY,x,y,z);
> FYP := unapply(FYP,x,y,z);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input left and right endpoints separated by blank.\n`);
> A := scanf(`%e`)[1];
> B := scanf(`%e`)[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(`%e`)[1];
> printf(`Input Y(%.10e).\n`, B);
> BETA := scanf(`%e`)[1];
> TK := (BETA-ALPHA)/(B-A);
> printf(`TK := %.8e\n`, TK);
> printf(`Input new TK? Enter Y or N.\n`);
> AA := scanf(`\n%c`)[1];
> if AA = "Y" or AA = "y" then
> printf(`input new TK\n`);
> TK := scanf(`%f`)[1];
> fi;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input an integer > 1 for the number of subintervals.\n`);
> N := scanf(`%d`)[1];
> if N <= 1 then
> printf(`Number must exceed 1.\n`);
> else
> OK := TRUE;
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input Tolerance.\n`);
> TOL := scanf(`%f`)[1];
> if TOL <= 0 then
> printf(`Tolerance must be positive.\n`);
> else
> OK := TRUE;
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input maximum number of iterations.\n`);
> NN := scanf(`%d`)[1];
> if NN <= 0 then
> printf(`Must be 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, `NONLINEAR SHOOTING METHOD\n\n`);
> fprintf(OUP, `  I    X(I)         W1(I)        W2(I)\n`);
> # Step 1
> H := (B-A)/N;
> K := 1;
> # TK is already computed.
> OK := FALSE;
> # Step 2
> while K <= NN and OK = FALSE do
> # Step 3
> W1[0] := ALPHA;
> W2[0] := TK;
> U1 := 0 ;
> U2 := 1;
> # Step 4
> # Runge-Kutta method for systems is used in Steps 5 and 6
> for I from 1 to N do
> # Step 5
> X := A+(I-1)*H;
> T := X+0.5*H;
> # Step 6
> K11 := H*W2[I-1];
> K12 := H*F(X,W1[I-1],W2[I-1]);
> K21 := H*(W2[I-1]+0.5*K12);
> K22 := H*F(T,W1[I-1]+0.5*K11,W2[I-1]+0.5*K12);
> K31 := H*(W2[I-1]+0.5*K22);
> K32 := H*F(T,W1[I-1]+0.5*K21,W2[I-1]+0.5*K22);
> K41 := H*(W2[I-1]+K32);
> K42 := H*F(X+H,W1[I-1]+K31,W2[I-1]+K32);
> W1[I] := W1[I-1]+(K11+2*(K21+K31)+K41)/6;
> W2[I] := W2[I-1]+(K12+2*(K22+K32)+K42)/6;
> K11 := H*U2;
> K12 := H*(FY(X,W1[I-1],W2[I-1])*U1+FYP(X,W1[I-1],W2[I-1])*U2);
> K21 := H*(U2+0.5*K12);
> K22 := H*(FY(T,W1[I-1],W2[I-1])*(U1+0.5*K11)+FYP(T,W1[I-1],W2[I-1])*(U2+0.5*K21));
> K31 := H*(U2+0.5*K22);
> K32 := H*(FY(T,W1[I-1],W2[I-1])*(U1+0.5*K21)+FYP(T,W1[I-1],W2[I-1])*(U2+0.5*K22));
> K41 := H*(U2+K32);
> K42 := H*(FY(X+H,W1[I-1],W2[I-1])*(U1+K31)+FYP(X+H,W1[I-1],W2[I-1])*(U2+K32));
> U1 := U1+(K11+2*(K21+K31)+K41)/6;
> U2 := U2+(K12+2*(K22+K32)+K42)/6;
> od;
> # Step 7
> # Test for accuracy
> if abs(W1[N]-BETA) < TOL then
> # Step 8
> I := 0;
> fprintf(OUP, `%3d %13.8f %13.8f %13.8f\n`, I, A, ALPHA, TK);
> for I from 1 to N do
> J := I+1;
> X := A+I*H;
> fprintf(OUP, `%3d %13.8f %13.8f %13.8f\n`, I, X, W1[J-1], W2[J-1]);
> od;
> fprintf(OUP, `Convergence in %d iterations\n`, K);
> fprintf(OUP, ` t = %14.7e\n`, TK);
> # Step 9
> OK := TRUE;
> else
> # Step 10
> # Newton's method applied to improve TK
> TK := TK-(W1[N]-BETA)/U1;
> K := K+1;
> fi;
> od;
> # Step 11
> # Procedure completed unsuccessfully.
> if OK = FALSE then
> fprintf(OUP, `Method failed after %d iterations\n`, NN);
> fi;
> fi;
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> RETURN(0);
> end;
> alg112();