How solve a PDE system with Specific Initial Condition?

Question

I'm trying 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 boundary condition

The boundary conditions are:

And the initial conditions are:

or by scheme

The code is written below

(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[-4.2 <= x <= -4.7 && -4.2 <= y <= -4.7 && 
     4.2 <= x <= 4.7 && -4.2 <= y <= -4.7 && -4.2 <= x <= -4.7 && 
      4.2 <= y <= 4.7 && -0.5 <= x <= 0.5 && -0.5 <= y <= 0.5 && 
     4.2 <= x <= 4.7 && 4.2 <= y <= 4.7, 1, 0], 
  N2[0, x, y] == 
   If[-4.2 <= x <= -4.7 && -4.2 <= y <= -4.7 && 
     4.2 <= x <= 4.7 && -4.2 <= y <= -4.7 && -4.2 <= x <= -4.7 && 
      4.2 <= y <= 4.7 && -0.5 <= x <= 0.5 && -0.5 <= y <= 0.5 && 
     4.2 <= x <= 4.7 && 4.2 <= y <= 4.7, 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"}}]

I believe the problem arises from the way the initial condition was implemented. Besides trying to implement these conditions using If, I tried to use Piecewise, but I was not successful. Can anybody help me?

Answer 1

I don't have time right now to fully solve the problem, but there is an issue with the way you are describing the conditionals. You need to flip some signs and use some Ors instead of Ands

The new system is given by:

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[(-4.2 >= x >= -4.7 && -4.2 >= y >= -4.7) || (4.2 <= x <= 
        4.7 && -4.2 >= y >= -4.7) || (-4.2 >= x >= -4.7 && 
       4.2 <= y <= 4.7) || (-0.5 <= x <= 0.5 && -0.5 <= y <= 
        0.5) || (4.2 <= x <= 4.7 && 4.2 <= y <= 4.7), 1, 0], 
  N2[0, x, y] == 
   If[(-4.2 >= x >= -4.7 && -4.2 >= y >= -4.7) || (4.2 <= x <= 
        4.7 && -4.2 >= y >= -4.7) || (-4.2 >= x >= -4.7 && 
       4.2 <= y <= 4.7) || (-0.5 <= x <= 0.5 && -0.5 <= y <= 
        0.5) || (4.2 <= x <= 4.7 && 4.2 <= y <= 4.7), 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"}}]

Right now, NDSolve is throwing a fit, but, at first glance, discontinuities in the ICs are likely the culprits.

Let me know if you want me to push this any further.

Update

You seem to get a pretty good approximation of the initial conditions by upping the grid points in the spatial discretization. Changing the MaxPoints to 200 there isn't too big of a hit with respect to computation time, but visible gains in accuracy. So, finding the solution to the PDE

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

Plotting:

n1plots = 
  Table[Plot3D[N1[t, x, y] /. sol[[1]], {x, -L, L}, {y, -L, L}, 
    AspectRatio -> Automatic, 
    PlotRange -> {Automatic, Automatic, {0, 1}}], {t, 
    Subdivide[0, T, 25]}];

n2plots = 
  Table[Plot3D[N2[t, x, y] /. sol[[1]], {x, -L, L}, {y, -L, L}, 
    AspectRatio -> Automatic, 
    PlotRange -> {Automatic, Automatic, {0, 1}}], {t, 
    Subdivide[0, T, 25]}];

N1:

N2: