Computing Gauss Legendre quadrature for large $N$

Question

I've been scanning across the web, and haven't found a good method to compute the Gauss Legendre abscissas and weights $\{ x_j, w^j \} _{j=1}^N$ for large $N\in\mathbb{N}$. My question is how to do it, and why should it work?

To those who need some background:

The goal is to approximate an integral by a discrete interpolating sum:

$$\int\limits_{(-1,1)} f(x) dx = \sum\limits_{j=-N}^N f(x_j)w^j $$

The question is, how to choose $\{ x_j, w^j \} _j$ appropriately.The Gauss Legendre quadrature tells you (for good reasons) to choose $x_j $ to be the roots of the $n$-th Legendre polynomial.

Problem : A straightforward computation of the Legendre polynomial for high $N$ is highly unstable, as it involves "big" coefficients of alternating signs.

EDIT: I've found a very simple code that computes weight and abscissas using eigenvalues of a symmetric matrix, but doesn't seem use Golub Welsch. The matrix is

$$\forall 1\leq j\leq N-1 \, A_{j,j+1} = A_{j+1,j} = \frac{j}{\sqrt{4j^2 -1}}$$ with all other entries are zero. The discussion about it was split to another post.

Answer 1

There are asymptotic methods that essentially give you $N$ nodes and weights in $O(N)$ time if the precision is assumed to be fixed (e.g. at double precision).

See Nicholas Hale and Alex Townsend, "Fast and Accurate Computation of Gauss-Legendre and Gauss-Jacobi Quadrature Nodes and Weights", SIAM J. Sci. Comput., 35(2) (a PDF is available at http://eprints.maths.ox.ac.uk/1629/1/finalOR79.pdf, Wayback Machine).

They claim that their algorithm achieves double precision accuracy for $N \ge 100$.

For $N < 100$, you may as well precompute a big table with perfect accuracy using a computer algebra system or arbitrary precision library of your choice (or look up tables that others have published).

As to computing Legendre polynomials in a numerically stable way, use the three-term recurrence $(n+1) P_{n+1}(x) = (2n+1) x P_n(x) - n P_{n-1}(x)$ to evaluate $P(x)$ directly instead of computing the coefficients of the polynomial and using Horner's rule (similarly for $P'(x)$).

Update (2019): the issue of efficient computation for large N with arbitrary precision (and also with rigorous error bounds) is addressed in new work by myself and Marc Mezzarobba: SIAM Journal on Scientific Computing, 2018, Vol. 40, No. 6 : pp. C726-C747 https://doi.org/10.1137/18M1170133 (https://arxiv.org/abs/1802.03948).

Answer 2

To add to Fredrik Johansson's answer: A nice history of algorithms for computing Gauss quadrature rules can be found in this SIAM News article by Alex Townsend. Therein, it is stated that the "final chapter" was written by Ignace Bogaert in this SISC paper, which gives an algorithm that is even faster and more accurate than the algorithm of Hale & Townsend.

There is a free open-source implementation at https://sourceforge.net/projects/fastgausslegendrequadrature/