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.

print(`This is the Runge-Kutta-Fehlberg Method.`):

print(`Input the function f(t,y) in terms of t and y`):

print(`For example: y-t^2+1`):

F := scanf(`%a`)[1]:print(`f(t,y) = `):print(F);

F := unapply(F,t,y):

OK := FALSE:

while OK = FALSE do

print(`Input left and right endpoints separated by blank`):

A := scanf(`%f`)[1]:print(`a = `):print(A):

B := scanf(`%f`)[1]:print(`b = `):print(B):

if A >= B then

print(`Left endpoint must be less than right endpoint`):

else

OK := TRUE:

fi:

od:

print(`Input the initial condition`):

ALPHA := scanf(`%f`)[1]:print(`y(a)= `):print(ALPHA):

OK := FALSE:

while OK = FALSE do

print(`Input tolerance`):

TOL := scanf(`%f`)[1]:print(`Tolerance= `):print(TOL):

if TOL <= 0 then

print(`Tolerance must be a positive.`):

else

OK := TRUE:

fi:

od:

OK := FALSE:

while OK = FALSE do

print(`Input minimum and maximum mesh spacing separated by a space`):

HMIN := scanf(`%f`)[1]:print(`Minimum spacing = `):print(HMIN):

HMAX := scanf(`%f`)[1]:print(`Maximum spacing = `):print(HMAX):

if HMIN < HMAX and HMIN > 0 then

OK := TRUE:

else

print(`Minimum mesh spacing must be a positive real number`):

print(`and less than the maximum mesh spacing`):

fi:

od:

if OK = TRUE then

print(`Choice of output method:`):

print(`1. Output to screen`):

print(`2. Output to text file`):

print(`Please enter 1 or 2`):

FLAG := scanf(`%d`)[1]:print(`Your choice= `):print(FLAG):

if FLAG = 2 then

print(`Input the file name in the form - drive:\134\134name.ext`):

print(`For example A:\134\134OUTPUT.DTA`):

NAME := scanf(`%s`)[1]:print(`Output file is `):print(NAME):

OUP := fopen(NAME,WRITE,TEXT):

else

OUP := default:

fi:

fprintf(OUP, `RUNGE-KUTTA-FEHLBERG METHOD\134n\134n`):

fprintf(OUP, ` T(I) W(I) H R\134n\134n`):

# Step 1

H := HMAX:

T := A:

W := ALPHA:

fprintf(OUP, `%12.7f %11.7f 0 0\134n`, 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\134n`, 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:

# 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:

od:

if OK = FALSE then

fprintf(OUP, `Minimal H exceeded\134n`):

fi:

# Step 12

# Process is complete

if OUP <> default then

fclose(OUP):

print(`Output file`,NAME,` created successfully`):

fi:

SUlUaGlzfmlzfnRoZX5SdW5nZS1LdXR0YS1GZWhsYmVyZ35NZXRob2QuRzYi

SU5JbnB1dH50aGV+ZnVuY3Rpb25+Zih0LHkpfmlufnRlcm1zfm9mfnR+YW5kfnlHNiI=

STVGb3J+ZXhhbXBsZTp+eS10XjIrMUc2Ig==

SSpmKHQseSl+PX5HNiI=

LChJInlHNiIiIiIqJEkidEdGJCIiIyEiIkYlRiU=

SVJJbnB1dH5sZWZ0fmFuZH5yaWdodH5lbmRwb2ludHN+c2VwYXJhdGVkfmJ5fmJsYW5rRzYi

SSVhfj1+RzYi

JCIiIUYj

SSVifj1+RzYi

JCIiIyIiIQ==

STxJbnB1dH50aGV+aW5pdGlhbH5jb25kaXRpb25HNiI=

SSd5KGEpPX5HNiI=

JCIiJiEiIg==

STBJbnB1dH50b2xlcmFuY2VHNiI=

SSxUb2xlcmFuY2U9fkc2Ig==

JCIiIiEiJg==

SWZuSW5wdXR+bWluaW11bX5hbmR+bWF4aW11bX5tZXNofnNwYWNpbmd+c2VwYXJhdGVkfmJ5fmF+c3BhY2VHNiI=

STNNaW5pbXVtfnNwYWNpbmd+PX5HNiI=

JCIiIiEiIw==

STNNYXhpbXVtfnNwYWNpbmd+PX5HNiI=

JCIjRCEiIw==

STlDaG9pY2V+b2Z+b3V0cHV0fm1ldGhvZDpHNiI=

STQxLn5PdXRwdXR+dG9+c2NyZWVuRzYi

STcyLn5PdXRwdXR+dG9+dGV4dH5maWxlRzYi

STRQbGVhc2V+ZW50ZXJ+MX5vcn4yRzYi

SS5Zb3VyfmNob2ljZT1+RzYi

IiIi

RUNGE-KUTTA-FEHLBERG METHOD

T(I) W(I) H R 0.0000000 0.5000000 0 0 0.2500000 0.9204886 0.2500000 0.0000062 0.4865525 1.3964916 0.2365525 0.0000045 0.7293336 1.9537496 0.2427811 0.0000043 0.9793336 2.5864269 0.2500000 0.0000038 1.2293336 3.2604615 0.2500000 0.0000024 1.4793336 3.9520965 0.2500000 0.0000007 1.7293336 4.6308277 0.2500000 0.0000015 1.9793336 5.2574869 0.2500000 0.0000043 2.0000000 5.3054896 0.0206664 0.0000000

JSFH