Turing instability of the Brusselator (Prigogine & Lefever 1968)¶
Family: patterns · Script: examples/patterns/brusselator.py
Obtained |
b_c to ≤ 9e-16 relative (3 sets) |
Reference |
Prigogine & Lefever 1968; Peña & Pérez-García 2001, PRE 64 |
dt u = Du dx^2 u + a - (b + 1) u + u^2 v
dt v = Dv dx^2 v + b u - u^2 v
The uniform state (u, v) = (a, b/a) loses stability to a stationary pattern of wavenumber k_c at
b_c = (1 + a eta)^2, k_c^2 = a / sqrt(Du Dv), eta = sqrt(Du / Dv)
(Peña & Pérez-García 2001, Phys. Rev. E 64, 056213, Sec. II), provided b_c < 1 + a^2 (the Hopf
threshold of the k = 0 mode), i.e. Dv / Du large enough. The nonlinear IVP is written as in the
article and linearized by the IR around the uniform state; a perturbation exp(i k x + sigma t)
lives in one Fourier group (2x2 pencil). b enters linearly: the parametric pencil gives the neutral
b(k) by Newton on Re(sigma), minimized over k (the domain length 2 pi / k is rebuilt per k).
Run: python examples/patterns/brusselator.py
Run it: python examples/patterns/brusselator.py (in the conatus environment).
1"""Turing instability of the Brusselator (Prigogine & Lefever 1968).
2
3 dt u = Du dx^2 u + a - (b + 1) u + u^2 v
4 dt v = Dv dx^2 v + b u - u^2 v
5
6The uniform state (u, v) = (a, b/a) loses stability to a stationary pattern of wavenumber k_c at
7
8 b_c = (1 + a eta)^2, k_c^2 = a / sqrt(Du Dv), eta = sqrt(Du / Dv)
9
10(Peña & Pérez-García 2001, Phys. Rev. E 64, 056213, Sec. II), provided b_c < 1 + a^2 (the Hopf
11threshold of the k = 0 mode), i.e. Dv / Du large enough. The nonlinear IVP is written as in the
12article and linearized by the IR around the uniform state; a perturbation exp(i k x + sigma t)
13lives in one Fourier group (2x2 pencil). b enters linearly: the parametric pencil gives the neutral
14b(k) by Newton on Re(sigma), minimized over k (the domain length 2 pi / k is rebuilt per k).
15
16Run: python examples/patterns/brusselator.py
17"""
18import numpy as np
19import scipy.optimize as so
20
21from conatus.expert import (IVP, CartesianCoordinates, ComplexFourier, Domain, Field, Parameter,
22 ParametricEVP, critical, dt, lap, linearize)
23
24
25def brusselator(k: float, a: float, Du: float, Dv: float, b0: float):
26 """Parametric pencil in b of the Brusselator linearized around (a, b/a), wavenumber k."""
27 cs = CartesianCoordinates("x")
28 x = cs["x"]
29 D = Domain(cs, [ComplexFourier(x, 4, (0, 2 * np.pi / k), dealias=1)]) # modes 0, +-k
30 u, v = Field("u", D), Field("v", D)
31 A, B = Parameter("a", a), Parameter("b", b0)
32 ivp = IVP([u, v], name="Brusselator")
33 ivp.add_equation((dt(u) - Du * lap(u) + (B + 1) * u, A + u * u * v))
34 ivp.add_equation((dt(v) - Dv * lap(v) - B * u, -u * u * v))
35 evp = linearize(ivp, around={u: A, v: B / A}, eigenvalue=Parameter("sigma", dtype=complex))
36 return ParametricEVP(evp, [B], group=lambda g: g[0].k > 0), evp
37
38
39def turing_threshold(a: float, Du: float, Dv: float):
40 """(b_c, k_c): neutral b(k) of the least stable mode, minimized over k."""
41 b_guess = (1 + a * np.sqrt(Du / Dv)) ** 2
42 crit = {}
43
44 def neutral_b(k):
45 P, _ = brusselator(k, a, Du, Dv, b_guess)
46 crit[k] = critical(P, "b", bracket=(0.2 * b_guess, 1 + a * a), growth="real", xtol=1e-9)
47 return crit[k].value
48
49 kg = np.sqrt(a / np.sqrt(Du * Dv))
50 opt = so.minimize_scalar(neutral_b, bounds=(0.7 * kg, 1.3 * kg), method="bounded",
51 options={"xatol": 1e-10})
52 return neutral_b(float(opt.x)), float(opt.x), crit[float(opt.x)]
53
54
55if __name__ == "__main__":
56 for a, Du, Dv in ((2.0, 1.0, 8.0), (3.0, 1.0, 20.0), (1.5, 2.0, 16.0)):
57 bc, kc, cp = turing_threshold(a, Du, Dv)
58 bth, kth = (1 + a * np.sqrt(Du / Dv)) ** 2, np.sqrt(a / np.sqrt(Du * Dv))
59 print(f"a={a}, Du={Du}, Dv={Dv}: b_c = {bc:.12f} (theory {bth:.12f}, gap {bc - bth:.1e}), "
60 f"k_c = {kc:.8f} (theory {kth:.8f}), sigma = {cp.mode.sigma:.1e}, residual {cp.mode.residual:.1e}")
61 assert abs(bc - bth) < 1e-12 and abs(kc - kth) < 1e-6, (a, bc, bth, kc, kth) # closed form
62 print(cp.explain())
63 print(brusselator(kc, a, Du, Dv, bc)[1].explain())