I would like to solve a PDE system reaction-diffusion type (2D spatial + 1 temporal) coupled as described below. Another question of this same system was solved here:

System of nonlinear PDE 2D (Reaction-Diffusion type) with periodic (in laterals) and null flux (top and bottom) boundary condition

Also, initial conditions is assumed with a noise.

The system is shown below:

enter image description here

The code:

L = 5;(*length of square*)
pts = 100;
T = 300;(*Time integration*)
α = 1.5;
k = 0.1;
ρ = 18.5;
a = 92;
b = 64;
d = 10;
γ = 9;
(*system of nonlinear PDE*)
pde = {D[u[t, x, y], 
     t] == (D[u[t, x, y], x, x] + 
       D[u[t, x, y], y, y]) + γ (α - 
        u[t, x, y] - (ρ u[t, x, y] v[t, x, y])/(
        1 + u[t, x, y] + k u[t, x, y] u[t, x, y])), 
   D[v[t, x, y], t] == 
    d (D[v[t, x, y], x, x] + 
        D[v[t, x, y], y, 
         y]) + γ (α (b - v[t, x, y]) - (ρ u [t, x, 
          y] v[t, x, y])/(1 + u[t, x, y] + k u[t, x, y] u[t, x, y]))};
(*Newman boundary condition*)
bc = {(D[u[t, x, y], x] /. x -> 0) == 
    0, (D[u[t, x, y], x] /. x -> 2 L) == 0, 
   u[t, x, 0] == u[t, x, 2 L] , (D[v[t, x, y], x] /. x -> 0) == 
    0, (D[v[t, x, y], x] /. x -> 2 L) == 0, 
   v[t, x, 0] == v[t, x, 2 L]};
(*initial condition*)
ic = {u[0, x, y] == 
      Flatten[Table[{x, y, RandomReal[{0.01, 1}]}, {x, 0, 2 L, 
         2/pts}, {y, 0, 2 L, 2/pts}], 1]][x, y], 
   v[0, x, y] == 
      Flatten[Table[{x, y, RandomReal[{0.01, 1}]}, {x, 0, 2 L, 
         2/pts}, {y, 0, 2 L, 2/pts}], 1]][x, y]};
eqns = Flatten@{pde, bc, ic};

sol = NDSolve[eqns, {u, v}, {t, 0, T}, {x, 0, 2 L}, {y, 0, 2 L}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> pts, "MaxPoints" -> pts}}];

But, there are errors:

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent. NDSolve::eerri: Warning: estimated initial error on the specified spatial grid in the direction of independent variable x exceeds prescribed error tolerance. NDSolve::eerri: Warning: estimated initial error on the specified spatial grid in the direction of independent variable y exceeds prescribed error tolerance.

Can anybody help me?

  • $\begingroup$ its clear from the picture u[t,x,y] depend on t, but from where you getting this x & y ? Also the gradient of U with respect to what !? $\endgroup$
    – nufaie
    Dec 9, 2022 at 14:06
  • 1
    $\begingroup$ Please notice, as indicated by the message, these are warnings i.e. they only indicate that something may be wrong. eerri warning is often a benign warning. (Just search in this site. ) ibcinc is indeed a problem in your case, for more info check this: mathematica.stackexchange.com/a/127411/1871 $\endgroup$
    – xzczd
    Dec 9, 2022 at 14:10
  • $\begingroup$ @xzczd I think the ibcinc warning is also benign in this case, because any inconsistencies will be immediately smoothed out by diffusion at t==0. $\endgroup$
    – Chris K
    Dec 9, 2022 at 18:56
  • 1
    $\begingroup$ I think your code is basically working, but may be quite slow for such high resolution and long time. The initial conditions might also be problematic for a different reason -- I recommend starting near the spatially homogeneous equilibrium given by NSolve[{0 == \[Gamma] (\[Alpha] - u - (\[Rho] u v)/(1 + u + k u u)), 0 == (\[Alpha] (b - v) - (\[Rho] u v)/(1 + u + k u u))}, {u, v}] which is {u -> 0.00128753, v -> 63.0009}. $\endgroup$
    – Chris K
    Dec 9, 2022 at 18:57
  • $\begingroup$ @ChrisK If the zero Neumann b.c. isn't important, the warning will be benign, otherwise we need to adjust the ScaleFactor sub-option a bit. (See the discussion in the post linked in my last comment. ) $\endgroup$
    – xzczd
    Dec 10, 2022 at 1:39

1 Answer 1


I must admit, I haven't tried to debug what's wrong with your implementation (a first guess could be ensuring periodicity in your initial conditions, see below) - but here's some working FEM code to get you started:


h[\[Rho]_,k_][u_,v_]=(\[Rho] u v)/(1+u+k u^2);

Derivative[1,0,0][u][t,x,y]-\[Gamma] f[a,\[Rho],k][u[t,x,y],v[t,x,y]]-Inactive[Div][( Inactive[Grad][u[t,x,y],{x,y}]),{x,y}]==0,
Derivative[1,0,0][v][t,x,y]-\[Gamma] g[\[Alpha],b,\[Rho],k][u[t,x,y],v[t,x,y]]-Inactive[Div][(d Inactive[Grad][v[t,x,y],{x,y}]),{x,y}]==0

Initial Conditions



Boundary Conditions



EvaluationMonitor:>(currentTime=Row[{"t = ",CForm[t]}])];,currentTime];


This seems to work, but it's quite slow to converge with the parameters and initial conditions you want. FWIW, here's the output of $v$ for 10 evenly-spaced time frames between $t=(0,1]$ (Note the periodic boundary conditions are correctly respected along x):



enter image description here

Spectral Methods (?)

Seems like you only have periodic boundary conditions along the horizontal direction, but if this isn't crucial to your problem - you can consider using FFT-based spectral methods. Here's an example walk-through (click on the Wolfram buttons on the top to see the Wolfram code) with code using the Gray-Scott and Cahn-Hilliard reaction-diffusion equations.

  • 2
    $\begingroup$ 2 PeriodicBoundaryCondition aren't enough for periodicity of derivative, you need 2*2, and ghost layer or triangular mesh, see discussion here: mathematica.stackexchange.com/q/204546/1871 $\endgroup$
    – xzczd
    Dec 10, 2022 at 1:49

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.