Function integration with singularities using scipy quad procedure

I am using a function quad

from scipy.integrate v0.19.1

to integrate functions with a square root singularity at each end of the integration interval, like for example

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

      

        
        
        
      

    

(I'm using the function sqrt

from numpy v1.12.0

) which immediately gives the correct pi result:

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

      

        
        
        
      

    

According to the function's documentation, the quad

keyword points

should be used to indicate the locations of singularities or breaks of the integrand, but if I specify the points [1, -1]

where the specified integrand is singluar, we get a Warning and nan

as a 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 elaborate on why this is quad

causing these problems if the special features of the integrand are specified and just works fine if these points are not specified?

EDIT: I think I figured out the correct solution to this problem. For the case when someone is facing similar problems, I quickly want to share my findings:

I want to integrate a function of the kind 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 limits of integration, and g(x)

has features within these limits if alpha, beta < 0

, then you should simply integrate f(x)

using g(x)

as the weighting function. For a subroutine, quad

this is possible using the arguments weight

and wvar

. With these arguments, you can also deal with different types of features and problematic oscillatory behavior. The ordered function g(x)

defined above can be used by setup weight='alg'

and use wvar=(alpha, beta)

to specify metrics in g(x)

.

Since 1/sqrt(1-x**2) = (x+1)**(-1/2) * (1-x)**(-1/2)

now I can handle the integral like this:

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

      

        
        
        
      

    

which gives the correct answer pi

with very high precision, whether I use the argument points=(-1, 1)

or not (which, as I understand it now, should only be used if the features / breaks might not be handled by choosing the appropriate weighting function).