Implementation, Turing Space, and Mode Calculations
The key change in the PDE session is spatial structure. The state is no longer a short vector of variables. It is a field sampled on a grid.
u(x, t)
The unknown is now a field over space, not a short ODE state vector.
Diffusion enters through the Laplacian, which links each location to its neighbors.
Turn the theory they already saw into concrete calculations for Gierer-Meinhardt and the later Gray-Scott transition.
Implement Neumann, Dirichlet, and periodic conditions on Cartesian grids.
Write the centered finite-difference stencil with vectorized shifts.
Advance the full field with a direct time-stepping loop.
Compute the equilibrium, the Jacobian, the mode matrix, and the unstable-mode list.
The same implementation logic works in 1D and on rectangles.
Zero flux means copy the nearest interior value.
u[0] = u[1], u[-1] = u[-2]
Fix the value at the boundary after every Euler step.
u[0] = u_left, u[-1] = u_right
Wrap one edge onto the other. With np.roll(), the stencil already does it.
On a rectangle, the zero-flux pattern is
\[ \Delta u(x) \approx \frac{u(x+\Delta x) - 2u(x) + u(x-\Delta x)}{\Delta x^2} \]
\[ \Delta u \approx u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4u_{i,j} \]
The same np.roll() idea extends later to rectangular grids.
Once the right-hand side is available, the PDE update is just
\[ U^{n+1} = U^n + \Delta t\,G(U^n). \]
du_dt, dv_dt = gierer_meinhardt_pde(0, uv, dx=dx)
uv[0] = uv[0] + dt * du_dt
uv[1] = uv[1] + dt * dv_dt
uv[:, 0] = uv[:, 1]
uv[:, -1] = uv[:, -2]The order matters. First apply Euler, then reimpose the boundary data.
\[ \frac{\partial u}{\partial t} = \Delta u + \gamma \left(a - b u + \frac{u^2}{v}\right), \qquad \frac{\partial v}{\partial t} = d\,\Delta v + \gamma (u^2 - v) \]
This is the model we use to turn the Turing theory into concrete calculations.
Set the reaction terms to zero:
\[ f(u, v) = a - bu + \frac{u^2}{v} = 0, \qquad g(u, v) = u^2 - v = 0. \]
From \(g(u_*, v_*) = 0\),
\[ v_* = u_*^2. \]
Substitute into \(f = 0\):
\[ a - bu_* + 1 = 0 \quad \Longrightarrow \quad u_* = \frac{a+1}{b}, \qquad v_* = \left(\frac{a+1}{b}\right)^2. \]
For \(a = 0.4\) and \(b = 1\),
\[ u_* = 1.4, \qquad v_* = 1.96. \]
At a general point \((u, v)\),
\[ J(u, v) = \begin{pmatrix} -b + 2u/v & -u^2/v^2 \\ 2u & -1 \end{pmatrix}. \]
At \((u_*, v_*)\) this becomes
\[ J_* = \begin{pmatrix} 2b/(a+1) - b & -\left(b/(a+1)\right)^2 \\ 2(a+1)/b & -1 \end{pmatrix}. \]
For \(a = 0.4\) and \(b = 1\),
\[ J_* = \begin{pmatrix} 0.4286 & -0.5102 \\ 2.8 & -1 \end{pmatrix}. \]
\[ f_u + g_v < 0, \qquad f_u g_v - f_v g_u > 0, \qquad g_v + d f_u > 2\sqrt{d(f_u g_v - f_v g_u)} \]
The homogeneous equilibrium is stable without spatial variation, but some spatial mode can grow.
At \(a = 0.4\), \(b = 1\), the point \(d = 30\) is inside the Turing window, while \(d = 20\) is outside.
For Neumann boundary conditions on \([0, L]\),
\[ \phi_n(x) = \cos\left(\frac{n\pi x}{L}\right), \qquad \lambda_n = -\left(\frac{n\pi}{L}\right)^2, \qquad n = 0, 1, 2, \dots \]
This is the step the professor emphasized. The spatial part is solved by hand, and those \(\lambda_n\) are what enter the temporal calculation.
For each mode \(n\),
\[ A_n = J + \lambda_n D, \qquad D = \operatorname{diag}(1, d). \]
Equivalently,
\[ A_n = J - k_n^2 D, \qquad k_n = \frac{n\pi}{L}. \]
The temporal eigenvalues of mode \(n\) are the roots of
\[ \det\left(A_n - \sigma^{(n)} I\right) = 0. \]
Positive \(\sigma^{(n)}\) means that spatial mode is unstable.
unstable_modes = []
for n in range(1, num_modes):
lambda_n = -(n * np.pi / length) ** 2
A_n = J + lambda_n * D
growth_rate = np.max(np.linalg.eigvals(A_n).real)
if growth_rate > 0:
unstable_modes.append((n, growth_rate))
unstable_modes.sort(key=lambda item: item[1], reverse=True)This list does two things.
Take \(\gamma = 1\), \(a = 0.4\), \(b = 1\).
\(d = 30\)
\(d = 20\)
The pipeline does not change.
The conceptual change is the boundary condition.
At the usual Gray-Scott base state \((u, v) = (1, 0)\), the linearized system is damped, so the classical Turing calculation becomes interesting only at a nontrivial homogeneous equilibrium when one exists.
Applied Math Lab