1
$\begingroup$

Suppose you want to approximately solve a system of ODEs, using some numerical method (Euler, RK, BDF, whatever):

$\frac{du}{dt} = f(u)$

To do this you need to select time steps which solve the ODEs to some tolerable precision. One approach, for selecting timesteps, is to look at timescales of change of each variable, and select a timestep related to the minimum timescale:

$\Delta t = \alpha\min\left(\frac{u}{f(u)}\right)$

Here, $\alpha$ is some constant and $\Delta t$ is the timestep. I’ve seen many people in astronomy and planetary science use this approach.

Question: is this a "bad" approach?

It seems not the best to me. It seems much better to use things like embedded runge-kutta for stepsize control, which selects the biggest timestep which admits less than some specified error (rtol or atol).

$\endgroup$
2
  • $\begingroup$ The right approach is to go in segments, first using one time step and a then a smaller time step (say 1/2 size), and verify the level of agreement; if above the tolerance then repeat reducing the time step. I believe Numerical Recipes has some illustrative implementation of this approach. A mature ODE package like SUNDIALS would have this implemented in a pretty robust way. $\endgroup$ Jul 5 at 22:50
  • $\begingroup$ For a direct comparison of CFL-based step size selection (similar to what is discussed here) versus an error-estimation based approach, see some of the experiments in this paper: arxiv.org/abs/2104.06836 $\endgroup$ Jul 7 at 10:18
1
$\begingroup$

As stated there's no re-rejection mechanism, i.e. ability to decrease the stepsize after a step has potentially failed. This is required for implicit methods which have Newton steps since there's a chance the $\Delta t$ is large enough that the (quasi-)Newton is unstable, in which case it needs to pullback on time. This instability can sometimes be seen via convergence rates, but in other cases it's caught by a posteriori error estimators which are absent here. Additionally, it's relying on a linear error estimate, while a posteriori estimators are inherently nonlinear using multiple values from the step, so this is likely to not be correct on very nonlinear problems. These facts together mean that this method likely is only applicable to non-stiff or at most semi-stiff (PDE) problems.

However, one of the bigger factors is that it's just suboptimal. Embedded error estimates are able to get a nonlinear error estimate with 0 extra $f$ calculations. This gives a linear error estimate with 1 extra $f$ calculation, meaning its likely more expensive while being less robust than a method which is developed to have an internal error estimator. For this reason I would be surprised to see it in any non-stiff ODE solver that is used in practice. That said, such a scheme can be a quick (albeit inefficient) way to change a method without embedded error estimates into something that has a decent adaptivity algorithm, if $\alpha$ is chosen sufficiently small for stability and one is willing to calculate $f$ quite a few more times.

$\endgroup$
2
  • $\begingroup$ I've read your article on why not to use forward Euler. It would be useful to have an analogous article which explains why using the "time scale approach" I described isn't the best. Maybe the article would have some benchmarks. Most everyone around me in the world of planetary science uses home-brewed ODE solvers with the timescale approach. I mostly worry that they are not actually solving the ODEs they are trying to solve to any reasonable precision. $\endgroup$ Jul 7 at 16:58
  • $\begingroup$ Yeah, it wouldn't be hard to show a few examples where it actually gets unbounded error. I'll add it to the queue. $\endgroup$ Jul 10 at 19:15
1
$\begingroup$

I do not think this is a bad approach, but it is not a very precise way to select timesteps either.

Admittedly, I have not come across this sort of timestep heuristic before, but looking at a linear test problem provides some insight as to why this is reasonable. For $y' = \lambda y$, the conditions becomes $\Delta t = \frac{\alpha}{\lambda}$. This looks exactly like a timestep restriction from the linear stability analysis, where $\alpha$ would be related to the size of the stability region. For a linear PDEs that has been discrtized in space, the condition looks an aweful lot like a CFL condition (if we generalize your condition to vector-valued $f$), as you can get expressions of the form $\Delta t = \alpha \, O(\Delta x^k)$. So to me, this heuristic seems to be guided by (linear) stabililty and serves as a reasonable upper bound on the timestep.

Conversely, more standard and popular timestep controllers based on embedded methods are guided by an approximation of the local truncation error committed at each step. Stability is not really a direct consideration here. The standard approach also gives more control via tolerances and can change from step to step. Sometimes embedded methods are not readily available, though, so this does give some credence to the "timescales of change" approach.

Edit: On second thought, I have come across something roughly resembling this with singular perturbation problems: $y' = \varepsilon^{-1} f(y)$, $0 < \varepsilon \ll 1$. During an initial transient phase, timesteps need to be on the order of $O(1/\varepsilon)$ for stability and accuracy. Past this transient phases, the solution is smooth and stiff. For some methods, e.g. B-convergent Runge-Kutta, choosing the timestep purely based on the $\varepsilon$ timescale would force you to take much smaller steps than required.

Edit 2: Just to be clear, I agree with your assessment. Standard error controllers are superior in most ways, and I would only consider "timescales of change" if I had no way to estimate the local truncation error of a step.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.