Gauss–Legendre quadrature

In numerical analysis, Gauss–Legendre quadrature is a form of Gaussian quadrature for approximating the definite integral of a function. For integrating over the interval [−1, 1], the rule takes the form:

${\displaystyle \int _{-1}^{1}f(x)\,dx\approx \sum _{i=1}^{n}w_{i}f(x_{i})}$

where

• n is the number of sample points used,
• w_i are quadrature weights, and
• x_i are the roots of the nth Legendre polynomial.

This choice of quadrature weights w_i and quadrature nodes x_i is the unique choice that allows the quadrature rule to integrate degree 2n − 1 polynomials exactly.

Many algorithms have been developed for computing Gauss–Legendre quadrature rules. The Golub–Welsch algorithm presented in 1969 reduces the computation of the nodes and weights to an eigenvalue problem which is solved by the QR algorithm.^[1] This algorithm was popular, but significantly more efficient algorithms exist. Algorithms based on the Newton–Raphson method are able to compute quadrature rules for significantly larger problem sizes. In 2014, Ignace Bogaert presented explicit asymptotic formulas for the Gauss–Legendre quadrature weights and nodes, which are accurate to within double-precision machine epsilon for any choice of n ≥ 21.^[2] This allows for computation of nodes and weights for values of n exceeding one billion.^[3]

History

Carl Friedrich Gauss was the first to derive the Gauss–Legendre quadrature rule, doing so by a calculation with continued fractions in 1814.^[4] He calculated the nodes and weights to 16 digits up to order n=7 by hand. Carl Gustav Jacob Jacobi discovered the connection between the quadrature rule and the orthogonal family of Legendre polynomials. As there is no closed-form formula for the quadrature weights and nodes, for many decades people were only able to hand-compute them for small n, and computing the quadrature had to be done by referencing tables containing the weight and node values. By 1942 these values were only known for up to n=16.

Definition

Graphs of Legendre polynomials (up to n = 5)

For integrating f over ${\displaystyle [-1,1]}$ with Gauss–Legendre quadrature, the associated orthogonal polynomials are Legendre polynomials, denoted by P_n(x). With the n-th polynomial normalized so that P_n(1) = 1, the i-th Gauss node, x_i, is the i-th root of P_n and the weights are given by the formula^[5]

