Significantly higher roundtrip errors with new IAS15 adaptive timestepping scheme

[rahil-makadia]

I am seeing significantly higher roundtrip integration errors if I use the new adaptive stepping scheme (mode 2, see attached figures). The figures are generated from the Figure 2 section of the publication_figures.ipynb notebook by just setting the sim.ri_ias15.adaptive_mode to either 1 or 2. I have also implemented the new timestepping scheme in my small-body OD/propagation library and am seeing worse behavior with long-term integration of (65803) Didymos (I suspect this is just because (3666) Holman is a main-belt asteroid vs. Didymos being an NEA).

Is this an expected consequence of the new adaptive timestep algorithm? If so, @hannorein do you have recommendations for when to use each scheme? Because I am seeing great stability for integrations with close approaches (e.g., Apophis) when I slightly modify the minimum timestep with the new algorithm (the previous one would oscillate significantly in those scenarios). So, my current intuition dictates that if I am expecting to propagate deep encounters/impacts, I should use the new timestep algorithm, but not in other cases such as long integrations (~100 years). Although I'm not sure what to prefer if I don't have this information beforehand...

[rahil-makadia]

Correction, I am seeing similar roundtrip errors with Didymos, not worse. I was misinterpreting my results but the issue still stands

[hannorein]

The old timestep algorithm might stall for very close encounters as we discuss in Pham et al. Both time stepping algorithms result in the same timestep for circular orbits (that's how we've normalized the new one). However, for eccentric or hyperbolic orbits, the new algorithm gives a timestep ~2.5 times larger. I suspect this explains the difference in accuracy.

Try the following: use the new algorithm, but reduce the epsilon value by a factor of 1000. The simulations should be ~2-3x slower, the accuracy should match what you obtained with the previous algorithm, and there should be no more (fewer) issues with stalling simulations for very close encounters.

(FWIW, there is an open pull request #102 to just use the old algorithm for now as a default)

[rahil-makadia]

Reducing the tolerance did make an improvement without a significant performance hit. Can you explain a bit more what the intuition was behind reducing it by a factor of 1000 just in case I run into something similar in the future? Seems like the new algorithm is enabling much lower epsilon values than the old one...

I'll keep in mind that the open pull request might get merged when I am comparing the GRSS propagator with ASSIST in the future!

[hannorein]

These are all great questions!

Yes! You should be able to use much lower epsilon values without getting stuck.

For the intuition about the factor 1000, see Equation 17. The timestep is proportional to $\epsilon^{1/7}$. Thus reducing epsilon by 1000 reduces the timestep by about 2.7.

[rahil-makadia]

Ah I see, so because we know the factor by which the new algorithm is increasing timesteps, we can figure out what the new epsilon should be with the updated algorithm to match the old one. Makes sense, thanks! Feel free to close the issue

[hannorein]

Correct. We could have chosen the new default timestep to always be smaller than the old one, but that would be way too pessimistic in most situations. This one here is an exception because very high precision is required during a close encounter.