Scipy.integrate 给出奇怪的结果

Question

我仍在与 scipy.integrate.quad 作斗争。

抛开所有细节，我有一个积分需要评估。该函数的形式是 x 中函数的乘积的积分形式，如下所示：

Z(k) = f(x) g(k/x) / abs(x)

我确定积分的范围是在两个正数之间。奇怪的是，当我选择一个我知道必须包含 x 的所有正值的宽范围时（例如从 1 到 10,000,000 的积分），它会快速积分并给出看起来正确的答案。但是，当我找出确切的限制（我知道 f(x) 在很多实数线上为零）并使用这些限制时，我得到了另一个不同的答案。它们并没有太大不同，尽管我知道第二个更准确。

经过多次摆弄，我让它工作正常，但随后需要添加一个除法 - 我至少得到了 z 的计算函数的“平滑”答案。在添加求幂（这是必需的）之前，我已经以正常的方式工作了，但是现在生成的函数（z）变得越来越振荡和奇特。

知道这里发生了什么吗？我知道这段代码来自一个旧的 Fortran 库，所以肯定存在一些已知问题，但我找不到参考文献。

核心代码如下：

def normal(x, mu, sigma) :
    return (1.0/((2.0*3.14159*sigma**2)**0.5)*exp(-(x-
                                      mu)**2/(2*sigma**2)))

def integrand(x, z, mu, sigma, f) : 
return np.exp(normal(z/x, mu, sigma)) * getP(x, f._x, f._y) / abs(x)




    for _z in range (int(z_min), int(z_max) + 1, 1000):
        z.append(_z)
        pResult = quad(integrand, lb, ub,
                       args=(float(_z), MU-SIGMA**2/2, SIGMA, X),
                       points = [100000.0],
                       epsabs = 1, epsrel = .01)    # drop error estimate of tuple 
        p.append(pResult[0])   # drop error estimate of tuple 

顺便说一下，getP() 返回一个线性插值、分段连续但非平滑的函数，以给出平滑地拟合直方图的离散“桶”之间的积分器值。

Answer 1

与许多数值方法一样，它可能对渐近线、零等非常敏感。唯一的选择是如果它接受的话，继续给它“提示”。