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_ode2_test(n,dom,kind,numeigval)
function [lam,phi,x]=eig_ode2_test(n,dom,kind,numeigval)
% Example: -y''+100cos^2(x)y=lambda y, x in [0,pi]
%          y(0)=y(pi), y'(0)=y'(pi)
% From: M. Plum, ZAMP 41 (1990), 205-226
% call: [lam,phi,x]=eig_ode2_test(128,[0 pi],2,7);
% exact eigenvalues:
% 9.7432204_0^6
% 28.6851_38^41
% 46.47783_4^7
% 62.9864_86^93
% 77.8052_37^44
% 91.80_07^11
% 98.975_4^8
% call: [lam,phi,x]=eig_ode2_test(64,[0 pi],2,7);
%
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;myBC;[V,L]=eig(full(A),full(B));
[LL,ind]=sort(diag(L));VV=V(:,ind);
lam=LL(1:numeigval);
for s=1:numeigval
phi(:,s)=t2x(VV(:,s),kind);
end
toc;myOUT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDE
% describes the differential operator
A=-D^2+100*funm(full(X),@cos)^2;B=speye(n);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myBC
% implements the homogeneous boundary conditions
T=cpv(n,dom,dom);
A(n-1,:)=T(1,:)-T(2,:);B(n-1,:)=A(n-1,:)/1.e8;
A(n,:)=T(1,:)*D-T(2,:)*D;B(n,:)=A(n,:)/1.e8;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
display('Numerical eigenvalues         Plum enclosure');
plum=[9.74322040,9.74322046;28.685138,28.685141;46.477834,46.477837;...
62.986486,62.986493;77.805237,77.805244;91.8007,91.8011;98.9754,98.9758];
for k=1:7
fprintf('  %12.8f,       [%11.8f, %12.8f]\n',lam(k),plum(k,1),plum(k,2));
end
xx=linspace(dom(1),dom(2),1000);
for k=1:numeigval
fx = barycheb(xx,phi(:,k),x,kind);
figure(k);plot(xx,fx);grid;
title(['\lambda_{',num2str(k),'} = ',num2str(lam(k))]);
xlabel('x');ylabel(['\phi_{',num2str(k),'}(x)']);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end