Code covered by the BSD License  

Highlights from
Chebpack

image thumbnail

Chebpack

by

 

15 Jul 2011 (Updated )

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

[LL,evec, x]=Greenberg6(n,dom,kind,numeigval)
function [LL,evec, x]=Greenberg6(n,dom,kind,numeigval)
% From: Leon Greenberg, Marco Marletta, Algorithm 775: the code SLEUTH for 
% solving fourth-order Sturm-Liouville problems, ACM Transactions on 
% Mathematical Software, 23, 4(1997), 453-493, Example 6 
%
% (Py'')''-(Sy')'=lambda y, P=(1-x^2)^2, S=12-4x^2
% -8y'(-1)=lambda y(-1), 8y'(1)=lambda y(1)
% -(Py'')'+Sy'=lambda/x y-lambda/x/8Py'', x=-1,1
%
% Call [lam,phi,x]=Greenberg6(96,[-1,1],2,10)
D=deriv(n,dom);X=mult(n,dom);E=speye(n);P=(E-X^2)^2;S=12*E-4*X^2;
[x,w]=pd(n,dom,kind);
A=D^2*(P*D^2)-D*(S*D);B=E;
T=cpv(n,dom,dom);
A(n-3,:)=-8*T(1,:)*D;B(n-3,:)=T(1,:);
A(n-2,:)=8*T(2,:)*D;B(n-2,:)=T(2,:);
A(n-1,:)=T(1,:)*(-D*(P*D^2)+S*D);B(n-1,:)=T(1,:)*(-E+P*D^2/8);
A(n,:)=T(2,:)*(-D*(P*D^2)+S*D);B(n,:)=T(2,:)*(E-P*D^2/8);
[V,L]=eig(full(A),full(B));LD=diag(L);
[LL,ind]=sort(LD,'ascend');VV=V(:,ind);evec=t2x(VV,kind);
k=(0:numeigval-1)';
lex=8*k+16*k.*(k-1)+8*k.*(k-1).*(k-2)+k.*(k-1).*(k-2).*(k-3);
display([LL(1:10),lex]);
for k=1:5,figure(k);plot(x,evec(:,k));title(num2str(k-1));pause(0.01);end;

Contact us