Editing Gaussian quadrature (section)

==== The Golub-Welsch algorithm ====
The three-term recurrence relation can be written in matrix form <math>J\tilde{P} = x\tilde{P} - p_n(x) \mathbf{e}_n</math> where <math>\tilde{P} = \begin{bmatrix} p_0(x) & p_1(x) & \cdots & p_{n-1}(x) \end{bmatrix}^\mathsf{T}</math>, <math>\mathbf{e}_n</math> is the <math>n</math>th standard basis vector, i.e., <math>\mathbf{e}_n = \begin{bmatrix} 0 & \cdots & 0  & 1 \end{bmatrix}^\mathsf{T}</math>, and {{mvar|J}} is the following [[tridiagonal matrix]], called the Jacobi matrix:
<math display="block">\mathbf{J} = \begin{bmatrix}
       a_0 &      1 &      0 &  \cdots  &       0 \\
       b_1 &    a_1 &      1 &  \ddots  &  \vdots \\
         0 &    b_2 & \ddots &  \ddots  &       0 \\
    \vdots & \ddots & \ddots &  a_{n-2} &       1 \\
         0 & \cdots &       0 & b_{n-1} & a_{n-1}
\end{bmatrix}.</math>

The zeros <math>x_j</math> of the polynomials up to degree {{mvar|n}}, which are used as nodes for the Gaussian quadrature can be found by computing the eigenvalues of this matrix. This procedure is known as ''Golub–Welsch algorithm''.

For computing the weights and nodes, it is preferable to consider the [[Symmetric matrix|symmetric]] tridiagonal matrix <math>\mathcal{J}</math> with elements
<math display="block">\begin{align}
                                            \mathcal{J}_{k,i} = J_{k,i} &= a_{k-1}        & k &= 1,2,\ldots,n \\[2.1ex]
  \mathcal{J}_{k-1,i} = \mathcal{J}_{k,k-1} = \sqrt{J_{k,k-1}J_{k-1,k}} &= \sqrt{b_{k-1}} & k &= \hphantom{1,\,}2,\ldots,n.
\end{align}</math>

That is,

<math display="block">\mathcal{J} = \begin{bmatrix}
       a_0 & \sqrt{b_1} &          0 &        \cdots  &              0 \\
\sqrt{b_1} &        a_1 & \sqrt{b_2} &        \ddots  &         \vdots \\
         0 & \sqrt{b_2} &     \ddots &        \ddots  &              0 \\
    \vdots &     \ddots &     \ddots &        a_{n-2} & \sqrt{b_{n-1}} \\
         0 &     \cdots &          0 & \sqrt{b_{n-1}} &        a_{n-1}
\end{bmatrix}.</math>

{{math|'''J'''}} and <math>\mathcal{J}</math> are [[similar matrices]] and therefore have the same eigenvalues (the nodes). The weights can be computed from the corresponding eigenvectors: If <math>\phi^{(j)}</math> is a normalized eigenvector (i.e., an eigenvector with euclidean norm equal to one) associated with the eigenvalue {{mvar|x<sub>j</sub>}}, the corresponding weight can be computed from the first component of this eigenvector, namely:
<math display="block">w_j = \mu_0 \left(\phi_1^{(j)}\right)^2</math>

where <math>\mu_0</math> is the integral of the weight function
<math display="block">\mu_0 = \int_a^b \omega(x) dx.</math>

See, for instance, {{harv|Gil|Segura|Temme|2007}} for further details.