---
jupyter:
  description: 'Turing patterns with planex2d and NGSolve: Schnakenberg''s reaction-diffusion model, the growth rate of
    every wavelength from the linear analysis checked against the waves that grow on a square, the spots they settle
    into, and a tapering tail unrolled flat, whose tip carries rings where it is narrower than the shortest growing
    wavelength.'
---

# Turing patterns

Why can a spotted animal have a ringed tail? Alan Turing showed that two substances that react with each other and
diffuse at different rates can turn a uniform mixture into a pattern. The uniform state is unstable to waves of one band
of wavelengths, and the waves that grow fastest set the spacing of the spots or stripes that form. James Murray took
this as a model of the markings on an animal's coat and found that the shape of the skin decides the pattern. A tail so
thin that no wave around it can grow carries only rings. A spotted animal can therefore have a ringed tail, as the
cheetah and the genet have, but a striped animal cannot have a spotted one.

The example uses Schnakenberg's model, one of the standard two-substance reactions. Its linear analysis gives the
growth rate of every wavelength. NGSolve and SciPy then let a pattern grow from a little noise on a square, where the
growth rates of the waves are checked against the analysis, and on the skin of a tapering tail, unrolled flat, where
the spots give way to rings at the circumference the analysis predicts.

This page is part of a collection of [field and simulation examples](fields-and-simulation.md) using planex2d.

![A pattern growing out of noise on the skin of a tapering tail, unrolled and drawn eight times around, with spots on the wide part and rings at the tip](../_static/gallery/turing-patterns.gif){.result}

```python
import math

import matplotlib.pyplot as plt
import ngsolve as ngs
import numpy as np
import planex as px
import scipy.sparse as sp
from scipy.fft import dctn
from scipy.optimize import minimize_scalar
from scipy.sparse.linalg import splu
```

```python nbsphinx="hidden"
from docs_figures import MAP, Chart, Row
```

## The growth rates

The two concentrations $u$ and $v$ react and diffuse as

$$\frac{\partial u}{\partial t} = a - u + u^2 v + \nabla^2 u, \qquad \frac{\partial v}{\partial t} = b - u^2 v + d \, \nabla^2 v,$$

in the model's own units of length and time. With $a$ = 0.1, $b$ = 0.9 and $v$ diffusing $d$ = 20 times as fast as $u$,
the uniform state $u_0 = a + b$, $v_0 = b / (a + b)^2$ is stable without diffusion. A small wave of wavenumber $k$ on it
grows at the rate of the larger eigenvalue of $J - k^2 \operatorname{diag}(1, d)$, where $J$ is the Jacobian of the
reactions. It grows where the determinant of that matrix is negative, between the two roots
$k^2 = \left(p \pm \sqrt{p^2 - 4 d \det J}\right) / 2d$ with $p = d J_{11} + J_{22}$:

```python
a, b, d = 0.1, 0.9, 20
u0, v0 = a + b, b / (a + b) ** 2
J = np.array([[2 * u0 * v0 - 1, u0**2], [-2 * u0 * v0, -(u0**2)]])


def growth(k):
    """the growth rate of a wave with the wavenumber k on the uniform state"""
    return np.linalg.eigvals(J - k**2 * np.diag([1, d])).real.max()


kmax = minimize_scalar(lambda k: -growth(k), bounds=(0.1, 1), method="bounded", options={"xatol": 1e-10}).x
p = d * J[0, 0] + J[1, 1]
k1, k2 = (math.sqrt((p + s * math.sqrt(p**2 - 4 * d * np.linalg.det(J))) / (2 * d)) for s in (-1, 1))
wavelength = 2 * math.pi / kmax
print("fastest-growing wavenumber", kmax, "| rate", growth(kmax), "| wavelength", wavelength)
print("growing from k =", k1, "to", k2, "| shortest growing wavelength", 2 * math.pi / k2)
```

## Spots on a square

The square is ten wavelengths across, meshed by netgen with triangle sides of at most a twelfth of a wavelength. Each
time step takes the diffusion implicitly, with the mass matrix lumped onto the vertices, and the reactions explicitly,
so that the two matrices are factored once. A face with identified edges gets NGSolve's `Periodic` space, which the tail
below needs. The run starts from the uniform state with a noise of a millionth and returns $u$ at the times asked for:

