clear all
Nmax=40;
counter=1;
for i=3:Nmax
Npoints=i;
Npointsarray(counter)=Npoints;
xpoints=linspace(-1,1,Npoints)’;
%[xpoints,w]=lgwt(Npoints,-1,1);
%[xpoints,w,P]=lglnodes(Npoints);
ypoints=exp(xpoints).*sin(5*xpoints);
n=length(xpoints)-1; %order of polynomial
D=DerivMatrix(xpoints,n);
dfexact=exp(xpoints).*(sin(5*xpoints)+5*cos(5*xpoints));
dfapprox=D*ypoints;
l2norm(counter)=(1./Npoints)*sum(abs(dfexact-dfapprox));
counter=counter+1;
end
plot(Npointsarray,log10(l2norm),’k-‘,’LineWidth’,3)
axis([0 Nmax -15 2])
xlabel(‘N’)
ylabel(‘log_{10}(error)’)
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 dyint=DerivLagrangePolynomial(x,y,n,xint)
dyint=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
temp=zeros(size(xint));
for j=1:n+1
if j ~= i
temp=temp+1./(xint-x(j));
end
end
dyint=dyint+y(i).*(Lix.*temp);
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