The timeline below gives some idea of the development of the key ideas.
Open heat_equation.vti (click on this link) and set it running. Scroll down in the Info Pane to see the formula:
This says that the rate of change of chemical a is equal to the Laplacian of chemical a. This means that any location that is lower than its neighbors will increase, and any location that is higher than its neighbors will decrease. This equation describes how heat will spread out in a solid or how ink would diffuse through a liquid.
|
We compute the Laplacian of each chemical at each location by finite differencing with a stencil of radius 1. In 1D this means that for a string of values:
the Laplacian at B will be given by: A + C - 2B. In 2D, for a cell C surround by 4 neighbors N, S, E, W the five-point stencil gives a Laplacian of: N + S + E + W - 4C. In general, we use a graph Laplacian, which is equal to the sum of (P-centre) for each neighbor P. On meshes, the neighbors are defined to be any cell that shares a vertex with the central cell. |
Having computed the rate of change of the chemical at each location, we then apply the Euler method to estimate the new value at some point in the future. Knowing the current value a and the rate of change delta_a gives:
Here, timestep is a value that determines how accurate the resulting simulation will be - smaller values will give more accurate solutions but take longer to compute because more steps are needed to get to the same place. You can change the timestep in our heat equation by scrolling down in the Info Pane to 'timestep' and editing the value. If you set it to much larger than 0.2 then the values quickly diverge, due to numerical instabilities. In general it is important to check that any witnessed behaviour is genuine and not an artefact of a timestep that is too large.
There are alternatives to the Euler method that allow for larger timesteps (e.g. Runge-Kutta) but we haven't implemented them yet.
|
Have a look at heat_equation_interpolation.vti. This is the same heat equation as we saw above, but with an extra chemical b. The heat equation applies to chemical a but only in those places where the pattern was originally zero. This has the effect of 'locking-in' certain values and forcing the other values to reach an equilibrium with them. |
Now open Turing1952/spots.vti. This system has two chemicals (a and b) which are shown side by side. The formula is more complicated:
delta_a = k * (16.0f - a * b) + D_a * laplacian_a;
delta_b = max( -b, k * (a * b - b - 12.0f) + D_b * laplacian_b );
Here k, D_a and D_b are parameters, and can be edited in the Info Pane to give different behavior.
Set the system running. After a while, spots appear.
To see the same system in one dimension, scroll down in the Info Pane to "Dimensions" and edit to set it to 64 x 1 x 1. The line graph shows chemical a in a brighter color. Soon a pattern of peaks and troughs appears, in the same way as the 2D system. In 3D (try 32 x 32 x 32) we get blobs appearing.
Turing suggested that this behavior of forming a regular pattern of spots is the sort of thing that might be used in nature by plants and animals wanting to direct their own growth (morphogenesis) or for the patterns on the skins of leopards, for example. While plausible, no evidence of this was found. I think it is fair to say that for many years the theory was poorly regarded amongst biologists, especially after a different mechanism was found to drive certain well-known examples of morphogenesis. It was only in 2006 that the first evidence for Turing's reaction-diffusion in biology was found, in the even placement of hair follicles on the skin of mice. Other systems have since been found.
Also look at Turing1952/spots_noisy.vti. Here we've used a third chemical c to add noise to the system. This shows how the patterns are robust to such interference.
After the concept of reaction-diffusion was identified, many other systems have been discovered that fit into the general scheme and show interesting phenomena. A partial list is shown below.