I have a problem involving numerical integration of the form:
$$I = \int_0^\infty \!dx \, w(x) f(x)$$
where the weighting function is a chi distribution of degree 2, i.e.,
$$w(x) = x \, e^{\frac{-x^2}{2}}$$
Does someone know the relevant orthogonal polynomial or how it could be derived? Is there a general method for deriving orthogonal polynomials for arbitrary weights and intervals of integration?
Thanks!
The theory is explained quite well in the Wikipedia article on Gaussian quadrature. You need to read the section on the general formula for the weights and on the recurrence relation to generate the orthogonal polynomials. I personally prefer to use a different method to derive quadrature rules which can more easily deal with additional information such as information about function values or (higher order) derivatives at specific points or known integral evaluations of the function multiplied by other weight functions or infinite summations of the function etc.. This works as follows.
Cauchy's integral formula allows us to write an integral of an analytic function $f(x)$ multiplied by an analytic weight function $w(x)$ as
$$\int_a^b w(x) f(x) dx = \frac{1}{2\pi i}\oint_C f(z) g(z)dz \tag{1}$$
where $C$ is a counterclockwise contour encircling the interval $[a,b]$, and
$$g(z) = \int_a^b \frac{w(x) dx}{z-x} \tag{2}$$
Expanding $g(z)$ around $z = \infty$ yields:
$$g(z) = \sum_{k=0}^{\infty}\frac{c_k}{z^{k+1}}\tag{3}$$
where the expansion coefficients $c_k$ are the moments of the weight function:
$$c_k = \int_a^b w(x) x^k dx$$
Conventional quadrature rules can then be obtained from (3) using the Padé approximant method. This allows us to find a rational function $Q_n(z)$ with numerator and denominator of degree $n$ in $\dfrac{1}{z}$ that matches the series (3) up to order $2n$. The function $Q_n(z)$ will in general have $n$ simple poles and can therefore be written as:
$$Q_n(z) = \sum_{j=1}^n \frac{w_j}{z-x_j}\tag{4}$$
If we approximate $g(z)$ in (1) by this expression and apply Cauchy's integral formula, we get a quadrature rule:
$$\int_a^b w(x) f(x) dx \approx \sum_{j=1}^n w_j f(x_j)$$
Since for large $z$, $g(z)$ differs from $Q(z)$ by terms of order $2n+1$ or larger in $\dfrac{1}{z}$, if follows that the quadrature rue is exact for $f(z)$ a polynomial of order $2n-1$, as $f(z) g(z)$ and $f(z) Q_n(z)$ will then differ from each other by terms of order $2$ or higher in $\dfrac{1}{z}$ which don't contribute to the contour integral.
For $w(x) = x \exp\left(-\frac{1}{2} x^2\right)$ we have:
$$c_k = \int_0^\infty x^{k+1} \exp\left(-\frac{1}{2} x^2\right)dx = 2^{\frac{k}{2}}\left(\frac{k}{2}\right)! $$
To get a 10 point quadrature rule we need to construct the [10/10] Padé approximant to $\sum_{k=0}^{19}\dfrac{c_k}{z^{k+1}}$. This takes just a few seconds of typing commands in Mathematica. The abscissa are the zeros of the denominator, they are at the points
$$ \begin{split} x_1 &= 0.12851261647070163641977064064860341108585877838908\\ x_2 &= 0.41572412829781859808211088609524755876388678431860\\ x_3 &= 0.83285826988445711612170229447704512679103176241683 \\ x_4 &= 1.35203760832396670771080034445078357894777166770144 \\ x_5 &= 1.95247254854000705863961082251262520444757933349651 \\ x_6 &= 2.62247785450006256126816957605308062657357748919288 \\ x_7 &= 3.35981743051373842409516754701476134720041887362005\\ x_8 &= 4.17354153848492249357376499539115142382485188180738 \\ x_9 &= 5.09226305313159606471611471587701418927804036939713\\ x_{10} &= 6.20248114919776450555405845591268240887766132895859 \end{split} $$
The weights are the residues of the Padé approximant at these poles:
$$ \begin{split} w_1 &= 0.02706653728875616663430922411999220730385427989898\\ w_2 &= 0.13619130027820014133524095926466149691350609506380\\ w_3 &= 0.27810052085034958748288747332588952454574290605869 \\ w_4 &= 0.3049206954472171684651353227872738721713569641388 \\ w_5 &= 0.1847200419742236837090651175552161912003223283323 \\ w_6 &= 0.059210304776936908820384872373088408001861428585 \\ w_7 &= 0.009187069314768842373513866382962725029786819745\\ w_8 &= 0.000591695015446369243719887139810902494014093397 \\ w_9 &= 0.0000117998889029225848877292645219561178319426849\\ w_{10} &=3.516519820935085554778658271622172314209555\times 10^{-8} \end{split} $$