Recently, I got a request how one can find the quadrature and weights of a Gauss-Legendre quadrature rule for large n. It seems that the internet has these points available free of charge only up to n=12. Below is the
MATLAB program that finds these values for any n. I tried the program for n=25 and it gave results in a minute or so. The results output up to 32 significant digits.
_______________________________________________________
% Program to get the quadrature points
% and weight for Gauss-Legendre Quadrature
% Rule
clc
clear all
syms x
% Input n: Quad pt rule
n=14;
% Calculating the Pn(x)
% Legendre Polynomial
% Using recursive relationship
% P(order of polynomial, value of x)
% P(0,x)=1; P(1,x)=0;
% (i+1)*P(i+1,x)=(2*i+1)*x*P(i,x)-i*P(i-1,x)
m=n-1;
P0=1;
P1=x;
for i=1:1:m
Pn=((2.0*i+1)*x*P1-i*P0)/(i+1.0);
P0=P1;
P1=Pn;
end
if n==1
Pn=P1;
end
Pn=expand(Pn);
quadpts=solve(vpa(Pn,32));
quadpts=sort(quadpts);
% Finding the weights
% Formula for weights is given at
% http://mathworld.wolfram.com/Legendre-GaussQuadrature.html
% Equation (13)
for k=1:1:n
P0=1;
P1=x;
m=n;
% Calculating P(n+1,x)
for i=1:1:m
Pn=((2.0*i+1)*x*P1-i*P0)/(i+1.0);
P0=P1;
P1=Pn;
end
Pn=P1;
weights(k)=vpa(2*(1-quadpts(k)^2)/(n+1)^2/ …
subs(Pn,x,quadpts(k))^2,32);
end
fprintf(‘Quad point rule for n=%g \n’,n)
disp(‘ ‘)
disp(‘Abscissas’)
disp(quadpts)
disp(‘ ‘)
disp(‘Weights’)
disp(weights’)_______________________________________________________
