python - Use scipy.integrate.quad to integrate complex numbers -

i'm using right scipy.integrate.quad integrate real integrands. situation appeared need integrate complex integrand. quad seems not able it, other scipy.integrate routines, ask: there way integrate complex integrand using scipy.integrate, without having separate integral in real , imaginary parts?

what's wrong separating out real , imaginary parts? scipy.integrate.quad requires integrated function return floats (aka real numbers) algorithm uses.

import scipy scipy.integrate import quad  def complex_quadrature(func, a, b, **kwargs):     def real_func(x):         return scipy.real(func(x))     def imag_func(x):         return scipy.imag(func(x))     real_integral = quad(real_func, a, b, **kwargs)     imag_integral = quad(imag_func, a, b, **kwargs)     return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:]) 

e.g.,

>>> complex_quadrature(lambda x: (scipy.exp(1j*x)), 0,scipy.pi/2) ((0.99999999999999989+0.99999999999999989j),  (1.1102230246251564e-14,),  (1.1102230246251564e-14,)) 

which expect rounding error - integral of exp(i x) 0, pi/2 (1/i)(e^i pi/2 - e^0) = -i(i - 1) = 1 + ~ (0.99999999999999989+0.99999999999999989j).

and record in case isn't 100% clear everyone, integration linear functional, meaning ∫ { f(x) + k g(x) } dx = ∫ f(x) dx + k ∫ g(x) dx (where k constant respect x). or our specific case ∫ z(x) dx = ∫ re z(x) dx + ∫ im z(x) dx z(x) = re z(x) + im z(x).

if trying integration on path in complex plane (other along real axis) or region in complex plane, you'll need more sophisticated algorithm.

note: scipy.integrate not directly handle complex integration. why? heavy lifting in fortran quadpack library, in qagse.f explicitly requires functions/variables real before doing "global adaptive quadrature based on 21-point gauss–kronrod quadrature within each subinterval, acceleration peter wynn's epsilon algorithm." unless want try , modify underlying fortran handle complex numbers, compile new library, aren't going work.

if want gauss-kronrod method complex numbers in 1 integration, @ wikipedias page , implement directly done below (using 15-pt, 7-pt rule). note, memoize'd function repeat common calls common variables (assuming function calls slow if function complex). did 7-pt , 15-pt rule, since didn't feel calculating nodes/weights myself , ones listed on wikipedia, getting reasonable errors test cases (~1e-14)

import scipy scipy import array  def quad_routine(func, a, b, x_list, w_list):     c_1 = (b-a)/2.0     c_2 = (b+a)/2.0     eval_points = map(lambda x: c_1*x+c_2, x_list)     func_evals = map(func, eval_points)     return c_1 * sum(array(func_evals) * array(w_list))  def quad_gauss_7(func, a, b):     x_gauss = [-0.949107912342759, -0.741531185599394, -0.405845151377397, 0, 0.405845151377397, 0.741531185599394, 0.949107912342759]     w_gauss = array([0.129484966168870, 0.279705391489277, 0.381830050505119, 0.417959183673469, 0.381830050505119, 0.279705391489277,0.129484966168870])     return quad_routine(func,a,b,x_gauss, w_gauss)  def quad_kronrod_15(func, a, b):     x_kr = [-0.991455371120813,-0.949107912342759, -0.864864423359769, -0.741531185599394, -0.586087235467691,-0.405845151377397, -0.207784955007898, 0.0, 0.207784955007898,0.405845151377397, 0.586087235467691, 0.741531185599394, 0.864864423359769, 0.949107912342759, 0.991455371120813]     w_kr = [0.022935322010529, 0.063092092629979, 0.104790010322250, 0.140653259715525, 0.169004726639267, 0.190350578064785, 0.204432940075298, 0.209482141084728, 0.204432940075298, 0.190350578064785, 0.169004726639267, 0.140653259715525,  0.104790010322250, 0.063092092629979, 0.022935322010529]     return quad_routine(func,a,b,x_kr, w_kr)  class memoize(object):     def __init__(self, func):         self.func = func         self.eval_points = {}     def __call__(self, *args):         if args not in self.eval_points:             self.eval_points[args] = self.func(*args)         return self.eval_points[args]  def quad(func,a,b):     ''' output 15 point estimate; , estimated error '''     func = memoize(func) #  memoize function skip repeated function calls.     g7 = quad_gauss_7(func,a,b)     k15 = quad_kronrod_15(func,a,b)     # don't have faith in error estimate taken wikipedia     # without incorporating how should scale changing limits     return [k15, (200*scipy.absolute(g7-k15))**1.5] 

test case:

>>> quad(lambda x: scipy.exp(1j*x), 0,scipy.pi/2.0) [(0.99999999999999711+0.99999999999999689j), 9.6120083407040365e-19] 

i don't trust error estimate -- took wiki recommended error estimate when integrating [-1 1] , values don't seem reasonable me. e.g., error above compared truth ~5e-15 not ~1e-19. i'm sure if consulted num recipes, more accurate estimate. (probably have multiple (a-b)/2 power or similar).

recall, python version less accurate calling scipy's quadpack based integration twice. (you improve upon if desired).