Programming the SIR model with Euler’s method

We now combine the SIR equations, initial conditions, time arrays and forward Euler method to simulate a complete epidemic. The program repeatedly recalculates all three rates because the biological flows change with the epidemic state.

Biological scenario

A closed population contains 1,000 people. Initially, 990 are susceptible, 10 are infectious and none have recovered. We use \(\beta=0.3\) per day and \(\gamma=0.1\) per day.

We want to estimate the epidemic trajectory for 160 days and determine the infectious peak, its timing and the final compartment sizes.

Ideas brought together

LESSON 3SIR differential equations and the rate function

LESSON 4Valid initial conditions and population conservation

LESSON 5Time grids and NumPy storage arrays

LESSON 6Repeated forward Euler updates

The SIR equations

\[\frac{dS}{dt}=-\beta\frac{SI}{N},\]

\[\frac{dI}{dt}=\beta\frac{SI}{N}-\gamma I,\]

\[\frac{dR}{dt}=\gamma I.\]

At each time point, the current values of \(S\), \(I\) and \(R\) determine the infection and recovery flows. Therefore, the rates must be recalculated inside the loop.

The three Euler updates

\[S_{n+1}=S_n+\left(\frac{dS}{dt}\right)_n\Delta t\]

\[I_{n+1}=I_n+\left(\frac{dI}{dt}\right)_n\Delta t\]

\[R_{n+1}=R_n+\left(\frac{dR}{dt}\right)_n\Delta t\]

Essential rule: use one common current state

Calculate all three derivatives from \(S_n,I_n,R_n\) before storing any next value. Do not calculate \(I_{n+1}\) using the newly updated \(S_{n+1}\). That would no longer be the standard forward Euler method.

Step 1: Write a reusable rate function

def sir_rates(S, I, R, beta, gamma):
    N = S + I + R
    infection_flow = beta * S * I / N
    recovery_flow = gamma * I

    dS_dt = -infection_flow
    dI_dt = infection_flow - recovery_flow
    dR_dt = recovery_flow

    return dS_dt, dI_dt, dR_dt

The function receives one epidemic state and returns its three instantaneous rates. It does not update the compartments.

Step 2: Prepare time and storage

number_of_steps = int(round(end_time / dt))
time = np.linspace(0, end_time, number_of_steps + 1)

S = np.zeros(number_of_steps + 1)
I = np.zeros(number_of_steps + 1)
R = np.zeros(number_of_steps + 1)

S[0], I[0], R[0] = S0, I0, R0

np.zeros() creates fixed-size numerical arrays. The initial conditions are placed at index 0. The later zeros are storage positions and will be replaced by the loop.

Step 3: Apply Euler repeatedly

for n in range(number_of_steps):
    dS_dt, dI_dt, dR_dt = sir_rates(
        S[n], I[n], R[n], beta, gamma
    )

    S[n + 1] = S[n] + dS_dt * dt
    I[n + 1] = I[n] + dI_dt * dt
    R[n + 1] = R[n] + dR_dt * dt

The function call uses only values at index n. After all derivatives are returned, the three next values are stored at index n + 1.

Step 4: Extract epidemic outcomes

peak_index = np.argmax(I)
peak_infectious = I[peak_index]
peak_time = time[peak_index]

np.argmax(I) returns the index at which \(I\) is largest.
The same index retrieves the peak value from I.
It also retrieves the corresponding peak time from time.

Run the complete SIR simulation

Run the program to see the numerical summary, selected-day table and compartment graph. You can change \(\beta\), \(\gamma\) or dt.

Interactive PythonComplete SIR Euler simulation

import numpy as np
import matplotlib.pyplot as plt

def sir_rates(S, I, R, beta, gamma):
    """Return the deterministic SIR rates at one state."""
    
    N = S + I + R
    infection_flow = beta * S * I / N
    recovery_flow = gamma * I
    
    dS_dt = -infection_flow
    dI_dt = infection_flow - recovery_flow
    dR_dt = recovery_flow
    
    return dS_dt, dI_dt, dR_dt

# Initial conditions
S0 = 990
I0 = 10
R0 = 0
N = S0 + I0 + R0

# Model parameters
beta = 0.3
gamma = 0.1

# Numerical settings
end_time = 160
dt = 0.1

# Validate settings
if min(S0, I0, R0, beta, gamma) < 0:
    raise ValueError("Compartments and parameters cannot be negative.")
