Search code examples
pythonsimulationnumerical-methodsnumerical-integration

Why should I use normalised units in numerical integration?


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.)


Solution

  • 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.