# Harvested prey–predator model incorporating a prey refuge

I want to plot the phase diagram of prey predator versus prey refuge to see how the prey refuge influences the population of prey and predator. And this is the system

$x'=\alpha x(1-x/k)-\beta\frac{(1-m)xy}{1+a(1-m)x}-q_1E_1x$

$y'=-\gamma y+c\beta\frac{(1-m)xy}{1+a(1-m)x}-q_2E_2y$

The prey predator with Holling type II model is incorporating a prey refuge, $mx$ and $k$, $\alpha$, $\gamma$, $c$ and $\beta/\alpha$ are the carrying capacity, growth rate of prey, death rate of predator, conversion factor denoting the number of newly born predators for each captured prey and maximum number of prey that can be eaten by each predator in unit time respectively.

I have the numerical value for $a=0.02$, $k=100$, $\alpha=10$, $\beta=0.6$, $\gamma=0.09$, $c=0.02$.

Thanks so much!

• I see few demos here on this subject demonstrations.wolfram.com/search.html?query=predator%20prey Commented Apr 12, 2014 at 6:26
• An article about this system can be found here. Commented Apr 12, 2014 at 6:53
• I think NDSolve is the function you are looking for. Commented Apr 12, 2014 at 9:44
• A posts with 1000s of constants must miss a few ;-) Your m?? Commented Apr 12, 2014 at 9:59

## 2 Answers

I solved more general system linked by @SjoerdC.deVries in the comments reproducing figure 3 and 4 - to prove it is correct. You can simplify this to version you need.

Clear["Global*"]

al = 2;
a = 2/1000;
k = 600;
b = 1/10;
g = 46/10^5;
c = 1/100;
m = 1/100;
E1 = 1;
q1 = 2/10;
E2 = 813/1000;
q2 = 2/100;
Tf = 300;

eqs = {
x'[t] ==
al x[t] (1 - x[t]/k) - b (1 - m) x[t] y[t]/(1 + a (1 - m) x[t]) -
q1 E1 x[t],
y'[t] == -g y[t] + c b (1 - m) x[t] y[t]/(1 + a (1 - m) x[t]) -
q2 E2 y[t],
x[0] == 2, y[0] == 8};

s = NDSolve[eqs, {x, y}, {t, Tf}];

Plot[Evaluate[{x[t], y[t]} /. s], {t, 0, Tf}, PlotStyle -> Automatic,
ImageSize -> 300, PlotRange -> All, Frame -> True]
ParametricPlot[Evaluate[{x[t], y[t]} /. s], {t, 0, Tf},
AspectRatio -> 1, PlotRange -> All, ImageSize -> 300, Frame -> True]


• interesting, your default plot colors are not the v9 defaults... :) Commented Apr 12, 2014 at 18:05
• Are these V 10 plots? They look different from V9 screen shot: !Mathematica graphics The lines are thicker also. Commented Apr 12, 2014 at 19:35
• Well "it is not V9" says it all ;-) You can actually get same thing now on Raspberry Pi. Commented Apr 12, 2014 at 20:15

Some modifications for the code in Vitaliy Kaurov's answer:

<<MaTeX
<< c:\CurvesGraphics6\CurvesGraphics6.m


Initial condition outside the limit cycle:

