File Exchange

image thumbnail

Legendre Laguerre and Hermite - Gauss Quadrature

version 1.2 (3.25 KB) by

Nodes and weights for Legendre Laguerre and Hermite - Gauss Quadrature

4.85714
10 Ratings

20 Downloads

Updated

View License

This .zip file contains 3 mfiles for computing the nodes and weights for Legendre Laguerre and Hermite - Gauss Quadrature of any order n.

Contrary to most of the files in use, the included files are all based on a symmetrical companion matrix, the eigenvalues of which are always real and correspond to the roots of the respective polynomials.

Hence, there is a 100% certainty of avoiding complex roots at high orders. It is known that the latter issue often causes huge numerical troubles.

Comments and Ratings (15)

Jakob Ameres

very useful!

Leo Lee

Can I ask a question? Can I use this for multivariate normal case?

brabbit

Osvaldo

Very clever and precise.

Is it me or these routines are not working correctly in Matlab 2013b?
PS. I have been using them since 2010.

Simple and fast.

afu2007 ??

nice!

David Holdaway

David Holdaway (view profile)

I added this on the end to refine the eigenvalues via newtons method

if mod(n,2)==1
x(ceil(n/2))=0;
end

z=x(x>=0);
success = false;
for its=1:50 %maxit=50, usually will take 2
p1 = pi^(-1/4);
p2=0;
for j=1:n %make hermite we need
p3=p2;
p2=p1;
p1=z.*sqrt(2/j).*p2-sqrt((j-1)/j).*p3;
end
pp = sqrt(2*n).*p2;
z1=z;
z=z1-p1./pp;
if all(abs(z-z1)< 20*eps)
success = true;
break
end
end
if ~success
warning('failed to converge to desired accuracy')
end
w(x>=0) = 2./pp.^2;
w(x<=0) = flipud(w(x>=0));
x(x>=0) = z;
x(x<=0) = -flipud(z);

David Holdaway

David Holdaway (view profile)

Oh and I forgot to add this required you insist x = 0 is the middle value for odd via
if mod(n,2)==1
x(ceil(n/2))=0;
end

David Holdaway

David Holdaway (view profile)

Very nice submission, simple and fast

One addition I did is this

% Linear map from[-1,1] to [a,b]
x=(a*(1-x)+b*(1+x))/2;
w=(b-a).*w/2;

a and b are user provided, and I think now that you can shift from [1 1] to [a b]

Very useful program.

How can I calculate the GaussLegendre in the [0 1] space?

Updates

1.2

Some minor changes to the code's comments.

1.1

No updates were made. Only the description was changed: it should read "the eigenvalues of the companion matrix" instead of "the roots of the companion matrix". Obviously, a matrix has no roots.

MATLAB Release
MATLAB 7.0.4 (R14SP2)

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video