Python中的复杂集成答案

【问题标题】：Complex integration in PythonPython中的复杂集成
【发布时间】：2015-12-08 22:07:55
【问题描述】：

有一个 MATLAB 函数 quadgk 可以计算复积分，或者至少可以计算具有极点和奇点的函数。在 Python 中，有一个通用的scipy.integrate.quad，可以方便地沿实轴进行集成。是否有与 MATLAB 的quadgk 等效的 Python？我只能在另一个 SO 问题中找到dr jimbob's code，这似乎不再适用于 Python 3.4。

【问题讨论】：

我建议你强调你需要一个复杂的contour。对于真实轮廓，您仍然可以使用scipy.integrate.quad 分别计算实部和虚部。如果所有其他方法都失败了，您仍然可以手动完成：定义一个前端函数，为两个复杂限制定义线性路径 a1+i*a2 和 b1+i*b2，它们之间有一个实参数，然后使用 quad产生的实部和虚部......我承认它远非漂亮，我希望有一个更简单、现成的解决方案:)

标签： python matlab function numpy scipy

【解决方案1】：

我不认为 SciPy 确实提供了与 MATLAB 的 quadgk 等效的功能，但您在 this question 中链接的代码只需稍作改动即可在 Python 3 中运行：

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 = list(map(func, eval_points))    # Python 3: make a list here
    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:                     # Python 3: no need to inherit from 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]

例如，

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

【讨论】：

我认为问题出在复杂的轮廓上，比如从-i 整合到i。该解决方案是否也适用于这种情况？
如果我的集成极限达到无穷大，比如print(quadgk(lambda x: scipy.exp(1j*x), -np.inf,0))，怎么办？然后我得到[(nan+nan*j), nan]。