Differential Equations - Gray Scott Equations

Introduction to Gray Scott Equations

The Gray-Scott equations are a set of partial differential equations used to model reaction-diffusion systems, particularly in chemistry and biology. They describe the interaction between two chemical species, typically referred to as and , that diffuse and react with each other over time and space.

The Gray-Scott model is known for producing a variety of complex, spatially structured patterns, such as spots, stripes, and waves. The nature of the patterns depends on the parameters , , , and . The model can exhibit both steady-state and oscillatory behavior.

Formulation of Gray-Scott Equations

The Gray-Scott model is defined by the following system of PDEs:

where:

• and are the concentrations of the two chemical species.
• and are the diffusion coefficients of and , respectively.
• is the feed rate of the reactant .
• is the rate at which decays.
• is the Laplacian operator, representing diffusion in space.

Interpretation of Terms

• Diffusion Terms: and represent the spatial diffusion of species and . These terms cause the chemicals to spread out over time.
• Reaction Terms: and represent the interaction between and . The term in the equation for indicates that is consumed in the reaction to produce , while the term in the equation for indicates that is produced in the reaction.
• Feed and Decay Terms: represents the addition of species into the system at a constant rate, while represents the decay of species .

Derivation of Numerical Solutions

We try to keep everything simple. The Gray-Scott equations are discretized in space using a finite difference method, and the resulting system of ODEs is solved using explicit Euler time-stepping.

To give the formulation, we use to denote the time step size, and to denote the spatial step sizes. For the time derivative, we approximate it using the forward difference:

The Laplacian in 2D cartesian coordinates is approximated using the five-point stencil:

is updated in a similar way. New values of and are calculated using the following update rules:

denote at the next time step. The update rules are applied to all grid points in the domain. After updating all grid points, we set and repeat the process for the next time step.

Here we formally present the algorithm for solving the Gray-Scott equations:

1. Use a rectangular grid to discretize the spatial-temporal domain.
2. Use forward difference to approximate the time derivative, and five-point central difference to approximate the Laplacian to obtain at each grid point. The term form the truncation error of the finite difference method. The same applies to . The equation applies for time steps 。
3. Drop the truncation error terms and rearrange the equation to obtain the update rule .
4. The initial condition (IC) is set by assigning given values to all at . The periodic boundary condition (BC) is set by copying the values of at the opposite boundary, i.e., and the same for . The boundary condition is applied at each time step .

The update rule can be efficiently implemented on fragment shaders, which are widely supported by modern GPUs. We use a straightforward implementation in GLSL to demonstrate the algorithm.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

#version 330

in vec2 uv0;
uniform sampler2D texture0;
out vec4 fragColor;

uniform float c_du;
uniform float c_dv;
uniform float c_f;
uniform float c_k;

uniform float dx;
uniform float dy;
uniform float hx;
uniform float hy;
uniform float dt;

void main() {
    vec2 center = texture(texture0, uv0).xy;
    vec2 left = texture(texture0, uv0 + vec2(-dx, 0)).xy;
    vec2 right = texture(texture0, uv0 + vec2(dx, 0)).xy;
    vec2 bottom = texture(texture0, uv0 + vec2(0, -dy)).xy;
    vec2 top = texture(texture0, uv0 + vec2(0, dy)).xy;

    float lu = (left.x + right.x - 2.0 * center.x) / (hx * hx)
        + (top.x + bottom.x - 2.0 * center.x) / (hy * hy);
    float lv = (left.y + right.y - 2.0 * center.y) / (hx * hx)
        + (top.y + bottom.y - 2.0 * center.y) / (hy * hy);
    float du = c_du * lu - center.x * center.y * center.y + c_f * (1.0 - center.x);
    float dv = c_dv * lv + center.x * center.y * center.y - (c_f + c_k) * center.y;
    vec2 result = center + vec2(du, dv) * dt;

    fragColor = vec4(result, 0.0, 1.0);
}

The periodic boundary condition is not explicitly implemented in the shader. Instead, we use the GL_REPEAT texture wrapping mode to handle the boundary condition.

Test of Accuracy and Stability

Accuracy

It is difficult to obtain an analytical solution for the Gray-Scott equations. From the truncation error term, we could expect the method to have a first-order accuracy in time and second-order accuracy in space.

Convergence

We could empirically show that the method converges under the condition

First we test the convergence of the method by increasing the spatial resolution while keeping the time step size constant. We set , , and . We use a fixed time step size and change the spatial resolution. The result indicate that the system fail to start evolving when the spatial resolution is too low. However, when the spatial resolution is too high, the system is numerically unstable and diverges in a few time steps. It is only possible to obtain stable results when the spatial resolution is within a certain range. Note that since different random noise is added to the initial condition, the exact pattern of the solution may vary in each simulation.

Simulation result at T=20000 with different spatial resolutions. The solution diverges when N > 291.

We also tried to change the time step size while keeping the spatial resolution constant. The result shows that lower time step size leads to consistent results, and the system diverges when the time step size exceeds a certain value.

Simulation result at T=20000 and N=256 with different t'. The solution diverges when t' > 1.17

We finally demonstrate that by decreasing the time step size and increasing the spatial resolution simultaneously, the method could be conditionally convergent. Following the condition , we set to

in these simulations. The result shows that the method is convergent under this condition.

Simulation result when the convergence condition is satified.

Stability

We test the stability of the method by extending final simulation time while keeping the spatial resolution and time step size constant. The result shows that the system is stable and the solution converges to a steady state.

Simulation result at different time steps

Simulation Results

We wrote a GPU-accelerated simulation program using OpenGL and GLSL. The program solves the Gray-Scott equations on a 2D grid and visualizes the results in real-time.

Screenshot of the simulation program

The program contains a GUI that allows users to adjust the parameters of the Gray-Scott equations, such as , , , and $k. Users can observe the patterns generated by the model in real-time and explore the dynamics of reaction-diffusion systems. To run the program, simply run python grayscott.py in the terminal.

We present the results of the simulation with different parameters. The patterns generated by the Gray-Scott equations are highly sensitive to the parameters. By adjusting these parameters, we can observe a wide range of patterns, including spots, stripes, and waves. The following images show the phase diagram of the Gray-Scott equation and the results obtained at different locations in the domain. We set , , , , and .

The phase diagram of the Gray-Scott equation (Pearson 1993)

Results obtained at different conditions

Conclusion

The Gray-Scott equations are a powerful tool for modeling reaction-diffusion systems and generating complex spatial patterns. By discretizing the equations and solving them numerically, we can observe a wide range of patterns, including spots, stripes, and waves. The numerical method we presented is conditionally convergent and stable under certain conditions. By adjusting the parameters of the model, we can explore the rich dynamics of reaction-diffusion systems and gain insights into the underlying mechanisms that give rise to complex patterns in nature.

许可协议

本文采用 署名-相同方式共享 4.0 国际 许可协议，转载请注明出处。

分享文章