Admin: Create Numerics/Integration of Functions: comprehensive coverage of Gaussian quadrature, numerical integration, method of moments, adaptive quadrature, Monte Carlo, and more — Python-centric with code examples. (via create-page on MediaWiki MCP Server)

2026-06-22T04:46:30Z

Create Numerics/Integration of Functions: comprehensive coverage of Gaussian quadrature, numerical integration, method of moments, adaptive quadrature, Monte Carlo, and more — Python-centric with code examples. (via create-page on MediaWiki MCP Server)

New page

==Philosophy==

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

<math>
I = \int_{a}^{b} f(x) \, dx
</math>

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 <code>scipy.integrate.quad</code> 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:

<math>
I \approx \sum_{i=0}^{n} w_i \, f(x_i)
</math>

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

A quadrature rule of degree <math>d</math> integrates all polynomials up to order <math>d</math> ''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 <math>f(a)</math> and <math>f(b)</math>:

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

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

Degree: 1 (exact for linear polynomials).

===Simpson's Rule (n=2)===

Fit a parabola through three equally spaced points:

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

Error: <math>-\frac{(b-a)^5}{2880} f^{(4)}(\xi)</math>.

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 <math>x^3</math> error term cancels.

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

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

where <math>h = (b-a)/3</math>. 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:

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

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

Composite Simpson's: error <math>\mathcal{O}(h^4)</math>. 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 <math>x_i</math> and weights <math>w_i</math> to maximize the polynomial degree of exactness. With <math>n</math> points, a Gaussian rule has degree <math>2n - 1</math> — nearly double Newton-Cotes (which has degree <math>n-1</math> for even <math>n</math>).

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 <math>[-1, 1]</math> with unit weight function <math>w(x) = 1</math>:

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

where <math>x_i</math> are the roots of the Legendre polynomial <math>P_n(x)</math>, and

<math>
w_i = \frac{2}{(1 - x_i^2) [P_n'(x_i)]^2}
</math>

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

<syntaxhighlight lang="python">
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))
</syntaxhighlight>

===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 <math>n</math>-point Gauss rule (degree <math>2n-1</math>), and a <math>(2n+1)</math>-point Kronrod extension that reuses the <math>n</math> Gauss nodes and adds <math>n+1</math> 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 <code>scipy.integrate.quad</code>.

<syntaxhighlight lang="python">
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
</syntaxhighlight>

===Other Gaussian Flavors===

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

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

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

==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:

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

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

<syntaxhighlight lang="python">
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)
</syntaxhighlight>

==Clenshaw-Curtis (The FFT-Based Alternative)==

Clenshaw-Curtis quadrature uses Chebyshev nodes <math>x_k = \cos(k\pi/N)</math> (the extrema of <math>T_N(x)</math>) and weights computed via FFT. With <math>n+1</math> points, it has degree <math>n</math> — half that of Gauss-Legendre — but the weights are ''all positive'' and the nodes are nested (the <math>n</math>-point set is a subset of the <math>2n</math>-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 <math>\mathcal{O}(N \log N)</math> 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 <math>n</math> 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 <math>\mu_k = \int x^k f(x) \, dx</math> for <math>k = 0, 1, \dots, m</math>, approximate <math>\int g(x) f(x) \, dx</math> by constructing a quadrature rule that exactly integrates the first <math>m</math> monomials against <math>f(x)</math>.

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

===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:

<math>
n(x) \approx \sum_{i=1}^{N} w_i \, \delta(x - x_i)
</math>

The <math>2N</math> parameters <math>\{w_i, x_i\}</math> are determined by matching the first <math>2N</math> moments:

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

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).

<syntaxhighlight lang="python">
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])
</syntaxhighlight>

===Why QMOM Wins===

For tracking particle size distributions in CFD, QMOM with <math>N=3</math> or <math>N=4</math> 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 <math>n</math> points per axis in <math>d</math> dimensions uses <math>n^d</math> evaluations. For <math>d > 4</math>, you need Monte Carlo.

