Simulating a partial differential equation - reaction-diffusion systems and Turing patterns

I want simulate a reaction-diffusion system described by a PDE called the FitzHugh–Nagumo equation.

The system that has been proposed by Alan Turing as a model of animal coat pattern formation and is exhibited by,

subject to Newman Boundary Conditions and Random Initial Conditions. I'm following the steps in this tutorial https://ipython-books.github.io/124-simulating-a-partial-differential-equation-reaction-diffusion-systems-and-turing-patterns/, but I was not successful. I thought of this code below.

L = 100; (*length of square*)
T = 50; (*Time integration*)
\[Tau] = 0.1; (*parameter*)
a = 0.00028; (*parameter*)
b = 0.005; (*parameter*)
k = -0.005; (*parameter*)

(*system of nonlinear PDE*)
pde = {D[u[t, x, y], t] ==
a*(D[u[t, x, y], x, x] + D[u[t, x, y], y, y]) + u[t, x, y] -
u[t, x, y]^3 - v[t, x, y] + k,
D[v[t, x, y], t] ==
b/\[Tau]*(D[v[t, x, y], x, x] + D[v[t, x, y], y, y]) +
1/\[Tau] u[t, x, y] - 1/\[Tau] v[t, x, y]};
(*Newman boundary condition*)
bc = {(D[u[t, x, y], x] /. x -> -L) ==
0, (D[u[t, x, y], x] /. x -> L) ==
0, (D[u[t, x, y], y] /. y -> -L) ==
0, (D[u[t, x, y], y] /. y -> L) ==
0, (D[v[t, x, y], x] /. x -> -L) ==
0, (D[v[t, x, y], x] /. x -> L) ==
0, (D[v[t, x, y], y] /. y -> -L) ==
0, (D[v[t, x, y], y] /. y -> L) == 0};
(*initial condition*)
ic = {u[0, x, y] == RandomReal[], v[0, x, y] == RandomReal[]};
eqns = Flatten@{pde, bc, ic};

sol = NDSolve[eqns, {u, v}, {t, 0, T}, {x, -L, L}, {y, -L, L}];

Table[DensityPlot[{u[t, x, y]/.sol,v[tx,y]/.sol}, {x, -L, L}, {y, -L, L},
AspectRatio -> Automatic
, PlotRange -> All
, ColorFunction -> "SunsetColors"
, {t, 0, T, 0.1}]

Can someone help me? Thanks in advance.

You're really close. Just a couple things: 1) your initial conditions are spatially uniform random numbers, which will prevent pattern formation, and 2) according to the linked post, the domain should be from -1 to 1, not -100 to 100 (those are 100 grid points they used).

Here's a version that addresses those two issues, plus uses Method -> {"IDA", "ImplicitSolver" -> {"GMRES"}} to greatly speed up NDSolve.

pts = 100;
L = 1;(*length of square*)
T = 8;(*Time integration*)
τ = 0.1;(*parameter*)
a = 0.00028;(*parameter*)
b = 0.005;(*parameter*)
k = -0.005;(*parameter*)

(*system of nonlinear PDE*)
pde = {D[u[t, x, y], t] == a*(D[u[t, x, y], x, x] + D[u[t, x, y], y, y])
+ u[t, x, y] - u[t, x, y]^3 - v[t, x, y] + k,
D[v[t, x, y], t] == b/τ*(D[v[t, x, y], x, x] + D[v[t, x, y], y, y])
+ 1/τ u[t, x, y] - 1/τ v[t, x, y]};

(*Newman boundary condition*)

bc = {(D[u[t, x, y], x] /. x -> -L) == 0,
(D[u[t, x, y], x] /. x -> L) == 0,
(D[u[t, x, y], y] /. y -> -L) == 0,
(D[u[t, x, y], y] /. y -> L) == 0,
(D[v[t, x, y], x] /. x -> -L) == 0,
(D[v[t, x, y], x] /. x -> L) == 0,
(D[v[t, x, y], y] /. y -> -L) == 0,
(D[v[t, x, y], y] /. y -> L) == 0};

(*initial condition*)

ic = {u[0, x, y] == Interpolation[Flatten[
Table[{x, y, RandomReal[]}, {x, -L, L, 2/pts}, {y, -L, L, 2/pts}], 1]][x, y],
v[0, x, y] == Interpolation[Flatten[
Table[{x, y, RandomReal[]}, {x, -L, L, 2/pts}, {y, -L, L, 2/pts}], 1]][x, y]};

eqns = Flatten@{pde, bc, ic};

sol = NDSolve[eqns, {u, v}, {t, 0, T}, {x, -L, L}, {y, -L, L},
Method -> {"MethodOfLines", "SpatialDiscretization" ->
{"TensorProductGrid", "MinPoints" -> pts, "MaxPoints" -> pts},
Method -> {"IDA", "ImplicitSolver" -> {"GMRES"}}}];

GraphicsGrid[Table[
time = t0 + 3*t1;
DensityPlot[u[time, x, y] /. sol, {x, -L, L}, {y, -L, L}, ColorFunction -> "SunsetColors",
PlotLabel -> "t=" <> ToString[time], Ticks -> False]
, {t1, 0, 2}, {t0, 0, 2}], ImageSize -> 600] As an alternative, maybe you could discretize the PDE to ODEs yourself, either manually or using @xzczd's functions here.

I'm curious about your "ecology" tag. Got another application in mind?

• Ah, making Turing patterns was one of many things I wanted to do, but never found the time for... thanks for this! – J. M. will be back soon Oct 6 '18 at 16:28
• I did some work on this here: physicsforums.com/threads/… – Dominic Apr 7 at 10:26
• @J.M.willbebacksoon thanks for the bonus :) – Chris K Sep 19 at 18:05
• Consider it a very belated show of thanks. ;) – J. M. will be back soon Sep 21 at 5:18