Turing patterns

I am new in learning Turing patterns. Is there any sample code available to generate such patterns in ecology model (Lotka–Volterra model)?  The above figure is taken from this paper, and is based on the following equations:  • There is a Wolfram Demonstrations Project demo here – ciao Jun 2 at 9:02
• Since this site is about Mathematica, the software, it is appreciated if questions include the theory (in this case, the differential equations) so that answers can focus on the coding. – C. E. Jun 2 at 9:13
• Can someone guide me to some texts also where several types of patterns are discussed in detail? For example, what are cold spots, hot spots, labyrinthine pattern etc. – Sankha Jun 2 at 14:15
• Have you seen this – user21 Jun 3 at 5:34
• @Sankha This site is for questions related to Mathematica. Your question (in the comment, not the main post) is a physics/biology question that has nothing to do with Mathematica. It would be off-topic here. – Szabolcs Jun 3 at 8:03

I developed a reaction-diffusion-advection model of pattern formation in semi-arid vegetation (tiger bush) 20 years ago, which shows a type of Turing instability. Plants ($$n$$) consume water ($$w$$) and facilitate each other by increasing water infiltration ($$wn^2$$ term). The model is set on a hillside so water advects downhill at speed $$v$$ and plants disperse as a diffusion term. $${\partial n \over \partial t}=wn^2-mn+\left({\partial^2 \over \partial x^2}+{\partial^2 \over \partial y^2}\right)n$$ $${\partial w \over \partial t}=a-w-wn^2+v{\partial w \over \partial x}$$

Here's a Mathematica implementation using NDSolve's MethodOfLines.

a = 0.3; (* nondimensional rainfall *)
m = 0.1; (* nondimensional plant mortality *)
v = 182.5; (* nondimensional water speed *)

tmax = 1000; (* max time *)
l = 200; (* nondimensional size of domain *)
pts = 40; (* numerical spatial resolution *)

(* random initial condition for plants *)
n0 = Interpolation[Flatten[Table[
{x, y, RandomReal[{0.99, 1.01}]}, {x, 0, l, l/pts}, {y, 0, l, l/pts}]
, 1], InterpolationOrder -> 0];

(* solve it *)
sol = NDSolve[{
D[n[x, y, t], t] == w[x, y, t] n[x, y, t]^2 - m n[x, y, t]
+ D[n[x, y, t], {x, 2}] + D[n[x, y, t], {y, 2}],
D[w[x, y, t], t] == a - w[x, y, t] - w[x, y, t] n[x, y, t]^2
- v D[w[x, y, t], x],
(* initial conditions *)
n[x, y, 0] == n0[x, y], w[x, y, 0] == a,
(* periodic boundary conditions *)
n[0, y, t] == n[l, y, t], w[0, y, t] == w[l, y, t],
n[x, 0, t] == n[x, l, t], w[x, 0, t] == w[x, l, t]
}, {w, n}, {t, 0, tmax}, {x, 0, l}, {y, 0, l},
Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> pts, "MaxPoints" -> pts}}
][];

(* look at final distribution *)
DensityPlot[Evaluate[n[x, y, tmax] /. sol], {x, 0, l}, {y, 0, l},
FrameLabel -> {"x", "y"}, PlotPoints -> pts,
ColorFunctionScaling -> False] Animated: Reference:

• Thanks for editing the question significantly and for the nice code. I will be greatful if you can share some texts dealing with several types of patterns and their physical meaning. – Sankha Jun 2 at 14:18
• I think those edits were mostly @Szabolcs's. One classic text to look at is Mathematical Biology II: Spatial Models and Biomedical Applications by James D. Murray. – Chris K Jun 2 at 14:26
• Ok, then thanks to both of you. – Sankha Jun 2 at 15:36
• Great answer! Your answer here is also relevant. – C. E. Jun 2 at 16:06
• Very lovely pictures! – J. M. will be back soon Jul 28 at 2:37

