8
$\begingroup$

The predator-prey model is governed by the following system of ode's. \begin{eqnarray} &&\displaystyle{\frac{dx}{dt}=r x\left(1 - \frac{x}{K}\right) - \frac{s y x}{1 + s \tau x}},\\[0.1cm] &&\displaystyle{\frac{dy}{dt}=-c y +d \frac{s x y}{1 + s \tau x}}, \end{eqnarray} where $r,K,c,d>0$, $s$ is the predator search rate and $\tau$ is time. The background of this model can be found here.

Using appropriate scaling as discussed here, the above the system can be written in the dimensionless form,

\begin{eqnarray} &&\displaystyle{\frac{dx}{dt}=x\left(1 - \frac{x}{k}\right) - \frac{m y x}{1 + x}},\\[0.1cm] &&\displaystyle{\frac{dy}{dt}=-c y + \frac{m x y}{1 + x}}. \end{eqnarray}

How to get solution of the above system?

How one can get phase portrait of the system?

Any help is highly appreciated.

$\endgroup$
15
  • $\begingroup$ the second equation should begin with $dy(x)/dt$ ? $\endgroup$ Commented Feb 18, 2017 at 14:25
  • $\begingroup$ Please write the equations in Mathematica syntax as well, and show what you tried. I assume you tried at least DSolve, and that if it couldn't solve it, you tried a numerical solution using NDSolve. What went wrong when you tried these? Explain where the difficulty is. $\endgroup$
    – Szabolcs
    Commented Feb 18, 2017 at 14:36
  • $\begingroup$ yes, I am sorry, edited. $\endgroup$ Commented Feb 18, 2017 at 14:37
  • $\begingroup$ What is the origin of this model? $\endgroup$
    – zhk
    Commented Feb 18, 2017 at 15:04
  • 4
    $\begingroup$ @SkSarifHassan Here it is mathematica.stackexchange.com/questions/45922/… $\endgroup$
    – zhk
    Commented Feb 18, 2017 at 15:12

1 Answer 1

15
$\begingroup$

Note

In my attempt to answer the OP's question, I have presented all most all the visuals/graphs which are important for the analysis of such models. If there is something missing or physically incorrect then please feel free to edit and correct?

Credit goes to @ChrisK for pointing me in the right direction, which made me able to carry out correct graphical analysis of the model.

The equations and values for the different parameters are take from the document cited in the question.

Solution & Plot

To solve such nonlinear systems, the best choice, in almost all cases is to use NDSolve to get numerical solutions.

c = 1; m = 3;
sol = With[{k = 3}, First@NDSolve[{x'[t] == x[t] (1 - x[t]/k) - m*x[t]*y[t]/(1 + x[t]), 
     y'[t] == m*x[t]*y[t]/(1 + x[t]) - c*y[t], x[0] == 1, 
     y[0] == 1}, {x, y}, {t, 0, 50}]]

Plot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 50}, 
PlotStyle -> {Red, Directive[Green, Dashed]}, Frame -> True,
PlotLegends -> LineLegend[{Red, Directive[Green, Dashed]}, {"prey", "predator"}]]

enter image description here

Nullclines

nk1 = With[{k = 3}, ContourPlot[{((1 - x/k) - m*y/(1 + x)), (m*x/(1 + x) - c), 
     y}, {x, -1.5, 3.5}, {y, -0.5, 1}, 
    Epilog -> {Text[Style["k>(m+c)/(m-c)", 12], Scaled[{0.7, 0.9}]]}]];

nk2 = With[{k = 1}, ContourPlot[{((1 - x/k) - m*y/(1 + x)), (m*x/(1 + x) - c), 
     y}, {x, -1.5, 3.5}, {y, -0.5, 1}, 
    Epilog -> {Text[Style["k<(m+c)/(m-c)", 12], Scaled[{0.7, 0.9}]]}]];

GraphicsGrid[{{nk1, nk2}}, ImageSize -> Large]

enter image description here

Phase Portrait

sol1[k_?NumericQ, x0_?NumericQ] := 
  First@NDSolve[{x'[t] == x[t] (1 - x[t]/k) - m*x[t]*y[t]/(1 + x[t]), 
     y'[t] == m*x[t]*y[t]/(1 + x[t]) - c*y[t], x[0] == x0, 
     y[0] == x0}, {x, y}, {t, 0, 200}];

ppk1 = ParametricPlot[
   Evaluate[{x[t], y[t]} /. sol1[1, #] & /@ Range[0, 1, 0.2]], {t, 0, 
    200}, Frame -> True, 
   Epilog -> {Text[Style["k=1", 20], Scaled[{0.8, 0.8}]]}, 
   ImageSize -> 200, PlotRange -> {{0, 2.5}, {0, 1.6}}];

ppk3 = ParametricPlot[
   Evaluate[{x[t], y[t]} /. sol1[3, #] & /@ Range[0, 1, 0.2]], {t, 0, 
    200}, Frame -> True, 
   Epilog -> {Text[Style["k=3", 20], Scaled[{0.8, 0.8}]]}, 
   ImageSize -> 200, PlotRange -> {{0, 2.5}, {0, 1.6}}];

GraphicsGrid[{{ppk1, ppk3}}]

enter image description here

sp1 = With[{k = 1}, 
   StreamPlot[{x*(1 - x/k) - m*x*y/(1 + x), m*x*y/(1 + x) - c*y}, {x, 
     0, 2.5}, {y, 0, 1.5}]];

sp2 = With[{k = 3}, 
   StreamPlot[{x*(1 - x/k) - m*x*y/(1 + x), m*x*y/(1 + x) - c*y}, {x, 
     0, 2.5}, {y, 0, 1.5}]];

GraphicsGrid[{{Show[ppk1, sp1, nk2], Show[ppk3, sp2, nk1]}}]

enter image description here

$\endgroup$
7
  • 2
    $\begingroup$ This is a nice start. For more interesting results, you should increase c relative to d (otherwise the predators can't persist no matter how much prey they have) and then use prey carrying capacity k as a bifurcation parameter. This is the well-known Rosenzweig-MacArthur predator-prey model. As you increase k there is a Hopf bifurcation where a predatory-prey limit cycle is born. $\endgroup$
    – Chris K
    Commented Feb 18, 2017 at 16:36
  • $\begingroup$ @MMM, how to achieve the bifurcation diagram of the model with respect to the parameter k? $\endgroup$ Commented Feb 19, 2017 at 4:08
  • 1
    $\begingroup$ @SkSarifHassan Check my edited response. $\endgroup$
    – zhk
    Commented Feb 19, 2017 at 6:39
  • $\begingroup$ @SkSarifHassan Is this what you wanted? or I am missing something? $\endgroup$
    – zhk
    Commented Feb 19, 2017 at 13:23
  • $\begingroup$ @MMM, yes it is. I was thinking from complex variable point of view. What about bifurcation diagram for k? $\endgroup$ Commented Feb 19, 2017 at 14:35

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