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.print(`This is Gaussian Quadrature for triple integrals. `): print(`Input the function F(x,y,z) in terms of x, y, and z `):print(`For example: z*sqrt(x^2+y^2) `):F := scanf(`%a`)[1]:print(`F(x,y,z) = `):print(F):F := unapply(F,x,y,z):print(`Input the functions alpha(x,y) and beta(x,y) in `):print(`terms of x and y, separated by a space `): alpha := scanf(`%a`)[1]:print(`alpha(x,y) = `):print(alpha):alpha := unapply(alpha,x,y):beta := scanf(`%a`)[1]:print(`beta(x,y) = `):print(beta):beta := unapply(beta,x,y):print(`Input the functions C(x), and D(x) in terms of x `):print(`separated by a space `):print(`For example: cos(x) sin(x) `):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 doprint(`Input lower limit of integration and upper limit of `): print(`integration separated by a blank space `):A := scanf(`%f`)[1]:B := scanf(`%f`)[1]:print(`a = `):print(A):print(`b = `):print(B):if A > B thenprint(`Lower limit must be less than upper limit `):elseOK := TRUE:fi:od: OK := FALSE:while OK = FALSE doprint(`Input three integers M > 1, N > 1, P > 1. This `):print(`implementation of Gaussian quadrature requires `):print(`that they are less than or equal to 5. `):print(`They will be used in `): print(`first, second, and third dimensions, resp. `):print(` - separate with blank. `):M := scanf(`%d`)[1]:N := scanf(`%d`)[1]:P := scanf(`%d`)[1]:print(`M = `):print(M):print(`N = `):print(N):print(`P = `):print(P):if M <= 1 or N <= 1 or P <= 1 then print(`Integers must be larger than 1. `):elseif M > 5 or N > 5 or P > 5 thenprint(`Integers must be less than or equal to 5. `):elseOK := TRUE:fi:fi:od:if OK = TRUE thenr[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 1H1 := (B-A)/2:H2 := (B+A)/2:# Use AJ in place of JAJ := 0:# Step 2for I1 from 1 to M do# Step 3X := H1*r[M-1,I1-1]+H2:JX := 0:C1 := C(X):D2 := D1(X):K1 := (D2-C1)/2:K2 := (D2+C1)/2:# Step 4for J from 1 to N do# Step 5Y := K1*r[N-1,J-1]+K2:JY := 0:# Use Z1 for BETA and Z2 for ALPHAZ1 := beta(X,Y):Z2 := alpha(X,Y):L1 := (Z1-Z2)/2:L2 := (Z1+Z2)/2:# Step 6for K from 1 to P doZ := L1*r[P-1,K-1]+L2:Q := F(X,Y,Z):JY := JY+co[P-1,K-1]*Q:od:# Step 7JX := JX+co[N-1,J-1]*L1*JY:od:# Step 8AJ := AJ+co[M-1,I1-1]*K1*JX:od:# Step 9AJ := AJ*H1:# Step 10OUP := default;fprintf(OUP,`Triple Integral Algorithm using Gaussian Quadrature \134n`):fprintf(OUP,`\134nThe triple integral of F from %12.8f to %12.8f is `, A, B):fprintf(OUP,`%16.11f\134n`, AJ):fprintf(OUP,` obtained with M = %3d , N = %3d and P = %3d\134n`, M, N, P):fi:JSFH