I did some work with the Brusselator some time ago. This is the reaction-diffusion equations which generate Turing patterns. There are some things you need to know:
(1) The non-linear PDEs have periodic boundary conditions. That means when you solve the system over a grid and you get to the end on the right side, the next point is on the left side. Same for the top and bottom. This is equivalent to solving the system over a torus.
(2) There was at the time some problems solving the system using NDSolve. Perhaps that has been resolved.
(3) The Laplacian in the system is sensitive to step size and is due to what I recall is von Neumann stability. Therefore, the step size is usually taken to be unity.

Below is a simple example not using NDSolve for these reasons and computing the Laplacian manually. And here is a reference for some of the work:

n = 64;
a = 4.5;
b = 7.5;
du = 2;
dv = 16;
dt = 0.01;
totaliter = 10000;
u = a + 0.3 RandomReal[{-0.5, 0.5}, {n, n}];
v = b/a + 0.3 RandomReal[{-0.5, 0.5}, {n, n}];

cf = Compile[{{uIn, _Real, 2}, {vIn, _Real,
2}, {aIn, _Real}, {bIn, _Real}, {duIn, _Real},
{dvIn, _Real},{dtIn, _Real}, {iterationsIn,
_Integer}},
Block[{u = uIn, v = vIn, lap, dt = dtIn, k = bIn +
1,kern = N[{{0, 1, 0}, {1, -4, 1}, {0, 1, 0}}], du =
duIn,
dv = dvIn},
Do[lap =
RotateLeft[u, {1, 0}] + RotateLeft[u, {0, 1}] +
RotateRight[u, {1, 0}] + RotateRight[u, {0, 1}] -
4*u;
u = u + dt (du lap + a - u (k - v u));
lap =
RotateLeft[v, {1, 0}] + RotateLeft[v, {0, 1}] +
RotateRight[v, {1, 0}] + RotateRight[v, {0, 1}] -
4*v;
v = v + dt (dv lap + u (b - v u));
, {iterationsIn}];
{u, v}]];

Timing[c1 = cf[u, v, a, b, du, dv, dt,
totaliter];]

ListDensityPlot[c1[]] Update: Wanted to update the recommendation below by Halirutan regarding global variables. Doing this reduced the execution time by 1/2. And also wanted to be more thorough and post the classical Turing patterns of stripes (b=7.5) and spots (b=7.0):

cf2 = With[{a = a, b = b},
Compile[{{uIn, _Real, 2}, {vIn, _Real,
2}, {aIn, _Real}, {bIn, _Real}, {duIn, _Real}, {dvIn, _Real}, \
{dtIn, _Real}, {iterationsIn, _Integer}},
Block[{u = uIn, v = vIn, lap, dt = dtIn, k = bIn + 1,
kern = N[{{0, 1, 0}, {1, -4, 1}, {0, 1, 0}}], du = duIn,
dv = dvIn},
Do[lap =
RotateLeft[u, {1, 0}] + RotateLeft[u, {0, 1}] +
RotateRight[u, {1, 0}] + RotateRight[u, {0, 1}] - 4*u;
u = u + dt (du lap + a - u (k - v u));
lap =
RotateLeft[v, {1, 0}] + RotateLeft[v, {0, 1}] +
RotateRight[v, {1, 0}] + RotateRight[v, {0, 1}] - 4*v;
v = v + dt (dv lap + u (b - v u));, {iterationsIn}];
{u, v}]]]; • Nice, +1. You need to be careful when using global variables like a and b inside compiled. They will trigger that the compiled function calls the main kernel again if it needs the values for those variables. You can wrap With[{a=a,b=b},...] around your definition. It makes a difference of 18 seconds vs. 25 seconds when using n=256. – halirutan Jun 3 at 20:23
• Thanks. I updated the code above. It cut execution time by 1/2. – Dominic Jun 4 at 10:46