alg055.txt

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

> restart;
> # ADAMS VARIABLE STEP-SIZE PREDICTOR-CORRECTOR ALGORITHM 5.5
> #
> # To approximate the solution of the initial value problem
> #        y' = f( t, y ), a <= t <= b, y(a) = ALPHA,
> #
> # with local truncation error within a given tolerance:
> #
> # INPUT:   endpoints a, b; initial condition ALPHA; tolerance TOL;
> #          maximum step size HMAX; minimum step size HMIN.
> #
> # OUTPUT:  I, T(I), W(I), H where at the Ith step W(I) approximates
> #          y(T(I)) and step size H was used or a message that the
> #          minimum step size was exceeded.
> alg055 := proc() local RK4, OK, A, B, ALPHA, TOL, HMIN, HMAX, FLAG, NAME, OUP, H, DONE, KK, NFLAG, I, TT, WP, WC, SIG, K, J, Q; global W,T,F,Part1,Part2;
> # Step 1
> RK4 := proc(H,X,Y,KK) local K1, K2, K3, K4; global W,T;
> K1 := H*F(X,Y);
> K2 := H*F(X+0.5*H,Y+0.5*K1);
> K3 := H*F(X+0.5*H,Y+0.5*K2);
> K4 := H*F(X+H,Y+K3);
> W[KK] := Y+(K1+2.0*(K2+K3)+K4)/6.0;
> T[KK] := X+H;
> end;
> printf(`This is the Adams Variable Step-size Predictor-Corrector Method\n`);
> printf(`Input the function F(t,y) in terms of t and y\n`);
> printf(`For example: y-t^2+1\n`);
> F := scanf(`%a`)[1];
> F := unapply(F,t,y);
> 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;
> OK := FALSE;
> printf(`Input the initial condition.\n`);
> ALPHA := scanf(`%f`)[1];
> 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 minimum and maximum mesh spacing  `);
> printf(`separated by a blank.\n`);
> HMIN := scanf(`%f`)[1];
> HMAX := scanf(`%f`)[1];
> if HMIN < HMAX and HMIN > 0 then
> OK := TRUE;
> else
> printf(`Minimum mesh spacing must be a  `);
> printf(`positive real number and less than\n`);
> printf(`the maximum mesh spacing.\n`);
> 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, `ADAMS VARIABLE STEP-SIZE PREDICTOR CORRECTOR METHOD\n\n`);
> fprintf(OUP, `           t           w           h       sigma\n`);
> # Step 2
> T[0] := A;
> W[0] := ALPHA;
> H := HMAX;
> # OK is used in place of FLAG to exit the loop in Step 4
> OK := TRUE;
> # DONE is used in place of LAST to indicate when the last
> # value is calculated
> DONE := FALSE;
> # Step 3
> for KK from 1 to 3 do
> RK4(H,T[KK-1],W[KK-1],KK);
> od;
> # NFLAG indicates calculation from RK4
> NFLAG := 1;
> I := 5;
> # Use TT in place of t
> TT := T[3] + H;
> # Step 4
> while DONE = FALSE do
> # Step 5
> # Predict W(I)
> Part1 := 55.0*F(T[I-2],W[I-2])-59.0*F(T[I-3],W[I-3]);
> Part2 := 37.0*F(T[I-4],W[I-4])-9.0*F(T[I-5],W[I-5]);
> WP := W[I-2]+H*(Part1+Part2)/24.0;
> # Correct W(I)
> Part1 := 9.0*F(TT,WP)+19.0*F(T[I-2],W[I-2])-5.0*F(T[I-3],W[I-3]);
> Part2 := F(T[I-4],W[I-4]);
> WC := W[I-2]+H*(Part1+Part2)/24.0;
> SIG := 19.0*abs(WC-WP)/(270.0*H);
> # Step 6
> if SIG <= TOL then
> # Step 7
> # Result accepted
> W[I-1] := WC;
> T[I-1] := TT;
> # Step 8
> if NFLAG = 1 then
> K := I-3;
> KK := I-1;
> # Previous results are also accepted
> for J from K to KK do
> fprintf(OUP, `%12.8f %11.8f %11.8f %11.8f\n`, T[J-1], W[J-1], H, SIG);
> od;
> fprintf(OUP, `%12.8f %11.8f %11.8f %11.8f\n`, T[I-1], W[I-1], H, SIG);
> else
> # Previous results were already accepted
> fprintf(OUP, `%12.8f %11.8f %11.8f %11.8f\n`, T[I-1], W[I-1], H, SIG);
> fi;
> # Step 9
> if OK = FALSE then
> DONE := TRUE;
> else
> # Step 10
> I := I+1;
> NFLAG := 0;
> # Step 11
> if SIG <= 0.1*TOL or T[I-2]+H > B then
> # Increase H if more accuracy has been obtained than is required, or
> # decrease H to include B as a mesh point
> # Step 12
> # To avoid underflow
> if SIG <= 1.0e-20 then
> Q := 4.0;
> else
> Q := (0.5*TOL/SIG)^(1/4);
> fi;
> # Step 13
> if Q > 4.0 then
> H := 4.0*H;
> else
> H := Q * H;
> fi;
> # Step 14
> if H > HMAX then
> H := HMAX;
> fi;
> # Step 15
> if T[I-2]+4.0*H > B then
> H := 0.25*(B-T[I-2]);
> if H < TOL then
> DONE := TRUE;
> fi;
> OK := FALSE;
> fi;
> # Step 16
> for KK from I-1 to I+2 do
> RK4(H,T[KK-1],W[KK-1],KK);
> od;
> NFLAG := 1;
> I := I+3;
> fi;
> fi;
> else
> # FALSE branch for Step 6 - result rejected
> # Step 17
> Q := (0.5*TOL/SIG)^(1/4);
> # Step 18
> if Q < 0.1 then
> H := 0.1 * H;
> else
> H := Q * H;
> fi;
> # Step 19
> if H < HMIN then
> fprintf(OUP, `HMIN exceeded\n`);
> DONE := TRUE;
> else
> if T[I-2]+4.0*H > B then
> H := 0.25*(B-T[I-2]);
> fi;
> if NFLAG = 1 then
> I := I-3;
> fi;
> for KK from I-1 to I+2 do
> RK4(H,T[KK-1],W[KK-1],KK);
> od;
> I := I+3;
> NFLAG := 1;
> fi;
> fi;
> # Step 20
> TT := T[I-2] + H;
> od;
> # Step 21
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg055();