Tutorials 
Numerical Integration with SciPy
Numerical integration computes the area under a curve — or the solution to a differential equation — when an exact analytical answer is unavailable or inconvenient. You'd use it to evaluate integrals of functions with no closed form, to integrate noisy experimental data, or to simulate a system that evolves over time according to known rules. SciPy's `integrate` module handles all three cases: `quad` for 1D integrals, `dblquad` for 2D integrals, and `solve_ivp` for ODEs. ### 1D Integration with quad `quad` computes a definite integral to high accuracy using adaptive Gaussian quadrature. It takes a function and integration limits, and returns both the result and an error estimate. It handles smooth functions, singularities, and even infinite integration bounds.
from scipy.integrate import quad
import numpy as np
# ∫ sin(x) dx from 0 to π = 2
result, error = quad(np.sin, 0, np.pi)
print(f"∫sin(x) dx from 0 to π = {result:.6f} (error estimate: {error:.2e})")
# Gaussian integral: ∫ exp(-x²) dx from -∞ to ∞ = √π
result2, error2 = quad(lambda x: np.exp(-x**2), -np.inf, np.inf)
print(f"∫exp(-x²) dx = {result2:.6f}")
print(f"√π = {np.sqrt(np.pi):.6f}")- `quad` returns a tuple `(result, error)` — always unpack both so you can check the accuracy. - Passing `np.inf` and `-np.inf` as limits handles improper integrals over infinite ranges — SciPy maps them to a finite domain internally. - The error estimate is an upper bound on the absolute error, typically much smaller in practice. ### Visualizing What quad Computes It helps to see the area that `quad` is computing. Here's an example that shades the region under a curve alongside the numerical result.
from scipy.integrate import quad
import numpy as np
import matplotlib.pyplot as plt
f = lambda x: np.exp(-0.5 * x**2) / np.sqrt(2 * np.pi) # standard normal PDF
a, b = -1.96, 1.96
area, _ = quad(f, a, b)
print(f"P(-1.96 ≤ X ≤ 1.96) = {area:.4f}") # ≈ 0.95
x = np.linspace(-4, 4, 400)
x_fill = np.linspace(a, b, 200)
plt.figure(figsize=(7, 4))
plt.plot(x, f(x), color="steelblue", linewidth=2)
plt.fill_between(x_fill, f(x_fill), alpha=0.3, color="steelblue",
label=f"Area = {area:.4f}")
plt.axvline(a, color="gray", linestyle="--", linewidth=1)
plt.axvline(b, color="gray", linestyle="--", linewidth=1)
plt.title("∫ Standard Normal PDF from −1.96 to 1.96")
plt.xlabel("x")
plt.ylabel("Density")
plt.legend()
plt.tight_layout()
plt.show()P(-1.96 ≤ X ≤ 1.96) = 0.9500
- The shaded area represents the probability that a standard normal variable falls within ±1.96 — the basis of 95% confidence intervals. - `fill_between(x_fill, f(x_fill), ...)` fills between the curve and y=0, which is exactly what `quad` computes. ### 2D Integration with dblquad `dblquad` computes a double integral ∫∫ f(y, x) dy dx over a rectangular (or variable-boundary) region. Note the argument order: the inner integrand function receives `(y, x)`, and the bounds are given as `x_low, x_high, y_low_func, y_high_func`.
from scipy.integrate import dblquad
import numpy as np
# ∫₀¹ ∫₀¹ x·y dy dx = 0.25
result, error = dblquad(lambda y, x: x * y, 0, 1, 0, 1)
print(f"∫∫ x·y dy dx over [0,1]×[0,1] = {result:.4f} (error: {error:.2e})")
# Integrate over a triangular region: 0 ≤ x ≤ 1, 0 ≤ y ≤ x
result2, error2 = dblquad(lambda y, x: x + y, 0, 1, 0, lambda x: x)
print(f"∫∫ (x+y) dy dx over triangle = {result2:.4f} (error: {error2:.2e})")
- In `dblquad`, the function signature is `f(y, x)` — inner variable first, outer second. This is the opposite of what you might expect.
- `y_low=0, y_high=lambda x: x` makes the y upper limit a function of x, integrating over the triangular region where 0 ≤ y ≤ x.
- For higher-dimensional integrals, use `scipy.integrate.nquad`.
### Solving ODEs with solve_ivp
An ordinary differential equation (ODE) describes how a quantity changes over time: dy/dt = f(t, y). `solve_ivp` ("solve initial value problem") integrates the equation forward from a known starting point, using an adaptive Runge-Kutta method by default.
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
# Predator-prey (Lotka-Volterra) model:
# dx/dt = αx − βxy (prey grows, gets eaten)
# dy/dt = δxy − γy (predators eat prey, die off)
alpha, beta, delta, gamma = 1.0, 0.1, 0.075, 1.5
def lotka_volterra(t, state):
x, y = state
return [alpha*x - beta*x*y, delta*x*y - gamma*y]
sol = solve_ivp(lotka_volterra, t_span=[0, 40], y0=[10, 5],
dense_output=True, max_step=0.1)
t = np.linspace(0, 40, 1000)
y = sol.sol(t)
plt.figure(figsize=(9, 4))
plt.plot(t, y[0], color="steelblue", linewidth=2, label="Prey (x)")
plt.plot(t, y[1], color="crimson", linewidth=2, label="Predators (y)")
plt.title("Lotka-Volterra Predator-Prey Model")
plt.xlabel("Time")
plt.ylabel("Population")
plt.legend()
plt.tight_layout()
plt.show()- `t_span=[0, 40]` is the integration window; `y0=[10, 5]` sets the initial populations of prey and predators. - `dense_output=True` enables `sol.sol(t)` — a continuous interpolant you can evaluate at any time point, not just the solver's internal steps. - `max_step=0.1` limits the step size to prevent the adaptive solver from taking steps so large it misses rapid oscillations. `scipy.integrate` covers everything from a single number to a full dynamical system simulation. For finding where a function equals zero, see [root finding with SciPy](/tutorials/root-finding-with-scipy), and for optimization, see [SciPy optimize basics](/tutorials/scipy-optimize-basics).