alg053.txt

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

> restart;
> # RUNGE-KUTTA-FEHLBERG ALGORITHM 5.3
> #
> # 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 STEPSIZE HMAX; MINIMUM STEPSIZE HMIN.
> #
> # OUTPUT:  T, W, H WHERE W APPROXIMATES Y(T) AND STEPSIZE H WAS 
> # USED OR A MESSAGE THAT THE MINIMUM STEPSIZE WAS EXCEEDED. 
> alg053 := proc() local F, OK, A, B, ALPHA, TOL, HMIN, HMAX, FLAG, NAME, OUP, H, T, W, K1, K2, K3, K4, K5, K6, R, DELTA;
> printf(`This is the Runge-Kutta-Fehlberg 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 tolerance\n`);
> TOL := scanf(`%f`)[1];
> if TOL <= 0 then
> printf(`Tolerance must be a positive.\n`);
> else
> OK := TRUE;
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input minimum and maximum mesh spacing separated by a space\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 positive real number\n`); 
> printf(`and less than 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, `RUNGE-KUTTA-FEHLBERG METHOD\n\n`);
> fprintf(OUP, `    T(I)           W(I)          H           R\n\n`);
> # Step 1
> H := HMAX;
> T := A;
> W := ALPHA;
> fprintf(OUP, `%12.7f %11.7f    0           0\n`, T, W);
> OK := TRUE;
> # Step 2
> while T < B and OK = TRUE do
> # Step 3
> K1 := H*F(T,W);
> K2 := H*F(T+H/4,W+K1/4);K3 := H*F(T+3*H/8,W+(3*K1+9*K2)/32);
> K4 := H*F(T+12*H/13,W+(1932*K1-7200*K2+7296*K3)/2197);
> K5 := H*F(T+H,W+439*K1/216-8*K2+3680*K3/513-845*K4/4104);
> K6 := H*F(T+H/2,W-8*K1/27+2*K2-3544*K3/2565+1859*K4/4104-11*K5/40);
> # Step 4
> R := abs(K1/360-128*K3/4275-2197*K4/75240.0+K5/50+2*K6/55)/H;
> # Step 5
> if R <= TOL then
> # Step 6
> # Approximation accepted
> T := T+H;
> W := W+25*K1/216+1408*K3/2565+2197*K4/4104-K5/5;
> # Step 7
> fprintf(OUP, `%12.7f %11.7f %11.7f %11.7f\n`, T, W, H, R);
> fi;
> # Step 8
> # To avoid underflow
> if R > 1.0E-20 then
> DELTA := 0.84 * exp(0.25 * log(TOL / R));
> else
> DELTA := 10.0;
> fi;
> # Step 9
> # Calculate new H
> if DELTA <= 0.1 then
> H := 0.1*H;
> else
> if DELTA >= 4 then
> H := 4.0 * H;
> else
> H := DELTA * H;
> fi;
> fi;
> # Step 10
> if H > HMAX then
> H := HMAX;
> fi;
> # Step 11
> if H < HMIN then
> OK := FALSE;
> else
> if T+H > B then
> if abs(B-T) < TOL then
> T := B;
> else
> H := B-T;
> fi;
> fi;
> fi;
> od;
> if OK = FALSE then
> fprintf(OUP, `Minimal H exceeded\n`);
> fi;
> # Step 12
> # Process is complete
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg053();