Gaussian Quadrature: Smart Numerical Integration Explained

Gaussian quadrature approximates a definite integral with very few function evaluations. The trapezoidal and Simpson rules fix their sample points and only optimize the weights — Gauss asked: what if we choose points and weights both optimally? With just n points, the rule then integrates every polynomial up to degree 2n − 1 exactly, double what n fixed points could manage.

Function evaluations are often expensive — each might run a simulation or solve an equation — so Simpson-level accuracy from a handful of points saves real computing time. Gaussian quadrature is the default integrator in finite-element software, physics codes, and statistics libraries for exactly this reason.

Key Concepts

• Gauss-Legendre rule: the standard form on the interval [−1, 1]. The nodes are the roots of the n-th Legendre polynomial, and the weights are chosen to make the rule exact for polynomials up to degree 2n − 1.
• Nodes are not evenly spaced: they cluster toward the ends of the interval, which is precisely what buys the extra accuracy.
• Change of interval: any [a, b] is mapped to [−1, 1] with a linear substitution before applying the rule.
• Error order: the error involves the (2n)-th derivative of f, so the rule is spectacular for smooth functions and ordinary for rough ones.

Formulas

• Gauss-Legendre rule: ∫_−1^1 f(x) dx ≈ Σ_i=1..n w_i f(x_i).
• Two-point rule: nodes x = ±1/√3 ≈ ±0.5773503, weights w₁ = w₂ = 1. Exact for all polynomials up to degree 3.
• Three-point rule: nodes 0 and ±√(3/5) ≈ ±0.7745967, weights 8/9 and 5/9. Exact up to degree 5.
• Interval transformation: ∫_a^b f(x) dx = (b − a)/2 · ∫_−1^1 f(((b − a)/2)t + (a + b)/2) dt. The error is proportional to f^(2n)(ξ), assuming f is 2n times differentiable.

Worked Example

Compute ∫₀¹ x³ dx with the two-point Gauss rule. The exact answer is 1/4 = 0.25, and since x³ has degree 3 = 2·2 − 1, the rule should be exact. Map [0, 1] to [−1, 1]: x = (t + 1)/2, dx = dt/2. Nodes t = ±0.5773503 become x₁ = 0.2113249, x₂ = 0.7886751. Estimate: 0.5 × (0.2113249³ + 0.7886751³) = 0.5 × (0.0094374 + 0.4905626) = 0.5 × 0.5 = 0.25 — exactly the true integral, from only two function values where Simpson's rule would need three.

Common Misconceptions

• "The sample points are evenly spaced." They are deliberately not — the nodes are roots of Legendre polynomials, denser near the endpoints. Even spacing would destroy the 2n − 1 exactness.
• "More points are always worth it." For smooth functions yes, but for functions with kinks or singularities the high derivatives in the error term are large, and splitting the interval (composite Gauss rules) beats raising n.
• "It works only on [−1, 1]." The linear change of variables maps any finite interval onto [−1, 1]; variants (Gauss-Laguerre, Gauss-Hermite) handle infinite intervals.

Real-World Applications

• Finite element analysis: every element stiffness matrix is assembled by Gaussian quadrature.
• Quantum physics: expectation values of wavefunctions need high-accuracy integrals from few points.
• Statistics: Bayesian computations integrate posterior densities with Gauss-Hermite rules.
• Finance: pricing options under complex models integrates payoff functions numerically.

Gaussian quadrature chooses both nodes and weights optimally, so n points integrate every polynomial up to degree 2n − 1 exactly — twice the power of fixed-point rules. With the two-point rule (±1/√3, weights 1) already exact through cubics, it delivers maximum accuracy per function evaluation for smooth problems.