3

I have a problem with scipy.integrate.quad. In short, I have a really long and complicated set of nested functions and integrals which include an integral of a decreasing function which must be integrated over the specific range 10^2 < x < 10^20.

To demonstrate this problem simply, consider the integral of y=x^(-2) between these values using scipy.integrate.quad:

import numpy as np
from scipy.integrate import quad

def func(x,z):
    "decreasing function to test"
    # print x, x**-2.
    return x**(-2.)

# Test quad:
print('Quad', quad(lambda x: func(x, 0.), 1e2, 1e20)[0])

The true answer to this integral is 0.01, however quad returns a very small value, and by uncommenting the print line in the function you can see that it only integrates over the largest x values (or is biased that way due to the log scale of the x-axis). I need to find some way of fixing this problem.

I know that you can get the correct answer by using other methods such as the Simpsons or trapezoidal rules:

from scipy.integrate import simpson, trapezoid

xarray=np.logspace(2, 20, 1000)

print('Simpson logspace x', simpson(func(xarray, 0.), xarray))
print('Trapezoid logspace x', trapezoid(func(xarray, 0.), xarray))

However, these include passing an array into the function, which is not possible with my actual code. Thus I'd have to include a for loop to generate a y array to integrate with, which slows the whole program down unacceptably.

Is there any trick to making quad work for this sort of x range, or anything which works like quad which I can use instead?

5
  • 2
    quad(lambda x: func(x,0.),1e2,np.inf) will work.
    – cel
    Oct 25, 2015 at 9:04
  • And do you really need a for loop to generate y for your x? Why don't your functions support element-wise array operations? (Anyway the quad approach should be more accurate if you get it to work properly) Oct 25, 2015 at 10:57
  • Thanks guys for replying. Unfortunately cel, for my function quad(lambda x: func(x,0.),1e2,np.inf) doesn't help. It returns a negative (wrong!) number which doesn't match the simpson or trapezium rules between 10^2 and 10^20. Andras, the reason my function doesn't support element-wise operation at the moment is because it loads in arrays at the base level and I haven't programmed it to handle taking arrays at the top. I may have to rewrite it if needs be, but it'll take ages, and I hoped there was a simpler fix to quad.
    – SeaWalk
    Oct 25, 2015 at 18:41
  • Does that negative return value refer to the 1/x^2 example? If the negative is for your actual problem: are you sure it's wrong? And have you tried, just for proof-of-concept, using a list comprehension instead of a for loop to generate your data for trapz (to see whether that's marginally faster)? And next time you might want to ping cel by adding @cel to your comment, otherwise they won't get a notification (I just left this question open just for such an occasion). Oct 25, 2015 at 18:59
  • 1
    There's also the points argument that you can use if you have some feeling about how the function behaves.
    – cel
    Oct 25, 2015 at 19:32

1 Answer 1

1

Try the variable substitution x = exp(t).

import numpy as np
from scipy.integrate import quad

# original
def f(x):
    return x**-2.

a, b = 10**2., 10**20.
quad(f, a, b)
# (2.754484024437256e-17, 5.2212478622117524e-17)

# with substitution
def g(t):
    expt = np.exp(t)
    return f(expt) * expt

quad(g, np.log(a), np.log(b))
# (0.00999999999999998, 3.0000387625246216e-09)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.