On Computing the Points and Weights for Gauss--Legendre Quadrature

Three methods are compared for computing n points x i and weights w(x i ) for Gauss--Legendre quadrature. The eigensystem method computes the points as the eigenvalues of a symmetric tridiagonal matrix, and the weights are computed from its eigenvectors. The resulting error is O(n) for the points and O(n²) for the weights. Using an inverse Taylor series computation for the points, Yakimiw [J. Comput. Phys., 129 (1996), pp. 406--430] reduced the point error growth to O(1). He also reduced the error in the weights to O(n) and obtained uniform pointwise relative accuracy. However, the points x i cluster quadratically near the end points x = ± 1, where for some n = n max they coalesce or even exchange to within machine accuracy, and computation fails. Although this has not been a problem in the past, it could soon become so because weather models are now approaching resolutions with n > 1000. Here we compare these methods to the Fourier--Newton method, which consists of Newton's method for computing the transformed points θᵢ = xᵢ as the zeros of the Fourier representations of the Legendre polynomials Pₙ(θ). The θᵢ are near equally spaced, particularly in the neighborhood of x = 0. For a fixed arithmetic precision, this extends the range of computable points and weights to n=nₘₐₓ², where n max is the limit for methods based on x i . Both the accuracy and speed of the three methods are tabulated.