This is a featured recipe from the IPython Cookbook, the definitive guide to high-performance scientific computing and data science in Python.
Partial Differential Equations (PDEs) describe the evolution of dynamical systems involving both time and space. Examples in physics include sound, heat, electromagnetism, fluid flow, elasticity, among others. Examples in biology include tumor growth, population dynamics, and epidemic propagations.
PDEs are hard to solve analytically. Therefore, PDEs are often studied via numerical simulations.
In this featured recipe, we will illustrate how to simulate a reaction-diffusion system described by a PDE called the Fitzhugh–Nagumo equation. A reaction-diffusion system models the evolution of one or several variables subject to two processes: reaction (transformation of the variables into each other) and diffusion (spreading across a spatial region). Some chemical reactions may be described by this type of model, but there are other applications in physics, biology, ecology, and other disciplines.
Here, we simulate a system that has been proposed by Alan Turing as a model of animal coat pattern formation. Two chemical substances influencing skin pigmentation interact according to a reaction-diffusion model. This system is responsible for the formation of patterns that are reminiscent of the pelage of zebras, jaguars, and giraffes.
We will simulate this system with the finite difference method. This method consists of discretizing time and space, and replacing the derivatives by their discrete equivalents.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
The variable $u$ represents the concentration of a substance favoring skin pigmentation, whereas $v$ represents another substance that reacts with the first and impedes pigmentation.
At initialization time, we assume that $u$ and $v$ contain independent random numbers on every grid point. Besides, we take Neumann boundary conditions: we require the spatial derivatives of the variables with respect to the normal vectors to be null on the boundaries of the domain $E$.
Let's define the four parameters of the model.
a = 2.8e-4
b = 5e-3
tau = .1
k = -.005
size = 100 # size of the 2D grid
dx = 2./size # space step
T = 10.0 # total time
dt = .9 * dx**2/2 # time step
n = int(T/dt)
U = np.random.rand(size, size)
V = np.random.rand(size, size)
We can compute the values of this operator on the grid using vectorized matrix operations. Because of side effects on the edges of the matrix, we need to remove the borders of the grid in the computation.
def laplacian(Z):
Ztop = Z[0:-2,1:-1]
Zleft = Z[1:-1,0:-2]
Zbottom = Z[2:,1:-1]
Zright = Z[1:-1,2:]
Zcenter = Z[1:-1,1:-1]
return (Ztop + Zleft + Zbottom + Zright - 4 * Zcenter) / dx**2
# We simulate the PDE with the finite difference method.
for i in range(n):
# We compute the Laplacian of u and v.
deltaU = laplacian(U)
deltaV = laplacian(V)
# We take the values of u and v inside the grid.
Uc = U[1:-1,1:-1]
Vc = V[1:-1,1:-1]
# We update the variables.
U[1:-1,1:-1], V[1:-1,1:-1] = \
Uc + dt * (a * deltaU + Uc - Uc**3 - Vc + k), \
Vc + dt * (b * deltaV + Uc - Vc) / tau
# Neumann conditions: derivatives at the edges
# are null.
for Z in (U, V):
Z[0,:] = Z[1,:]
Z[-1,:] = Z[-2,:]
Z[:,0] = Z[:,1]
Z[:,-1] = Z[:,-2]
plt.imshow(U, cmap=plt.cm.copper, extent=[-1,1,-1,1]);
plt.xticks([]); plt.yticks([]);
Whereas the variables when completely random at initialization time, we observe the formation of patterns after a sufficiently long simulation time.
Let's explain how the finite difference method allowed us to implement the update step. We start from the following system of equations:
We first use the following scheme for the discrete Laplace operator:
We also use this scheme for the time derivative of $u$ and $v$:
We end up with the following iterative update step:
Here, our Neumann boundary conditions state that the spatial derivatives with respect to the normal vectors are null on the boundaries of the domain E:
We implement these boundary conditions by duplicating values in matrices $U$ and $V$ on the edges (see code).
In order to ensure that our numerical scheme converges to a numerical solution that is close to the actual (unknown) mathematical solution, the stability of the scheme needs to be ascertained. One can show that, here, a sufficient condition for the stability is:
Here are further references on partial differential equations, reaction-diffusion systems, and numerical simulations of those systems.
This was a featured recipe from the IPython Cookbook, by Cyrille Rossant, Packt Publishing, 2014.