restart:

# EXTRAPOLATION ALGORITHM 5.6

# 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 minimum stepsize was exceeded.

print(`This is Gragg Extrapolation`):

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(`File name is `):print(NAME):

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

else

OUP := default:

fi:

fprintf(OUP, `GRAGG EXTRAPOLATION\134n\134n`):

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

# STEP 1

NK[0] := 2:

NK[1] := 4:

for J from 1 to 3 do

I1 := 2*J:

NK[I1] := 3*NK[I1-1]/2:

NK[I1+1] := 2*NK[I1-1]:

od:

# STEP 2

T0 := A:

W0 := ALPHA:

H := HMAX:

# DONE is used in place of FLAG to exit the loop in Step 4

DONE := FALSE:

# STEP 3

for I1 from 1 to 7 do

for J from 1 to I1 do

Q[I1-1,J-1] := (NK[I1]*1/NK[J-1])*(NK[I1]*1/NK[J-1]):

od:

# STEP 4

while DONE = FALSE do

# STEP 5

K := 1:

# when desired accuracy achieved, NFLAG is set to 1

NFLAG := 0:

# STEP 6

while K <= 8 and NFLAG = 0 do

# STEP 7

HK := H/NK[K-1]:

T := T0:

W2 := W0:

# Euler first step

W3 := W2+HK*F(T, W2):

T := T0+HK:

# STEP 8

M := NK[K-1]-1:

for J from 1 to M do

W1 := W2:

W2 := W3:

# midpoint method

W3 := W1+2*HK*F(T, W2):

T := T0+(J+1)*HK:

od:

# STEP 9

# endpoint correction to compute Y(K,1)

Y[K] := (W3+W2+HK*F(T, W3))/2:

# STEP 10

# NOTE: Y(K-1)=Y(K-1,1),Y(K-2)=Y(K-2,2),...,

# Y(1)=Y(K-1,K-1) since only previous row of table

# is saved

if K >= 2 then

# STEP 11

J := K:

# save Y(K-1,K-1)

V := Y[1]:

# STEP 12

while J >= 2 do

# extrapolation to compute

# Y(J-1) := Y(K,K-J+2)

Y[J-1] := Y[J]+(Y[J]-Y[J-1])/(Q[K-2,J-2]-1):

J := J-1:

od:

# STEP 13

if abs(Y[1] - V) <= TOL then

NFLAG := 1:

fi:

# Y(1) accepted as new w

fi:

# STEP 14

K := K+1:

od:

# STEP 15

K := K-1:

# STEP 16

if NFLAG = 0 then

# STEP 17

# new value for w rejected, decrease H

H := H/2:

# STEP 18

if H < HMIN then

fprintf(OUP, `HMIN exceeded\134n`):

DONE := TRUE:

fi:

else

# STEP 19

# new value for w accepted

W0 := Y[1]:

T0 := T0 + H:

fprintf(OUP, `%12.8f %11.8f %11.8f %6d\134n`, T0, W0, H, K):

# STEP 20

# increase H if possible

if T0 >= B then

DONE := TRUE:

# Procedure completed successfully.

else

if T0 + H > B then

H := B - T0:

# Terminate at t = B.

else

if K <= 3 then

H := 2*H:

fi:

if H > HMAX then

H := H/2:

fi:

od:

# STEP 21

if OUP <> default then

fclose(OUP):

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

fi:

STxUaGlzfmlzfkdyYWdnfkV4dHJhcG9sYXRpb25HNiI=

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=

JCIiIiEiKg==

SUlJbnB1dH5taW5pbXVtfmFuZH5tYXhpbXVtfm1lc2h+c3BhY2luZ35+RzYi

STZzZXBhcmF0ZWR+Ynl+YX5ibGFuay5HNiI=

STNNaW5pbXVtfnNwYWNpbmd+PX5HNiI=

JCIiJiEiIw==

STNNYXhpbXVtfnNwYWNpbmd+PX5HNiI=

JCIjRCEiIw==

STlDaG9pY2V+b2Z+b3V0cHV0fm1ldGhvZDpHNiI=

STQxLn5PdXRwdXR+dG9+c2NyZWVuRzYi

STcyLn5PdXRwdXR+dG9+dGV4dH5maWxlRzYi

STRQbGVhc2V+ZW50ZXJ+MX5vcn4yRzYi

STBZb3VyfmNob2ljZX5pc35HNiI=

IiIi

GRAGG EXTRAPOLATION T W H K 0.25000000 0.92048729 0.25000000 5 0.50000000 1.42563937 0.25000000 5 0.75000000 2.00399999 0.25000000 6 0.87500000 2.31618735 0.12500000 8 0.93750000 2.47711152 0.06250000 3 1.06250000 2.80710828 0.12500000 4 1.18750000 3.14571936 0.12500000 8 1.31250000 3.48993088 0.12500000 3 1.56250000 4.18103966 0.25000000 6 1.81250000 4.84728492 0.25000000 7 2.00000000 5.30547195 0.18750000 4

JSFH