Software 14.4.1 (Exact and linear waves). A symbolic check that \(1/(1 + e^{z/\sqrt 2})\) is a bistable front, solving the residual for the speed; then the roots of the linearisation at \(u = 0\) for the Fisher equation at speeds \(1\), \(2\), \(3\):

>>> import sympy as sp
>>> x, t, a, c = sp.symbols('x t a c', real=True)
>>> lam = sp.symbols('lambda')
>>> z = x - c*t  # bistable: u_t = u_xx + u(1-u)(u-a)
>>> U = 1 / (1 + sp.exp(z / sp.sqrt(2)))
>>> res = U.diff(t) - U.diff(x, 2) - U * (1 - U) * (U - a)
>>> sp.solve(sp.simplify(res * (1 + sp.exp(z/sp.sqrt(2)))**3
...                      / sp.exp(z/sp.sqrt(2))), c)        # the wave speed
[sqrt(2)*(1/2 - a)]
>>> for cc in (1, 2, 3):                   # Fisher: linearised at u = 0
...     print(cc, sp.solve(lam**2 + cc*lam + 1, lam))
...
1 [-1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2]
2 [-1]
3 [-3/2 - sqrt(5)/2, -3/2 + sqrt(5)/2]

Reading the output. The residual vanishes exactly when \(c = \sqrt 2(\frac 12 - a)\). For the Fisher front, \(c = 1\) gives complex roots, so a profile with that speed would oscillate about \(0\); \(c = 2\) gives the double root \(-1\), the boundary case; \(c = 3\) gives two negative roots, a node.

MATLAB/Octave: roots([1 c 1]) for each c.

Software 14.4.2 (The Fisher front). The Fisher–KPP equation solved by explicit differences on \([0, 200]\) with \(\Delta x = 0.05\) and \(r = 0.2\), from \(u = 1\) on \(x < 5\); the position where \(u = \frac 12\) at four times, the average speeds over the intervals between them, and Bramson’s prediction \(2 - 3/2t\) at the mid-times:

>>> import numpy as np
>>> L, N = 200.0, 4000                     # Fisher-KPP: u_t = u_xx + u(1 - u)
>>> x = np.linspace(0, L, N + 1); dx = x[1]; dt = 0.2 * dx**2
>>> u = np.where(x < 5, 1.0, 0.0)          # a population released on the left
>>> def front(u): return x[np.argmax(u < 0.5)]          # where u crosses 1/2
>>> pos = {}; T = 0.0
>>> for target in (20, 40, 60, 80):
...     while T < target - 1e-12:
...         lap = np.r_[0, u[2:] - 2*u[1:-1] + u[:-2], 0] / dx**2
...         u = u + dt * (lap + u * (1 - u)); u[0] = 1.0; u[-1] = 0.0
...         T += dt
...     pos[target] = round(front(u), 2)
...
>>> pos
{20: 38.3, 40: 77.2, 60: 116.6, 80: 156.2}
>>> [round((pos[b] - pos[a]) / 20, 4)
...  for a, b in ((20, 40), (40, 60), (60, 80))]   # front speeds
[1.945, 1.97, 1.98]
>>> [round(2 - 3 / (2*tm), 4) for tm in (30, 50, 70)]   # Bramson
[1.95, 1.97, 1.9786]

Reading the output. The front advances \(38.9\), \(39.4\) and \(39.6\) in successive intervals of \(20\), speeds \(1.945\), \(1.970\), \(1.980\), approaching \(2\) from below; Bramson’s formula gives \(1.950\), \(1.970\), \(1.979\), in agreement to about \(0.005\). The slow logarithmic approach to the minimal speed is a hallmark of pulled fronts.

MATLAB/Octave: the same explicit loop; find(u < 0.5, 1) for the front.

Software 14.4.3 (A Turing pattern). For the Schnakenberg model, the largest growth rate over \(0 \leq k \leq 2\) and the wavenumber attaining it for \(d = 1\), \(10\), \(40\); then a simulation with \(d = 40\) on a periodic line of length \(120\), by the Fourier method with diffusion treated implicitly, and the dominant wavenumber of the final state, its wavelength, and the range of \(u\):

>>> A, B = 0.1, 0.9  # Schnakenberg kinetics, steady state
>>> us, vs = A + B, B / (A + B)**2
>>> J = np.array([[-1 + 2*us*vs, us**2], [-2*us*vs, -us**2]])
>>> def growth(k, d):  # largest Re eigenvalue of J - k^2 D
...     return np.linalg.eigvals(J - k**2 * np.diag([1.0, d])).real.max()
...
>>> ks = np.linspace(0, 2, 2001)
>>> for d in (1.0, 10.0, 40.0):
...     g = np.array([growth(k, d) for k in ks])
...     print(d, round(g.max(), 4), round(ks[g.argmax()], 3))
...
1.0 -0.1 0.0
10.0 0.0572 0.564
40.0 0.411 0.421
>>> Lp, M, d = 120.0, 1024, 40.0  # simulate d = 40 on a periodic line
>>> xp = np.arange(M) * Lp / M; kk = 2 * np.pi * np.fft.fftfreq(M, d=Lp/M)
>>> rng = np.random.default_rng(0)
>>> U = us + 0.01 * rng.standard_normal(M)
>>> V = vs + 0.01 * rng.standard_normal(M)
>>> h = 0.05  # diffusion implicit, reactions explicit
>>> for _ in range(int(400 / h)):
...     f = A - U + U**2 * V; g = B - U**2 * V
...     U = np.fft.ifft(np.fft.fft(U + h*f) / (1 + h*kk**2)).real
...     V = np.fft.ifft(np.fft.fft(V + h*g) / (1 + h*d*kk**2)).real
...
>>> spec = np.abs(np.fft.fft(U - U.mean()))[:M//2]
>>> kdom = kk[np.argmax(spec)]
>>> kdom, round(2*np.pi / kdom, 2), round(U.min(), 3), round(U.max(), 3)
(0.41887902047863906, 15.0, 0.149, 2.696)

Reading the output. With equal diffusivities every mode decays (the largest rate is \(-0.1\), at \(k = 0\)). At \(d = 10\), just above the threshold \(8.57\), a narrow band grows slowly, at rate \(0.057\); at \(d = 40\) the fastest mode has \(k = 0.421\). The simulation settles on \(k = 0.419\), the allowed wavenumber \(2\pi \cdot 8/120\) nearest to it: eight peaks of wavelength \(15\), with \(u\) between \(0.15\) and \(2.70\).

MATLAB/Octave: max(real(eig(J - k^2*diag([1 d])))) over a grid of k.