# How to solve a reaction-diffusion?

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: 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] ==
Interpolation[
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] ==
Interpolation[
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?

• 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 !? Dec 9, 2022 at 14:06
• 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 Dec 9, 2022 at 14:10
• @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. Dec 9, 2022 at 18:56
• 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}. Dec 9, 2022 at 18:57
• @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. ) Dec 10, 2022 at 1:39

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:

### Equations

h[\[Rho]_,k_][u_,v_]=(\[Rho] u v)/(1+u+k u^2);
f[a_,\[Rho]_,k_][u_,v_]=a-u-h[\[Rho],k][u,v];
g[\[Alpha]_,b_,\[Rho]_,k_][u_,v_]=\[Alpha](b-v)-h[\[Rho],k][u,v];

eqns[d_,\[Gamma]_,\[Alpha]_,\[Rho]_,k_,a_,b_]={
};


### Initial Conditions

horizontalPeriodicRandomField[{min_,max_},pts_,L_]:=Block[
{rand=RandomReal[{min,max},{pts,pts}]},
rand[[-1]]=rand[];
ListInterpolation[rand,{{-L,L},{-L,L}},PeriodicInterpolation->{True,False}]]

ics={
u[0,x,y]==horizontalPeriodicRandomField[{0.01,1},100,5][x,y],
v[0,x,y]==horizontalPeriodicRandomField[{0.01,1},100,5][x,y]
};


### Boundary Conditions

bcs={
PeriodicBoundaryCondition[u[t,x,y],x==-5,TranslationTransform[{10,0}]],
PeriodicBoundaryCondition[v[t,x,y],x==-5,TranslationTransform[{10,0}]]
};


### Solution

Monitor[{ufun,vfun}=NDSolveValue[
{eqns[10,9,3/2,37/2,1/10,92,64],bcs,ics},
{u,v},{x,y}\[Element]Rectangle[{-5,-5},{5,5}],{t,0,1},
Method->{"PDEDiscretization"->{"MethodOfLines","SpatialDiscretization"->{"FiniteElement","MeshOptions"->{"MaxCellMeasure"->0.01}}}},
EvaluationMonitor:>(currentTime=Row[{"t = ",CForm[t]}])];,currentTime];


### Visualization

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):


vFrames=With[{sol=vfun,reg=Rectangle[{-5,-5},{5,5}]},
Table[ContourPlot[sol[t,x,y],{x,y}\[Element]reg,PlotRange->All,Frame->None,Axes->None,PlotPoints->25],{t,Rest@Subdivide}]];

Rasterize[Multicolumn[vFrames,5,Appearance->"Horizontal"],ImageSize->800] ## 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 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 Dec 10, 2022 at 1:49