Python+Scipy+Integration: dealing with precision errors in functions with spikes
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).
>>>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)
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).
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 )