使用 scipy 的 quad 例程对具有奇点的函数进行积分
Integrating a function with singularities using scipy's quad routine
我正在使用 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)
这会以非常高的准确性产生正确答案 pi，无论我是否使用参数 points=(-1, 1) （据我现在了解，应该只使用，如果 singularities/discontinuities 无法通过选择合适的加权函数来处理）。
参数points是针对在积分区间内出现的singularities/discontinuities。它不适用于间隔的端点。因此，在您的示例中，没有 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 则输出毫无价值。
我正在使用 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)
这会以非常高的准确性产生正确答案 pi，无论我是否使用参数 points=(-1, 1) （据我现在了解，应该只使用，如果 singularities/discontinuities 无法通过选择合适的加权函数来处理）。
参数points是针对在积分区间内出现的singularities/discontinuities。它不适用于间隔的端点。因此，在您的示例中，没有 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 则输出毫无价值。