Solving a Reaction-diffusion system

Question

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.

Answer 1

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]

Answer 2

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