> 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. 
> alg044 := proc() local F, C, D, OK, A, B, N, M, NN, MM, H, AN, AE, AO, I, X, YA, YB, HX, BN, BE, BO, J, Y, Z, A1, AC;
> printf(`This is Simpsons Method for double integrals.\n\n`);
> printf(`Input the functions F(X,Y), C(X), and D(X) in terms of x and y separated by a space.\n`);
> printf(`For example: cos(x+y) x^3 x \n`);
> F := scanf(`%a`)[1];
> F := unapply(F,x,y);
> C := scanf(`%a`)[1];
> C := unapply(C,x);
> D := scanf(`%a`)[1];
> D := unapply(D,x);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input lower limit of integration and `);
> printf(`upper limit of integration\n`);
> printf(`separated by a blank\n`);
> A := scanf(`%f`)[1];
> B := scanf(`%f`)[1];
> if A > B then
> printf(`Lower limit must be less than upper limit\n`);
> else
> OK := TRUE;
> fi;
> od; 
> OK := FALSE;
> while OK = FALSE do
> printf(`Input two even positive integers N, M ; there will`);
> printf(` be N subintervals for outer\n`);
> printf(`integral and M subintervals for inner `);
> printf(`integral - separate with blank\n`);
> N := scanf(`%d`)[1];
> M := scanf(`%d`)[1];
> if M > 0 and N > 0 and M mod 2 = 0 and N mod 2 = 0 then
> OK := TRUE;
> else
> printf(`Integers must be even and positive\n`);
> 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 I from 0 to N do
> # Step 3
> # Composite Simpson's method for X
> X := A+I*H;
> YA := C(X);
> YB := D(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 I = 0 or I = N then
> AN := AN+A1;
> else
> if I mod 2 = 0 then
> AE := AE + A1;
> else
> AO := AO + A1;
> fi;
> fi;
> od;
> # Step 9
> # Use AC of J
> AC := (AN + 2 * AE + 4 * AO) * H /3;
> # Step 10
> printf(`\nThe double integral of F from %12.8f to %12.8f is\n`, A, B);
> printf(`%12.8f`, AC);
> printf(` obtained with N := %3d and M := %3d\n`, N, M);
> fi;
> RETURN(0);
> end;
> alg044();