Mathematical foundations of Gauss-Legendre Quadrature.
While Newton-Cotes rules like the trapezoidal and Simpson formulas restrict evaluation nodes to equally spaced grids, Carl Friedrich Gauss posed a profound question in 1814: what if we treat the evaluation nodes themselves as free parameters? By choosing the n nodes to be the roots of the degree-n Legendre orthogonal polynomial, Gauss proved that one can integrate every polynomial of degree up to 2n - 1 exactly—doubling the algebraic precision of standard formulas using the exact same number of function evaluations. Furthermore, all Gauss-Legendre weights are unconditionally positive, eliminating the subtractive cancellation that plagues high-order Newton-Cotes rules. In this computational laboratory, we compute the roots of Legendre polynomials, verify degree-5 exactness for n = 3 points, and demonstrate the sixth-order O(h⁶) super-convergence of composite Gauss-Legendre quadrature.
Viva voce questions & analytical theory
Why is Gauss-Legendre quadrature called 'optimal'?
In standard Newton-Cotes rules, the n node locations are fixed a priori, leaving n weights to satisfy n equations (degree n - 1 exactness). Carl Friedrich Gauss treated both the n nodes and n weights as 2n free parameters, choosing nodes as roots of orthogonal Legendre polynomials to achieve exactness for all polynomials of degree ≤ 2n - 1.
How are Legendre polynomials defined and why are their roots ideal for quadrature?
Legendre polynomials P_n(x) are orthogonal on [-1, 1] with respect to weight w(x) = 1: ∫_{-1}^1 P_n(x) P_m(x) dx = (2 / (2n+1)) δ_{nm}. Because P_n(x) is orthogonal to all polynomials of degree < n, any polynomial f of degree ≤ 2n - 1 can be divided as f = q · P_n + r (with deg(q), deg(r) < n), making ∫ f = ∫ r exact at the roots of P_n.
What is the Golub-Welsch algorithm?
The Golub-Welsch algorithm (1969) computes Gauss quadrature nodes and weights by framing the orthogonal polynomial recurrence as a symmetric tridiagonal Jacobi matrix. The eigenvalues of this matrix are the quadrature nodes, and the squared first components of normalized eigenvectors yield the weights.
Why are all Gauss-Legendre quadrature weights strictly positive?
Because the cardinal basis polynomials L_{n,j}(x) satisfy ∫_{-1}^1 (L_{n,j}(x))² dx > 0. Since (L_{n,j}(x))² has degree 2n - 2 ≤ 2n - 1, the quadrature formula integrates it exactly: w_j = ∑ w_k (L_{n,j}(x_k))² = ∫ (L_{n,j}(x))² dx > 0, ensuring unconditional numerical stability.
What is the asymptotic error convergence rate for composite n-point Gauss quadrature?
For an n-point rule applied over sub-panels of width h, the error scales as O(h^{2n}). For n = 3, the composite rule achieves sixth-order convergence O(h⁶), halving step size reduces error by 2⁶ = 64.
How does Gauss-Legendre quadrature differ from Gauss-Chebyshev or Gauss-Laguerre quadrature?
Gauss quadrature generalizes to any weight function w(x) on interval [a, b]. Gauss-Legendre uses w(x)=1 on [-1, 1]; Gauss-Chebyshev uses w(x)=1/√(1-x²) on [-1, 1]; Gauss-Laguerre uses w(x)=e^{-x} on [0, ∞); and Gauss-Hermite uses w(x)=e^{-x²} on (-∞, ∞).
Why are Gauss-Kronrod extensions used in adaptive integration packages (like QUADPACK / scipy.integrate.quad)?
Gauss-Legendre nodes cannot be reused when increasing order from n to 2n+1. Alexander Kronrod developed optimal extensions that add n + 1 new nodes to existing n Gauss nodes, allowing reuse of all n function evaluations to compute an embedded error estimate.
Algorithm precautions and stability rules
- Always perform the affine transformation x = (b-a)/2 · ξ + (a+b)/2 with full 64-bit precision when mapping canonical nodes from [-1, 1] to [a, b].
- Ensure precomputed Legendre roots and weights are stored with at least 16 significant decimal digits to prevent precision loss.
- Verify that the integrand has at least 2n continuous derivatives to achieve the theoretical O(h^{2n}) convergence rate.
- When integrating over unbounded intervals, do not use Gauss-Legendre; switch to Gauss-Laguerre or Gauss-Hermite quadrature.
- Sum positive weighted contributions using sorted or compensated summation to prevent floating point accumulation bias.
- For functions with interior discontinuities, subdivide the domain at the singular points so each sub-panel remains smooth.
Theoretical limitations
- Standard Gauss-Legendre nodes are non-nested: doubling n requires discarding all previous function evaluations unless Gauss-Kronrod extensions are used.
- Gauss quadrature does not automatically evaluate endpoints, making it inconvenient if boundary values are needed for other purposes.
- Like all 1D quadrature rules, tensor-product extension to d dimensions requires n^d points, suffering from exponential curse of dimensionality.