12
$\begingroup$

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,

enter image description here

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.

$\endgroup$
17
+200
$\begingroup$

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]

Mathematica graphics

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?

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.