5

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 ?

jss
  • 129
  • 4
reg
  • 51
  • 2
  • 1
    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
  • 3
    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

2 Answers2

4

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.

Fredrik Johansson
  • 375
  • 2
  • 7
2

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.

Aurelius
  • 2,365
  • 1
  • 15
  • 19
  • 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
  • 1
    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