alg054.txt

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

> restart;
> # ADAMS-FOURTH ORDER PREDICTOR-CORRECTOR ALGORITHM 5.4
> #
> # To approximate the solution of the initial value problem
> #        y' = f(t,y), a <= t <= b, y(a) = alpha,
> # at N+1 equally spaced points in the interval [a,b].
> #
> # INPUT:   endpoints a,b; initial condition alpha; integer N.
> #
> # OUTPUT:  approximation w to y at the (N+1) values of t.
> alg054 := proc() local F, OK, A, B, ALPHA, N, FLAG, NAME, OUP, H, T, W, I, K1, K2, K3, K4, T0, W0, J;
> printf(`This is Adams-Bashforth 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;
> printf(`Input the initial condition\n`);
> ALPHA := scanf(`%f`)[1];
> OK := FALSE;
> while OK = FALSE do
> printf(`Input an integer > 3 for the number of subintervals\n`);
> N := scanf(`%d`)[1];
> if N <= 3 then
> printf(`Number must be greater than 3\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, `ADAMS-BASHFORTH FOURTH ORDER PREDICTOR CORRECTOR METHOD\n\n`);
> fprintf(OUP, `    t           w\n`);
> # Step 1
> H := (B-A)/N;
> T[0] := A;
> W[0] := ALPHA;
> fprintf(OUP, `%5.3f %11.7f\n`, T[0], W[0]);
> # Step 2
> for I from 1 to 3 do
> # Steps 3 and 4
> # Compute starting values using Runger-Kutta 4
> T[I] := T[I-1]+H;
> K1 := H*F(T[I-1], W[I-1]);
> K2 := H*F(T[I-1]+0.5*H, W[I-1]+0.5*K1);
> K3 := H*F(T[I-1]+0.5*H, W[I-1]+0.5*K2);
> K4 := H*F(T[I], W[I-1]+K3);
> W[I] := W[I-1]+(K1+2.0*(K2+K3)+K4)/6.0;
> # Step 5
> fprintf(OUP, `%5.3f %11.7f\n`, T[I], W[I]);
> od;
> # Step 6
> for I from 4 to N do
> # Step 7
> # T0 and W0 are used in place of t, w resp.
> T0 := A+I*H;
> # Predict W(I)
> Part1 := 55.0*F(T[3],W[3])-59.0*F(T[2],W[2])+37.0*F(T[1],W[1]);
> Part2 :=-9.0*F(T[0],W[0]);
> W0 := W[3]+H*(Part1+Part2)/24.0;
> # Correct W(I)
> Part1 := 9.0*F(T0,W0)+19.0*F(T[3],W[3])-5.0*F(T[2],W[2]);
> Part2 := F(T[1],W[1]);
> W0 := W[3]+H*(Part1+Part2)/24.0;
> # Step 8
> fprintf(OUP, `%5.3f %11.7f\n`, T0, W0);
> # Step 9
> # Prepare for next iteration
> for J from 1 to 3 do
> T[J-1] := T[J];
> W[J-1] := W[J];
> od;
> # Step 10
> T[3] := T0;
> W[3] := W0;
> od;
> fi;
> # Step 11
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> RETURN(0);
> end;
> alg054();