```python
def simulate(face, times, dt=0.1, seed=1):
    """the mesh and u at the given times, grown from the uniform state with a little noise"""
    mesh = ngs.Mesh(px.OCCGeometry(face).GenerateMesh(maxh=wavelength / 12))
    V = ngs.H1(mesh, order=1)
    if face.identifications:
        V = ngs.Periodic(V)
    trial, test = V.TnT()
    free = np.flatnonzero(V.FreeDofs())

    def matrix(form):
        m = ngs.BilinearForm(form).Assemble().mat
        return sp.csr_matrix(m.CSR(), shape=(V.ndof, V.ndof))[free][:, free]

    M = np.asarray(matrix(trial * test * ngs.dx).sum(axis=1)).ravel()
    K = matrix(ngs.grad(trial) * ngs.grad(test) * ngs.dx)
    solve_u, solve_v = (splu((sp.diags(M) + dt * D * K).tocsc()).solve for D in (1, d))
    rng = np.random.default_rng(seed)
    u, v = (c * (1 + 1e-6 * rng.standard_normal(len(free))) for c in (u0, v0))
    found, t = [], 0.0
    for T in times:
        while t < T - dt / 2:
            uuv = u * u * v
            u, v = solve_u(M * (u + dt * (a - u + uuv))), solve_v(M * (v + dt * (b - uuv)))
            t += dt
        found.append(ngs.GridFunction(V))
        found[-1].vec.FV().NumPy()[free] = u
    return mesh, found


side = round(10 * wavelength)
square = px.QuickSketch().Rectangle(side, side).Face()
mesh, (early, later, emerging, settled) = simulate(square, (5, 15, 50, 800))
print(mesh.ne, "triangles")
```

On a square with insulated sides, every pattern is a sum of cosine waves, and while the pattern is small each wave
grows on its own. The cosine coefficients of $u$, taken on a grid of 256 by 256 points, give each wave's growth rate
between the times 5 and 15. The waves along one side lie on the curve of the analysis (Figure 1):

```python
def cosines(g, n=256):
    """the squared cosine coefficients of u on the square and their wavenumbers"""
    x = (np.arange(n) + 0.5) * side / n
    gx, gy = np.meshgrid(x, x)
    values = np.asarray(g(mesh(gx.ravel(), gy.ravel()))).reshape(n, n)
    m = np.arange(n) * math.pi / side
    return np.hypot(*np.meshgrid(m, m)), dctn(values, norm="ortho") ** 2


k, before = cosines(early)
_, after = cosines(later)
measured = np.log(after / before) / (2 * (15 - 5))
waves = (k[0] > k1) & (k[0] < k2)
predicted = np.array([growth(q) for q in k[0][waves]])
fast = predicted > 0.1
print(waves.sum(), "growing waves along one side; where faster than 0.1, largest difference",
      np.abs(measured[0][waves] - predicted)[fast].max())
# plt.plot(k[0][waves], measured[0][waves], "o"); plt.show()
```

```python nbsphinx="hide-input"
ks = np.linspace(0, 1, 400)
Chart("Wavenumber k", "Growth rate", lines={"Linear analysis": (ks, [growth(q) for q in ks])},
      dots={"Waves along one side of the square": (k[0][waves], measured[0][waves])},
      marks={"Fastest-growing wave": (kmax, growth(kmax))}, title="The growth rate of every wavelength")
```

Once the waves are large, the reactions couple them, and they settle into a hexagonal pattern of spots (Figure 2). The
settled spots lie closer together than the wavelength that grew fastest, though still inside the band of growing
wavenumbers:

```python
k, power = cosines(settled)
power[0, 0] = 0
print("strongest wavenumber at the end", k.ravel()[power.argmax()], "| fastest-growing", kmax)
# px.PlotField(settled, mesh, square)
```

```python nbsphinx="hide-input"
def sampled(mesh, g, turn=None):
    """u for Row and PlotField, at points of the face; turn: the angle of the tail's sector, into which points of a
    drawing around the tail are turned back"""
    def field(face, plot=True, **_):
        (x0, y0), (x1, y1) = face.bounding_box
        gx, gy = np.meshgrid(np.linspace(x0, x1, 500), np.linspace(y0, y1, 500))
        xy = np.c_[gx.ravel(), gy.ravel()]
        xy = xy[face.Distance(xy, signed=True) < -0.05]
        x, y = xy.T
        if turn is not None:
            r, angle = np.hypot(x, y), np.clip(np.arctan2(y, x) % turn, 1e-9, turn - 1e-9)
            x, y = r * np.cos(angle), r * np.sin(angle)
        return {"points": [tuple(p) for p in xy], "values": np.asarray(g(mesh(x, y))).ravel(), "label": "u"}
    return field


Row((square, "At the time 50", {"field": sampled(mesh, emerging), "cmap": MAP.reversed()}),
    (square, "At the time 800", {"field": sampled(mesh, settled), "cmap": MAP.reversed()}),
    caption="Spots growing out of noise on the square, and the settled pattern")
```

## A tapering tail

