Gray-Scott Parameters

Module 4: Space and Pattern Formation

The feed rate $f$ and kill rate $k$ decide what the Gray-Scott solver produces. Pearson (1993) mapped this $(f, k)$ plane with forward Euler on a periodic $256 \times 256$ grid of physical size $2.5 \times 2.5$, time step $1$, and diffusion coefficients $d_1 = 2 \times 10^{-5}$, $d_2 = 10^{-5}$. He started from $(1, 0)$ with a central $20 \times 20$ square set to $(1/2, 1/4)$ plus 1% noise. The course values $d_1 = 0.1$, $d_2 = 0.05$ in grid units ($dx = 1$) keep his ratio $d_1/d_2 = 2$ but are not the same scaling, so compare patterns qualitatively.

Pearson (1993) also reports three facts about the kinetics that you should keep in mind:

The trivial state $(1, 0)$ is linearly stable for all positive $f$ and $k$.
The patterns arise from finite-amplitude perturbations, not from a small-amplitude instability.
With $d_1/d_2 = 2$, Turing patterns occur only in a narrow region near $f = k = 0.0625$.

Five Parameter Sets

Use $d_1 = 0.1$ and $d_2 = 0.05$:

Type	$f$	$k$
A	0.040	0.060
B	0.014	0.047
C	0.062	0.065
D	0.078	0.061
E	0.082	0.059

Figure 1: $v(x, y)$ at $t = 6000$ for the five types, on a periodic $96 \times 96$ grid, $dx = 1$, $dt = 2$, starting from one noisy square. On this small grid and time window, Type B decays back to $(1, 0)$.

Run each set longer and on the full $250 \times 250$ grid of the assignment. Which sets give isolated spots, which give stripes or labyrinths, and which fill the domain with a moving front?

Connection Back to Turing Analysis

The steady-state, Jacobian, and mode-matrix pipeline from Turing Numerics still applies. Only the spectrum changes: with periodic boundaries the wavenumbers are $2\pi n/L$.

The trivial state $(u, v) = (1, 0)$ has Jacobian $\operatorname{diag}(-f, -(f+k))$, so it is linearly stable for every mode. The patterns in Figure 1 therefore start from the finite-amplitude square, not from a Turing instability of $(1, 0)$ (Pearson 1993). Most of Pearson’s patterns occur where no nontrivial uniform state exists at all.

The nontrivial homogeneous equilibria satisfy $u_* v_* = f + k$. Substituting into $u_* v_*^2 = f(1 - u_*)$ gives

\[ u_* = \frac{1}{2}\left(1 \pm \sqrt{1 - \frac{4(f+k)^2}{f}}\right), \qquad v_* = \frac{f+k}{u_*}, \tag{1}\]

which exist when $4(f+k)^2 \le f$. The Jacobian at any homogeneous equilibrium is

\[ J_{\mathrm{GS}} = \begin{pmatrix} -v_*^2 - f & -2u_*v_* \\ v_*^2 & 2u_*v_* - (f+k) \end{pmatrix}. \tag{2}\]

Exercise

Implement both helpers. Keep only equilibria with $u_* > 0$.

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 (u_star, v_star)
    return equilibria


def gray_scott_jacobian(u_star, v_star, f=0.040, k=0.060):
    # TODO: return the 2x2 Jacobian array
    pass

Solution

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)],
        ]
    )

Exercise

For Type A ($f = 0.040$, $k = 0.060$), do nontrivial equilibria exist? What about Type C?

Solution

Type A: $4(f+k)^2/f = 4 \cdot 0.01 / 0.04 = 1$, so the discriminant is $0$ and the two nontrivial equilibria merge at $u_* = 0.5$, $v_* = 0.2$, a saddle-node point. In floating point the discriminant comes out as about $-2 \times 10^{-16}$, so the code above returns only $(1, 0)$. Type C: $4 \cdot 0.127^2 / 0.062 \approx 1.04 > 1$, so only $(1, 0)$ exists. In both cases the pattern cannot come from a Turing instability of a nontrivial uniform state that is stable without diffusion.

Once both helpers work, reuse find_unstable_modes_2d from Unstable Modes in 2D. Swap in gray_scott_jacobian, use $D = \operatorname{diag}(d_1, d_2)$ with no $\gamma$, and replace $\pi/L$ by $2\pi/L$ for the periodic spectrum.

Exercise

Type D ($f = 0.078$, $k = 0.061$) has nontrivial equilibria. List them, compute $\operatorname{tr} J$ and $\det J$ at each one, and scan the periodic modes on $L = 250$ with $d_1 = 0.1$, $d_2 = 0.05$. Is any of them Turing unstable?

def find_unstable_modes_gs(u_star, v_star, f, k, d1=0.1, d2=0.05, length=250.0, num_modes=60):
    jac = gray_scott_jacobian(u_star, v_star, f, k)
    diffusion = np.diag([d1, d2])
    unstable_modes = []
    # TODO: loop over (n, m), including (0, 0), with k2 = (2 pi / L)**2 (n**2 + m**2)
    # TODO: store ((n, m), growth_rate) when the growth rate is positive
    return unstable_modes

