MATH3104 - Mathematical Biology
Primer - Finite Differences Method
Dietmar Oelz
May 10, 2021
Here we focus on straightforward ways to treat numerically reaction-drift-diffusion equations
in 1D. Typically we are looking for solutions u = u(t, x) where t ≥ 0 is time and x ∈ [a, b] ⊂ R.
In the two assignments for this part of MATH3104 we’ll restrict to periodic boundary con-
ditions, i.e. ρ(t, a) = ρ(t, b) and ∂xρ(t, a) = ∂xρ(t, b), therefore all the differential operators and
equations we will discuss in this primer are implicitly meant to be coupled to periodic b.c. .
1 Discrete approximation of functions
Numerically, i.e. on the computer, we cannot represent functions on intervals, but only vectors.
To approximate a function u = u(x) by a vector u¯ = (u1, u2, ..., uN ) we may use an equi-distant
grid with grid-size ∆x = (b− a)/(N − 1) and N grid points x1 = a, ..., xi = a+ i∆x,...,xN = b
such that ui ≈ u(xi).
For time-dependent problems we write the solution u(t, x) as u¯n = (un1 , ..., u
n
N ) where n =
0, 1, 2... is an index of time such that uni ≈ u(n∆t, i∆x). Here ∆t > 0 is the size of the time-steps.
2 Finite differences
The most straightforward way to discretise derivatives is ∂xu(xi) ≈ ui+1−ui∆x (forward difference
operator), ∂xu(xi) ≈ ui−ui−1∆x (backward difference operator) or ∂xu(xi) ≈ ui+1−ui−12∆x (central
difference operator)
For second derivatives the easiest approach is to combine a forward and a backward difference:
∂xxu(xi) ≈
ui+1−ui
∆x − ui−ui−1∆x
∆x
=
ui+1 − 2ui + ui−1
(∆x)2
3 Discrete diffusion operator in 1D
The linear diffusion operator u 7→ D∂xx assigns the 2nd derivative (times D) to a given function.
Analogously the discrete diffusion operator is a linear map - a matrix - applied to a vector of
grid values u¯ = (u1, ..., uN ), namely
u1...
un
 7→ 1
(∆x)2

u2 − 2u1 + uN
u3 − 2u2 + u1
...
uN − 2uN−1 + uN−2
u1 − 2uN + uN−1
 .
1
In matrix notation we obtain u¯ = (u1, ..., uN ) 7→ DMu¯, where
M =
1
(∆x)2

−2 1 0 0 0 · · · 0 1
1 −2 1 0 0 · · · 0 0
0 1 −2 1 0 · · · 0 0
...
. . .
...
0 0 0 · · · 0 1 −2 1
1 0 0 · · · 0 0 1 −2

