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. print(`This is Gaussian Quadrature for double integrals. `): print(`Input the function F(x,y) in terms of x and y `):print(`For example: sqrt(x^2+y^2) `):F := scanf(`%a`)[1]:print(`F(x,y) = `):print(F):F := unapply(F,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 two integers M > 1 and N > 1. This `):print(`implementation of Gaussian quadrature requires `):print(`both to be less than or equal to 5. `):print(`M is used for the outer integral and N for the inner integral - separated by a space. `):M := scanf(`%d`)[1]:N := scanf(`%d`)[1]:print(`M = `):print(M):print(`N = `):print(N):if M <= 1 or N <= 1 then print(`Integers must be larger than 1. `):elseif M > 5 or N > 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 1 H1 := (B-A)/2:H2 := (B+A)/2:# use AJ in place of J AJ := 0:# STEP 2 for I1 from 1 to M do# STEP 3 X := H1*r[M-1,I1-1]+H2:JX := 0:C1 := C(X):D2 := D1(X):K1 := (D2-C1)/2:K2 := (D2+C1)/2:# STEP 4 for J from 1 to N doY := 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,I1-1]*K1*JX:od:# STEP 6 AJ := AJ*H1:# STEP 7 OUP := default:fprintf(OUOP,`Double Integral Algorithm using Gaussian Quadrature \134n`):fprintf(OUP,`\134nThe double integral of F from %12.8f to %12.8f is\134n`, A, B):fprintf(OUP,` %16.11e`, AJ):fprintf(OUP,` obtained with M = %3d and N = %3d\134n`, M, N):fi:JSFH