${\displaystyle w_{i}={\frac {2}{\left(1-x_{i}^{2}\right)\left[P'_{n}(x_{i})\right]^{2}}}.}$

Some low-order quadrature rules are tabulated below for integrating over ${\displaystyle [-1,1]}$.

Number of points, n	Points, xi		Weights, wi
1	0		2
2	${\displaystyle \pm {\frac {1}{\sqrt {3}}}}$	±0.57735…	1
3	0		${\displaystyle {\frac {8}{9}}}$	0.888889…
	${\displaystyle \pm {\sqrt {\frac {3}{5}}}}$	±0.774597…	${\displaystyle {\frac {5}{9}}}$	0.555556…
4	${\displaystyle \pm {\sqrt {{\frac {3}{7}}-{\frac {2}{7}}{\sqrt {\frac {6}{5}}}}}}$	±0.339981…	${\displaystyle {\frac {18+{\sqrt {30}}}{36}}}$	0.652145…
	${\displaystyle \pm {\sqrt {{\frac {3}{7}}+{\frac {2}{7}}{\sqrt {\frac {6}{5}}}}}}$	±0.861136…	${\displaystyle {\frac {18-{\sqrt {30}}}{36}}}$	0.347855…
5	0		${\displaystyle {\frac {128}{225}}}$	0.568889…
	${\displaystyle \pm {\frac {1}{3}}{\sqrt {5-2{\sqrt {\frac {10}{7}}}}}}$	±0.538469…	${\displaystyle {\frac {322+13{\sqrt {70}}}{900}}}$	0.478629…
	${\displaystyle \pm {\frac {1}{3}}{\sqrt {5+2{\sqrt {\frac {10}{7}}}}}}$	±0.90618…	${\displaystyle {\frac {322-13{\sqrt {70}}}{900}}}$	0.236927…

For integrating over a general real interval ${\displaystyle [a,b]}$, a change of interval can be applied to convert the problem to one of integrating over ${\displaystyle [-1,1]}$.

Algorithms

Newton–Raphson methods

Several researchers have developed algorithms for computing Gauss–Legendre quadrature nodes and weights based on the Newton–Raphson method for finding roots of functions. Various optimizations for this specific problem are used. For instance, some methods compute ${\displaystyle \theta _{i}=\arccos x_{i}}$ to avoid issues associated with clustering of the roots ${\displaystyle x_{i}}$ near the ends of the interval ${\displaystyle -1}$ and ${\displaystyle 1}$.^[6][7] Some methods utilize formulas to approximate the weights and then use a few iterations of Newton-Raphson to lower the error of the approximation to below machine precision.^[8][6]

Golub–Welsch

In 1969, Golub and Welsch published their method for computing Gaussian quadrature rules given the three term recurrence relation that the underlying orthogonal polynomials satisfy.^[1] They reduce the problem of computing the nodes of a Gaussian quadrature rule to the problem of finding the eigenvalues of a particular symmetric tridiagonal matrix. The QR algorithm is used to find the eigenvalues of this matrix. By taking advantage of the symmetric tridiagonal structure, the eigenvalues can be computed in ${\displaystyle {\mathcal {O}}(n^{2})}$ time, as opposed to the ${\displaystyle {\mathcal {O}}(n^{3})}$ time expected for a generic eigenvalue problem.

Asymptotic formulas

Various methods have been developed that use approximate closed-form expressions to compute the nodes. As mentioned above, in some methods formulas are used as approximations to the nodes, after which some Newton-Raphson iterations are performed to refine the approximation. In a 2014 paper, Ignace Bogaert derives asymptotic formulas for the nodes that are exact up to machine precision for ${\displaystyle n\geq 21}$ and for the weights that are exact up to machine precision for ${\displaystyle n\geq 30}$.^[2] This method does not require any Newton-Raphson iterations or evaluations of Bessel functions as other methods do. As shown in the paper, the method was able to compute the nodes at a problem size of one billion in 11 seconds. Since the nodes near the endpoints of ${\displaystyle [-1,1]}$ become very close to each other at this problem size, this method of computing the nodes is sufficient for essentially any practical application in double-precision floating point.

Higher precision

Johansson and Mezzarobba^[9] describe a strategy to compute Gauss–Legendre quadrature rules in arbitrary-precision arithmetic, allowing both small and large ${\displaystyle n}$. A rule of order ${\displaystyle n=1000}$ with 1000 digits of precision can be calculated in around one second. Their method uses Newton–Raphson iteration together with several different techniques for evaluating Legendre polynomials. The algorithm also provides a certified error bound.

Gil, Segura and Temme^[10] describe iterative methods with fourth order convergence for the computation of Gauss–Jacobi quadratures (and, in particular, Gauss–Legendre). The methods do not require a priori estimations of the nodes to guarantee its fourth-order convergence. Computations of quadrature rules with even millions of nodes and thousands of digits become possible in a typical laptop.

Comparison with other quadrature rules

Gauss–Legendre quadrature is optimal in a very narrow sense for computing integrals of a function f over [−1, 1], since no other quadrature rule integrates all degree 2n − 1 polynomials exactly when using n sample points. However, this measure of accuracy is not generally a very useful one---polynomials are very simple to integrate and this argument does not by itself guarantee better accuracy on integrating other functions.

Clenshaw–Curtis quadrature is based on approximating f by a polynomial interpolant at Chebyshev nodes and integrates polynomials of degree up to n exactly when given n samples. Even though it does not integrate polynomials or other functions that are analytic in a large neighborhood of [−1, 1] as well as Gauss–Legendre quadrature, Clenshaw–Curtis converges at approximately the same rate as Gauss–Legendre quadrature for most (non-analytic) integrands.^[11] Also, Clenshaw–Curtis shares the properties that Gauss–Legendre quadrature enjoys of convergence for all continuous integrands f and robustness to numerical rounding errors.^[12] Clenshaw–Curtis is straightforward to implement in ${\displaystyle {\mathcal {O}}(n\log n)}$ time by FFT-based methods.

Newton–Cotes quadrature is based on approximating f by a polynomial interpolant at equally-spaced points in [−1, 1], and like Clenshaw–Curtis also integrates polynomials of degree up to n exactly when given n samples. However, it suffers from Runge's phenomenon as n increases; Newton–Cotes does not converge for some continuous integrands f, and when it does converge it may suffer from numerical rounding errors.^[12] Thus, both Clenshaw–Curtis and Gauss–Legendre enjoy substantial benefits over Newton–Cotes for moderate to large n.

References

1. Golub, Gene H.; Welsch, John H. (1969). "Calculation of Gauss quadrature rules". Math. Comp. 23 (106): 221–230. doi:10.1090/S0025-5718-69-99647-1. JSTOR 2004418.
2. Bogaert, I. (2014). "Iteration-free computation of Gauss–Legendre quadrature nodes and weights". SIAM J. Sci. Comput. 36 (3): C1008–C1026. doi:10.1137/140954969. hdl:1854/LU-5683230. S2CID 690538.
3. Townsend, A. (2015). "The race for high order Gauss–Legendre quadrature". SIAM News. 48: 1–3.
4. C. F. Gauss, Methodus nova integralium valores per approximationem inveniendi, Comment. Soc. Reg. Scient. Gotting. Recent., (1814).
5. (Abramowitz & Stegun 1983, p. 887)
6. Hale, N.; Townsend, A. (2013). "Fast and accurate computation of Gauss–Legendre and Gauss–Jacobi quadrature nodes and weights". SIAM J. Sci. Comput. 35 (2): A652–A674. doi:10.1137/120889873.
7. Swarztrauber, P. N. (2003). "On computing the points and weights for Gauss–Legendre quadrature". SIAM J. Sci. Comput. 24: 945–954. doi:10.1137/S2064827500379690 (inactive 2024-10-15).{{cite journal}}: CS1 maint: DOI inactive as of October 2024 (link)
8. I. Bogaert, B. Michiels, and J. Fostier, O(1) computation of Legendre polynomials and Gauss–Legendre nodes and weights for parallel computing, SIAM J. Sci. Comput., 34 (2012), pp. 83–101.
9. Johansson, F.; Mezzarobba, M. (2018). "Fast and rigorous arbitrary-precision computation of Gauss–Legendre quadrature nodes and weights". SIAM J. Sci. Comput. 40 (6): C726–C747. arXiv:1802.03948. doi:10.1137/18M1170133.
10. Gil, A.; Segura, J.; Temme, M. M. (2021). "Fast and reliable high accuracy computation of Gauss–Jacobi quadrature". Numer. Algorithms. 87: 1391–1419. arXiv:2008.08641. doi:10.1007/s00211-019-01066-2. S2CID 189762478.
11. Lloyd N. Trefethen. 2012. Approximation Theory and Approximation Practice. Society for Industrial and Applied Mathematics, USA.
12. Trefethen, L.N. (2008). "Is Gauss Quadrature Better than Clenshaw–Curtis?". SIAM Rev. 50 (1): 67–87. doi:10.1137/060659831.

Sources

• Abramowitz, Milton; Stegun, Irene Ann, eds. (1983) [June 1964]. "Chapter 25.4, Integration". Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Applied Mathematics Series. Vol. 55 (Ninth reprint with additional corrections of tenth original printing with corrections (December 1972); first ed.). Washington D.C.; New York: United States Department of Commerce, National Bureau of Standards; Dover Publications. ISBN 978-0-486-61272-0. LCCN 64-60036. MR 0167642. LCCN 65-12253.