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.

Example 1: SIR epidemic spread over 10 days

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.11111 986.319 12.4396 1.24184
2.22222 981.761 15.4529 2.78566
3.33333 976.135 19.1636 4.7017
4.44444 969.208 23.7171 7.07526
5.55556 960.713 29.2776 10.0095
6.66667 950.348 36.0267 13.6257
7.77778 937.769 44.1636 18.0671
8.88889 922.609 53.8916 23.4991
10 904.494 65.3967 30.1088
Example 2: SIR integration using RK23 method

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.11111 986.319 12.4393 1.24167
2.22222 981.764 15.4515 2.78472
3.33333 976.14 19.1601 4.69945
4.44444 969.223 23.7073 7.07004
5.55556 960.74 29.2599 9.99991
6.66667 950.388 36.0007 13.6109
7.77778 937.836 44.122 18.0423
8.88889 922.7 53.8348 23.4655
10 904.614 65.3221 30.0638
Example 3: Full epidemic curve showing peak and decline

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.55556 960.713 29.2776 10.0095
11.1111 883.036 78.8615 38.1023
16.6667 719.319 174.251 106.429
22.2222 491.286 275.392 233.323
27.7778 299.172 302.787 398.041
33.3333 186.119 258.033 555.848
38.8889 127.446 191.642 680.912
44.4444 97.3906 132.333 770.276
50 81.2082 87.9965 830.795
Example 4: SIR model using BDF implicit method

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.11111 986.318 12.4398 1.24196
2.22222 981.761 15.4535 2.78584
3.33333 976.133 19.165 4.70242
4.44444 969.205 23.7182 7.07715
5.55556 960.71 29.2777 10.0122
6.66667 950.344 36.0259 13.6297
7.77778 937.77 44.1586 18.0717
8.88889 922.619 53.8785 23.5029
10 904.515 65.3749 30.1105

Python Code

Show 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"Error: Invalid input: All initial values, parameters, and timesteps must be numbers. Details: {type(e).__name__}"

    # Validate ranges
    if susceptible_init < 0:
        return "Error: Invalid input: s_initial must be non-negative."
    if infected_init < 0:
        return "Error: Invalid input: i_initial must be non-negative."
    if recovered_init < 0:
        return "Error: Invalid input: r_initial must be non-negative."
    if infection_rate < 0:
        return "Error: Invalid input: beta must be non-negative."
    if recovery_rate < 0:
        return "Error: Invalid input: gamma must be non-negative."
    if time_end <= time_start:
        return "Error: Invalid input: t_end must be greater than t_start."
    if num_timesteps <= 0:
        return "Error: 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"Error: 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"Error: Integration error: {type(e).__name__}: {str(e)}"

    if not solution.success:
        return f"Error: 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

Initial number of susceptible individuals
Initial number of infected individuals
Initial number of recovered individuals
Infection rate parameter (per contact per time)
Recovery rate parameter (per time)
Start time of integration
End time of integration
Number of timesteps to evaluate the solution
Integration method for solving the ODE system