.
Note that the periodic boundary condition implies that xN is taken as the “left” neighbour of
x1, which itself is taken as the “right” neighbour of xN .
4 Linear Elliptic equation (Poisson, etc)
To solve a stationary (not time-dependent) equation such as 0 = D∆u− αu+ f solve
0 = DMu¯− αu¯+ f¯
i.e. u¯ = (DM − α1)−1 (−f¯), where f¯ is a vector corresponding to the discretised version of the
right hand side function f and 1 is the identity matrix.
5 Diffusion equation
If we add dependence on time then the canonical 2nd order PDE is the diffusion equation
∂tρ = D∂xxρ
with initial datum ρ(t = 0, .) = ρI(.) and diffusion coefficient D > 0.
1. Start to iterate in time with the vector u¯0 which is a discretised version of ρI . For ρI =
δ(x− c) (Dirac-δ, i.e. one particle located at x = c), take
u¯0i =
{
1
∆x i = i0
0 otherwise
where i0 is the index which corresponds to the position of the δ-peak given by c. Note
that this guarantees that integration of u¯0 using a straightforward quadrature formulate
such as
∑
i u¯
0
i∆x (in Matlab: sum(ubar0)*deltax) returns 1 just as for the Dirac delta.
2. To compute u¯1 given u¯0, and then u¯2 given u¯1, etc. (time-stepping) use schemes such as
(a) Explicit Euler: u¯
n−u¯n−1
∆t = DMu¯
n−1, hence u¯n = u¯n−1 + ∆tDMu¯n−1 Attention:
make sure to satisfy the diffusive CFL-condition, D∆t/∆x2 ≤ 1/2, otherwise the
scheme is unstable.
Note that if you solve the diffusion equation starting from a delta-Distribution with
this explicit scheme and the optimal value of ∆t according to the diffusive CFL-
condition, you may obtain zeros in every second grid point. In this case, to avoid
that, just set ∆t to a slightly smaller value.
(b) Implicit euler u¯
n−u¯n−1
∆t = DMu¯
n, hence u¯n = (1−∆tDM)−1 u¯n−1. Note that im-
plicit schemes are typically more stable than explicit ones, less restrictive when it
comes to step sizes, but harder to implement.
(c) or a midpoint rule such as u¯
n−u¯n−1
∆t = D
1
2M(u¯
n + u¯n−1), etc.
2
Use may use the following Matlab code as a skeleton for your numerical codes. It approx-
imates the solution of 
∂tρ = D∆∂xxρ ,
ρ(0) = ρ(10) ,
ρ(t = 0, x) = 100× δ(x− 5) .
a=0; b=10;
Dcoef=0.1;
N=100; deltax=(b-a)/N; deltat=0.1;
xval=linspace(a, b-deltax, N);
M=diff(diff([[zeros(1, N-1), 1]; eye(N); [1, zeros(1, N-1)]]))/deltax^2;
h=figure;
rho=zeros(N, 1); rho(N/2)=100/deltax;
t=0;
while t<100
t=t+deltat;
rhoold=rho;
rho=(eye(N)-Dcoef*deltat*M)\rhoold;
figure(h);
plot(xval, rho);
yl=ylim;
ylim([0, 100]);
drawnow
end
6 Scalar conservation laws (1D)
∂tu+ ∂x(f(u)) = 0 .
The simplest scheme (1st order) is the Lax-Friedrichs method, i.e.
u¯ni − 12
(
u¯n−1i+1 + u¯
n−1
i−1
)
∆t
+
f(un−1i+1 )− f(un−1i−1 )
2∆x
= 0 ,
which implies that
u¯ni =
1
2
(
u¯n−1i+1 + u¯
n−1
i−1
)− ∆t
2∆x
(
f(un−1i+1 )− f(un−1i−1 )
)
. (1)
Attention: The following CFL-condition,
max
i=1..N
∣∣∣∣f(u¯n−1i )u¯n−1i
∣∣∣∣∆t∆x ≤ 1/2 , (2)
guarantees the stability of the scheme. Also bear in mind that this scheme generates significant
numerical diffusion, unless the left hand side in (2) is close to its upper bound and unless the
grid is sufficiently fine (∆t, ∆x small). Therefore it might be a good idea to adapt the step-size
∆t in every time-step.
Note that as before the periodic boundary condition is implemented identifying xN+1 with x1
and x0 with xN . To obtain the shifted vectors (u¯i±1), you may want to use MATLAB-
circshift which implements the periodic boundary “automatically”.
3
circshift([1.0; 2.0; 4.0], 1)=[4.0; 1.0; 2.0]
circshift([1.0; 2.0; 4.0], -1)=[2.0; 4.0; 1.0]
and so u¯i±1 = circshift(u¯i,∓1).
E.g. if the flux is given by f(u) = u∂xS, then, starting with a discrete version of S, use
circshift to obtain a finite difference approximation of ∂xS, here written as (∆S)i, and then set
f(un−1i ) = ui(∆S)i.
7 Putting things together: reaction-drift-diffusion
Imagine you want to combine diffusion and something else such as reaction and/or drift,
∂tu+ ∂x(f(u)) = D∂xxu+ g(u) .
1. One approach is to compute an explicit time stepping of the entire equation,
u¯ni =
1
2
(
u¯n−1i+1 + u¯
n−1
i−1
)− ∆t
2∆x
(
f(un−1i+1 )− f(un−1i−1 )
)
+ (∆tDMu¯n−1)i + ∆t g(un−1i ) .
Depending on whether drift or diffusion dominates, you will need to satisfy either the
diffusive or the hyperbolic CFL conditions (and be aware of the magnitude of numerical
diffusion).
2. Fully implicit scheme such as u¯
n−u¯n−1
∆t = DMu¯
n + g(u¯n, x, t) which needs to be solved
either exactly using Newton iteration (if f is non-linear) or approximately (only the first
or a few Newton iterations).
3. You might also want to convert some u¯n−1 into u¯n to obtain (semi-)implicit schemes with
better stability. You might try (here for the reaction-diffusion equation ∂tu = ∆u+ u(1−
u)),
u¯n − u¯n−1
∆t
= Mu¯n + u¯n−1(1− u¯n)
or
u¯n − u¯n−1
∆t
= Mu¯n + u¯n(1− u¯n−1) ,
making sure that the resulting equation is linear in u¯n.
4. Another approach - the one I recommend for its simplicity - is to use a splitting
method where - within one time-step - you compute e.g.
(a) one time-step of a scheme for the diffusion equation (see section 5) such as u¯n =
(1−∆tDM)−1 u¯n−1 (fully implicit).
(b) then rename (u¯n−1 = u¯n!) and compute one time-step of an ODE-scheme for the
reaction part of the equation such as
u¯ni = u¯
n−1
i + ∆t u¯
n−1
i (1− u¯ni ) ,
i.e. u¯ni =
(
u¯n−1i (1 + ∆t)
)
/
(
1 + ∆t u¯n−1i
)
.
(c) then rename (u¯n−1 = u¯n!) and compute one time-step for the drift part (assignment
8) such as Lax-Friedrichs Method (1).
4
8 Systems of equations
Again the splitting scheme is the easiest way to deal with systems of equations. For the parabolic-
elliptic Keller-Segel system, in every time step you may want to
1. obtain the vector (S
n− 1
2
i )i solving the elliptic equation for the chemoattractant given the
right hand side u¯n−1 (see section 4),
2. and then compute u¯n as time-update of the drift-diffusion equation given the chemoat-
tractant concentration (S
n− 1
2
i )i (see section 7).
Likewise, for the Turing system in every time-step you may want to
1. find (u
n− 1
2
i ) computing one Euler step with the reaction term f(u
n−1
i , v
n−1
i ) and then
compute (uni ) by solving diffusion equation.
2. and then apply the same procedure (vn−1i )i to obtain (v
n
i )i.
5











































































































































































































































































学霸联盟