Solution

There are three equilibria. $(1, 0)$ is stable, with $\det J > 0$ and $\operatorname{tr} J < 0$, and no mode grows. The upper state $u_* \approx 0.548$, $v_* \approx 0.254$ has $\det J < 0$: it is a saddle, unstable already without diffusion, so the flat mode $(0, 0)$ grows too. The lower state $u_* \approx 0.452$, $v_* \approx 0.307$ is stable without diffusion and no mode grows for $d_1/d_2 = 2$. So none of them shows a Turing instability, in line with Pearson (1993).

def find_unstable_modes_gs(u_star, v_star, f, k, d1=0.1, d2=0.05, length=250.0, num_modes=60):
    jac = gray_scott_jacobian(u_star, v_star, f, k)
    diffusion = np.diag([d1, d2])
    unstable_modes = []
    for n in range(num_modes):
        for m in range(num_modes):
            k2 = (2 * np.pi / length) ** 2 * (n**2 + m**2)
            growth_rate = np.max(np.linalg.eigvals(jac - k2 * diffusion).real)
            if growth_rate > 0:
                unstable_modes.append(((n, m), growth_rate))
    unstable_modes.sort(key=lambda item: item[1], reverse=True)
    return unstable_modes


f, k = 0.078, 0.061
for u_star, v_star in gray_scott_equilibria(f, k):
    jac = gray_scott_jacobian(u_star, v_star, f, k)
    modes = find_unstable_modes_gs(u_star, v_star, f, k)
    print(f"({u_star:.3f}, {v_star:.3f})  tr = {np.trace(jac):.4f}  det = {np.linalg.det(jac):.5f}  unstable modes: {len(modes)}")

References

Pearson, John E. 1993. “Complex Patterns in a Simple System.” Science 261 (5118): 189–92. https://doi.org/10.1126/science.261.5118.189.

