# Solving a Reaction-diffusion system

I'm just a beginner in Mathematica software. I have version 12.

I am looking for a program that will solve a prey-predator system with diffusion. My problem is that I tried a previously posted program, but the code is not compatible with Mathematica 12. Another problem is that I need to plot the solution as an animation to show the periodic motion. The program is next:

 pts = 100;
t0 = 0;
t1 = 2;
tmax = 3000;
(*length of square*)  L = 1;(*Time integration*) T = 8;(* Diffusion \
parameter for the prey*)d1 = 0.00028;(* Diffusion parameter for the \
predator*)d2 = 0.00028;(* Fertility parameter for the prey*)a = 1;(* \
Mortality parameter of the prey in the presence of predator*)b = 1;(* \
Fertility parameter of the predator*) c = 1;(* Fertility parameter of \
the predator in the presence of the prey*) d = 1;
(* Mortality parameter of the predator*) e = 1;

(*system of nonlinear PDE*)
pde = {D[u[t, x, y], t] ==
d1 (D[u[t, x, y], x, x] + D[u[t, x, y], y, y]) + a u[t, x, y] -
b u[t, x, y]*w[t, x, y],
D[w[t, x, y], t] ==
d2 (D[w[t, x, y], x, x] + D[w[t, x, y], y, y]) + c w[t, x, y] +
d u[t, x, y]*w[t, x, y] - e w[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[w[t, x, y], x] /. x -> -L) ==
0, (D[w[t, x, y], x] /. x -> L) ==
0, (D[w[t, x, y], y] /. y -> -L) ==
0, (D[w[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],
w[0, x, y] ==
Interpolation[
Flatten[Table[{x, y, RandomReal[]}, {x, -L, L, 2/pts}, {y, -L,
L, 2/pts}], 1]][x, y]};
(*ic={u[0,x,y]\[Equal]0,w[0,x,y]\[Equal]0};*)

eqns = Flatten@{pde, bc, ic};
sol = NDSolve[eqns, {u, w}, {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 + 2*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]
(Animate[DensityPlot[
Evaluate[u[x, y, time] /. sol], {x, 0, L}, {y, 0, L},
FrameLabel -> {"x", "y"}, PlotPoints -> pts,
ColorFunctionScaling -> False], {time, 0, 3000}]


It appears the following warnings: NDSolve: Warning: boundary and initial conditions are inconsistent NDSolve: Warning: estimated initial error on the specified spatial grid in the direction of independent variable x exceeds prescribed error tolerance NDSolve: "Further output of

 \!$$\*StyleBox[RowBox[{\"InterpolatingFunction\", \ \"::\", \"dmval\"}], \"MessageName\"]$$ will be suppressed during \


this calculation."

Also the plot anime doesn't out.

If someone will help me, I will be very grateful.

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, [by clicking the checkmark sign](tinyurl.com/4srwe26 – Dunlop Jan 7 '20 at 19:46
• Have you seen this? – Henrik Schumacher Jan 7 '20 at 19:47
• Perhaps to get more assistance you can show the code (plus the link to the code) that you used. Furthermore there are quite some examples of similar codes that you can find if you search for Reaction-Diffusion equations – Dunlop Jan 7 '20 at 19:48
• Does this answer your question? Simulating a partial differential equation - reaction-diffusion systems and Turing patterns – Dunlop Jan 7 '20 at 19:50
• Kamal, the example form (178559) does run in version 12. It just throws warnings because several parameters (for example the boundary conditions and the time step size) have been chosen suboptimally. – Henrik Schumacher Jan 8 '20 at 9:37

## 2 Answers

In version 12, the code can be simplified. Firstly, one should normalize the equations by d1 to cover all stages of the process by replacing t->d1 t. Then we have {t,0,1} in the new variables, which corresponds to {t,0,1/d1} in the original equations. Secondly, change pts to pts=10 in the initial data. Thirdly, we use NeumannValue[] in the boundary conditions, then we get

pts = 10;
tmax = 1;
(*length of square*)L = 1;(*Time integration*)T = 1;(*Diffusion \
parameter for the prey*)d1 = 0.00028;(*Diffusion parameter for the \
predator*)d2 = 0.00028;(*Fertility parameter for the prey*)a = \
1;(*Mortality parameter of the prey in the presence of predator*)b = \
1;(*Fertility parameter of the predator*)c = 1;(*Fertility parameter \
of the predator in the presence of the prey*)d = 1;
(*Mortality parameter of the predator*)e = 1;

(*system of nonlinear PDE*)
pde = {D[u[t, x, y],
t] - (d1 (D[u[t, x, y], x, x] + D[u[t, x, y], y, y]) +
a u[t, x, y] - b u[t, x, y]*w[t, x, y])/d1,
D[w[t, x, y],
t] - (d2 (D[w[t, x, y], x, x] + D[w[t, x, y], y, y]) +
c w[t, x, y] + d u[t, x, y]*w[t, x, y] - e w[t, x, y])/d1};

u0 = Interpolation[
Flatten[Table[{x, y, RandomReal[]}, {x, -L, L, 2/pts}, {y, -L, L,
2/pts}], 1]]; w0 =
Interpolation[
Flatten[Table[{x, y, RandomReal[]}, {x, -L, L, 2/pts}, {y, -L, L,
2/pts}], 1]];
reg = Rectangle[{-L, -L}, {L, L}];

ic = {u[0, x, y] == u0[x, y], w[0, x, y] == w0[x, y]};
(*Newman boundary condition*)
bc = NeumannValue[0, True];
eqns = {pde == {bc, bc}, ic};
Monitor[sol =
NDSolve[eqns, {u, w}, {t, 0, T}, {x, y} \[Element] reg,
EvaluationMonitor :> (monitor = Row[{"t=", t}])], monitor];


Solution visualization

Table[
DensityPlot[
Evaluate[u[t, x, y] /. First[sol]], {x, -L, L}, {y, -L, L},
ColorFunction -> Hue, PlotLabel -> Row[{"t=", t}], Frame -> False,
PlotRange -> All], {t, 0.002, .01, .002 }]

frames = Table[
DensityPlot[
Evaluate[u[t, x, y] /. First[sol]], {x, -L, 0}, {y, -L, 0},
ColorFunction -> "SunsetColors", PlotLabel -> t, Frame -> False,
PlotRange -> All], {t, 0, tmax, .02 tmax}];

ListAnimate[frames]  • Thank you very much for the answer, this is really what I was really asking. Thank you for the effort again. Kind regards, – Kamal Khalil Feb 28 '20 at 15:16
• @KamalKhalil You are welcome! – Alex Trounev Feb 28 '20 at 15:43

The Helmholtz equation is a reaction-diffusion equation. There are many examples in the documentation.

Maybe something like this:

region = RegionDifference[Rectangle[{0, 0}, {5, 10}], Disk[{5, 5}, 3]];
op = -Laplacian[u[x, y], {x, y}] - u[x, y];
bcs = {DirichletCondition[u[x, y] == 0, x == 0 && 8 <= y <= 10],
DirichletCondition[u[x, y] == 100, (x - 5)^2 + (y - 5)^2 == 3^2]};
ufun = NDSolveValue[{op == 0, bcs}, u, Element[{x, y}, region]];
ContourPlot[ufun[x, y], Element[{x, y}, region],
ColorFunction -> "Temperature", AspectRatio -> Automatic] For more information have a look at Solving PDE with FEM

You will receive better answers if your question is more specific. What have you tried (code) and what was the output (error messages).

• @KamalKhalil, can you add that information to your post. There is an edit button and then you should add have you have tried so far. – user21 Jan 8 '20 at 15:39
• user21. Thank you very much. I posted my program in the edited version. – Kamal Khalil Jan 8 '20 at 15:50