SIR

Overview

The SIR function solves the classic SIR (Susceptible-Infectious-Recovered) compartmental model, a foundational framework in mathematical epidemiology for simulating infectious disease dynamics. First proposed by Kermack and McKendrick in 1927, the SIR model partitions a population into three compartments and tracks how individuals transition between them over time.

The model is governed by a system of ordinary differential equations (ODEs):

\frac{dS}{dt} = -\frac{\beta S I}{N}, \quad \frac{dI}{dt} = \frac{\beta S I}{N} - \gamma I, \quad \frac{dR}{dt} = \gamma I

where S, I, and R represent the susceptible, infectious, and recovered populations respectively, N is the total population, \beta is the transmission rate (contacts per unit time multiplied by transmission probability), and \gamma is the recovery rate (the reciprocal of the average infectious period). The basic reproduction number R_0 = \beta / \gamma determines whether an epidemic will occur: if R_0 > 1, the infection spreads; otherwise, it dies out.

This implementation uses scipy.integrate.solve_ivp, SciPy’s general-purpose ODE solver, which provides 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 background on compartmental models in epidemiology, see the Wikipedia article on compartmental models.

The SIR model is appropriate for diseases where recovery confers lasting immunity, such as measles, mumps, and rubella. Variants include SEIR (adding an exposed compartment), SIS (no lasting immunity), and SIRV (vaccination dynamics). The function returns time-series data for all three compartments, enabling visualization of epidemic curves and analysis of outbreak timing and magnitude.

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

Excel Usage

=SIR(s_initial, i_initial, r_initial, beta, gamma, t_start, t_end, timesteps, solve_ivp_method)

• s_initial (float, required): Initial number of susceptible individuals
• i_initial (float, required): Initial number of infected individuals
• r_initial (float, required): Initial number of recovered individuals
• beta (float, required): Infection rate parameter (per contact per time)
• gamma (float, required): Recovery rate parameter (per time)
• 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 evaluate the solution
• solve_ivp_method (str, optional, default: “RK45”): Integration method for solving the ODE system

Returns (list[list]): 2D list with header [t, S, I, R], or error message string.

Examples

Example 1: Demo case 1

Inputs:

s_initial	i_initial	r_initial	beta	gamma	t_start	t_end	timesteps
990	10	0	0.3	0.1	0	10	10

Excel formula:

=SIR(990, 10, 0, 0.3, 0.1, 0, 10, 10)

Expected output:

t	S	I	R
0	990	10	0
1.111	986.3	12.44	1.242
2.222	981.8	15.45	2.786
3.333	976.1	19.16	4.702
4.444	969.2	23.72	7.075
5.556	960.7	29.28	10.01
6.667	950.3	36.03	13.63
7.778	937.8	44.16	18.07
8.889	922.6	53.89	23.5
10	904.5	65.4	30.11

Example 2: Demo case 2

Inputs:

s_initial	i_initial	r_initial	beta	gamma	t_start	t_end	timesteps	solve_ivp_method
990	10	0	0.3	0.1	0	10	10	RK23

Excel formula:

=SIR(990, 10, 0, 0.3, 0.1, 0, 10, 10, "RK23")

Expected output:

t	S	I	R
0	990	10	0
1.111	986.3	12.44	1.242
2.222	981.8	15.45	2.785
3.333	976.1	19.16	4.699
4.444	969.2	23.71	7.07
5.556	960.7	29.26	10
6.667	950.4	36	13.61
7.778	937.8	44.12	18.04
8.889	922.7	53.83	23.47
10	904.6	65.32	30.06

Example 3: Demo case 3

Inputs:

s_initial	i_initial	r_initial	beta	gamma	t_start	t_end	timesteps
990	10	0	0.3	0.1	0	50	10

Excel formula:

=SIR(990, 10, 0, 0.3, 0.1, 0, 50, 10)

Expected output:

t	S	I	R
0	990	10	0
5.556	960.7	29.28	10.01
11.11	883	78.86	38.1
16.67	719.3	174.3	106.4
22.22	491.3	275.4	233.3
27.78	299.2	302.8	398
33.33	186.1	258	555.8
38.89	127.4	191.6	680.9
44.44	97.39	132.3	770.3
50	81.21	87.99	830.8

Example 4: Demo case 4

Inputs:

