cumulative simpson integration with scipy

Question

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"?

Answer 1

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.

Answer 2

Your question has been answered a long time ago, but I came across the same problem recently. I wrote some functions to compute such cumulative integrals for equally spaced points; the code can be found on GitHub. The order of the interpolating polynomials ranges from 1 (trapezoidal rule) to 7. As Daniel pointed out in the previous answer, some choices have to be made on how the intervals are summed, especially at the borders; results may thus be sightly different depending on the package you use. Be also aware that the numerical integration may suffer from Runge's phenomenon (unexpected oscillations) for high orders of polynomials.

Here is an example:

import numpy as np
from scipy import integrate as sp_integrate
from gradiompy import integrate as gp_integrate

# Definition of the function (polynomial of degree 7)
x = np.linspace(-3,3,num=15)
dx = x[1]-x[0]
y = 8*x + 3*x**2 + x**3 - 2*x**5 + x**6 - 1/5*x**7
y_int = 4*x**2 + x**3 + 1/4*x**4 - 1/3*x**6 + 1/7*x**7 - 1/40*x**8

# Cumulative integral using scipy
y_int_trapz = y_int [0] + sp_integrate.cumulative_trapezoid(y,dx=dx,initial=0)
print('Integration error using scipy.integrate:')
print('  trapezoid = %9.5f' % np.linalg.norm(y_int_trapz-y_int))

# Cumulative integral using gradiompy
y_int_trapz = gp_integrate.cumulative_trapezoid(y,dx=dx,initial=y_int[0])
y_int_simps = gp_integrate.cumulative_simpson(y,dx=dx,initial=y_int[0])
print('\nIntegration error using gradiompy.integrate:')
print('  trapezoid = %9.5f' % np.linalg.norm(y_int_trapz-y_int))
print('  simpson   = %9.5f' % np.linalg.norm(y_int_simps-y_int))

# Higher order cumulative integrals
for order in range(5,8,2):
    y_int_composite = gp_integrate.cumulative_composite(y,dx,order=order,initial=y_int[0])
    print('  order %i   = %9.5f' % (order,np.linalg.norm(y_int_composite-y_int)))
    
# Display the values of the cumulative integral
print('\nCumulative integral (with initial offset):\n',y_int_composite)

You should get the following result:

'''
Integration error using scipy.integrate:
  trapezoid = 176.10502

Integration error using gradiompy.integrate:
  trapezoid = 176.10502
  simpson   =   2.52551
  order 5   =   0.48758
  order 7   =   0.00000

Cumulative integral (with initial offset):
 [-6.90203571e+02 -2.29979407e+02 -5.92267425e+01 -7.66415188e+00
  2.64794452e+00  2.25594840e+00  6.61937372e-01  1.14797061e-13
  8.20130517e-01  3.61254267e+00  8.55804341e+00  1.48428883e+01
  1.97293221e+01  1.64257877e+01 -1.13464286e+01]
'''

Answer 3

I would go with Daniel's solution. But you need to be careful if the function that you are integrating is itself subject to fluctuations. Simpson's requires the function to be well-behaved (meaning in this case, one that is continuous).

There are techniques for making a moderately badly behaved function look like it is better behaved than it really is (really forms of approximation of your function) but in that case you have to be sure that the function "adequately" approximates yours. In that case you might make the intervals may be non-uniform to handle the problem.

An example might be in considering the flow of a field that, over longer time scales, is approximated by a well-behaved function but which over shorter periods is subject to limited random fluctuations in its density.