5
$\begingroup$

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!

$\endgroup$
9
$\begingroup$

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]

enter image description here

$\endgroup$
  • $\begingroup$ interesting, your default plot colors are not the v9 defaults... :) $\endgroup$ – Sjoerd C. de Vries Apr 12 '14 at 18:05
  • $\begingroup$ Are these V 10 plots? They look different from V9 screen shot: !Mathematica graphics The lines are thicker also. $\endgroup$ – Nasser Apr 12 '14 at 19:35
  • $\begingroup$ Well "it is not V9" says it all ;-) You can actually get same thing now on Raspberry Pi. $\endgroup$ – Vitaliy Kaurov Apr 12 '14 at 20:15
2
$\begingroup$

Some modifications for the code in Vitaliy Kaurov's answer

Initial condition outside the limit cycle

\[Alpha] = 200/207; 
a = 4; 
k = 3; 
\[Beta] = 200/69 + 1/50; 
m = 1/2; 
c = 4761/20000; 
\[Gamma] = 1/10; 
E1 = 500/621; 
q1 = 2/10; 
E2 = 3/4; 
q2 = 2/100; 
Tf = 1000; 
eqs1 = {Derivative[1][x][t] == (-E1)*q1*x[t] + x[t]*(1 - x[t]/k)*\[Alpha] - 
      ((1 - m)*x[t]*y[t]*\[Beta])/(1 + a*(1 - m)*x[t]), 
    Derivative[1][y][t] == (-E2)*q2*y[t] + (c*(1 - m)*x[t]*y[t]*\[Beta])/
       (1 + a*(1 - m)*x[t]) - y[t]*\[Gamma], x[0] == 1.1, y[0] == 1.1}; 
s1 = NDSolve[eqs1, {x, y}, {t, Tf}]; 
Plot[Evaluate[{x[t], y[t]} /. First[s1]], {t, 0, Tf}, 
  PlotRange -> Full, PlotStyle -> {{Blue, Thickness[0.005]}, 
    {Red, Thickness[0.005]}}, Frame -> True, 
  FrameLabel -> {Style["Time", 18], Style["Population densities", 18]}, 
  LabelStyle -> Directive[Black, Bold, Medium], RotateLabel -> True, 
  PlotPoints -> 500, ImageSize -> Large, DisplayFunction -> Identity]
ParametricPlot[Evaluate[{x[t], y[t]} /. s1], {t, 0, Tf}, 
  PlotStyle -> {Black, Thickness[0.003]}, AspectRatio -> 1, 
  PlotRange -> All, Frame -> True, RotateLabel -> True, 
  FrameLabel -> {Style["Prey population density", 20], 
    Style["Predator population density", 20]}, RotateLabel -> False, 
  LabelStyle -> Directive[Black, Bold, Medium], ImageSize -> Large]

enter image description here

enter image description here

Initial condition inside the limit cycle

\[Alpha] = 200/207; 
a = 4; 
k = 3; 
\[Beta] = 200/69 + 1/50; 
m = 1/2; 
c = 4761/20000; 
\[Gamma] = 1/10; 
E1 = 500/621; 
q1 = 2/10; 
E2 = 3/4; 
q2 = 2/100; 
Tf = 1000; 
eqs2 = {Derivative[1][x][t] == (-E1)*q1*x[t] + x[t]*(1 - x[t]/k)*\[Alpha] - 
      ((1 - m)*x[t]*y[t]*\[Beta])/(1 + a*(1 - m)*x[t]), 
    Derivative[1][y][t] == (-E2)*q2*y[t] + (c*(1 - m)*x[t]*y[t]*\[Beta])/
       (1 + a*(1 - m)*x[t]) - y[t]*\[Gamma], x[0] == 1.02, y[0] == 1.02}; 
s2 = NDSolve[eqs2, {x, y}, {t, Tf}]; 
Plot[Evaluate[{x[t], y[t]} /. First[s2]], {t, 0, Tf}, 
  PlotRange -> Full, PlotStyle -> {{Blue, Thickness[0.005]}, 
    {Red, Thickness[0.005]}}, Frame -> True, 
  FrameLabel -> 
   {Style["\!\(\*\nStyleBox[\"Time\",\nFontColor->GrayLevel[0]]\)", 
     18], Style["\!\(\*\nStyleBox[\"Population\",\nFontColor->GrayLevel\
[0]]\)\!\(\*\nStyleBox[\" \
\",\nFontColor->GrayLevel[0]]\)\!\(\*\nStyleBox[\"densities\",\nFontCol\
or->GrayLevel[0]]\)", 18]}, LabelStyle -> Directive[Black, Bold, 
    Medium], RotateLabel -> True, PlotPoints -> 500, 
  ImageSize -> Large, DisplayFunction -> Identity]
