Search code examples
pythonscipyphysicsnumerical-integration

Integration from a set of acceleration data to position


I'm trying for a project to integrate acceleration data in order to have an approximation of the position. I used a real simple set of data to start with, with a constant acceleration.

from scipy.integrate import cumtrapz

t = [0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7]
a = [-9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8]
v = cumtrapz(a, t, initial=0)
z = cumtrapz(v, t, initial=5)

The result is quite satisfying, apart from the fact that the initial condition for the position is only respected for the first value, and I don't understand how I can change this ? enter image description here


Solution

  • First of all, scipy.integrate.cumtrapz is left for backward compatibility, you should use the newer scipy.integrate.cumulative_trapezoid function. Secondly, if you read the documentation you'll see that initial is not an initial condition but simply a value prepended onto the array, which would normally be one element shorter than the original data.

    initial : scalar, optional

    If given, insert this value at the beginning of the returned result. 0 or None are the only values accepted. Default is None, which means res has one element less than y along the axis of integration.

    There is also a deprecation warning that, starting version 1.12.0, providing anything but 0 or None will result in a warning

    This shows the sizes of the result:

    from scipy.integrate import cumulative_trapezoid
    
    t = [0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7]
    a = [-9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8, -9.8]
    print(len(t), len(cumulative_trapezoid(a, t)), len(cumulative_trapezoid(a, t, initial=0)))  # 15, 14, 15
    

    To enforce your initial condition, you should set initial=0 and then add the initial condition to the result.

    v0 = 0.
    z0 = 5.
    v = cumulative_trapezoid(a, t, initial=0) + v0
    z = cumulative_trapezoid(v, t, initial=0) + z0