19
$\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$
1

1 Answer 1

24
+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$
8
  • 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$ Oct 6, 2018 at 16:28
  • $\begingroup$ I did some work on this here: physicsforums.com/threads/… $\endgroup$
    – Dominic
    Apr 7, 2019 at 10:26
  • $\begingroup$ @J.M. thanks for the bonus :) $\endgroup$
    – Chris K
    Sep 19, 2019 at 18:05
  • $\begingroup$ Consider it a very belated show of thanks. ;) $\endgroup$ Sep 21, 2019 at 5:18
  • $\begingroup$ @ChrisK thank you for introducing IDA, I had a code that wasn't even compiling but now it creates a result :D Few questions: The reason why IDA is required to be in MethodOfLines is that IDA changes the solvers when required, right? And GMRES is the name of the method to be changed? $\endgroup$
    – confused
    Mar 17, 2021 at 11:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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