3

我正在尝试使用 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)
4

2 回答 2

3

您可能想尝试其他集成方法，例如integrate.romberg()方法。

或者，您可以使用 获取函数较大的点的位置，然后使用一些启发式方法将积分间隔缩短到函数最大值附近（weighted_ftag_2(x_samples).argmax()位于对于您的问题：它必须始终包含在您的功能最大的区域中的点。x_samples[….argmax()]x_samples

更一般地说，关于要集成的功能的任何特定信息都可以帮助您获得其积分的良好价值。我会结合一个对你的函数很好的方法（Scipy 提供的众多方法之一）和积分区间的合理分割（例如按照上面建议的线）。

于 2010-07-06T13:26:34.480 回答
1

如何在每个整数范围 [x, x+1) 上评估您的函数 f()，然后将 eg 相加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 )
于 2010-07-07T10:35:46.713 回答