Numerical Integration with scipy.integrate

Numerical integration is one of the most practically important areas of scientific computing. In many real-world problems, integrals cannot be evaluated in closed form, or the underlying function is known only through measured data rather than a formula. scipy.integrate provides a comprehensive suite of tools that address both situations: high-accuracy adaptive quadrature for callable functions, rule-based integration for tabular data, and powerful solvers for ordinary differential equations. Understanding how and when to use each tool — and how to interpret their outputs — is essential for anyone doing quantitative work in science, engineering, or data analysis.

Numerical Quadrature with quad and dblquad

The workhorse of one-dimensional numerical integration in SciPy is scipy.integrate.quad. It implements an adaptive Gaussian quadrature algorithm (specifically QUADPACK under the hood), meaning it automatically refines its sampling of the integrand in regions where the function varies rapidly. This makes it far more reliable than fixed-step rules for general analytic functions.

To use quad, you pass it three required arguments: the integrand as a Python callable, the lower limit, and the upper limit. It always returns a tuple of two values: the estimated integral and an estimate of the absolute error in that result. This built-in accuracy information is a key advantage — you are never left guessing how trustworthy the answer is.

import numpy as np
from scipy.integrate import quad

# Integrate the Gaussian function from 0 to 1
def gaussian(x):
    return np.exp(-x**2)

result, error = quad(gaussian, 0, 1)
print(f"Integral: {result:.10f}")
print(f"Estimated absolute error: {error:.2e}")
# Integral: 0.7468241328
# Estimated absolute error: 8.29e-15

The error estimate of roughly 10⁻¹⁵ confirms the result is accurate to essentially machine precision. In practice you should always inspect this second element; if it is unexpectedly large, the integration may be struggling with a singularity, a highly oscillatory integrand, or a very wide domain.

One extremely useful feature is support for infinite limits. By passing np.inf or -np.inf as a limit, quad internally applies a change of variables to map the infinite domain to a finite one, then integrates on that transformed domain. This means you can compute integrals like the total probability under a normal distribution without any manual tricks.

# Integrate exp(-x^2) from -inf to +inf (should equal sqrt(pi))
result, error = quad(gaussian, -np.inf, np.inf)
print(f"Result:    {result:.10f}")
print(f"sqrt(pi):  {np.sqrt(np.pi):.10f}")
# Result:    1.7724538509
# sqrt(pi):  1.7724538509

When your integrand depends on extra parameters beyond x, pass them through the args keyword argument as a tuple rather than using global variables, which keeps code clean and reusable.

def power_law(x, n):
    return x**n

result, error = quad(power_law, 0, 1, args=(3,))
print(result)  # 0.25, i.e., integral of x^3 from 0 to 1

Double integrals over two-dimensional domains are handled by scipy.integrate.dblquad. Its signature is slightly more involved because the limits of the inner integral can themselves be functions of the outer variable, enabling non-rectangular domains. For a rectangular domain the inner limits are simply constants, which you can pass as lambda functions returning those constants.

from scipy.integrate import dblquad

# Integrate f(x, y) = x * y^2 over [0,1] x [0,2]
def f(y, x):          # NOTE: inner variable y comes first
    return x * y**2

result, error = dblquad(f, 0, 1,          # x limits
                           0, 2)           # y limits (constants)
print(result)   # 2/3 * 1 = 0.6667...
# Exact: integral_0^1 x dx * integral_0^2 y^2 dy = 0.5 * 8/3 = 4/3
print(result)   # 1.3333...

A common source of confusion: in dblquad the integrand is written as f(y, x) — the inner variable of integration comes first. The outer limits (for x here) are scalars; the inner limits for y can be callables of x. For example, integrating over the triangle where 0 ≤ y ≤ x within [0,1] looks like:

result, error = dblquad(f, 0, 1,
                           0, lambda x: x)   # inner y goes from 0 to x

Integrating Tabular Data with simpson and trapezoid

Adaptive quadrature requires a callable function. But in many real situations — sensor readings, experimental measurements, simulation outputs sampled at fixed time steps — you have only a table of (x, y) pairs and no analytical formula. SciPy provides two classical rules for this case: the trapezoidal rule and Simpson's rule.

