MATLAB: Numerical integration gauss legendre

gausslegendre

I want a code which compute double integrals with gauss-legendre method.

Best Answer

There is a function in matlab's Symbolic Toolbox Mupad called "numeric::gldata" described at:

 http://www.mathworks.com/help/symbolic/mupad_ref/numeric-gldata.html

as well as a number of entries in the File Exchange which provide for the Gauss-Legendre method for single integrals. I would suppose you already know of these.

If your double integral is to be taken over a rectangular area, you can regard your problem as a single integral with respect to one variable of single integrals with respect to a second variable for some function of these two variables. The weights and abscissae of the Gauss-Legendre nodes in the one direction and those in the other direction would be selected in a grid pattern of points within the rectangular area. There is a pdf article at the site:

http://portal.tugraz.at/portal/page/portal/Files/i2020/files/Abschlussarbeiten/fischnaller09m.pdf

which explains this method as applied to double integrals better than I can here in their equation 2.7 on page 9. As explained in this article, your computation would be a sum of sums of products. Each term would the the product of the integrand value at a node multiplied by two weight values, one for each of the two directions. As you are no doubt aware, the Gauss-Legnedre method dictates the spacing of nodes where the integrand function is to be evaluated. This seems to me to be the correct approach to your problem.

Note that the above "sum of sums" can be expressed nicely in matlab by the product of three arrays and a scale factor:

 A/4*W1'*F*W2

where W1 and W2 are column vectors of the node weights in the two direction, where F is a matrix of appropriate values of the integrand at the nodes, and where the scalar A is the area of the rectangle.

(Corrected: The division by 4 above is needed because the Gauss-Legendre weights and abscissae are intended for the interval from -1 to +1 and the weights sum to 2, so for a double integral we need to divide the area by 4.)

Related Solutions

MATLAB: How to integrate a functional of integrals numerically

The integrals in the numerator and denominator of φ(z,j) are actually double integrals of ψ w.r.t. x and y, so I calculated them that way in the loop, substituting in the required values for z, and then computed φ(z,j) from them and beta. (The entire calculation required about 214 seconds on my machine.)

This code:

zv = 0.1:0.1:1;
jv = 0:31;
fpsixyz = @(x,y,z)  exp(-(((x-y).*z).^2)/4)./(1+y.^2);      % Function to integrate
B = @(j) 1E-5 * 2.^j;                                       % Beta
for k1 = 1:size(jv,2)
    for k2 = 1:size(zv,2)
        z = zv(k2);
        psiz = integral2(@(x,y) fpsixyz(x,y,z), 0, Inf, -Inf, Inf); % Integrate over [x,y]
        j = jv(k1);
        phizj(k1,k2) = psiz ./ (B(j) + psiz);       % Compute phi(z,j)
    end
end
figure(1)
surfc(phizj)
xlabel('\itz\rm')
ylabel('\itj\rm')
zlabel('\phi(\itz, j\rm)')
set(gca, 'XTick', 1:2:10, 'XTickLabel', zv(1:2:end))
axis([xlim    0  31    0  1])

produces:

MATLAB: Matlab Code for the Gauss Legendre Quadrature

Have your vectors of weights be the same size as your vectors of locations, ordered so that corresponding weights and locations are in the same position:

%% parameters
a = -1; % lower bound
b = 1; % upper bound 
n = [-1/sqrt(3) 1/sqrt(3)]; %location values for n=2
n1 = [-0.86113 -0.33998 0.33998 0.86113]; %location values for n=4
n2 = [-0.93246 -0.66120 -0.23861 0.23861 0.66120 0.993246]; % location values for n=6
w = [1 1]; %weight for n=2
w1 = [0.34785 0.625214 0.625214 0.34785]; %weights for n=4
w2 = [0.17132 0.36076 0.6791 0.6791 0.36076 0.17132]; % weights for n=6
f = @(x) (1)./(x.^2 .* (sqrt(x.^2 +1)));%function 
%% Gauss_Legendre Quadrature 
% n=2
gaussleg = sum(w.*f(n));
%n=4
gaussleg1 = sum(w1.*f(n1));
%n=6
gaussleg2 = sum(w2.*f(n2));

Higher precision in your weights and locations will give more accurate results.

(edit) I also changed f slightly so that you could pass in the entire vector of locations at once.

Best Answer

Related Solutions

MATLAB: How to integrate a functional of integrals numerically

MATLAB: Matlab Code for the Gauss Legendre Quadrature

Related Question