# How solve a PDE system with Specific Initial Condition?

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?

• 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! Dec 4, 2017 at 23:02

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[], {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[], {x, -L, L}, {y, -L, L},
AspectRatio -> Automatic,
PlotRange -> {Automatic, Automatic, {0, 1}}], {t,
Subdivide[0, T, 25]}];


N1: N2: • When you are available I would like you to help me solve this problem.
– SAC
Dec 5, 2017 at 12:24
• @SAC I've updated the answer with an improved solution and some plotting. Dec 7, 2017 at 17:55