Preface

This section studies some first order nonlinear ordinary differential equations describing the time evolution (or “motion”) of those hamiltonian systems provided with a first integral linking implicitly both variables to a motion constant. An application has been performed on the Lotka--Volterra predator-prey system, turning to a strongly nonlinear differential equation in the phase variables.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part III of the course APMA0340
Introduction to Linear Algebra with Mathematica

Glossary

Predator-Prey Equations

Some situations require more than one differential equation to model a particular situation. We might use a system of differential equations to model two interacting species, say where one species preys on the other. For example, we can model how the population of Canadian lynx (lynx Canadians) interacts with a the population of snowshoe hare (lepus americanis) (see https://www.youtube.com/watch?v=ZWucOrSOdCs).

A good historical data are known for the populations of the lynx and snowshoe hare from the Hudson Bay Company. This large Canadian retail company, which owns and operates a large number of retail stores in North America and Europe, including Galeria Kaufhof, Hudson's Bay, Home Outfitters, Lord & Taylor, and Saks Fifth Avenue, was originally founded in 1670 as a fur trading company. The Hudson Bay Company kept accurate records on the number of lynx pelts that were bought from trappers from 1821 to 1940. The company noticed that the number of pelts varied from year to year and that the number of lynx pelts reached a peak about every ten years. The ten year cycle for lynx can be best understood using a system of differential equations.

The primary prey for the Canadian lynx is the snowshoe hare. We will denote the population of hares by H(t) and the population of lynx by L(t), where t is the time measured in years. We will make the following assumptions for our predator-prey model.

We now have a system of differential equations that describe how the two populations interact,
\[ \dot{H} = a\,H -b\, H\, L ,\qquad \dot{L} = c\,H\,L -\delta\, L , \]
where we use the dot notation to represent the derivatives with respect to time t.
Alfred James Lotka
The system of differential equation used to model predator-prey interactions
\[ \begin{split} \dot{x} = x \left( a - \alpha y \right) , \\ \dot{y} = y \left( -\delta + bx \right) , \end{split} \]
is usually referred to as Lotka--Volterra model. Here x is the number of prey (for example, rabbits); y is the number of some predator (for example, foxes); \( \dot{x} \ \mbox{ and} \ \dot{y} \) represent the instantaneous growth rates of the two populations; t represents time; a, α, δ, and b are positive real parameters describing the interaction of the two species.

The Lotka--Volterra system of equations is an example of a Kolmogorov model, which is a more general framework that can model the dynamics of ecological systems with predator-prey interactions, competition, disease, and mutualism. The equations which model the struggle for existence of two species (prey and predators) bear the name of two scientists: Lotka (1880--1949) and Volterra (1860--1940). They lived in different countries, had distinct professional and life trajectories, but they are linked together by their interest and results in mathematical modeling.

The predator–prey model was initially proposed by Alfred J. Lotka in the theory of autocatalytic chemical reactions in 1910. Lotka was born in Lemberg, Austria-Hungary, but his parents immigrated to the US. In 1925, he utilized the equations to analyze predator-prey interactions. Lotka published almost a hundred articles on various themes in chemistry, physics, epidemiology or biology, about half of them being devoted to population issues. He also wrote six books.

Vito Volterra

The same set of equations was published in 1926 by Vito Volterra, a mathematician and physicist, who had become interested in mathematical biology because of the impact by the marine biologist Umberto D'Ancona (1896--1964). Vito Volterra was born in Ancona, then part of the Papal States, into a very poor Jewish family. He attended the University of Pisa, where he became professor of rational mechanics in 1883. His most famous work was done on integral equations. In 1892, he became professor of mechanics at the University of Turin and then, in 1900, professor of mathematical physics at the University of Rome La Sapienza. His daughter Luisa married Umberto D’Ancona.

The predator-prey system of equations was later extended by many researchers, including C. S. Holling, Arditi--Ginzburg model, Rosenzweig-McArthur model, and some others.

The critical points of the Lotka--Volterra system of equations are the solutions of the algebraic equations

\[ x \left( a - \alpha y \right) =0 , \qquad y \left( -\delta + bx \right) = 0, \]
namely (0,0) and (δ/b, 𝑎/α). Next we try to examine the local behavior of solutions near each critical point using linearization approach. Neglecting the nonlinear terms, we obtain the corresponding linear system near the origin:
\[ \frac{\text d}{{\text d}t} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{bmatrix} a & 0 \\ 0 & -\delta \end{bmatrix}\begin{pmatrix} x \\ y \end{pmatrix} . \]
The eigenvalues of this system of equations are λ1 = 𝑎 > 0 and λ2 = - δ < 0. So the origin is a saddle point, unstable. At another critical point, linearization yields
\[ \frac{\text d}{{\text d}t} \begin{pmatrix} x- \delta / b \\ y - a/\alpha \end{pmatrix} = \begin{bmatrix} 0 & -\frac{\alpha \delta}{b} \\ \frac{ab}{\alpha} & 0 \end{bmatrix}\begin{pmatrix} x- \delta / b \\ y- a/\alpha \end{pmatrix} . \]
The eigenvalues of the above system of equations are pure complex conjugate:
\[ \lambda = \pm {\bf j} \sqrt{a\delta} . \]
Therefore, the Lotka--Volterra system of equation has not hyperbolic equilibrium point (δ/b, 𝑎/α) and linearization is inconclusive. So we need to find another approach. It is not hard to verify that the energy function
\[ E = bx + \alpha y - \delta \,\ln x - a\,\ln y + \mbox{constant} \]
is a Lyapunov function for the predator-prey model. Indeed, its derivative with respect to the system of differential equation is
\[ \dot{E} = b\,\dot{x} + \alpha\,\dot{y} - \frac{\delta}{x}\,\dot{x} - \frac{a}{y}\,\dot{y} = \left( a - \alpha\,y \right) \left( b\,x - \delta \right) + \left( -\delta + b\,x \right) \left( -a + \alpha\,y \right) \equiv 0. \]
Therefore, the stationary point (δ/b, a/α) is a stable critical point.

We may try to find the general solution of the Lotka--Volterra system of equations. From both equations, we get

\[ \frac{\dot{x}}{x} = a- \alpha \,y \qquad \mbox{and} \qquad \frac{\dot{y}}{y} = - \delta + b\,x . \]
Multiplying the first relation by (-δ + bx) and the second one by (𝑎 - αy), we obtain
\[ \frac{\dot{x}}{x} \left( - \delta + b\,x \right) = \left( a- \alpha \,y \right) \left( - \delta + b\,x \right) = \frac{\dot{y}}{y} \left( a- \alpha \,y \right) . \]
So we separated variables
\[ \frac{\dot{x}}{x} \left( - \delta + b\,x \right) = \frac{\dot{y}}{y} \left( a- \alpha \,y \right) , \]
which upon integration yields the energy integral.
Example 1:
s = NDSolve[{H'[t] == 0.4*H[t] - 0.01*H[t]*L[t],
L'[t] == -0.3*L[t] + 0.005*H[t]*L[t], H[0] == 70, L[0] == 50}, {H, L}, {t, 0, 100}]
a = Plot[Evaluate[H[t] /. s], {t, 0, 100}, PlotStyle -> {Thick, Red},
PlotRange -> {20, 85}, PlotLabels -> Placed[{"L(t)"}, Above], PlotLegends -> {"L(t)"}]
b = Plot[Evaluate[L[t] /. s], {t, 0, 100}, PlotStyle -> {Thick, Blue},
PlotRange -> {20, 80}, PlotLabels -> Placed[{"H(t)"}, Below], PlotLegends -> {"H(t)"}]
Show[{a, b}, PlotRange -> {20, 85}]

Notice that the predator population, L, begins to grow and reaches a peak after the prey population, H reaches its peak. As the prey population declines, the predator population also declines. Once the predator population is smaller, the prey population has a chance to recover that the cycle begins again.

f[x_, y_] = {x - 2 x y, -.5 y + x y};
Print["X' = ", MatrixForm[f[x, y]]]
Out[4]=
X' = (x-2 x y
-0.5 y+x y )
Print["x' = ", f[x, y][[1]], ", y' = ", f[x, y][[2]]]
Out[5]=
x' = x-2 x y, y' = -0.5 y+x y

Show[VectorPlot[f[x, y], {x, 0, 2}, {y, 0, 2}], Frame -> True,
BaseStyle -> {FontFamily -> "Times", FontSize -> 14}]
Show[VectorPlot[
f[x, y]/(10^-8 + Norm[f[x, y]]), {x, 0, 2}, {y, 0, 2}],
Frame -> True, BaseStyle -> {FontFamily -> "Times", FontSize -> 14}]


p2 = Show[
ParametricPlot[
Table[{x[t], y[t]} /. sol[[j]], {j, Length[ic]}], {t, 0, T},
PlotRange -> {{0, 4}, {0, 3}}],
Table[Graphics[{Arrowheads[.03],
Arrow[{ic[[j]], ic[[j]] + .01 f[ic[[j, 1]], ic[[j, 2]] ] }]}],
{j, Length[ic]}],
AxesLabel -> {"x", "y"},
BaseStyle -> {FontFamily -> "Times", FontSize -> 14},
PlotLabel -> "Phase portrait"
]


Plot[{x[t] /. sol[[1]], y[t] /. sol[[1]]}, {t, 0, T},
PlotStyle -> {{Thickness[0.01], RGBColor[1, 0, 0]}, {Thickness[0.01],
RGBColor[0, 1, 0]}}, AxesOrigin -> {0, 0},
AxesLabel -> {"t", "x(t),y(t)"},
BaseStyle -> {FontFamily -> "Times", FontSize -> 14},
PlotLabel -> "Time courses: Red - x(t), Green - y(t)"]

we can place two plots sude by side:

a = 2.2; b = 0.03; c = 1.4; d = 0.02;
eq1 = r'[t] == a*r[t] - b*r[t]*f[t]
eq2 = f'[t] == -c*f[t] + d*f[t]*r[t]
sol = NDSolve[{eq1, eq2, r[0] == 150, f[0] == 10}, {r, f}, {t, 0, 30}]
gr1 = Plot[Evaluate[{r[t], f[t]} /. sol], {t, 0, 30}, PlotStyle -> {{Thick, Hue[0.7]}, { Dashing[{0.03}], Hue[1]}}, Ticks -> {{{0.01, 0}, {10, "time"}, 20}, {500, {700, "number"}, 1000}}, LabelStyle -> {FontFamily -> "Times", FontSize -> 16}, PlotLabel -> "solid blue: rabbits, dashed red: foxes"];
gr2 = ParametricPlot[Evaluate[{r[t], f[t]} /. sol], {t, 0, 30}, PlotStyle -> {Thick, Hue[0.3]}, Ticks -> {{{0.01`, 0}, 500, {700, "rabbits"}, 1000}, {500, {700, "foxes"}, 1000}}, LabelStyle -> {FontFamily -> "Times", FontSize -> 16}];
sp = StreamPlot[{2.2` r - 0.03 f r, -1.4 f + 0.02` f r}, {r, 0, 200}, {f, 0, 150}, StreamScale -> Large]
GraphicsRow[{gr1, gr2, sp}]
    ■
Example 2: Consider the system of autonomous equations

\[ \dot{x} = x^2 -3\,xy, \qquad \dot{y} = 2xy-y^2 . \]
First, we plot the direction field using StreamPLot command:


StreamPlot[{x^2 - 3*x*y, 2*x*y - y^2}, {x, -3, 3}, {y, -3, 3}]
Then we can use NDSolve function to solve the given system numerically subject to some initial conditions:


deq1 = x'[t] == x[t]^2 - 3*x[t]*y[t];
deq2 = y'[t] == 2*x[t]*y[t] - y[t]^2;
soln = NDSolve[{deq1, deq2, x[0] == 1, y[0] == 1}, {x[t], y[t]}, {t, -10, 10}];
Next, we plot the corresponding solutions curve with the command
x = soln[[1, 1, 2]]; y = soln[[1, 2, 2]]; ParametricPlot[{x, y}, {t, -10, 10}, PlotStyle -> Thick]
We can plot as many such solutions curves as we want, and then display them simultaneously. For example, the following command creates a list of solution curves corresponding to the initial conditions \( x(0) = n/10, \quad y(0) =n/10 \) for \( n=-4,-3,\ldots , 3,4. \)


curve = Table[n, {n, -4, 4}]; curvep = Table[n, {n, 0, 4}]; Do[Clear[x, y];
soln = NDSolve[{deq1, deq2, x[0] == n/10, y[0] == n/10}, {x[t], y[t]}, {t, -10, 10}];
x = soln[[1, 1, 2]];
y = soln[[1, 2, 2]];
curve[[n]] = ParametricPlot[{x, y}, {t, -10, 10}, PlotStyle -> Thick], {n, -4, 4}];
Show[curve[[-4]], curve[[-3]], curve[[-2]], curve[[-1]], curve[[1]], curve[[2]],
curve[[3]], curve[[4]], AspectRatio -> 1]
To display solutions in the first quadrant, we type:



Do[Clear[x, y]; soln = NDSolve[{deq1, deq2, x[0] == n/10, y[0] == n/10}, {x[t], y[t]}, {t, -10, 10}];
x = soln[[1, 1, 2]];
y = soln[[1, 2, 2]];
curvep[[n]] = ParametricPlot[{x, y}, {t, -10, 10}, PlotStyle -> Thick], {n, 0, 4}];
Show[curvep[[0]], curvep[[1]], curvep[[2]], curvep[[3]], curvep[[4]], AspectRatio -> 1]
    ■
Example 3: Consider the following competition of two species for a single nutrient:
\[ \begin{split} \dot{S} = \gamma\,S(t) \left[ 1 - S(t) /K \right] - \frac{m_1}{y_1} \,\frac{x_1 S(t)}{a_1 + S(t)} - \frac{m_2}{y_2} \,\frac{x_2 S(t)}{a_2 + S(t)} , \\ \dot{x}_1 = \frac{m_1 x_1 (t)\,S(t)}{a_1 + S(t)} - \delta_1 x_1 (t) , \\ \dot{x}_2 = \frac{m_2 x_2 (t)\,S(t)}{a_2 + S(t)} - \delta_2 x_1 (t) , \end{split} \]
under the assumptions
\[ S(0) = S_0 > 0, \quad x_1 (0) = x_{10} > 0, \quad x_2 (0) = x_{20} > 0. \]
Here xi is the number of the i-th predator at time t, S(t) is the number of the prey at time t, mi is the maximum growth (birth) rate of the i-th predator, δi is the death rate for the i-th predator, yi is the yield factor for the i-th predator, which is the prey density at which the functional response of the predator is half maximal. The parameters γ and K are the intrinsic rate of increase and the carrying capacity or the prey population, respectively.
Theorem: Let
  1. \( \displaystyle 0 < a_1 \delta_1 < K \left( m_1 - \delta_1 \right) \) and
  2. \( \displaystyle a_2 \delta_2 > K \left( m_2 - \delta_2 \right) \) or \( b_2 \le 1 . \)
If \( \displaystyle K < a_1 + \frac{2a_1 \delta_1}{m_1 - \delta_1} , \) then
\[ \begin{split} \lim_{t\to\infty} S(t) &= S^{\ast} = \frac{a_1 \delta_1}{m_1 - \delta_1} , \\ \lim_{t\to\infty} x_1 (t) &= x_1^{\ast} = \frac{\gamma \left( 1-S^{\ast} /K \right) \left( a_1 + S^{\ast} \right)}{m_1 / y_1} , \\ \lim_{t\to\infty} x_2 (t) &= 0. \end{split} \]
The coexistence of two predators is possible although the problem remains open analytically.     ■

 

Non Kolmogorov models

Though straightforward, the Lotka-Volterra model makes various assumptions regarding the environment and evolution of predator and prey populations, limiting its capabilities when applying it to real-world examples. Five key assumptions are as follows:
  1. The population of prey is sustainable, as they are able to find and have ample food at all times.
  2. The food supply of predators depends solely on the population of the prey, meaning that the predator eats nothing else but the prey.
  3. The rate of change of population is proportional to its size.
  4. The environment is static. It does not change in favor of one species over the other, nor does genetic adaptation play a role.
  5. The predators’ appetites are limitless.

Thus, we can see that the Lotka-Volterra model fails to impose restrictions on certain parameters, in turn failing to produce realistic projections for predator and prey population data from the field. Two main problems with the Lotka-Volterra model that deter accurate predictions are as follows:

  1. In the absence of predators (y = 0), the prey population grows exponentially; there is nothing to limit their growth.
  2. The rate at which predators remove prey is βxy. By dividing this term by y, we see that the rate an individual predator consumes prey is βx. This rate increases as x increases, without any limit for x. Thus, in relation to assumption 5, this implies that an individual predator can never get full.

In order to address the limitations of the Lotka-Volterra model, Rosenzweig and MacArthur proposed their own model in 1963, which more realistically portrays predator-prey interactions.

Example 4: Consider the following predator-prey system with Beddington--DeAngelis functional response (see complete mathematical analysis in the article of D.T.Dimitrov and H.V. Kojouharov, Complete mathematical analysis of predator–prey models with linear prey growth and Beddington–DeAngelis functional response, Applied Mathematics and Computation, 162, (2), 523--538, 2005)
\[ \dot{x} = x - \frac{a\,xy}{1+x+y} , \qquad \dot{y} = \frac{b\,xy}{1+x+y} -c\,y , \]
where the constants are a = 6, b = 7.5, and c = 5.

Mathematical analysis of the Beddington--DeAngelis system shows that there exist two equilibria (0,0) and \( \left( \frac{ac}{ab - b -ac} , \frac{b}{ab - b -ac} \right) = (4,1) , \) being globally stable in the interior of the first quadrant. The eigenvalues of the Jacobian matrix at the origin are \( \lambda_1 =1 \quad\mbox{and}\quad \lambda_2 =5 , \) and the eigenvalues of the Jacobian at the point (4,1) are given by \( - \frac{1}{12} \pm {\bf j} \, \frac{\sqrt{119}}{12} . \)

StreamDensityPlot[{x-x*y/(1+x+y), 2*x*y/(1+x+y) -3*y}, {x,-3,3}, {y,-3,3}, StreamColorFunction->"Rainbow"]
    ■
Example 5: The Rosenzweig-MacArthur model assumes that in the absence of the predator, the prey population grows according to the logistic growth equation, not exponentially:
\[ \frac{{\text d}x}{{\text d}t} = \alpha x \left( 1 - \frac{x}{K} \right) , \]
where K represents the carrying capacity of the population. The solution to this logistic growth equation is:
\[ x(t) = \frac{K x_0 e^{\alpha t}}{K + x_0 \left( e^{\alpha t} -1 \right)} \]
where \( \lim_{t\to\infty} x(t) = K . \)

The Rosenzweig-MacArthur model assumes that the rate at which an individual predator consumes prey has a maximum value:

\[ \delta \,\frac{x}{b+x} . \]
As x increasingly gets larger, this ratio approaches δ:
\[ \lim_{x\to\infty} \,\delta \,\frac{x}{b+x} = \delta . \]
The parameter b controls how rapidly an individual predator becomes full, placing a limit over their appetite, and δ is the maximum of this rate. Thus, the Rosenzweig--MacArthur model addresses key problems as a result of the assumptions made by the Lotka--Volterra model, rendering a more realistic project of real-world predator-prey interactions. Putting the equations together, it can be seen that the Rosenzweig-MacArthur model is as follows:
\[ \begin{split} \frac{{\text d}x}{{\text d}t} = \alpha x \left( 1 - \frac{x\beta y}{b+x} \right) \\ \frac{{\text d}y}{{\text d}t} = -y \left( y - \frac{\delta x}{b+x} \right) . \end{split} \tag{5.1} \]
However, this model is the result upon scaling of the general form. The general form of the Rosenzweig-MacArthur model can be written as follows
\[ \begin{split} \frac{{\text d}x}{{\text d}t} = r x \left( 1 - \frac{x}{K} \right) -y\,h(x) , \\ \frac{{\text d}y}{{\text d}t} = -y \left( -c + d\,h(x) \right) . \end{split} \tag{5.2} \]
where
\[ h(x) = \frac{ax}{b+x} . \]
    ■

 

  1. M. Rosenzweig and R. MacArthur, Graphical representation and stability con-ditions of predator-prey interaction, American Naturalist (1963) 97, 209-2. https://www.jstor.org/stable/2458702
  2. Barabas, Gyorgy, D'Andrea, Rafael, and Stump, Simon Maccracken, Chesson's coexistence theory, Ecological Monographs, 2018, https://doi.org/10.1002/ecm.1302
  3. Batiha, K., Numerical Solutions of the Multispecies Predator-Prey Model by Variational Iteration Method, Journal of Computer Science, 2007, Vol. 3 (7): 523-527, 2007
  4. Bayliss, A., Nepomnyashchy, A.A., Volpert, V.A., Mathematical modeling of cyclic population dynamics, Physica D: Nonlinear Phenomena, 2019, Volume 394, doi: 10.1016/j.physd.2019.01.010
  5. Dellal, M., Lakrib, M., Sari, T., The operating diagram of a model of two competitors in a chemostat with an external inhibitor, Mathematical Biosciences, 302, No 8, 2018, 27--45.
  6. Dimitrov, D.T. and Kojouharov, H.V., Complete mathematical analysis of predator–prey models with linear prey growth and Beddington–DeAngelis functional response, Applied Mathematics and Computation, 162, (2), 523--538, 2005.
  7. Goel, N.S., Maitra, S.C., and Montroll, E., On the Volterra and other nonlinear models of interacting populations, Reviews of Modern Physics, 1971, Vol. 43, pp. 231-- https://doi.org/10.1103/RevModPhys.43.231
  8. Hsu, S.B., Hubbell, S.P., and Paul Waltman, Competing predators, SIAM Journal on Applied Mathematics, 35, No 4, 1978, 617--625.
  9. Hsu, S.B., Hubbell, S.P., and Paul Waltman, A mathematical theory for single-nutrient competition in continuous cultures of micro-organisms, SIAM Journal on Applied Mathematics, Vol 32, No 2, 1977, 366--383.
  10. Hsu, S.B., Hubbell, S.P., and Paul Waltman, A contribution to the theory of competing predators, Ecological Monographs, Vol 48, No 3, 1978, 337--349.
  11. May, R.M., Leonard, W.J., Nonlinear Aspects of Competition Between Three Species, SIAM Journal on Applied Mathematics, Vol. 29, Issue 2, pp. https://doi.org/10.1137/0129022
  12. Molla, H., Rahman, M.S., Sarwardi, S., Dynamics of a Predator–Prey Model with Holling Type II Functional Response Incorporating a Prey Refuge Depending on Both the Species, International Journal of Nonlinear Sciences and Numerical Simulation, 2018, Vol. 20, No. 1, https://doi.org/10.1515/ijnsns-2017-0224
  13. Olek. S. An accurate solution to the multispecies Lotka-Volterra equations, SIAM Review (Society for Industrial and Applied Mathematics), 1994, Vol. 36(3) (1994), 480–488; https://doi.org/10.1137/1036104
  14. Ruan, J., Tan, Y., Zhang, C., A Modified Algorithm for Approximate Solutions of Lotka-Volterra systems, Procedia Engineering, 2011, Volume 15, 2011, Pages 1493-1497; https://doi.org/10.1016/j.proeng.2011.08.277
  15. Scarpello, G.M. and Ritelli, D., A new method for the explicit integration of Lotka--Volterra equations, 2003, Divulgaciones Matematicas, Vol. 11, No. 1, pp. 1--17.
  16. Seo, Gunog and Wolkowicz, Gail S. K., Sensitivity of the dynamics of the general Rosenzweig--MacArthur model to the mathematical form of the functional response: a bifurcation theory approach. Journal of Mathematical Biology, 76, No 7, 2018, 1873--1906.

 

Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions