使用scipy.integrate.quad对复数进行积分

Question

我现在正在使用scipy.integrate.quad这个工具，成功地对一些真实的函数进行积分。现在出现了一种情况，我需要对一个复杂的函数进行积分。quad似乎无法处理这个问题，其他的scipy.integrate工具也是如此。所以我想问：有没有办法使用scipy.integrate对复杂的函数进行积分，而不需要把积分分成实部和虚部来处理呢？

Answer 1

我知道我来得有点晚，但也许我的一个项目 quadpy 能帮上忙。这个

import quadpy
import numpy

val, err = quadpy.quad(lambda x: numpy.exp(1j * x), 0, 1)
print(val)

可以正确地给出

(0.8414709848078964+0.4596976941318605j)

Answer 2

把它分成实部和虚部有什么问题呢？scipy.integrate.quad这个函数要求被积分的函数返回的是浮点数（也就是实数），这是它使用的算法所需要的。

import scipy
from 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:])

比如，

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

你会期待到的舍入误差 - 从0到π/2的exp(i x)的积分是(1/i)(e^i π/2 - e^0) = -i(i - 1) = 1 + i ~ (0.99999999999999989+0.99999999999999989j)。

为了确保大家都明白，积分是一个线性运算，这意味着∫ { f(x) + k g(x) } dx = ∫ f(x) dx + k ∫ g(x) dx（这里的k是一个常数，与x无关）。在我们的具体情况下，∫ z(x) dx = ∫ Re z(x) dx + i ∫ Im z(x) dx，因为z(x) = Re z(x) + i Im z(x)。

如果你想在复平面上（除了沿着实轴）进行路径或区域的积分，你需要一个更复杂的算法。

注意：Scipy.integrate不会直接处理复数积分。为什么呢？因为它在FORTRAN的QUADPACK库中完成了繁重的工作，特别是在qagse.f中，这个函数明确要求在进行“基于21点高斯-克朗罗德积分的全局自适应求积”之前，函数/变量必须是实数。所以，除非你想修改底层的FORTRAN代码来处理复数，并将其编译成一个新的库，否则你是无法让它正常工作的。

如果你真的想用复数进行一次积分的高斯-克朗罗德方法，可以查看维基百科页面，并直接实现，如下所示（使用15点和7点规则）。注意，我对函数进行了记忆化处理，以便重复调用常用变量（假设函数调用很慢，因为函数可能非常复杂）。另外，我只用了7点和15点规则，因为我不想自己计算节点/权重，而这些正是维基百科上列出的，但对于测试用例得到了合理的误差（大约是~1e-14）。

import scipy
from 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 is the 15 point estimate; and the estimated error '''
    func = Memoize(func) #  Memoize function to skip repeated function calls.
    g7 = quad_gauss_7(func,a,b)
    k15 = quad_kronrod_15(func,a,b)
    # I don't have much faith in this error estimate taken from wikipedia
    # without incorporating how it should scale with changing limits
    return [k15, (200*scipy.absolute(g7-k15))**1.5]

测试用例：

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

我对误差估计不太信任——我从维基百科上拿了一些推荐的误差估计，用于从[-1到1]的积分，但这些值对我来说似乎不太合理。例如，上面的误差与真实值相比大约是~5e-15，而不是~1e-19。如果有人查阅数值计算的书籍，可能会得到更准确的估计。（可能需要乘以(a-b)/2的某个幂或类似的东西）。

记住，Python版本的准确性不如直接调用scipy的QUADPACK积分两次。（如果需要，你可以对此进行改进）。