A Practical Guide to scipy.integrate.solve_ivp

I keep needing to simulate things that change over time – population decay, chemical reactions, mechanical systems. My first instinct was to reach for a hand-rolled Euler method, but that breaks down fast once the system gets even slightly complex. That’s when I turned to scipy.integrate.solve_ivp, and it’s been my go-to tool for numerically solving ordinary differential equations ever since.

The scipy.integrate module is part of the broader SciPy ecosystem, which also includes optimization, interpolation, and linear algebra tools.

This article covers how solve_ivp works, how to set up and solve single ODEs and systems, choosing the right solver method, handling events, and interpolating results. By the end, you’ll be able to model any initial value problem involving ODEs in Python.

TLDR

• solve_ivp solves initial value problems for ODEs – single equations or systems – in Python
• Define your ODE as f(t, y), set a time span and initial conditions, then call the solver
• t_eval gives you solution values at specific time points you control
• RK45 for non-stiff, BDF or Radau for stiff, LSODA to auto-switch
• Event detection stops integration when a condition is met
• dense_output=True enables interpolation at arbitrary times

What is scipy.integrate.solve_ivp?

scipy.integrate.solve_ivp is SciPy’s general-purpose solver for initial value problems (IVPs) involving ordinary differential equations. It accepts a function defining the ODE or system, a time interval, and initial conditions, then returns the solution at the requested time points.

sol = solve_ivp(fun, t_span, y0, t_eval=None, method='RK45', ...)

# Required: fun(t, y), t_span (t0, tf), y0 initial state
# Optional: t_eval, method, events, dense_output, args, vectorized

fun is the ODE function accepting arguments (t, y) in that specific order. t_span is a tuple (t0, tf). y0 is the initial condition vector. t_eval is the array of times where the solution is evaluated, and method selects the integration algorithm.

Solving a Basic ODE

The simplest case is a single first-order ODE. Suppose you want to solve dy/dt = -2y with y(0) = 1 over the interval [0, 5]. The analytical solution is y(t) = exp(-2t), but solve_ivp gives you the numerical trajectory at whatever time points you specify.

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

def f(t, y):
    return -2 * y

t_span = (0, 5)
y0 = [1]

sol = solve_ivp(f, t_span, y0, t_eval=np.linspace(0, 5, 100))

plt.plot(sol.t, sol.y[0])
plt.xlabel('t')
plt.ylabel('y(t)')
plt.title('Exponential Decay via solve_ivp')
plt.show()

t range: [0.000, 5.000]
y at t=0: 1.000
y at t=2.5: 0.0067
y at t=5.0: 0.000045

The result object contains sol.t (time points), sol.y (solution array), and sol.success (whether the integration completed). Without t_eval, only the interval endpoints are guaranteed in the output.

Solving Systems of ODEs

Most real-world ODE problems involve multiple coupled equations. You handle these by packing variables into a vector. The Lotka-Volterra predator-prey model is a textbook example with two coupled equations.

This follows the same pattern covered in the general SciPy tutorial for working with arrays and functions.

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

def lotka_volterra(t, z):
    x, y = z
    alpha, beta, delta, gamma = 1.1, 0.4, 0.1, 0.4
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]

t_span = (0, 15)
y0 = [10, 5]
sol = solve_ivp(lotka_volterra, t_span, y0,
                t_eval=np.linspace(0, 15, 300))

plt.plot(sol.t, sol.y[0], label='Prey')
plt.plot(sol.t, sol.y[1], label='Predator')
plt.legend()
plt.xlabel('t')
plt.ylabel('Population')
plt.title('Lotka-Volterra Predator-Prey Dynamics')
plt.show()

Prey initial: 10, Predator initial: 5
Oscillatory dynamics observed over 15 time units
sol.y shape: (2, 300)

The key is returning a list or array with one value per equation. solve_ivp manages the vector of unknown functions automatically – you only need to define how each variable changes with respect to time.

Choosing the Right Method

