PYTHON / SCIPY
Integration and differential equations
Use quad and the sampled-data rules for definite integrals, and solve_ivp for ODE initial-value problems, including tolerances, events and stiff methods.
What you will learn
- Read quad's return value as (integral, estimated absolute error), not just a number
- Convert an nth-order ODE into an n-component first-order system for solve_ivp
- Set accuracy with rtol/atol; t_eval only chooses where the answer is reported
- Stop an integration on a condition with a terminal event function plus direction
Understanding Integration and differential equations
SciPy separates the two things people mean by "integration". scipy.integrate.quad takes a Python callable and two limits and evaluates the integral itself, adaptively subdividing the interval until an internal error estimate falls below tolerance, which is why it returns a pair: the value and an estimated absolute error. If you only have samples, such as a measured signal on a fixed grid, you cannot ask for more evaluations, so you use a fixed rule instead: trapezoid, simpson, or cumulative_trapezoid for a running total. Their accuracy is decided by the spacing of your data, not by a tolerance you request.
solve_ivp solves initial-value problems written as a single explicit first-order system, y'(t) = f(t, y) with y(t0) = y0, where y is a vector. Higher-order equations go in by hand: for y'' = -y you choose the state [y, y'] and return [y', -y], keeping the component order identical to the order in y0. The solver then takes its own adaptive steps; t_eval only says where you want output, and those values come from a polynomial interpolant (dense output) built on the steps the solver already took. That means extra output points cost almost nothing and buy no accuracy, because rtol and atol are the accuracy knobs.
The default method, RK45, is an explicit Runge-Kutta pair and the right first choice for smooth, non-stiff problems. When a system mixes rates differing by orders of magnitude, stability forces an explicit method into microscopic steps and the run appears to hang; switching to method='Radau' or 'BDF' (implicit, optionally with an analytic jac) fixes that. Events are the other feature worth learning: hand solve_ivp a function whose zero crossing matters, mark it terminal, and the solver brackets the crossing between two steps and refines it with a root finder instead of you scanning the output grid.
import numpy as np
from scipy.integrate import quad, solve_ivp
# quad: adaptive definite integral of a callable; exact value here is 2
area, err = quad(np.sin, 0.0, np.pi)
print(f"quad area = {area:.12f} error estimate < 1e-10: {err < 1e-10}")
# solve_ivp: dy/dt = -0.5*y, y(0) = 2 -> y(t) = 2*exp(-0.5*t)
def decay(t, y):
return -0.5 * y
sol = solve_ivp(decay, (0.0, 10.0), [2.0],
t_eval=[0.0, 2.5, 5.0, 7.5, 10.0],
rtol=1e-10, atol=1e-12)
print("solver finished:", sol.success)
for t, y in zip(sol.t, sol.y[0]):
print(f"t={t:4.1f} y={y:.8f} exact={2 * np.exp(-0.5 * t):.8f}")solve_ivp integrates a first-order vector field using its own adaptive steps, so you supply f(t, y), y0 and tolerances rather than a step size or an output grid.
Worked examples
Second-order ODE as a two-component system
Turns y'' = -y into a first-order system and checks the result against cos(t).
import numpy as np
from scipy.integrate import solve_ivp
def sho(t, u):
y, v = u # state = [position, velocity]
return [v, -y] # dy/dt = v, dv/dt = -y
t_eval = np.array([0.0, np.pi / 3, 2 * np.pi / 3, np.pi])
sol = solve_ivp(sho, (0.0, np.pi), [1.0, 0.0], t_eval=t_eval,
rtol=1e-10, atol=1e-12)
print("sol.y shape:", sol.y.shape)
for t, y in zip(sol.t, sol.y[0]):
print(f"t/pi = {t / np.pi:.3f} y = {y:7.4f} cos t = {np.cos(t):7.4f}")
print("agrees to 1e-8:", bool(np.max(np.abs(sol.y[0] - np.cos(t_eval))) < 1e-8))Example explained
Line 1sho unpacks the state and returns the derivative of each component in exactly the order used in y0.
Line 2y0 = [1.0, 0.0] encodes y(0)=1 and y'(0)=0, whose closed form is cos(t).
Line 3sol.y has one row per state component and one column per requested time, hence shape (2, 4).
Line 4rtol=1e-10 with atol=1e-12 pushes RK45's error far below the 1e-8 check; loosening them widens the gap.
Integrating data you cannot re-evaluate
Compares trapezoid, simpson and a running cumulative integral on 11 fixed samples of x**3.
import numpy as np
from scipy.integrate import trapezoid, simpson, cumulative_trapezoid
x = np.linspace(0.0, 1.0, 11) # spacing h = 0.1, samples only
y = x ** 3 # exact integral over [0, 1] is 0.25
print(f"trapezoid: {trapezoid(y, x):.6f}")
print(f"simpson: {simpson(y, x=x):.6f}")
c = cumulative_trapezoid(y, x, initial=0.0)
print(f"running integral at x=0.5 and x=1.0: {c[5]:.6f} {c[-1]:.6f}")Example explained
Line 1trapezoid overshoots by h**2/12 * (f'(1) - f'(0)) = 0.0025, the leading error term of the trapezoid rule.
Line 2simpson uses the same 11 samples but is exact for cubics, so it returns 0.25.
Line 3cumulative_trapezoid with initial=0.0 returns an array as long as x, so c[5] is the area up to x=0.5 and c[-1] matches trapezoid.
Line 4No tolerance argument exists here: refining the answer means collecting more samples.
Stopping at a condition with a terminal event
Finds when a thrown ball returns to height zero using an event function instead of scanning output times.
import numpy as np
from scipy.integrate import solve_ivp
g = 9.81
def ball(t, s):
height, vel = s
return [vel, -g]
def hits_ground(t, s):
return s[0] # event fires where the height is zero
hits_ground.terminal = True
hits_ground.direction = -1 # only downward crossings
sol = solve_ivp(ball, (0.0, 100.0), [0.0, 20.0], events=hits_ground,
rtol=1e-10, atol=1e-12)
t_hit = sol.t_events[0][0]
print("status:", sol.status, "(1 means a terminal event stopped it)")
print(f"landing time = {t_hit:.6f} analytic 2*v0/g = {2 * 20.0 / g:.6f}")
print(f"last t in solution = {sol.t[-1]:.6f}")Example explained
Line 1hits_ground returns the height, and solve_ivp detects a sign change between two steps, then refines the root on the interpolant.
Line 2direction = -1 is essential here: the height is exactly 0 at t=0, so without it the terminal event fires immediately at t=0.
Line 3terminal = True makes the integration stop at the root, so sol.t[-1] equals the event time and sol.status becomes 1.
Line 4t_events is a list with one array per event function, hence the double indexing in sol.t_events[0][0].
Important notes
quad's second return value is an estimated absolute error; if SciPy emits IntegrationWarning that estimate is unreliable, so split the interval at kinks or singularities, or pass points=... for finite discontinuities.
odeint is the older LSODA interface with reversed f(y, t) order and no event support; prefer solve_ivp, using method='LSODA' if you want the same automatic stiff/non-stiff switching.
Common mistakes
Writing the right-hand side as f(y, t), the odeint argument order, and passing it to solve_ivp: it runs without error but integrates -0.5*t instead of -0.5*y, silently giving the wrong curve.
Passing a scalar initial condition, solve_ivp(f, (0, 10), 2.0), which raises ValueError: `y0` must be 1-dimensional; a one-equation system still needs [2.0].
Treating t_eval as an accuracy setting: adding 10000 output times changes nothing because they are interpolated from the same adaptive steps, and points outside t_span raise ValueError.
Try it yourself
Change, predict, then run
Solve dy/dt = 0.8*y*(1 - y/50) with y(0) = 1 on [0, 20] using solve_ivp with rtol=1e-8, printing y at t = 0, 5, 10, 15, 20 beside the closed form 50/(1 + 49*np.exp(-0.8*t)). Then add a terminal event that fires when y reaches 25 and compare the crossing time with np.log(49)/0.8.
Open the Python workspaceCheck your understanding
A solve_ivp call with rtol=1e-3 gives a visibly wrong solution. You change t_eval from 10 points to 10000 points and the values at the original 10 times are identical. Why?
- t_eval only selects where the internal solution is sampled; step size and error are controlled by rtol and atol
- The solver caches its result between calls, so the second run reused the first solution
- 10000 output points exceed the internal step limit, so the extra points were discarded
- RK45 always uses fixed steps of (tf - t0)/10, so output density cannot affect anything
Show answer
solve_ivp chooses its steps from a local error estimate compared against rtol/atol, then evaluates a dense-output interpolant at whatever t_eval you asked for, so output density cannot change accuracy; option 3 is tempting but wrong because RK45 in solve_ivp is adaptive, not fixed-step, and it never derives its step size from the number of output points.