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

Question

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.

Answer 1

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?