2.8

2.8 | 6 ratings Rate this file 16 Downloads (last 30 days) File Size: 1.51 KB File ID: #8067

Gauss-Laguerre

by Jordi Soler Penades

 

19 Jul 2005 (Updated 19 Jul 2005)

Integrates a functio using Gauss-Laguerre cuadrature method

| Watch this File

File Information
Description

One function produces the Laguerre polynomial and the other integrates

MATLAB release MATLAB 6.5 (R13)
Tags for This File  
Everyone's Tags
Tags I've Applied
Add New Tags Please login to tag files.
Comments and Ratings (9)
25 Jan 2006 hanna mermelstein

The idea is quite nice, but the procedure does not work properly. From n >35 the zeros of the Laguerre polynomals start to be complex, what they should not be. Seems, that the roots-routine fails at this point. Up to n=32 the zeros are correctly calculated (see Krylov, Approximate calculation of integrals, Macmillan company, new york : 1962). But the weights (called coefficients in gausslaguerr.m) are completely wrong, thus the integration does not compute the right values at all.

25 Jan 2006 hanna mermelstein

The idea is quite nice, but the procedure does not work properly. From n >35 the zeros of the Laguerre polynomals start to be complex, what they should not be. Seems, that the roots-routine fails at this point. Up to n=32 the zeros are correctly calculated (see Krylov, Approximate calculation of integrals, Macmillan company, new york : 1962). But the weights (called coefficients in gausslaguerr.m) are completely wrong, thus the integration does not compute the right values at all.

11 Jun 2006 zarelly orellana

hola deseo mucha informacion sobre algoritmos de laguerre enmatlab

28 Sep 2007 Teddy Qian

I fixed some mistakes,now it's weights are right,just n<32.
if n>32,the problem is still there.
gausslaguerr.m
function [a,b]=gausslaguerr(n)
p=polelague(n);
x=roots(p(n+1,:));
x=x(x~=0);
for k=1:n
    q(k)=(factorial(n))^2*x(k)*exp(x(k))/(polyval(p(n+2,:),x(k)))^2;
end
a=x;
b=q;

polelague.m
function p=polelague(n)
p(1,1)=1;
p(2,1:2)=[-1 1];
for k=2:n+1
   p(k+1,1:k+1)=((2*(k-2)*[0 p(k,1:k)]+3*[0 p(k,1:k)]-[p(k,1:k) 0]-(k-1).^2*[0 0 p(k-1,1:k-1)]));
end

28 Jul 2008 Geert V.D.

The problem with the complex roots is caused by the the 'roots' function of Matlab. This determines the roots as the eigenvalues of the polynomials companion matrix 'CM', wich is such that det(x*I-CM)=P_n(x). To this aim, the general form of the companion matrix is used.
However, when the degree of the polynomial becomes 'too' large, the companion matrix is badly conditioned and yields inaccurate (in this case complex) roots.
The solution to this problem is hence to write your own 'roots' function, with a feasible companion matrix. Since we want the roots to be real, we would like the companion matrix to be symmetrical, as we know that all eigenvalues of a symmetrical matrix are real.
Below is the code for the Gauss-Laguerre, Gauss-Legendre and Gauss-Hermite quadratures. As you can see, the code consists of 2 blocks: first construct a symmetrical companion matrix and then determine the (real) eigenvalues (i.e. the roots of the polynomial). All three functions produce the correct abscissas and weights, for any value n>=2.

function [x, w] = GaussLegendre_2(n)

i = 1:n-1;
a = i./sqrt(4*i.^2-1);
CM = diag(a,1) + diag(a,-1);

[V L] = eig(CM);
[x ind] = sort(diag(L));
V = V(:,ind)';
w = 2 * V(:,1).^2;

function [x, w] = GaussLaguerre_2(n, alpha)

i = 1:n;
a = (2*i-1) + alpha;
b = sqrt( i(1:n-1) .* ((1:n-1) + alpha) );
CM = diag(a) + diag(b,1) + diag(b,-1);

[V L] = eig(CM);
[x ind] = sort(diag(L));
V = V(:,ind)';
w = gamma(alpha+1) .* V(:,1).^2;

function [x, w] = GaussHermite_2(n)

i = 1:n-1;
a = sqrt(i/2);
CM = diag(a,1) + diag(a,-1);

[V L] = eig(CM);
[x ind] = sort(diag(L));
V = V(:,ind)';
w = sqrt(pi) * V(:,1).^2;

23 Aug 2008 K C

Geert VD's Code is efficient using the companion matrix.
As a rule of thumb with matlab, you should always avoid using the loops if you can use vectorized methods.

17 May 2009 Johannes Ganze

Seems that all the Gauss type quadrature methods fail at large n ,for gauss-legendre quadrature the critical point is approximately n=50.

13 Jun 2011 Süleyman Ördek

götverenler

23 Aug 2011 Elleesse

What is the input 'alpha' in Geert V.D.'s approach to computing the Gauss-Laguerre quadrature?

function [x, w] = GaussLaguerre_2(n, alpha)

i = 1:n;
a = (2*i-1) + alpha;
b = sqrt( i(1:n-1) .* ((1:n-1) + alpha) );
CM = diag(a) + diag(b,1) + diag(b,-1);

[V L] = eig(CM);
[x ind] = sort(diag(L));
V = V(:,ind)';
w = gamma(alpha+1) .* V(:,1).^2;

Your response would be greatly appreciated.

Please login to add a comment or rating.
Tag Activity for this File
Tag Applied By Date/Time
integration Jordi Soler Penades 22 Oct 2008 07:53:45
gausslaguerre Jordi Soler Penades 22 Oct 2008 07:53:45
mathematics Jordi Soler Penades 22 Oct 2008 07:53:45
method Jordi Soler Penades 22 Oct 2008 07:53:45
integration Morison Cohen 13 Mar 2011 07:35:18
gausslaguerre younes 23 May 2011 12:17:57
gausslaguerre Süleyman Ördek 13 Jun 2011 15:32:32
integration Süleyman Ördek 13 Jun 2011 15:32:36
gausslaguerre Isaac 15 Dec 2011 22:13:14
integration Isaac 15 Dec 2011 22:13:17

Contact us at files@mathworks.com