# The Lotka-Volterra predator-prey model

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.

• 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. Commented Feb 17, 2015 at 0:18
• Be sure that all your constants have been assigned numerical values before plotting the solution. Commented Feb 17, 2015 at 0:19
• 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. Commented Feb 17, 2015 at 0:55
• It strongly looks as a homework. Commented Feb 16, 2017 at 19:11

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]


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}},
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}},
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}},
Epilog -> {Text[Style["r=2", 20], Scaled[{0.6, 0.9}]]}]},
ImageSize -> 500, Spacings -> -50]


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


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


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


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


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


• 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." Commented Feb 17, 2015 at 1:25
• @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. Commented Feb 17, 2015 at 1:30
• 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! Commented Feb 17, 2015 at 1:32

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' :) ).