Integrating a function with singularities using scipy's quad routine

Question

I am using the quad function from scipy.integrate v0.19.1 to integrate functions with a square-root like singularity at each end of the integration interval such as for example

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

(I use the sqrt function from numpy v1.12.0) which immediately yields the correct result pi:

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

According to the documentation of the quad function the keyword points should be used to indicate the locations of singularities or discontinuities of the integrand, but if I indicate the points [1, -1] where the above integrand is singluar I get a Warning and nan as result:

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)

Can someone clarify, why quad produces these problems if the singularities of the integrand are specified and just runs fine if these points are not indicated?

EDIT: I think I figured out the correct solution for this problem. For the case someone else encounters similar problems I quickly want to share my findings:

I you want to integrate a function of the form f(x)*g(x) with a smooth function f(x) and g(x) = (x-a)**alpha * (b-x)**beta, where a and b are the the integration limits and g(x) has singularities at these limits if alpha, beta < 0, then you are supposed to just integrate f(x) using g(x) as weighting function. For the quad routine, this is possible using the weight and wvar arguments. With these arguments you are also able to handle singularities of different kinds and problematic oscillatory behavior. The weighting function g(x) defined above can be used by setting weight='alg' and using wvar=(alpha, beta) to specify the exponents in g(x).

Since 1/sqrt(1-x**2) = (x+1)**(-1/2) * (1-x)**(-1/2) I can now handle the integral as follows:

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

which yields the correct answer pi with a very high accuracy, no matter if I use the argument points=(-1, 1) or not (which, as far as I understand now, should only be used, if the singularities/discontinuities can not be handled by choosing an appropriate weighting function).

Answer 1

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

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

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