LOTKA_VOLTERRA

Overview

The LOTKA_VOLTERRA function numerically solves the Lotka–Volterra predator-prey equations, a foundational model in mathematical biology that describes the dynamics of two interacting species—one as predator and one as prey. Developed independently by Alfred J. Lotka (1925) and Vito Volterra (1926), this system of coupled first-order nonlinear ordinary differential equations remains widely used in ecology, economics, and systems dynamics.

The model is defined by the following equations:

\frac{dx}{dt} = \alpha x - \beta xy

\frac{dy}{dt} = \delta xy - \gamma y

where x represents the prey population density, y represents the predator population density, \alpha is the prey birth rate, \beta is the predation rate coefficient, \delta is the predator reproduction rate per prey consumed, and \gamma is the predator death rate. The model produces characteristic oscillatory dynamics: prey populations grow exponentially in the absence of predators, predators increase when prey is abundant, and both populations cycle through peaks and troughs over time.

This implementation uses SciPy’s solve_ivp function to perform the numerical integration. The function supports multiple integration methods including explicit Runge-Kutta schemes (RK45, RK23, DOP853) for non-stiff problems and implicit methods (Radau, BDF, LSODA) for stiff systems. The default RK45 method uses a fifth-order Runge-Kutta formula with fourth-order error control. For more details, see the SciPy integration documentation and the SciPy GitHub repository.

The classic predator-prey oscillations predicted by this model have been observed in natural populations, most notably in the historical Hudson’s Bay Company fur trading records showing coupled fluctuations between lynx and snowshoe hare populations. For a comprehensive mathematical treatment, see the Wikipedia article on Lotka–Volterra equations.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=LOTKA_VOLTERRA(prey_initial, predator_initial, alpha, beta, delta, gamma, t_start, t_end, timesteps, lv_method)
  • prey_initial (float, required): Initial number of prey.
  • predator_initial (float, required): Initial number of predators.
  • alpha (float, required): Prey birth rate.
  • beta (float, required): Predation rate coefficient.
  • delta (float, required): Predator reproduction rate per prey eaten.
  • gamma (float, required): Predator death rate.
  • t_start (float, required): Start time of integration.
  • t_end (float, required): End time of integration.
  • timesteps (int, optional, default: 10): Number of timesteps to solve for.
  • lv_method (str, optional, default: “RK45”): Integration method.

Returns (list[list]): 2D list with header [t, Prey, Predator], or error message string.

Example 1: Predator-prey oscillations over 20 time units

Inputs:

prey_initial predator_initial alpha beta delta gamma t_start t_end timesteps
40 9 0.1 0.02 0.01 0.1 0 20 10

Excel formula:

=LOTKA_VOLTERRA(40, 9, 0.1, 0.02, 0.01, 0.1, 0, 20, 10)

Expected output:

t Prey Predator
0 40 9
2 30.2896 15.0284
4 18.2074 19.9381
6 9.61768 21.3943
8 5.08375 20.1818
10 2.89689 17.849
12 1.82351 15.3012
14 1.26782 12.9134
16 0.964522 10.8081
18 0.793433 9.00443
20 0.697514 7.48223
Example 2: Lotka-Volterra integration using RK23 method

Inputs:

prey_initial predator_initial alpha beta delta gamma t_start t_end timesteps lv_method
40 9 0.1 0.02 0.01 0.1 0 20 10 RK23

Excel formula:

=LOTKA_VOLTERRA(40, 9, 0.1, 0.02, 0.01, 0.1, 0, 20, 10, "RK23")

Expected output:

t Prey Predator
0 40 9
2 30.2897 15.0284
4 18.2074 19.9381
6 9.61769 21.3943
8 5.08376 20.1818
10 2.89689 17.849
12 1.82351 15.3012
14 1.26782 12.9134
16 0.964525 10.8081
18 0.793433 9.00444
20 0.697513 7.48224
Example 3: Extended simulation showing population decline phase

Inputs:

prey_initial predator_initial alpha beta delta gamma t_start t_end timesteps
50 10 0.1 0.02 0.01 0.1 0 40 10

Excel formula:

=LOTKA_VOLTERRA(50, 10, 0.1, 0.02, 0.01, 0.1, 0, 40, 10)

