Solving SIR and SEIR with solve_ivp

Writing Euler’s method from scratch reveals how numerical approximation works. For practical modelling, SciPy provides tested differential-equation solvers. Here we use solve_ivp while keeping the biological equations visible.

Biological scenario

We want to solve deterministic SIR and SEIR epidemics in populations of 1,000 people using the same transmission and recovery parameters. The SEIR model additionally uses a five-day mean exposed period.

We will compare the infectious peaks and inspect how introducing the exposed stage changes epidemic timing.

What solve_ivp does

The name means “solve an initial-value problem.” An initial-value problem contains:

1 EQUATIONSThe derivative function

2 TIME INTERVALStarting and ending times

3 INITIAL STATECompartment values at the start

4 PARAMETERSBiological rate values

5 OUTPUT TIMESTimes at which results are requested

solve_ivp repeatedly evaluates the derivative function, chooses internal step sizes and returns an approximate solution.

Import the solver

from scipy.integrate import solve_ivp

scipy is a scientific-computing library.
integrate is its integration and differential-equation module.
solve_ivp is the function being imported.

The derivative function has a required form

The solver calls a function using the pattern function(t, y, ...).

def sir_rhs(t, y, beta, gamma):
    S, I, R = y
    N = S + I + R

    infection_flow = beta * S * I / N
    recovery_flow = gamma * I

    return [
        -infection_flow,
        infection_flow - recovery_flow,
        recovery_flow
    ]

Part	Meaning
`t`	The current time supplied by the solver. The basic SIR equations do not use it explicitly, but the required function signature includes it.
`y`	The current state array. For SIR, it contains \([S,I,R]\).
`S, I, R = y`	Unpacks the state array into meaningful compartment names.
`beta, gamma`	Additional biological parameters passed through args.
`return [...]`	Returns one derivative for every entry in y, in the same order.

Call solve_ivp

sir_solution = solve_ivp(
    sir_rhs,
    t_span=(0, 200),
    y0=[990, 10, 0],
    t_eval=reporting_times,
    args=(beta, gamma),
    method="RK45",
    rtol=1e-8,
    atol=1e-10
)

Argument	Definition and use
`sir_rhs`	The derivative function. It is passed without parentheses because the solver must call it repeatedly.
`t_span=(0, 200)`	The integration interval: start at day 0 and finish at day 200.
`y0=[990, 10, 0]`	The initial state in the same order expected by the derivative function.
`t_eval=reporting_times`	The times at which solution values should be returned.
`args=(beta, gamma)`	A tuple of extra arguments passed to the derivative function.
`method="RK45"`	An adaptive Runge–Kutta method. It is also the default method.
`rtol`	Relative error tolerance, scaled in relation to solution size.
`atol`	Absolute error tolerance, important when solution values are close to zero.

Output times are not necessarily solver steps

t_eval specifies where results are reported. It does not force the adaptive solver to use those intervals as its internal steps. The solver may take smaller or larger internal steps while controlling estimated error.

This differs from the manual Euler program, where dt was both the update step and the spacing of stored values.

Understand the returned solution

solution.tone-dimensional array of returned times

solution.ytwo-dimensional array of compartment results

solution.successBoolean indicating whether integration finished successfully

For SIR, solution.y has three rows:

S, I, R = sir_solution.y

Each row has one value for every entry of solution.t. For SEIR, the solution has four rows.

Run both models

Interactive PythonSIR and SEIR with solve_ivp

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def sir_rhs(t, y, beta, gamma):
    """Return SIR derivatives in the order S, I, R."""
    
    S, I, R = y
    N = S + I + R
    
    infection_flow = beta * S * I / N
    recovery_flow = gamma * I
    
    return [
        -infection_flow,
        infection_flow - recovery_flow,
        recovery_flow
    ]

def seir_rhs(t, y, beta, sigma, gamma):
    """Return SEIR derivatives in the order S, E, I, R."""
    
    S, E, I, R = y
    N = S + E + I + R
    
    infection_flow = beta * S * I / N
    progression_flow = sigma * E
    recovery_flow = gamma * I
    
    return [
        -infection_flow,
        infection_flow - progression_flow,
        progression_flow - recovery_flow,
        recovery_flow
    ]

# Shared parameters
beta = 0.3
gamma = 0.1
sigma = 0.2

# Time interval and requested output times
start_time = 0
end_time = 200
reporting_times = np.linspace(start_time, end_time, 2001)

# Solve SIR
sir_solution = solve_ivp(
    sir_rhs,
    t_span=(start_time, end_time),
    y0=[990, 10, 0],
    t_eval=reporting_times,
    args=(beta, gamma),
    method="RK45",
    rtol=1e-8,
    atol=1e-10
)

