Solution of an integral with strange imprecision of gamma functions

Trying to solve the following integral, with $n,m \in \mathbb{Z^+}$, $\alpha>1$, $0 < \epsilon < 1$, and $\Gamma(.)$ and $\Gamma(.,.)$ the gamma and incomplete gamma functions, respectively: $$I(\alpha,m,n) = \frac{(\alpha (n-1))^n }{\Gamma (n)-\Gamma \left(n,\frac{(n-1) \alpha }{\epsilon +1}\right)}\int_1^{1/\epsilon}e^{\frac{(y-1) (\alpha -\alpha n)}{y}}(y-1)^{n-1} y^{m-n-1} dy$$ My approach has been as follows, leading to expressing the integral with gamma functions, which causes precision problems. Since $(y-1)^{n-1}= \sum _{i=0}^{n-1} (-1)^i \binom{n-1}{i} y^{n-1-i}$, and since $$\int_1^{1/\epsilon } y^{m-2-i} e^{\frac{(y-1) }{y}(\alpha -\alpha n)} \, dy = \frac{e^{\alpha (1- n)} \left(\frac{1}{\alpha -\alpha n}\right)^{i-m} (\Gamma (i-m+1,\alpha -n \alpha )-\Gamma (i-m+1,-(n-1) \alpha \epsilon ))}{\alpha (n-1)},$$ we get as a possible solution a summation showing ratios of differences of various gamma functions:

$I(\alpha, m,n)=e^{\alpha (1- n)} (\alpha (n-1))^{n-1}$

$$\sum _{i=0}^{n-1} (-1)^i \binom{n-1}{i}\frac{ \left(\frac{1}{\alpha -\alpha n}\right)^{i-m} (\Gamma (i-m+1,\alpha -n \alpha )-\Gamma (i-m+1,-(n-1) \alpha \epsilon ))}{\Gamma (n)-\Gamma \left(n,\frac{(n-1) \alpha }{\epsilon +1}\right)}$$

Now the hitch. I am computing (for control) the expression of $I(\alpha,m, n)$ integral using high precision numerical integration. With n=10, the numerical integration matches the gamma function. Beyond, it, the gamma functions (using Wolfram's Mathematica) produce different results (and from familiarity with the problem, I know that the gamma is wrong compared to the numerical). Yet I need $n$ around $10^4$. Further, other special functions such as the exponential integral give even worse results.

The results are in the attached picture. With $\epsilon=\frac{1}{100}$, $I(\frac{3}{2},1,n)$:

The questions are:

1) Is there a possible closed form for $I$not entailing special functions?

2) Is there a better way to solve the integral by avoiding the summation?

With gratitude.

• When calculating the values in the table, what value of $\epsilon$ did you use ? – Johannes Trost Sep 18 '15 at 15:11
• .01. I will add. Thanks – Nero Sep 18 '15 at 15:26
• Have you checked the difference between the asymptotic expansions of the two Gamma functions in the denominator for large n? Your results suggest that the denominator rapidly becomes very small as you increase n. – Biswajit Banerjee Sep 18 '15 at 22:55
• How many bits/digits of precision is your numerical integration using? I'm not familiar with Mathematica but trying to reproduce your numbers in sympy. – Ahmed Fasih Nov 6 '15 at 14:54
• Have you tried the usual strategies, like trying to solve (or have MMA solve) simpler versions first: $\epsilon -> 0$, $\alpha = 1$, $m=n$, etc? – Mitch Nov 6 '15 at 15:17

I don't think closed form is possible, but here is (possibly) a more controlled series formula for the function: $$I(\alpha, m, n) = \sum_{j=0}^{\infty} {}^{m+j-1}C_{j} \,\,(\alpha(n-1))^{-j} \left[\frac{\Gamma (j+n) - \Gamma (j+n, \alpha(n-1)(1-\epsilon))}{\Gamma(n) - \Gamma(n, (n-1)\alpha/(1+\epsilon))}\right]$$
I got this using the variable substitution $z=(y-1)/y$ and carrying out a Taylor expansion in powers of $z$.
• Had the same result. This probably scales a lot better for high $n$. All terms are strictly positive. I wonder if it is possible to find good upper and lower bounds. – StijnDeVuyst Nov 7 '15 at 0:32