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.

KS_evol.m
function [lam,phi,phip,t,C,kodmax]=KS_evol(n,numeigval,dom,m,dt,K)
% 
% From: B. FORNBERG, A PSEUDOSPECTRAL FICTITIOUS POINT METHOD FOR HIGH
% ORDER INITIAL-BOUNDARY VALUE PROBLEMS, SIAM J. SCI. COMPUT., 
% Vol. 28, No. 5, pp. 17161729
%
% u_t+u_xxxx+u_xx+uu_x=0, [-16*pi,16*pi]
% u(x,0)=0.1sin^2(x/2) in (-12pi,-10pi) and 0 otherwise
% u(-16*pi,t)=u(16*pi,t)=u_x(-16*pi,t)=u_x(16*pi,t)=0
%
% Call: [lam,phi,phip,t,C]=KS_evol(256,100,[-16*pi,16*pi],0,0.1,1500);
% Kuramoto-Sivashinski problem, Lyapunov-Schmidt method
% 
disp('Linear part');tic;
phi=zeros(n,numeigval);phip=zeros(n,numeigval);
E=speye(n);D=deriv(n,dom);X=mult(n,dom);F=(dom(2)^2*E-X^2)^2;
[xleg,wleg]=pd(n,dom,3);rleg=ones(n,1);Tleg=cpv(n,xleg,dom);
A=F*(D^4+D^2)*F;B=F*F;AA=A(1:n-4,1:n-4);BB=B(1:n-4,1:n-4);
[V,L]=eig(full(AA),full(BB));
[LL,ind]=sort(diag(L),'ascend');VV=V(:,ind);
lam=LL(1:numeigval);
for s=1:numeigval
    phi(:,s)=Tleg*F*[VV(:,s);0;0;0;0];nphi=sqrt(wleg*phi(:,s).^2);
    phi(:,s)=phi(:,s)/nphi;
    phip(:,s)=Tleg*D*F*[VV(:,s);0;0;0;0]/nphi;
end;
disp(['Calculated ',num2str(numeigval),' eigenvalues and eigenfunctions']);
I=wip(phi,phi,rleg,wleg);orthogonal=max(max(abs(I-eye(numeigval))));
disp(['test orthogonality: ',num2str(orthogonal)]);toc;
disp('Evolution');tic;
U=ICK(xleg);
u0=U;Uo=wip(phi,u0,rleg,wleg);C=zeros(numeigval,K+1);C(:,1)=Uo;
Lam=1+dt/2*lam;
p=plot(xleg,u0,'EraseMode','none');axis([dom(1) dom(2) -5 5]);
grid;xlabel('x');ylabel('u(x,t)');
for k=1:K % time steps
    c0=Uo;uo=phi*Uo;uop=phip*Uo;
    [Nuo,~,~]=feval(@NK,xleg,uo,uop);
    CONST=uo-dt/2*phi*(Uo.*lam)+dt/2*Nuo;kodmax=1;
    [c,~,kod]=lisc(@Nwork,Lam,phi,phip,xleg,rleg,wleg,1.e-6,m,c0);
    set(p,'color','w');set(p,'Ydata',phi*c,'color','b');
    title(['t = ',num2str(k*dt,'%3.1f')]);drawnow;grid;pause(0.01);
    C(:,k+1)=c;if kod>kodmax, kodmax=kod;end
    Uo=c;
end;
% Output
solv=phi*C;t=(1:K+1)'*dt;toc;
% Output
surf((xleg),t,(solv)'), shading interp, lighting phong, axis tight
view([30 70]), colormap(autumn); set(gca,'zlim',[-5 5])
light('color',[1 1 0],'position',[-1,2,2])
material([0.30 0.60 0.60 40.00 1.00]);
xlabel('x');ylabel('t');zlabel('u(x,t)');
function U=ICK(x)
U=zeros(size(x));
ind=find(-12*pi<x & x<-10*pi);
U(ind)=0.1*sin(x(ind)/2).^2;
end
function [Nu,Npu,Npup]=NK(x,u,up)
Nu=-u.*up;
Npu=-up;
Npup=-u;
end
function [Nu,Npu,Npup]=Nwork(xleg,u,up)
    [Nu,Npu,Npup]=feval(@NK,xleg,u,up); 
    Nu=dt/2*Nu+CONST;
    Npu=dt/2*Npu;
    Npup=dt/2*Npup;
end
end

Contact us