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,I1,I2,I3,time]=KdV(n,dom,kind,dt,K)
function [x,t,solnum,I1,I2,I3,time]=KdV(n,dom,kind,dt,K)
%
% From: K. Brauer, The Korteweg-de Vries Equation: History, 
% exact Solutions, and graphycal Representation, 
% http://www.usf.uni-osnabrueck.de/~kbrauer/solitons.html
%
% u_t=-u_xxx-3(u^2)_x, u(-50,t)=u(50,t)=u_x(50,t)=0, t in [-20,20]
% u(x,-20)=s1+s2+s3, si=bi*sech(sqrt(bi/2)*(x-2bit))^2, 
% b1=0.4,b2=0.7, b3=1
%
% Test:  I1=int_-20^20 u(x,t) dx = const,
%        I2=int_-20^20 u^2(x,t) dx = const,
%        I3=int_-50^50 (u_x^2(x,t)-2*u^3(x,t)) dx = const.
% Call: [x,t,solnum,I1,I2,I3,time]=KdV(256,[-50,50],2,0.05,800);
% 
[x,w]=pd(n,dom,kind);Tbc=cpv(n,dom,dom);D=deriv(n,dom);t(1)=-20;
solnum(:,1)=0.4*sech(sqrt(0.4/2)*(x-2*0.4*t)).^2+...
    0.7*sech(sqrt(0.7/2)*(x-2*0.7*t)).^2+1*sech(sqrt(1/2)*(x-2*1*t)).^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([dom(1)-0.1 dom(2)+0.1 -0.1 2.1]);grid;
title(['t= ',num2str(0)]);xlabel('x');ylabel('u(x,t)');hold on;pause(0.01)
for k=2:K+1
    time(k)=cputime;
    t(k)=-20+(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);
    time(k)=cputime-time(k);
    set(p,'color','w');set(p,'Ydata',solnum(:,k),'color','b');
    title(['t = ',num2str(t(k),'%3.1f'),'   iter = ',num2str(j,'%2.0f')]);
    hold on;pause(0.01)
end
figure(2);I1=w*solnum;I2=w*solnum.^2;I3=w*(t2x(D*sol,kind).^2-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');
figure(4); waterfall(x,t(1:25:750),solnum(:,1:25:750)');
xlabel x;ylabel t;zlabel u(x,t);title('The solution');
end 

Contact us