CODEWITHSUNDEEP
pythonscipyquad
TMueller83
TMueller83

Reputation: 612

Integrating a function with singularities using scipy's quad routine

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).

Upvotes: 9

Views: 5330

Answers (1)

user6655984
user6655984

Reputation:

The parameter points is meant for singularities/discontinuities that occur within the interval of integration. It is not intended for endpoints of the interval. So, in your example the version without points is the correct approach. (I can't pinpoint what goes wrong when endpoints are included in points without diving into FORTRAN code that SciPy wraps.)

Compare with the following example, where singularities occur within the interval of integration: 

>>> 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)

Here the inclusion of points is appropriate, and yields the correct result, while without points the output is worthless.

Upvotes: 4

Related Questions