α = 200/207;
a = 4;
k = 3;
β = 200/69 + 1/50;
m = 1/2;
c = 4761/20000;
γ = 1/10;
E1 = 500/621;
q1 = 2/10;
E2 = 3/4;
q2 = 2/100;
Tf = 1000;
eqs1 = {x'[t] == -E1 q1 x[t] + x[t] (1 - x[t]/k) α - ((1 - m) x[t]*y[t] β)/(1 + a (1 - m) x[t]),
y'[t] == -E2 q2 y[t] + (c (1 - m) x[t]*y[t] β)/(1 + a (1 - m) x[t]) - y[t]γ, x[0] == 1.1, y[0] == 1.1};
s1 = NDSolve[eqs1, {x, y}, {t, Tf}];

ts1 = Plot[Evaluate[{x[t], y[t]} /. First[s1]], {t, 0, Tf},
PlotRange -> All,
PlotStyle -> {{Blue, Thickness[0.003]}, {Red, Thickness[0.003]}},
Frame -> True,
FrameLabel -> { MaTeX["\\text{t}", Magnification -> 1],
MaTeX["\\text{Population densities}", Magnification -> 1]},
LabelStyle -> Directive[Black, Tiny], RotateLabel -> True,
PlotPoints -> 500, ImageSize -> {300, 200},
DisplayFunction -> Identity]

pp1 = ParametricPlot[Evaluate[{x[t], y[t]} /. s1], {t, 0, 100},
Axes -> False, Oriented -> True, ArrowPositions -> {0.05, 0.985},
PlotStyle -> {Black, Thickness[0.002], Arrowheads[0.035]},
AspectRatio -> 1, PlotRange -> All, Frame -> True,
RotateLabel -> True,
FrameLabel -> {
MaTeX["\\text{Prey population density}", Magnification -> 1],
MaTeX["\\text{Predator population density}", Magnification -> 1]},
RotateLabel -> False, LabelStyle -> Directive[Black, Tiny],
ImageSize -> {250, 250}]


Initial condition inside the limit cycle:

α = 200/207;
a = 4;
k = 3;
β = 200/69 + 1/50;
m = 1/2;
c = 4761/20000;
γ = 1/10;
E1 = 500/621;
q1 = 2/10;
E2 = 3/4;
q2 = 2/100;
Tf = 1000;
eqs2 = {x'[t] == -E1 q1 x[t] + x[t] (1 - x[t]/k) α - ((1 - m) x[t]*y[t] β)/(1 + a (1 - m) x[t]),
y'[t] == -E2 q2 y[t] + (c (1 - m) x[t]*y[t] β)/(1 + a (1 - m) x[t]) - y[t]γ, x[0] == 1.02, y[0] == 1.02};
s2 = NDSolve[eqs2, {x, y}, {t, Tf}];

ts2 = Plot[Evaluate[{x[t], y[t]} /. First[s2]], {t, 0, Tf},
PlotRange -> All,
PlotStyle -> {{Blue, Thickness[0.003]}, {Red, Thickness[0.003]}},
Frame -> True,
FrameLabel -> { MaTeX["\\text{t}", Magnification -> 1],
MaTeX["\\text{Population densities}", Magnification -> 1]},
LabelStyle -> Directive[Black, Tiny], RotateLabel -> True,
PlotPoints -> 500, ImageSize -> {300, 300},
DisplayFunction -> Identity]

pp2 = ParametricPlot[Evaluate[{x[t], y[t]} /. s2], {t, 900, Tf},
PlotStyle -> {Black, Thickness[0.008]}, AspectRatio -> 1,
PlotRange -> All, Frame -> True, RotateLabel -> True,
FrameLabel -> {
MaTeX["\\text{Prey population density}", Magnification -> 1],
MaTeX["\\text{Predator population density}", Magnification -> 1]},
RotateLabel -> False, LabelStyle -> Directive[Black, Tiny],
ImageSize -> {250, 250}]


The initial conditions and equilibrium point:

p0 = Part[{x, y} /. NSolve[{-E1 q1 x + x (1 - x/k) α - ((1 - m) x*y β)/(1 + a (1 - m) x), -E2 q2 y + (c (1 - m) x*y β)/(1 + a (1 - m) x) - y γ} == {0, 0} && x > 0 && y > 0, {x, y}], 1];
(*{0.97972, 0.992966}*)
IC1 = {1.02, 1.02};
IC2 = {1.1, 1.1};

pointIC1 =
Graphics[{PointSize[0.015], Darker[Green], Point[IC1]},
IC1 - {-0.08, -0.011}, Axes -> True, PlotRange -> All,
ImageSize -> {250, 250}];

pointIC2 =
Graphics[{PointSize[0.015], Darker[Green], Point[IC2]},
IC2 - {-0.08, -0.011}, Axes -> True, PlotRange -> All,
ImageSize -> {250, 250}];


Finally, when applying the command Show produces

Show[pp2, pp21, pp1, pointIC1, pointIC2,
Epilog -> {{PointSize[0.018], Black, Point[p0]}, {PointSize[0.013],
White, Point[p0]}}, Axes -> False, ImageSize -> {350, 350}]


• Please post code in InputForm if you want to contribute anything with code. Commented Nov 10, 2017 at 8:14
• Ok, I appreciate the suggestion. Commented Nov 10, 2017 at 8:23
• This looks like it would be a valuable contribution if it were posted properly. Please read the instructions on formatting code given here and here Commented Nov 10, 2017 at 11:32
• Thanks for your support M. Goldberg. Commented Nov 10, 2017 at 20:16