Gaussian quadrature

Deep dive

How it works

In: integrand f, interval [a, b], point count n or tolerance · Out: integral estimate (plus error estimate if adaptive) · State: node/weight table per n, cacheable · Cost: n evaluations of f; setup O(n²) by Newton, O(n) by asymptotics

1. Nodes are the roots of Pₙ, from the recurrence k·P_k = (2k-1)·x·P_{k-1} − (k-1)·P_{k-2}.
2. Newton (Newton root finding) from the guess cos(π(i − ¼)/(n + ½)), with Pₙ′ = n(xPₙ − P_{n-1})/(x² − 1).
3. Weights wᵢ = 2/((1 − xᵢ²)·Pₙ′(xᵢ)²); they sum to 2.
4. Map to [a, b]: x = (b−a)/2·t + (a+b)/2, so ∫f ≈ (b−a)/2·Σ wᵢ f(x(tᵢ)).
5. Adaptive: compare the n- and 2n-point results, accept below the tolerance, otherwise bisect with half the tolerance each.

Why 2n-1: write a polynomial of degree ≤ 2n-1 as Pₙ·q + r with deg q, r < n. ∫Pₙ·q = 0 by orthogonality, the rule gives 0 because Pₙ vanishes at the nodes, and r is integrated exactly by interpolation.

Nodes by eigenvalues. Golub-Welsch (1969): nodes are the eigenvalues of the symmetric tridiagonal Jacobi matrix of the recurrence; weights come from the first eigenvector components. Hale-Townsend and Bogaert give O(n) iteration-free methods for large n.

Other weights. Hermite (e^(-x²) on ℝ), Laguerre (e^(-x) on [0, ∞)), Jacobi ((1-x)^α(1+x)^β). Radau fixes one endpoint as a node (exact to 2n-2), Lobatto both (2n-3).

Gauss-Kronrod. Add n+1 nodes to the n Gauss nodes so the 2n+1 points form a higher-order rule: G7/K15 reuses the 7 Gauss evaluations, and |K − G| estimates the error. QUADPACK's QAGS uses 21-point Gauss-Kronrod with global bisection and extrapolation for endpoint singularities.

Code

use std::f64::consts::PI;

/// P_n(z) and P_n'(z) from the three-term recurrence k P_k = (2k-1) z P_{k-1} - (k-1) P_{k-2}.
fn legendre(n: usize, z: f64) -> (f64, f64) {
    let (mut p0, mut p1) = (1.0, z);
    for k in 2..=n {
        let kf = k as f64;
        (p0, p1) = (p1, ((2.0 * kf - 1.0) * z * p1 - (kf - 1.0) * p0) / kf);
    }
    (p1, n as f64 * (z * p1 - p0) / (z * z - 1.0))
}

/// Gauss-Legendre rule on [-1, 1]: (nodes ascending, weights).
pub fn gauss_legendre(n: usize) -> (Vec<f64>, Vec<f64>) {
    let (mut x, mut w) = (vec![0.0; n], vec![0.0; n]);
    for i in 0..(n + 1) / 2 {
        let mut z = (PI * (i as f64 + 0.75) / (n as f64 + 0.5)).cos(); // guess for the i-th largest root
        for _ in 0..100 {
            let (p, dp) = legendre(n, z);
            let dz = p / dp;
            z -= dz; // Newton step
            if dz.abs() < 1e-15 { break; }
        }
        let dp = legendre(n, z).1;
        (x[i], x[n - 1 - i]) = (-z, z);
        w[i] = 2.0 / ((1.0 - z * z) * dp * dp);
        w[n - 1 - i] = w[i];
    }
    (x, w)
}

/// Fixed n-point rule on [a, b]: x = (b-a)/2 t + (a+b)/2, dx = (b-a)/2 dt.
pub fn integrate(f: &dyn Fn(f64) -> f64, a: f64, b: f64, rule: &(Vec<f64>, Vec<f64>)) -> f64 {
    let (h, m) = (0.5 * (b - a), 0.5 * (a + b));
    h * rule.0.iter().zip(&rule.1).map(|(t, w)| w * f(m + h * t)).sum::<f64>()
}

/// Adaptive bisection: error estimate = |G_2n - G_n|; tolerance halves per split.
pub fn adaptive(f: &dyn Fn(f64) -> f64, a: f64, b: f64, tol: f64, n: usize, evals: &mut usize) -> f64 {
    let (lo, hi) = (gauss_legendre(n), gauss_legendre(2 * n));
    fn rec(f: &dyn Fn(f64) -> f64, a: f64, b: f64, tol: f64, lo: &(Vec<f64>, Vec<f64>),
           hi: &(Vec<f64>, Vec<f64>), depth: u32, evals: &mut usize) -> f64 {
        let (g, g2) = (integrate(f, a, b, lo), integrate(f, a, b, hi));
        *evals += lo.0.len() + hi.0.len();
        if (g2 - g).abs() <= tol || depth == 50 { return g2; }
        let m = 0.5 * (a + b);
        rec(f, a, m, 0.5 * tol, lo, hi, depth + 1, evals) + rec(f, m, b, 0.5 * tol, lo, hi, depth + 1, evals)
    }
    rec(f, a, b, tol, &lo, &hi, 0, evals)
}

