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.22222 28.9261 15.6977
4.44444 15.8964 20.5712
6.66667 7.73692 21.2014
8.88889 3.91387 19.2132
11.1111 2.21194 16.4277
13.3333 1.41609 13.681
15.5556 1.01761 11.2497
17.7778 0.808012 9.19041
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.22222 28.9261 15.6978
4.44444 15.8964 20.5712
6.66667 7.73693 21.2014
8.88889 3.91388 19.2132
11.1111 2.21195 16.4277
13.3333 1.41609 13.681
15.5556 1.01761 11.2497
17.7778 0.808013 9.19041
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.44444 13.6696 26.5662
8.88889 2.22936 22.3564
13.3333 0.664603 15.1378
17.7778 0.345661 9.91073
22.2222 0.263517 6.43845
26.6667 0.258283 4.17556
31.1111 0.298029 2.71017
35.5556 0.382165 1.76385
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.11111 29.9605 7.23091
2.22222 27.8038 10.2737
3.33333 23.5138 13.8244
4.44444 18.0448 17.1168
5.55556 12.8068 19.3494
6.66667 8.70582 20.2219
7.77778 5.87382 19.9516
8.88889 4.03785 18.9429
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.
    """
    # Validate input types
    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)
    except Exception:
        return "Error: Invalid input: all parameters must be numbers."
    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."
    # Validate lv_method
    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}."
    # Create time array for evaluation
    t_eval = np.linspace(t0, t1, ntp)

    # Lotka-Volterra ODE system
    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]
    # Integrate
    try:
        sol = solve_ivp(
            lv_ode,
            [t0, t1],
            [prey0, predator0],
            method=lv_method,
            t_eval=t_eval,
            vectorized=False,
            rtol=1e-6,
            atol=1e-8
        )
    except Exception as e:
        return f"Error: solve_ivp error: {e}"
    if not sol.success:
        return f"Error: Integration failed: {sol.message}"
    # Prepare output: header row + data rows
    result = [["t", "Prey", "Predator"]]
    for i in range(len(sol.t)):
        t_val = float(sol.t[i])
        prey_val = float(sol.y[0][i])
        predator_val = float(sol.y[1][i])
        # Disallow nan/inf
        if any([
            t_val != t_val or t_val == float('inf') or t_val == float('-inf'),
            prey_val != prey_val or prey_val == float('inf') or prey_val == float('-inf'),
            predator_val != predator_val or predator_val == float('inf') or predator_val == float('-inf')
        ]):
            return "Error: Invalid output: nan or inf encountered in results."
        result.append([t_val, prey_val, predator_val])
    return result

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.