0
$\begingroup$

Suppose we have the model:

tmax = 2000;
\[Beta]1 = 0.0001;
\[Beta]2 = 0.006;
p = 0.9;
q = 0.8;
\[Xi]1 = 0.8;
\[Xi]2 = 0.9;
\[Epsilon] = 0.002;
p1 = 0.01;
p2 = 0.03;
\[Alpha] = 0.01;
\[Mu] = 0.01;
\[Nu] = 0.55;
SIIJA = NDSolveValue[{
    S'[t] == \[Nu] - \[Beta]1*S[t]*I2[t] - \[Beta]2*S[t]*J[t] - \[Mu]*
       S[t],
    I1'[t] == 
     p*\[Beta]1*S[t]*I2[t] + 
      q*\[Beta]2*S[t]*J[t] + \[Xi]1 *J[t] - (\[Epsilon] + \[Mu])*I1[t],
    I2'[t] == (1 - p) \[Beta]1*S[t]*I2[t] + (1 - q)*\[Beta]2*S[t]*
       J[t] + \[Epsilon]*I1[t] + \[Xi]2*J[t] - (p1 + \[Mu])*I2[t],
    J'[t] == p1*I2[t] - (\[Xi]1 + \[Xi]2 + p2 + \[Mu])*J[t],
    A'[t] == p2*J[t] - (\[Alpha] + \[Mu])*A[t] ,
    S[0] == 50,
    I1[0] == 50,
    I2[0] == 50,
    J[0] == 20,
    A[0] == 20},
   {S, I1, I2, J, A},
   {t, 0, tmax}];
{f1, f2, f3, f4, f5} = SIIJA;

