使用 scipy 的 quad 例程将函数与奇点集成

Question

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

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

（我使用来自numpy v1.12.0 的sqrt 函数）立即产生正确的结果pi：

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

根据quad 函数的文档，关键字points 应用于指示被积函数的奇点或不连续点的位置，但如果我指出上述被积函数为单点的点[1, -1]，我会得到结果是警告和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，那么您应该使用g(x) 作为权重函数 来整合f(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)

无论我是否使用参数points=(-1, 1)，它都会以非常高的准确性产生正确的答案pi（据我所知，只有在奇点/不连续点可以不能通过选择适当的加权函数来处理）。

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 则输出毫无价值。