System of nonlinear PDE 2D (Reaction-Diffusion type) with periodic boundary condition

Question

I want to solve a system of Pde (2D) reaction diffusion type using NDSolve

whose boundary conditions are

and the initial conditions are

or

I thought of the following code

(*parameters*)

L = 5;
T = 10;

(*system of nonlinear PDE*)

pde = {D[N1[t, x, y], t] == 
    D[N1[t, x, y], x, x] + 
     D[N1[t, x, y], y, 
      y] + (1 - N1[t, x, y] - 0.5 N2[t, x, y]) N1[t, x, y], 
   D[N2[t, x, y], t] == 
    D[N2[t, x, y], x, x] + 
     D[N2[t, x, y], y, 
      y] + (1 - N2[t, x, y] - 0.5 N1[t, x, y]) N2[t, x, y]};

(*periodic boundary condition*)

bc = {N1[t, -L, y] == N1[t, L, y], N1[t, x, -L] == N1[t, x, L], 
   N2[t, -L, y] == N2[t, L, y], N2[t, x, -L] == N2[t, x, L]};

(*initial condition*)

ic = {N1[0, x, y] == If[-0.5 <= x <= 1 && -0.5 <= y <= 1, 1, 0], 
   N2[0, x, y] == If[-0.5 <= x <= 1 && -0.5 <= y <= 1, 0, 1]};
eqns = Flatten@{pde, bc, ic};

{N1, N2} = 
 NDSolve[eqns, {N1, N2}, {t, 0, T}, {x, -L, L}, {y, -L, L}, 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid"}}]

However, the following errors appear.

Also, when constructing the plot (by DensityPlot), I verify that there is a failure in the initial condition, as shows

Can someone help me?

Answer 1

As @xzczd alluded to, your code is almost correct. The message NDSolve::mxsst is more of a note than an error or even a warning.

Here's a working version because I love reaction-diffusion equations (this appears to be a diffusive Lotka-Volterra competition model, parameterized so that the species can stably coexist). I used a diffusion coefficient d=0.002, to see the traveling wave better, increased the number of grid points, and fixed your initial conditions to match the description.

(*parameters*)
L = 5;
T = 100;
d = 0.002;

(*system of nonlinear PDE*)

pde = {
D[N1[t, x, y], t] == d*(D[N1[t, x, y], x, x] + D[N1[t, x, y], y, y])
+ (1 - N1[t, x, y] - 0.5 N2[t, x, y]) N1[t, x, y], 
D[N2[t, x, y], t] == d*(D[N2[t, x, y], x, x] + D[N2[t, x, y], y, y])
+ (1 - N2[t, x, y] - 0.5 N1[t, x, y]) N2[t, x, y]
};

(*periodic boundary condition*)

bc = {N1[t, -L, y] == N1[t, L, y], N1[t, x, -L] == N1[t, x, L], 
N2[t, -L, y] == N2[t, L, y], N2[t, x, -L] == N2[t, x, L]};

(*initial condition*)

ic = {N1[0, x, y] == If[-0.5 <= x <= 0.5 && -0.5 <= y <= 0.5, 1, 0], 
N2[0, x, y] == If[-0.5 <= x <= 0.5 && -0.5 <= y <= 0.5, 0, 1]};

eqns = Flatten@{pde, bc, ic};

sol = NDSolve[eqns, {N1, N2}, {t, 0, T}, {x, -L, L}, {y, -L, L}, 
Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 200}}][[1]];

Here's a movie: