Gauss–Legendre Quadrature

Gauss–Legendre Quadrature is a high-order numerical integration algorithm in computation and algorithms and numerical analysis that approximates the definite integral of a continuous function f(x) over the canonical interval [−1, 1] (or any finite interval [a, b] via affine coordinate transformation) as a weighted sum of n function evaluations at non-equispaced interior nodes x₁, x₂, ..., x_n. Introduced by Carl Friedrich Gauss in an 1814 memoir using continued-fraction expansions of hypergeometric series and reformulated in terms of orthogonal Legendre polynomials P_n(x) by Carl Gustav Jacob Jacobi in 1826, it is the archetype of the broader family of Gaussian quadrature rules.

Whereas Newton–Cotes formulas such as the trapezoidal rule, Simpson's rule, and Romberg's method fix the n evaluation nodes at predetermined equally spaced grid points—leaving only the n weights w_i free to match polynomials of degree at most n − 1 (or n for odd n)—Gauss–Legendre quadrature optimizes both the n node positions x_i and the n weights w_i simultaneously. Exploiting all 2n degrees of freedom achieves the maximal possible algebraic degree of precision 2n − 1: an n-point Gauss–Legendre rule integrates every polynomial of degree ≤ 2n − 1 with zero truncation error, avoids the Runge phenomenon of high-degree equispaced interpolation, and converges geometrically for analytic integrands.

Maximal Degree of Precision via Legendre Polynomial Roots

A general n-point quadrature rule with unit weight function ω(x) = 1 on [−1, 1] takes the form:

Let P_n(x) denote the n-th Legendre polynomial, defined by Rodrigues' formula P_n(x) = (2ⁿ n!)⁻¹ (dⁿ/dxⁿ)[(x² − 1)ⁿ] or the three-term recurrence (k + 1)P_k+1(x) = (2k + 1)*xP_k*(x) − kP_k−1(x) with P₀(x) = 1 and P₁(x) = x. The Legendre polynomials form an orthogonal family in L²([−1, 1]):

Consequently, P_n(x) is orthogonal to every polynomial q(x) of degree ≤ n − 1: ∫₋₁¹ P_n(x)q(x) dx = 0. Furthermore, all n roots of P_n(x) are real, simple, symmetric about the origin, and lie inside the open interval (−1, 1).

To see why placing the nodes x₁, ..., x_n at the roots of P_n(x) (P_n(x_i) = 0) yields precision 2n − 1, let p(x) be any polynomial of degree ≤ 2n − 1. Dividing p(x) by P_n(x) using polynomial Euclidean division gives:

Integrating p(x) over [−1, 1], the first term vanishes by orthogonality of P_n(x) to q(x). Evaluating p(x) at the roots x_i of P_n(x) gives p(x_i) = q(x_i) · 0 + r(x_i) = r(x_i). Choosing the n interpolatory Lagrange weights w_i = ∫₋₁¹ ℓ_i(x) dx (which integrate degree-(n − 1) polynomials r(x) with zero error) yields:

Conversely, no n-point rule can attain degree of precision 2n, because the degree-2n non-negative polynomial π_n²(x) = ∏_i=1..n (x − x_i)² has a positive integral ∫₋₁¹ π_n²(x) dx > 0 while evaluating to zero at all n nodes.

Weight Formula, Positivity, and Affine Interval Transformation

Using the Christoffel–Darboux identity for Legendre polynomials, each quadrature weight w_i has the closed-form expression:

Because every weight w_i is positive (in contrast to high-order Newton–Cotes rules, which develop large alternating negative weights for n ≥ 9), Gauss–Legendre quadrature is unconditionally stable against roundoff amplification in IEEE 754 arithmetic and converges for every Riemann-integrable function as n → ∞ (Steklov's theorem). To integrate over an arbitrary finite interval [a, b], the affine change of variables t = ((b − a) / 2)x + (a + b) / 2 maps [−1, 1] onto [a, b]:

Peano Kernel Error Formula and Geometric Convergence

For a function f ∈ C²ⁿ([a, b]), Hermite osculatory interpolation yields the classical truncation error remainder for some ξ ∈ (a, b):

By Stirling's approximation, the factorial coefficient decays as (π / 4ⁿ)², and for any function f(z) analytic inside a Bernstein ellipse E_ρ in the complex plane with foci at ±1 and semiaxes summing to ρ > 1, the quadrature error decays geometrically as O(ρ⁻²ⁿ).

Computing Nodes and Weights: Golub–Welsch Versus Glaser–Liu–Rokhlin

To compute (x_i, w_i) numerically for arbitrary n, Gene Golub and John Welsch (1969) formulated the Golub–Welsch algorithm, which constructs the symmetric tridiagonal Jacobi matrix J_n ∈ ℝ^n×n from the three-term recurrence coefficients (zero diagonal, super/subdiagonal entries β_k = k / √(4k² − 1) for k = 1, ..., n − 1). The eigenvalues of J_n (computed via the symmetric tridiagonal QR algorithm) are the roots x_i, and the weights equal twice the square of the first component of each normalized eigenvector (w_i = 2 q_1,i²) in O(n²) operations. For large n ≥ 100, asymptotic Newton's method root-finders initialized by Tricomi–Bessel expansions (Glaser, Liu, & Rokhlin, 2007; Hale & Townsend, 2013) compute all n nodes and weights in O(n) linear time to full machine precision.

Because the roots of P_n(x) and P_2n(x) do not overlap (except x = 0 for odd n), doubling n in a naive adaptive integrator discards all previous function evaluations. Gauss–Kronrod quadrature (such as the G7–K15 and G10–K21 pairs in QUADPACK) remedies this by interleaving n + 1 Kronrod nodes (roots of the Stieltjes polynomial E_n+1(x)) between the n Gauss nodes, reusing all n existing evaluations to produce an embedded (2n + 1)-point rule of degree 3n + 1 (or 3n + 2) for a-posteriori error control. In multidimensional scientific computing, tensor-product Gauss–Legendre points serve as the standard stiffness-matrix integration quadrature across isoparametric finite elements and spectral-element fluid dynamics solvers.

Why are Gauss–Legendre nodes clustered near the endpoints ±1 rather than equally spaced?

The roots of P_n(x) satisfy the asymptotic angular relation x_i ≈ −cos((4i − 1)π / (4n + 2)), projecting uniform angles onto the horizontal diameter of the unit circle. This Chebyshev-like O(1/n²) clustering near ±1 counteracts the growth of the polynomial nodal product ∏(x − x_i) near the boundaries, preventing the wild boundary oscillations (Runge's phenomenon) that destabilize high-order equispaced Newton–Cotes formulas.

When can Romberg's method or the trapezoidal rule outperform Gauss–Legendre quadrature?

Gauss–Legendre quadrature is optimal for general smooth non-periodic functions on finite intervals, achieving degree 2n − 1. However, for smooth periodic functions integrated over a full period (or functions decaying doubly exponentially after a tanh-sinh transformation), the simple equispaced trapezoidal rule and Romberg's method enjoy Euler–Maclaurin boundary cancellation of all odd derivatives and allow complete 100% node reuse when halving the step size.

How does Gauss–Legendre quadrature handle integrands with endpoint singularities?

Because all Gauss–Legendre nodes x_i lie strictly inside the open interval (−1, 1), the rule never evaluates f(x) at x = ±1 and can be applied directly to integrable endpoint singularities such as f(x) = (1 − x)^−1/2 without dividing by zero. However, because f⁽²ⁿ⁾(x) diverges at the endpoint, convergence slows from geometric to algebraic unless the algebraic singularity is absorbed into a weighted Gauss–Jacobi quadrature rule with weight ω(x) = (1 − x)^α(1 + x)^β.

1. Gauss, C. F. (1814). "Methodus nova integralium valores per approximationem inveniendi." Commentationes Societatis Regiae Scientiarum Gottingensis Recentiores, 3, 39–76. https://gdz.sub.uni-goettingen.de/id/PPN35283028X_0003
2. Golub, G. H., & Welsch, J. H. (1969). "Calculation of Gauss quadrature rules." Mathematics of Computation, 23(106), 221–230. https://doi.org/10.1090/S0025-5718-69-99647-1
3. Laurie, D. P. (1997). "Calculation of Gauss-Kronrod quadrature rules." Mathematics of Computation, 66(219), 1133–1145. https://doi.org/10.1090/S0025-5718-97-00861-2
4. Trefethen, L. N. (2008). "Is Gauss quadrature better than Clenshaw–Curtis?" SIAM Review, 50(1), 67–87. https://doi.org/10.1137/060659831
5. Hale, N., & Townsend, A. (2013). "Fast and accurate computation of Gauss–Legendre and Gauss–Jacobi quadrature nodes and weights." SIAM Journal on Scientific Computing, 35(2), A652–A674. https://doi.org/10.1137/120889873