2

I'm calculating the upper incomplete gamma function using Python's mpmath.gammainc:

import numpy as np
from mpmath import gammainc    

z = 0 # define power of t as t^(-1)
a = 0.5+0.4j # integral lower limit
b = np.inf # integral upper limit

myfun = np.array([gammainc(z,a,b,regularized=False)], dtype=complex)

It is a one-dimensional integral as defined in the mpmath docs. I want to compare this result myfun using scipy's quad function:

myfun2 = scipy.integrate.quad(exp(-t)/t, a, inf)[0]

However I don't think quad accepts complex arguments for the upper/lower bounds of integration. I don't know if it's possible to separate the problem into real/imaginary parts. Any ideas?

Medulla Oblongata
  • 3,771
  • 8
  • 36
  • 75
  • Could you gives more details about the mathematical definition of the integration with complex bounds? I don't think the formula `exp(-t)/t` is actually used for the complex case, I am trying to make sense of the [Incomplete_gamma_function](https://en.wikipedia.org/wiki/Incomplete_gamma_function). `quad` will not work with complex argument: an integration path have to be defined (or an area) on the complex plane... – xdze2 Aug 08 '18 at 13:37

1 Answers1

4

The integral should be taken over the horizontal half-line going from a to the right. (This will fail if a is a negative real number, but that is a branch cut territory anyway). This halfline is parameterized by a+t where t is real and goes from 0 to infinity. So, integrate exp(-(a+t))/(a+t) from 0 to infinity.

Also, quad from SciPy requires real-valued functions, so separate it into a real and imaginary part. And don't forget that the function must be passed as a callable, like lambda t: np.exp(-t), not just exp(-t)

from scipy import integrate
myfun2_re = integrate.quad(lambda t: np.real(np.exp(-(a+t))/(a+t)), 0, np.inf)[0]
myfun2_im = integrate.quad(lambda t: np.imag(np.exp(-(a+t))/(a+t)), 0, np.inf)[0]
myfun2 = myfun2_re + 1j*myfun2_im
print(myfun2)

This prints

(0.3411120086192922-0.36240971724285814j)

compared to myfun, 

array([ 0.34111201-0.36240972j])

By the way, there is no need to wrap myfun into an array if you just want to convert a single number: complex(gammainc(...)) would do.

  • "The integral should be taken over the horizontal half-line going from a to the right." is it somehow general or specific to the incomplete gamma function? – xdze2 Aug 09 '18 at 11:44
  • Specific to this function, which is defined as an integral from a to positive infinity when a is real. That sets the general direction where the path of integration should go. There is a lot of freedom in choosing its shape because such an integral is independent of path as long as we don't go around a singularity (there is one at 0) or encounter large values at infinity (there are none in the right half plane). Horizontal line is the simplest path available. Bottom line: one needs to know some complex analysis to work with special functions. –  Aug 09 '18 at 12:38