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 986.723 12.1717 1.10512
2 982.752 14.7987 2.44955
3 977.948 17.9689 4.08307
4 972.151 21.7845 6.06473
5 965.174 26.3605 8.46515
6 956.809 31.824 11.3667
7 946.823 38.3125 14.8644
8 934.956 45.9757 19.0685
9 920.937 54.9593 24.1038
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 986.724 12.1714 1.10501
2 982.753 14.7978 2.44891
3 977.953 17.9658 4.08101
4 972.162 21.7774 6.06073
5 965.193 26.3482 8.45872
6 956.844 31.8017 11.3545
7 946.87 38.2832 14.847
8 935.028 45.9306 19.0419
9 921.03 54.9005 24.069
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 965.174 26.3605 8.46515
10 904.499 65.3953 30.1054
15 778.311 141.556 80.1327
20 583.953 240.28 175.767
25 385.927 300.161 313.912
30 244.974 290.906 464.121
35 164.389 238.396 597.215
40 119.716 178.782 701.502
45 95.3084 127.205 777.487
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 986.723 12.1718 1.10523
2 982.751 14.7992 2.44972
3 977.946 17.9701 4.08357
4 972.148 21.7858 6.066
5 965.171 26.3611 8.46758
6 956.807 31.8234 11.37
7 946.82 38.3118 14.8684
8 934.958 45.9692 19.0732
9 920.948 54.9449 24.1074
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.
    """
    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)

      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."

      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)}."

      t_eval = np.linspace(time_start, time_end, num_timesteps + 1)

      def sir_ode(t, y):
        susceptible, infected, recovered = y
        total_population = susceptible + infected + recovered
        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]

      solution = scipy_solve_ivp(
        sir_ode,
        [time_start, time_end],
        [susceptible_init, infected_init, recovered_init],
        method=solve_ivp_method,
        t_eval=t_eval
      )
      if not solution.success:
        return f"Error: Integration failed: {solution.message}"

      result = [["t", "S", "I", "R"]]
      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])]
        if not all(np.isfinite(v) for v in row):
          return "Error: Invalid output: nan or inf detected."
        result.append(row)
      return result
    except Exception as e:
      return f"Error: {str(e)}"

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