Diffusion.m

Presentation
Boundary Conditions
Initial Condition

Presentation

Here the code for the diffusion equation. We took the Vibrating String code and modified it to show the diffusion phenomenon. The equation used is the well-known diffusion equation also the Heat Equation without source : \displaystyle f_t=\mu f_{xx} The resolution is the same that the vibrating string, a march in time for the time derivative and a differentiation matrix for the spatial derivative.

clear all; close all;

% parameters
N=200; % number of gridpoints
L=50; % domain length
mu=0.5; % thermic conductivity
dt=0.05; % time step
tmax=20; % final time
x0=L/2;
t0=0.5;

% differentiation matrices
scale=-2/L;
[x,DM] = chebdif(N,2); 
dx=DM(:,:,1)*scale;    
dxx=DM(:,:,2)*scale^2;    
x=(x-1)/scale; 
I=eye(N); 

% system matrices
E=I;
A=mu*dxx;

Boundary Conditions

We know that at the boundary conditions, the solution is zero, thus we put this value in the differentiation matrix like this :

% boundary conditions
E([1 N],:)=0; 
A([1 N],:)=I([1 N],:);

% march in time matrix 
M=(E-A*dt/2)\(E+A*dt/2);

Initial Condition

For the initial condition we use the Gaussian function that represent well the diffusion at the initial time : \displaystyle T=\frac{1}{\sqrt{2\pi t_0}}exp\left(\frac{-(x-x_0)^2}{2t_0}\right)

% initial condition
q=1/sqrt(2*pi*t0)*exp(-(x-x0).^2/2/t0);

% marching loop
tvec=0:dt:tmax; 
Nt=length(tvec);
filename = 'Diffusion.gif';
for ind=1:Nt     
    q=M*q; % one step forward
    t=ind*dt;
    
    Tt=1/sqrt(2*pi*(t+t0))*exp(-(x-x0).^2/2/(t+t0));
   
    % plotting
    subplot(1,2,1);    
    plot(x,q,'k+',x,Tt,'g-'); 
    title('Diffusion');
    xlabel('x'); ylabel('T');
    axis([L/4,3*L/4,-0.1,0.5])
    
    subplot(1,2,2);
    plot(t,norm((q-Tt),2),'b.'); grid on; hold on
    drawnow
    title('Error')
    xlabel('t'); ylabel('T');
    
end