- Let's import the packages.

In [ ]:

```
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
```

- We will simulate the following system of partial differential equations on the domain $E=[-1,1]^2$:

\begin{align*}
\frac{\partial u}{\partial t} &= a \Delta u + u - u^3 - v + k\\
\tau\frac{\partial v}{\partial t} &= b \Delta v + u - v\\
\end{align*}

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.

In [ ]:

```
a = 2.8e-4
b = 5e-3
tau = .1
k = -.005
```

- We discretize time and space. The following condition ensures that the discretization scheme we use here is stable:

In [ ]:

```
size = 80 # size of the 2D grid
dx = 2./size # space step
```

In [ ]:

```
T = 10.0 # total time
dt = .9 * dx**2/2 # time step
n = int(T/dt)
```

In [ ]:

```
U = np.random.rand(size, size)
V = np.random.rand(size, size)
```

- Now, we define a function that computes the discrete Laplace operator of a 2D variable on the grid, using a five-point stencil finite difference method. This operator is defined by:

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.

In [ ]:

```
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
```

In [ ]:

```
# 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]
```

- Finally, we display the variable $u$ after a time $T$ of simulation.

In [ ]:

```
plt.imshow(U, cmap=plt.cm.copper, extent=[-1,1,-1,1]);
```

You'll find all the explanations, figures, references, and much more in the book (to be released later this summer).

IPython Cookbook, by Cyrille Rossant, Packt Publishing, 2014 (500 pages).