sp.m

LiScNLS is a Matlab application for the numerical study of some nonlinear differential equations o

function [lam,phi,phip,x,wint,wipr,dom]=sp(fun,n,errtol)
%
% Calculates the eigenvalues lam, the eigenfunctions phi and their derivatives phip of the problem
% -p*y"-dp*y'+q*y=lambda*r*y, x in dom=(a,b)
% alp(11)*y(a)+alp(12)*y'(a)=0
% alp(21)*y(b)+alp(22)*y'(b)=0
% using the Chebyshev-tau method 
% fun gives p,dp,q,r,alp,dom
% 
e=0:n;
D=2*spdiags(repmat(e',1,floor(n/2)),1:2:n-1,n,n);
tm=(-1).^(0:n-1);tp=ones(1,n);tm(1)=0.5;tp(1)=0.5;
time=cputime;
[P,DP,Q,R,alp,dom,x,wint,wipr]=feval(fun,n);% the problem
if isa(alp,'char') 
    lam=P;phi=DP;phip=Q;
    return;
end
alpha=(dom(2)-dom(1))/2;
D=D/alpha;
% the discretization
A=-P*D^2-DP*D+Q*eye(n);B=R*eye(n);
AA(3:n,:)=A(1:n-2,:);BB(3:n,:)=B(1:n-2,:);
AA(1,:)=tm*(alp(1,1)*eye(n)+alp(1,2)*D);
AA(2,:)=tp*(alp(2,1)*eye(n)+alp(2,2)*D);
BB(1,:)=1/5.e10*AA(1,:);
BB(2,:)=1/6.e10*AA(2,:);
[V,L]=eig(AA,BB);
DD=diag(L);
% filter the "good" eigenvalues (real and ordered)
ind1=find(abs(imag(DD))<errtol);
eigenval=DD(ind1);eigenvec=V(:,ind1);
eigenval=real(eigenval);eigenvec=real(eigenvec);% eigenvectors in the spectral space
[eigl,ind]=sort(eigenval);
% 
t=x2t(n,x,dom);
eigv=t'*eigenvec(:,ind); % eigenvectors in the physical space
flag=0;
% test the orthogonality of the filtered eigenvectors
for i=1:length(ind)
    I=wip(eigv(:,i),eigv(:,1:i),wipr,wint,dom);
    eigv(:,i)=eigv(:,i)/sqrt(I(i));I=I/I(i);
    I(i)=0;ind2=find(abs(I)>errtol);
    if length(ind2)>0
        % testing if  multiple (or clusters of) eigenvalues exist
        if max(abs(diff(eigl([ind2;i]))))>10*errtol break;end
        % Gram-Schmidt orthogonalization of corresponding eigenvectors 
        indm=fliplr([ind2;i]');flag=1;
        eigvm=eigv(:,indm);
        psi=zeros(size(eigvm));
        for ii=1:length(indm)
            v=eigvm(:,ii);
            for jj=1:ii
                v=v-wip(v,psi(:,jj),wipr,wint,dom)*psi(:,jj);
            end
            psi(:,ii)=v/sqrt(wip(v,v,wipr,wint,dom));
        end
        eigv(:,indm)=psi;
    end
end;
e=cputime-time;
if length(ind)==0 i=1;end;
disp([num2str(i-1),' eigenvalues computed by sp in ',num2str(e),' sec']);
lam=eigl(1:i-1);phi=eigv(:,1:i-1);
if flag==1
    phicof=inv(t')*phi;
    test=max(max(abs(AA*phicof-BB*phicof*diag(lam))));
    disp('Multiple (cluster) eigenvalues detected');
    disp(['Eigenfunctions test:   err=',num2str(test)]);
end
phip=t'*(D*(t'\phi));
for k=1:i-1
    if phi(1,k)<0 phi(:,k)=-phi(:,k);phip(:,k)=-phip(:,k);end
end