Warning: Mega-nerd stuff ahead.
Way back, when I hadn’t derailed into programming yet (or did I?) I wrote a Object Oriented Visual C++ program to perform “identification” of mechanical systems. Identification is the process of determining certain parameters of a system, such as directly measurable weight and dimensions, but also of un-measurables like friction. You just can’t put a robot arm on some quit of measuring device (e.g. a scale) and read the friction from it.
One of the ways to perform identification is to conjure up a mathematical/mechanical description of your system. These consist of a complex set of equations, which may remind you of the Newtonian laws from your high school days. You know: F=m*a and the like. They are more complex because there are multiple dimensions and the differential equations might be coupled. The equations contain both the known and unknown parameters.
After that you will perform a test on your system, say the aforementioned arm, and record both the input signal (the voltages on the individual motors of the arm in time) and the resulting position, velocities and accelerating, i.e. the output signals.
Then you take the same set of input signals and feed these into your mathematical model. For this model the known parameters have been entered. The unknown parameters act as your “knobs“, that you can “turn” until the calculated output signal (from the model) has the best fit to the recorded output signals.
The turning of the knobs is done with a new set of differential equations, that is even more complex than the mathematical model that we started with. On average it contains around 10 coupled equations that are stiff. Stiff means that they are hard to solve with the simple mechanism for integration, such as Euler. You will need one hell of a differential equation solver to get that done.
So, what I did for my thesis on system identification was to write a multithreaded differential equation solver. I took an existing Fortran implementation called VODE, that implements a variable step ordinary differential equation (ODE) solver, and rewrote it to C++. I ran into some numerical boundaries that were very hard to track. Did you know that the x86 code for x*x might result in different results than x pow 2? Well I didn’t. I took me two weeks to find that in the millions of iterations that the solver performs. I did my fair share of power debugging there. The results of the calculations were well within numerical precision range for the algorithms used, but I wanted to get the results exactly the same as for the Fortran implementation. In the end they did.
To make a long story short, about two years ago I rewrote it to C#. The source code for the VariableOde integrator solution and the original Fortran code is included below. Everybody needs to do his bit of solving differential equations, don’t you? Go ahead and download it already :-). OK, but if you are a mathematician you might just appreciate it.integrator.zip (60.22 KB)