Search code examples
pythonnumerical-methodsdifferential-equationsnumericaldifferentiation

Calculating acceleration by numerical differentiation in Python


From stylus movement measurements I got two measurements: the time and the distance. So I have two arrays, first an array of timestamps (in milliseconds) and then an array of the same size of distance measurements. For example the two arrays could look like:

distance = [1.4142,1.0000,1,0,1.0000,1.0000,0,0,1.0000,1.0000,0,1.0000,2.0000,2.2361,0,3.0000,3.6056,3.1623,3.1623,0,3.6056,3.1623,3.1623,0,1.4142,2.2361,1.0000,0,0]

timestamps = [1563203.5,1563208,1563210.5,1563213.5,1563218.5,1563223.5,1563226.5,1563229,1563233.5,1563238.5,1563242.5,1563245,1563248.5,1563253.5,1563258,1563260.5,1563263.5,1563268.5,1563273.5,1563276.5,1563279,1563283.5,1563288.5,1563292.5,1563295,1563298.5,1563303.5,1563307,1563317.5]

I think the first derivative gives me the speed and the second derivative gives me the acceleration.

I'm interested in calculating the acceleration using numerical differentiation. How can this be done in python?


Solution

  • The first derivative can be calculated as (f(x+h) - f(x-h)) / (2h). This gives an estimated error on the order of h^2.

    The second derivative can be calculated as (f(x+h) - 2f(x) + f(x-h)) / h^2. This also gives an estimated error on the order of h^2.

    Note: you can use (f(x+h) - f(x))/h or (f(x) - f(x-h))/h to estimate the derivative, but these give estimated errors on the order of h, which is larger than h^2.

    I haven't looked at your data in detail, but these formulas don't quite work out if you don't have a constant h interval between your data points. In that case I would simply calculate the acceleration in two steps:

    v(x_i) = (p(x_(i+1) - p(x_(i-1))) / (x_(i+1) - x_(i-1))
    a(x_i) = (v(v_(i+1) - v(x_(i-1))) / (x_(i+1) - x_(i-1))
    

    You basically just feed in your data for the first derivative back into your derivation algorithm and you get the second derivative.