使用不同正交积分切比雪夫权重函数时的不同结果答案

【问题标题】：Different results while integrating Chebyshev weight function using different quadratures使用不同正交积分切比雪夫权重函数时的不同结果
【发布时间】：2014-10-07 00:50:15
【问题描述】：

有人可以解释我在使用 3 个不同的例程和 2 个不同的指数表示形式集成 Chebyshev weight function 时观察到的以下行为吗？在每种情况下，预期的答案都是 Pi：

from scipy.integrate import quadrature, quad, fixed_quad

print fixed_quad(lambda x: 1/(1 - x**2)**(1/2), -1, 1)
print fixed_quad(lambda x: 1/(1 - x**2)**(0.5), -1, 1)

print quadrature(lambda x: 1/(1 - x**2)**(1/2), -1, 1) 
print quadrature(lambda x: 1/(1 - x**2)**(0.5), -1, 1)

print quad(lambda x: 1/(1 - x**2)**(1/2), -1, 1)
print quad(lambda x: 1/(1 - x**2)**(0.5), -1, 1)

给了

(2.0000000000000009, None)
(2.8254100794589787, None)

(1.9999999999999996, 4.4408920985006262e-16)
/usr/lib/python2.7/dist-packages/scipy/integrate/quadrature.py:168: AccuracyWarning:    maxiter (50) exceeded. Latest difference = 6.965869e-04
 AccuracyWarning)
(3.107110439388189, 0.00069658693569163432)

(2.0, 2.220446049250313e-14)
(3.141592653589564, 6.200200353134733e-10)

首先，可以看出，答案取决于指数中给出的是 1/2 还是 0.5：为什么会这样？

其次，结果取决于选择的例程；有人可以解释为什么 FORTRAN 的 QUADPACK 得到正确的答案，但高斯求积完全错过了结果吗？

注意：here 解决了类似的问题，但没有解决上述两个具体问题。

【问题讨论】：

标签： python scipy integration numerical-methods

【解决方案1】：

在 python 2.7（您正在使用）中，整数除法是默认设置。这意味着 1/2 将计算为 0。如果要使用浮点除法，请将 from __future__ import division 添加到代码顶部。

【讨论】：

谢谢@gjdanis，这确实回答了我的第一个问题。
@MaviPranav：乐于助人！一直发生在我身上！

【解决方案2】：

正如@gjdanis 指出的那样，在python 2.7 中，1/2 是0（除非您在代码中包含from __future__ import division）。

你的被积函数在 1 和 -1 处有奇点。 fixed_quad 和 quadrature 使用权重函数 w(x) = 1 进行高斯求积，所以这些奇点处理不好。

fixed_quad 不是自适应的（因此得名）。默认顺序是 5。您必须增加顺序（很多）才能获得合理的近似值：

In [179]: print fixed_quad(lambda x: 1/(1 - x**2)**(0.5), -1, 1, n=100)
(3.124265558250825, None)

In [180]: print fixed_quad(lambda x: 1/(1 - x**2)**(0.5), -1, 1, n=2000)
(3.1407221810853478, None)

quadrature 只需按递增顺序调用fixed_quad（直到maxiter 参数给出的最大值），直到连续积分估计之间的差异足够小。打印的警告告诉您已达到最大订单但未满足所需的容错。 maxiter 的默认值为 50；您需要增加maxiter 以获得更好的结果。例如，这是maxiter=200 的结果：

In [2]: print quadrature(lambda x: 1/(1 - x**2)**(0.5), -1, 1, maxiter=200)
/Users/warren/local_scipy/lib/python2.7/site-packages/scipy/integrate/quadrature.py:183: AccuracyWarning: maxiter (200) exceeded. Latest difference = 4.353464e-05
  AccuracyWarning)
(3.1329074742380407, 4.3534643496823122e-05)

如果你使用maxiter，你也应该明智地使用miniter。 quadrature 天真地从 miniter 开始，并将阶数增加 1，直到误差估计足够小或达到 maxiter。

要详细了解fixed_quad 和quadrature 的工作原理，请查看源代码：https://github.com/scipy/scipy/blob/master/scipy/integrate/quadrature.py

正如您所指出的，quad 是 Fortran 库 QUADPACK 的包装器。这段代码比fixed_quad 和quadrature 使用的简单高斯求积要复杂得多。

【讨论】：