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_test(n,dom, kind,dt,K)
function [x,t,solnum]=pde_nonlin_test(n,dom, kind,dt,K)
% Example: u_t=niu u_xx - [(u^2)/2]_x, x in [-1,1],t>0, niu=0.01/pi
%          u(x,0)=-sin(pi x), u(-1,t)=u(1,t)=0 (Burgers equation)
% From: C. Basdevant et. al., Comput. Fluids, 14, (1986), pp. 23-41
% Test:  max |u_x| = 152.00516 at the instant t_max=1.6037/pi
% call: [x,t,solnum]=pde_nonlin_test(513,[-1,1],2,5.e-4,1200);
% 
x=pd(n,dom,kind);D=deriv(n,dom);T=cpv(n,dom,dom);niu=0.01/pi;
myDElin;A=speye(n)-dt*A;myINIT;uo=x2t(f,kind);solnum(:,1)=f;sol(:,1)=uo;
h = waitbar(0,'Please wait...');set(h,'name','Burgers equation');t=zeros(K+1);
for k=2:K+1
    waitbar((k-1)/K,h);
    myDEnonlin;myBC;un=A\b;t(k)=(k-1)*dt;
    uo=un;solnum(:,k)=t2x(uo,kind);sol(:,k)=uo;
end
close(h);myOUT;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myINIT
    % describes the initial condition in physical form
    f=-sin(pi*x);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDElin
    % describes the linear spatial part of the pde
    A=niu*(D^2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDEnonlin
    % describes the linear spatial part of the pde
    b=uo-dt*D*x2t((t2x(uo,kind).^2)/2,kind);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myBC
    % describes the boundary conditions
    A(end-1,:)=T(1,:);A(end,:)=T(2,:);
    b(n-1)=0;b(n)=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
    p=plot(x,f,'EraseMode','none');xlabel('x');ylabel('u(x,t)');
    axis([-1.2 1.2 -1.2 1.2]);title('t= 0');grid;hold on;pause(1);
    for j=1:K
        if rem(j,20)==0,set(p,'color','w');
            set(p,'Ydata',solnum(:,j+1),'color','b');
            title(['t= ',num2str(t(j+1))]);hold on;pause(0.1)
        end
    end
    m=abs(cpv(n,0,dom)*(D*sol));[M,ind]=max(m);
    display(['maxdiff= ',num2str(M),', tmax= ',num2str(ind*dt)]);
    display(['maxdiff_analytical= ',num2str(152.00516),...
        ',tmax_analytical= ',num2str(1.6037/pi)]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

Contact us