Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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, 
       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

  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

  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?

share|improve this question

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.