> restart;
> # GAUSSIAN DOUBLE INTEGRAL ALGORITHM 4.5 
> #
> # 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; positive integers m, n. (Assume that the
> #           roots r(i,j) and coefficients c(i,j) are available for 
> #           i equals m and n for 1<= j <= i.
> # 
> # OUTPUT:   approximation J to I. 
> alg045 := proc() local F, C, D, OK, A, B, M, N, r, co, H1, H2, AJ, I, X, JX, C1, D1, K1, K2, J, Y, Q;
> printf(`This is Gaussian Quadrature for double integrals.\n`);  
> printf(`Input the function F(x,y) in terms of x and y\n`);
> printf(`For example: sqrt(x^2+y^2)\n`);
> F := scanf(`%a`)[1];
> F := unapply(F,x,y);
> printf(`Input the functions C(x), and D(x) in terms of x `);
> printf(`separated by a space\n`);
> printf(`For example: cos(x) sin(x)\n`);
> 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 upper limit of `); 
> printf(`integration separated by a blank space\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 integers M > 1 and N > 1. This\n`);
> printf(`implementation of Gaussian quadrature requires\n`);
> printf(`both to be less than or equal to 5.\n`);
> printf(`M is used for the outer integral and N for the inner integral - separated by a space.\n`);
> M := scanf(`%d`)[1];
> N := scanf(`%d`)[1];
> if M <= 1 or N <= 1 then 
> printf(`Integers must be larger than 1.\n`);
> else
> if M > 5 or N > 5 then
> printf(`Integers must be less than or equal to 5.\n`);
> else
> OK := TRUE;
> fi;
> fi;
> od;
> if OK = TRUE then
> r[1,0] := 0.5773502692; r[1,1] := -r[1,0]; co[1,0] := 1.0;
> co[1,1] := 1.0; r[2,0] := 0.7745966692; r[2,1] := 0.0;
> r[2,2] := -r[2,0]; co[2,0] := 0.5555555556; co[2,1] := 0.8888888889;
> co[2,2] := co[2,0]; r[3,0] := 0.8611363116; r[3,1] := 0.3399810436;
> r[3,2] := -r[3,1]; r[3,3] := -r[3,0]; co[3,0] := 0.3478548451;
> co[3,1] := 0.6521451549; co[3,2] := co[3,1]; co[3,3] := co[3,0];
> r[4,0] := 0.9061798459; r[4,1] := 0.5384693101; r[4,2] := 0.0;
> r[4,3] := -r[4,1]; r[4,4] := -r[4,0]; co[4,0] := 0.2369268850;
> co[4,1] := 0.4786286705; co[4,2] := 0.5688888889; co[4,3] := co[4,1];
> co[4,4] := co[4,0];
> # STEP 1 
> H1 := (B-A)/2;
> H2 := (B+A)/2;
> # use AJ in place of J 
> AJ := 0;
> # STEP 2 
> for I from 1 to M do
> # STEP 3 
> X := H1*r[M-1,I-1]+H2;
> JX := 0;
> C1 := C(X);
> D1 := D(X);
> K1 := (D1-C1)/2;
> K2 := (D1+C1)/2;
> # STEP 4 
> for J from 1 to N do
> Y := K1 * r[N-1,J-1]+K2;
> Q := F(X, Y);
> JX := JX + co[N-1,J-1]*Q;
> od;
> # STEP 5 
> AJ := AJ+co[M-1,I-1]*K1*JX;
> od;
> # STEP 6 
> AJ := AJ*H1;
> # STEP 7 
> printf(`\nThe double integral of F from %12.8f to %12.8f is\n`, A, B);
> printf(`  %.10e`, AJ);
> printf(` obtained with M = %3d and N = %3d\n`, M, N);
> fi;
> RETURN(0);
> end;
> alg045();