1
$\begingroup$

I have the following system:

enter image description here

For this system I have to calculate the three equilibria (critical points). Here are the equations in Mathematica:

eqH = HH'[t] == (bH/NN)*(NN - CC[t] - HH[t])*HH[t] - (sH*aH)*HH[t]
eqC = CC'[t] == (bC/NN)*(NN - CC[t] - HH[t])*CC[t] - (sC + aC)*CC[t]
eqS = SS[t] == NN - HH[t] - CC[t]

L have the following values for the constants:

NN := 400
sS := 1/5
sC := 1/7
sH := 1/5
bC := 0.45
bH := 0.4
aC := 0.1
aH := 0.1

But now I don't know what to do next. For critical points I don't know how to solve the equations H'[t] == 0, C'[t] == 0, S[t] == 0

$\endgroup$
  • 1
    $\begingroup$ Have you seen Solve? $\endgroup$ – Michael E2 Feb 4 '18 at 23:18
  • $\begingroup$ @MichaelE2 I tried to solve the system like this: system = {eqH, eqC} , solSystem = DSolve[system, {hh, cc}, t] but it's not ok $\endgroup$ – Darius Ionut Feb 5 '18 at 0:05
  • $\begingroup$ You don't want to set S[t]==0 but you do need to set eqH and eqC equal to zero in Solve. Don't use DSolve since you're solving algebraic equations. $\endgroup$ – Chris K Feb 5 '18 at 2:40
3
$\begingroup$
Block[{Nc, Sc, dC, dH, bC, bH, aC, aH, St, Ht, Ct},

 With[{$Nc = 400, $Sc = 0.2, $dC = 1./7, $dH = 0.2, $bC = 0.45, $bH = 0.4, $aC = 0.1, $aH = 0.1},
  With[{paramsRl = Thread[{Nc, Sc, dC, dH, bC, bH, aC, aH} -> {$Nc, $Sc, $dC, $dH, $bC, $bH, $aC, $aH}]},

   With[{xRHS = (bH/Nc) (Nc - Ct - Ht) Ht - (dH + aH) Ht},
    With[{yRHS = (bC/Nc) (Nc - Ct - Ht) Ct - (dC + aC) Ct},
     With[{CST = Sc - Nc + Ht + Ct},

      Block[{horIsoclRl, vertIsoclRl, critPtsRl, expr},
       (* horizontal isoclines ie dC/dH \[Equal] 0 *)
       horIsoclRl = Solve[yRHS == 0, Ct];

       (* vertical isoclines ie dC/dH \[Equal] Infintiy *)
       vertIsoclRl = Solve[xRHS == 0, Ct];

       (* critical points ie Ht' \[Equal] 0 && Ct' \[Equal] 0 *)
       critPtsRl = Solve[{xRHS == 0, yRHS == 0}, {Ht, Ct}];

       Block[{isoclGraph, critPtsGraph, cstGraph, critPtsTab, paramsTab, C, H},

        (* plot isoclines - horizontal: Gray, vertical: Blue *)
        With[{xMn = -150, xMx = 350},
         isoclGraph = MapIndexed[
           Graphics[
             If[
               #2[[1]] == 1,
               {Gray, Thin, Opacity[0.5], 
                Line[{{xMn, #1[[-1]] /. Ht -> xMn}, {xMx, #1[[-1]] /. Ht -> xMx}}]},
               {Lighter[Blue], Thin, Opacity[0.5], Line[{{xMn, #1[[-1]] /. Ht -> xMn}, {xMx, #1[[-1]] /. Ht -> xMx}}]}
               ] /. paramsRl] &,
           {{Ht, Ct} /. horIsoclRl, {Ht, Ct} /. vertIsoclRl},
           {2}
           ];

         (* plot critical points - Red *)
         critPtsGraph = Graphics[
             {PointSize[Large], Darker[Red], Point[#]}
             ] & /@ (({Ht, Ct} /. critPtsRl) /. paramsRl);

         (* plot constraint - Brown *)
         cstGraph = MapIndexed[
           Graphics[
             {Brown, Thin, Dashed, Opacity[0.7], 
               Line[{{xMn, #1[[-1]] /. Ht -> xMn}, {xMx, #1[[-1]] /. Ht -> xMx}}]} /. paramsRl] &, (({Ht, Ct} /. Solve[CST == 0, Ct]) /. paramsRl)];

         (* create critical points' table *)
         critPtsTab = Grid[

           Partition[
            Style[{Ht, Ct}, Darker[Blue]] /. critPtsRl /. paramsRl, 2,
             2, {1, 1}, {}],
           Alignment -> Left
           ];

         (* create parameters' table *)
         paramsTab = Grid[
           Prepend[
            Partition[
             Thread[

              Style[#, Gray, 12] & /@ {N, S, Subscript[δ, C], 
                 Subscript[δ, H], Subscript[β, C], Subscript[β, H], 
                 Subscript[α, C], Subscript[α, H]} -> {$Nc, $Sc, $dC, $dH, $bC, $bH, $aC, $aH}
          ], 3, 3, {1, 1}, {}], {Style["Parameters", 22], 
         SpanFromLeft}], Alignment -> Center];

         (* assemble output *)
         Labeled[
          Legended[
           Show[isoclGraph, critPtsGraph, cstGraph,
            Frame -> True,
            ImageSize -> Medium,
            FrameLabel -> {
              {Style[C, Bold, 16], None},
              {Style[H, Bold, 16], None}
              },
            RotateLabel -> False
            ],
           Column[{
             LineLegend[
              {Directive[{Gray, Thick, Opacity[0.9]}],
               Directive[{Lighter[Blue], Thick, Opacity[0.9]}],
               Directive[{Brown, Thick, Dashed, Opacity[0.9]}]
               },
              {Style["Horizontal isoclines", 18],
               Style["Vertical isoclines", 18],
               Style["Constraint (?)", 18]
               }],
             Labeled[
              PointLegend[
               {Directive[{PointSize[Large], Darker[Red]}]},
               {Style["Critical points", 18]}
               ],
              critPtsTab
              ]
            }]
           ],
          paramsTab,
          Left
          ]

         ]

        ]


       ]

      ]
     ]
    ]
   ]
  ]
 ]

enter image description here

$\endgroup$
4
$\begingroup$

Are you looking for something like this?

NN = 400; sS = 1/5; sC = 1/7; sH = 1/5; bC = 0.45; bH = 0.4; aC = 0.1; aH = 0.1;

eqH = (bH/NN)*(NN - CC[t] - HH[t])*HH[t] - (sH*aH)*HH[t] == 0;

eqC = (bC/NN)*(NN - CC[t] - HH[t])*CC[t] - (sC + aC)*CC[t] == 0;

eqS = NN - HH[t] - CC[t] - SS[t] == 0

NSolve[{eqH, eqC, eqS}, {CC[t], HH[t], SS[t]}]
$\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.