# Solve SEIR
seir_solution = solve_ivp(
    seir_rhs,
    t_span=(start_time, end_time),
    y0=[990, 5, 5, 0],
    t_eval=reporting_times,
    args=(beta, sigma, gamma),
    method="RK45",
    rtol=1e-8,
    atol=1e-10
)

# Stop if either solver failed
if not sir_solution.success:
    raise RuntimeError(sir_solution.message)
if not seir_solution.success:
    raise RuntimeError(seir_solution.message)

# Unpack results
S_sir, I_sir, R_sir = sir_solution.y
S_seir, E_seir, I_seir, R_seir = seir_solution.y

# Calculate peaks
sir_peak_index = np.argmax(I_sir)
seir_peak_index = np.argmax(I_seir)

sir_peak = I_sir[sir_peak_index]
seir_peak = I_seir[seir_peak_index]

sir_peak_time = sir_solution.t[sir_peak_index]
seir_peak_time = seir_solution.t[seir_peak_index]

# Population-conservation errors
sir_total = S_sir + I_sir + R_sir
seir_total = S_seir + E_seir + I_seir + R_seir

sir_error = np.max(np.abs(sir_total - sir_total[0]))
seir_error = np.max(np.abs(seir_total - seir_total[0]))

# Display solver information
print("Solver information")
print("--------------------------------------------")
print(f"SIR success:       {sir_solution.success}")
print(f"SIR evaluations:   {sir_solution.nfev}")
print(f"SEIR success:      {seir_solution.success}")
print(f"SEIR evaluations:  {seir_solution.nfev}")

# Display model outcomes
print("\nModel   Peak I   Peak day   Final R   Max total error")
print("-------------------------------------------------------")
print(
    f"SIR    {sir_peak:7.2f}   {sir_peak_time:8.1f}  "
    f"{R_sir[-1]:8.2f}   {sir_error:.2e}"
)
print(
    f"SEIR   {seir_peak:7.2f}   {seir_peak_time:8.1f}  "
    f"{R_seir[-1]:8.2f}   {seir_error:.2e}"
)

# Plot both solutions
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))

axes[0].plot(sir_solution.t, S_sir, label="Susceptible")
axes[0].plot(sir_solution.t, I_sir, label="Infectious")
axes[0].plot(sir_solution.t, R_sir, label="Recovered")
axes[0].set_title("SIR")
axes[0].set_xlabel("Time (days)")
axes[0].set_ylabel("Number of people")
axes[0].grid(alpha=0.25)
axes[0].legend()

axes[1].plot(seir_solution.t, S_seir, label="Susceptible")
axes[1].plot(seir_solution.t, E_seir, label="Exposed")
axes[1].plot(seir_solution.t, I_seir, label="Infectious")
axes[1].plot(seir_solution.t, R_seir, label="Recovered")
axes[1].set_title("SEIR")
axes[1].set_xlabel("Time (days)")
axes[1].grid(alpha=0.25)
axes[1].legend()

fig.suptitle("Adaptive numerical solutions from solve_ivp")
fig.tight_layout()

Output

Run the code to see the result.

Additional returned information

Attribute	Meaning
`solution.success`	True when the solver reaches the end successfully.
`solution.message`	A textual explanation of termination or failure.
`solution.nfev`	Number of derivative-function evaluations performed.
`solution.t`	Returned time values.
`solution.y`	Returned state values, organised as variables by times.

Euler and solve_ivp

Manual Euler method

One fixed update formula
Step size selected directly
Easy to inspect and learn
Accuracy requires careful step-size testing

solve_ivp

Established numerical methods
Adaptive internal step sizes
Error tolerances supplied directly
Solver success and diagnostic information returned

Using a library solver does not replace mathematical understanding. The modeller must still define the correct equations, initial conditions, parameters, time interval and biological interpretation.

Biological interpretation

The SIR model moves newly infected people directly into the infectious compartment. The SEIR model delays infectiousness through the exposed compartment. This usually shifts and reshapes the infectious curve.

The displayed comparison is illustrative rather than a controlled scientific comparison because the two models have different compartment structures and initial-state definitions.

Important conditions and common errors

The order of y0, unpacked variables and returned derivatives must match.
The derivative function must return the same number of values as the initial-state vector.
All requested t_eval values must lie inside t_span and be ordered in the integration direction.
Always inspect solution.success before interpreting results.
Tighter tolerances usually require more computation; they do not correct an incorrect model.
Population conservation and non-negativity should still be checked.

Try changing the tolerances

First run the models with rtol=1e-8 and atol=1e-10. Then use rtol=1e-4 and atol=1e-6.

Compare nfev for each run.
Compare the reported peak values and times.
Compare the population-conservation errors.
Explain the trade-off between requested accuracy and computation.