Solving two-component two-dimensional differential Equation with NDSolve (Brusselator Model)

Question

I'm trying to solve a two-component two-dimensional reaction-diffusion differential equation system with Mathematica. The background of the model is the so called "Brusselator Model" where one can find a nice outline online from the Institute of Theoretical Physics in Münster, Germany. The model is given by the differential equation system:

$$u_t=D_u\Delta u+a-(b+1)u+u^2 v$$ $$v_t=D_v\Delta v+b u-u^2 v$$

with $u,v$ as the solution variables depending on $x$, $y$ and $t$, $\Delta$ as the Laplace operator $D_u$ and $D_v$ as the diffusion constants and $a,b$ as positive reaction rate constants.

This system exhibits a transition from small fluctuations around the uniform steady state solution of the model into a high-amplitude stripe pattern as shown in Fig. 8.6 of the online reference, when the parameterset $D_u$, $D_v$, $a$ and $b$ is in the instability region. Based on stability analysis (see reference) the system is unstable when e.g. using the parameter set $D_u=5$, $D_v=12$, $a=3$ and $b=9$.

The steady state solution of the system is: $$u_0=a=3, v_0=\frac{b}{a}=1/3$$

In Mathematica we define a small fluctuation (e.g. $\pm2\%$) on the steady state solutions $u_0$ and $v_0$ as initial values (at $t=0$) of the differential equations

u0 = Interpolation[Flatten[Table[{{x, y}, 3 + If[x == 0 || y == 0 || x == 50 || y == 50, 
       0, 
       RandomReal[{-a/50, a/50} /. a -> 3]], {0, 0}}, {x, 0,50}, {y, 0, 50}], 1], 
       InterpolationOrder -> 2];
v0 = Interpolation[Flatten[Table[{{x, y}, 9/u0[x, y], {0, 0}}, {x, 0, 50}, {y, 0, 50}], 1], 
       InterpolationOrder -> 2];

here I restricted the initial value boundaries to the steady state without random fluctuation and enforced that the first derivate of the interpolation function is zero at the boundaries too, to keep the boundary and initial conditions in the differential equation consistent. One can plot the initial condition $u_0$ with the following statements

GraphicsRow[{
  ContourPlot[u0[x, y], {x, 0, 50}, {y, 0, 50}, 
    PlotLegends -> Automatic, 
    ColorFunction -> "SolarColors", 
    FrameLabel -> {"x", "y"}, 
    ImageSize -> Medium], 
  Plot3D[u0[x, y], {x, 0, 50}, {y, 0, 50}, PlotRange -> {2.9, 3.1}]}]

My approach to solve this system for a $x$-$y$ space interval of $50\times 50$ units and the mentioned parameters, I tried out the following statement

NDSolve[Evaluate[{
  D[u[x, y, t], t] == Du (Derivative[2, 0, 0][u][x, y, t] + Derivative[0, 2, 0][u][x, y, t]) + a - (b + 1) u[x, y, t] + v[x, y, t] (u[x, y, t])^2,
  D[v[x, y, t], t] == Dv (Derivative[2, 0, 0][v][x, y, t] + Derivative[0, 2, 0][v][x, y, t]) + b u[x, y, t] - v[x, y, t] (u[x, y, t])^2,
  u[x, y, 0] == u0[x, y],
  v[x, y, 0] == v0[x, y],
  Derivative[1, 0, 0][u][0, y, t] == 0, Derivative[1, 0, 0][u][L, y, t] == 0,
  Derivative[0, 1, 0][u][x, 0, t] == 0, Derivative[0, 1, 0][u][x, L, t] == 0,
  Derivative[1, 0, 0][v][0, y, t] == 0, Derivative[1, 0, 0][v][L, y, t] == 0,
  Derivative[0, 1, 0][v][x, 0, t] == 0, Derivative[0, 1, 0][v][x, L, t] == 0,
} /. {L -> 50, Du -> 5, Dv -> 12, a -> 3, b -> 9}], 
{u, v}, {x, 0, 50}, {y, 0, 50}, {t, 0, 0.05 4000}, MaxSteps -> {15, 15, Infinity}]

but get the following error:

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent.

which is strange since I enforced the Neumann boundary conditions through the interpolation function derivate. Finally NDSolve runs forever, yielding no result even with the rather low setting of 15 for MaxSteps in x and y direction. Any proposal for improvements?

Answer 1

I tried your excellent idea of using an Interpolation over a "noisy" point set as the initial condition to a reaction-diffusion system (Gierer-Meinhardt in one spatial dimension). It seems that one cannot really enforce the boundary condition df/dx=0 at x=0 or x=L. It would be possible with a "clamped spline" implementation which uses exactly those boundary conditions, but unfortunately this is not available in Mathematica (yet). So what I did was to make a list of {x,y} coordinates, x runs from 0 to 1 in increments of 0.01, y is some const value. Then I picked a few indices randomly between 10 and 90 and perturbed the corresponding y coordinates, then fitted an Interpolation over them. This makes sure that the "flanks" are practically constant, and there is some noise "in the middle". If I feed that as initial conditions to NDSolve, it stops complaining about inconsistent initial/boundary conditions, and produces the expected result. Of course your problem is more complex in 2 spatial dimensions, but maybe this trick helps.

Here is a function that generates such "noisy almost-constant" functions (note that it contains some ugly hard-coded stuff but the idea should be understandable I hope):

noisefunc[const_, nperturb_: 10, noisefract_: 0.1] := Module[
  {ipoints, ridx},
  ipoints = Table[{xi, const}, {xi, 0.0, 1.0, 0.01}];

  ridx = Table[RandomInteger[{10, 90}], {nperturb}];

  Do[ipoints[[ridx[[ri]]]][[2]] += 
    noisefract*RandomReal[{-1, 1}], {ri, nperturb}
   ];
  Return[Interpolation[ipoints, Method -> "Spline"]]
  ]

And here is an example usage:

nf=noisefunc[1.0]
Plot[nf[x],{x,0,1},PlotRange->All]

with the result:

Of course you have to play with the settings. If there is "too much noise", then NDSolve may tell you that the max. number of grid points is not enough etc.