scipy.integrate.trapezoid approximates the integral by drawing straight lines between consecutive data points and summing the resulting trapezoid areas. For each interval [x_i, x_{i+1}] the contribution is (y_i + y_{i+1}) / 2 * (x_{i+1} - x_i). It is simple, robust, and works correctly for both equally and unequally spaced x values.

import numpy as np
from scipy.integrate import trapezoid

# Approximate integral of sin(x) from 0 to pi using 10 points
x = np.linspace(0, np.pi, 10)
y = np.sin(x)

result = trapezoid(y, x)
print(f"Trapezoidal result: {result:.6f}")  # ~1.979 (exact is 2.0)

The trapezoid rule introduces an error proportional to h² (where h is the step size), so doubling the number of points reduces the error by a factor of four. For smooth functions, however, you can do significantly better with the same number of points by using scipy.integrate.simpson.

Simpson's rule fits a quadratic polynomial through each trio of consecutive points, then integrates that parabola exactly. This reduces the error to order h⁴, a dramatic improvement for smooth data. The classical form of Simpson's rule requires an odd number of data points (an even number of intervals); SciPy handles even-numbered arrays gracefully by applying a corrected formula for the last interval.

from scipy.integrate import simpson

result_simp = simpson(y, x=x)
print(f"Simpson result:     {result_simp:.6f}")  # ~2.0001, much closer to 2.0

Both functions accept x as an optional keyword argument. If you omit it, SciPy assumes unit spacing (i.e., x = [0, 1, 2, ...]), which is occasionally useful but usually you should supply the actual x array to get physically meaningful results. Both methods also accept unequally spaced x values, which is important when working with real sensor data that may not arrive at perfectly regular intervals.

Consider a practical example: you have acceleration data from an accelerometer and want to find velocity by integration.

# Simulated accelerometer data at irregular time stamps
t = np.array([0.0, 0.1, 0.25, 0.4, 0.5, 0.7, 1.0])
a = np.array([0.0, 1.2, 2.3, 2.1, 1.8, 0.9, 0.0])  # m/s^2

velocity_change = trapezoid(a, t)
print(f"Velocity change (trap): {velocity_change:.4f} m/s")

velocity_change_s = simpson(a, x=t)
print(f"Velocity change (simp): {velocity_change_s:.4f} m/s")

