Deriving Gauss-Hermite weights

approximate integrationnumerical methodsorthogonal-polynomials

I'm currently dealing with Gauss quadrature and I'm having trouble deriving the formula for the Gauss-Hermite quadrature weights. For reference: in my course the Hermite polynomials are defined with the recurrence relation: $$ H_{n+1}(x)=xH_n(x)- \frac{n}{2}H_{n-1}(x)$$
So I set up the corresponding tridiagonal Jacobi Matrix with eigenvalues equal to the roots of $H_n$, defined as:\begin{pmatrix}0&\sqrt{\frac12}&&&\\\sqrt{\frac12}&0&\sqrt{\frac22}&&\\&\sqrt{\frac22}&\ddots&\ddots&\\&&\ddots&\ddots&\sqrt{\frac{n-1}2}\\&&&\sqrt{\frac{n-1}2}&0\end{pmatrix} Now for the weights: I've read that they are supposed to solve the following linear system of equations:$$\sum^n_{i=1}H_k(x_i)*\omega_i=\left\{
\begin{aligned}
(H_0,H_0), for\ k&=0 \\
0, for\ k&=1,2,…,n-1 \\
\end{aligned}
\right. $$ where $(-,-)$ denotes the scalar product of the orthogonal polynomials and $x_i$ are the roots of $H_n$. Here is where I'm unsure how to continue. I know this must be transformed into a linear system of equations but I am unsure how. I've seen a derivation for the Legendre Polynomials which includes finding the eigenvectors of the Jacobi Matrix, but I wouldn't even know to begin if I were forced to find the eigenvectors. Could somebody steer me in the right direction? Thank you in advance

Best Answer

You have started the Galub-Welsch algorithm, noting that the eigenvalues of the tridiagonal array are the zeros that you desire. The weights are equal to the squares entries of the first eigenvector (times the integral of the zeroth Hermite polynomials against the weighting function, which I believe is $\sqrt{\pi}$ by your three term recurrence relation). "First" as defined by the lowest eigenvalue (most negative root).

See Equation (2.6) in https://web.stanford.edu/class/cme335/S0025-5718-69-99647-1.pdf

The algorithm is:

  1. Create Tridiagonal matrix $T$
  2. Find the eigenvectors and eigenvectors $q_{i,j}$ and eigenvalues $x_i$ of $T$
  3. The eigenvalues are the zeros of $H_N(x)$ (as you know)
  4. Sort the eigenvectors by eigenvalues, take the lowest $q_1$ (or highest by symmetry) the weights are equal to the square of the entries multiplied by $(H_0,H_0)$ (in your notation), ie $$w_i = q_{1,i}^2 \int_{-\infty}^\infty H_0(x)^2\omega(x)dx=q_{1,i}^2\sqrt{\pi}$$

Alternatively, what you have written in the sum (they are known as the "tower equations" as per Stoer and Bulirsch) can be solved by using a linear algebra package -- given you have a package that can calculate the Hermite polynomials.

  1. Install a package that has the Hermite polynomials calculated, such as "SpecialPolynomials" in Julia (or equivalently use the Taylor expansion from the wiki)
  2. Construct the matrix $\mathcal{H}_{ij} = H_i(x_j)$ where $i,j\in\{0,N-1\}$
  3. The right hand side is equal to $r^t=(\sqrt{\pi},0,0,0,\cdots)$
  4. Solve $\mathcal{H}\cdot w = r$, ie $w=\mathcal{H} \backslash r$ in Julia
Related Question