7

I have a set of coupled ODEs that I'm solving numerically. The independent variable is time and runs from values of $10^{15}$ to $10^{17}$ in units of seconds. The state variables in their usual physical units have large numbers like $10^{10}$ and $10^{50}$. As these state variables evolve with time, their values systematically increase by ten to twenty orders of magnitude.

I am able to solve my equations with SciPy's solve_ivp (either RK23 or the implicit BDF method) using the default relative error tolerance (rtol=1e-3) and an arbitrary absolute tolerance atol=1. Amazingly, this is with all variables in their default physical units and huge orders of magnitude. However because of the huge scales of my variables, I am wondering if it is worth investing some time in Nondimensionalization.

What advantages would nondimensionalization offer? Would it make the ODE solver behave better in terms of finishing the integration faster but still being precise (i.e., giving the same solution I have now in the dimensionalized version)? Would it provide a natural choice for the atol (and rtol?) parameter which currently I don't really know how to set?

Wrzlprmft
  • 2,032
  • 12
  • 32
quantumflash
  • 193
  • 5

1 Answers1

6

Background

Let’s have a brief look at what parameters an integrator typically has and what they are used for. All of them govern step-size adaption:

  • Relative tolerance (rtol) is usually the most important parameter governing accuracy. After an integration step, the integrator estimates the magnitude of the error it makes and if it exceeds the relative tolerance, a smaller step is taken instead. Values between $10^{-5}$ and $10^{-3}$ are typical. If you have any way to estimate which relative error would affect your analysis, a few orders of magnitude below this is a good choice. However, usually it’s just how accurate you want your integration to be. Since this value is relative, nondimensionalisation should not affect it.

  • Absolute tolerance (atol) is like relative tolerance, just going by absolute values. It is conventionally cumulative with absolute tolerance, i.e., the error needs to be below atol+rtol*state. Unless you set it to a high value or atol to zero, it becomes only relevant if a dynamical variable crosses zero, as relative tolerance cannot cover this well.

  • Most integrators have other parameters like the maximum and minimum step size, the initial step size and further parameters governing the step-size adaption.

Why nondimensionalise?

For a normal integrator, the default values of these parameters are made for a setup where all your dynamical variables have the same order of magnitude, that order of magnitude is 1, and the smallest time scale of your dynamics also has an order of magnitude of 1.

If this does not apply, you have to consider all integration parameters (not just the tolerances) as to whether they need to be adjusted. For many integrators, you also have the problem that you cannot set the absolute tolerance per component. Unless you know what you are doing and can set all the parameters accordingly, it’s easier and safer to nondimensionalise.

Possible consequences of badly set integration parameters include a failing integration, bogus results, and an overly slow integration.

Your specific case

Some thoughts on your specific scenario:

  • Your time scales are very high. I would have to look into the details of the integrators you used, but it might very well be that you barely made use of the adaptive integration, and instead progressed with the default maximum time step (which in turn is still very small compared to the time scale of your dynamics) and thus wasted computation time.

  • Your absolute error was probably be too low. If your dynamical variables never become zero, it is not surprising if this had no effect. You may also just have been lucky. In general, if you have zero crossings and you have dynamical variables of different orders of magnitude, you cannot have an absolute error that is appropriate for all of them (see above).

  • Your relative error is okay.

  • If your dynamics varies on several orders of magnitude, consider performing the entire integration in the logarithmic domain.

Digression: Why can’t the integrator do this?

Adapted from this answer of mine

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 be spared from problems arising from overly large or small dynamical variables or time scales.
  • About fifty in a thousand users (including the above) miss an opportunity to learn some fundamentals about how integrators work and reading documentations.
  • About one in thousand users will run into a horrible problem with the new 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
  • 2,032
  • 12
  • 32
  • I'm kind of late but how do I decide whether to integrate in the logarithmic domain for both the state variables AND time (dlogX/dlogt) or just the state variables (dlogX/dt)? Integrating in the logarithmic domain also seems to be sensitive to the initial conditions because replacing dX/dt = XdlogX / (t dlogt) means multiplying the RHS by t/X. – quantumflash Jan 29 '23 at 01:12
  • Also, what do rtol and atol intuitively mean for dlogX/dlogt with logarithmic state variables, and should atol scale with the initial conditions somehow? This is especially confusing if the initial condition is X=1, so dlog(X)=0. What would a good atol be in this case...? – quantumflash Jan 29 '23 at 01:40
  • Logarithmic time is something that only makes sense if you have a special point in time such as the big bang (in cosmology) or the time of a mutation (in population genetics) and if your timescales span many orders of magnitude (which doesn’t apply to you). Also, you obviously cannot do this, if you want to start at $t=0$. If that is given, I don’t see how this would introduce sensitivity to initial conditions as the dynamics itself should not change. – Wrzlprmft Jan 29 '23 at 10:36
  • Also, what do rtol and atol intuitively mean for dlogX/dlogt with logarithmic state variables, and should atol scale with the initial conditions somehow? – When working logarithmically, you usually only want to work with atol as it effectively gives you relative tolerance. Any rtol would mean that you accept your actual relative error to be larger when the order of magnitude of your variable is larger, which usually makes little sense. – Wrzlprmft Jan 29 '23 at 10:40
  • This is especially confusing if the initial condition is X=1, so dlog(X)=0. What would a good atol be in this case...? – I don’t see a problem with atol here, as zero is not a special point for it. Interpreting rtol may be weird, but see my above comment for that. – Wrzlprmft Jan 29 '23 at 10:41
  • Thank you so much! I'm not starting at t=0 but in my units of seconds, I start a bit after the big bang at t~10^15 sec and go until 10^17 sec, so the integration time spans two orders of magnitude and I'm working with large numbers for the time variable. I also just generally like the idea of having both the numerator and denominator of the derivative be log so the whole derivative is dimensionless rather than just dlogX/dt. In the case of dlogX/dlogt, how does atol give you a relative tolerance? Is it because dlogX=dX/X and dlogt=dt/t, i.e., fractional changes in X and t? – quantumflash Jan 29 '23 at 15:49
  • So are you saying I should set rtol=0 and atol = the maximum fractional change in X (and t???) that I am willing to accept? I don't understand why we wouldn't want our relative error to indeed be larger as the order of magnitude of X increases? For example, say at t=10^17, X is 10^50, that means logX is 50 and logt=17. What would atol=1e-6 and rtol=1e-6 mean here? And similarly at an earlier t=10^15, if X is 10^10, logX is 10 and logt=15 -- what would the same atol=1e-6 and rtol=1e-6 mean here? It feels wrong to set rtol=0. – quantumflash Jan 29 '23 at 15:55
  • how does atol give you a relative tolerance. – In brief, it is because multiplications become additions in the logarithmic domain. Mathematically, suppose you want to limit your absolute error of the variable $x$ to $ε$. Then your margin of allowed results is $x·(1±ε)$. In the logarithmic domain, this becomes $$\log(x·(1±ε)) = \log(x) + \log(1±ε) ≈ \log(x) ± ε.$$ – Wrzlprmft Jan 30 '23 at 14:37