Diffusion of the dead - The maths of zombie invasions. Part 5, Time of first interaction.

Previously, we saw how to simulate the diffusion equation, which we're using to model zombie motion. Critically, the diffusion equation allows us to predict the density of zombies at all places and for all time. 

There are many questions we could answer with this equation; however, the most pressing question to any survivors is, "How long do we have before the first zombie arrives?". The mathematical formulation of this question is, "For what time, $t_z$, does $Z(L, t_z) = 1$?". Unfortunately, this does not have a nice solution that can be evaluated easily. However, we can calculate $t_z$ numerically to any degree of accuracy we choose by using a simple solution searching technique known as the "Bisection Search Algorithm"

The first thing to notice about $Z(L,t)$ is that it is monotonically increasing in time, thus, if $t_1>t_2$ then  $Z(L,t_1)>Z(L,t_2)$. This can be seen in a number of ways. For example, by watching the diffusion simulation from last time we see that as time increases the zombie population on the boundary only ever increases. This makes intuitive sense, because diffusion is causing the zombie population to spread out, so that it becomes uniform everywhere.

Using this knowledge then a simple method of solving $Z(L, t_z) = 1$ would be to substitute in a value of $t$. If $Z(L, t) < 1$ then we double $t$ and consider $Z(L, 2t)$, which will be greater than $Z(L,t)$, because of the monotonic property we discussed above. We keep doubling $t$ until we reach a value such that $Z(L, 2^nt) > 1$. Defining $t_0 = 2^{n-1}t$ and $t_1 = 2^nt$, then we know that there exists some $t_z$ between $t_0$ and $t_1$ such that $Z(L, t_z) = 1$.

To gain better approximations to the value of $t_z$, we start halving this domain and only keep the half that contains the solution. However, how do we know which half contains the solution? We once again rely on the monotonic property and evaluate the function at the midpoint of the domain. Explicitly, if  $Z(L,(t_0+t_1)/2) < 1$ then the solution is in the right half. Alternatively, if $Z(L,(t_0+t_1)/2) > 1$ then the solution is in the left half.

An example illustrating this concept is shown in Figure 1. In the initial setup, $Z(L,(t_0+t_1)/2)>1$, thus, we redefine $t_1 = (t_0 + t_1)/2$ and repeat the process.

Figure 1. Bisection technique. After each iteration, the domain, shown by the dashed lines,
becomes half as long as it was before.

After each iteration, we halve the size of the interval $[t_0, t_1]$, making it smaller and smaller. Thus, by design, since we always have $t_z$ within $[t_0, t_1]$ and by repeating the halving process, we can estimate $t_z$ to any accuracy we like.

The benefit of this method is its simplicity and reliability; it will always work. However, the cost of this reliability comes at the price of speed. If the initial searching region is very big, it may take a large number of iterations before the process produces an answer to an accuracy with which we are happy. There are quicker methods, but these are usually more complex and sometimes they may fail to find a solution altogether.

If the zombies are closing in on you and you need to compute their interaction times more quickly, we direct you to consider Newton-Raphson techniques rather than the bisection method. Although if the zombies really are that close may I suggest you focus on fighting them off before you read any further?

No matter what solution method you use you should end up with a graph like Figure 2, which illustrates the time in minutes we have before we meet a zombie depending on how far away we are and how fast the zombies are diffusing.

Figure 2. Time in minutes until the first zombie arrives for various rates of diffusion and distances.

At this point we can consider two general strategies. Either we run away, or we try and slow the zombie down. Both of these approaches will indeed increase the time it takes for the zombies to get to you. However, Figure 2 clearly shows that we gain much more time by running away compared to the strategy of slowing the zombies down. Next time I will expand on this point and show how to estimate the zombie interaction time much quicker.

________________________________________________________________

________________________________________

My amazing wife on her zombie themed hen do.
Note that the theme had nothing to do with me.

You may be wondering how I chose the variable ranges for Figure 2. Well, the Oxford maths department is approximately 90m away from a graveyard so the distances simply came from a matter of self preservation. The diffusion speed on the other hand was generated with the aid of my loving wife. I got her to wander randomly around, staggering like a zombie. I timed her and measured how far she had moved and, thus, calculated that her zombie diffusion rate was approximately 115m$^2$/minute. Don't worry we did the experiment three times and took an average, which ensures that this value is accurate.

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;

