Choosing an appropriate time step for a discrete & continuous dynamics simulation

Question

I have inherited of a flight dynamics simulation in C++ which represents a small drone with it's autopilot, actuator dynamics and a solid state IMU.

Hence, it is composed of a few models, some continuous (flight dynamics & actuators), integrated with a runge kutta 4 scheme, some discrete (the autopilot and IMU). I have full control over the physics timestep. The autopilot is supposed to run at 500hz, the solid state IMU at 2000Hz.

I am to find a "correct" time step value in order to minimize errors while maintaining a reasonnable conputationnal time.

I tried plotting the difference in mechanical energy between the highest frequency i ran the model at (16Khz) and the others simulations (500hz, 1000hz, 2000hz, 4000hz, 8000hz) I have selected these frequencies in order for the simulation to step precisely on the "activations" of the autopilot and IMU.

I did the same on a L2 norm composed of the flight dynamics state variables (speed, position, rotationnal speed and euler angles).

The results were very different from what i would have expected:

While the errors seem acceptable (errMax ~ 0.3% ), i do not understand why there is such an increase in error around the 2000 / 4000hz point. Also, the errors dont seem to decrease with the time step. This leaves me quite puzzled on the relevance of my approach.

Would anyone know any reason that could cause an increased error around a specific time step in a simulation that mixes continuous and discrete state models?

Also, is there any kind of relevant physical analysis that would allow one to get a norm that accurately represents a simulation state (in order to study the convergence of said simulation)

Answer 1

This is speculation, as I do not know the Lipschitz constant or the derivative scales of your simulation. Also, there might be some resonance effect in the interplay of the discrete and continuous parts. But what I would first draw attention to is that the error of RK4 (and any other method) has a V shape in a loglog plot vs. the step size. This is the result from two contributions to the error, one of size $\mu/h$ where $\mu$ represents the size of floating point noise per integration step and $1/h$ the number of integration steps. The other contribution is from the method error and is of size $Ch^4$. They both balance at around $h=\sqrt[5]\mu$, which for floating point is around $10^{-3}$.

This is for test models of small dimension with appreciable scales in values and derivatives. In general it should give the impression that RK4, like any higher order method, has its best behavior at relatively large step sizes. It could be that your reference solution is already far into the side of too small step sizes. Then the accumulation of floating point errors gives it a quasi-random distortion that should disqualify it from being a reference solution.

What you can do to test this hypothesis is to compare the lower frequency solutions to each other, double and half step size, using Richardson extrapolation, and check if the error scales like $h^4$ in most components of the error vector.

Alternatively you can compare the RK4 solution to a fixed-step "abuse" of a higher order method, for instance, the 5th order stage of Fehlberg or Dormand-Prince, or something of even higher order. For the lower frequencies this should give a valid reference solution.

I implemented this fixed-step for DoPri5 to test the order of the method in https://stackoverflow.com/a/54502790/3088138, the stages look more cluttered, but the principle remains the same.