Code covered by the BSD License

Chebpack

Damian Trif (view profile)

15 Jul 2011 (Updated )

The MATLAB package Chebpack solves specific problems for differential or integral equations.

eig42F.m
function [lam,phi,x]=eig42F(n,dom,kind,numeigval)
% From: Paul T. Dawkins et.al., The Origin and Nature of Spurious Eigenvalues
% in the Spectral Tau Method, J. COMP. PHYS., 147, 441462 (1998)
%
% y''''=s y'', x in [-1,1]
% y(-1)=y(1)=0=y'(-1)=y'(1),
%
% Call: [lam,phi,x]=eig42F(64,[-1 1],2,60);
% Calculates numeigval eigenvalues lam and eigenfunctions phi at the grid x
% by the Galerkin method
%
tic;phi=zeros(n,numeigval);fxpfex=zeros(1000,5);
x=pd(n,dom,kind);X=mult(n,dom);D=deriv(n,dom);E=speye(n);
A=(E-X^2)^2*D^4*(E-X^2)^2;B=(E-X^2)^2*D^2*(E-X^2)^2;
A=A(1:n-4,1:n-4);B=B(1:n-4,1:n-4);
[V,L]=eig(full(A),full(B));
[LL,ind]=sort(diag(L),'descend');VV=V(:,ind);
lam=LL(1:numeigval);toc;
for s=1:numeigval
phi(:,s)=t2x((eye(n)-X^2)^2*[VV(:,s);0;0;0;0],kind);
end
disp('Numerical eigenvalues');disp(lam(1:5));
xx=linspace(dom(1),dom(2),1000);
for k=1:5
if k==1, fex=1+cos(pi*xx);
elseif k==2,
fex=sin(sqrt(-lam(2))*xx)-xx*sqrt(-lam(2))*cos(sqrt(-lam(2)));
elseif k==3, fex=1-cos(2*pi*xx);
elseif k==4,
fex=sin(sqrt(-lam(4))*xx)-xx*sqrt(-lam(4))*cos(sqrt(-lam(4)));
else fex=1+cos(3*pi*xx);
end
figure(k);
fx = barycheb(xx,phi(:,k),x,kind);fxpfex(:,k)=fx./fex;
plot(xx,fx,'b',xx,fxpfex(:,k),'r');
grid;legend('numeric','numeric/exact');
title(['\lambda_{',num2str(k),'} = ',num2str(lam(k))]);
xlabel('x');ylabel(['\phi_{',num2str(k),'}(x)']);
end
disp('test eigenvalues');
disp([-lam/pi^2,sin(sqrt(-lam))-sqrt(-lam).*cos(sqrt(-lam))]);
end