clear all
a=-3;
b=2;
%xpoints=linspace(a,b,11)’;
%[xpoints,w]=lgwt(21,-1,1);
[xpoints,w,P]=lglnodes(51);
tpoints=((b-a)/2)*xpoints+(a+b)/2;
fpoints=exp(tpoints).*sin(5*tpoints);
n=length(xpoints)-1; %order of polynomial
xint=-1:0.01:1;
tint=((b-a)/2)*xint+(a+b)/2;
fint=LagrangePolynomial(tpoints,fpoints,n,tint);
freal=exp(tint).*sin(5*tint);
dfreal=exp(tint).*(sin(5*tint)+5*cos(5*tint));
figure(1)
hold off
plot(tpoints,fpoints,’ko’,’MarkerSize’,20,’MarkerFaceColor’,’r’)
hold on
plot(tint,fint,’b-‘,’Linewidth’,4)
plot(tint,freal,’k-‘,’Linewidth’,4)
xlabel(‘t’);ylabel(‘f(t)’)
D=DerivMatrix(xpoints,n);
derivfpoints=(2/(b-a))*D*fpoints;
figure(2)
hold off
plot(tint,dfreal,’k-‘,’Linewidth’,4)
hold on
plot(tpoints,derivfpoints,’ko’,’MarkerSize’,20,’MarkerFaceColor’,’b’)
xlabel(‘t’);ylabel(‘df(t)/dt’)
function D=DerivMatrix(x,n)
D=zeros(n+1,n+1);
for i=1:n+1
num(i)=1.0;
for k=1:n+1
if k ~= i
num(i)=num(i)*(x(i)-x(k));
end
end
end
for j=1:n+1
den(j)=1.0;
for k=1:n+1
if k ~= j
den(j)=den(j)*(x(j)-x(k));
end
end
end
for i=1:n+1
for j=1:n+1
if i ~= j
D(i,j)=num(i)/(den(j)*(x(i)-x(j)));
end
end
end
for i=1:n+1
for k=1:n+1
if i~=k
D(i,i)=D(i,i)+1./(x(i)-x(k));
end
end
end
end
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,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