st = Style[#, 15, Black] &;

Plot[{f1[t], f2[t], f3[t], f4[t], f5[t]}, {t, 0, tmax}, 
 PlotStyle -> {Blue, Purple, Yellow, Orange, Red}, Frame -> True, 
 FrameLabel -> st /@ {"Time", "Density"}, 
 PlotLegends -> 
  Placed[LineLegend[{Blue, Purple, Yellow, Orange, Red}, {"S(t)", 
     "I1(t)", "I2(t)", "J(t)", "A(t)" }, 
    LegendFunction -> Framed], {0.85, 0.45}], ImageSize -> 500]

Giving:

enter image description here

How would I find the value of say I1, I2 etc at time 1000 for example from the graph?

Edit: From the original model in the original post, we had total population converging to $55$, this was just as expected as $N$ will converge to $\frac{\nu}{\mu}$ when $t \rightarrow \infty$. However, consider the same model with different parameter values:

tmax = 20000;
\[Beta]1 = 0.0001;
\[Beta]2 = 0.006;
p = 0.3;
q = 0.4;
\[Xi]1 = 0.001;
\[Xi]2 = 0.003;
\[Epsilon] = 0.0002;
p1 = 0.01;
p2 = 0.03;
\[Alpha] = 0.01;
\[Mu] = 0.01;
\[Nu] = 0.55;
SIIJA = NDSolveValue[{
    S'[t] == \[Nu] - \[Beta]1*S[t]*I2[t] - \[Beta]2*S[t]*J[t] - \[Mu]*
       S[t],
    I1'[t] == 
     p*\[Beta]1*S[t]*I2[t] + 
      q*\[Beta]2*S[t]*J[t] + \[Xi]1 *J[t] - (\[Epsilon] + \[Mu])*I1[t],
    I2'[t] == (1 - p) \[Beta]1*S[t]*I2[t] + (1 - q)*\[Beta]2*S[t]*
       J[t] + \[Epsilon]*I1[t] + \[Xi]2*J[t] - (p1 + \[Mu])*I2[t],
    J'[t] == p1*I2[t] - (\[Xi]1 + \[Xi]2 + p2 + \[Mu])*J[t],
    A'[t] == p2*J[t] - (\[Alpha] + \[Mu])*A[t] ,
    S[0] == 50,
    I1[0] == 50,
    I2[0] == 50,
    J[0] == 20,
    A[0] == 20},
   {S, I1, I2, J, A},
   {t, 0, tmax}];
{f1, f2, f3, f4, f5} = SIIJA;

st = Style[#, 15, Black] &;

Plot[{f1[t], f2[t], f3[t], f4[t], f5[t]}, {t, 0, tmax}, 
 PlotStyle -> {Blue, Purple, Yellow, Orange, Red}, Frame -> True, 
 FrameLabel -> st /@ {"Time", "Density"}, 
 PlotLegends -> 
  Placed[LineLegend[{Blue, Purple, Yellow, Orange, Red}, {"S(t)", 
     "I1(t)", "I2(t)", "J(t)", "A(t)" }, 
    LegendFunction -> Framed], {0.85, 0.75}], ImageSize -> 500]

When we do f1[20000] + f2[20000] + f3[20000] + f4[20000] + f5[20000] it gives us 51.3631 which is not what we want.. why is it doing this?

$\endgroup$
11
  • 1
    $\begingroup$ Right klick with the mouse in the graphics and choose: "GetCoordinates" $\endgroup$ Nov 30 '21 at 16:56
  • $\begingroup$ @DanielHuber These are not accurate.. And since there's a overlap, is there a way to type something in Mathematica like N[S(1000)]? $\endgroup$
    – Math
    Nov 30 '21 at 16:57
  • 2
    $\begingroup$ Why you do not simply write e.g.: f1[1000]? $\endgroup$ Nov 30 '21 at 17:01
  • $\begingroup$ @DanielHuber Thank you, I completely forgot about this $\endgroup$
    – Math
    Nov 30 '21 at 17:08
  • $\begingroup$ If you add up the right hand sides of your system to get N'[t], you get \[Nu] - (\[Alpha] + \[Mu]) A[t] - \[Mu] I1[t] - \[Mu] I2[t] - \[Mu] J[t] - \[Mu] S[t]. That \[Alpha] loss term is the leak. $\endgroup$
    – Chris K
    Dec 1 '21 at 13:54
2
$\begingroup$

Use Manipulate to specify time

Clear["Global`*"]

$Version

(* "12.3.1 for Mac OS X x86 (64-bit) (June 19, 2021)" *)

Constants

tmax = 2000;
β1 = 10^-4;
β2 = 6*^-3;
p = 9/10;
q = 4/5;
ξ1 = 4/5;
ξ2 = 9/10;
ϵ = 2*^-3;
p1 = 10^-2;
p2 = 3*^-2;
α = 10^-2;
μ = 10^-2;
ν = 11/20;

Solution

SIIJA = NDSolveValue[{
    S'[t] == ν - β1*S[t]*I2[t] - β2*S[t]*J[t] - μ*S[t],
    I1'[t] == p*β1*S[t]*I2[t] + 
      q*β2*S[t]*J[t] + ξ1*J[t] - (ϵ + μ)*I1[t],
    I2'[t] == (1 - p) β1*S[t]*I2[t] + (1 - q)*β2*S[t]*
       J[t] + ϵ*I1[t] + ξ2*J[t] - (p1 + μ)*I2[t],
    J'[t] == p1*I2[t] - (ξ1 + ξ2 + p2 + μ)*J[t],
    A'[t] == p2*J[t] - (α + μ)*A[t],
    S[0] == 50, I1[0] == 50, I2[0] == 50, J[0] == 20, A[0] == 20},
   {S, I1, I2, J, A}, {t, 0, tmax}];

{f1, f2, f3, f4, f5} = SIIJA;

st = Style[#, 15, Black] &;

Display

Manipulate[
 Column[{
   Plot[{f1[t], f2[t], f3[t], f4[t], f5[t]}, {t, 0, tmax},
    PlotStyle -> {Blue, Purple, Yellow, Orange, Red},
    Frame -> True,
    FrameLabel -> st /@ {"Time", "Density"},
    PlotLegends ->
     Placed[
      LineLegend[
       {"S", "I1", "I2", "J", "A"},
       LegendFunction -> Framed,
       Background -> White],
      {0.85, 0.45}],
    ImageSize -> 500,
    MaxRecursion -> 5,
    Prolog -> {AbsoluteThickness[1], Gray, Dashed,
      InfiniteLine[{{tv, 0}, {tv, 60}}]}],
   Row[
    Riffle[
     StringForm["`` = ``", #[[1]],
        ScientificForm[#[[2]], 3]] & /@
      Transpose[{{S, I1, I2, J, A},
        {f1[tv], f2[tv], f3[tv], f4[tv], f5[tv]}}],
     "; "]]},
  ItemSize -> {35, Automatic}],
 {{tv, 1000, t}, 0, tmax, 10, Appearance -> "Labeled"}]

enter image description here

$\endgroup$
1
  • $\begingroup$ Very nice solution, can you please check my edit in the original post? $\endgroup$
    – Math
    Dec 1 '21 at 12:49

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