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,solnum]=pde_lin_matr_ex1(n,dom,kind,dt,K)
function [x,solnum]=pde_lin_matr_ex1(n,dom,kind,dt,K)
% Example u_t=u_xx+20u_x, x in [a,b],t>0
%         u(a,t)=alpha=1, u(b,t)=beta=2 ,u(x,0)=uo, 
% Steady solution: u(x)=2-exp(-20(x+4))
% Remark: uses matricial exponential method in spectral form
% call: [x,solnum]=pde_lin_matr_ex1(64,[-4,2],2,0.01,100);
% 
x=pd(n,dom,kind);D=deriv(n,dom);T=cpv(n,dom,dom);
alpha=[];beta=[];t=zeros(K+1);
myINIT;myDE;myBC;solnum(:,1)=uo;sol(:,1)=x2t(solnum(:,1),kind);
for k=1:K
    t(k+1)=k*dt;AE=expm(t(k+1)*At);BE=(AE-eye(n-2))/At;
    ub=AE*sol(1:n-2,1)+BE*Att*[alpha;beta];
    sol(:,k+1)=[ub;Tbb\([alpha;beta]-Tb*ub)];
    solnum(:,k+1)=t2x(sol(:,k+1),kind);
end
myOUT;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myINIT
    % describes the initial condition in physical form
    uo=max(0,1-abs(x));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDE
    % describes the linear spatial part of the pde
    A=D^2+20*D;
    Ab=A(1:n-2,1:n-2);Abb=A(1:n-2,n-1:n);
    Tb=T(:,1:n-2);Tbb=T(:,n-1:n);
    Att=Abb/Tbb;At=Ab-Att*Tb;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myBC
    % describes the boundary conditions
    alpha=1;beta=2;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
    % describes the output of the program
    figure(1);
    p=plot(x,uo,'EraseMode','none');axis([-4.1 2.1 -0.1 2.1]);grid;
    title(['t= ',num2str(0)]);xlabel('x');ylabel('u(x,t)');hold on;pause(1)
    for j=2:K+1
        set(p,'color','w');set(p,'Ydata',solnum(:,j),'color','b');
        title(['t= ',num2str(t(j))]);hold on;pause(0.1)
    end
    figure(2);semilogy(x,abs(solnum(:,end)-(2-exp(-20*(x+4)))));
    xlabel('x');ylabel('err');title('The error at the last time step');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

Contact us