Gaussian quadrature for arbitrary weight function except for the Method of Undetermined Coefficients The question is as follow:
For the integration  
$$\int_{-\infty}^\infty w(x)f(x)dx,$$
where $w(x)$ is any form of distribution, if I want to solve the integration with Gaussian-Hermite quadrature by using look-up table method to get the nodes and weights, how should I do it?  
I know the Method of Undetermined Coefficients and Moment-Matching method, but the two methods are troublesome to solve a group of non-linear equation sets. I have ever read that the integration above can be done by the way that  the Rackwitz-Fiessler transformation is first used to transform the variable to standard normal variable and then the nodes and weights are derived by using Gaussian-Hermite quadrature. However, I don't know the specific procedures in it.  
I really hope somebody can help me out. Thanks!
 A: The abbreviated explanation for implementation of Gaussian quadrature looks like this. I'll define:


*

*the inner product $(f,g)=\int_{-\infty}^\infty w(x) f(x) g(x) dx$ on suitable functions $f,g \in L^2_w$ (if you don't know what $L^2_w$ means, ignore it). 

*the corresponding norm $\| f \|=(f,f)^{1/2}$.


Next I'll define two sequences of polynomials $p_k,q_k$ with


*

*$p_k=q_k/\| q_k \|$

*$q_0=1$ 

*$q_1=x-(xp_0,p_0)p_0$

*$q_k=xp_{k-1}-(xp_{k-1},p_{k-1})p_{k-1}-(xp_{k-1},p_{k-2})p_{k-2},k \geq 2$
This last equation is called a "three term recurrence relation". In this setting, $p_k$ has degree $k$, has norm $1$, and is orthogonal to all the other $p_j$. 
Define $a_k=(xp_k,p_k)$ for $k \geq 0$, and $b_k=(xp_k,p_{k-1})$ for $k \geq 1$. (Note that these are exactly the integrals computed to generate the $p_k$.) Then we define the $n \times n$ matrix $T_{ij}$ with $T_{ii}=a_{i-1}$ for $i \geq 1$ and $T_{i-1,i}=T_{i,i-1}=(b_{i-1})^{1/2}$ for $i \geq 2$ and all other $T_{ij}=0$. Finally we define $c=(1,1)$ or in other words $c=\int_{-\infty}^\infty w(x) dx$. ($c$ had to come in somewhere, because the $p_k$ remain invariant under multiplication of $w$ by a constant.)
Then the nodes $x_i$ of the $n$ point Gaussian quadrature method are the eigenvalues of $T$. The weight $w_i = c (y^{(i)}_1)^2$ where $y^{(i)}$ is a unit eigenvector with eigenvalue $x_i$. Since this is a symmetric tridiagonal eigenproblem, it is pretty easy to solve, relative to its size. The task of computing the $a_i,b_i$ will probably take more time than solving the eigenproblem (unless they can be computed in closed form, which they can for certain simple $w$).
This procedure is called the Golub-Welsch algorithm. It works primarily because of the fundamental theorem of Gaussian quadrature, which says that the nodes of a $n$ point Gaussian quadrature are the roots of $p_n$. You can find more details at https://en.wikipedia.org/wiki/Gaussian_quadrature
I'll add that technically this is equivalent to "undetermined coefficients", i.e. to setting up the nonlinear system $\sum_{i=1}^n w_i x_i^k = \int_{-\infty}^\infty w(x) x^k dx$ for $k=0,1,\dots,2n-1$. But in practice Golub-Welsch will be much easier than solving that nonlinear system. That is because:


*

*the nonlinear system is poorly scaled. (It would be well-scaled if you used $p_k$ instead of $x^k$.)

*Casting the nonlinear system into a symmetric tridiagonal eigenproblem makes it far easier.

