2

I have a problem with scipy.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 numpy.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 simps
    from numpy import trapz

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

    print 'Simpson logspace x',simps(func(xarray,0.),xarray)
    print 'Trapezoid logspace x',trapz(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?

Will Vousden
  • 32,488
  • 9
  • 84
  • 95
SeaWalk
  • 63
  • 4
  • 2
    `quad(lambda x: func(x,0.),1e2,np.inf)` will work. – cel Oct 25 '15 at 09: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) – Andras Deak -- Слава Україні Oct 25 '15 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 '15 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](http://stackoverflow.com/a/22108640/5067311) 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). – Andras Deak -- Слава Україні Oct 25 '15 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 '15 at 19:32

0 Answers0