Numerics/Integration of Functions

Philosophy

Integration of functions is one of the three pillars of computational mathematics, alongside differentiation and root-finding. The central problem: given a function $ f(x) $ defined on $ [a,b] $, compute

$ I = \int_{a}^{b} f(x) \, dx $

with as few function evaluations as possible, to a prescribed accuracy. That's it. Everything else — adaptive quadrature, Gaussian rules, Monte Carlo, moments — is just strategy for getting the answer cheap and honest.

Opinionated stance: If you reach for scipy.integrate.quad without understanding what it does under the hood, you are doing it wrong. That function is an adaptive Gauss-Kronrod scheme — knowing why that matters is the difference between trusting your answer and being a passenger in someone else's code.

The Quadrature Mindset

All numerical integration reduces to a weighted sum:

$ I \approx \sum_{i=0}^{n} w_i \, f(x_i) $

where $ x_i $ are the abscissas (nodes) and $ w_i $ are the weights. Every method in this article — Newton-Cotes, Gaussian quadrature, Clenshaw-Curtis, Monte Carlo — is a particular choice of $ \{x_i, w_i\} $. The art is picking the right pair for your integrand.

A quadrature rule of degree $ d $ integrates all polynomials up to order $ d $ exactly. This is the benchmark. If your integrand is smooth (analytic), high-degree rules crush low-degree ones. If your integrand has kinks, corners, or singularities, no amount of polynomial exactness will save you — you need to handle the nasty part analytically, split the interval, or change variables.

Newton-Cotes: The Obvious (and Usually Wrong) Choice

Trapezoidal Rule (n=1)

The simplest closed Newton-Cotes formula. Fit a straight line through $ f(a) $ and $ f(b) $:

$ \int_a^b f(x) \, dx \approx \frac{b-a}{2} \left[ f(a) + f(b) \right] $

Error: $ -\frac{(b-a)^3}{12} f''(\xi) $ for some $ \xi \in [a,b] $.

Degree: 1 (exact for linear polynomials).

Simpson's Rule (n=2)

Fit a parabola through three equally spaced points:

$ \int_a^b f(x) \, dx \approx \frac{b-a}{6} \left[ f(a) + 4f\!\left(\frac{a+b}{2}\right) + f(b) \right] $

Error: $ -\frac{(b-a)^5}{2880} f^{(4)}(\xi) $.

Degree: 3. Simpson's rule is exact for cubics, which is why it overperforms relative to the naive expectation. This is a happy accident of symmetry — the $ x^3 $ error term cancels.

Simpson's 3/8 Rule (n=3)

$ \int_a^b f(x) \, dx \approx \frac{3h}{8} \left[ f(a) + 3f(a+h) + 3f(a+2h) + f(b) \right] $

