clear all
xpoints=linspace(-1,1,5);
%[xpoints,w]=lgwt(25,-1,1);
%[xpoints,w,P]=lglnodes(2);
ypoints=1./(1+25*xpoints.^2);
n=length(xpoints)-1; %order of polynomial
xint=-1:0.01:1;
yint=LagrangePolynomial(xpoints,ypoints,n,xint);
yreal=1./(1+25*xint.^2);
hold off
plot(xpoints,ypoints,’ko’,’MarkerSize’,20,’MarkerFaceColor’,’r’)
hold on
plot(xint,yint,’b-‘,’Linewidth’,4)
plot(xint,yreal,’k-‘,’Linewidth’,4)
xlabel(‘x’)
ylabel(‘y’)
function yint=LagrangePolynomial(x,y,n,xint)
yint=zeros(size(xint));
for i=1:n+1
Lix=1.0;
for j=1:n+1
if j ~= i
Lix=Lix.*(xint-x(j))/(x(i)-x(j));
end
end
yint=yint+y(i).*Lix;
end
end
function [x,w]=lgwt(N,a,b)
% lgwt.m
%
% This script is for computing definite integrals using Legendre-Gauss
% Quadrature. Computes the Legendre-Gauss nodes and weights on an interval
% [a,b] with truncation order N
%
% Suppose you have a continuous function f(x) which is defined on [a,b]
% which you can evaluate at any x in [a,b]. Simply evaluate it at all of
% the values contained in the x vector to obtain a vector f. Then compute
% the definite integral using sum(f.*w);
%
% Written by Greg von Winckel – 02/25/2004
% This MATLAB Function is taken from MATLAB Central File Exchange
N=N-1;
N1=N+1; N2=N+2;
xu=linspace(-1,1,N1)’;
% Initial guess
y=cos((2*(0:N)’+1)*pi/(2*N+2))+(0.27/N1)*sin(pi*xu*N/N2);
% Legendre-Gauss Vandermonde Matrix
L=zeros(N1,N2);
% Derivative of LGVM
Lp=zeros(N1,N2);
% Compute the zeros of the N+1 Legendre Polynomial
% using the recursion relation and the Newton-Raphson method
y0=2;
% Iterate until new points are uniformly within epsilon of old points
while max(abs(y-y0))>eps
L(:,1)=1;
Lp(:,1)=0;
L(:,2)=y;
Lp(:,2)=1;
for k=2:N1
L(:,k+1)=( (2*k-1)*y.*L(:,k)-(k-1)*L(:,k-1) )/k;
end
Lp=(N2)*( L(:,N1)-y.*L(:,N2) )./(1-y.^2);
y0=y;
y=y0-L(:,N2)./Lp;
end
% Linear map from[-1,1] to [a,b]
x=(a*(1-y)+b*(1+y))/2;
% Compute the weights
w=(b-a)./((1-y.^2).*Lp.^2)*(N2/N1)^2;
end
function [x,w,P]=lglnodes(N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% lglnodes.m
%
% Computes the Legendre-Gauss-Lobatto nodes, weights and the LGL Vandermonde
% matrix. The LGL nodes are the zeros of (1-x^2)*P’_N(x). Useful for numerical
% integration and spectral methods.
%
% Reference on LGL nodes and weights:
% C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Tang, “Spectral Methods
% in Fluid Dynamics,” Section 2.3. Springer-Verlag 1987
%
% Written by Greg von Winckel – 04/17/2004
% Contact: gregvw@chtm.unm.edu
% This MATLAB Function is taken from MATLAB Central File Exchange
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Truncation + 1
N1=N+1;
% Use the Chebyshev-Gauss-Lobatto nodes as the first guess
x=cos(pi*(0:N)/N)’;
% The Legendre Vandermonde Matrix
P=zeros(N1,N1);
% Compute P_(N) using the recursion relation
% Compute its first and second derivatives and
% update x using the Newton-Raphson method.
xold=2;
while max(abs(x-xold))>eps
xold=x;
P(:,1)=1; P(:,2)=x;
for k=2:N
P(:,k+1)=( (2*k-1)*x.*P(:,k)-(k-1)*P(:,k-1) )/k;
end
x=xold-( x.*P(:,N1)-P(:,N) )./( N1*P(:,N1) );
end
w=2./(N*N1*P(:,N1).^2);
end