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]=OrrSom(n,dom,kind,numeigval)
function [lam,phi,x]=OrrSom(n,dom,kind,numeigval)
% From: J. J. Dongarra, B. Straughan, D. W. Walker, Chebyshev tau -
% QZ Algorithm Methods for Calculating Spectra of Hydrodynamic Stability
% Problems, Applied Numerical Mathematics archive, 22, 4(1996), 399 - 434,
% Also from: McFadden et. al., Elimination of spurious egenvalues in the
% Chebyshev tau spectral method, J. COMP. PHYS., 91, 228-239 (1990)
%
% (D^2-a^2)^2/(-iaRe)y+U(D^2-a^2)y-U''y=\lambda(D^2-a^2)y, x in [-1,1]
% y(-1)=y(1)=0=y'(-1)=y'(1), a=1,Re=5772 or Re=1.e4,U=1-x^2
%
% Call: [lam,phi,x]=OrrSom(128,[-1 1],2,60);
%
tic;A=zeros(n);B=zeros(n);phi=zeros(n,numeigval);a=1;Re=5772;%Re=1.e4;
x=pd(n,dom,kind);X=mult(n,dom);D=deriv(n,dom);E=speye(n);
A=(E-X^2)^2*((D^2-a^2*E)^2/(-1i*Re*a)+(E-X^2)*(D^2-a^2*E)+2*E)*(E-X^2)^2;
B=(E-X^2)^2*(D^2-a^2*E)*(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));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;
display(lam);
plot(real(-1i*lam),imag(-1i*lam),'*k');grid;
xlabel('real(\lambda)');ylabel('imag(\lambda)');
end