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

问题描述 投票：0回答：1

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

我有一个包含多个变量的函数，并尝试对它进行积分，同时保持另一个变量不变。然而，对于次要变量的不同值，积分过程会产生完全不同的结果 - 尽管函数变化不大！

这是我的示例代码（不幸的是，不可重现）：

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
的输出是

tPrime 0.16 value 6.63310536371 erorr 0.000345564616621
tPrime 0.18 value 5.93645658492 erorr 0.00045956820487
tPrime 0.22 value 0.359208009237 erorr 3.98801002485e-15

错误都在我的容忍范围之内，但数值跳动很大。这是图表：

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

(np.array([myFunc(pp, 0.6, ttPrime, Param) for pp in pGrid])).sum()*(pGrid[1]-pGrid[0])
0.35538925924926557

这意味着较小的值实际上是正确的。因此，我还将在此处添加 

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.')
python numpy scipy numerical-integration
1个回答
2
投票
当

quad

 为 True 时，
full_output
对于返回值有一个有点时髦的约定。如果 
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.

这意味着数值积分失败。如果不了解更多关于被积数

myFunc
的信息，很难说它失败的原因。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.