--- title: "Gray-Scott Parameters" subtitle: "Module 4: Space and Pattern Formation" format: html --- The feed rate $f$ and kill rate $k$ decide what the Gray-Scott solver produces. @pearson1993complex mapped this $(f, k)$ plane with forward Euler on a periodic $256 \times 256$ grid of physical size $2.5 \times 2.5$, time step $1$, and diffusion coefficients $d_1 = 2 \times 10^{-5}$, $d_2 = 10^{-5}$. He started from $(1, 0)$ with a central $20 \times 20$ square set to $(1/2, 1/4)$ plus 1% noise. The course values $d_1 = 0.1$, $d_2 = 0.05$ in grid units ($dx = 1$) keep his ratio $d_1/d_2 = 2$ but are not the same scaling, so compare patterns qualitatively. @pearson1993complex also reports three facts about the kinetics that you should keep in mind: - The trivial state $(1, 0)$ is linearly stable for all positive $f$ and $k$. - The patterns arise from finite-amplitude perturbations, not from a small-amplitude instability. - With $d_1/d_2 = 2$, Turing patterns occur only in a narrow region near $f = k = 0.0625$. ## Five Parameter Sets Use $d_1 = 0.1$ and $d_2 = 0.05$: | Type | $f$ | $k$ | |---|---|---| | A | 0.040 | 0.060 | | B | 0.014 | 0.047 | | C | 0.062 | 0.065 | | D | 0.078 | 0.061 | | E | 0.082 | 0.059 | ```{python} #| label: fig-gsp-types #| fig-cap: '$v(x, y)$ at $t = 6000$ for the five types, on a periodic $96 \times 96$ grid, $dx = 1$, $dt = 2$, starting from one noisy square. On this small grid and time window, Type B decays back to $(1, 0)$.' #| fig-width: 10 #| fig-height: 2.4 #| echo: false import matplotlib.pyplot as plt import numpy as np from amlab.pdes_2d.gray_scott import gray_scott_pde sets = {"A": (0.040, 0.060), "B": (0.014, 0.047), "C": (0.062, 0.065), "D": (0.078, 0.061), "E": (0.082, 0.059)} fig, axs = plt.subplots(1, 5, figsize=(10, 2.4)) for ax, (label, (f, k)) in zip(axs, sets.items()): n, r = 96, 8 rng = np.random.default_rng(0) uv = np.ones((2, n, n)) uv[1] = 0.0 c = n // 2 uv[:, c - r : c + r, c - r : c + r] = 0.5 uv[:, c - r : c + r, c - r : c + r] *= 1 + 0.1 * rng.standard_normal((2, 2 * r, 2 * r)) for _ in range(3000): uv = uv + 2.0 * gray_scott_pde(0, uv, d1=0.1, d2=0.05, f=f, k=k) ax.imshow(uv[1], cmap="viridis", vmin=0, vmax=0.5) ax.set_title(f"{label}: f={f}, k={k}", fontsize=8) ax.axis("off") plt.tight_layout() plt.show() plt.close() ``` Run each set longer and on the full $250 \times 250$ grid of the assignment. Which sets give isolated spots, which give stripes or labyrinths, and which fill the domain with a moving front? ## Connection Back to Turing Analysis The steady-state, Jacobian, and mode-matrix pipeline from [Turing Numerics](turing-numerics.qmd) still applies. Only the spectrum changes: with periodic boundaries the wavenumbers are $2\pi n/L$. The trivial state $(u, v) = (1, 0)$ has Jacobian $\operatorname{diag}(-f, -(f+k))$, so it is linearly stable for every mode. The patterns in @fig-gsp-types therefore start from the finite-amplitude square, not from a Turing instability of $(1, 0)$ [@pearson1993complex]. Most of Pearson's patterns occur where no nontrivial uniform state exists at all. The nontrivial homogeneous equilibria satisfy $u_* v_* = f + k$. Substituting into $u_* v_*^2 = f(1 - u_*)$ gives $$ u_* = \frac{1}{2}\left(1 \pm \sqrt{1 - \frac{4(f+k)^2}{f}}\right), \qquad v_* = \frac{f+k}{u_*}, $$ {#eq-gsp-steady} which exist when $4(f+k)^2 \le f$. The 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-gsp-jacobian} ::: {.callout-note} ## Exercise Implement both helpers. Keep only equilibria with $u_* > 0$. ```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 (u_star, v_star) 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"} ## Solution ```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)], ] ) ``` ::: ::: {.callout-note} ## Exercise For Type A ($f = 0.040$, $k = 0.060$), do nontrivial equilibria exist? What about Type C? ::: ::: {.callout-tip collapse="true"} ## Solution Type A: $4(f+k)^2/f = 4 \cdot 0.01 / 0.04 = 1$, so the discriminant is $0$ and the two nontrivial equilibria merge at $u_* = 0.5$, $v_* = 0.2$, a saddle-node point. In floating point the discriminant comes out as about $-2 \times 10^{-16}$, so the code above returns only $(1, 0)$. Type C: $4 \cdot 0.127^2 / 0.062 \approx 1.04 > 1$, so only $(1, 0)$ exists. In both cases the pattern cannot come from a Turing instability of a nontrivial uniform state that is stable without diffusion. ::: Once both helpers work, reuse `find_unstable_modes_2d` from [Unstable Modes in 2D](gierer-meinhardt-2d-modes.qmd). Swap in `gray_scott_jacobian`, use $D = \operatorname{diag}(d_1, d_2)$ with no $\gamma$, and replace $\pi/L$ by $2\pi/L$ for the periodic spectrum. ::: {.callout-note} ## Exercise Type D ($f = 0.078$, $k = 0.061$) has nontrivial equilibria. List them, compute $\operatorname{tr} J$ and $\det J$ at each one, and scan the periodic modes on $L = 250$ with $d_1 = 0.1$, $d_2 = 0.05$. Is any of them Turing unstable? ```python def find_unstable_modes_gs(u_star, v_star, f, k, d1=0.1, d2=0.05, length=250.0, num_modes=60): jac = gray_scott_jacobian(u_star, v_star, f, k) diffusion = np.diag([d1, d2]) unstable_modes = [] # TODO: loop over (n, m), including (0, 0), with k2 = (2 pi / L)**2 (n**2 + m**2) # TODO: store ((n, m), growth_rate) when the growth rate is positive return unstable_modes ``` ::: ::: {.callout-tip collapse="true"} ## Solution There are three equilibria. $(1, 0)$ is stable, with $\det J > 0$ and $\operatorname{tr} J < 0$, and no mode grows. The upper state $u_* \approx 0.548$, $v_* \approx 0.254$ has $\det J < 0$: it is a saddle, unstable already without diffusion, so the flat mode $(0, 0)$ grows too. The lower state $u_* \approx 0.452$, $v_* \approx 0.307$ is stable without diffusion and no mode grows for $d_1/d_2 = 2$. So none of them shows a Turing instability, in line with @pearson1993complex. ```python def find_unstable_modes_gs(u_star, v_star, f, k, d1=0.1, d2=0.05, length=250.0, num_modes=60): jac = gray_scott_jacobian(u_star, v_star, f, k) diffusion = np.diag([d1, d2]) unstable_modes = [] for n in range(num_modes): for m in range(num_modes): k2 = (2 * np.pi / length) ** 2 * (n**2 + m**2) growth_rate = np.max(np.linalg.eigvals(jac - k2 * diffusion).real) if growth_rate > 0: unstable_modes.append(((n, m), growth_rate)) unstable_modes.sort(key=lambda item: item[1], reverse=True) return unstable_modes f, k = 0.078, 0.061 for u_star, v_star in gray_scott_equilibria(f, k): jac = gray_scott_jacobian(u_star, v_star, f, k) modes = find_unstable_modes_gs(u_star, v_star, f, k) print(f"({u_star:.3f}, {v_star:.3f}) tr = {np.trace(jac):.4f} det = {np.linalg.det(jac):.5f} unstable modes: {len(modes)}") ``` :::