While searching the web, I came across the following algorithm for the Gauss-Legendre quadrature. I wasn't able to find a reference or a proof of my own as to why it works. I'll present it, and the question is - Why does it work?

The Method:

  1. We construct the Matrix $A \in M_{N\times N} (\mathbb{R} )$ by $\forall 1\leq j\leq N-1 \, A_{j,j+1} = A_{j+1,j} = \frac{j}{\sqrt{4j^2 -1}}$, and all other entries are zero.

  2. For some reason, the Legendre polynomial of order $N$ is the characteristic polynomial $p_N(x) = {\rm det}(x{\rm Id - A})$. Therefore, the eigenvalues of $A$ are the abscissas of the quadrature.

  3. For the weights, we take $v^j$ to be the eigenvector of the $j-th$ smallest (negative) eigenvalue, and we have that $w_j =2(v^j _1) ^2, \, \forall 1\leq j \leq N$, with $v^j_1$ denoting the first coordinate.

I don't understand this method, nor can I find a proper reference to it in standard Numerical Integration books (Davis and Rebinowitz) or Orthogonal Polynomials books (Szego).

Meta - Note: in my old question I asked for a stable and efficient algorithm for the Gauß-Legendre quadrature, but none of the proposed methods was the one I've found. This is why I've opened a new question.


This is a particular implementation of a more general method, described in John Boyd's Why Eigenvalues Are Roots: A Derivation of the One-Dimensional Companion Matrix for General Orthogonal Polynomials (restricted access). Gauss-Legendre quadrature of order $n$ needs the roots of the Legendre polynomial $P_n(x)$ and a numerical root solver must guarantee that the roots are real. By reformulating the root finding problem into an eigenvalue problem for a symmetric matrix a numerical instability for large $n$ is avoided. The specific choice of symmetric companion matrix used here is derived on page 9 of this paper.

For the weights, I think the method described in the OP needs correction: It is not the $j$-th component of the eigenvector of the smallest eigenvalue that determines the weight $w_j$, but the first component of the $j$-th eigenvector. (I looked at the code linked by the OP, and it seems that is indeed what it does.) The relation is explained on page 223 of the Golub-Welsch paper, where all of this originates (or in more detail on page 6 of these notes).

  • $\begingroup$ Thanks, you are right, I've corrected the formula. I'll sit down and read the references. $\endgroup$ – Amir Sagiv Apr 13 '16 at 13:16
  • $\begingroup$ I can't seem to access Boyd's book or have it in our library. Is it available anywhere else? Is the same material covered somewhere else? $\endgroup$ – Amir Sagiv Apr 13 '16 at 13:32
  • $\begingroup$ added a reference that derives the specific companion matrix used in the code, I think that paper is open access. $\endgroup$ – Carlo Beenakker Apr 13 '16 at 14:11
  • $\begingroup$ Yep, this is exactly what I was looking for. Thanks! $\endgroup$ – Amir Sagiv Apr 13 '16 at 18:36

You are describing the Golub-Welsch algorithm. As described, it costs O(n^2) operations to compute n-point Gauss-Legendre quadrature rule.

If you are interested in large Gauss-Legendre quadrature rule, then better (faster and more accurate) algorithms now exist. See your previous question: Computing Gauss Legendre Quadrature for Large N


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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