使用 scipy.integrate.quad 对复数进行积分 [英] Use scipy.integrate.quad to integrate complex numbers

问题描述

我现在正在使用 scipy.integrate.quad 来成功集成一些真正的被积函数.现在出现了一种情况，我需要对一个复杂的被积函数进行积分.quad 似乎无法做到，因为其他 scipy.integrate 例程，所以我问:有没有办法使用 scipy.integrate 对复杂的被积函数进行积分，而不必将实部和虚部中的积分分开?

解决方案

将其分成实部和虚部有什么问题?scipy.integrate.quad 需要集成函数为其使用的算法返回浮点数(又名实数).

导入scipy从 scipy.integrate 导入四边形def complex_quadrature(func, a, b, **kwargs):def real_func(x):返回 scipy.real(func(x))def imag_func(x):返回 scipy.imag(func(x))real_integral = quad(real_func, a, b, **kwargs)imag_integral = quad(imag_func, a, b, **kwargs)返回 (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 + i ~ (0.99999999999999989+0.99999999999999989j).

记录，以防每个人都不是 100% 清楚，积分是一个线性函数，这意味着 ∫ { 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 点 Gauss-Kronrod 正交的全局自适应正交，并由 Peter Wynn 的 epsilon 算法加速"之前为实数."因此，除非您想尝试修改底层的 FORTRAN 以使其处理复数，否则将其编译到一个新的库中，否则您将无法使其工作.

如果您真的想在一次积分中使用复数进行 Gauss-Kronrod 方法，请查看 wikipedias page 并按照下面的方法直接实现(使用 15-pt, 7-pt 规则).注意，我记住了函数来重复对公共变量的公共调用(假设函数调用很慢，好像函数非常复杂).也只做了 7-pt 和 15-pt 规则，因为我不想自己计算节点/权重，那些是维基百科上列出的，但是测试用例会出现合理的错误(~1e-14)

导入scipy从 scipy 导入数组def quad_routine(func, a, b, x_list, w_list):c_1 = (b-a)/2.0c_2 = (b+a)/2.0eval_points = map(lambda x: c_1*x+c_2, x_list)func_evals = 地图(函数，eval_points)返回 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.7415907397, 0.741597397, 0.745597394949494w_gauss = 数组([0.129484966168870, 0.279705391489277, 0.381830050505119, 0.417959183673469, 0.381830059590.381830059590.780.597879879879879879800返回 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]返回 quad_routine(func,a,b,x_kr, w_kr)类记忆(对象):def __init__(self, func):self.func = funcself.eval_points = {}def __call__(self, *args):如果参数不在 self.eval_points 中:self.eval_points[args] = self.func(*args)返回 self.eval_points[args]def quad(func,a,b):''' 输出是 15 点估计值；和估计误差'''func = Memoize(func) # 记忆函数以跳过重复的函数调用.g7 = quad_gauss_7(func,a,b)k15 = quad_kronrod_15(func,a,b)# 我不太相信来自维基百科的这个错误估计# 没有考虑它应该如何随着不断变化的限制而扩展返回 [k15, (200*scipy.absolute(g7-k15))**1.5]

测试用例:

<预><代码>>>>四边形(lambda x:scipy.exp(1j*x)，0，scipy.pi/2.0)[(0.99999999999999711+0.9999999999999689j), 9.6120083407040365e-19]

我不相信误差估计——当从 [-1 到 1] 集成时，我从 wiki 中获取了一些推荐的误差估计，并且这些值对我来说似乎不合理.例如，上面与真值相比的错误是 ~5e-15 而不是 ~1e-19.我敢肯定，如果有人咨询过 num 食谱，您可以获得更准确的估计.(可能必须通过 (a-b)/2 乘以某种力量或类似的东西).

回想一下，python 版本不如仅调用 scipy 的基于 QUADPACK 的集成两次准确.(如果需要，您可以改进它).

I'm using right now the scipy.integrate.quad to successfully integrate some real integrands. Now a situation appeared that I need to integrate a complex integrand. quad seems not be able to do it, as the other scipy.integrate routines, so I ask: is there any way to integrate a complex integrand using scipy.integrate, without having to separate the integral in the real and the imaginary parts?

解决方案

What's wrong with just separating it out into real and imaginary parts? scipy.integrate.quad requires the integrated function return floats (aka real numbers) for the algorithm it uses.

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:])

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 is what you expect to rounding error - integral of exp(i x) from 0, pi/2 is (1/i)(e^i pi/2 - e^0) = -i(i - 1) = 1 + i ~ (0.99999999999999989+0.99999999999999989j).

And for the record in case it isn't 100% clear to everyone, integration is a linear functional, meaning that ∫ { f(x) + k g(x) } dx = ∫ f(x) dx + k ∫ g(x) dx (where k is a constant with respect to x). Or for our specific case ∫ z(x) dx = ∫ Re z(x) dx + i ∫ Im z(x) dx as z(x) = Re z(x) + i Im z(x).

If you are trying to do a integration over a path in the complex plane (other than along the real axis) or region in the complex plane, you'll need a more sophisticated algorithm.

Note: Scipy.integrate will not directly handle complex integration. Why? It does the heavy lifting in the FORTRAN QUADPACK library, specifically in qagse.f which explicitly requires the functions/variables to be real before doing its "global adaptive quadrature based on 21-point Gauss–Kronrod quadrature within each subinterval, with acceleration by Peter Wynn's epsilon algorithm." So unless you want to try and modify the underlying FORTRAN to get it to handle complex numbers, compile it into a new library, you aren't going to get it to work.

If you really want to do the Gauss-Kronrod method with complex numbers in exactly one integration, look at wikipedias page and implement directly as done below (using 15-pt, 7-pt rule). Note, I memoize'd function to repeat common calls to the common variables (assuming function calls are slow as if the function is very complex). Also only did 7-pt and 15-pt rule, since I didn't feel like calculating the nodes/weights myself and those were the ones listed on wikipedia, but getting reasonable errors for test cases (~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]

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 the error estimate -- I took something from wiki for recommended error estimate when integrating from [-1 to 1] and the values don't seem reasonable to me. E.g., the error above compared with truth is ~5e-15 not ~1e-19. I'm sure if someone consulted num recipes, you could get a more accurate estimate. (Probably have to multiple by (a-b)/2 to some power or something similar).

Recall, the python version is less accurate than just calling scipy's QUADPACK based integration twice. (You could improve upon it if desired).

这篇关于使用 scipy.integrate.quad 对复数进行积分的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！