Code covered by the BSD License

# Chebpack

### Damian Trif (view profile)

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,
% 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
```