Diffusion of the dead - The maths of zombie invasions. Part 3, Diffusive motion.

As discussed previously, we are going to model the zombie motion using the diffusion equation. In this post we introduce the gritty details. I've interpreted the mathematical symbols intuitively, so, if you stick with it, you should find yourself understanding more than you ever thought you could.

It is impossible to overstate the importance of the diffusion equation. Wherever the movement of a modelled species can be considered random and directionless, the diffusion equation will be found. This means that by understanding the diffusion equation we are able to describe a host of different systems such as heat conduction through solids, gases (e.g. smells) spreading out through a room, proteins moving round the body, molecule transportation in chemical reactions and rainwater seeping through soil, to name but a few of the great numbers of applications.

If you've never come across diffusion before, or want to know more about it's basic properties the video below is a very good primer, although feels very much like a "Look around you" episode.

The mathematical treatment of diffusion begins by defining the variables that we will need. Let the density of zombies at a point $x$ and at a time $t$ be $Z(x,t)$ then the density has to satisfy the diffusion equation,

\begin{equation}
\frac{\partial Z}{\partial t}(x,t)=D\frac{\partial^2 Z}{\partial x^2}(x,t).
\end{equation}
To some an equation can be scarier than any zombie, but fear not. I am going to break this equation down into bits so that you are able to see the reality behind the mathematics.

Notice that the equation is made up of two terms, the left-hand side and the right-hand side, which are defined to be equal. Explicitly, the left-hand side is known as the time derivative and it simply tell us how the zombie density is changing over time,
\begin{equation}
\frac{\partial Z}{\partial t}(x,t)=\text{rate of change of $Z$ over time at a point $x$}.
\end{equation}

Although the numerical value of this term is important, what is more important is if the term is positive or negative. Specifically, if $\partial Z/\partial t$ is positive then $Z$ is increasing at that point in time, and, vice-versa, if $\partial Z/\partial t$ is negative then $Z$ is decreasing. Thus, we use this term to tell us how the zombie population is changing over time.

The term on the right-hand side is known as the second spatial derivative and it is a little more complicated than the time derivative. Essentially it encapsulates the idea that the zombies move from areas of high density to areas of low density (i.e. they spread out). To aid your intuitive understanding of this term see Figure 1.

Figure 1. A typical initial zombie density graph. There are regions of high zombie activity, e.g. a graveyard, and there are regions of low zombie density, e.g. your local library.

In the figure, there are initially more zombies on the left of the space than the right. Just before the peak in density the arrow (which is the tangent to the curve known as the spatial derivative, or $\partial Z/\partial x$ at this point) is pointing upwards. This means that as $x$ increases, so does the zombie density, $Z$. At this point

\begin{equation}
\frac{\partial Z}{\partial x}=\textrm{rate of change of $Z$ as $x$ increases} > 0.
\end{equation}

Just after the peak the arrow is pointing down thus, at this point,

\begin{equation}
\frac{\partial Z}{\partial x}=\textrm{rate of change of $Z$ as $x$ increases} < 0.
\end{equation}

Thus, at the peak, the spatial derivative is decreasing, because it goes from positive to negative. This, in turn, means that the second derivative is negative at the peak, because a negative second derivative means the first derivative is decreasing. This is analogous to statements made above about the sign of the time derivative and the growth, or decay, of the zombie population.

In summary, our hand wavy argument tells us that at local maximum $\partial^2 Z/\partial x^2<0$. Using the equality of the diffusion equation, this means that at a local maximum the time derivative is negative and, thus, the density of zombies is decreasing. A similar argument shows that the population of zombies at a local minimum increases. In summary, we see that diffusion causes zombies to move from regions of high density to low density.

Finally, we mention the factor $D$, which is called the diffusion coefficient. $D$ is a positive constant that controls the rate of movement. Specifically, the larger $D$ is the faster the zombies spread out.

