# numerical integration with possible division by 'zero'

I am trying to integrate

$$\int^1_0 t^{2n+2}\exp\left({\frac{\alpha r_0}{t}}\right)dt$$

which is a simple transformation of

$$\int^{\infty}_1 x^{2n}\exp(-\alpha r_0 x)dx$$

using $t = \frac1{x}$ because it is difficult to numerically approximate improper integrals. This does, however, lead to the problem of evaluating the new integrand near zero. It will be very easy to get the proper number of quadrature nodes seeing as the interval is only of length 1 (so the comparable $dt$ can be made very small), but what sort of considerations should I make when integrating near zero?

On some level, I think that simply taking $\int^1_\epsilon t^{2n+2}\exp({\frac{\alpha r_0} {t}})dt$ is a good idea where $\epsilon$ is some small number. However, what number should I choose? Should it be machine epsilon? Is division by machine epsilon a well quantified number? Furthermore, if division my machine epsilon (or close to it) gives an incredibly large number, then taking $\exp(\frac{1}{\epsilon})$ will become even larger.

How should I account for this? Is there a way to have a well defined numerical integral of this function? If not, what is the best way of integrating the function?

• Have you considered using Monte Carlo? Dec 13, 2011 at 17:15
• I feel like it would not fix the problem. Monte Carlo integration is often reserved for high dimensional integrals. I would run into the exact same problems with Monte Carlo, I would simply have less control over where my function is being evaluated. Dec 13, 2011 at 17:16
• You might be right. Dec 13, 2011 at 17:21
• I think it would still be good to have an answer (perhaps to a separate, more general question) explaining how one does numerical integration when the function is divergent at one limit, for the general case where it's not possible to do the integral analytically. Then again, that could just as well be found in Numerical Recipes... Dec 13, 2011 at 19:03
• @Faheem: "Monte Carlo is an extremely bad method; it should be used only when all alternative methods are worse." - Alan Sokal Dec 14, 2011 at 7:31

This can be done by integration by parts: $$\int^\infty_1 x e^{-ax} = \frac{-1}{a} x e^{-ax}\mid^\infty_1 - \frac{-1}{a} \int^\infty_1 e^{-ax} = \frac{e^{-a}}{a} + \frac{e^{-a}}{a^2} = \frac{a+1}{a^2} e^{-a}$$ and continuing on by induction $$\int^\infty_1 x^k e^{-ax} = \frac{-1}{a} x^k e^{-ax}\mid^\infty_1 - \frac{-k}{a} \int^\infty_1 x^{k-1} e^{-ax} = \frac{e^{-a}}{a} + \frac{k}{a} \int^\infty_1 x^{k-1} e^{-ax}$$ so that $$I(k) = \frac{e^{-a}}{a} + \frac{k}{a} I(k-1)$$ and $I(0) = \frac{e^{-a}}{a}$.
• It's often a good idea to ask a computer algebra system whenever you have an integral like this. Maple evaluates "$\int_1^\infty x^{2n} \exp(-\alpha x) \mathrm{d}x \text{ assuming }n :: \text{nonnegint}, \alpha > 0$" immediately into $\Gamma(2n+1, \alpha) \alpha^{-2n-1}$; I'm sure Mathematica does the same. (Still a good idea to verify it numerically, of course, which these guys can typically also do.) Dec 15, 2011 at 4:06