# 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? – Faheem Mitha Dec 13 '11 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. – drjrm3 Dec 13 '11 at 17:16
• You might be right. – Faheem Mitha Dec 13 '11 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... – David Z Dec 13 '11 at 19:03
• @Faheem: "Monte Carlo is an extremely bad method; it should be used only when all alternative methods are worse." - Alan Sokal – J. M. Dec 14 '11 at 7:31

## 2 Answers

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}$.

• absolutely no idea how i overlooked this. thank you. – drjrm3 Dec 13 '11 at 18:09
• Clever substitutions and integration by parts should always be one of the first things you do with unruly integrals. – J. M. Dec 14 '11 at 4:30
• 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.) – Erik P. Dec 15 '11 at 4:06
• Actually Mathematica chooses to represent the answer as ExpIntegralE[-2 n, a r]. If you run FunctionExpand on it, then it gives the same answer as Maple. – Searke Jan 3 '12 at 21:36

Have a look at QUADPACK. It has routines for integration over (semi-)infinite domains.