A tail is close to a cone, and the skin of a cone unrolls flat into a sector of a ring without stretching, so that
diffusion on it is diffusion on the flat sector. The two straight sides of the sector are the same line on the tail,
and identifying them closes the skin. The tail is 15 wavelengths long, 1.5 wavelengths around at the root and a fifth of
one at the tip:

```python
around_root, around_tip, length = 1.5 * wavelength, 0.2 * wavelength, 15 * wavelength
turn = (around_root - around_tip) / length
r_tip = around_tip / turn


def sector(turns=1):
    """the skin of the tail unrolled, or several turns of it side by side"""
    sk = px.QuickSketch((r_tip, 0)).Line(length, name="side_a")
    sk.Arc(center=(0, 0), sweep=turns * math.degrees(turn), name="root").Turn(90).Line(length, name="side_b")
    return sk.Arc(center=(0, 0), sweep=-turns * math.degrees(turn), name="tip").Face()


tail = sector()
tail.Identify("side_a", "side_b")
print("sector of", math.degrees(turn), "degrees | area", tail.mass,
      "| cone's side", (around_root + around_tip) / 2 * length)
tail_mesh, (skin,) = simulate(tail, (800,))
```

At each distance from the tip, the circle around the tail shows whether it carries a ring, with the same value all the
way around, or spots. The rings end close to where the circumference reaches the shortest wavelength that can grow
(Figure 3):

```python
def spread(r):
    """how much u varies around the tail at the radius r of the sector"""
    angle = np.linspace(0, turn, 60)
    values = np.asarray(skin(tail_mesh(r * np.cos(angle), r * np.sin(angle)))).ravel()
    return values.max() - values.min()


radii = np.arange(r_tip + 0.5, r_tip + length / 2, 0.25)
first = next(r for r in radii if spread(r) > 0.1)
print("rings up to the circumference", (first - 0.25) * turn, "| spots from", first * turn,
      "| shortest growing wavelength", 2 * math.pi / k2)
# px.PlotField(skin, tail_mesh, tail)
```

```python nbsphinx="hide-input"
seams = [px.Line((r_tip * math.cos(j * turn), r_tip * math.sin(j * turn)),
                 ((r_tip + length) * math.cos(j * turn), (r_tip + length) * math.sin(j * turn))) for j in range(1, 8)]
Row((sector(8), "The skin of the tail, eight turns side by side",
     {"field": sampled(tail_mesh, skin, turn), "cmap": MAP.reversed(), "dashed": seams}),
    caption="The tail's skin unrolled, its sector drawn eight times around")
```

See [Identifications](../integrations/netgen/identifications.md) for identified edges and the rotational sector, and
[Chladni figures](chladni-figures.md) for patterns set by the shape of a plate.

```python nbsphinx="hidden"
import sys
from matplotlib.animation import FuncAnimation
sys.path.insert(0, "..")
from docs_figures import INK
from gallery import save_animation

# the pattern growing out of the noise on the tail, at times spaced evenly on a log scale: the waves grow
# exponentially until the spots form, then settle
times = np.geomspace(30, 800, 30)
grown_mesh, grown = simulate(tail, times)
fan = sector(8)
(x0, y0), (x1, y1) = fan.bounding_box
gx, gy = np.meshgrid(np.linspace(x0, x1, 520), np.linspace(y0, y1, 400))
xy = np.c_[gx.ravel(), gy.ravel()]
inside = fan.Distance(xy, signed=True) < -0.05
r, angle = np.hypot(xy[inside, 0], xy[inside, 1]), np.clip(np.arctan2(xy[inside, 1], xy[inside, 0]) % turn, 1e-9, turn - 1e-9)
where = grown_mesh(r * np.cos(angle), r * np.sin(angle))
height = (x1 - x0 + 8) / 1.5
fig, ax = plt.subplots(figsize=(6, 4))


def frame(k):
    ax.clear()
    image = np.full(gx.size, np.nan)
    image[inside] = np.asarray(grown[min(k, len(grown) - 1)](where)).ravel()
    ax.imshow(image.reshape(gx.shape), origin="lower", extent=(x0, x1, y0, y1), cmap=MAP.reversed(), vmin=0.2, vmax=3.1,
              interpolation="bilinear")
    fan.Plot(ax=ax, face_col=None, line_col=INK, arc_col=INK, legend=False)
    ax.set_aspect("equal", adjustable="box")
    ax.set_xlim(x0 - 4, x1 + 4)
    ax.set_ylim((y0 + y1) / 2 - height / 2, (y0 + y1) / 2 + height / 2)
    ax.axis("off")


animation = FuncAnimation(fig, frame, frames=len(grown) + 14, interval=80)
save_animation(animation, "turing-patterns", fps=12, bare=True)
plt.close(fig)
```