6
$\begingroup$

I have a Mathematica assignment where we need to numerically solve a Lotka-Volterra predator-prey population equations. $$\frac{dx}{dt}=ax-bxy,\quad\frac{dy}{dt}=-cx+dxy,$$ where $a,b,c,d\in\mathbb{R}_+$.

So far, I have managed to come up with this for the representation of the equation in Mathematica:

s = DSolve[{x'[t] == a x[t] - b x[t] y[t], y'[t] == -c y[t] + d x[t] y[t]}, {x[t], y[t]}, t]

However, I am having trouble getting the equation to plot. I have been trying to use:

Plot[Evaluate[{x[t], y[t]} /. s], {t, 0, 20}, WorkingPrecision -> 20]

as I have found online, but nothing shows on the plot.

Can anyone tell me what I am doing wrong? I am new to Mathematica.

Improve this question
$\endgroup$
4
  • $\begingroup$ Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$
    – bbgodfrey
    Commented Feb 17, 2015 at 0:18
  • $\begingroup$ Be sure that all your constants have been assigned numerical values before plotting the solution. $\endgroup$
    – bbgodfrey
    Commented Feb 17, 2015 at 0:19
  • $\begingroup$ Thanks for the input! I tried setting all my variables, just by saying a=1, b=1, c=1, d=1 (or 1, 2, 3, 4) to no avail. $\endgroup$
    – elykl33t
    Commented Feb 17, 2015 at 0:55
  • $\begingroup$ It strongly looks as a homework. $\endgroup$
    – Alexei Boulbitch
    Commented Feb 16, 2017 at 19:11

4 Answers 4

Reset to default
18
$\begingroup$

The dimensional phase-space plot

The analysis of the predator-prey model is incomplete without the phase-space plot. So, here it is,

a = 2/3; b = 4/3; c = 1; d = 1;

sol[x0_?NumericQ] := First@NDSolve[{x'[t] == a*x[t] - b*x[t]*y[t], y'[t] == c*x[t]*y[t] -
d*y[t], x[0] == x0, y[0] == x0}, {x, y}, {t, 0, 20}];

