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 )