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.

pde_nonlin_ex3.m
function [x,t,solnum]=pde_nonlin_ex3(n,dom, kind,dt,K)
% Example: u_t=i/2 u_xx +iu|u|^2, x in [-4,4],t>0, 
%          u(x,0)=4exp(2ix)sech(4x), u(-4,t)=0, u(4,t)=0
% From: T. A. Driscoll et.al., BIT Numerical Mathematics,(2008)48: 701723
% Remark: Cubic Schroedinger equation with reflection
%         uses  2-backward difference / 3-Adams-Bashforth        
% call: [x,t,solnum]=pde_nonlin_ex3(128,[-4,4],2,0.01,400);
%
x=pd(n,dom,kind);D=deriv(n,dom);T=cpv(n,dom,dom);
A=speye(n)-1i*dt/2*(D^2);A(end-1,:)=T(1,:);A(end,:)=T(2,:);
solnum(:,1)=4*exp(2i*x).*sech(4*x);t(1)=0;
sol(:,1)=x2t(solnum(:,1),kind);% initial condition
sol(n-1,1)=0;sol(n,1)=0;% boundary condition is satisfied by uo
figure(1)
p=plot(x,abs(solnum(:,1)),'EraseMode','none');
axis([-4.2 4.2 -0.2 5.2]);title('t= 0');grid;xlabel('x');ylabel('u(x,t)');
for k=2:3 % starting iterations
    ux=solnum(:,k-1);uo=sol(:,k-1);
    b=uo+dt*1i*x2t(ux.*(ux.*conj(ux)),kind);
    b(n-1)=0;b(n)=0; %boundary condition is satisfied by b
    un=A\b;t(k)=(k-1)*dt;
    solnum(:,k)=t2x(un,kind);sol(:,k)=un;
    set(p,'color','w');
    set(p,'Ydata',abs(solnum(:,k)),'color','b');
    title(['t= ',num2str(t(k))]);drawnow;
end
A=1.5*speye(n)-1i*dt/2*(D^2);A(end-1,:)=T(1,:);A(end,:)=T(2,:);
for k=4:K % 2-backward difference / 3-Adams-Bashforth
    uxc=solnum(:,k-1);uc=sol(:,k-1);Gc=1i*x2t(uxc.*(uxc.*conj(uxc)),kind);
    uxo=solnum(:,k-2);uo=sol(:,k-2);Go=1i*x2t(uxo.*(uxo.*conj(uxo)),kind);
    uxoo=solnum(:,k-3);uoo=sol(:,k-3);Goo=1i*x2t(uxoo.*(uxoo.*conj(uxoo)),kind);
    b=2*uc-0.5*uo+dt*(8*Gc-7*Go+2*Goo)/3;
    b(n-1)=0;b(n)=0; %boundary condition is satisfied by b
    un=A\b;t(k)=(k-1)*dt;
    solnum(:,k)=t2x(un,kind);sol(:,k)=un;
    set(p,'color','w');
    set(p,'Ydata',abs(solnum(:,k)),'color','b');
    title(['t= ',num2str(t(k))]);drawnow;
end
figure(2);mesh(x,t,abs(solnum)');xlabel x;ylabel t;zlabel u(x,t);view([20,50]);title('The solution');
end

Contact us