s_initial	i_initial	r_initial	beta	gamma	t_start	t_end	timesteps	solve_ivp_method
990	10	0	0.3	0.1	0	10	10	BDF

Excel formula:

=SIR(990, 10, 0, 0.3, 0.1, 0, 10, 10, "BDF")

Expected output:

t	S	I	R
0	990	10	0
1.111	986.3	12.44	1.242
2.222	981.8	15.45	2.786
3.333	976.1	19.17	4.702
4.444	969.2	23.72	7.077
5.556	960.7	29.28	10.01
6.667	950.3	36.03	13.63
7.778	937.8	44.16	18.07
8.889	922.6	53.88	23.5
10	904.5	65.37	30.11

Python Code

from scipy.integrate import solve_ivp as scipy_solve_ivp
import numpy as np

def sir(s_initial, i_initial, r_initial, beta, gamma, t_start, t_end, timesteps=10, solve_ivp_method='RK45'):
    """
    Solves the SIR system of ordinary differential equations for infection dynamics using scipy.integrate.solve_ivp (see [scipy.integrate.solve_ivp](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html)).

    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:
        s_initial (float): Initial number of susceptible individuals
        i_initial (float): Initial number of infected individuals
        r_initial (float): Initial number of recovered individuals
        beta (float): Infection rate parameter (per contact per time)
        gamma (float): Recovery rate parameter (per time)
        t_start (float): Start time of integration
        t_end (float): End time of integration
        timesteps (int, optional): Number of timesteps to evaluate the solution Default is 10.
        solve_ivp_method (str, optional): Integration method for solving the ODE system Valid options: RK45, RK23, DOP853, Radau, BDF, LSODA. Default is 'RK45'.

    Returns:
        list[list]: 2D list with header [t, S, I, R], or error message string.
    """
    # Validate and convert inputs
    try:
        susceptible_init = float(s_initial)
        infected_init = float(i_initial)
        recovered_init = float(r_initial)
        infection_rate = float(beta)
        recovery_rate = float(gamma)
        time_start = float(t_start)
        time_end = float(t_end)
        num_timesteps = int(timesteps)
    except (TypeError, ValueError) as e:
        return f"Invalid input: All initial values, parameters, and timesteps must be numbers. Error: {type(e).__name__}"

    # Validate ranges
    if susceptible_init < 0:
        return "Invalid input: s_initial must be non-negative."
    if infected_init < 0:
        return "Invalid input: i_initial must be non-negative."
    if recovered_init < 0:
        return "Invalid input: r_initial must be non-negative."
    if infection_rate < 0:
        return "Invalid input: beta must be non-negative."
    if recovery_rate < 0:
        return "Invalid input: gamma must be non-negative."
    if time_end <= time_start:
        return "Invalid input: t_end must be greater than t_start."
    if num_timesteps <= 0:
        return "Invalid input: timesteps must be a positive integer."

    # Validate method
    valid_methods = ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']
    if solve_ivp_method not in valid_methods:
        return f"Invalid input: solve_ivp_method must be one of {', '.join(valid_methods)}."
    # Create time array for evaluation
    t_eval = np.linspace(time_start, time_end, num_timesteps)

    # Define SIR ODE system
    def sir_ode(t, y):
        """SIR model differential equations."""
        susceptible, infected, recovered = y
        total_population = susceptible + infected + recovered

        # Avoid division by zero
        if total_population == 0:
            return [0, 0, 0]

        force_of_infection = infection_rate * susceptible * infected / total_population

        dS_dt = -force_of_infection
        dI_dt = force_of_infection - recovery_rate * infected
        dR_dt = recovery_rate * infected

        return [dS_dt, dI_dt, dR_dt]

    # Solve the ODE system
    try:
        solution = scipy_solve_ivp(
            sir_ode,
            [time_start, time_end],
            [susceptible_init, infected_init, recovered_init],
            method=solve_ivp_method,
            t_eval=t_eval
        )
    except Exception as e:
        return f"Integration error: {type(e).__name__}: {str(e)}"

    if not solution.success:
        return f"Integration failed: {solution.message}"

    # Build result as 2D list with header row
    header = ["t", "S", "I", "R"]
    result = [header]

    for i in range(len(solution.t)):
        row = [
            float(solution.t[i]),
            float(solution.y[0][i]),
            float(solution.y[1][i]),
            float(solution.y[2][i])
        ]
        result.append(row)

    return result

Online Calculator