Scipy的四边形在“相当平滑”功能上失败

Question

数据生成函数相当混乱，所以我会尽量弄清楚。如果做得不够好，请评论。

我有一个函数，有几个变量，试着积分，另一个不变。然而，对于第二变量的不同值，积分过程产生完全不同的结果-尽管函数没有太大改变！在

这是我的示例代码（unfort。不可复制）：

fig, ax = plt.subplots()
tPrimeGrid = [0.16, 0.18, 0.22]
from scipy import integrate
for ttPrime in tPrimeGrid:
    # compute grid
    lowerBar = continuousTime.getPLowerBarAnalytical(ttPrime, Param.S, Param)
    upperBar = Param.r

    # integrate function
    h = integrate.quad(myFunc, lowerBar, upperBar, args=(0.6, ttPrime, Param), full_output=True)
    print('tPrime {} value {} erorr {}'.format(ttPrime, h[0], h[1]))
    # plot function to be evaluated
    pGrid = np.linspace(lowerBar, upperBar, 100)
    plt.plot(pGrid, [myFunc(pp, 0.6, ttPrime, Param) for pp in pGrid], label='t {} (h: {})'.format(ttPrime, h[0]))
plt.legend()

quad的输出是

^{pr2}$

所有的错误都在我的容忍范围之内，但是值跳了很多。图表如下：

所以这些函数的形状非常相似。用肉眼看，积分的差别确实没有意义。于是我“手动”计算了它们：

^{3}$

意味着这个小值实际上是正确的。因此，我还将在这里添加quad()的完整输出：

integrate.quad(myFunc, lowerBar, upperBar, args=(0.6, 0.16, Param), full_output=True)
Out[19]: 
(6.634157704675579,
 0.004721148834250418,
 {'alist': array([ 1.99994895,  1.78826738,  1.86060326,  1.79090489,  1.93030163,
          1.72120652,  1.96515082,  1.75605571,  1.98257541,  1.7734803 ,
          1.9912877 ,  1.7821926 ,  1.99564385,  1.78872682,  1.99782193,
          1.78654874,  1.99940018,  1.78763778,  1.99945548,  1.99891096,
          1.78845456,  1.99972774,  1.99918322,  1.78831843,  1.99988939,
          1.99931935,  1.7881823 ,  1.99999149,  1.99942145,  1.7882844 ,
          1.9998979 ,  1.9999447 ,  1.99940443,  1.78825036,  1.99996597,
          1.99986387,  1.99991492,  1.99995746,  1.99938742,  1.78827589,
          1.99993194,  1.99990641,  1.99998298,  1.99988089,  1.99995321,
          1.99939592,  1.78827164,  1.99994044,  1.99995108,  1.99940231]),
  'blist': array([ 1.99995108,  1.78827164,  1.93030163,  1.86060326,  1.96515082,
          1.75605571,  1.98257541,  1.7734803 ,  1.9912877 ,  1.7821926 ,
          1.99564385,  1.78654874,  1.99782193,  1.79090489,  1.99891096,
          1.78763778,  1.99940231,  1.7881823 ,  1.99972774,  1.99918322,
          1.78872682,  1.99986387,  1.99931935,  1.78845456,  1.9998979 ,
          1.99938742,  1.78825036,  2.        ,  1.99945548,  1.78831843,
          1.99990641,  1.99994895,  1.99942145,  1.78826738,  1.99998298,
          1.99988089,  1.99993194,  1.99996597,  1.99939592,  1.7882844 ,
          1.99994044,  1.99991492,  1.99999149,  1.99988939,  1.99995746,
          1.99940018,  1.78827589,  1.9999447 ,  1.99995321,  1.99940443]),
  'elist': array([  4.61134496e-03,   4.14135579e-06,   1.27804393e-15,
           1.55808958e-15,   5.54429458e-16,   0.00000000e+00,
           2.90617202e-16,   0.00000000e+00,   2.27720143e-15,
           0.00000000e+00,   3.91514838e-15,   0.00000000e+00,
           2.15987477e-14,   5.40749179e-17,   1.19682144e-12,
           0.00000000e+00,   1.03521880e-04,   0.00000000e+00,
           2.48275100e-08,   6.85735811e-13,   6.78454062e-18,
           8.65204618e-09,   1.64268148e-13,   3.39437556e-18,
           4.65487056e-08,   1.12644353e-13,   0.00000000e+00,
           6.19443638e-08,   4.08066769e-09,   8.48813317e-19,
           5.99147851e-07,   1.21478039e-06,   9.39741081e-10,
           0.00000000e+00,   6.32087778e-14,   2.19912539e-10,
           6.97448944e-09,   1.85114515e-13,   1.28558440e-13,
           2.12217047e-19,   8.89451362e-08,   8.45167083e-09,
           1.25621961e-14,   2.59393721e-08,   2.45738052e-15,
           1.19106919e-12,   1.06110581e-19,   4.91812027e-08,
           2.72831861e-15,   3.06819516e-12]),
  'iord': array([ 1, 17, 50, 49, 48, 47, 46, 45, 44, 43, 42, 29, 40, 39, 38, 23, 35,
         11, 34,  3,  7, 21, 30, 18, 12,  8,  6,  0,  0,  0,  0,  0,  0,  0,
          0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0], dtype=int32),
  'last': 50,
  'neval': 2079,
  'rlist': array([  1.04205921e-01,   6.52426361e-06,   1.15115963e-01,
           1.40340233e-01,   4.99385660e-02,   0.00000000e+00,
           2.33230318e-02,   0.00000000e+00,   1.12779305e-02,
           0.00000000e+00,   5.54633015e-03,   0.00000000e+00,
           2.75040018e-03,   4.87063561e-03,   1.36955727e-03,
           0.00000000e+00,   8.85474806e-03,   0.00000000e+00,
           1.73731989e+00,   3.41803845e-04,   6.11097092e-04,
           1.73678061e+00,   1.70814248e-04,   3.05738170e-04,
           2.00530044e-01,   8.53852175e-05,   0.00000000e+00,
           5.29681716e-06,   1.51991445e-01,   7.64543067e-05,
           2.17985628e-01,   2.00512965e-01,   7.26773425e-02,
           0.00000000e+00,   1.06564581e-05,   3.34545761e-01,
           5.59012296e-01,   5.32838374e-06,   1.06721256e-05,
           1.91148122e-05,   3.34512038e-01,   2.38773007e-01,
           5.32807434e-06,   1.85664300e-01,   2.66423055e-06,
           5.33597722e-06,   9.55759147e-06,   1.85647516e-01,
           1.33212494e-06,   8.93844378e-03])},
 'The maximum number of subdivisions (50) has been achieved.\n  If increasing the limit yields no improvement it is advised to analyze \n  the integrand in order to determine the difficulties.  If the position of a \n  local difficulty can be determined (singularity, discontinuity) one will \n  probably gain from splitting up the interval and calling the integrator \n  on the subranges.  Perhaps a special-purpose integrator should be used.')

下面是我的相关问题：

1. 为什么在这样相似的形状下，有些参数起作用，有些参数不起作用？我如何“预测”未来的这些失败？在
2. 完整的输出（下面）告诉我已经实现了最大数量的细分。它并没有告诉我错误可能是未估计到的。那一定是暗示？5.9-0.0004不接近0.3。在
3. 由于增加限额没有帮助（见下文），有哪些潜在的替代方案？如何集成此功能？在

我尝试增加limit，但这只给出了以下（不是更好的）输出：

integrate.quad(expectedUtility, lowerBar, upperBar, args=(0.6, 0.16, Param), full_output=True, limit=500)
Out[24]: 
(6.633112814769514,
 4.74743687572826e-06,
[...]
 'The occurrence of roundoff error is detected, which prevents \n  the requested tolerance from being achieved.  The error may be \n  underestimated.')

Answer 1

当full_output为真时，quad对于返回值有一个有些奇怪的约定。如果quad没有检测到任何问题，则返回y, abserr, infodict。如果quad检测到某种形式的故障，则返回y, abserr, infodict, message。infodict不包含表示失败的简单status字段，因此您必须检查是否存在第四个返回值。如果存在，则是一个描述问题的字符串。（在代码中，您可能会添加类似if len(h) > 3: ...handle the failure...）在错误情况的完整输出中，您可以看到message是包含以下内容的字符串：

The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze
the integrand in order to determine the difficulties.  If the position of a
local difficulty can be determined (singularity, discontinuity) one will
probably gain from splitting up the interval and calling the integrator
on the subranges.  Perhaps a special-purpose integrator should be used.

这意味着数值积分失败了。为什么不知道被积函数很难。在