High Order Numerical Differentiation in C++

Markus Buchholz
3 min readJan 19, 2023


The following article depicts the simple numerical methods to compute derivatives in C++. I am going to focus on robotic applications and simply describe how to compute velocity, acceleration, and jerk while robots move along the path with certain time constraints.
The computational methods for obtaining the high-order derivatives are not limited to robotics applications only, therefore displayed below methods are of course applicable to any other domains.
We are going to compute the first, second, and third derivatives of the path since the velocity, acceleration, and jerk are often evaluated while the robot moves (mobile robot or manipulator).

The source code you can find on my GitHub.

Please note, I limit the article to depicting only the final equations without focusing on the approach step to derive final formulas. I do not focus also on numerical errors.

Numerical Differentiation

Numerical differentiation is a collection of methods to estimate the value
of the derivative of the function at a given point. Wherever the form of the function is not known analytically or its form is very complicated, and for functions that values are known only for a discrete set of arguments, differentiation numeric is applicable.

Newton’s method is the simplest way to calculate the differential of a function f(x).
The definition of the derivative of a function is used here, which is approximated by the quotient differential:

Assuming that when calculating f (x) the values for the previous arguments are known, it is worth considering approximating the derivative with a difference quotient in the form (backward differentiation):

In the following article, we have to apply the second and third numerical differential since we need (besides the velocity), acceleration and jerk.

The second derivative (acceleration) can be computed as follows,

The third derivative (jerk) can be computed as follows,

Using the specified above equations I applied them to compute discussed values of robot motion.

First I used a cubic polynomial which has four coefficients. Therefore, it can satisfy the position and velocity constraints at the initial and final points. I do not compute (plot) jerk since the acceleration for cubic path jumps at boundaries and affects jerk to be infinitive. (note plots are scaled)

motion curves (by author)

In order to estimate the jerk value I used the following polynomial (seven—order) which guarantees that the jerk is zero at the start and stop),

motion curves (by author)

Plotting requires incorporating the header file which has to be in the same folder as your cpp (a file you can clone from my repository).
Your program can be compiled as follows,

g++ my_prog.cpp -o my_prog -I/usr/include/python3.8 -lpython3.8//


//folder tree

├── my_prog
├── my_prog.cpp
├── matplotlibcpp.h

Thank you for reading.