# Numerical integration of a single integral

I am trying to evaluate this integral directly using numerical integration functions in Mathematica and Python.

$$\int_0^\infty {y^{-(n+1)}\prod_{i=1}^n \gamma(a_i+1,b_i y) }dy$$

where $\gamma(a,x)=\int_0^x dt\ t^{a-1}e^{-t}$ is the lower incomplete Gamma function. In my actual problem $n$ can take values from 10 to 100, where $a_i$ can take values from 100 to 1000.

Here is the Python code

import numpy as np
import scipy.integrate as integrate
import scipy.special as special

def integrand(a,b,y):
l=len(a)
x=pow(y,-(l+1))
for ai,bi in zip(a,b):
# x*=special.gammainc(ai+1,bi*y)
x*=1/bi*special.gammainc(ai+1,bi*y)
return x

ns=[2,4,5,6]

print integrate.quad(lambda x: integrand([n for n in ns],[0.1,0.8,0.3,0.4],x),0,np.inf)
print integrate.quad(lambda x: integrand([10*n for n in ns],[0.1,0.8,0.3,0.4],x),0,np.inf)
print integrate.quad(lambda x: integrand([100*n for n in ns],[0.1,0.8,0.3,0.4],x),0,np.inf)


the corresponding output is

$(8.231566916434924e-05, 2.7894916135933814e-10)$
$(5.884546215247849e-09, 1.1634128468162877e-08)$
$(0.0, 0.0)$

Which is disappointing considering the modest $n$=4. In the second output $a_i$ ~ $60$ and the error is greater for one order of magnitude than the resulting integral. Don't need to comment the third trial case with $a_i$ ~ $600$, with zero error!

Of course that is not magic, I am using the full integration interval $[0,\infty]$, and ask for numerical tricks like reducing integration interval and any other treatment to get reasonable results.

I made some easy progress by multiplying by a constant ($\prod_r 1/b_i$) and the numbers becomes more stable within all the range of my application.

For example, by plotting the integrand I see an unimodal curve which span the ranges $[5,60]$, $[50,600]$ and $[500,6000]$ respectively, vanishing outside. If I truncate the integration intervals respectively the results are quite acceptable and agree with Wolfram Mathematica results. I tried other larger sets, and find the intervals by trial/error. So I will be happy if I can find the proper intervals.

How can I find the interval for any sets of $a_i$ and $b_i$ and a multiplicative constant that made the numbers stable ?

• You can always do a change of variables, like $u = \frac{1}{1+y}$, to change the interval of the integral to something finite. – spektr Feb 16 '16 at 16:54
• Try calling quad with epsabs=0.0, maybe the tiny magnitude is what's confusing quad. – Kirill Feb 16 '16 at 17:12
• I want to point out just how badly behaved the integrand is: in the second example, the integrand goes from $(0,0)$ to $(205.6,6.2\times10^{200})$ then goes to zero as $y^{-5}$. Perhaps you can look at Laplace's method or similar to get an approximate value in another way? – Kirill Feb 16 '16 at 20:36

Improper integrals like this are often no match for double exponential quadrature. I can compute your integrals with mpmath. For larger n, it appears that a higher precision needs to be used.

from mpmath import mp, gammainc, quad, power, inf, nstr

def integrand(a,b,y):
l=len(a)
x=power(y,-(l+1))
for ai,bi in zip(a,b):
x*=gammainc(ai+1,0,bi*y,regularized=True) / bi
return x

ns=[2,4,5,6]

mp.pretty = True
mp.dps = 30
print quad(lambda x: integrand([n for n in ns],[0.1,0.8,0.3,0.4],x),[0,inf],error=True)
print quad(lambda x: integrand([10*n for n in ns],[0.1,0.8,0.3,0.4],x),[0,inf],error=True)
mp.dps = 100
print nstr(quad(lambda x: integrand([100*n for n in ns],[0.1,0.8,0.3,0.4],x),[0,inf],error=True),20)


The output is:

(0.0000823156691641082138376342855035, 1.0e-41)
(0.0000000156728839596991659411718340463, 1.0e-28)
(1.6721149546712801741e-12, 1.0e-32)


Splitting the interval at one point also does the trick:

>>> mp.dps = 30
>>> print quad(lambda x: integrand([100*n for n in ns],[0.1,0.8,0.3,0.4],x),[0,5000,inf],error=True)
(1.6721149546712801740899041739e-12, 1.0e-32)


If you don't know which parameters to use, you could make it automatic: check the error estimate from quad() and, if it's too large, either double the precision or subdivide the interval in two and recurse.

Quadrature on improper integrals can always be problematic. If you're finding large errors, I'd suggest one of two things:

• As choward suggested, use a change of variable trick to change the integral to have finite bounds, then use quadrature as usual.

• Use a quadrature method that uses basis functions appropriate for your problem. In this case, the Laguerre Polynomials (or Rational Chebyshev on semi-infinite intervals) work (ref Boyd: Chebyshev and Fourier Spectral Methods.) Laguerre Quadrature should work better. There's even a numpy module for it.

• But the weight function is $y^{-1-n}$, not $e^{-y}$ as for Gauss-Laguerre? – Kirill Feb 16 '16 at 20:22
• You can rewrite it as shown in the wiki link there - you'd write $y^{-1-n} e^y e^{-y}$. But good point on how badly behaved the integrand is, I didn't look at it that close. That's probably always going to be horrible. – Aurelius Feb 16 '16 at 20:55
• I think the problem with writing $e^{-x}(e^xf(x))$ is that the errors depend on high derivatives of the integrand $g(x) = e^xf(x)$, which will be unusually large due to the exponential present. For example, try $\frac{\pi}{2} = \int_0^\infty e^{-x} \frac{e^x}{1+x^2}\,dx$ using that suggestion, and compare its (very slow) convergence with $\int_0^\infty \frac{e^{-x}}{1+x^2}$. I don't think wikipedia is right about this at all. In fact, I remember a question I answered on MSE – Kirill Feb 16 '16 at 21:20
• Ah, fair points. – Aurelius Feb 16 '16 at 21:48