ParametricPlot[Evaluate[{x[t], y[t]} /. s2], {t, 0, Tf}, 
  PlotStyle -> {Black, Thickness[0.003]}, AspectRatio -> 1, 
  PlotRange -> All, Frame -> True, RotateLabel -> True, 
  FrameLabel -> {Style["\!\(\*\nStyleBox[\"Prey\",\nFontColor->GrayLeve\
l[0]]\)\!\(\*\nStyleBox[\" \
\",\nFontColor->GrayLevel[0]]\)\!\(\*\nStyleBox[\"population\",\nFontCo\
lor->GrayLevel[0]]\)\!\(\*\nStyleBox[\" \
\",\nFontColor->GrayLevel[0]]\)\!\(\*\nStyleBox[\"density\",\nFontColor\
->GrayLevel[0]]\)", 20], Style["\!\(\*\nStyleBox[\"Predator\",\nFontCol\
or->GrayLevel[0]]\)\!\(\*\nStyleBox[\" \
\",\nFontColor->GrayLevel[0]]\)\!\(\*\nStyleBox[\"population\",\nFontCo\
lor->GrayLevel[0]]\)\!\(\*\nStyleBox[\" \
\",\nFontColor->GrayLevel[0]]\)\!\(\*\nStyleBox[\"density\",\nFontColor\
->GrayLevel[0]]\)", 20]}, RotateLabel -> False, 
  LabelStyle -> Directive[Black, Bold, Medium], ImageSize -> Large]

enter image description here enter image description here

First phase portrait

g1 = ParametricPlot[Evaluate[{x[t], y[t]} /. s2], {t, 0, Tf}, 
  PlotStyle -> {Black, Thickness[0.003]}, AspectRatio -> 1, 
  PlotRange -> All, Frame -> True, RotateLabel -> True, 
  FrameLabel -> { Style["\!\(\*
StyleBox[\"Prey\",\nFontColor->GrayLevel[0]]\)\!\(\*
StyleBox[\" \",\nFontColor->GrayLevel[0]]\)\!\(\*
StyleBox[\"population\",\nFontColor->GrayLevel[0]]\)\!\(\*
StyleBox[\" \",\nFontColor->GrayLevel[0]]\)\!\(\*
StyleBox[\"density\",\nFontColor->GrayLevel[0]]\)", 20], Style["\!\(\*
StyleBox[\"Predator\",\nFontColor->GrayLevel[0]]\)\!\(\*
StyleBox[\" \",\nFontColor->GrayLevel[0]]\)\!\(\*
StyleBox[\"population\",\nFontColor->GrayLevel[0]]\)\!\(\*
StyleBox[\" \",\nFontColor->GrayLevel[0]]\)\!\(\*
StyleBox[\"density\",\nFontColor->GrayLevel[0]]\)", 20]}, 
  RotateLabel -> False, LabelStyle -> Directive[Black, Bold, Medium], 
  ImageSize -> Large]

enter image description here

Second phase portrait

g2 = ParametricPlot[Evaluate[{x[t], y[t]} /. s1], {t, 0, Tf}, 
  PlotStyle -> {Black, Thickness[0.003]}, AspectRatio -> 1, 
  PlotRange -> All, Frame -> True, RotateLabel -> True, 
  FrameLabel -> { Style["\!\(\*
StyleBox[\"Prey\",\nFontColor->GrayLevel[0]]\)\!\(\*
StyleBox[\" \",\nFontColor->GrayLevel[0]]\)\!\(\*
StyleBox[\"population\",\nFontColor->GrayLevel[0]]\)\!\(\*
StyleBox[\" \",\nFontColor->GrayLevel[0]]\)\!\(\*
StyleBox[\"density\",\nFontColor->GrayLevel[0]]\)", 20], Style["\!\(\*
StyleBox[\"Predator\",\nFontColor->GrayLevel[0]]\)\!\(\*
StyleBox[\" \",\nFontColor->GrayLevel[0]]\)\!\(\*
StyleBox[\"population\",\nFontColor->GrayLevel[0]]\)\!\(\*
StyleBox[\" \",\nFontColor->GrayLevel[0]]\)\!\(\*
StyleBox[\"density\",\nFontColor->GrayLevel[0]]\)", 20]}, 
  RotateLabel -> False, LabelStyle -> Directive[Black, Bold, Medium], 
  ImageSize -> Large]

enter image description here

Finally, when applying the command Show[g2,g1] produces

Show[g2, g1, FrameLabel -> {Style["Prey population density", 20], 
    Style["Predator population density", 20]}]

enter image description here

$\endgroup$
  • 2
    $\begingroup$ Please post code in InputForm if you want to contribute anything with code. $\endgroup$ – J. M. is away Nov 10 '17 at 8:14
  • $\begingroup$ Ok, I appreciate the suggestion. $\endgroup$ – E. Chan-López Nov 10 '17 at 8:23
  • 2
    $\begingroup$ 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 $\endgroup$ – m_goldberg Nov 10 '17 at 11:32
  • $\begingroup$ Thanks for your support M. Goldberg. $\endgroup$ – E. Chan-López Nov 10 '17 at 20:16

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.