3
$\begingroup$

I have the following ODE system

dx1 = -0.7 x1 y2 
dy1 = 0.7 x1 y2 + 1.7 (1-x1-y1) y2 - y1 
dx2 = -0.7 x2 y1 
dy2 = 0.7 x2 y1 + 1.7 (1-x2-y2) y1 - y2

I have found that, for initial conditions of the type $x_1=1-y_1, x_2=1-y_2$, solutions either:

  • reach a steady-state where $y_1=y_2=0$
  • or reach a steady-state where $y_1,y_2>0$.

Now, I want to plot the plane of initial conditions $(y_1(0),y_2(0))$ and color each point in the plane blue or red depending on whether the steady-state falls in the first or second bucket.

I know how to use streamplot to explore the phase-space of solutions but I am not sure how to go about plotting the regions depending on the steady-state value. Any ideas would be much appreciated.

$\endgroup$
  • $\begingroup$ By the way, what's the system a model of? $\endgroup$ – Chris K Feb 18 at 20:18
  • 1
    $\begingroup$ It's an extension of an epidemic model that exhibits bistability in the initial conditions. See this paper $\endgroup$ – rpa Feb 18 at 20:21
3
$\begingroup$

Here's an approach based on this answer by @ssch.

First, a function to solve the system and report the final values of y1[tmax] and y2[tmax]:

res[y10_?NumericQ, y20_?NumericQ] := Module[{sol},
  sol = NDSolve[{
    x1'[t] == -0.7 x1[t] y2[t],
    y1'[t] == 0.7 x1[t] y2[t] + 1.7 (1 - x1[t] - y1[t]) y2[t] - y1[t],
    x2'[t] == -0.7 x2[t] y1[t],
    y2'[t] == 0.7 x2[t] y1[t] + 1.7 (1 - x2[t] - y2[t]) y1[t] - y2[t],
    x1[0] == 1 - y10, y1[0] == y10, x2[0] == 1 - y20, y2[0] == y20},
    {x1, y1, x2, y2}, {t, 0, tmax}][[1]];
  Return[{y1[tmax], y2[tmax]} /. sol]
];

Then RegionPlot where the final values are close enough to {0, 0}:

tmax = 100;
tol = 10^-3;

RegionPlot[{Norm[res[y10, y20]] < tol, Norm[res[y10, y20]] > tol},
  {y10, 0, 0.2}, {y20, 0, 0.2}, 
  PlotPoints -> 20, MaxRecursion -> 2, PlotStyle -> {Blue, Red}, FrameLabel -> {y10, y20}]

Mathematica graphics

See also this question.

$\endgroup$
  • $\begingroup$ Thanks!! I ended up figuring out another way and I posted my solution as the same time as yours. However, I'm accepting your answer as I think it's better! $\endgroup$ – rpa Feb 18 at 20:18
2
$\begingroup$

Here's what I ended up doing in case anyone else has a similar question:

I used NDSolve to numerically solve the system. I then extracted the values of $y_1$ and $y_2$ at time $t=100$, which I know is sufficiently long for the system to reach steady-state. I then threshold the average of $y_1$ and $y_2$ and set the threshold close to zero. Then, I did this for multiple values of the initial conditions using Table and plotted the resulting array using arrayplot.

ArrayPlot[
 Table[(i1[100] + i2[100])/2 /. 
   NDSolve[{s1'[t] == -0.7 s1[t] i2[t], 
      i1'[t] ==  0.7 s1[t] i2[t] + 1.7 (1 - s1[t] - i1[t]) i2[t] - i1[t], 
      s2'[t] == -0.7 s2[t] i1[t], 
      i2'[t] == 0.7 s2[t] i1[t] + 1.7 (1 - s2[t] - i2[t]) i1[t] - i2[t], 
      i1[0] == p1, i2[0] == p2, s1[0] == 1 - p1, s2[0] == 1 - p2},
      {i1, i2}, {t, 0, 100}][[1]], {p1, 0, 1, 0.005}, {p2, 0, 1, 0.005}], 
   ColorFunction -> (If[# > 0.02, Red, Blue] &), DataReversed -> True]

enter image description here

$\endgroup$

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.