Numerical integration gauss legendre

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

Answer 1

[Roger Stafford]

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.)

[candas]

Thanks a lot for your answer. I tried to make a code for computing an integral in [a,b]x[c,d] but when I'm running the code, the answer is exceeding the real value of the integral. I used the m-file lgwt.m in my code for the Legendre-Gauss nodes and weights. Can you see where is the error of this code?

function [I]=gl2d(f,a,b,c,d,n,m)

% %

[w1,x1]=lgwt(n,a,b);

[w2,x2]=lgwt(m,c,d);

if a==b | c==d

    I=0;

    return

end

I=0;

for i=1:n

    for j=1:m

        I=I+(((b-a)*(d-c))/4)*w1(i)*w2(j)*f(x1(i),x2(j));

    end

end

[Roger Stafford]

The function 'lgwt' which you are using is presumably the File Exchange contribution by Greg von Winckel which he wrote in 2004. According to his documentation it is supposed to be called by

 [x,w] = lgwt(N,a,b)

with the returned node positions first, followed by the node weights, whereas you have them in reversed order in your code:

 [w1,x1] = lgwt(n,a,b)
 [w2,x2] = lgwt(m,c,d)

That would certainly lead to an erroneous answer.

The next problem is the scale factor. In the form I gave you my assumption was that the weights that would be furnished would be those corresponding to an interval from -1 to +1 which is where the roots of the Legendre polynomials are all located and that they would consequently add up to 2. However the weights which 'lgwt' returns are already rescaled to those appropriate to quadrature over the interval from a to b. Hence the corrective factor I gave, (b-a)*(c-d)/4, should not be applied here with 'lgwt'. You just need the sum of the following over all i and j:

 w1(i)*w2(j)*f(x1(i),x2(j))

As a final check you should verify that sum(w1) equals b-a and sum(w2) equals d-c (to within roundoff error of course,) and that the positions of x1 and x2 are located appropriately within the respective intervals [a,b] and [c,d].

[candas]

Thank you so much! I can do what you 've told me or just change

[x1,w1] = lgwt(n,a,b)

[x2,w2] = lgwt(m,c,d)

 to

[x1,w1] = lgwt(n,-1,1)

[x2,w2] = lgwt(m,-1,1)

and it's the same. Thanks again!