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

问题描述

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

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?

推荐答案

仅将其分为实部和虚部是怎么回事? scipy.integrate.quad要求集成函数为其使用的算法返回浮点数(即实数).

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

例如

>>> 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).

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).

出于记录的目的，如果不是所有人都清楚知道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).

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.

注意:Scipy.integrate将不会直接处理复杂的集成.为什么?它在FORTRAN QUADPACK 库(尤其是 qagse.f ，它明确要求函数/变量必须为实数，然后才能进行基于每个子间隔内21点高斯-克朗德德正交的全局自适应正交运算，并通过Peter Wynn的epsilon算法进行加速."因此，除非您想尝试修改底层的FORTRAN以使其能够处理复数，然后将其编译到新的库中，否则您将无法使其正常工作.

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.

如果您真的想在一个积分中对复数进行Gauss-Kronrod方法，请查看 Wikipedias页面，并按照以下步骤直接实施(使用15点，7点规则).注意，我记住函数会重复对公共变量的公共调用(假设函数调用很慢，就好像函数非常复杂).也只做7点和15点规则，因为我自己不想计算节点/权重，而这些是维基百科上列出的，但是对于测试用例却得到了合理的误差(〜1e-14)

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]

测试用例:

>>> 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才能达到某种幂或类似的强度).

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).

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

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屋！