4

I was simulating the solar system (Sun, Earth and Moon). When I first started working on the project, I used the base units: meters for distance, seconds for time, and metres per second for velocity. Because I was dealing with the solar system, the numbers were pretty big, for example the distance between the Earth and Sun is 150·10⁹ m.

When I numerically integrated the system with scipy.solve_ivp, the results were completely wrong. Here is an example of Earth and Moon trajectories. Earth and Moon

But then I got a suggestion from a friend that I should use standardised units: astronomical unit (AU) for distance and years for time. And the simulation started working flawlessly!

My question is: Why is this a generally valid advice for problems such as mine? (Mind that this is not about my specific problem which was already solved, but rather why the solution worked.)

Wrzlprmft
  • 4,234
  • 1
  • 28
  • 54
Edward Henry Brenner
  • 157
  • 1
  • 9
  • Please correct and update your plot. Are the units actually meters, or is this in astronomical units [AU]? What is the time span used for integration, your particle is barely moving, the distance is `4e-11` of the radius, corresponding to `0.0002`seconds for the Earth orbit. – Lutz Lehmann Jul 24 '21 at 06:46
  • 1
    I am voting to reopen this question as it is not asking for debugging help and therefore does not require code. It is covered by all but the first main points of [What topics can I ask about here?](/help/on-topic) (“a software algorithm”, ”software tools commonly used by programmers”, “a practical, answerable problem that is unique to software development”). Also, as I hopefully demonstrated with my answer, this question can be reasonably answered without code and further clarification. – Wrzlprmft Mar 19 '22 at 10:27
  • 2
    @MisterMiyagi: I think the added last sentence should have made that already clear, but I boiled it down even further. – Wrzlprmft Mar 19 '22 at 10:41
  • Similar question on [scicomp.se]: [How does non-dimensionalization improve the behavior of ODE solvers?](https://scicomp.stackexchange.com/questions/41600/how-does-non-dimensionalization-improve-the-behavior-of-ode-solvers/41606#41606) – Wrzlprmft Jul 13 '22 at 08:04

1 Answers1

8

Most, if not all integration modules work best out of the box if:

  • your dynamical variables have the same order of magnitude;
  • that order of magnitude is 1;
  • the smallest time scale of your dynamics also has the order of magnitude 1.

This typically fails for astronomical simulations where the orders of magnitude vary and values as well as time scales are often large in typical units.

The reason for the above behaviour of integrators is that they use step-size adaption, i.e., the integration step is adjusted to keep the estimated error at a defined level. The step-size adaption in turn is governed by a lot of parameters like absolute tolerance, relative tolerance, minimum time step, etc. You can usually tweak these parameters, but if you don’t, there need to be some default values and these default values are chosen with the above setup in mind.

Digression

You might ask yourself: Can these parameters not be chosen more dynamically? As a developer and maintainer of an integration module, I would roughly expect that introducing such automatisms has the following consequences:

  • About twenty in a thousand users will not run into problems like yours.
  • About fifty a thousand users (including the above) miss an opportunity to learn rudimentary knowledge about how integrators work and reading documentations.
  • About one in thousand users will run into a horrible problem with the automatisms that is much more difficult to solve than the above.
  • I need to introduce new parameters governing the automatisms that are even harder to grasp for the average user.
  • I spend a lot of time in devising and implementing the automatisms.
Wrzlprmft
  • 4,234
  • 1
  • 28
  • 54
  • 1
    I suppose this answer is right, but it doesn't explain the problem from the beginning. From the beginning: ODE solvers, like the ones in solve_ivp, must select time steps which do a "reasonable job" of solving the ODEs. How do they do this? Well at each step the integrator estimates the local error introduced with a give step size ([lots of ways to do this](https://en.wikipedia.org/wiki/Adaptive_step_size)). They then select the step size which introduces LESS local error than is prescribed by the user... I've run out of space. but might be worth adjusting this answer to explain from the start – nicholaswogan Jul 25 '21 at 05:20
  • 2
    @nicholaswogan: *I've run out of space.* – Exactly. I can write my own textbook chapter on step-size adaption here, but then plenty of these already exist (and are much easier to find than the answer to this question). I have to presume that the asker/reader has some background knowledge or is able to obtain it with the right keyword. Anyway, I gave a one-sentence summary. – Wrzlprmft Jul 25 '21 at 07:44