The GrayâScott model is a reactionâdiffusion system used to look at how simple chemical interactions can generate complex patterns. It describes the evolution of two chemical species, U and V, diffusing across a surface while undergoing an autocatalytic reaction. Unlike linear diffusion equations, which smooth out disturbances, the GrayâScott system is highly nonlinear and can produce stable spots, stripes, waves, and other very interesting Turing patterns.
Because of its non-linear nature, we must simulate the system numerically. Here, we simulate it on a 2D grid using finite difference approximations for the diffusion terms and explicit time-stepping for the reactions. Almost like a cellular automaton, but instead of having discrete states, each cell holds continuous concentration values for chemicals U and V.
The GrayâScott system models a simplified autocatalytic reaction:
U + 2V â 3V (autocatalysis)
V â P (decay)
Here, U is fed into the system at a constant rate, while V is the reactant that autocatalytically converts U into more of itself. Because of
these interactions and imbalances, small changes in the concentrations of U and V can grow and lead to complex patterns.
When diffusion is included and the system is treated as a continuous 2D medium, the concentrations
U(x,y,t) and V(x,y,t)—meaning 'the concentration of U or V at position (x,y) and time t'—evolve according to the following pair of partial differential equations:
The first term, \(\ D_\_ \nabla^2 U \ \), is what represents our diffusion, and is what causes the chemicals to spread out over time. The second term, \(\ \pm UV^2 \ \), represents the autocatalytic reaction where species U is consumed to produce more of species V. The third term in each equation represents the feed and removal rate of chemicals from the system. This maintains a non-equilibrium state, which is necessary for forming those interesting patterns.
Because the system has no closed-form solution, we approximate the PDEs to simulate the reaction-diffusion process over discrete time steps. We represent the 2D surface as a grid of cells, where each cell holds the information about the concentrations of U and V present at that location. The Laplacian operator \(\ \nabla^2\ \) is approximated using a finite difference stencil. This simulation uses a 9-point stencil for better accuracy, but other stencils (like 4-point cross or even a diagonal) can also be used.
To evolve the system in time, we step time by a small increment via the forward Euler method: