Problem when using non-integer step sizes in runge_kutta_fixed_step_size

[lrbusquets]

One year ago, working on question 1e of the Numerical Astrodynamics assignment 1, I reported the following issue in the course forum: when I used runge_kutta_fixed_step_size for different step sizes from an np. logspace array, I obtained an error figure that looked like this

whereas when I intriduced an np.round() to my step size array it turned to the expected

Although there seems to be no theoretical reason that prevents a Runge-Kutta integrator to work well with arbitrary step sizes, it looked as if some bug introduced significant extra error if a non-integer step was used. Dominic asked me to open an issue here, but it eventually went out of my mind I never did it.

Now, working on my MSc thesis, I encountered the same issue again. When doing a convergence plot just as we did in our first PnO assignment, I found that in order to fix it from this

to this

I just needed to apply an np.round() to the step size array I was working with. Note the scale of the vertical axis: in the former the errors are increased by a factor of several thousands from the nominal rounding error regime.

I think this might be a relevant issue for several reasons:

1. If there is some bug that causes this issue (e.g., some float precision reduced where it does not have to inside the runge kutta integrator), could it be also affecting the variable step size solutions?
2. In any ongoing research using TUDAT, could this error be affecting the accuracy of solutions believed to be valid?
3. If there is some theoretical limitation that prevents runge kutta integrators from working well when step sizes have too many decimals, it should be stated in the documentation!
4. For sure some students will find this issue again now that the Numerical Astrodynamics course is ongoing, and they might also end up spending hours before figuring out that a round() is all that needs to be applied.

[DominicDirkx]

Hi Lluc,

Thanks for posting, and reminding me of this. But, it should still be addressed in the underlying code. First, to answer your few questions:

1. I suppose it could, but since the variable step integrators will in general take steps that are not integers, we should have seen this issue arising there already, but we don't
2. In theory, yes. But, since this issue will always show up when doing a test of validity of integrator settings, people will catch it if they are doing their work properly
3. No, I suspect this is somehow related to floating point issues in the time representation, but I'm not yet sure where/how
4. We've slightly modified the Numerical Astrodynamics assignment, so that this issue no longer pops up.

Could you share a link to a repo with the code that generates your results? That will be helpful to debugging it,

Dominic

[lrbusquets]

Hi Dominic,

Thanks for your reply. If you believe the bug is in my code, please do not spend your time on it. In my view, however, I don't understand how the cause can be anywhere else than in the Tudat implementation of the integrator. Here's what I think:

The case I reported last year (assignment 3, not 1 as I mistakenly stated) was based on the comparison of a numerically propagated Keplerian orbit with respect to its corresponding analytical solution. What I came across this time, in contrast, is based on a numerical convergence analysis using perturbed dynamics, in which each error value is obtained from the difference between the trajectory propagated with the corresponding time step and the trajectory propagated with a timestep twice as small. In other words, these are two completely different strategies that suggest that something is going on when the step sizes used have several decimal figures, and the only thing the two strategies have in common is the usage of the runge_kutta_fixed_step_size function.

Code for the Ganymede plot: file "integrator_analysis_Q1e" in https://github.com/lrbusquets/NumericalAstro_2023_lrbusquets/tree/main/Assignment3
Code for my thesis: file "planets_and_asteroid_integrator_analysis.py" in https://github.com/lrbusquets/master_thesis/blob/master/planets_and_asteroid_integrator_analysis.py

My first thought was that this could be reasonably explained if there were an inconsistency in the float precision used in the C++ scripts that conduct the integration. For instance, if the evaluation of each state

is conducted maintaining the double precision of the step size h coming from python (where it is a float64), but in the intermediate evaluations

the state derivative function f reduces the number of decimal points of the step size h (for instance, because inside the C++ script it is declared as a float instead of a double), then the intermediate derivatives are not evaluated at the exact instants they should, and therefore it can be reasonable that an extra numerical error accumulates over time.

If this were the case, I could imagine that variable step size integrators do not show any misbehaviour because the C++ function that determines the next step size to be used stores it as a float variable in the first place. But I know I should look at the scripts instead of speculating.

I might do some more tests to check whether or not this explanation is supported and I will post here in case I find anything new. Of course I might well be absolutely wrong, so don't hesitate to mention any other possible explanation you have or if you want me tro do any test in particular.