6
$\begingroup$

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.

Improve this question
$\endgroup$
1

2 Answers 2

Reset to default
3
$\begingroup$

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

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks, you are right, I've corrected the formula. I'll sit down and read the references. $\endgroup$
    – Amir Sagiv
    Apr 13, 2016 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, 2016 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, 2016 at 14:11
  • $\begingroup$ Yep, this is exactly what I was looking for. Thanks! $\endgroup$
    – Amir Sagiv
    Apr 13, 2016 at 18:36
1
$\begingroup$

You are describing the Golub-Welsch algorithm. As described, it costs O(n^2) operations to compute n-point Gauss-Legendre quadrature rule.
http://www.ams.org/journals/mcom/1969-23-106/S0025-5718-69-99647-1/

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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.