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.

[x,t,solnum]=pde_nonlin_ex4(n,dom,kind,dt,K)
function [x,t,solnum]=pde_nonlin_ex4(n,dom,kind,dt,K)
% Example: u_t=-u_xxx-3(u^2)_x, u(-20,t)=u(20,t)=u_x(20,t)=0, 
%          u(x,0)=2sech(x)^2 (Korteweg-de Vries equation)
% Exact solution:  solex=2*sech(x-4t)^2,
% Test:  I1=int_-20^20 u(x,t) dx = const,
%        I2=int_-20^20 u^2(x,t) dx = const,
%        I3=int_-20^20 (u_x^2(x,t)-u^3(x,t)) dx = const.
% Remark: uses matricial method/fixed point iterations
% call: [x,t,solnum]=pde_nonlin_ex4(128,[-20,20],2,0.01,400);
% 
[x,w]=pd(n,dom,kind);Tbc=cpv(n,dom,dom);D=deriv(n,dom);t(1)=0;
solnum(:,1)=2*sech(x).^2;
sol(:,1)=x2t(solnum(:,1),kind);
L=-D^3;T(1,:)=Tbc(1,:);T(2,:)=Tbc(2,:);T(3,:)=Tbc(2,:)*D;
Lb=L(1:n-3,1:n-3);Lbb=L(1:n-3,n-2:n);
Tb=T(:,1:n-3);Tbb=T(:,n-2:n);
Ltt=Lbb/Tbb;Lt=Lb-Ltt*Tb;
E=expm(Lt*dt/2);Phi=(E-eye(n-3))/(Lt*dt/2);
Pu=[eye(n-3);-Tbb\Tb];Pd=zeros(n-3,n);Pd(1:n-3,1:n-3)=eye(n-3);
figure(1);
p=plot(x,solnum(:,1),'EraseMode','none');axis([-20.1 20.1 -0.1 2.1]);grid;
title(['t= ',num2str(0)]);xlabel('x');ylabel('u(x,t)');hold on;pause(0.3)
for k=2:K+1
     t(k)=(k-1)*dt;uo=Pd*sol(:,k-1);
     for j=1:50
     un=E*Pd*sol(:,k-1)-dt/2*Phi*Pd*(3*D*x2t(t2x(Pu*uo,kind).^2,kind));
     if norm(un-uo)<1.e-7, break;end
     uo=un;
     end
     sol(:,k)=Pu*(E*uo-dt/2*Phi*Pd*(3*D*x2t(t2x(Pu*uo,kind).^2,kind)));
     solnum(:,k)=t2x(sol(:,k),kind);
     set(p,'color','w');set(p,'Ydata',solnum(:,k),'color','b');
     title(['t= ',num2str(t(k))]);hold on;pause(0.1)
end
figure(2);I1=w*solnum;I2=w*solnum.^2;I3=w*(t2x(D*sol,kind).^2-solnum.^3);
plot(t,I1,'r',t,I2,'b',t,I3,'g');legend('I1','I2','I3');xlabel('t');ylabel('The invariants');
figure(3); mesh(x,t,solnum');xlabel x;ylabel t;zlabel u(x,t);title('The solution');
for j=1:K+1,solex(:,j)=2*sech(x-4*t(j)).^2;nor(j)=max(abs(solex(:,j)-solnum(:,j)));end;
figure(4);plot(1:K+1,nor);xlabel('t');ylabel('err');title('The errors');
end 

Contact us