scipy quad seems to provide inaccurate integral value

Question

I am puzzled by the fact that integrate.quad from scipy returns inaccurate results even for simple integrals. In the following, I will give a simple example for which it is easily possible to cross-check with analytical computation.

However, I would like to know a general resolution to the issue by knowing what is the best and the most reliable integral method in scipy/numpy. It is because, the actual integral I need to do is highly complicated, and an analytical solution and/or re-parametrization is not possible.

Here is a simple example:

from scipy.integrate import quad
xmax=1.0e23
xmin=0.5
def myFUN(x):
    return 1/(1+x)**(3.0)
print( quad( myFUN, xmin,xmax )[0] )

that returns:

2.3171047961935996e-40

Now, let's do it analytically:

def ana(x):
    return -0.5/(1+x)**2.0 
print( ana(xmax)-ana(xmin) )

one gets:

0.2222222222222222

Note that the numerical output differed from the analytical computation by almost 40-orders of magnitude.

Thank you in advance.

Answer 1

Imagine that you were doing the problem by, say, the trapezium rule. You break the range up into, say, 1000 intervals. If xmin = 0.5 and xmax = 1e23 then the first internal value is 1e20 and the value of the function is about 1e-60 or, to all intents and purposes, 0. Obviously scipy's quad routine will use something more advanced than trapezium rule, but the consequences will be the same.

So, if you choose to use a library routine with "near-infinity" as one of the arguments then you would have to take special measures (e.g. truncate the integral range to something where you can bound the part that is ignored below a tolerable level).

If you use a sensible integration range then the function is well-behaved.

from scipy.integrate import quad
xmax=10.0
xmin=0.5

def myFUN(x):
    return 1/(1+x)**(3.0)

def indefinite_integral( x ):
    return -0.5 / ( 1 + x ) ** 2

print( "Using library routine: ", quad( myFUN, xmin,xmax )[0] )
print( "Using your maths knowledge: ", indefinite_integral( xmax ) - indefinite_integral( xmin ) )

Output:

Using library routine:  0.21808999081726355
Using your maths knowledge:  0.21808999081726355

EDIT:

You can usually transform the integral so that the integral range is O(1). For example, in your case, you could substitute

u = 1/(1+x)³

Then f(x).dx transforms to (-1/3)u^-1/3.du.

In the code below I have gone back to your original integral limits, but also made this substitution. (I also flipped the integral limits over to make the integrand positive, but that's not vital.) Then it gives you the answer you desire on the full range.

from scipy.integrate import quad
xmax=1e23
xmin=0.5

def myFUN( x ):
    return 1.0 / ( 1.0 + x ) ** 3.0

def indefinite_integral( x ):
    return -0.5 / ( 1 + x ) ** 2

print( "Using library routine: ", quad( myFUN, xmin,xmax )[0] )
print( "Using your maths knowledge: ", indefinite_integral( xmax ) - indefinite_integral( xmin ) )

######################

def altFUN( u ):
    return 1.0 / ( 3.0 * u ** ( 1.0 / 3.0 ) )

umin = myFUN( xmax )
umax = myFUN( xmin )
print( "Using transformed integral: ", quad( altFUN, umin, umax )[0] )

Output:

Using library routine:  2.3171047961935996e-40
Using your maths knowledge:  0.2222222222222222
Using transformed integral:  0.2222222222222222

Answer 2

I wrote a Python integration package called lintegrate that works on the logarithm of functions, which you could potentially use instead. After installing lintegrate with:

pip install lintegrate

you could do:

import numpy as np

from lintegrate import lcquad

def myfun(x, *args):
    """
    Define the natural logarithm of the function you want to integrate.
    """

    # natural log of 1 / (1 + x)**3
    return -3.0 * np.log(1 + x)


xmin = 0.5
xmax = 1e23

# return the natural log of the original function
lint = lcquad(myfun, xmin, xmax)[0]

# exponentiate to get your required answer
print(f"Integral = {np.exp(lint)}")
Integral = 0.2222222222800475

You can see the documentation for lcquad here. Internally it use the GSL cquad function.