Apparent integer counter overflow bug detected in the use of the integrate_const method

[CFMaguire]

I am looking at n-body simulations of the Trappist-1 seven planet system orbiting a dwarf star. For the software I adapted the sophisticated solar_system.cpp example in the tutorials which uses the integrate_const method. Since the orbital period of the innermost planet is 1.5 days such simulations require time steps on the order of 0.1 day. So simulations on the million year time scale will require billions of calls in the integration solution. When I first tried such a long simulation all went well until about 587 thousand years at which the time reported by the simulation switched to an equally large negative value. Thereafter the integration observer reported decreasingly negative time evaluations as would be expected if the previous time value was getting a positive increment. This behavior appeared consistent with an integer counter overflowing the INT_MAX value (2147483647, 2.15 billion) after which it changes to the most negative value. In order to reproduce this effect I adapted the harmonic oscillator program to use only the integrate_const method for ever smaller time steps. That program version is the attached harmonic_oscillator_Debug.cpp file whose output file debugOutput13July2022. txt for a 3 nanosecond time step reproduces the behavior. In fact this is an infinite loop bug since the maximum time requested will never be achieved and the program cycles between positive and negative time values. A user workaround is available. The requested duration of the simulation can be split into successive segments such that the overflow integer counter is never encountered. The end of one segment becomes the initial state vector for the next segment. The files harmonic_oscillator_Fix.cpp and its output file fixOutput13July2022.txt show the use of this workaround. When I applied the workaround to my own n-body simulation code it appeared to work over a 1.5 million year calculation. The final orbits were stable with the expected orbital parameters. As a new user for the ODEINT system my apologies in advance if this is a spurious report.

harmonic_oscillator_Debug.cpp.gz

debugOutput13July2022.txt.gz

harmonic_oscillator_Fix.cpp.gz

fixOutput13July2022.txt.gz