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.

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

print(`This is the Adams Variable Step-size Predictor-Corrector 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:

OK := FALSE:

print(`Input the initial condition.`):

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

while OK = FALSE do

print(`Input tolerance.`):

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

if TOL <= 0 then

print(`Tolerance must be positive.`):

else

OK := TRUE:

fi:

od:

OK := FALSE:

while OK = FALSE do

print(`Input minimum and maximum mesh spacing `):

print(`separated by a blank.`):

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 `):

print(`positive real number and less than`):

print(`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 is `):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, `ADAMS VARIABLE STEP-SIZE PREDICTOR CORRECTOR METHOD\134n\134n`):

fprintf(OUP, ` t w h sigma\134n`):

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

I1 := 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[I1-2],W[I1-2])-59.0*F(T[I1-3],W[I1-3]):

Part2 := 37.0*F(T[I1-4],W[I1-4])-9.0*F(T[I1-5],W[I1-5]):

WP := W[I1-2]+H*(Part1+Part2)/24.0:

# Correct W(I)

Part1 := 9.0*F(TT,WP)+19.0*F(T[I1-2],W[I1-2])-5.0*F(T[I1-3],W[I1-3]):

Part2 := F(T[I1-4],W[I1-4]):

WC := W[I1-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[I1-1] := WC:

T[I1-1] := TT:

# Step 8

if NFLAG = 1 then

K := I1-3:

KK := I1-1:

# Previous results are also accepted

for J from K to KK do

fprintf(OUP, `%12.8f %11.8f %11.8f %11.8f\134n`, T[J-1], W[J-1], H, SIG):

od:

fprintf(OUP, `%12.8f %11.8f %11.8f %11.8f\134n`, T[I1-1], W[I1-1], H, SIG):

else

# Previous results were already accepted

fprintf(OUP, `%12.8f %11.8f %11.8f %11.8f\134n`, T[I1-1], W[I1-1], H, SIG):

fi:

# Step 9

if OK = FALSE then

DONE := TRUE:

else

# Step 10

I1 := I1+1:

NFLAG := 0:

# Step 11

if SIG <= 0.1*TOL or T[I1-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[I1-2]+4.0*H > B then

H := 0.25*(B-T[I1-2]):

if H < TOL then

DONE := TRUE:

fi:

OK := FALSE:

fi:

# Step 16

for KK from I1-1 to I1+2 do

RK4(H,T[KK-1],W[KK-1],KK):

od:

NFLAG := 1:

I1 := I1+3:

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\134n`):

DONE := TRUE:

else

if T[I1-2]+4.0*H > B then

H := 0.25*(B-T[I1-2]):

fi:

if NFLAG = 1 then

I1 := I1-3:

fi:

for KK from I1-1 to I1+2 do

RK4(H,T[KK-1],W[KK-1],KK):

od:

I1 := I1+3:

NFLAG := 1:

fi:

# Step 20

TT := T[I1-2] + H:

od:

# Step 21

if OUP <> default then

fclose(OUP):

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

fi:

SWpuVGhpc35pc350aGV+QWRhbXN+VmFyaWFibGV+U3RlcC1zaXplflByZWRpY3Rvci1Db3JyZWN0b3J+TWV0aG9kRzYi

SU5JbnB1dH50aGV+ZnVuY3Rpb25+Zih0LHkpfmlufnRlcm1zfm9mfnR+YW5kfnlHNiI=

STVGb3J+ZXhhbXBsZTp+eS10XjIrMUc2Ig==

SSpmKHQseSl+PX5HNiI=

LChJInlHNiIiIiIqJEkidEdGJCIiIyEiIkYlRiU=

SVNJbnB1dH5sZWZ0fmFuZH5yaWdodH5lbmRwb2ludHN+c2VwYXJhdGVkfmJ5fmJsYW5rLkc2Ig==

SSVhfj1+RzYi

JCIiIUYj

SSVifj1+RzYi

JCIiIyIiIQ==

ST1JbnB1dH50aGV+aW5pdGlhbH5jb25kaXRpb24uRzYi

SSh5KGEpfj1+RzYi

JCIiJiEiIg==

STFJbnB1dH50b2xlcmFuY2UuRzYi

SS1Ub2xlcmFuY2V+PX5HNiI=

JCIiIiEiJg==

SUlJbnB1dH5taW5pbXVtfmFuZH5tYXhpbXVtfm1lc2h+c3BhY2luZ35+RzYi

STZzZXBhcmF0ZWR+Ynl+YX5ibGFuay5HNiI=

STNNaW5pbXVtfnNwYWNpbmd+PX5HNiI=

JCIiIiEiIw==

STNNYXhpbXVtfnNwYWNpbmd+PX5HNiI=

JCIjRCEiIw==

STlDaG9pY2V+b2Z+b3V0cHV0fm1ldGhvZDpHNiI=

STQxLn5PdXRwdXR+dG9+c2NyZWVuRzYi

STcyLn5PdXRwdXR+dG9+dGV4dH5maWxlRzYi

STRQbGVhc2V+ZW50ZXJ+MX5vcn4yRzYi

STBZb3VyfmNob2ljZX5pc35HNiI=

IiIi

ADAMS VARIABLE STEP-SIZE PREDICTOR CORRECTOR METHOD t w h sigma

0.12570162 0.70023174 0.12570162 0.00000405 0.25140325 0.92309475 0.12570162 0.00000405 0.37710487 1.16738749 0.12570162 0.00000405 0.50280649 1.43174769 0.12570162 0.00000405 0.62850811 1.71463013 0.12570162 0.00000461 0.75420974 2.01428289 0.12570162 0.00000521 0.87991136 2.32871936 0.12570162 0.00000591 1.00561298 2.65568688 0.12570162 0.00000671 1.13131460 2.99263101 0.12570162 0.00000760 1.25701623 3.33665521 0.12570162 0.00000862 1.38271785 3.68447492 0.12570162 0.00000978 1.48572739 3.96974089 0.10300954 0.00000703 1.58873692 4.25276738 0.10300954 0.00000703 1.69174646 4.53100881 0.10300954 0.00000703 1.79475600 4.80164289 0.10300954 0.00000703 1.89776554 5.06154198 0.10300954 0.00000776 1.92332415 5.12397135 0.02555862 0.00000004 1.94888277 5.18547178 0.02555862 0.00000004 1.97444138 5.24598540 0.02555862 0.00000004 2.00000000 5.30545284 0.02555862 0.00000004

JSFH