Gauss-Legendre

Motivation

Instead of evaluating the integrand at equally spaced nodes as Newton-Cotes methods do, Gauss-Legendre quadrature places nodes at the zeros of Legendre polynomials — positions that are optimal in the sense that an \(n\)-point rule integrates all polynomials of degree up to \(2n-1\) exactly. This means it achieves the same accuracy as a \(2n\)-point trapezoidal rule with only half the function evaluations, making it the preferred method for smooth integrands on finite intervals.

Example

use integrate::gauss_quadrature::legendre_rule;

let square = |x: f64| x * x;

let a = 0.0;
let b = 1.0;

let n: usize = 100;

let integral = legendre_rule(square, a, b, n);
println!("{}", integral);

Understanding Gauss-Legendre rule

Legendre polynomials

With respect to the inner product

\[ \langle f, g \rangle = \int_{-1}^{1} f(x) \cdot g(x) \cdot w(x) \, dx \]

the Legendre polynomials, defined by Rodrigues’ formula

\[ P_n(x) = \frac{1}{2^n n!} \frac{\partial^{n} (x^2 - 1)^n}{\partial x^n} \quad \text{for} \quad n > 0 \]

and \(P_0(x) = 1\), form an orthogonal family with weight function \(w(x) = 1\) on \([-1, 1]\).

Quadrature formula on \([-1, 1]\)

The \(n\)-point Gauss-Legendre quadrature formula \(GL_n(f)\), for approximating \(\int_{-1}^{1} f(x)\,dx\), is given by

\[ GL_n(f) = A_1 f(x_1) + \cdots + A_n f(x_n) \]

where \(x_1, \ldots, x_n\) are the zeros of \(P_n\) and the weights are

\[ A_i = \frac{2}{1 - x_i^2} \cdot \frac{1}{n^2 P_{n-1}(x_i)^2} \quad \text{for} \quad i = 1, \ldots, n \]

All weights \(A_i\) are positive, which makes the rule numerically stable.

Change of variables to \([a, b]\)

To integrate over a general interval \([a, b]\), the substitution \(t = \dfrac{2(x - a)}{b - a} - 1\) gives

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

so the \(n\)-point Gauss-Legendre approximation on \([a, b]\) is

\[ GL_n(f, a, b) = A_1’ f(x_1’) + \cdots + A_n’ f(x_n’) \]

where the mapped nodes and weights are

\[ x_i’ = \frac{b-a}{2} \cdot x_i + \frac{b+a}{2}, \qquad A_i’ = \frac{b-a}{2} \cdot A_i \quad \text{for} \quad i = 1, \ldots, n \]

Truncation error

The truncation error for the rule on \([-1, 1]\) is

\[ \int_{-1}^{1} f(x) \, dx - GL_n(f) = K \cdot \frac{f^{(2n)}(c)}{(2n)!} \]

where \(c \in (-1, 1)\) is unknown and \(K\) is a constant determined by

\[ K = \int_{-1}^{1} x^{2n} \, dx - GL_n(x^{2n}) \]

The same form holds on \([a, b]\) with \(c \in (a, b)\). A corollary is that if \(f^{(2n)}(x) = 0\) for all \(x \in [a, b]\) — i.e. if \(f\) is a polynomial of degree at most \(2n - 1\) — then the rule is exact.

Node and weight computation

Most Gauss quadrature rules compute nodes and weights by building a tridiagonal Jacobi matrix and finding its eigenvalues (Golub-Welsch algorithm), which costs \(O(n^2)\) in total. The Legendre rule in this crate uses a different strategy that delivers each node and weight independently in \(O(1)\) time, regardless of \(n\):

• \(n \leq 100\): each \((x_i, A_i)\) pair is read directly from a precomputed table of \(\theta\)-coordinates and weights stored in the source. The lookup is a single array index with no arithmetic at all.

• \(n > 100\): Bogaert’s asymptotic formula expresses each node as \(x_i = \cos\theta_i\), where \(\theta_i\) is obtained from the \(k\)-th zero of the Bessel function \(J_0\) via a fixed-degree Chebyshev polynomial correction (6 terms for the node, 8 terms for the weight). The entire computation is a handful of fused-multiply-add operations — no Newton iteration, no eigenvalue solver.

The result is that the per-node cost is \(O(1)\) and the error on every node and weight is within a few ulps of the true value, even for very large \(n\).

Limitations

Gauss-Legendre quadrature is best suited for smooth functions on finite intervals. It is not appropriate for functions with singularities, discontinuities, or functions defined on unbounded domains. Rapidly oscillating functions may also require a very large number of nodes \(n\). For semi-infinite or doubly-infinite domains use the Gauss-Laguerre or Gauss-Hermite rules instead.

References

Node and weight computation based on the original C++ implementation by Ignace Bogaert, featuring \(O(1)\) per-node complexity and errors within a few ulps:

\[ [1]: \text{Bogaert I., } \textit{Iteration-free computation of Gauss-Legendre quadrature nodes and weights,} \\ \text{SIAM Journal on Scientific Computing, Volume 36, Number 3, 2014, pages A1008\text{–}A1026.} \]

For more details, here is a link to the article.