Solving ODEs with solve_ivp

SciPy · Practical · 12 min

Population models, chemical kinetics, epidemics and mechanics are all ODEs you integrate numerically.

An ordinary differential equation gives the rate of change of a state: dy/dt = f(t, y). integrate.solve_ivp(fun, t_span, y0) integrates it from t_span[0] to t_span[1] starting at y0.

  • fun(t, y) takes time first and the state vector, and returns the derivative vector of the same length. Higher-order equations are rewritten as first-order systems: for x'' = -x use the state [x, v] with derivative [v, -x].
  • t_eval=np.linspace(...) chooses output times; otherwise you get the solver's own steps. dense_output=True gives sol.sol(t) for any t.
  • args=(...) passes parameters to fun.

The result has t, y (shape (n_states, n_times), one row per variable), success and message.

Methods: the default "RK45" suits most non-stiff problems. If it crawls with tiny steps, the problem is likely stiff (very different time scales, as in chemical kinetics); switch to "LSODA" (switches automatically), "BDF" or "Radau".

Accuracy is controlled by rtol (default 1e-3) and atol (1e-6). The defaults are loose; tighten them for long simulations.

Events are functions of (t, y) whose zero crossings are recorded in sol.t_events; set event.terminal = True to stop there, for example when a ball hits the ground.

Example

import numpy as np
from scipy.integrate import solve_ivp

def falling(t, state, drag):          # state = [height, velocity]
    h, v = state
    return [v, -9.81 - drag * v]

def hit_ground(t, state, drag):
    return state[0]                   # event when height crosses zero
hit_ground.terminal = True
hit_ground.direction = -1

sol = solve_ivp(falling, (0, 10), [100.0, 0.0], args=(0.1,),
                events=hit_ground, rtol=1e-8, atol=1e-8)
print(sol.success, sol.y.shape[0])          # True 2
print(round(sol.t_events[0][0], 3))         # 4.882 (seconds until impact)
print(round(sol.y_events[0][0][1], 2))      # -37.89 (impact velocity, m/s)

Remember

  • fun(t, y) returns dy/dt; rewrite higher-order equations as systems
  • Switch to LSODA, BDF or Radau for stiff problems
  • Tighten rtol/atol and use events to stop at conditions

Watch out for

  • Writing the derivative function as fun(y, t), the old odeint order.

Check yourself

  1. RK45 takes millions of tiny steps on a chemical kinetics model. What is the likely fix?

Sign in to answer, see the explanations, and get this lesson back as a quick review later.

Official docs: docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html