And with that you now understand one of the most important partial differential equations in all of mathematics. That wasn't too hard was it? Next time we discuss the solution of the diffusion equation including some simulations and Matlab code for you to try yourself.

Diffusion of the dead - The maths of zombie invasions. Part 2, Important questions you need to ask in a zombie outbreak.

We begin modelling a zombie population in the same way that a mathematician would approach the modelling of any subject. We, first, consider what questions we want to ask, as the questions will direct which techniques we use to solve the problem. Secondly, we consider what has been done before and what factors were missing in order to achieve the answers we desire. This set of blog posts will consider three questions:

1. How long will it take for the zombies to reach us?
2. Can we stop the infection?
3. Can we survive?

In order to answer these questions I, Ruth Baker, Eamonn Gaffney and Philip Maini focused on the motion of the zombies as their speed and directionality would have huge effects on these three questions. Explicitly, we used a mathematical description of diffusion as a way to model the zombies motion. This was discussed in a previous post, but I recap the main points here.

• The original zombie infection article by Robert Smith? did not include zombie, or human movement.
• Zombies are well known for is their slow, shuffling, random motion. The end of Dawn of the Dead (shown in the YouTube clip below) gives some great footage of zombies just going about their daily business.
• This random motion is perfectly captured through the mathematics of diffusion.

Of course, there is plenty of evidence to suggest that zombies are attracted to human beings, as they are the predator to our prey. However, as we will see, we are going to be over run on a time scale of minutes! Thus, although mathematicians can model directed motion, and chasing, these additional components complicate matters. Further, random motion leads to some nice simple scaling formulas that can be used to quickly calculate how long you approximately have left before you meet a zombie.

Another simplifying assumption that we make is that we can model the zombie (and human) populations as continuous quantities. Again, this is incorrect as zombies are discrete units (even if they are missing body parts). Since we are making an assumption we will create an error in our solution. But how big is this error? In particular, if the error in the assumption is smaller than the errors in our observable data set then we do not have to worry too much. The error introduced by this assumption is actually dependent on the size of the population we are considering. The more individuals you have, you more the population will act like a continuous quantity. Since there are a lot of corpses out there, we do not think this assumption is too bad.

Note that we could model the motion of each zombie individually, however,  the computing power needed by such a simulation is much larger than the continuum description, which can be solved completely analytically. This is particularly important in the case of the zombie apocalypse, where time spent coding a simulation, may be better spent scavenging.

These are the basic assumptions we made when modelling a zombie population. Although I have tried to justify them you may have reservations about their validity. That is the very nature of mathematical modelling; try the simplest thing, first, and compare it to data. If you reproduce the phenomena that you are interested in then you have done your job well. However, if there is a discrepancy between the data and your maths then you have to revisit your assumptions and adapt them to make them more realistic.

Next time we contend with the equations and model the motion of the zombie as a random walker.

Diffusion of the dead - The maths of zombie invasions. Part 1, For those who can't wait.

Many years ago I posted a blog post about an academic article I had written about zombies. Finally, the article was published in Mathematical Modelling of Zombies. The book is designed to be readable by anyone with an interest in mathematics. However, those with an numerical background will find that it pushes them further as it does not shy away from clearly displaying the mathematics, whilst explaining the methods behind the madness.

As always, thank you to Martin Berube
for the use of his zombie image.

The articles range over a number of fields and are simply a means to dress up our everyday techniques in a way that is more palatable for a non-mathematical audience. Of course not all of you will want to shell out for the book. Thus, I've decided to essentially serialize the chapter in the next few posts, thus, we will be looking at all of the results of our paper and, hopefully, I'll be explaining the mathematics more clearly, so that any one is able to follow it. I may even throw in a matlab code or two, so that anyone is able to reproduce the results.

For those of you who are just interested in a quick review, you can read The Times article, or my brief version on the University of Oxford's Mathematical Institute website. Alternatively, if you are too tired to read you can always watch, or listen to a recorded version from the Athens Science Festival, Cambridge Science Festival, or on The Science of Fiction radio show. Finally, you could always see me live, when I'm giving one of my talks.

Next time, we will start at the beginning by modelling the zombie motion.