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.

[lam,phi,x]=eig_ode4_2_test(n,dom,kind,numeigval)
function [lam,phi,x]=eig_ode4_2_test(n,dom,kind,numeigval)
% Example: y^(4)=lambda y'', x in [-1,1]
%          y(-1)=y'(-1)=y(1)=y'(1)=0
% From: M. Charalambides, F.Waleffe, SIAM J. Numer. Anal. 47, pp. 48-68
% Exact eigenvalues: lam=-n^2 pi^2 or -qn^2, qn=tan(qn)
% call:  [lam,phi,x]=eig_ode4_2_test(80,[-1 1],2,32);
%
tic;A=zeros(n);B=zeros(n);phi=zeros(n,numeigval);
x=pd(n,dom,kind);X=mult(n,dom);D=deriv(n,dom);
myDE;
[V,L]=eig(full(A),full(B));
[LL,ind]=sort(abs(diag(L)));VV=V(:,ind);
lam=-LL(1:numeigval);
for s=1:numeigval
phi(:,s)=t2x((eye(n)-X^2)^2*[VV(:,s);0;0;0;0],kind);
end
toc;
myOUT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDE
% describes the differential operator
E=eye(n);A=D^4*(E-X^2)^2;B=D^2*(E-X^2)^2;
A=A(1:n-4,1:n-4);B=B(1:n-4,1:n-4);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
I=lam<0;slam=sqrt(-lam(I));
slam1=slam;slam1(2:2:end)=0;slam2=slam;slam2(1:2:end)=0;
display(['lambda                ','sqrt(-lambda)/pi=n','   sqrt(-lambda)-tan(sqrt(-lambda))=0']);
display([lam(I),slam1/pi,slam2.*cos(slam2)-sin(slam2)]);
xx=linspace(dom(1),dom(2),1000);
for k=1:5
figure(k);
fx = barycheb(xx,phi(:,k),x,kind);
plot(xx,fx);grid;
title(['\lambda_{',num2str(k),'} = ',num2str(lam(k))]);
xlabel('x');ylabel(['\phi_{',num2str(k),'}(x)']);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end