13.1 Gaussian Quadrature Theory

Gaussian Quadrature for Integration

Fundamental Principles and Advantages

Gaussian quadrature approximates definite integrals by choosing evaluation points (nodes) and weights that maximize accuracy for a given number of function evaluations. Unlike Newton-Cotes methods that use equally spaced points, Gaussian quadrature places nodes strategically to extract more information from each evaluation.

The central result is this: with $n$ nodes, Gaussian quadrature exactly integrates polynomials of degree $2n - 1$ or less. That's roughly twice the polynomial degree you'd expect from $n$ points, and it's the reason the method is so efficient. You get both the node locations $x_i$ and the weights $w_i$ as free parameters, so you have $2n$ unknowns to satisfy $2n$ polynomial exactness conditions.

For smooth, well-behaved functions, this translates to high accuracy with relatively few function evaluations. The tradeoff is that the nodes are no longer equally spaced, so you can't reuse function values if you increase $n$.

Efficiency and Applications

Gaussian quadrature is commonly used across scientific computing:

• Computational physics: quantum mechanics calculations, scattering problems
• Engineering: structural analysis via finite element methods, heat transfer
• Finance: option pricing models requiring numerical integration
• Statistics: evaluating probability distributions that lack closed-form CDFs

The method adapts to different problems by choosing an appropriate weight function. Each weight function corresponds to a different family of Gaussian quadrature (covered in detail below).

Gaussian Quadrature Formula and Orthogonal Polynomials

General Formula and Components

The general Gaussian quadrature formula is:

$\int_{a}^{b} w(x)f(x)\,dx \approx \sum_{i=1}^{n} w_i f(x_i)$

where:

• $w(x)$ is the weight function (nonnegative, integrable over $[a,b]$)
• $f(x)$ is the integrand you want to approximate
• $x_i$ are the nodes (abscissas), determined as roots of the degree-$n$ orthogonal polynomial associated with $w(x)$
• $w_i$ are the weights, computed so the formula is exact for polynomials up to degree $2n - 1$

The weights can be expressed through Lagrange basis polynomials. Specifically, $w_i = \int_a^b w(x) L_i(x)\,dx$, where $L_i(x)$ is the $i$-th Lagrange interpolating polynomial built on the nodes. This construction guarantees the maximum degree of precision.

Role of Orthogonal Polynomials

Orthogonal polynomials are what make Gaussian quadrature work. A family of polynomials $\{P_0, P_1, P_2, \ldots\}$ is orthogonal with respect to $w(x)$ on $[a,b]$ if:

$\int_a^b w(x) P_m(x) P_k(x)\,dx = 0 \quad \text{for } m \neq k$

The nodes of $n$-point Gaussian quadrature are the roots of $P_n(x)$. This choice is what pushes the exactness from degree $n - 1$ (which any interpolatory rule achieves) up to degree $2n - 1$.

Different families correspond to different quadrature types:

$P_{k+1}(x) = (A_k x + B_k)P_k(x) - C_k P_{k-1}(x)$

This recurrence is essential in practice because it lets you compute nodes and weights efficiently and with good numerical stability, rather than finding polynomial roots from explicit coefficient formulas.

Gaussian Quadrature Accuracy vs Other Methods

Convergence Properties

For analytic (infinitely smooth) functions, Gaussian quadrature exhibits exponential convergence: the error decreases exponentially as you increase $n$. This is dramatically faster than the algebraic convergence of Newton-Cotes rules.

The error for $n$-point Gaussian quadrature applied to a function $f$ has the form:

$E_n(f) = \frac{f^{(2n)}(\xi)}{(2n)!} \int_a^b w(x)\left[P_n(x)\right]^2 dx$

for some $\xi \in (a,b)$. Notice the $2n$-th derivative appears, confirming exactness for polynomials of degree up to $2n - 1$.

Comparative Analysis

Here's how the error scaling compares for a fixed step size $h$ or number of nodes:

