0
$\begingroup$

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

enter image description here

The boundary conditions are:

enter image description here

And the initial conditions are:

enter image description here

or by scheme

enter image description here

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?

$\endgroup$
  • $\begingroup$ 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! $\endgroup$ – Chris K Dec 4 '17 at 23:02
3
$\begingroup$

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:

enter image description here

N2:

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.