Diffusion of the dead - The maths of zombie invasions. Part 4, Simulating zombie movement.

Last time we presented the diffusion equation
\begin{equation}
\frac{\partial Z}{\partial t}(x,t)=D\frac{\partial^2 Z}{\partial x^2}(x,t).
\end{equation}
and demonstrated that it had the right properties to model zombie motion. However, stating the equation is not enough. We must add additional information to the system before we can solve the problem uniquely. Specifically, we need to define: the initial state of the system, where the boundaries of the system are, and, finally, what happens to the zombies at the boundaries.

For the boundary condition we assume that the zombies cannot move out of the region $0\leq x\leq L$. This creates theoretical boundaries which the population cannot cross; the zombies will simply bounce off these boundaries and be reflected back into the domain. Since no zombies can cross the boundaries at $x=0$ and $x=L$, the `flux of zombies' across these points is zero,
\begin{equation}
\frac{\partial Z}{\partial x}(0,t)=0=\frac{\partial Z}{\partial x}(L,t) \text{ (the zero flux boundary conditions).}
\end{equation}

For the initial condition, we assume that the zombies are all localised in one place, a graveyard for example. Thus the zombies have a density of $Z_0$ zombies/metre in the region $0\leq x \leq 1$,
\begin{equation}
Z(x,0)=\left\{\begin{array}{cc}
Z_0&\text{for }0\leq x\leq 1,\\
0&\text{for } x>1,
\end{array}
\right. \text{ (the initial condition).}
\end{equation}

The diffusion with the given initial and boundary conditions can be solved exactly and has the form,
\begin{equation}
Z(x,t)=\frac{Z_0}{L}+\sum^\infty_{n=1}\frac{2Z_0}{n\pi}\sin\left(\frac{n\pi}{L}\right) \cos\left(\frac{n\pi}{L}x\right)\exp\left({-\left( \frac{n\pi}{L} \right)^2Dt}\right).\label{Solution}
\end{equation}
Although I have made this solution magically appear from nowhere, the solution can be rigorously produced using methods called "separation of variables" and "Fourier series".

Separating the variables means that we assume space and time components of the solution are not coupled together in a complicated manner. Namely, the solution can be written as a spatial component, multiplied by a time component, which can be seen above because the space variable, $x$, only appears in the cosine function, whilst the time variable, $t$, only appears in the exponential function. Fourier series allows us to write a function as an infinite summation of sines and cosines. Although this may seem cumbersome the Fourier series usually have very nice properties, which allow them to be manipulated with ease.

The `$\sin$' and `$\cos$' functions are those that the reader may remember from their trigonometry courses. The exponential function, `$\exp$', is one of the fundamental operators of mathematics, but for now the only property that we are going to make use of is that if $a>0$ then $\exp(-at)\rightarrow 0$ as $t\rightarrow \infty$. Using this fact we can see that as $t$ becomes large most of the right-hand side of the solution becomes very small, approximately zero. Hence, for large values of $t$, we can approximate
\begin{equation}
Z(x,t)\approx\frac{Z_0}{L}.\label{Longtime_approximation}
\end{equation}
This means that as time increases, zombies spread out evenly across the available space, with average density $Z_0/L$ everywhere.

The simulation illustrates two solutions to the diffusion equation for the zombie density. Namely, the video above compares the  analytical solution, given in equation \eqref{Solution}, with a direct numerical solution of the diffusion equation. If you stop the video right at the start, the initial condition can be seen and it shows that there is a large density of zombies between $0\leq x \leq 1$. The movie then illustrates how diffusion causes the initial peak to spread out filling the whole domain. After 500 time units the density of zombies has become uniform throughout the domain. The Matlab codes for plotting this solutions can be found below.

There are a couple of things to note in this simulation. Firstly, we say that the analytical solution is only a approximation because we can not simulate an infinite number of cosine terms. The movie shows only the first 1000 terms and, as you can see, the comparison between the two results is pretty good. Secondly we have not given units for the density, space or time, thus, are we using seconds and metres or minutes and miles? Well, the answer is that in some sense it doesn't matter, whilst in other cases it matters a great deal. Since we are only interested in visualising the dynamics of diffusion, we can keep the units arbitrary. However, if we had a specific application, with data, then we would have to ensure that all our units were consistent.

This finishes our look at solving the diffusion equation. Next time we actually use these solutions to provide us with the first of ours answers. These answers will then provide us with a strategy to deal with the eventual zombie apocalypse.

________________________________________________________________

________________________________________

Below you will find the Matlab code used to generate the above movie. If you have Matlab then you can simply copy and paste the text into a script and run it immediately.

function Diffusion_solution
clear
close all
clc

Z0=100; %Initial density.
L=10; %Length of domain.
D=0.1; %Diffusion coefficient.
dx=0.1; %Space step.
x=[0:dx:L]; %Spatial discretisation.
dt=1; %Time step.
final_time=500; %Final time.
times=0:dt:final_time; %Time discretisation.
iter=1; %Parameter initialisation.

%% Analytical solution.
Z=ones(1+final_time/dt,length(x))*Z0/L;
for t=0:dt:final_time
for n=1:1000
Z(iter,:)=Z(iter,:)+2*Z0/(n*pi)*sin(n*pi/L)*cos(x*n*pi/L)*...
exp(-(n*pi/L)^2*t*D);
end
iter=iter+1;
end

%% Direct numerical solution.
sol = pdepe(0,@(x,t,u,DuDx)Zombies(x,t,u,DuDx,D),@(x)ICs(x,Z0),@BCs,x,times);

%% Plotting
for i=1:final_time+1
plot(x,Z(i,:),'b','linewidth',3)
hold on
plot(x,sol(i,:),'g--','linewidth',3)
set(gca,'yTick',[0:20:100])
xlabel('Distance, x','fontsize',20);
ylabel('Density, Z','fontsize',20);
set(gca,'fontsize',20)
grid off
axis([0 10 0 100])
title(['Time=',num2str((i-1)*dt)])
legend('Approximate analytic solution','Approximate numerical solution','location','northoutside')
drawnow
hold off
end

function value = ICs(x,Z0)
% Setting the initial condition, which is a step function. Namely, the
% density of zombies is Z0 when x is less than 1 and 0 everywhere else.
if x < 1
value=Z0;
else
value=0;
end
function [c,b,s] = Zombies(~,~,~,DuDx,D)
% Matlab syntax for the diffusion equation.
c = 1;
b = D*DuDx(1);
s = 0;
function [pl,ql,pr,qr] = BCs(~,~,~,~,~)
% Matlab syntax for the zero flux boundary conditions.
pl = 0;
ql = 1;
pr = 0;
qr = 1;