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.

  • I know I already got some answers, but I'm really intrigued whether or not anyone understands the solution I've quoted in the "edit" section – Amir Sagiv May 23 '15 at 9:37
up vote 4 down vote accepted

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 http://sourceforge.net/projects/fastgausslegendrequadrature/.

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

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)$).

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.