Numerical solution of the Lotka-Volterra predator-prey equation¶

In this notebook (which covers roughly the same ground as this video) we use scipy.integrate.solve_ivp() to solve the Lotka-Volterra equations: \begin{align*} \dot{x} &= ax-bxy \\ \dot{y} &= -cy+dxy. \end{align*} Here $a,b,c$ and $d$ are positive constants; we will take them all to be one by default. These equations can be taken as a model of the dynamics of a population of predators and a population of prey. For example, $x$ could be the number of fish in a lagoon, and $y$ could be the number of sharks. Over time, both sharks and fish are born, the sharks eat some fish, and some sharks die, especially if there are too few fish: the equations are a crude model of these processes.

The above equation can be written in terms of the vector $u=\begin{bmatrix} x\\ y\end{bmatrix}$ as $$ \frac{du}{dt} = \begin{bmatrix} \dot{x} \\ \dot{y} \end{bmatrix} = \begin{bmatrix} ax - bxy \\ -cy+dxy \end{bmatrix} = \begin{bmatrix} au_0 - bu_0u_1 \\ -cu_1+du_0u_1 \end{bmatrix}, $$ and that is how we will formulate it in Python.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

In the cell below we call solve_ivp() to solve the Lotka-Volterra equations for $0\leq t\leq 24$ with $x=2.9$ and $y=1.0$ at $t=0$. The return value is an object sol with various attributes. The most important of these is sol.y, which is a numpy array of shape (2,4000). The first row sol.y[0] contains the values of $x$ at equally spaced times between $t=0$ and $t=24$, and the second row sol.y[1] contains the corresponding values of $y$. We plot $x$ and $y$ against $t$.

In [2]:

def du_dt(t, u, a = 1, b = 1, c = 1, d = 1):
    """ Return du/dt for the Lotka-Volterra model. """
    x, y = u
    return np.array([a * x - b * x * y, - c * y + d * x * y])

t = np.linspace(0, 24, 4000)
u0 = [2.9, 1.0]

sol = solve_ivp(du_dt, [t[0], t[-1]], u0, t_eval = t, method = 'DOP853', atol = 1e-9, args = (1, 1, 1, 1))
fish = sol.y[0]
sharks = sol.y[1]
U = np.log(sharks) + np.log(fish) - sharks - fish

plt.figure(figsize = (12, 6))
plt.plot(t, fish, color = 'red', label = 'Fish (thousands)')
plt.plot(t, sharks, color = 'blue', label = 'Sharks (tens)')
plt.xlabel('Time (years)')
plt.ylabel('Population')
plt.legend(loc='upper right')

Out[2]:

<matplotlib.legend.Legend at 0x21bf0db1d10>

No description has been provided for this image

We now plot $y$ against $x$ instead of plotting $x$ and $y$ against $t$. We plot a number of different curves corresponding to different initial conditions at $t=0$.

In [3]:

plt.figure(figsize = (8, 8))
t = np.linspace(0, 10, 1000)
ics = np.arange(1.0, 3.0, 0.1)
for f0 in ics:
    sol = solve_ivp(du_dt, [t[0], t[-1]], [f0, 1.0], t_eval = t, method = 'DOP853', atol = 1e-9, args = (1, 1, 1, 1))
    plt.plot(sol.y[0], sol.y[1])
plt.xlabel('Fish (thousands)')
plt.ylabel('Sharks (tens)')

Out[3]:

Text(0, 0.5, 'Sharks (tens)')

In [ ]: