Python+Scipy+Integration: dealing with precision errors in functions with spikes

Question

Hi,

I am trying to use scipy.integrate.quad to integrate a function over a very large range (0..10,000). The function is zero over most of its range but has a spike in a very small range (e.g. 1,602..1,618).

When integrating, I would expect the output to be positive, but I guess that somehow quad's guessing algorithm is getting confused and outputting zero. What I would like to know is, is there a way to overcome this (e.g. by using a different algorithm, some other parameter, etc.)? I don't usually know where the spike is going to be, so I can't just split the integration range and sum the parts (unless somebody has a good idea on how to do that).

Thanks!

Sample output:

>>>scipy.integrate.quad(weighted_ftag_2, 0, 10000)
(0.0, 0.0)
>>>scipy.integrate.quad(weighted_ftag_2, 0, 1602)
(0.0, 0.0)
>>>scipy.integrate.quad(weighted_ftag_2, 1602, 1618)
(3.2710994652983256, 3.6297354011338712e-014)
>>>scipy.integrate.quad(weighted_ftag_2, 1618, 10000)
(0.0, 0.0)

Answer 1

You might want to try other integration methods, such as the integrate.romberg() method.

Alternatively, you can get the location of the point where your function is large, with weighted_ftag_2(x_samples).argmax(), and then use some heuristics to cut the integration interval around the maximum of your function (which is located at x_samples[….argmax()]. You must taylor the list of sampled abscissas (x_samples) to your problem: it must always contain points that are in the region where your function is maximum.

More generally, any specific information about the function to be integrated can help you get a good value for its integral. I would combine a method that works well for your function (one of the many methods offered by Scipy) with a reasonable splitting of the integration interval (for instance along the lines suggested above).

Answer 2

How about evaluating your function f() over each integer range [x, x+1), and adding up e.g. romb(), as EOL suggests, where it's > 0:

from __future__ import division
import numpy as np
from scipy.integrate import romb

def romb_non0( f, a=0, b=10000, nromb=2**6+1, verbose=1 ):
    """ sum romb() over the [x, x+1) where f != 0 """
    sum_romb = 0
    for x in xrange( a, b ):
        y = f( np.arange( x, x+1, 1./nromb ))
        if y.any():
            r = romb( y, 1./nromb )
            sum_romb += r
            if verbose:
                print "info romb_non0: %d %.3g" % (x, r)  # , y
    return sum_romb

#...........................................................................
if __name__ == "__main__":
    np.set_printoptions( 2, threshold=100, suppress=True )  # .2f

    def f(x):
        return x if (10 <= x).all() and (x <= 12).all() \
            else np.zeros_like(x)

    romb_non0( f, verbose=1 )