restart:

# ADAPTIVE QUADRATURE ALGORITM 4.3 #

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

# a given tolerance TOL:

# INPUT: endpoints a, b: tolerance TOL: limit N to number of levels

# OUTPUT: approximation APP or message that N is exceeded.

print(`This is Adaptive Quadrature with Simpson's Method.`):

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

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

BB := scanf(`%f`)[1]:print(`a = `):print(AA):print(`b = `):print(BB):

if AA > BB then

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

else

OK := TRUE:

fi:

od:

OK := FALSE:

while OK = FALSE do

print(`Input tolerance.`):

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

if EPS > 0 then

OK := TRUE:

else

print(`Tolerance must be positive.`):

fi:

od:

OK := FALSE:

while OK = FALSE do

print(`Input the maximum number of levels.`):

N := scanf(`%d`)[1]:print(`Maximum number of levels = `):print(N):

if N > 0 then

OK := TRUE:

else

print(`Number must be positive`):

fi:

od:

if OK = TRUE then

CNT := 0:

OK := TRUE:

# Step 1

APP := 0:

I1 := 1:

TOL[I1] := 10*EPS:

A[I1] := AA:

H[I1] := 0.5*(BB-AA):

FA[I1] := F(AA):

CNT := CNT+1:

FC[I1] := F((AA+H[I1])):

CNT := CNT+1:

FB[I1] := F(BB):

CNT := CNT+1:

# Approximation from Simpson's method for entire interval

S[I1] := H[I1]*(FA[I1]+4*FC[I1]+FB[I1])/3:

L[I1] := 1:

# Step 2

while I1 > 0 and OK = TRUE do

# Step 3

FD := F((A[I1]+0.5*H[I1])):

CNT := CNT+1:

FE := F((A[I1]+1.5*H[I1])):

CNT := CNT+1:

# Approximations from Simpson's method for halves of intervals

S1 := H[I1]*(FA[I1]+4*FD+FC[I1])/6:

S2 := H[I1]*(FC[I1]+4*FE+FB[I1])/6:

# Save data at this level

V[0] := A[I1]:

V[1] := FA[I1]:

V[2] := FC[I1]:

V[3] := FB[I1]:

V[4] := H[I1]:

V[5] := TOL[I1]:

V[6] := S[I1]:

LEV := L[I1]:

# Step 4

# Delete the level

I1 := I1-1:

# Step 5

if abs(S1+S2-V[6]) < V[5] then

APP := APP+(S1+S2):

else

if LEV >= N then

OK := FALSE:

# Procedure fails

else

# Add one level

# Data for right half subinterval

I1 := I1+1:

A[I1] := V[0]+V[4]:

FA[I1] := V[2]:

FC[I1] := FE:

FB[I1] := V[3]:

H[I1] := 0.5*V[4]:

TOL[I1] := 0.5*V[5]:

S[I1] := S2:

L[I1] := LEV+1:

I1 := I1+1:

# Data for left half subinterval

A[I1] := V[0]:

FA[I1] := V[1]:

FC[I1] := FD:

FB[I1] := V[2]:

H[I1] := H[I1-1]:

TOL[I1] := TOL[I1-1]:

S[I1] := S1:

L[I1] := L[I1-1]:

fi:

od:

if OK = FALSE then

print(`Level exceeded. Method failed to produce accurate approximation.`):

else

OUP := default:

fprintf(OUP,`Adaptive Quadrature\134n`):

fprintf(OUP,`\134nThe integral of F from %12.8f to %12.8f is\134n`, AA, BB):

fprintf(OUP,`%12.8f to within %14.8e\134n`, APP, EPS):

fprintf(OUP,`The number of function evaluations is: %d\134n`, CNT):

fi:

SVNUaGlzfmlzfkFkYXB0aXZlflF1YWRyYXR1cmV+d2l0aH5TaW1wc29uJ3N+TWV0aG9kLkc2Ig==

SUZJbnB1dH50aGV+ZnVuY3Rpb25+Rih4KX5pbn50ZXJtc35vZn54RzYi

STRGb3J+ZXhhbXBsZTp+Y29zKHgpRzYi

SShGKHgpfj1+RzYi

LCQqJi1JJHNpbkc2JCUqcHJvdGVjdGVkR0koX3N5c2xpYkc2IjYjLCQqJEkieEdGKSEiIiIjNSIiIkYtISIjIiQrIg==

SUZJbnB1dH5sb3dlcn5saW1pdH5vZn5pbnRlZ3JhdGlvbn5hbmR+RzYi

STt1cHBlcn5saW1pdH5vZn5pbnRlZ3JhdGlvbkc2Ig==

STVzZXBhcmF0ZWR+Ynl+YX5ibGFua0c2Ig==

SSVhfj1+RzYi

JCIiIiIiIQ==

SSVifj1+RzYi

JCIiJCIiIQ==

STFJbnB1dH50b2xlcmFuY2UuRzYi

SS1Ub2xlcmFuY2V+PX5HNiI=

JCIiIiEiJQ==

SURJbnB1dH50aGV+bWF4aW11bX5udW1iZXJ+b2Z+bGV2ZWxzLkc2Ig==

STxNYXhpbXVtfm51bWJlcn5vZn5sZXZlbHN+PX5HNiI=

IiM2

Adaptive Quadrature The integral of F from 1.00000000 to 3.00000000 is -1.42601481 to within 1.00000000e-04 The number of function evaluations is: 93

JSFH