These methods are indispensable in experimental science. You often cannot re-run an experiment to get more data points, so choosing the more accurate rule (Simpson's) costs nothing extra computationally and can meaningfully improve your result.

Introduction to Ordinary Differential Equations (ODEs)

An ordinary differential equation describes how one or more quantities change with respect to a single independent variable (typically time). The simplest and most fundamental form is the first-order ODE:

dy/dt = f(t, y)

This says: the rate of change of y at time t is determined by some function f that may depend on both the current time and the current value of y. To pin down a unique solution you also need an initial condition — the value of y at some starting time t₀. Together, the ODE and its initial condition form an initial value problem (IVP).

A classic example is exponential decay: dy/dt = -k·y, with y(0) = y₀. The analytical solution is y(t) = y₀·exp(-k·t), but for more complex ODEs such closed-form solutions do not exist, and we must rely on numerical methods.

Higher-order ODEs (those involving second or higher derivatives) must be rewritten as systems of first-order equations before they can be handed to a numerical solver. The key insight is that a second derivative can always be treated as a new variable. For example, Newton's second law for a damped harmonic oscillator is:

m·x'' + c·x' + k·x = 0

Dividing through by m and introducing the state variables y₀ = x (position) and y₁ = x' (velocity), we get the equivalent first-order system:

• dy₀/dt = y₁
• dy₁/dt = -(c/m)·y₁ - (k/m)·y₀

The initial condition is now a vector: [x(0), x'(0)]. This vectorization is the universal pattern for translating any physical problem into a form SciPy can solve.

Solving ODEs with solve_ivp

scipy.integrate.solve_ivp is the modern, recommended interface for solving IVPs in SciPy. Its core signature is:

from scipy.integrate import solve_ivp

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

• fun: the derivative function with signature f(t, y)
• t_span: a two-element tuple (t_start, t_end)
• y0: a list or array of initial values, one per state variable
• t_eval: optional array of times at which to record the solution

The returned object is a Bunch (a dict-like object). The solution is stored in:

• result.t: a 1-D array of time points
• result.y: a 2-D array of shape (n_variables, n_time_points), where each row is one state variable's trajectory

Here is a complete example solving the damped oscillator system defined above:

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

# System parameters
m, c, k = 1.0, 0.5, 4.0   # mass, damping, spring constant

def oscillator(t, y):
    dydt = [y[1],                              # dy0/dt = y1 (velocity)
            -(c/m)*y[1] - (k/m)*y[0]]         # dy1/dt = acceleration
    return dydt

y0 = [1.0, 0.0]            # start at x=1, v=0
t_span = (0, 20)
t_eval = np.linspace(0, 20, 500)

result = solve_ivp(oscillator, t_span, y0, t_eval=t_eval)

print(result.success)      # True if integration succeeded
plt.plot(result.t, result.y[0], label='Position')
plt.plot(result.t, result.y[1], label='Velocity')
plt.legend(); plt.xlabel('Time'); plt.show()

The method argument controls which algorithm is used. The default 'RK45' is an explicit Runge-Kutta method of order 4/5 with adaptive step size control — an excellent general-purpose choice. However, some ODEs are stiff: they contain dynamics on vastly different time scales simultaneously, forcing explicit methods to take extremely small steps to remain stable. For stiff problems, implicit methods are necessary:

• 'Radau': implicit Runge-Kutta of order 5, very robust for stiff problems
• 'BDF': Backward Differentiation Formula, a variable-order method well-suited to very stiff systems
• 'LSODA': automatically switches between non-stiff (Adams) and stiff (BDF) methods

# Stiff example: Robertson chemical kinetics
def robertson(t, y):
    return [-0.04*y[0] + 1e4*y[1]*y[2],
             0.04*y[0] - 1e4*y[1]*y[2] - 3e7*y[1]**2,
             3e7*y[1]**2]

result = solve_ivp(robertson, (0, 1e11), [1, 0, 0],
                   method='Radau', rtol=1e-6, atol=1e-9)
print(result.success)

Attempting to solve the Robertson problem with 'RK45' would be either extremely slow (requiring millions of tiny steps) or outright unstable. Recognizing stiffness and selecting an appropriate method is one of the most important practical skills in numerical ODE solving.

The dense_output=True option instructs the solver to build a continuous interpolant of the solution as it goes, not just sample at discrete steps. After the solve, you can call result.sol(t) with any array of time values — even ones not used during integration — to get smooth, interpolated state values. This is useful for finding events, making high-resolution plots, or evaluating the solution at arbitrary points without re-running the solver.

result = solve_ivp(oscillator, t_span, y0, dense_output=True)
t_fine = np.linspace(0, 20, 5000)
y_fine = result.sol(t_fine)       # shape: (2, 5000)
plt.plot(t_fine, y_fine[0])
plt.show()

Defining the Derivative Function for solve_ivp

Getting the derivative function right is the most critical step in solving any ODE numerically. The function must have the signature f(t, y) where t is a scalar (the current time) and y is a NumPy array whose length equals the number of state variables. The function must return an array-like of the same length, where each element is the time derivative of the corresponding state variable.

# Example: three coupled first-order ODEs
def lorenz(t, y, sigma=10, rho=28, beta=8/3):
    dydt = [sigma * (y[1] - y[0]),
            y[0] * (rho - y[2]) - y[1],
            y[0] * y[1] - beta * y[2]]
    return dydt

Notice that the Lorenz system above has parameters (sigma, rho, beta) in addition to the state. Because solve_ivp only calls your function with (t, y), there are two clean ways to pass extra parameters:

• Default arguments (as above): simple and readable when parameters are fixed.
• Closure: define an outer function that captures parameters as local variables and returns the derivative function. This is especially useful when parameters may vary or are computed at runtime.

def make_oscillator(m, c, k):
    def oscillator(t, y):
        return [y[1], -(c/m)*y[1] - (k/m)*y[0]]
    return oscillator

osc = make_oscillator(m=2.0, c=0.3, k=5.0)
result = solve_ivp(osc, (0, 30), [1.0, 0.0], t_eval=np.linspace(0, 30, 600))

You can also use a lambda for concise one-liners, though this becomes hard to read for anything more than trivial equations:

# Exponential decay: dy/dt = -0.3*y
result = solve_ivp(lambda t, y: [-0.3*y[0]], (0, 10), [1.0])

A crucial practical point: units must be consistent throughout the derivative function. If t is in seconds but you accidentally use a decay constant in units of per-minute, your results will be off by a factor of 60. Always write down the physical units of each state variable and each parameter before coding the derivative function, and verify them by dimensional analysis.

When computing derivatives, use y[0], y[1], etc. to index into the state vector. For large systems it is often clearer to unpack at the top of the function:

def pendulum(t, y, g=9.81, L=1.0):
    theta, omega = y           # unpack for clarity
    dtheta_dt = omega
    domega_dt = -(g/L) * np.sin(theta)
    return [dtheta_dt, domega_dt]

This style makes the mapping between code variables and physical quantities transparent and reduces the risk of index errors in systems with many state variables.

Interpreting and Visualizing ODE Solutions

Once solve_ivp returns, your first action should always be to check result.success. This boolean is True if the integrator reached t_end successfully, and False if it terminated early — for example because the step size became too small to make progress (often a sign of a singularity) or the maximum number of steps was exceeded. If result.success is False, result.message provides a human-readable explanation.

result = solve_ivp(oscillator, (0, 100), [1.0, 0.0])
if not result.success:
    print("Integration failed:", result.message)
else:
    print("Integration succeeded.")

The solution arrays result.t and result.y are straightforward to work with. Remember that result.y is indexed by row, so result.y[0] is the entire time series of the first state variable, result.y[1] is the second, and so on.

# Time series plot
plt.figure(figsize=(10, 4))
plt.plot(result.t, result.y[0], label='Position x(t)')
plt.plot(result.t, result.y[1], label='Velocity v(t)')
plt.xlabel('Time (s)')
plt.ylabel('State')
plt.legend()
plt.tight_layout()
plt.show()

Beyond time-series plots, phase portraits are one of the most powerful visualization tools for dynamical systems. A phase portrait plots one state variable against another (without showing time explicitly), revealing the qualitative structure of the system's behavior — whether trajectories spiral inward to a stable equilibrium, orbit a limit cycle, or diverge chaotically.

# Phase portrait: position vs velocity
plt.figure(figsize=(5, 5))
plt.plot(result.y[0], result.y[1])
plt.xlabel('Position x')
plt.ylabel('Velocity v')
plt.title('Phase Portrait of Damped Oscillator')
plt.grid(True)
plt.show()

For a damped oscillator, the phase portrait is an inward spiral converging to the origin, visually confirming that the system loses energy over time. For an undamped oscillator (conservative system) it would be a closed ellipse. For the Lorenz system, the famous butterfly attractor emerges when you plot three state variables against each other.

Controlling numerical accuracy is done through two tolerances:

• rtol (relative tolerance): controls the error relative to the current magnitude of each state variable. Default is 1e-3.
• atol (absolute tolerance): controls the error floor for state variables near zero. Default is 1e-6.

Reducing either tolerance forces the solver to take smaller steps and perform more function evaluations, increasing accuracy at the cost of computation time. A common practical guideline is to tighten tolerances until your result stops changing meaningfully, then use one step tighter to add safety margin.

# High-accuracy integration
result_fine = solve_ivp(oscillator, (0, 20), [1.0, 0.0],
                        t_eval=t_eval,
                        rtol=1e-8, atol=1e-10)

result_coarse = solve_ivp(oscillator, (0, 20), [1.0, 0.0],
                          t_eval=t_eval,
                          rtol=1e-3, atol=1e-6)

# Compare the two solutions to assess accuracy impact
max_diff = np.max(np.abs(result_fine.y[0] - result_coarse.y[0]))
print(f"Max difference between fine and coarse: {max_diff:.2e}")

A practical workflow for any ODE problem: start with the default tolerances and RK45; check result.success; plot time series and inspect for non-physical behavior (e.g., energy growing without bound in a conservative system); if accuracy seems insufficient or the problem is stiff, switch method or tighten tolerances; finally, verify results against known special cases or analytical solutions where available. Following this process systematically makes numerical ODE solving both reliable and efficient.