The Monte Carlo estimator:

<math>
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)
</math>

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

Variance reduction techniques:
* '''Importance sampling:''' Draw from <math>p(\mathbf{x}) \propto |f(\mathbf{x})|</math> instead of uniform. Reweight by <math>1/p(\mathbf{x})</math>.
* '''Stratified sampling:''' Partition <math>\Omega</math> 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 <math>\mathcal{O}((\log N)^d / N)</math>, which for moderate <math>d</math> and large <math>N</math> beats <math>1/\sqrt{N}</math> substantially.

<syntaxhighlight lang="python">
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
</syntaxhighlight>

==Dealing with Singularities==

===Subtract the Singularity===

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

<math>
\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
</math>

The first term is analytic; the second has a removable singularity (the numerator now vanishes at <math>x=0</math>) 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 <math>[0,1]</math> with <math>\sqrt{x}</math> singularities:

<math>
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
</math>

The integrand is now analytic at <math>t=0</math>. Gaussian quadrature or even Simpson will work perfectly.

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

<math>
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
</math>

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 <math>f(z)</math> for complex <math>z</math>, 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 <math>\int_0^\infty f(x) dx</math> with algebraic decay (<math>f(x) \sim x^{-\alpha}</math>), map to <math>[0,1]</math>:

<math>
t = \frac{x}{1+x}, \quad x = \frac{t}{1-t}, \quad dx = \frac{dt}{(1-t)^2}
</math>

Then <math>\int_0^\infty f(x) dx = \int_0^1 f\!\left(\frac{t}{1-t}\right) \frac{dt}{(1-t)^2}</math>. Gaussian quadrature on <math>[0,1]</math> 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:

# '''1D, smooth, finite interval:''' Gauss-Legendre, <math>n=16</math> to <math>n=64</math>. Done.
# '''1D, unknown smoothness:''' Adaptive Gauss-Kronrod (G7/K15). <code>scipy.integrate.quad</code>.
# '''1D, known endpoint singularities:''' Gauss-Jacobi with matching <math>\alpha, \beta</math>, ''or'' variable transformation + Gauss-Legendre.
# '''1D, periodic integrand:''' Composite trapezoidal on uniform grid — it converges ''exponentially'' for smooth periodic functions (Euler-Maclaurin cancellation).
# '''1D, infinite interval with known decay:''' Gauss-Laguerre or Gauss-Hermite.
# '''1D, infinite interval, algebraic decay:''' Map to <math>[0,1]</math> + Gauss-Legendre.
# '''2D-4D:''' Product Gauss-Legendre, or sparse grids (Smolyak) for high order.
# '''>4D:''' Quasi-Monte Carlo (Sobol sequence).
# '''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 <math>h</math> has an error expansion in even powers of <math>h</math> (Euler-Maclaurin formula). By computing the trapezoidal approximation at successively halved step sizes <math>h, h/2, h/4, \dots</math>, you can extrapolate away the <math>h^2, h^4, h^6, \dots</math> error terms.

<syntaxhighlight lang="python">
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
</syntaxhighlight>

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

==Sparse Grids (Smolyak Quadrature)==

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

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

<math>
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}
</math>

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

==Sinc Quadrature==

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

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

The error is <math>\mathcal{O}(e^{-c\sqrt{N}})</math>. Sinc quadrature is closely related to the trapezoidal rule on <math>\mathbb{R}</math>, 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:

<math>
\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
</math>

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.

<syntaxhighlight lang="python">
# 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
</syntaxhighlight>

==Summary Table==

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

[[Category:Numerics]]
[[Category:Numerical Recipes]]
[[Category:Python]]

Numerics/Integration of Functions - Revision history

Admin: Create Numerics/Integration of Functions: comprehensive coverage of Gaussian quadrature, numerical integration, method of moments, adaptive quadrature, Monte Carlo, and more — Python-centric with code examples. (via create-page on MediaWiki MCP Server)