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_ode.m
function [lam,phi,phip,t,C]=KS_ode(n,numeigval,dom,dt,tfin)
% 
% 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_ode(256,100,[-16*pi,16*pi],0.1,150);
% Kuramoto-Sivashinski problem, direct 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;
u0=ICK(xleg);c0=wip(phi,u0,rleg,wleg);
[t,C]=ode15s(@Kfun,0:dt:tfin,c0);C=C';
solv=phi*C;toc;
% Output
figure(1);p=plot(xleg,u0,'EraseMode','none');axis([-50 50 -5 5]);
grid;xlabel('x');ylabel('u(x,t)');
for k=1:length(t)
set(p,'color','w');set(p,'Ydata',solv(:,k),'color','b');
drawnow;grid;title(['t = ',num2str(t(k),'%3.1f')]);pause(0.01)
end
figure(2)
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 dcdt=Kfun(~,c)
u=phi*c;up=phip*c;
Nu=-u.*up;
coef=wip(phi,Nu,rleg,wleg);
dcdt=-diag(lam)*c+coef;
end
end

Contact us