alg052.txt

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

> restart;
> # RUNGE-KUTTA (ORDER 4) ALGORITHM 5.2
> #
> # TO APPROXIMATE THE SOLUTION TO THE INITIAL VALUE PROBLEM:
> #            Y' = F(T,Y), A<=T<=B, Y(A) = ALPHA,
> # AT (N+1) EQUALLY SPACED NUMBERS 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.
> alg052 := proc() local F, OK, A, B, ALPHA, N, FLAG, NAME, OUP, H, T, W, I, K1, K2, K3, K4;
> printf(`This is the Runge-Kutta Order Four 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 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, `RUNGE-KUTTA FOURTH ORDER METHOD\n\n`);
> fprintf(OUP, `    t           w\n\n`);
> # Step 1
> H := (B-A)/N;
> T := A;
> W := ALPHA;
> fprintf(OUP, `%5.3f %11.7f\n`, T, W);
> # Step 2
> for I from 1 to N do
> # Step 3
> # Use K1, K2, K3, K4 for K(1), K(2), K(3), K(4) resp.
> K1 := H*F(T,W);
> K2 := H*F(T+H/2.0, W+K1/2.0);
> K3 := H*F(T+H/2.0, W+K2/2.0);
> K4 := H*F(T+H,W+K3);
> # Step 4
> # Compute W(I)
> W := W+(K1+2.0*(K2+K3)+K4)/6.0;
> # Compute T(I)
> T := A+I*H;
> # Step 5
> fprintf(OUP, `%5.3f %11.7f\n`, T, W);
> od;
> # Step 6
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg052();