fn main() {
    // weights sum to 2, nodes symmetric
    for n in [1, 2, 3, 5, 8, 20, 64] {
        let (x, w) = gauss_legendre(n);
        assert!((w.iter().sum::<f64>() - 2.0).abs() < 1e-14, "sum n={n}");
        for i in 0..n { assert!((x[i] + x[n - 1 - i]).abs() < 1e-15); }
    }
    // exact for x^k, k <= 2n-1; not for x^(2n)
    for n in [2, 4, 8, 16] {
        let rule = gauss_legendre(n);
        let mut worst: f64 = 0.0;
        for k in 0..2 * n {
            let ex = if k % 2 == 0 { 2.0 / (k as f64 + 1.0) } else { 0.0 };
            worst = worst.max((integrate(&|x| x.powi(k as i32), -1.0, 1.0, &rule) - ex).abs());
        }
        assert!(worst < 1e-14, "n={n} worst={worst}");
        let k = 2 * n;
        let err = (integrate(&|x| x.powi(k as i32), -1.0, 1.0, &rule) - 2.0 / (k as f64 + 1.0)).abs();
        assert!(err > 1e-12);
        println!("n={n:2}: worst error k<=2n-1 = {worst:.1e}; error at x^{k} = {err:.2e}");
    }
    // sin on [0, pi]
    for n in [2, 5, 7, 8] {
        let e = (integrate(&|x: f64| x.sin(), 0.0, PI, &gauss_legendre(n)) - 2.0).abs();
        println!("sin n={n}: error {e:.2e}");
        if n == 8 { assert!(e < 1e-12); }
    }
    // sqrt on [0, 1]
    for tol in [1e-6, 1e-10] {
        let mut ev = 0;
        let v = adaptive(&|x: f64| x.sqrt(), 0.0, 1.0, tol, 5, &mut ev);
        let e = (v - 2.0 / 3.0).abs();
        assert!(e <= tol);
        let r = integrate(&|x: f64| x.sqrt(), 0.0, 1.0, &gauss_legendre(15));
        println!("sqrt tol={tol:.0e}: error {e:.2e}, evals {ev}; fixed 15-pt error {:.2e}", (r - 2.0 / 3.0).abs());
    }
}

Complexity and constants

• Newton: O(n) per recurrence, a few iterations per root, n/2 roots: O(n²) per rule.
• Measured above (monomials on [-1, 1]), n = 2, 4, 8, 16: worst error over k ≤ 2n-1 is at most 4.4·10⁻¹⁶; x^(2n) misses by 1.8·10⁻¹, 1.2·10⁻², 4.7·10⁻⁵, 7.2·10⁻¹⁰.
• ∫₀^π sin = 2: error 6.4·10⁻² (n = 2), 1.1·10⁻⁷ (5), 1.8·10⁻¹² (7), 4.2·10⁻¹⁵ (8).
• ∫₀¹ √x = 2/3: 15 fixed points leave 2.8·10⁻⁵. Adaptive 5/10: tolerance 10⁻⁶ used 585 evaluations, ended 2.4·10⁻¹³ off; tolerance 10⁻¹⁰ used 1,455, ended 1.1·10⁻¹⁶ off. The estimate bounds the lower rule's error, so the returned G₂ₙ beats the tolerance.

Efficiency

• Cost: n evaluations of f plus a dot product; a cached table is free. Bound: evaluations of f.
• Levers:
  • cache nodes and weights per n — setup O(n²) → 0 — repeated calls, element loops
  • Gauss-Kronrod pair — reuses the n Gauss evaluations in the error estimate (QUADPACK design) — adaptive codes
  • split at known kinks and singularities — restores geometric rate — piecewise-smooth f
  • O(n) node algorithms — FastGaussQuadrature.jl README shows 100,000 nodes in about 2 ms — n in the thousands
• Naive → standard → state of the art:
  • fixed Gauss-Legendre, nodes recomputed per call · O(n²) setup + n evaluations
  • cached table, adaptive n/2n bisection · 585 evaluations for √x at 10⁻⁶ (measured)
  • adaptive Gauss-Kronrod with extrapolation (QUADPACK QAGS, GSL, Boost) · singularities handled, error bound returned

Implementation notes

• The starting guess is what sends Newton to the right root; keep it. Stop on |dz| < 1e-15.
• Gauss-Hermite for E[g(Z)], Z ~ N(μ, σ²): x = μ + √2·σ·t, weights wᵢ/√π.
• Positive weights keep round-off near ε·Σ|f| (floating-point error).

Pitfalls

• Kinks and endpoint singularities (√x, log x) make convergence algebraic: split or substitute.
• Oscillatory f (sin(100x)): needs several points per period; split into panels or use Filon/Levin-type methods.
• n versus 2n Gauss reuses no points; Kronrod does.
• An n versus 2n difference can be small by luck on non-analytic f.

Variants and alternatives

In production

Engineering: element Gauss rules in FEM codes. Libraries: QUADPACK inside SciPy quad and GSL, Boost.Math gauss_kronrod, FastGaussQuadrature.jl. Statistics: lme4 glmer uses adaptive Gauss-Hermite with nAGQ points per axis. Finance: Gauss-Hermite expectations in pricing.

Implementations

Sources

• Gene H. Golub, John H. Welsch, Calculation of Gauss Quadrature Rules, Mathematics of Computation 23, 1969. paper
• Nicholas Hale, Alex Townsend, Fast and Accurate Computation of Gauss-Legendre and Gauss-Jacobi Quadrature Nodes and Weights, SIAM J. Sci. Comput. 35, 2013. paper
• Ignace Bogaert, Iteration-Free Computation of Gauss-Legendre Quadrature Nodes and Weights, SIAM J. Sci. Comput. 36, 2014. paper
• Lloyd N. Trefethen, Is Gauss Quadrature Better than Clenshaw-Curtis?, SIAM Review 50, 2008. paper
• R. Piessens, E. de Doncker-Kapenga, C. W. Überhuber, D. K. Kahaner, QUADPACK: A Subroutine Package for Automatic Integration, Springer, 1983. book · code
• Boost.Math, Gauss-Kronrod Quadrature. docs · GSL, Numerical Integration. docs