restart:

# DOUBLE INTEGAL ALGORITHM 4.4

# To approximate I = double integral ( ( f(x,y) dy dx ) ) with limits

# of integration from a to b for x and from c(x) to d(x) for y:

# INPUT: endpoints a, b: even positive integers m, n.

# OUTPUT: approximation J to I.

print(`This is Simpsons Method for double integrals. `):

print(`Input the functions F(X,Y), C(X), and D(X) in terms of x and y separated by a space. `):

print(`For example: cos(x+y) x^3 x `):

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

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

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

C := unapply(C,x):

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

D1 := unapply(D1,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 two even positive integers N, M : there will`):

print(` be N subintervals for outer `):

print(`integral and M subintervals for inner `):

print(`integral - separate with blank `):

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

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

if M > 0 and N > 0 and M mod 2 = 0 and N mod 2 = 0 then

OK := TRUE:

else

print(`Integers must be even and positive `):

fi:

od:

if OK = TRUE then

# Step 1

H := (B-A)/N:

# Use AN, AE, AO for J(1), J(2), J(3) resp.

# End terms

AN := 0:

# Even terms

AE := 0:

# Odd terms

AO := 0:

# Step 2

for I1 from 0 to N do

# Step 3

# Composite Simpson's method for X

X := A+I1*H:

YA := C(X):

YB := D1(X):

HX := (YB-YA)/M:

# Use BN, BE, BO for K(1), K(2), K(3) resp.

# End terms

BN := F(X,YA)+F(X,YB):

# Even terms

BE := 0:

# Odd terms

BO := 0:

# Step 4

for J from 1 to M-1 do

# Step 5

Y := YA+J*HX:

Z := F(X, Y):

# Step 6

if J mod 2 = 0 then

BE := BE+Z:

else

BO := BO+Z:

fi:

od:

# Step 7

# Use A1 for L, which is the integral of F(X(I),Y) from

# C(X(I)) to D(X(I)) by the Composite Simpson's method

A1 := (BN+2*BE+4*BO)*HX/3:

# Step 8

if I1 = 0 or I1 = N then

AN := AN+A1:

else

if I1 mod 2 = 0 then

AE := AE + A1:

else

AO := AO + A1:

fi:

od:

# Step 9

# Use AC of J

AC := (AN + 2 * AE + 4 * AO) * H /3:

# Step 10

OUP := default:

fprintf(OUP,`Double Integrals Using Simpson's Composite Method \134n`):

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

fprintf(OUP,`%12.8f`, AC):

fprintf(OUP,` obtained with N := %3d and M := %3d\134n`, N, M):

fi:

SVBUaGlzfmlzflNpbXBzb25zfk1ldGhvZH5mb3J+ZG91YmxlfmludGVncmFscy5+fkc2Ig==

SWBwSW5wdXR+dGhlfmZ1bmN0aW9uc35GKFgsWSksfkMoWCksfmFuZH5EKFgpfmlufnRlcm1zfm9mfnh+YW5kfnl+c2VwYXJhdGVkfmJ5fmF+c3BhY2Uufkc2Ig==

ST5Gb3J+ZXhhbXBsZTp+Y29zKHgreSl+eF4zfnh+fkc2Ig==

SSpGKHgseSl+PX5HNiI=

LUkkZXhwRzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliRzYiNiMqJkkieUdGJyIiIkkieEdGJyEiIg==

SShDKHgpfj1+RzYi

KiRJInhHNiIiIiQ=

SShEKHgpfj1+RzYi

KiRJInhHNiIiIiM=

SUZJbnB1dH5sb3dlcn5saW1pdH5vZn5pbnRlZ3JhdGlvbn5hbmR+RzYi

STx1cHBlcn5saW1pdH5vZn5pbnRlZ3JhdGlvbn5HNiI=

STZzZXBhcmF0ZWR+Ynl+YX5ibGFua35HNiI=

SSVhfj1+RzYi

JCIiIiEiIg==

SSVifj1+RzYi

JCIiJiEiIg==

SVNJbnB1dH50d29+ZXZlbn5wb3NpdGl2ZX5pbnRlZ2Vyc35OLH5Nfjp+dGhlcmV+d2lsbEc2Ig==

ST5+YmV+Tn5zdWJpbnRlcnZhbHN+Zm9yfm91dGVyfkc2Ig==

SUdpbnRlZ3JhbH5hbmR+TX5zdWJpbnRlcnZhbHN+Zm9yfmlubmVyfkc2Ig==

SUBpbnRlZ3JhbH4tfnNlcGFyYXRlfndpdGh+Ymxhbmt+RzYi

SSVOfj1+RzYi

IiM1

SSVNfj1+RzYi

IiM1

Double Integrals Using Simpson's Composite Method The double integral of F from 0.10000000 to 0.50000000 is 0.03330546 obtained with N := 10 and M := 10

JSFH