Python+Scipy+Integration：处理带有尖峰的函数中的精度误差

Question

我正在尝试使用 scipy.integrate.quad 在非常大的范围（0..10,000）内集成函数。该函数在其大部分范围内为零，但在非常小的范围内有一个尖峰（例如 1,602..1,618）。

积分时，我希望输出为正，但我猜想四边形的猜测算法不知何故变得混乱并输出零。我想知道的是，有没有办法克服这个问题（例如通过使用不同的算法、其他一些参数等）？我通常不知道峰值会在哪里，所以我不能只分割积分范围并对各部分求和（除非有人对如何做到这一点有很好的想法）。

谢谢！

输出示例：

>>>scipy.integrate.quad(weighted_ftag_2, 0, 10000)
(0.0, 0.0)
>>>scipy.integrate.quad(weighted_ftag_2, 0, 1602)
(0.0, 0.0)
>>>scipy.integrate.quad(weighted_ftag_2, 1602, 1618)
(3.2710994652983256, 3.6297354011338712e-014)
>>>scipy.integrate.quad(weighted_ftag_2, 1618, 10000)
(0.0, 0.0)

Answer 1

您可能想尝试其他集成方法，例如

integrate.romberg()

方法。

或者，您可以使用

weighted_ftag_2(x_samples).argmax()

获取函数较大点的位置，然后使用一些启发式方法来缩短函数最大值周围的积分区间（位于

x_samples[….argmax()]

处。您必须泰勒您的问题的采样横坐标列表 (

x_samples

)：它必须始终包含位于函数最大值区域的点。

更一般地说，有关要集成的功能的任何具体信息都可以帮助您获得其积分的良好价值。我会将一种适合您的函数的方法（Scipy 提供的许多方法之一）与积分区间的合理分割（例如沿着上面建议的路线）结合起来。

Answer 2

如何在每个整数范围 [x, x+1) 上评估函数 f()， 并相加，例如

romb()

，正如 EOL 所建议的，当它 > 0 时：

from __future__ import division
import numpy as np
from scipy.integrate import romb

def romb_non0( f, a=0, b=10000, nromb=2**6+1, verbose=1 ):
    """ sum romb() over the [x, x+1) where f != 0 """
    sum_romb = 0
    for x in xrange( a, b ):
        y = f( np.arange( x, x+1, 1./nromb ))
        if y.any():
            r = romb( y, 1./nromb )
            sum_romb += r
            if verbose:
                print "info romb_non0: %d %.3g" % (x, r)  # , y
    return sum_romb

#...........................................................................
if __name__ == "__main__":
    np.set_printoptions( 2, threshold=100, suppress=True )  # .2f

    def f(x):
        return x if (10 <= x).all() and (x <= 12).all() \
            else np.zeros_like(x)

    romb_non0( f, verbose=1 )