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?
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.