Add interactive simulator links and expand unstable modes discussion in 2D Gierer-Meinhardt and Gray-Scott models

Author: Daniel Precioso, PhD

modules/pde-1d/turing-instability.qmd

@@ -10,6 +10,12 @@ Patterns do not appear by accident. Turing instability explains when a uniform s
 
 The previous page gave the intuition: local activation can reinforce a bump, inhibition can suppress nearby growth, and diffusion can either smooth or amplify spatial differences depending on the parameter regime. This page turns that picture into a linear stability test and asks a precise question, which parameter values and which spatial modes make the flat steady state lose stability?
 
+::: {.callout-note}
+## Interactive Simulator
+
+You can follow the theory on this page while running the [interactive 1D Gierer-Meinhardt simulator](https://bam-applied-math-lab.streamlit.app/Gierer_Meinhardt_1D). The simulator lets you pick a boundary condition, adjust domain length and reaction parameters, build a custom perturbation, and check which modes the simulation selects against the prediction from linear theory.
+:::
+
 ```{python}
 #| label: fig-pde-turing-space
 #| fig-cap: 'Turing space for the Gierer-Meinhardt model in the $(a, d)$ plane.'
@@ -236,7 +242,7 @@ plt.show()
 
 Mark the point $(a, d) = (0.40, 20)$ on your plot. Is it inside the instability region?
 
-This is the same Turing-space diagram that appears in the left panel of the interactive Gierer-Meinhardt code, so the algebra and the simulation refer to the same parameter plane.
+This is the same Turing-space diagram that appears in the [interactive 1D Gierer-Meinhardt simulator](https://bam-applied-math-lab.streamlit.app/Gierer_Meinhardt_1D), so the algebra and the simulation refer to the same parameter plane.
 
 ## Compute $F_n$ and the Unstable Mode List in 1D
 
@@ -294,154 +300,6 @@ def find_unstable_modes_1d(a=0.40, b=1.00, d=30.0, length=40.0, num_modes=10):
 
 The first entry in this list is the dominant linear mode, which predicts the earliest visible wavelength in the simulation.
 
-## Extend the Scan to 2D
-
-In two dimensions, the modes are pairs $(n_x, n_y)$. On a rectangle with side lengths $L_x$ and $L_y$ and Neumann boundaries,
-
-$$
-\lambda_{n_x,n_y} = -\left(\frac{n_x\pi}{L_x}\right)^2 - \left(\frac{n_y\pi}{L_y}\right)^2.
-$$ {#eq-pde-turing-2d-modes}
-
-The only real change is that you now scan a grid of mode pairs and build
-
-$$
-F_{n_x,n_y} = J + \lambda_{n_x,n_y}\,\mathrm{diag}(1, d).
-$$
-
-```python
-import numpy as np
-
-def find_unstable_modes_2d(a=0.40, b=1.00, d=30.0, length_x=20.0, length_y=50.0, num_modes=8):
-    fu, fv, gu, gv = gierer_meinhardt_jacobian(a, b)
-    jac = np.array([[fu, fv], [gu, gv]])
-    diffusion = np.diag([1.0, d])
-    unstable_modes = []
-
-    for nx in range(num_modes):
-        for ny in range(num_modes):
-            if nx == 0 and ny == 0:
-                continue
-
-            # TODO: compute lambda_nm and F_nm
-            # TODO: store mode pairs with positive growth rate
-
-    return unstable_modes
-```
-
-::: {.callout-tip collapse="true"}
-## Hint: 2D mode scan (click to expand)
-
-```python
-def find_unstable_modes_2d(a=0.40, b=1.00, d=30.0, length_x=20.0, length_y=50.0, num_modes=8):
-    fu, fv, gu, gv = gierer_meinhardt_jacobian(a, b)
-    jac = np.array([[fu, fv], [gu, gv]])
-    diffusion = np.diag([1.0, d])
-    unstable_modes = []
-
-    for nx in range(num_modes):
-        for ny in range(num_modes):
-            if nx == 0 and ny == 0:
-                continue
-
-            lambda_nm = -((nx * np.pi / length_x) ** 2 + (ny * np.pi / length_y) ** 2)
-            f_nm = jac + lambda_nm * diffusion
-            growth_rate = np.max(np.linalg.eigvals(f_nm).real)
-            if growth_rate > 0:
-                unstable_modes.append(((nx, ny), growth_rate))
-
-    unstable_modes.sort(key=lambda item: item[1], reverse=True)
-    return unstable_modes
-```
-:::
-
-For periodic boundaries, the same loop works after replacing each factor of $\pi / L$ by $2\pi / L$.
-
-## Apply the Same Workflow to Gray-Scott
-
-The theory on this page does not depend on the specific reaction law. Gray-Scott uses the same equilibrium to Jacobian to mode-matrix pipeline [@gray1984autocatalytic; @pearson1993complex].
-
-Write the reaction terms as
-
-$$
-R_u(u, v) = -u v^2 + f(1-u), \qquad R_v(u, v) = u v^2 - (f+k)v.
-$$
-
-One homogeneous equilibrium is always
-
-$$
-(u_*, v_*) = (1, 0).
-$$
-
-If $v_* \neq 0$, then $u_* v_* = f + k$ and the nontrivial equilibria satisfy
-
-$$
-u_* = \frac{1}{2}\left(1 \pm \sqrt{1 - \frac{4(f+k)^2}{f}}\right), \qquad v_* = \frac{f+k}{u_*},
-$$ {#eq-pde-turing-gs-steady}
-
-whenever $1 - 4(f+k)^2/f \ge 0$.
-
-The Gray-Scott Jacobian at a homogeneous equilibrium is
-
-$$
-J_{\mathrm{GS}} =
-\begin{pmatrix}
--v_*^2 - f & -2u_*v_* \\
-v_*^2 & 2u_*v_* - (f+k)
-\end{pmatrix}.
-$$ {#eq-pde-turing-gs-jacobian}
-
-At the trivial state $(1, 0)$, this Jacobian is diagonal and every linear mode decays. So the classical Turing test is only interesting at a nontrivial homogeneous equilibrium when one exists.
-
-For periodic boundaries on a rectangle,
-
-$$
-\lambda_{m,n} = -\left(\frac{2\pi m}{L_x}\right)^2 - \left(\frac{2\pi n}{L_y}\right)^2,
-$$
-
-and the mode matrix is
-
-$$
-F^{\mathrm{GS}}_{m,n} = J_{\mathrm{GS}} + \lambda_{m,n}\,\mathrm{diag}(d_1, d_2).
-$$
-
-```python
-import numpy as np
-
-def gray_scott_equilibria(f=0.040, k=0.060):
-    equilibria = [(1.0, 0.0)]
-    discriminant = 1 - 4 * (f + k) ** 2 / f
-
-    if discriminant >= 0:
-        root = np.sqrt(discriminant)
-        for sign in (+1, -1):
-            u_star = 0.5 * (1 + sign * root)
-            if u_star > 0:
-                equilibria.append((u_star, (f + k) / u_star))
-
-    return equilibria
-
-
-def gray_scott_jacobian(u_star, v_star, f=0.040, k=0.060):
-    return np.array(
-        [
-            [-v_star**2 - f, -2 * u_star * v_star],
-            [v_star**2, 2 * u_star * v_star - (f + k)],
-        ]
-    )
-```
-
-That is the whole bridge from Gierer-Meinhardt to Gray-Scott. Keep the eigenvalue scan, replace the Jacobian, replace the boundary spectrum, then sort the unstable modes and read off the dominant one.
-
-## Compare Two Cases
-
-Use your functions to compare $d=20$ and $d=30$ with $a=0.40$ and $b=1.00$.
-
-1. Does each case satisfy the Turing conditions?
-2. How many unstable modes appear in 1D?
-3. What is the leading mode in 2D?
-
-These answers will help you interpret the 1D assignment now and the 2D simulations in the next session.
-
 ## What's Next?
 
 Now that you know how to predict pattern-forming parameter regimes, close the 1D loop by completing the assignment. After that, carry the same ideas to a 2D grid in the next session.
@@ -450,5 +308,5 @@ Now that you know how to predict pattern-forming parameter regimes, close the 1D
 [PDEs in 2D](../pde-2d/index.qmd){.btn .btn-primary}
 
 ::: {.callout-tip}
-Reference code is available in [amlab/pdes/gierer_meinhardt_1d.py](https://github.com/daniprec/BAM-Applied-Math-Lab/blob/main/amlab/pdes/gierer_meinhardt_1d.py) and [amlab/pdes/gierer_meinhardt_2d.py](https://github.com/daniprec/BAM-Applied-Math-Lab/blob/main/amlab/pdes/gierer_meinhardt_2d.py).
+Reference code is available in [amlab/pdes/gierer_meinhardt_1d.py](https://github.com/daniprec/BAM-Applied-Math-Lab/blob/main/amlab/pdes/gierer_meinhardt_1d.py).
 :::

modules/pde-2d/gierer-meinhardt-2d.qmd

@@ -326,9 +326,79 @@ def update_frame(_):
 
 The 5-point stencil is the standard choice, but you can improve the rotational symmetry of the Laplacian with a 9-point stencil. That is a good extension if you want to compare numerical accuracy.
 
-## Explore
+## Unstable Modes on a Rectangle
 
-Use the Turing page to compare the cases $d=20$ and $d=30$. Does the pattern that forms on your grid match the change suggested by the unstable mode analysis?
+The Turing instability page showed how to list unstable modes in 1D by scanning the Neumann eigenvalues $\lambda_n = -(n\pi/L)^2$. In two dimensions, the modes are pairs $(n_x, n_y)$. On a rectangle with side lengths $L_x$ and $L_y$ and Neumann boundaries,
+
+$$
+\lambda_{n_x,n_y} = -\left(\frac{n_x\pi}{L_x}\right)^2 - \left(\frac{n_y\pi}{L_y}\right)^2.
+$$ {#eq-pde-turing-2d-modes}
+
+The mode matrix is the same as before:
+
+$$
+F_{n_x,n_y} = J + \lambda_{n_x,n_y}\,\mathrm{diag}(1, d).
+$$
+
+Extend your 1D mode scan to two dimensions.
+
+```python
+import numpy as np
+
+def find_unstable_modes_2d(a=0.40, b=1.00, d=30.0, length_x=20.0, length_y=50.0, num_modes=8):
+    fu, fv, gu, gv = gierer_meinhardt_jacobian(a, b)
+    jac = np.array([[fu, fv], [gu, gv]])
+    diffusion = np.diag([1.0, d])
+    unstable_modes = []
+
+    for nx in range(num_modes):
+        for ny in range(num_modes):
+            if nx == 0 and ny == 0:
+                continue
+
+            # TODO: compute lambda_nm and F_nm
+            # TODO: store mode pairs with positive growth rate
+
+    return unstable_modes
+```
+
+::: {.callout-tip collapse="true"}
+## Hint: 2D mode scan (click to expand)
+
+```python
+def find_unstable_modes_2d(a=0.40, b=1.00, d=30.0, length_x=20.0, length_y=50.0, num_modes=8):
+    fu, fv, gu, gv = gierer_meinhardt_jacobian(a, b)
+    jac = np.array([[fu, fv], [gu, gv]])
+    diffusion = np.diag([1.0, d])
+    unstable_modes = []
+
+    for nx in range(num_modes):
+        for ny in range(num_modes):
+            if nx == 0 and ny == 0:
+                continue
+
+            lambda_nm = -((nx * np.pi / length_x) ** 2 + (ny * np.pi / length_y) ** 2)
+            f_nm = jac + lambda_nm * diffusion
+            growth_rate = np.max(np.linalg.eigvals(f_nm).real)
+            if growth_rate > 0:
+                unstable_modes.append(((nx, ny), growth_rate))
+
+    unstable_modes.sort(key=lambda item: item[1], reverse=True)
+    return unstable_modes
+```
+:::
+
+For periodic boundaries, replace each factor of $\pi / L$ by $2\pi / L$.
+
+## Compare $d = 20$ and $d = 30$
+
+Use the functions above to compare $d = 20$ and $d = 30$ with $a = 0.40$ and $b = 1.00$.
+
+1. Does each case satisfy the Turing conditions?
+2. How many unstable mode pairs $(n_x, n_y)$ appear on your rectangle?
+3. What is the dominant mode pair?
+
+Then run the simulation for each case and check whether the observed stripe or spot pattern matches the dominant mode prediction.
 
 ## What's Next?
 

modules/pde-2d/gray-scott.qmd

@@ -210,7 +210,76 @@ The same steady-state to Jacobian to mode-matrix pipeline from the Turing page s
 - For Neumann boundaries, the admissible wavenumbers use factors of $\pi / L$.
 - For periodic boundaries, the admissible wavenumbers use factors of $2\pi / L$.
 
-For the usual base state $(u, v) = (1, 0)$, the linearized system is damped, so the classical Turing calculation only becomes interesting at a nontrivial homogeneous equilibrium when one exists. In that regime, the mode matrices are built exactly as before, but with the periodic Laplacian eigenvalues. That is the clean conceptual bridge from Gierer-Meinhardt to Gray-Scott.
+For the usual base state $(u, v) = (1, 0)$, the linearized system is damped, so the classical Turing calculation only becomes interesting at a nontrivial homogeneous equilibrium when one exists.
+
+The nontrivial homogeneous equilibria satisfy $u_* v_* = f + k$, and one can show that
+
+$$
+u_* = \frac{1}{2}\left(1 \pm \sqrt{1 - \frac{4(f+k)^2}{f}}\right), \qquad v_* = \frac{f+k}{u_*},
+$$ {#eq-pde-gs-steady}
+
+whenever $1 - 4(f+k)^2/f \ge 0$. The Gray-Scott Jacobian at any homogeneous equilibrium is
+
+$$
+J_{\mathrm{GS}} =
+\begin{pmatrix}
+-v_*^2 - f & -2u_*v_* \\
+v_*^2 & 2u_*v_* - (f+k)
+\end{pmatrix}.
+$$ {#eq-pde-gs-jacobian}
+
+Implement both helpers and run the mode-matrix eigenvalue scan with the periodic spectrum to find unstable modes for a given $(f, k, d_1, d_2)$.
+
+```python
+import numpy as np
+
+def gray_scott_equilibria(f=0.040, k=0.060):
+    equilibria = [(1.0, 0.0)]
+    discriminant = 1 - 4 * (f + k) ** 2 / f
+
+    if discriminant >= 0:
+        root = np.sqrt(discriminant)
+        for sign in (+1, -1):
+            u_star = 0.5 * (1 + sign * root)
+            # TODO: append the nontrivial equilibrium if u_star > 0
+
+    return equilibria
+
+
+def gray_scott_jacobian(u_star, v_star, f=0.040, k=0.060):
+    # TODO: return the 2x2 Jacobian array
+    pass
+```
+
+::: {.callout-tip collapse="true"}
+## Hint: Gray-Scott equilibria and Jacobian (click to expand)
+
+```python
+def gray_scott_equilibria(f=0.040, k=0.060):
+    equilibria = [(1.0, 0.0)]
+    discriminant = 1 - 4 * (f + k) ** 2 / f
+
+    if discriminant >= 0:
+        root = np.sqrt(discriminant)
+        for sign in (+1, -1):
+            u_star = 0.5 * (1 + sign * root)
+            if u_star > 0:
+                equilibria.append((u_star, (f + k) / u_star))
+
+    return equilibria
+
+
+def gray_scott_jacobian(u_star, v_star, f=0.040, k=0.060):
+    return np.array(
+        [
+            [-v_star**2 - f, -2 * u_star * v_star],
+            [v_star**2, 2 * u_star * v_star - (f + k)],
+        ]
+    )
+```
+:::
+
+Once both helpers work, reuse `find_unstable_modes_2d` from the Gierer-Meinhardt 2D page, swap in `gray_scott_jacobian`, and replace $\pi / L$ with $2\pi / L$ to get the periodic mode spectrum. That is the direct conceptual bridge from Gierer-Meinhardt to Gray-Scott [@gray1984autocatalytic; @pearson1993complex].
 
 ## Optional: 9-Point Stencil