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.
Best Answer
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).