ParametricPlot[Evaluate[{x[t], y[t]} /. sol[#] & /@ Range[0.9, 1.8, 0.1]], {t, 0, 20}, 
PlotRange -> All, Frame -> True]

enter image description here

The dimensionless Lotka–Volterra equations

Following the lecture notes for Mathematical Biology by Jeffrey R. Chasnov ( section 1.4 page no. 11), we can have the dimensionless equations in the following form,

$$\frac{d U}{dt}=r(U-UV),\quad \frac{dV}{dt}=\frac{1}{r}(UV-V),$$

where $r=\sqrt{\alpha/\beta}$, a nondimensional parameter. In the notes, the author has solved the above system using Matlab numerical solver ode45. Here, I will reproduce his results using Mathematica.

The solution to the above system is found using ParametricNDSolveValue for different values of $r$.

sol = ParametricNDSolveValue[{U'[t] == r*(U[t] - U[t]*V[t]), 
   V'[t] == 1/r*(U[t]*V[t] - V[t]), U[0] == 1.1, V[0] == 1}, {U,V}, {t, 0, 20}, {r}];

GraphicsColumn[{Plot[{Evaluate[sol[0.5][[1]]][t], 
    Evaluate[sol[0.5][[2]]][t]}, {t, 0, 19}, AspectRatio -> .4, 
   ImageSize -> 500, 
   PlotLegends -> 
    Placed[LineLegend[{Red, Directive[Green, Dashed]}, {"prey", 
       "predator"}], {0.89, 0.89}], PlotRange -> {{0, 19}, {0.8, 1.25}}, 
   PlotStyle -> {Red, Directive[Green, Dashed]}, Frame -> True, 
   FrameStyle -> AbsoluteThickness[.003], 
   FrameTicks -> {{Automatic, Automatic}, {Automatic, Automatic}}, 
   PlotRangePadding -> None, ImagePadding -> {{20, 1}, {20, 30}}, 
   Epilog -> {Text[Style["r=0.5", 20], Scaled[{0.6, 0.9}]]}], 
  Plot[{Evaluate[sol[1][[1]]][t], Evaluate[sol[1][[2]]][t]}, {t, 0, 
    19}, AspectRatio -> .4, PlotRange -> {{0, 19}, {0.8, 1.25}}, 
   PlotStyle -> {Red, Directive[Green, Dashed]}, Frame -> True, 
   FrameStyle -> AbsoluteThickness[.003], 
   FrameTicks -> {{Automatic, Automatic}, {Automatic, Automatic}}, 
   PlotRangePadding -> None, ImagePadding -> {{20, 1}, {20, 30}}, 
   Epilog -> {Text[Style["r=1", 20], Scaled[{0.6, 0.9}]]}], 
  Plot[{Evaluate[sol[2][[1]]][t], Evaluate[sol[2][[2]]][t]}, {t, 0, 
    19}, AspectRatio -> .4, PlotRange -> {{0, 19}, {0.8, 1.25}}, 
   PlotStyle -> {Red, Directive[Green, Dashed]}, Frame -> True, 
   FrameStyle -> AbsoluteThickness[.003], 
   FrameTicks -> {{Automatic, Automatic}, {Automatic, Automatic}}, 
   PlotRangePadding -> None, ImagePadding -> {{20, 1}, {20, 30}}, 
   Epilog -> {Text[Style["r=2", 20], Scaled[{0.6, 0.9}]]}]}, 
 ImageSize -> 500, Spacings -> -50]

enter image description here

Finally, plotting the phase-space plot for different values of $r$. For this I have utilized once again ?NumericQ.

sol[r_?NumericQ, U0_?NumericQ] :=First@NDSolve[{U'[t] == r*(U[t] - U[t]*V[t]),
V'[t] == 1/r*(U[t]*V[t] - V[t]), U[0] == U0, V[0] == U0}, {U,V}, {t, 0, 20}];

Grid[{{ParametricPlot[
    Evaluate[{U[t], V[t]} /. sol[0.5, #] & /@ 
      Range[0.9, 1.8, 0.1]], {t, 0, 20}, Frame -> True, 
    Epilog -> {Text[Style["r=0.5", 20], Scaled[{0.8, 0.8}]]}, 
    ImageSize -> 200, PlotRange -> {{0, 4}, {0, 4}}], 
   ParametricPlot[
    Evaluate[{U[t], V[t]} /. sol[1, #] & /@ Range[0.9, 1.8, 0.1]], {t,
      0, 20}, Frame -> True, 
    Epilog -> {Text[Style["r=1", 20], Scaled[{0.8, 0.8}]]}, 
    ImageSize -> 200, PlotRange -> {{0, 4}, {0, 4}}], 
   ParametricPlot[
    Evaluate[{U[t], V[t]} /. sol[2, #] & /@ Range[0.9, 1.8, 0.1]], {t,
      0, 20}, Frame -> True, 
    Epilog -> {Text[Style["r=2", 20], Scaled[{0.8, 0.8}]]}, 
    ImageSize -> 200, PlotRange -> {{0, 4}, {0, 4}}]}}]

enter image description here

The simply way to plot the phase planes of system is to use StreamPlot

Grid[{{With[{r = 0.5}, 
    StreamPlot[{r*(U - U*V), 1/r*(U*V - V)}, {U, 0, 5}, {V, 0, 5}, 
     Epilog -> {Text[Style["r=0.5", 20], Scaled[{0.8, 0.8}]]}]], 
   With[{r = 1}, 
    StreamPlot[{r*(U - U*V), 1/r*(U*V - V)}, {U, 0, 5}, {V, 0, 5}, 
     Epilog -> {Text[Style["r=1", 20], Scaled[{0.8, 0.8}]]}]], 
   With[{r = 2}, 
    StreamPlot[{r*(U - U*V), 1/r*(U*V - V)}, {U, 0, 5}, {V, 0, 5}, 
     Epilog -> {Text[Style["r=2", 20], Scaled[{0.8, 0.8}]]}]]}}]

enter image description here

Improve this answer
$\endgroup$
14
$\begingroup$

This code presented in this answer is essentially the same as the code given by DumpsterDoofus. I just want to point out the importance of having meaningful values for the parameters and initial conditions. I also want to show one might explore the parameter space with a Manipulate expression.

One of the phenomena demonstrated by the Lotka-Volterra model is that, under certain conditions, the predator and prey populations are cyclic with a phase shift between them. Here is a demonstration of this effect.

lv = 
  With[{a = 1.5, b = 1., c = 3., d = 1.},
    With[{x0 = 10., y0 = 6.}, 
      NDSolve[
        {x'[t] == a x[t] - b x[t] y[t], y'[t] == -c y[t] + d x[t] y[t], 
         x[0] == x0, y[0] == y0}, 
        {x[t], y[t]}, {t, 0, 20}]]]

Plot[Evaluate[{x[t], y[t]} /. lv], {t, 0, 20}]

plot

The static formulation is easily transformed into a dynamic one.

Manipulate[
  lv = NDSolve[
    {x'[t] == a x[t] - b x[t] y[t], y'[t] == -c y[t] + d x[t] y[t], 
    x[0] == x0, y[0] == y0}, {x[t], y[t]}, {t, 0, 20}]; 
  Plot[Evaluate[{x[t], y[t]} /. lv], {t, 0, 20}, PlotRange -> All, AxesOrigin -> {0,0}],
  {{lv, {}}, None},
  {{x0, 10}, 5, 30, 1},
  {{y0, 6}, 2, 20, 1},
  {{a, 1.5}, 1., 4.},
  {{b, 1.}, 1., 4.},
  {{c, 3.}, 1., 4.},
  {{d, 1.}, .1, 2.}]

demo

Improve this answer
$\endgroup$
10
$\begingroup$

You forgot to specify initial conditions on $x$ and $y$, which makes plotting impossible. However, even when initial conditions are specified, DSolve seems to choke on some sort of inverse-function garbage, so I'll sidestep it by using NDSolve instead (I'll let others determine why the analytic solution doesn't work):

{a, b, c, d} = {1, 2, 3, 4};
{X, Y} = NDSolveValue[{x'[t] == a x[t] - b x[t] y[t], 
    y'[t] == -c y[t] + d x[t] y[t], x[0] == 1, y[0] == 1}, {x, y}, {t,
     0, 1}];
Plot[{X[t], Y[t]}, {t, 0, 1}]

enter image description here

Improve this answer
$\endgroup$
3
  • $\begingroup$ Hey, thanks so much for this reply! I had tried setting initial conditions on x and y before, but I didn't seem to help so I had been adding and removing them in and out. This looks like exactly what I need! However, I tried to type it up myself, and am getting an error "There are fewer dependent variables, {y[t]}, than equations, so the system is overdetermined." $\endgroup$
    – elykl33t
    Commented Feb 17, 2015 at 1:25
  • $\begingroup$ @elykl33t: You can just copy and paste the code in my answer and execute it (no need to manually type it). Restart the kernel to make sure you don't have any stray definitions floating around, as that's a very common mistake. $\endgroup$
    – DumpsterDoofus
    Commented Feb 17, 2015 at 1:30
  • $\begingroup$ that was it! I initially did copy paste, and only tried the typing when it didn't work. I forget about the kernel restart! Thanks! $\endgroup$
    – elykl33t
    Commented Feb 17, 2015 at 1:32
2
$\begingroup$

There is no close-form answer because the solution is implicit. An equation of the phase plot can be obtained, however:

$$\frac{dy} { dx} =\frac{ (- c y + d x y) }{ (a x - b x y)} $$

which can rewritten as

$$\left(\frac{a}{y} - b\right) dy = \left(d -\frac{ c}{x}\right) dx$$

Hence the implicit equation of the phase plot that you can draw with, say,

ContourPlot:

a Log[y] - b y - d x + c Log[x] = Cte.

NDSolve must be used instead of DSolve which cannot invert this equation, wherein time information is absent.

I recommend you study it by hand, beginning with the variations of a y -> Log[y] - b y and x -> - d x + c Log[x], showing that a value is usually reached twice (this explains why the curve is a closed one, a.k.a. 'potato' :) ).

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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