python - 使用 scipy 的 quad 例程集成具有奇点的函数

Question

我正在使用 quad来自 scipy.integrate v0.19.1 的函数在积分区间的每一端积分具有平方根奇点的函数，例如

In [1]: quad(lambda x: 1/sqrt(1-x**2), -1, 1)

(我使用 sqrt 中的 numpy v1.12.0 函数)立即产生正确的结果 pi:

Out[1]: (3.141592653589591, 6.200897573194197e-10)

根据 quad 的文档函数关键字points应该用于指示被积函数的奇点或不连续点的位置，但如果我指示点 [1, -1]上面的被积函数是 singluar 我得到警告和 nan结果:

In [2]: quad(lambda x: 1/sqrt(1-x**2), -1, 1, points=[-1, 1])

RuntimeWarning: divide by zero encountered in double_scalars
IntegrationWarning: Extremely bad integrand behavior occurs at some
points of the integration interval.

Out[2]: (nan, nan)

有人可以澄清一下，为什么quad如果指定了被积函数的奇点，是否会产生这些问题，如果未指定这些点，则运行正常？

编辑:我想我找到了解决这个问题的正确方法。对于其他人遇到类似问题的情况，我想很快分享我的发现:

我想集成 f(x)*g(x) 形式的函数具有平滑功能f(x)和 g(x) = (x-a)**alpha * (b-x)**beta , 其中a和 b是集成限制和 g(x)如果 alpha, beta < 0 在这些极限处有奇点, 那么你应该只整合 f(x)使用 g(x)作为权重函数。对于 quad例程，这可以使用 weight和 wvar争论。有了这些参数，您还可以处理不同种类的奇点和有问题的振荡行为。加权函数g(x)上面定义的可以通过设置weight='alg'来使用并使用 wvar=(alpha, beta)在 g(x) 中指定指数.

自 1/sqrt(1-x**2) = (x+1)**(-1/2) * (1-x)**(-1/2)我现在可以按如下方式处理积分:

In [1]: quad(lambda x: 1, -1, 1, weight='alg', wvar=(-1/2, -1/2))
Out[1]: (3.1415926535897927, 9.860180640534107e-14)

产生正确答案 pi具有非常高的准确性，无论我是否使用参数 points=(-1, 1)或不(据我所知，只有在无法通过选择适当的加权函数来处理奇点/不连续点时才应使用)。

Answer 1

参数 points 表示出现在 积分区间内的奇点/不连续点。它不适用于间隔的端点。因此，在您的示例中，没有 points 的版本是正确的方法。 (如果不深入研究 SciPy 包装的 FORTRAN 代码，我无法确定当端点包含在 points 中时出了什么问题。)

与以下示例比较，其中奇点出现在积分区间内:

>>> quad(lambda x: 1/sqrt(abs(1-x**2)), -2, 2)
(inf, inf)
>>> quad(lambda x: 1/sqrt(abs(1-x**2)), -2, 2, points = [-1, 1])
(5.775508447436837, 7.264979728915932e-10)

这里包含 points 是合适的，并且会产生正确的结果，而没有 points 则输出毫无值(value)。