clear all
Npoints=50;
alpha = 5; %BC1 y(x=-1)=5
beta = 3; %BC2 y(x=1)=3
%xpoints=linspace(-1,1,Npoints)’;
[xpoints,w,P]=lglnodes(Npoints);
xpoints=-xpoints;
tpoints=(xpoints+3)/2;
y = @(t) t+4./t.^2;
yexact = y(tpoints);
n=length(xpoints)-1; %order of polynomial
D=DerivMatrix(xpoints,n);
D2 = D*D;
A=D2+diag(2./(xpoints+3))*D-diag(2./(xpoints+3).^2);
Amatrix=A(2:end-1,2:end-1);
Cvector=-alpha*A(2:end-1,1)-beta*A(2:end-1,end);
sol = Amatrix\Cvector;
yapprox = [alpha;sol;beta];
hold off
plot(tpoints,yapprox,’bo’)
hold on
plot(tpoints,yexact,’k-‘)
xlabel(‘t’);ylabel(‘y’);
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