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

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)

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.

• Thank you very much, i will investigate these points ASAP. However, i am not authorized to implement other integrations functions that what is already inside the simulation framework. (basically, i have access to forward euler, RK2 and RK4). Do you think that, by switching to Euler if i am able to distinguish an error scaling like h, that could be enough to say that my reference solution is already using a step size that's too small? – J.M Oct 7 '20 at 14:16
• Also, the program only uses double precision floating point variables. Correct me if i'm wrong but i belive this prevents rounding errors to occur as soon as h ~10-3 – J.M Oct 7 '20 at 14:22
• No, I think not that RK2 can be useful in the estimation of the RK4 error. You still have the use of Richardson extrapolation, that is, comparing the result of step $h$ with the result of step $2h$. It is also worth testing if the full simulation actually does have error order 4, the graphs look more like order 2. In that regard comparing with RK2 should be useful. (I do not understand the restriction, while it is a qualitative jump to implement a method with adaptive step size, implementing other explicit RK methods with fixed step size is practically just switching the Butcher tableau.) – Lutz Lehmann Oct 7 '20 at 14:23
• I'll go for that then. In terms of why i cannot implement any higher order method is because the framework i am using is a rather old tool, and is quite cluttered (implementing a new integration method would represent a few weeks of work, to say the least), plus, i dont have the "right" to do it. – J.M Oct 7 '20 at 14:29
• As said, the masking of floating point errors by the truncation error of the method occurs for $h>10^{-3}$ only for ODE that have scales $\sim 1$ in values and derivatives of the ODE function. If one rescales the ODE then the scale of the time can also change, so that the optimal $h$ can also have a drastically different value. You can, as example, do a solar system simulation in the km-sec setup or in the AU-day setup, physically it is the same, but the numerical values of corresponding step sizes have drastically different scales. – Lutz Lehmann Oct 7 '20 at 14:29