• Trapezoidal rule: error $O(h^2)$, exact for polynomials up to degree 1
• Simpson's rule: error $O(h^4)$, exact for polynomials up to degree 3
• $n$-point Gaussian quadrature: exact for polynomials up to degree $2n - 1$

For smooth functions, Gaussian quadrature needs far fewer evaluations to reach a target accuracy. However, the method can struggle with:

• Singularities within or at the endpoints of the integration interval
• Rapid oscillations that aren't well-captured by a polynomial approximation
• Discontinuities in the integrand or its derivatives

When these issues arise, two practical strategies help:

1. Composite Gaussian quadrature: subdivide the interval into smaller subintervals and apply Gaussian quadrature on each one
2. Adaptive Gaussian quadrature: automatically refine the subintervals where the integrand is difficult, concentrating effort where it's needed

The Lebesgue constant, which measures how the interpolation error amplifies, grows slowly for Gaussian nodes. This contributes to the method's numerical stability compared to equally spaced interpolation, where the Lebesgue constant grows exponentially (the Runge phenomenon).

Applying Gaussian Quadrature to Integrals

Specific Gaussian Quadrature Types

Gauss-Legendre quadrature is the most commonly used variant. It has weight function $w(x) = 1$ on $[-1, 1]$. To apply it to an arbitrary interval $[a, b]$, use the linear change of variables:

$x = \frac{b - a}{2}t + \frac{b + a}{2}, \quad t \in [-1, 1]$

which transforms the integral as:

$\int_a^b f(x)\,dx = \frac{b-a}{2}\int_{-1}^{1} f\!\left(\frac{b-a}{2}t + \frac{b+a}{2}\right)dt$

You then apply the standard Gauss-Legendre nodes and weights to the right-hand side.

Gauss-Chebyshev quadrature uses weight function $w(x) = 1/\sqrt{1 - x^2}$ on $[-1, 1]$. It's well-suited for integrands of the form $f(x)/\sqrt{1 - x^2}$, where $f$ is smooth. The nodes and weights have closed-form expressions, which is a nice computational advantage. Example: computing $\int_{-1}^{1} \frac{\cos(5x)}{\sqrt{1-x^2}}\,dx$.

Gauss-Hermite quadrature uses $w(x) = e^{-x^2}$ on $(-\infty, \infty)$. This is natural for problems involving Gaussian-type integrands, such as quantum mechanics expectation values. Example: evaluating $\int_{-\infty}^{\infty} x^2 e^{-x^2}\,dx$.

Advanced Techniques and Applications

Gauss-Laguerre quadrature uses $w(x) = e^{-x}$ on $[0, \infty)$, making it appropriate for integrands with exponential decay. Example: $\int_{0}^{\infty} x\, e^{-x}\,dx = 1$.

Gauss-Jacobi quadrature generalizes several of the above. Its weight function is $w(x) = (1-x)^\alpha(1+x)^\beta$ on $[-1, 1]$, with parameters $\alpha > -1$ and $\beta > -1$. Legendre corresponds to $\alpha = \beta = 0$, and Chebyshev to $\alpha = \beta = -1/2$. Adjusting $\alpha$ and $\beta$ lets you handle endpoint singularities of algebraic type directly within the quadrature framework.

Composite Gaussian quadrature subdivides a large interval into smaller subintervals and applies Gaussian quadrature on each. This is useful when:

• The interval is large and the integrand varies significantly across it
• The function has localized features (sharp peaks, near-singularities)

For example, to evaluate $\int_{0}^{10} \sin(x^2)\,dx$, you'd split $[0, 10]$ into subintervals and apply Gauss-Legendre on each, summing the results.

Adaptive Gaussian quadrature takes this further by automatically choosing where to refine. It estimates the error on each subinterval and subdivides only where the error is large. This is especially valuable for integrands like $\int_{0}^{1} x^{-1/2}\,dx$, where the singularity at $x = 0$ demands more nodes nearby but the rest of the interval is smooth.