% 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.
syms('OK', 'A', 'B', 'M', 'N', 'r', 'co', 'H1', 'H2', 'AJ', 'I', 'X');
syms('JX', 'C1', 'D1', 'K1', 'K2', 'J', 'Y', 'Q','x','s','y');
TRUE = 1;
FALSE = 0;
fprintf(1,'This is Gaussian Quadrature for double integrals.\n');
fprintf(1,'Input the function F(x,y) in terms of x and y\n');
fprintf(1,'For example: sqrt(x^2+y^2)\n');
s = input(' ','s');
F = inline(s,'x','y');
fprintf(1,'Input the functions C(x), and D(x) in terms of x ');
fprintf(1,'on separate lines\n');
fprintf(1,'For example: cos(x) \n');
fprintf(1,' sin(x) \n');
s = input(' ','s');
C = inline(s,'x');
s = input(' ','s');
D = inline(s,'x');
OK = FALSE;
while OK == FALSE
fprintf(1,'Input lower limit of integration and upper limit of\n ');
fprintf(1,'integration on separate lines\n');
A = input(' ');
B = input(' ');
if A > B
fprintf(1,'Lower limit must be less than upper limit\n');
else
OK = TRUE;
end
end
OK = FALSE;
while OK == FALSE
fprintf(1,'Input two integers M > 1 and N > 1 on separate lines.\n');
fprintf(1,'This implementation of Gaussian quadrature requires\n');
fprintf(1,'both to be less than or equal to 5.\n');
fprintf(1,'M is used for the outer integral and N for the inner\n');
fprintf(1,'integral.\n');
M = input(' ');
N = input(' ');
if M <= 1 | N <= 1
fprintf(1,'Integers must be greater than 1.\n');
else
if M > 5 | N > 5
fprintf(1,'Integers must be less than or equal to 5.\n');
else
docsity.com