if N <= 0:
    raise ValueError("The population must be positive.")
if dt <= 0:
    raise ValueError("dt must be positive.")

# Prepare time and storage
number_of_steps = int(round(end_time / dt))
time = np.linspace(0, end_time, number_of_steps + 1)

S = np.zeros(number_of_steps + 1)
I = np.zeros(number_of_steps + 1)
R = np.zeros(number_of_steps + 1)

S[0], I[0], R[0] = S0, I0, R0

# Forward Euler simulation
for n in range(number_of_steps):
    dS_dt, dI_dt, dR_dt = sir_rates(
        S[n], I[n], R[n], beta, gamma
    )
    
    S[n + 1] = S[n] + dS_dt * dt
    I[n + 1] = I[n] + dI_dt * dt
    R[n + 1] = R[n] + dR_dt * dt
    
    if min(S[n + 1], I[n + 1], R[n + 1]) < -1e-10:
        raise ValueError("A negative compartment occurred. Reduce dt.")

# Calculate outcomes
peak_index = np.argmax(I)
peak_infectious = I[peak_index]
peak_time = time[peak_index]
final_size = R[-1]

population = S + I + R
conservation_error = np.max(np.abs(population - N))

# Display summary
print("SIR epidemic summary")
print("-------------------------------------")
print(f"Peak infectious:       {peak_infectious:.2f}")
print(f"Time of peak:          {peak_time:.1f} days")
print(f"Final recovered:       {final_size:.2f}")
print(f"Maximum total error:   {conservation_error:.2e}")

# Display selected days
print("\nDay       S         I         R")
print("-------------------------------------")
for day in [0, 20, 40, 60, 80, 120, 160]:
    index = int(round(day / dt))
    print(
        f"{day:3d}  {S[index]:8.2f}  "
        f"{I[index]:8.2f}  {R[index]:8.2f}"
    )

# Plot the epidemic
plt.figure(figsize=(8, 4.8))
plt.plot(time, S, label="Susceptible")
plt.plot(time, I, label="Infectious")
plt.plot(time, R, label="Recovered")
plt.axvline(
    peak_time,
    color="grey",
    linestyle=":",
    label=f"Peak at day {peak_time:.1f}"
)
plt.xlabel("Time (days)")
plt.ylabel("Number of people")
plt.title("Deterministic SIR model: forward Euler method")
plt.grid(alpha=0.25)
plt.legend()
plt.tight_layout()

Output

Run the code to see the result.

How the outcome is calculated

Code	Technical purpose	Biological purpose
`R[-1]`	Uses index −1 to select the final array entry.	Gives the recovered population at the end of the simulation.
`np.argmax(I)`	Returns the position of the maximum value.	Locates the epidemic peak.
`np.abs(population - N)`	Calculates the size of each deviation from \(N\).	Measures population-conservation error.
`np.max(...)`	Selects the largest deviation.	Reports the worst conservation error across all times.
`plt.plot(time, I)`	Plots paired horizontal and vertical array values.	Shows how the infectious population changes through time.
`plt.legend()`	Displays the labels supplied to plotted lines.	Identifies the three compartments and peak marker.

Expected epidemic pattern

Susceptibledecreases continuously

Infectiousrises, reaches one peak, then falls

Recoveredincreases continuously

Initially, infection occurs faster than recovery, so the infectious population grows. As susceptible people are depleted, the infection flow weakens. The infectious peak occurs when the infection and recovery flows are approximately equal. After that point, recovery exceeds new infection and \(I\) declines.

The deterministic result is one smooth population-level trajectory. It represents average model behaviour and does not describe random individual infection and recovery events.

Important numerical and model conditions

The initial compartments and parameters must be non-negative.
The total population must be positive.
\(\Delta t>0\) and should be sufficiently small.
All rates in one Euler step must be calculated from the same current state.
A negative compartment signals an invalid numerical result, often caused by a time step that is too large.
Population conservation should be monitored across the entire trajectory.
The final simulation time must be long enough to capture the outcome of interest.

Try changing transmission

Run the model with beta = 0.15 and then with beta = 0.4, keeping all other values unchanged.

Predict which model will have the larger infectious peak.
Compare the peak times.
Compare the final recovered populations.
Explain the differences using the infection flow \(\beta SI/N\).