LOTKA_VOLTERRA

=LOTKA_VOLTERRA(prey_initial, predator_initial, alpha, beta, delta, gamma, t_start, t_end, [timesteps], [solve_ivp_method])

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

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

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

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

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, solve_ivp_method='RK45'):
    """
    Numerically solves the Lotka-Volterra predator-prey system of ordinary differential equations.
 
    Wraps: scipy.integrate.solve_ivp (SciPy ODE integration). See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html
 
    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.
        timesteps: Number of timesteps to solve for. Default is 10.
        solve_ivp_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)
        ntp = int(timesteps)
    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 ntp <= 0:
        return "Invalid input: timesteps must be a positive integer."
    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 solve_ivp_method
    allowed_methods = ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']
    if solve_ivp_method not in allowed_methods:
        return f"Invalid input: solve_ivp_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=solve_ivp_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)):
        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 "Invalid output: nan or inf encountered in results."
        result.append([t_val, prey_val, predator_val])
    return result