使用scipy.integrate.quad积分复数

Question

我现在正在使用scipy.integrate.quad成功地集成一些实际的被积物。现在出现了一种情况，我需要集成一个复杂的被积物。和其他scipy.integrate例程一样，quad似乎无法做到这一点，所以我问：有没有办法使用scipy.integrate集成复杂的被积函数，而不必将实部和虚部分开？

Answer 1

仅仅将其分为实部和虚构部分有什么问题？ 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,))

这就是您期望的舍入误差-exp（ix）的整数从0开始，pi / 2为（1 / i）（e ^ i pi / 2-e ^ 0）= -i（i-1）= 1
+我〜（0.99999999999999989 + 0.99999999999999989j）。

记录不是所有人都清楚地知道，积分是线性函数，这意味着∫{f（x）+ kg（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点高斯-
科隆德积分的全局自适应正交运算，并由Peter进行加速”。永利的epsilon算法。”
因此，除非您想尝试修改基础的FORTRAN以使其能够处理复数，然后将其编译到新的库中，否则您将无法使其工作。

如果您真的想在一次积分中使用复数进行Gauss-
Kronrod方法，请查看Wikipedias页面并按照下面的方法直接实现（使用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]集成时，我从Wiki上获取了一些建议的错误估计，而这些值对我来说似乎并不合理。例如，与真值相比，上述误差是〜5e-15而不是〜1e-19。我敢肯定，如果有人咨询了num食谱，那么您可以获得更准确的估计。（可能必须乘以(a-b)/2某种力量或类似的东西）。

回想一下，python版本的准确性不如仅仅两次调用scipy的基于QUADPACK的集成。（如果需要，您可以对其进行改进）。