Rosenzweig-MacArthur predator-prey model

Question

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.

Answer 1

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"}]]

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]

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}}]

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]}}]