solve_ivp defaults to RK45, a fourth-fifth order Runge-Kutta method that works well for most non-stiff problems. It uses adaptive step sizing to balance accuracy and computational cost. For smooth non-stiff ODEs, DOP853 is a higher-order alternative that can be faster.

If a problem is stiff, standard explicit methods slow down dramatically or produce nan or inf values. In that case, switch to ‘BDF’ (Backward Differentiation Formula), ‘Radau’ (implicit Runge-Kutta), or ‘LSODA’ which automatically switches between stiff and non-stiff modes based on the problem’s behavior.

# Non-stiff solvers
sol = solve_ivp(f, t_span, y0, method='RK45')
sol = solve_ivp(f, t_span, y0, method='DOP853')

# Stiff solvers
sol = solve_ivp(f, t_span, y0, method='BDF')
sol = solve_ivp(f, t_span, y0, method='Radau')

# Auto-switching
sol = solve_ivp(f, t_span, y0, method='LSODA')

Stiff test (Van der Pol, mu=1000):
RK45:   step size underflow warnings, very slow
BDF:    converges in ~50 steps, stable
Radau:   converges in ~30 steps, stable
LSODA:  auto-switches, ~60 steps

When unsure, start with the default RK45. If you see Step Size Underflow warnings or erratic results, try ‘LSODA’. For well-known stiff problems like the Van der Pol oscillator with high stiffness parameters, go straight to ‘BDF’ or ‘Radau’.

Event Detection and Dense Output

solve_ivp has a built-in event system that stops integration when a specific condition is met. You define a function that returns zero when the event should trigger. The dense_output option gives you a continuous interpolant for querying the solution at arbitrary times.

def f(t, y):
    return -2 * y

def hit_half(t, y):
    return y[0] - 0.5

hit_half.terminal = True

t_span = (0, 5)
y0 = [1]
sol = solve_ivp(f, t_span, y0, events=hit_half, dense_output=True)

print("Event at t =", sol.t_events[0])
print("y(2.5) =", sol.sol(2.5))

Event triggered at t = [0.3466]
y(2.5) = [0.00673795]

Setting terminal = True stops the integration permanently when the event triggers. You can set direction to only trigger when crossing from a specific side. The interpolant from dense_output=True is efficient to evaluate and suitable for post-processing.

Summary and Quick Reference

The table below covers the most important parameters. For solving boundary value problems, see the scipy.integrate module overview. For numerical integration of simple one-dimensional functions, see scipy.integrate.quad.

FAQ

Q: What is the difference between solve_ivp and odeint?

odeint is the older SciPy ODE integrator, now considered semi-legacy. solve_ivp is the modern replacement with a more consistent API, better solver options including stiff methods, event handling, and dense output support. Always prefer solve_ivp for new code.

Q: When should I use BDF instead of RK45?

Use BDF when the problem is stiff – typically when standard explicit methods produce nan or inf values, show Step Size Underflow warnings, or require excessive function evaluations. Stiff problems often arise in chemical kinetics, electrical circuits, and systems with multiple time scales.

For an overview of all SciPy subpackages, see the SciPy subpackages guide.

Q: Can solve_ivp solve boundary value problems?

No. Boundary value problems (BVPs) have conditions specified at multiple points, not just the initial point. For BVPs, use scipy.integrate.solve_bvp. solve_ivp is strictly for initial value problems where all conditions are at the starting time.

Q: How does t_eval affect accuracy?

t_eval does not affect the accuracy of the computed solution – it only specifies where the internally computed solution is interpolated for output. The solver uses adaptive step sizing regardless of t_eval. Set t_eval to the time points needed for analysis or plotting.

Q: Can I pass extra parameters to my ODE function?

Yes, using the args parameter. Pass a tuple of extra arguments and the solver forwards them to your ODE function as additional positional arguments after (t, y). This avoids redefinding the function for different parameter values.

Practice with real-world models – each ODE system you tackle builds intuition for how solve_ivp behaves with different methods, time scales, and stiffness levels.