Python+Scipy+积分：处理带尖峰函数的精度误差

Question

我正在尝试使用scipy.integrate.quad这个工具来计算一个函数在一个非常大的范围内的积分（从0到10,000）。这个函数在大部分范围内都是零，但在一个很小的范围内有一个尖峰（比如在1,602到1,618之间）。

在进行积分时，我本来期待结果是正数，但我猜这个quad的猜测算法可能搞混了，结果输出了零。我想知道有没有办法解决这个问题（比如使用不同的算法，或者调整一些参数等等）？我通常不知道尖峰会在哪里，所以不能简单地把积分范围分开然后把各部分相加（除非有人有好的主意来做到这一点）。

谢谢！

示例输出：

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

你可以试着在每个整数范围 [x, x+1) 上评估你的函数 f()，然后把那些大于 0 的结果加起来，比如用 romb()，就像 EOL 提到的那样。

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 )

Answer 2

你可以尝试其他的积分方法，比如 integrate.romberg() 这个方法。

另外，你可以通过 weighted_ftag_2(x_samples).argmax() 找到你函数值比较大的那个点的位置，然后用一些简单的技巧来缩小积分的区间，围绕着你函数的最大值进行积分（这个最大值的位置是 x_samples[….argmax()]）。你需要根据自己的问题调整样本点列表（x_samples），确保它总是包含在你函数值最大的区域里的点。

更一般来说，关于要积分的函数的任何具体信息都能帮助你得到一个好的积分值。我建议结合一个适合你函数的方法（可以参考 Scipy 提供的多种方法）和合理的积分区间划分（比如上面提到的那种方式）。