where $ h = (b-a)/3 $. Degree: 3 (same as Simpson's, but with one extra function evaluation — almost never worth it).

Composite Rules

Applying a low-order rule over many subintervals is the workhorse approach:

$ \int_a^b f(x) \, dx = \sum_{k=0}^{N-1} \int_{x_k}^{x_{k+1}} f(x) \, dx $

Composite trapezoidal with $ N $ panels: error $ \mathcal{O}(h^2) $ where $ h = (b-a)/N $.

Composite Simpson's: error $ \mathcal{O}(h^4) $. This is the sweet spot for many engineering applications.

Why Newton-Cotes Is Usually Wrong

High-order Newton-Cotes rules (n > 8) suffer from Runge's phenomenon: the weights oscillate and include negative values, making them numerically unstable. The weights grow exponentially in magnitude, and cancellation destroys precision. If you need high accuracy, you don't use high-order Newton-Cotes — you use Gaussian quadrature.

Gaussian Quadrature: The Right Way

The Big Idea

Gaussian quadrature chooses both the nodes $ x_i $ and weights $ w_i $ to maximize the polynomial degree of exactness. With $ n $ points, a Gaussian rule has degree $ 2n - 1 $ — nearly double Newton-Cotes (which has degree $ n-1 $ for even $ n $).

How? The nodes are the roots of orthogonal polynomials. This is not a coincidence — it's a deep and beautiful connection between approximation theory and spectral methods.

Gauss-Legendre (The Default)

For integrals on $ [-1, 1] $ with unit weight function $ w(x) = 1 $:

$ \int_{-1}^{1} f(x) \, dx \approx \sum_{i=1}^{n} w_i f(x_i) $

where $ x_i $ are the roots of the Legendre polynomial $ P_n(x) $, and

$ w_i = \frac{2}{(1 - x_i^2) [P_n'(x_i)]^2} $

Python implementation (opinionated: use NumPy's polynomial tools, not tables):

import numpy as np
from numpy.polynomial.legendre import leggauss

def gauss_legendre(f, a, b, n):
    """
    Integrate f(x) from a to b using n-point Gauss-Legendre.

    This is the workhorse. For smooth integrands, n=16 is usually
    more than enough for machine precision on [0,1].
    """
    x_ref, w_ref = leggauss(n)
    # Map from [-1, 1] to [a, b]
    x_mapped = 0.5 * (b - a) * x_ref + 0.5 * (a + b)
    w_mapped = 0.5 * (b - a) * w_ref
    return np.dot(w_mapped, f(x_mapped))

Gauss-Kronrod: The Error Estimator You Actually Need

Plain Gaussian quadrature has no built-in error estimate. Gauss-Kronrod solves this by computing two rules simultaneously: an $ n $-point Gauss rule (degree $ 2n-1 $), and a $ (2n+1) $-point Kronrod extension that reuses the $ n $ Gauss nodes and adds $ n+1 $ new ones. The difference between them estimates the error.

The classic pair is G7/K15 (7 Gauss nodes, 15 Kronrod nodes). This is the engine inside scipy.integrate.quad.

import numpy as np

# Tabulated G7/K15 nodes and weights on [-1, 1]
# (Generated via high-precision root-finding; these are standard.)
GK15_NODES = np.array([
    -0.9914553711208126, -0.9491079123427585, -0.8648644233597691,
    -0.7415311855993945, -0.5860872354676911, -0.4058451513773972,
    -0.2077849550078985,  0.0,
     0.2077849550078985,  0.4058451513773972,  0.5860872354676911,
     0.7415311855993945,  0.8648644233597691,  0.9491079123427585,
     0.9914553711208126
])

GK15_WEIGHTS_KRONROD = np.array([
    0.022935322010529, 0.063092092629979, 0.104790010322250,
    0.140653259715525, 0.169004726639267, 0.190350578064785,
    0.204432940075298, 0.209482141084728,
    0.204432940075298, 0.190350578064785, 0.169004726639267,
    0.140653259715525, 0.104790010322250, 0.063092092629979,
    0.022935322010529
])

def gauss_kronrod_15(f, a, b):
    """15-point Gauss-Kronrod quadrature on [a, b], returns (estimate, error)."""
    x = 0.5 * (b - a) * GK15_NODES + 0.5 * (a + b)
    w = 0.5 * (b - a) * GK15_WEIGHTS_KRONROD
    estimate = np.dot(w, f(x))

    # Gauss 7-point rule for comparison (reuses subset of nodes)
    g7_idx = [1, 3, 5, 7, 9, 11, 13]
    g7_weights = np.array([0.129484966168870, 0.279705391489277,
                           0.381830050505119, 0.417959183673469,
                           0.381830050505119, 0.279705391489277,
                           0.129484966168870])
    g7_x = x[g7_idx]
    g7_w = 0.5 * (b - a) * g7_weights
    g7_est = np.dot(g7_w, f(g7_x))

    error = np.abs(estimate - g7_est)
    return estimate, error

Other Gaussian Flavors

Quadrature	Weight Function	Interval	Orthogonal Polynomial
Gauss-Legendre	$ w(x)=1 $	$ [-1,1] $	Legendre $ P_n $
Gauss-Chebyshev	$ w(x)=1/\sqrt{1-x^2} $	$ [-1,1] $	Chebyshev $ T_n $
Gauss-Laguerre	$ w(x)=e^{-x} $	$ [0,\infty) $	Laguerre $ L_n $
Gauss-Hermite	$ w(x)=e^{-x^2} $	$ (-\infty,\infty) $	Hermite $ H_n $
Gauss-Jacobi	$ w(x)=(1-x)^\alpha(1+x)^\beta $	$ [-1,1] $	Jacobi $ P_n^{(\alpha,\beta)} $

When to use each:

• Gauss-Laguerre — Integrals over $ [0,\infty) $ with exponential decay. Ideal for Laplace transform inversion and radial integrals in quantum mechanics.
• Gauss-Hermite — Integrals over $ (-\infty,\infty) $ with Gaussian decay. The natural choice for expectation values under a normal distribution.
• Gauss-Chebyshev — Integrands with endpoint singularities of type $ 1/\sqrt{1-x^2} $. Great for functions with square-root branch cuts.
• Gauss-Jacobi — Endpoint singularities of type $ (1-x)^\alpha(1+x)^\beta $. When you know the singularity structure, Jacabi rules are surgical.

Pro tip: You can always map $ [a,b] $ to $ [-1,1] $ with a linear change of variables: $ x = \frac{b-a}{2}t + \frac{a+b}{2} $

Adaptive Quadrature: Let the Computer Do the Work

Fixed-order quadrature is naive. Adaptive quadrature recursively subdivides the interval where the error estimate is large and coarsens where it's small. The algorithm:

1. Compute $ I_{[a,b]} $ with rule of order $ p $
2. Compute $ I_{[a,\frac{a+b}{2}]} + I_{[\frac{a+b}{2},b]} $ with the same rule
3. If $ |I_{[a,b]} - (I_{left} + I_{right})| < \varepsilon $, accept and return
4. Otherwise, recurse on each half

This is the algorithm behind scipy.integrate.quad (which uses adaptive Gauss-Kronrod with the G7/K15 pair).

def adaptive_quad(f, a, b, tol=1e-12, max_depth=50):
    """
    Adaptive quadrature using the 15-point Gauss-Kronrod rule.
    Recurses on intervals where error exceeds tolerance.
    """
    def _recurse(a, b, depth):
        I, err = gauss_kronrod_15(f, a, b)
        if err < tol * (b - a) or depth >= max_depth:
            return I, err

        mid = 0.5 * (a + b)
        I_left, err_left = _recurse(a, mid, depth + 1)
        I_right, err_right = _recurse(mid, b, depth + 1)
        return I_left + I_right, err_left + err_right

    return _recurse(a, b, 0)

Clenshaw-Curtis (The FFT-Based Alternative)

Clenshaw-Curtis quadrature uses Chebyshev nodes $ x_k = \cos(k\pi/N) $ (the extrema of $ T_N(x) $) and weights computed via FFT. With $ n+1 $ points, it has degree $ n $ — half that of Gauss-Legendre — but the weights are all positive and the nodes are nested (the $ n $-point set is a subset of the $ 2n $-point set), enabling efficient adaptive schemes.

For smooth, non-pathological integrands, Gauss-Legendre converges roughly twice as fast. But Clenshaw-Curtis is trivially adaptive (reuse old evaluations), and the FFT weight computation makes it $ \mathcal{O}(N \log N) $ to set up. In practice, the difference is often negligible — both converge exponentially for analytic integrands.

Opinion: Use Gauss-Legendre if you can afford to pick $ n $ up front. Use Clenshaw-Curtis if you're building an adaptive scheme and nested nodes matter.

Method of Moments

The Statistical View

"Method of moments" in integration refers to two distinct things. The first is statistical: estimate an integral by matching moments of an approximating distribution to sample moments. For integration, this means:

Given moments $ \mu_k = \int x^k f(x) \, dx $ for $ k = 0, 1, \dots, m $, approximate $ \int g(x) f(x) \, dx $ by constructing a quadrature rule that exactly integrates the first $ m $ monomials against $ f(x) $.

This is exactly what Gaussian quadrature does — the $ n $-point Gauss rule with weight $ f(x) $ is the unique rule that matches the first $ 2n $ moments. The nodes of Gaussian quadrature are the roots of the polynomial orthogonal with respect to the weight $ f(x) $.

The Quadrature Method of Moments (QMOM)

In engineering — especially reacting flow, aerosol dynamics, and population balance modeling — QMOM approximates a continuous distribution by a small number of delta functions:

$ n(x) \approx \sum_{i=1}^{N} w_i \, \delta(x - x_i) $

The $ 2N $ parameters $ \{w_i, x_i\} $ are determined by matching the first $ 2N $ moments:

$ m_k = \int x^k n(x) \, dx = \sum_{i=1}^{N} w_i \, x_i^k \qquad k = 0, 1, \dots, 2N-1 $

This is solved via the Product-Difference (PD) algorithm (Gordon, 1968) or the Wheeler algorithm, which computes the recurrence coefficients of the orthogonal polynomials from the moments, then finds the nodes and weights via Golub-Welsch (eigenvalues of the Jacobi matrix).

import numpy as np

def golub_welsch(alpha, beta):
    """
    Compute Gaussian quadrature nodes and weights from recurrence
    coefficients alpha (diagonal) and sqrt(beta) (off-diagonal)
    of the Jacobi matrix J. Nodes = eigenvalues of J.
    Weights = w_0 * (first component of eigenvectors)^2.
    """
    n = len(alpha)
    J = np.diag(alpha)
    for i in range(n - 1):
        J[i, i + 1] = np.sqrt(beta[i + 1])
        J[i + 1, i] = np.sqrt(beta[i + 1])

    eigenvalues, eigenvectors = np.linalg.eigh(J)
    weights = beta[0] * eigenvectors[0, :]**2
    # Sort by eigenvalue
    idx = np.argsort(eigenvalues)
    return eigenvalues[idx], weights[idx]

def moments_to_quadrature(moments, n_nodes):
    """
    Given the first 2*n_nodes moments, compute the n_nodes-point
    Gaussian quadrature rule via the Wheeler algorithm.
    moments[k] = ∫ x^k w(x) dx for k = 0..2n-1.
    """
    assert len(moments) >= 2 * n_nodes
    M = 2 * n_nodes
    # Modified Chebyshev algorithm / Wheeler
    alpha = np.zeros(n_nodes)
    beta = np.zeros(n_nodes)
    beta[0] = moments[0]

    sigma = np.zeros((M, M))
    sigma[0, :] = moments[:M]

    for k in range(1, M):
        for ell in range(k, M):
            sigma[k, ell] = sigma[k - 1, ell + 1] \
                          - alpha[k - 1] * sigma[k - 1, ell] \
                          - beta[k - 1] * sigma[k - 2, ell] \
                              if k > 1 else sigma[k - 1, ell + 1] \
                                         - alpha[k - 1] * sigma[k - 1, ell]

        if k < n_nodes:
            alpha[k] = sigma[k, k + 1] / sigma[k, k] - sigma[k - 1, k] / sigma[k - 1, k - 1]
            beta[k] = sigma[k, k] / sigma[k - 1, k - 1]

    return golub_welsch(alpha[:n_nodes], beta[:n_nodes])

Why QMOM Wins

For tracking particle size distributions in CFD, QMOM with $ N=3 $ or $ N=4 $ nodes (6-8 moments) often matches the accuracy of sectional methods with hundreds of bins. The cost: solving a small eigenvalue problem per cell per timestep, versus transporting hundreds of scalars. The tradeoff is heavily in QMOM's favor for reacting multiphase flows.

Monte Carlo Integration: When Dimension Kills You

All the methods above scale exponentially with dimension — a product rule with $ n $ points per axis in $ d $ dimensions uses $ n^d $ evaluations. For $ d > 4 $, you need Monte Carlo.

The Monte Carlo estimator:

$ I = \int_{\Omega} f(\mathbf{x}) \, d\mathbf{x} \approx \frac{V}{N} \sum_{i=1}^{N} f(\mathbf{x}_i), \quad \mathbf{x}_i \sim \text{Uniform}(\Omega) $

Error: $ \mathcal{O}(1/\sqrt{N}) $, independent of dimension. This is the killer feature. For $ d=1 $, MC is terrible compared to Gaussian quadrature. For $ d=10 $, MC is your only option.

Variance reduction techniques:

• Importance sampling: Draw from $ p(\mathbf{x}) \propto |f(\mathbf{x})| $ instead of uniform. Reweight by $ 1/p(\mathbf{x}) $.
• Stratified sampling: Partition $ \Omega $ into strata, sample each, combine. Guarantees variance ≤ crude MC.
• Quasi-Monte Carlo (QMC): Use low-discrepancy sequences (Sobol, Halton) instead of pseudorandom. Error drops to $ \mathcal{O}((\log N)^d / N) $, which for moderate $ d $ and large $ N $ beats $ 1/\sqrt{N} $ substantially.

import numpy as np
from scipy.stats import qmc  # Quasi-Monte Carlo

def monte_carlo_integrate(f, bounds, n=100_000, method='sobol'):
    """
    Integrate f over a d-dimensional box defined by bounds = [(a1,b1), ...].
    Uses Sobol QMC by default. Crude MC if method='random'.
    """
    dim = len(bounds)
    if method == 'sobol':
        sampler = qmc.Sobol(d=dim, scramble=True)
        points = sampler.random(n)
    else:
        points = np.random.uniform(size=(n, dim))

    # Map to bounds
    volume = 1.0
    for i, (lo, hi) in enumerate(bounds):
        points[:, i] = lo + (hi - lo) * points[:, i]
        volume *= (hi - lo)

    vals = np.array([f(p) for p in points])
    estimate = volume * np.mean(vals)
    stderr = volume * np.std(vals) / np.sqrt(n)
    return estimate, stderr

Dealing with Singularities

Subtract the Singularity

If $ f(x) = g(x) / \sqrt{x} $ on $ [0, 1] $, separate out the nasty part:

$ \int_0^1 \frac{g(x)}{\sqrt{x}} dx = \int_0^1 \frac{g(0)}{\sqrt{x}} dx + \int_0^1 \frac{g(x) - g(0)}{\sqrt{x}} dx = 2g(0) + \int_0^1 \frac{g(x) - g(0)}{\sqrt{x}} dx $

The first term is analytic; the second has a removable singularity (the numerator now vanishes at $ x=0 $) and can be integrated with standard methods.

Change of Variables

For endpoint singularities, a nonlinear change of variables can remove them entirely. The classic transformation for $ [0,1] $ with $ \sqrt{x} $ singularities:

$ x = t^2, \quad dx = 2t \, dt, \quad \int_0^1 \frac{g(x)}{\sqrt{x}} dx = 2 \int_0^1 g(t^2) \, dt $

The integrand is now analytic at $ t=0 $. Gaussian quadrature or even Simpson will work perfectly.

For a $ 1/x $ singularity, the double-exponential (tanh-sinh) transformation maps the interval to $ (-\infty, \infty) $ with exponential decay, making the integrand amenable to trapezoidal rule (which converges exponentially for analytic functions on $ \mathbb{R} $):

$ x = \tanh\left( \frac{\pi}{2} \sinh t \right), \quad dx = \frac{\pi}{2} \frac{\cosh t}{\cosh^2\left( \frac{\pi}{2} \sinh t \right)} dt $

This is the basis of tanh-sinh quadrature, arguably the most robust black-box integrator for 1D problems with endpoint singularities.

Tellegen's Integration (Complex Plane)

If you can evaluate $ f(z) $ for complex $ z $, integrating along a contour in the complex plane that avoids a real-axis singularity is a clean way to handle integrable singularities without subtraction or variable transformations.

Infinite Intervals

For $ \int_0^\infty f(x) dx $ with algebraic decay ($ f(x) \sim x^{-\alpha} $), map to $ [0,1] $:

$ t = \frac{x}{1+x}, \quad x = \frac{t}{1-t}, \quad dx = \frac{dt}{(1-t)^2} $

Then $ \int_0^\infty f(x) dx = \int_0^1 f\!\left(\frac{t}{1-t}\right) \frac{dt}{(1-t)^2} $. Gaussian quadrature on $ [0,1] $ handles the transformed integrand.

For exponential decay, Gauss-Laguerre is purpose-built. For Gaussian decay, Gauss-Hermite.

Practical Decision Tree

Here is the opinionated flowchart for picking a method:

1. 1D, smooth, finite interval: Gauss-Legendre, $ n=16 $ to $ n=64 $. Done.
2. 1D, unknown smoothness: Adaptive Gauss-Kronrod (G7/K15). scipy.integrate.quad.
3. 1D, known endpoint singularities: Gauss-Jacobi with matching $ \alpha, \beta $, or variable transformation + Gauss-Legendre.
4. 1D, periodic integrand: Composite trapezoidal on uniform grid — it converges exponentially for smooth periodic functions (Euler-Maclaurin cancellation).
5. 1D, infinite interval with known decay: Gauss-Laguerre or Gauss-Hermite.
6. 1D, infinite interval, algebraic decay: Map to $ [0,1] $ + Gauss-Legendre.
7. 2D-4D: Product Gauss-Legendre, or sparse grids (Smolyak) for high order.
8. >4D: Quasi-Monte Carlo (Sobol sequence).
9. Integral over a distribution defined by moments: QMOM (method of moments).

Romberg Integration: Richardson Extrapolation on Steroids

Romberg integration applies Richardson extrapolation to the composite trapezoidal rule. The trapezoidal rule with $ h $ has an error expansion in even powers of $ h $ (Euler-Maclaurin formula). By computing the trapezoidal approximation at successively halved step sizes $ h, h/2, h/4, \dots $, you can extrapolate away the $ h^2, h^4, h^6, \dots $ error terms.

def romberg(f, a, b, n_levels=6):
    """
    Romberg integration. Returns the extrapolated integral and the
    full Romberg tableau.
    """
    R = np.zeros((n_levels, n_levels))
    h = b - a
    R[0, 0] = 0.5 * h * (f(a) + f(b))

    for i in range(1, n_levels):
        h *= 0.5
        # Composite trapezoidal with 2^i panels
        n_panels = 2**i
        x_new = a + h * np.arange(1, n_panels, 2)  # Odd indices only (new points)
        R[i, 0] = 0.5 * R[i-1, 0] + h * np.sum(f(x_new))

        # Richardson extrapolation
        for k in range(1, i + 1):
            R[i, k] = R[i, k-1] + (R[i, k-1] - R[i-1, k-1]) / (4**k - 1)

    return R[n_levels-1, n_levels-1], R

Romberg is elegant but has been largely superseded by adaptive Gaussian quadrature. Its one remaining niche: when you can only sample $ f $ on a uniform grid and Richardson extrapolation is the only way to boost accuracy without additional evaluations.

Sparse Grids (Smolyak Quadrature)

For moderate dimensions ($ d \leq 20 $), full tensor-product Gaussian quadrature is intractable. Smolyak sparse grids combine 1D rules of varying orders into a $ d $-dimensional rule with far fewer points — $ \mathcal{O}(N (\log N)^{d-1}) $ instead of $ \mathcal{O}(N^d) $ — while preserving almost the same polynomial exactness.

The construction: for a 1D quadrature sequence $ Q_k $ (e.g., Gauss-Legendre with $ 2^k-1 $ points), the level-$ L $ Smolyak rule is:

$ Q_L^{(d)} = \sum_{\substack{\mathbf{k} \in \mathbb{N}^d \\ |\mathbf{k}|_1 \leq L + d - 1}} (-1)^{L+d-1-|\mathbf{k}|_1} \binom{d-1}{|\mathbf{k}|_1 - L} \bigotimes_{j=1}^{d} Q_{k_j} $

In practice, use sparsegrids or call into TASMANIAN (Oak Ridge's sparse grid library).

Sinc Quadrature

For analytic functions on $ \mathbb{R} $ with exponential decay, sinc quadrature (using $ \operatorname{sinc}(x) = \sin(\pi x)/(\pi x) $ as the cardinal function) converges exponentially:

$ \int_{-\infty}^{\infty} f(x) \, dx \approx h \sum_{k=-\infty}^{\infty} f(kh), \qquad h = \mathcal{O}(1/\sqrt{N}) $

The error is $ \mathcal{O}(e^{-c\sqrt{N}}) $. Sinc quadrature is closely related to the trapezoidal rule on $ \mathbb{R} $, which also converges exponentially for analytic, decaying functions. If you find yourself reaching for sinc quadrature, double-check whether tanh-sinh (double-exponential) quadrature would do the job instead — it handles endpoint singularities and is more robust for finite-interval problems.

Numerical Differentiation of Integrals (Leibniz Rule)

When you need the derivative of an integral with respect to a parameter:

$ \frac{d}{d\alpha} \int_{a(\alpha)}^{b(\alpha)} f(x, \alpha) \, dx = f(b(\alpha), \alpha) \, b'(\alpha) - f(a(\alpha), \alpha) \, a'(\alpha) + \int_{a(\alpha)}^{b(\alpha)} \frac{\partial f}{\partial \alpha} \, dx $

The first two terms come from the moving boundaries; the third from differentiating inside. In practice, use automatic differentiation (JAX, PyTorch) rather than finite differences on the integral — the gradient can flow through the quadrature nodes and weights directly.

# With JAX: automatic differentiation through quadrature is trivial
import jax.numpy as jnp
from jax import grad

def integral(param):
    """∫_0^1 sin(param * x) dx, computed via Gauss-Legendre."""
    x, w = leggauss(32)  # precomputed
    x_mapped = 0.5 * (x + 1)  # map to [0,1]
    return 0.5 * jnp.dot(w, jnp.sin(param * x_mapped))

dI_dparam = grad(integral)
print(dI_dparam(2.0))  # Exact derivative, no finite differences

Summary Table

Method	Degrees of Freedom	Degree/Convergence	Best For
Trapezoidal	$ N $ points	$ \mathcal{O}(h^2) $	Periodic smooth functions (exponential)
Simpson's	$ N $ panels	$ \mathcal{O}(h^4) $	Quick-and-dirty 1D
Gauss-Legendre	$ n $ nodes	degree $ 2n-1 $	1D smooth, finite interval — gold standard
Gauss-Kronrod	$ 2n+1 $ nodes	degree $ 3n+1 $	Adaptive 1D with error estimate
Clenshaw-Curtis	$ n+1 $ nodes	degree $ n $	Adaptive with nested nodes
QMOM	$ 2N $ moments	matches $ 2N $ moments	Distributions defined by moments
Monte Carlo	$ N $ samples	$ \mathcal{O}(1/\sqrt{N}) $	High dimension ($ d > 4 $)
Quasi-Monte Carlo	$ N $ samples	$ \mathcal{O}((\log N)^d/N) $	Moderate-high dimension
Romberg	$ 2^L+1 $ evals	Extrapolated trapezoidal	Uniform-grid data
Sparse Grids	$ \mathcal{O}(N(\log N)^{d-1}) $	Near-full tensor product	$ d \leq 20 $
Tanh-Sinh	$ N $ nodes	Exponential	Endpoint singularities

References

• Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P. Numerical Recipes in C, 2nd ed. Cambridge University Press, 1992. (Ch. 4: Integration of Functions)
• Golub, G.H. and Welsch, J.H. "Calculation of Gauss Quadrature Rules." Mathematics of Computation, 23(106):221-230, 1969.
• Gordon, R.G. "Error Bounds in Equilibrium Statistical Mechanics." J. Math. Phys., 9(5):655-663, 1968. (Product-Difference algorithm for QMOM)
• McGraw, R. "Description of Aerosol Dynamics by the Quadrature Method of Moments." Aerosol Sci. Tech., 27(2):255-265, 1997.
• Trefethen, L.N. "Is Gauss Quadrature Better than Clenshaw-Curtis?" SIAM Review, 50(1):67-87, 2008.
• Takahasi, H. and Mori, M. "Double Exponential Formulas for Numerical Integration." Publ. RIMS, Kyoto Univ., 9:721-741, 1974.
• Smolyak, S.A. "Quadrature and Interpolation Formulas for Tensor Products of Certain Classes of Functions." Dokl. Akad. Nauk SSSR, 148(5):1042-1045, 1963.
• Stroud, A.H. Approximate Calculation of Multiple Integrals. Prentice-Hall, 1971.