A reaction diffusion system contains 2 pde, $xc, xp$, and one ode, xnn. A variable in my coupled pde system, $xc(t,x,y)$, has a piecewice initial condition for a duration. I am trying to redo the simulations in an old article.

"$xc(t,x)$ in a rectangular region between (-20,+20) is raised to x0h for 6 seconds"

I don't have a problem with it's "piecewiseness" ("$xc(t,x,y)$ in a rectangular region between (-20,+20) is raised to x0h", which I express as xc[0,x]==If[-mix < y < +mix && -mix< x <+mix, xc0h, xc0]).

But I don't know how can I write an initial condition saying "$xc(t,x,y)$ in a square region between (-mix,+mix) is raised to xc0h for 6 seconds".

I tried to add it as an additional initial condition but it didn't worked. Do you have any idea? For a fuller expression, this is the whole code:

mx = 600;
mix = mx/10;
mt = 600; 
xc0h  = 0.9;
 xc0 = 0.130;
 xp0  = 0.095;
xnn0  = 1;

 s = NDSolve[
   D[ xc[t, x, y], t] == 
    20 (D [xc[t, x, y], t] + D[ xc[t, x, y], t]) + 
     8.1 (0.567 + (0.433 xp[t, x, y] )/(xp[t, x, y] + 4.0)) xnn[t, x, 
       y] (0.111 + (0.899 xc[t, x, y])/(0.7 + xc[t, x, y])) - (
     2.0 xc[t, x, y])/(0.1 + xc[t, x, y]) + 0.02,
   D[ xp[t, x, y], t] == 
    300 (D[ xp[t, x, y], x] + D[ xp[t, x, y], y]) - 0.17 xp[t, x, y],
   D[ xnn[t, x, y], 
     t] == (1 - xc[t, x, y]^2/(xc[t, x, y]^2 + 0.7^2)) - 
     xnn[t, x, y],
   xc[0, x, y] == If[-mix < y < +mix&&-mix < x < +mix, xc0h, xc0],
   xc[t, 0, y] == xc0,
   xc[t, mx, y] == xc0,
   xc[t, x, 0] == xc0,
   xc[t, x, mx] == xc0,
   xp[0, x, y] == xp0,
   xp[t, 0, y] == xp0,
   xp[t, mx, y] == xp0,
   xp[t, x, 0] == xp0,
   xp[t, x, mx] == xp0,
   xnn[0, x, y] == xnn0
  {xc, xp},
  {t, 0, mt},
  {x, 0, mx},
  {y, 0, mx}]
  • $\begingroup$ Gram per second looks like speed, i.e. first derivative of the function, not its value. $\endgroup$
    – yarchik
    Jan 25, 2021 at 19:04
  • $\begingroup$ @yarchik you mean D[X1[0,x],t]==If[-20<x<20,10,5] ? Then how can I indicate the duration? There should be a difference between keeping $(-20,20)$ region in this value for 6 seconds and 60 seconds right? $\endgroup$
    – confused
    Jan 25, 2021 at 19:11
  • $\begingroup$ Hmm, I do not know, I am confused my the rectangular region, which is definitely a 2D region, not 1D. $\endgroup$
    – yarchik
    Jan 25, 2021 at 19:22
  • 2
    $\begingroup$ Could you give a fuller description of the model? What's X1 and are there other variables? What's the PDE? $\endgroup$
    – Chris K
    Jan 25, 2021 at 20:49
  • 2
    $\begingroup$ Could you add a verbal description, because a few things are still not clear: 1) there is no X1 in your code, 2) the same time derivative shows up in three places in the first two equations, is this a typo? and 3) what are the other variables supposed to be doing when the one is fixed? $\endgroup$
    – Chris K
    Jan 26, 2021 at 15:46

1 Answer 1


Since the only $x$-dependence is in the initial conditions, and there is no $y$-dependence at all, we can remove all references to $x$ and $y$. The resulting equations will be valid at all $(x,y)$.

We treat the initial conditions as parameters and obtain a solution with this code:

mt = 20;
soln = ParametricNDSolve[{
    D[xc[t], t] == 20 (D[xc[t], t] + D[xc[t], t]) + 
        8.1 (0.567 + (0.433 xp[t])/(xp[t] + 4.0)) *
        xnn[t] (0.111 + (0.899 xc[t])/(0.7 + xc[t])) - 
        (2.0 xc[t])/(0.1 + xc[t]) + 0.02,
   D[xp[t], t] == 300 (D[xp[t], t] + D[xp[t], t]) - 0.17 xp[t],
   D[xnn[t], t] == (1 - xc[t]^2/(xc[t]^2 + 0.7^2)) - xnn[t],   
   xc[0] == xc0,
   xp[0] == xp0,
   xnn[0] == xnn0},
  {xc, xp}, {t, 0, mt}, {xc0, xp0, xnn0}];

We can integrate out to a greater time, but after plotting the resulting (exponential) solutions we choose to plot only the shorter interval. Plots can be produced with the following code. The first plot is valid for positions $(x,y)$ where $x \notin (0,30)$ and the second plot is valid for positions $(x,y)$ where $0 < x < 30$. The third plot is valid for all $(x,y)$.

  xc0h = 0.9,
  xc0 = 0.130,
  xp0 = 0.095,
  xnn0 = 1},
   Plot[xc[xc0, xp0, xnn0][t] /. soln,
    {t, 0, mt}, Frame -> True, 
    FrameLabel -> {"t", "\!\(\*SubscriptBox[\(X\), \(c\)]\)", 
      "Initial condition is \!\(\*SubscriptBox[\(X\), \(c0\)]\)"}],
   Plot[xc[xc0h, xp0, xnn0][t] /. soln,
    {t, 0, mt}, GridLines -> Automatic,
    Frame -> True, 
    FrameLabel -> {"t", "\!\(\*SubscriptBox[\(X\), \(c\)]\)", 
      "Initial condition is \!\(\*SubscriptBox[\(X\), \(c0h\)]\)"}],
   Plot[xp[xc0h, xp0, xnn0][t] /. soln,
    {t, 0, mt}, GridLines -> Automatic,
    Frame -> True, 
    FrameLabel -> {"t", "\!\(\*SubscriptBox[\(X\), \(p\)]\)", 
      "Initial condition is \!\(\*SubscriptBox[\(X\), \(c0h\)]\)"}]}

enter image description here

Note that the DE for $x_p$ is not coupled to $x_c$ or to $x_{nn}$, so its solution is independent of the initial conditions $x_c(0)$ and $x_{nn}(0)$. It may be a little misleading to label the third plot with the same title as others. It is also misleading to not put the other ICs in the plot titles as well.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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