restart:

# ROMBERG ALGORITHM 4.2

# To approximate I = integral ( f(x) dx ) from a to b:

# INPUT: endpoints a, b: integer n.

# OUTPUT: an array R. ( R(2,n) is the approximation to I. )

# R is computed by rows: only 2 rows saved in storage

print(`This is Romberg Integration. \134n`):

print(`Input the function F(x) in terms of x `):

print(`For example: cos(x) `):

F := scanf(`%a`)[1]:print(`F(x) = `):print(F):

F := unapply(F,x):

OK := FALSE:

while OK = FALSE do

print(`Input lower limit of integration and `):

print(`upper limit of integration `):

print(`separated by a blank `):

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

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

if A > B then

print(`Lower limit must be less than upper limit `):

else

OK := TRUE:

fi:

od:

OK := FALSE:

while OK = FALSE do

print(`Input number of rows - no decimal point. `):

N := scanf(`%d`)[1]:print(`N = `):print(N):

if N > 0 then

OK := TRUE:

else

print(`Number must be a positive integer `):

fi:

od:

if OK = TRUE then

# STEP 1

H := B-A:

R[0,0] := (F(A)+F(B))/2*H:

# STEP 2

OUP := default:

fprintf(OUP,`Romberg Integration \134n`):

fprintf(OUP,`Initial Data:\134n`);

fprintf(OUP,`Limits of integration = [%12.8f, %12.8f]\134n`, A, B);

fprintf(OUP,`Number of rows = %3d\134n`, N);

fprintf(OUP,`\134nRomberg Integration Table:\134n`);

fprintf(OUP,`\134n%12.8f\134n\134n`, R[0,0]);

# STEP 3

for I1 from 2 to N do

# STEP 4

# approximation from Trapezoidal method

SUM := 0:

M := 2^(I1-2):

for K from 1 to M do

SUM := SUM+F(A+(K-0.5)*H):

od:

R[1,0] := (R[0,0]+H*SUM)/2:

# STEP 5

# extrapolation

for J from 2 to I1 do

L := 2^(2*(J-1)):

R[1,J-1] := R[1,J-2]+(R[1,J-2]-R[0,J-2])/(L-1):

od:

# STEP 6

for K from 1 to I1 do

fprintf(OUP,` %11.8f`,R[1,K-1]):

od:

fprintf(OUP,`\134n\134n`);

# STEP 7

H := H/2:

# since only two rows are kept in storage, this step is

# to prepare for the next row

# update row 1 of R

# STEP 8

for J from 1 to I1 do

R[0,J-1] := R[1,J-1]:

od:

fi:

ST9UaGlzfmlzflJvbWJlcmd+SW50ZWdyYXRpb24ufnwrRzYi

SUdJbnB1dH50aGV+ZnVuY3Rpb25+Rih4KX5pbn50ZXJtc35vZn54fkc2Ig==

STVGb3J+ZXhhbXBsZTp+Y29zKHgpfkc2Ig==

SShGKHgpfj1+RzYi

LUkkc2luRzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliRzYiNiNJInhHRic=

SUZJbnB1dH5sb3dlcn5saW1pdH5vZn5pbnRlZ3JhdGlvbn5hbmR+RzYi

STx1cHBlcn5saW1pdH5vZn5pbnRlZ3JhdGlvbn5HNiI=

STZzZXBhcmF0ZWR+Ynl+YX5ibGFua35HNiI=

SSVhfj1+RzYi

JCIiIUYj

SSVifj1+RzYi

JCIrYUVmVEohIio=

SUpJbnB1dH5udW1iZXJ+b2Z+cm93c34tfm5vfmRlY2ltYWx+cG9pbnQufkc2Ig==

SSVOfj1+RzYi

IiIn

Romberg Integration Initial Data: Limits of integration = [ 0.00000000, 3.14159265] Number of rows = 6 Romberg Integration Table:

-0.00000000 1.57079633 2.09439510 1.89611890 2.00455976 1.99857073 1.97423160 2.00026917 1.99998313 2.00000555 1.99357034 2.00001659 1.99999975 2.00000002 1.99999999 1.99839336 2.00000104 2.00000000 2.00000000 2.00000000 2.00000000

JSFH