cumulative simpson integration with scipy

numerical-integration, numpy, scipy

Solution

You can always write your own:

def cumsimp(func,a,b,num):
    #Integrate func from a to b using num intervals.

    num*=2
    a=float(a)
    b=float(b)
    h=(b-a)/num

    output=4*func(a+h*np.arange(1,num,2))
    tmp=func(a+h*np.arange(2,num-1,2))
    output[1:]+=tmp
    output[:-1]+=tmp 
    output[0]+=func(a)
    output[-1]+=func(b)
    return np.cumsum(output*h/3)

def integ1(x):
    return x

def integ2(x):
    return x**2

def integ0(x):
    return np.ones(np.asarray(x).shape)*5

First look at the sum and derivative of a constant function.

print cumsimp(integ0,0,10,5)
[ 10.  20.  30.  40.  50.]

print np.diff(cumsimp(integ0,0,10,5))
[ 10.  10.  10.  10.]

Now check for a few trivial examples:

print cumsimp(integ1,0,10,5)
[  2.   8.  18.  32.  50.]

print cumsimp(integ2,0,10,5)
[   2.66666667   21.33333333   72.          170.66666667  333.33333333]

Writing your integrand explicitly is much easier here then reproducing the simpson's rule function of scipy in this context. Picking intervals will be difficult to do when provided a single array, do you either:

- Use every other value for the edges of simpson's rule and the remaining values as centers?

- Use the array as edges and interpolate values of centers?

There are also a few options for how you want the intervals summed. These complications could be why its not coded in scipy.

Problem

I have some code which uses scipy.integration.cumtrapz to compute the antiderivative of a sampled signal. I would like to use Simpson's rule instead of Trapezoid. However scipy.integration.simps seems not to have a cumulative counterpart... Am I missing something? Is there a simple way to get a cumulative integration with "scipy.integration.simps"?

Original source