30

我现在正在使用 scipy.integrate.quad 来成功集成一些真正的被积函数。现在出现了一种情况,我需要整合一个复杂的被积函数。与其他 scipy.integrate 例程一样,quad 似乎无法做到这一点,所以我问:有没有办法使用 scipy.integrate 积分一个复杂的被积函数,而不必分离实部和虚部的积分?

4

2 回答 2

57

把它分成实部和虚部有什么问题? 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)。

并且记录一下,以防每个人都不是 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 as z(x) = Re z(x) + i Im z(x)。

如果您尝试对复平面中的路径(而不是沿实轴)或复平面中的区域进行积分,则需要更复杂的算法。

注意: Scipy.integrate 不会直接处理复杂的集成。为什么?它在 FORTRAN QUADPACK库中完成了繁重的工作,特别是在qagse.f中,它明确要求函数/变量在执行其“基于每个子区间内的 21 点 Gauss-Kronrod 求积的全局自适应求积,由 Peter 加速”之前为实数永利的 epsilon 算法。” 因此,除非您想尝试修改底层 FORTRAN 以使其处理复数,将其编译到新库中,否则您将无法使其正常工作。

如果您真的想在一个集成中使用复数执行 Gauss-Kronrod 方法,请查看wikipedias 页面并按如下方式直接实现(使用 15-pt、7-pt 规则)。注意,我记住了函数来重复对公共变量的常用调用(假设函数调用很慢,好像函数非常复杂)。也只做了 7-pt 和 15-pt 规则,因为我不想自己计算节点/权重,那些是维基百科上列出的,但在测试用例中得到了合理的错误(~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 的集成两次准确。(如果需要,您可以对其进行改进)。

于 2011-05-11T14:40:58.117 回答
9

我意识到我迟到了,但也许quadpy(我的一个项目)可以提供帮助。这

import quadpy
import numpy

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

正确给出

(0.8414709848078964+0.4596976941318605j)
于 2017-03-17T20:27:51.367 回答