solver Module#

ODE solver wrapper for compartmental models.

This module wraps scipy.integrate.solve_ivp() with epidemic-specific features: population conservation checks, stiff detection, and fallback to Radau for stiff systems.

Functions#

episia.models.solver.solve_model(derivatives, y0, t_span, t_eval=None, method='RK45', rtol=1e-06, atol=1e-08, max_step=inf, check_conservation=True, conservation_tol=0.001)[source]#

Solve an epidemic ODE system.

Parameters:

  • derivatives (Callable[[float, ndarray], ndarray]) – f(t, y) → dy/dt callable.

  • y0 (ndarray) – Initial state vector.

  • t_span (Tuple[float, float]) – (t_start, t_end).

  • t_eval (ndarray | None) – Output time points. If None uses 1000 points.

  • method (str) – scipy method: ‘RK45’ (default), ‘RK23’, ‘DOP853’, ‘Radau’, ‘BDF’, ‘LSODA’.

  • rtol (float) – Relative tolerance.

  • atol (float) – Absolute tolerance.

  • max_step (float) – Maximum step size (useful for stiff systems).

  • check_conservation (bool) – Raise if total population drifts > tol.

  • conservation_tol (float) – Fractional tolerance for conservation check.

Returns:

(t, solution) where solution has shape (n_compartments, len(t)).

Raises:

RuntimeError – Solver failure or population not conserved.

Return type:

Tuple[ndarray, ndarray]

episia.models.solver.estimate_herd_immunity(r0)[source]#

Herd immunity threshold: h = 1 - 1/R₀.

Parameters:

r0 (float) – Basic reproduction number.

Returns:

Fraction of population that needs immunity.

Raises:

ValueError – r0 <= 0.

Return type:

float

episia.models.solver.doubling_time(beta, gamma)[source]#

Early exponential doubling time T₂ = ln(2) / (β - γ).

Valid only during the initial exponential growth phase (S ≈ N).

Parameters:

  • beta (float) – Transmission rate.

  • gamma (float) – Recovery rate.

Returns:

Doubling time in the same units as beta/gamma.

Raises:

ValueError – beta <= gamma (no growth).

Return type:

float

Examples#

Basic solving:

import numpy as np
from episia.models.solver import solve_model

def sir_derivatives(t, y):
    S, I, R = y
    N = 1000
    beta, gamma = 0.3, 0.1
    dS = -beta * S * I / N
    dI =  beta * S * I / N - gamma * I
    dR =  gamma * I
    return np.array([dS, dI, dR])

y0 = [990, 10, 0]
t, sol = solve_model(
    derivatives=sir_derivatives,
    y0=y0,
    t_span=(0, 200),
    t_eval=np.linspace(0, 200, 1000)
)

Herd immunity:

from episia.models.solver import estimate_herd_immunity

for r0 in [1.5, 2.5, 3.5]:
    hit = estimate_herd_immunity(r0)
    print(f"R₀={r0:.1f}: need {hit:.1%} immune")

Doubling time:

from episia.models.solver import doubling_time

# Early epidemic growth
t_double = doubling_time(beta=0.3, gamma=0.1)
print(f"Cases double every {t_double:.1f} days")