> restart;
> # GAUSSIAN TRIPLE INTEGRAL ALGORITHM 4.6
> # 
> # To approximate I = triple integral ( ( f(x,y,z) dz dy dx ) ) with       
> # limits of integration from a to b for x, from c(x) to d(x) for y, 
> # and from alpha(x,y) to beta(x,y) for z.
> #
> # INPUT:   endpoints a, b; positive integers m, n, p. (Assume that the
> #          roots r(i,j) and coefficients c(i,j) are available for i
> #          equals m, n, and p and for 1 <= j <= i.
> alg046 := proc() local F, C, D, alpha, beta, OK, A, B, M, N, P, r, co, H1, H2, AJ, I, X, JX, C1, D1, K1, K2, J, Y, JY, Z1, Z2, L1, L2, K, Z, Q;
> printf(`This is Gaussian Quadrature for triple integrals.\n`);  
> printf(`Input the function F(x,y,z) in terms of x, y, and z\n`);
> printf(`For example: z*sqrt(x^2+y^2)\n`);
> F := scanf(`%a`)[1];
> F := unapply(F,x,y,z);
> 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);
> printf(`Input the functions alpha(x,y) and beta(x,y) in `);
> printf(`terms of x and y, separated by a space\n`); 
> alpha := scanf(`%a`)[1];
> alpha := unapply(alpha,x,y);
> beta := scanf(`%a`)[1];
> beta := unapply(beta,x,y);
> 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 three integers M > 1, N > 1, P > 1. This\n`);
> printf(`implementation of Gaussian quadrature requires\n`);
> printf(`that they are less than or equal to 5.\n`);
> printf(`They will be used in `);  
> printf(`first, second, and third dimensions, resp.\n`);
> printf(` - separate with blank.\n`);
> M := scanf(`%d`)[1];
> N := scanf(`%d`)[1];
> P := scanf(`%d`)[1];
> if M <= 1 or N <= 1 or P <= 1 then 
> printf(`Integers must be larger than 1.\n`);
> else
> if M > 5 or N > 5 or P > 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
> # Step 5
> Y := K1*r[N-1,J-1]+K2;
> JY := 0;
> # Use Z1 for BETA and Z2 for ALPHA
> Z1 := beta(X,Y);
> Z2 := alpha(X,Y);
> L1 := (Z1-Z2)/2;
> L2 := (Z1+Z2)/2;
> # Step 6
> for K from 1 to P do
> Z := L1*r[P-1,K-1]+L2;
> Q := F(X,Y,Z);
> JY := JY+co[P-1,K-1]*Q;
> od;
> # Step 7
> JX := JX+co[N-1,J-1]*L1*JY;
> od;
> # Step 8
> AJ := AJ+co[M-1,I-1]*K1*JX;
> od;
> # Step 9
> AJ := AJ*H1;
> # Step 10
> printf(`\nThe triple integral of F from %12.8f to %12.8f is   `, A, B);
> printf(`%.10e\n`, AJ);
> printf(` obtained with M = %3d , N = %3d and P = %3d\n`, M, N, P);
> fi;
> RETURN(0);
> end;
> alg046();