Expected output:

t Prey Predator
0 50 10
4 16.5187 25.9759
8 3.07668 23.8713
12 0.893843 17.1202
16 0.424973 11.7589
20 0.290423 7.99164
24 0.25511 5.41464
28 0.26586 3.66711
32 0.311038 2.48638
36 0.393447 1.69009
40 0.524585 1.15367
Example 4: Different rate parameters showing faster dynamics

Inputs:

prey_initial predator_initial alpha beta delta gamma t_start t_end timesteps lv_method
30 5 0.15 0.025 0.015 0.12 0 10 10 RK45

Excel formula:

=LOTKA_VOLTERRA(30, 5, 0.15, 0.025, 0.015, 0.12, 0, 10, 10, "RK45")

Expected output:

t Prey Predator
0 30 5
1 30.0522 6.97046
2 28.4189 9.60727
3 24.9856 12.7449
4 20.2814 15.8892
5 15.3205 18.3993
6 11.008 19.8578
7 7.73165 20.243
8 5.43802 19.7961
9 3.89594 18.816
10 2.87194 17.5481

Python Code

Show Code 
from scipy.integrate import solve_ivp
import numpy as np

def lotka_volterra(prey_initial, predator_initial, alpha, beta, delta, gamma, t_start, t_end, timesteps=10, lv_method='RK45'):
    """
    Numerically solves the Lotka-Volterra predator-prey system of ordinary differential equations.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        prey_initial (float): Initial number of prey.
        predator_initial (float): Initial number of predators.
        alpha (float): Prey birth rate.
        beta (float): Predation rate coefficient.
        delta (float): Predator reproduction rate per prey eaten.
        gamma (float): Predator death rate.
        t_start (float): Start time of integration.
        t_end (float): End time of integration.
        timesteps (int, optional): Number of timesteps to solve for. Default is 10.
        lv_method (str, optional): Integration method. Valid options: RK45, RK23, DOP853, Radau, BDF, LSODA. Default is 'RK45'.

    Returns:
        list[list]: 2D list with header [t, Prey, Predator], or error message string.
    """
    try:
      prey0 = float(prey_initial)
      predator0 = float(predator_initial)
      a = float(alpha)
      b = float(beta)
      d = float(delta)
      g = float(gamma)
      t0 = float(t_start)
      t1 = float(t_end)
      ntp = int(timesteps)

      if t1 <= t0:
        return "Error: Invalid input: t_end must be greater than t_start."
      if ntp <= 0:
        return "Error: Invalid input: timesteps must be a positive integer."
      if prey0 < 0 or predator0 < 0:
        return "Error: Invalid input: initial populations must be non-negative."
      if a < 0 or b < 0 or d < 0 or g < 0:
        return "Error: Invalid input: rate parameters must be non-negative."

      allowed_methods = ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']
      if lv_method not in allowed_methods:
        return f"Error: Invalid input: lv_method must be one of {allowed_methods}."

      t_eval = np.linspace(t0, t1, ntp + 1)

      def lv_ode(t, y):
        prey, predator = y
        dprey_dt = a * prey - b * prey * predator
        dpredator_dt = d * prey * predator - g * predator
        return [dprey_dt, dpredator_dt]

      sol = solve_ivp(
        lv_ode,
        [t0, t1],
        [prey0, predator0],
        method=lv_method,
        t_eval=t_eval,
        vectorized=False,
        rtol=1e-6,
        atol=1e-8
      )
      if not sol.success:
        return f"Error: Integration failed: {sol.message}"

      result = [["t", "Prey", "Predator"]]
      for i in range(len(sol.t)):
        row = [float(sol.t[i]), float(sol.y[0][i]), float(sol.y[1][i])]
        if not all(np.isfinite(v) for v in row):
          return "Error: Invalid output: nan or inf encountered in results."
        result.append(row)
      return result
    except Exception as e:
      return f"Error: {str(e)}"

Online Calculator

Initial number of prey.
Initial number of predators.
Prey birth rate.
Predation rate coefficient.
Predator reproduction rate per prey eaten.
Predator death rate.
Start time of integration.
End time of integration.
Number of timesteps to solve for.
Integration method.