LOTKA_VOLTERRA

Overview

The LOTKA_VOLTERRA function numerically solves the classic Lotka-Volterra predator-prey system of ordinary differential equations, modeling the interaction between two species: prey and predator. The equations describe how the populations of prey ( $x$ ) and predator ( $y$ ) change over time:

\begin{align*} \frac{dx}{dt} &= \alpha x - \beta x y \\ \frac{dy}{dt} &= \delta x y - \gamma y \end{align*}

where $\alpha$ is the prey birth rate, $\beta$ is the predation rate coefficient, $\delta$ is the predator reproduction rate per prey eaten, and $\gamma$ is the predator death rate. The function uses SciPy’s ODE solver (scipy.integrate.solve_ivp ) to integrate the system from t_start to t_end. Only the basic two-species model is exposed; extensions for additional species, harvesting, or stochastic effects are not supported. The integration method can be selected from several common solvers, but only the most used options are exposed. This example function is provided as-is without any representation of accuracy.

Usage

To use the function in Excel:


=LOTKA_VOLTERRA(prey_initial, predator_initial, alpha, beta, delta, gamma, t_start, t_end, [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 for integration.
t_end (float, required): End time for integration.
method (str, optional, default=‘RK45’): Integration method. One of RK45, RK23, DOP853, Radau, BDF, LSODA.

The function returns a 2D array (table) with columns: t, Prey, Predator, representing time and the two populations at each step. If the input is invalid or integration fails, a string error message is returned.

Examples

Example 1: Basic Case (Default Method)

Inputs:

prey_initial	predator_initial	alpha	beta	delta	gamma	t_start	t_end	method
40.0	9.0	0.1	0.02	0.01	0.1	0.0	20.0	RK45

Excel formula:


=LOTKA_VOLTERRA(40.0, 9.0, 0.1, 0.02, 0.01, 0.1, 0.0, 20.0)

Expected output:

t	Prey	Predator
0.000	40.000	9.000
20.000	0.698	7.482

Example 2: RK23 Method

Inputs:

prey_initial	predator_initial	alpha	beta	delta	gamma	t_start	t_end	method
40.0	9.0	0.1	0.02	0.01	0.1	0.0	20.0	RK23

Excel formula:


=LOTKA_VOLTERRA(40.0, 9.0, 0.1, 0.02, 0.01, 0.1, 0.0, 20.0, "RK23")

Expected output:

t	Prey	Predator
0.000	40.000	9.000
20.000	0.698	7.482

Example 3: Longer Time Span

Inputs:

prey_initial	predator_initial	alpha	beta	delta	gamma	t_start	t_end	method
50.0	10.0	0.1	0.02	0.01	0.1	0.0	40.0	RK45

Excel formula:


=LOTKA_VOLTERRA(50.0, 10.0, 0.1, 0.02, 0.01, 0.1, 0.0, 40.0)

Expected output:

t	Prey	Predator
0.000	50.000	10.000
40.000	0.525	1.154

Example 4: All Arguments Specified

Inputs:

prey_initial	predator_initial	alpha	beta	delta	gamma	t_start	t_end	method
30.0	5.0	0.15	0.025	0.015	0.12	0.0	10.0	RK45

Excel formula:


=LOTKA_VOLTERRA(30.0, 5.0, 0.15, 0.025, 0.015, 0.12, 0.0, 10.0, "RK45")

Expected output:

t	Prey	Predator
0.000	30.000	5.000
10.000	2.872	17.548

Python Code


from scipy.integrate import solve_ivp
from typing import List, Union
 
def lotka_volterra(
    prey_initial: float,
    predator_initial: float,
    alpha: float,
    beta: float,
    delta: float,
    gamma: float,
    t_start: float,
    t_end: float,
    method: str = 'RK45'
) -> Union[List[List[float]], str]:
    """
    Numerically solves the Lotka-Volterra predator-prey system of ordinary differential equations.
 
    Args:
        prey_initial: Initial number of prey.
        predator_initial: Initial number of predators.
        alpha: Prey birth rate.
        beta: Predation rate coefficient.
        delta: Predator reproduction rate per prey eaten.
        gamma: Predator death rate.
        t_start: Start time of integration.
        t_end: End time of integration.
        method: Integration method ('RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA'). Default is 'RK45'.
 
    Returns:
        2D list with header row: t, Prey, Predator. Each row contains time and population values, or an error message (str) if input is invalid.
 
    This example function is provided as-is without any representation of accuracy.
    """
    # 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)
    except Exception:
        return "Invalid input: all parameters must be numbers."
    if t1 <= t0:
        return "Invalid input: t_end must be greater than t_start."
    if prey0 < 0 or predator0 < 0:
        return "Invalid input: initial populations must be non-negative."
    if a < 0 or b < 0 or d < 0 or g < 0:
        return "Invalid input: rate parameters must be non-negative."
    # Validate method
    allowed_methods = ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']
    if method not in allowed_methods:
        return f"Invalid input: method must be one of {allowed_methods}."
    # 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]
    # Time points for output (100 points)
    t_eval = [t0 + (t1 - t0) * i / 99 for i in range(100)]
    try:
        sol = solve_ivp(
            lv_ode,
            [t0, t1],
            [prey0, predator0],
            method=method,
            t_eval=t_eval,
            vectorized=False,
            rtol=1e-6,
            atol=1e-8
        )
    except Exception as e:
        return f"solve_ivp error: {e}"
    if not sol.success:
        return f"Integration failed: {sol.message}"
    # Prepare output: header row + data rows
    result = [["t", "Prey", "Predator"]]
    for i in range(len(sol.t)):
        prey_val = sol.y[0][i]
        predator_val = sol.y[1][i]
        # Disallow nan/inf
        if any([
            isinstance(prey_val, float) and (prey_val != prey_val or prey_val == float('inf') or prey_val == float('-inf')),
            isinstance(predator_val, float) and (predator_val != predator_val or predator_val == float('inf') or predator_val == float('-inf'))
        ]):
            return "Invalid output: nan or inf encountered in results."
        result.append([sol.t[i], prey_val, predator_val])
    return result

Example Workbook

Link to Workbook