# Integration of an indefinite integral: matlab precision problem

The integral I need to evaluate is:

$$\int_x^{\infty} \frac{t^n}{e^{t} -1} dt$$

After some research I found a paper saying,

The numerical values of the two integrals [...] are easily calculated either directly or through series expansion.

It's easy to evaluate, right? So I try plugging this into my 32bit 7.7.0.471 (R2008b) MATLAB program and it cannot integrate it. Fortunately they actually cited this and I was lead to a mathematics table which stated, $$\int_x^{\infty} \frac{t^n dt}{e^t -1} = \sum_{k=1}^{\infty} \exp^{-kx} ( \frac{x^n}{k} + \frac{n x^{n-1}}{k^2} + \frac{ (n)(n-1) x^{n-2}}{k^3} + \cdots +\frac{n!}{k^{n+1}} )$$

I'm using $n=2$ so my actual expression is $$\int_x^{\infty} \frac{t^2 dt}{e^t -1} = \sum_{k=1}^{\infty} \exp(-kx) ( \frac{x^2}{k} + \frac{2 x^{1}}{k^2} + \frac{ (2)(1) x^{0}}{k^3})$$

The problem is, however, that my $x$ is very large for the environment, for example $10^{19}$, and the precision for $\exp$ causes it to be evaluated as 0. If I try to scale this on a $\log$ basis then I can't seem to get it converge even after several thousand summations. My software is also too old to be compatible with newer contributions which would solve my problem by allowing me to arbitrarily increase the precision.

This paper is from the mid 90's, why did they say it is so easy? Is there a better way, mathematically, to attack this that a novice like myself can understand and implement? How does one even find this series? It seems crazy.

I did find a paper that seems to be doing what I want, although I am not sure how to implement it.

• When $x$ is very large, $e^{-kx}$ is very very small. Are you sure it evaluates to 1? For $x=10^{19}$, the first term of the sum is on the order of $e^{-10^{19}}\approx 1/10^{4\times 10^{18}}$, which would underflow a double-precision floating point variable, and further terms are even smaller. It should converge very quickly for small $n$ and large $x$. Did you make a mistake? May 15 '13 at 20:53
• Sorry, meant 0. May 15 '13 at 20:57
• So when you say it doesn't converge, does it actually just converge to 0? A number that small would not be representable in double-precision floating point arithmetic, so it would just be 0. For $n=2$ and $x=10^{19}$ you could even just take the first term in the sum and it would be a good approximation of the actual value. The second term would be about $e^{-10^{19}}$ times the first term, and so negligible. May 15 '13 at 21:02
• You're missing the point of my problem: it's $\exp(- 10^{19} ) * ( (10^{19})^2/k + 2 * (10^{19}/k^2 + 1/(k^3) )$, that expression will not give me the right answer here with this limited precision. May 15 '13 at 21:05
• I expect the answers to be in the range of $10^{0}$ to $10^{-70}$ for my full range of $x$. May 15 '13 at 21:11

For $n>0, x=0$, this evaluates to $\Gamma(n+1)\ \mathrm{Li}_{n+1}(1)$, and for $n=2$ this is $2\ \zeta (3)$.
For the first half, maybe this helps regarding more references. Here, is an expression for the other half and even a closed (well, in terms of polylogs) for for $n=2$.
• I cannot assume $x=0$ or my work would be done. My calculations have an actual $x$ dependency and apparently it has been done before. (I just don't know how.) I tried using the polylogarithm, however, I cannot get the precision high enough even when using the VPA function in MATLAB May 15 '13 at 20:21
• @sciencenewbie: "I cannot assume $x=0$"... you don't have to assume anything. The point I make is that the result $2\ \zeta(3)$ for $\int_0^\infty$ enables you to switch from computing $\int_x^\infty$ to computing $\int_0^x$, for which we also know the series expansion. May 15 '13 at 23:20