Complex Ginzburg–Landau equation: Benjamin–Feir–Newell and Eckhaus instabilities of plane waves¶
Family: patterns · Script: examples/patterns/ginzburg_landau.py
Obtained |
complex form (conj() on the LHS, doubled system declared); BFN threshold: 1 + c₁c₃ < 1e-12 (Richardson k → 0); Eckhaus q_E² to 1e-7 |
Reference |
Newell 1974; Aranson & Kramer 2002, RMP 74 |
dt A = A + (1 + i c1) dx^2 A - (1 + i c3) |A|^2 A
(Aranson & Kramer 2002, Rev. Mod. Phys. 74, 99, eq. (1), with b = c1, c = c3). The plane wave
A = R exp(i (q x - omega t)), R^2 = 1 - q^2, omega = c3 R^2 + c1 q^2, is linearly unstable to long
waves (Benjamin–Feir–Newell) iff 1 + c1 c3 < 0 (Newell 1974; Aranson & Kramer 2002, Sec. II.B);
for 1 + c1 c3 > 0 it is Eckhaus-unstable iff q^2 > q_E^2 = (1 + c1 c3) / (3 + c1 c3 + 2 c3^2).
Written in the co-moving, co-rotating variable B = A exp(-i (q x - omega t)), the plane wave is the
constant state B = R and the equation keeps constant coefficients:
dt B = (1 + i omega) B + (1 + i c1)(dx^2 B + 2 i q dx B - q^2 B) - (1 + i c3) |B|^2 B.
The IR linearizes the complex form: conj(dB) appears (|B|^2 B = B B conj(B)), so the linear
operator is R-linear, not C-linear, and couples the Fourier modes k and -k. The discretizer solves
it as the equivalent doubled system in (dB, conj dB), every equation followed by its conjugate, the
eigenvalue standing for d/dt (declared in explain(); before conj() in an LHS was refused
and this example was written in the real variables B = X + i Y: same spectrum, checked in
tests/discretize/test_conj_lhs.py). A perturbation exp(i k x + sigma t) lives in one Fourier group
(2x2 pencil); the long-wave threshold is the k -> 0 limit.
Run: python examples/patterns/ginzburg_landau.py
Run it: python examples/patterns/ginzburg_landau.py (in the conatus environment).
1"""Complex Ginzburg–Landau equation: Benjamin–Feir–Newell and Eckhaus instabilities of plane waves.
2
3 dt A = A + (1 + i c1) dx^2 A - (1 + i c3) |A|^2 A
4
5(Aranson & Kramer 2002, Rev. Mod. Phys. 74, 99, eq. (1), with b = c1, c = c3). The plane wave
6A = R exp(i (q x - omega t)), R^2 = 1 - q^2, omega = c3 R^2 + c1 q^2, is linearly unstable to long
7waves (Benjamin–Feir–Newell) iff 1 + c1 c3 < 0 (Newell 1974; Aranson & Kramer 2002, Sec. II.B);
8for 1 + c1 c3 > 0 it is Eckhaus-unstable iff q^2 > q_E^2 = (1 + c1 c3) / (3 + c1 c3 + 2 c3^2).
9
10Written in the co-moving, co-rotating variable B = A exp(-i (q x - omega t)), the plane wave is the
11constant state B = R and the equation keeps constant coefficients:
12
13 dt B = (1 + i omega) B + (1 + i c1)(dx^2 B + 2 i q dx B - q^2 B) - (1 + i c3) |B|^2 B.
14
15The IR linearizes the complex form: conj(dB) appears (|B|^2 B = B B conj(B)), so the linear
16operator is R-linear, not C-linear, and couples the Fourier modes k and -k. The discretizer solves
17it as the equivalent doubled system in (dB, conj dB), every equation followed by its conjugate, the
18eigenvalue standing for d/dt (declared in explain(); before conj() in an LHS was refused
19and this example was written in the real variables B = X + i Y: same spectrum, checked in
20tests/discretize/test_conj_lhs.py). A perturbation exp(i k x + sigma t) lives in one Fourier group
21(2x2 pencil); the long-wave threshold is the k -> 0 limit.
22
23Run: python examples/patterns/ginzburg_landau.py
24"""
25import numpy as np
26import scipy.optimize as so
27
28from conatus.expert import (IVP, CartesianCoordinates, ComplexFourier, Domain, Field, Parameter,
29 ParametricEVP, conj, critical, diff, discretize, dt, eigenvalues, linearize,
30 sqrt)
31
32
33def cgle(k: float, c1: float, c3: float, q: float):
34 """Linearized CGLE around the plane wave of wavenumber q, one perturbation wavenumber k."""
35 cs = CartesianCoordinates("x")
36 x = cs["x"]
37 D = Domain(cs, [ComplexFourier(x, 4, (0, 2 * np.pi / k), dealias=1)]) # modes 0, +-k
38 B = Field("B", D, dtype=complex)
39 C1, C3, Q = Parameter("c1", c1), Parameter("c3", c3), Parameter("q", q)
40 R2 = 1 - Q ** 2 # plane-wave amplitude squared
41 w = C3 * R2 + C1 * Q ** 2 # plane-wave frequency
42
43 def disp(f): # (1 + i c1)(f'' + 2 i q f' - q^2 f)
44 return (1 + 1j * C1) * (diff(diff(f, x), x) + 2j * Q * diff(f, x) - Q ** 2 * f)
45 ivp = IVP([B], name="CGLE (co-moving frame)")
46 ivp.add_equation((dt(B) - (1 + 1j * w) * B - disp(B), -(1 + 1j * C3) * B * B * conj(B)))
47 sigma = Parameter("sigma", dtype=complex)
48 return linearize(ivp, around={B: sqrt(R2)}, eigenvalue=sigma)
49
50
51def growth(k: float, c1: float, c3: float, q: float) -> tuple[float, float]:
52 """Largest Re(sigma) over the groups +-k, and the largest eigen-residual."""
53 plan = discretize(cgle(k, c1, c3, q))
54 res = [eigenvalues(plan, g) for g in plan.groups if g[0].k != 0]
55 return (max(float(r.values.real.max()) for r in res), max(float(r.residuals.max()) for r in res))
56
57
58def bfn_threshold(k: float, c1: float) -> float:
59 """c3 at which the plane wave q = 0 becomes unstable at wavenumber k (parametric pencil in c3)."""
60 P = ParametricEVP(cgle(k, c1, -0.5 / c1, 0.0), [Parameter("c3", -0.5 / c1)],
61 group=lambda g: g[0].k > 0)
62 return critical(P, "c3", bracket=(-3.0 / c1, -0.2 / c1), growth="real", xtol=1e-9).value
63
64
65def eckhaus_threshold(k: float, c1: float, c3: float) -> float:
66 """q at which the plane wave becomes unstable at small perturbation wavenumber k."""
67 return so.brentq(lambda q: growth(k, c1, c3, q)[0], 1e-3, 0.99, xtol=1e-13)
68
69
70if __name__ == "__main__":
71 c1 = 2.0
72 print(f"Benjamin–Feir–Newell, c1 = {c1}: predicted c3 = -1/c1 = {-1 / c1}")
73 for k in (0.1, 0.05, 0.025):
74 c3k = bfn_threshold(k, c1)
75 print(f" k = {k:5.3f}: c3_threshold = {c3k:.12f}, 1 + c1 c3 = {1 + c1 * c3k:.3e}"
76 f" (exact 2x2: -(2 + k^2 (1 + c1^2)) / (2 c1) = {-(2 + k * k * (1 + c1 * c1)) / (2 * c1):.12f})")
77 assert abs(c3k + (2 + k * k * (1 + c1 * c1)) / (2 * c1)) < 1e-12, (k, c3k) # -> -1/c1 (Newell 1974)
78 for c3 in (-0.4, -0.6):
79 g, r = growth(0.05, c1, c3, 0.0)
80 print(f" c3 = {c3}: 1 + c1 c3 = {1 + c1 * c3:+.1f}, max Re(sigma) at k=0.05 = {g:+.3e} (residual {r:.1e})")
81 print("Eckhaus (1 + c1 c3 > 0): q_E^2 = (1 + c1 c3) / (3 + c1 c3 + 2 c3^2)")
82 for c1, c3 in ((0.0, 0.0), (0.5, 0.5), (1.0, -0.5)):
83 qE2 = (1 + c1 * c3) / (3 + c1 * c3 + 2 * c3 ** 2)
84 for k in (0.02, 0.01):
85 q = eckhaus_threshold(k, c1, c3)
86 print(f" c1={c1}, c3={c3}, k={k}: q^2 = {q * q:.10f} (theory {qE2:.10f}, gap {q * q - qE2:.2e})")
87 assert abs(q * q - qE2) < 1e-4, (c1, c3, q * q, qE2) # Eckhaus, at k = 0.01
88 print(cgle(0.05, 2.0, -0.6, 0.0).explain())