Seems that all the Gauss type quadrature methods fail at large n ,for gauss-legendre quadrature the critical point is approximately n=50.
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.
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;
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
zarelly orellana
hola deseo mucha informacion sobre algoritmos de laguerre enmatlab
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.
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.
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Jose (view profile)
Hello everyone, can someone explain me how can I run this algorithm. I need help please I'm new in matlab
Elleesse (view profile)
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.
Süleyman Ördek (view profile)
götverenler
Johannes Ganze (view profile)
Seems that all the Gauss type quadrature methods fail at large n ,for gauss-legendre quadrature the critical point is approximately n=50.
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.
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;
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
hola deseo mucha informacion sobre algoritmos de laguerre enmatlab
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.
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.