scipy.integrate.quad当要集成的函数也是整数时失败（有时）

Question

我用sqipy.integrate.quad来计算二重积分。基本上，我试图计算exp[-mu_wx_par]上的积分，其中mu_wx_par也是一个积分。在

我的代码基本上是有效的。但是，它返回的值不正确。在

import numpy as np
from scipy import integrate


def mu_wx_par(x, year, par):
    """ First function to be integrated """
    m = max(par['alfa'], 0) + par['beta'] * 10 ** (par['gamma'] * x)
    w = np.minimum(par['frem_a'] + par['frem_b'] * x + par['frem_c'] * x**2, 0)
    return m * (1 + w/100)**(year - par['innf_aar'])


def tpx_wx(x, t, par, year):
    """ Second function to be integrated (which contains an integral itself)"""
    result, _ = integrate.quad(lambda s: mu_wx_par(x + s, year + s, par), 0, t)
    return np.exp(-result)


def est_lifetime(x, year, par):
    """ Integral of second function. """
    result, _ = integrate.quad(lambda s: tpx_wx(x, s, par, year), 0, 125 - x)

    return result


# Test variables
par = {'alfa': 0.00019244401470947973,
       'beta': 2.420260552210541e-06,
       'gamma': 0.0525500987420195,
       'frem_a': 0.3244611019518985,
       'frem_b': -0.12382978382606026,
       'frem_c': 0.0011901237463116591,
       'innf_aar': 2018
       }


year = 2018
estimate_42 = est_lifetime(42, year, par)
estimate_43 = est_lifetime(43, year, par)
rough_estimate_42 = sum([tpx_wx(42, t, par, year) for t in range(0, 100)])

print(estimate_42)
print(estimate_43)
print(rough_estimate_42)
3.1184634065887544
46.25925442287578
47.86323490659588

estimate_42的值不正确。它应该与rough_estimate_42的值相同。但是请注意，estimate_43看起来不错。这是怎么回事？在

我使用的是scipyv1.1.0和numpyv1.15.1和Windows。在

有人认为函数几乎在任何地方都接近于零，这篇文章scipy integrate.quad return an incorrect value。事实并非如此，简单的tpx_wx从a=0到{}的简单绘图清楚地显示了这一点

^{pr2}$

Answer 1

这似乎与a known issue中某些Fortran代码在quad中编译的方式有关：递归地调用它在某些情况下会导致失败。另请参见Big integration error with integrate.nquad。在

除非用更好的标志重新编译SciPy，在Windows上似乎应该避免与quad进行嵌套集成。一种解决方法是在其中一个集成步骤中使用^{}方法。将est_lifetime中的quad替换为

integrate.romberg(lambda s: tpx_wx(x, s, par, year), 0, 125 - x, divmax=20)

结果为estimate_42生成{}，这与quad在Linux上返回的结果一致。在

---

可视化集成过程的一种方法是声明一个全局列表eval_points，并将eval_points.append(t)插入tpx_wx。使用相同版本的SciPy（本测试中为0.19.1），结果看起来不同。在

在Linux上：

在Windows上：

在窗口上，一个复杂点的邻域的迭代二分法过早终止，结果似乎是在子区间上的某个部分积分。在