In a posteriori error control for solving ODEs, one typically computes two different approximate solutions, one of which being "more accurate" and one of which being "less accurate". If $y_q^{n+1}$ is the more accurate approximation to the solution at the next time step $y^{n+1}$, and $y_p^{n+1}$ is the less accurate approximation, one can estimate the error for the less accurate approximation, e.g., $$\text{error}(y_p^{n+1}) \approx y_q^{n+1} - y_p^{n+1}.$$

This error estimate can then be used to determine if, for example, the time-step size should change according to some prescribed accuracy requirements by the user.

Q: Could one use "sequence acceleration" techniques and "error control" at the same time using just $y_q^{n+1}$ and $y_p^{n+1}$? Or would this be "having my cake and eating it"?

Some cases I had in mind:

  • $y_q^{n+1}$ could come from a formally higher-order approximation and $y_p^{n+1}$ from a lower-order approximation as in embedded Runge-Kutta schemes.
  • $y_q^{n+1}$ could come using a small time-step size when compared to $y_p^{n+1}$. Should I do "Richardson Extrapolation" in this case?
  • $y_q^{n+1}$ and $y_p^{n+1}$ could be "Milne pairs" which have the same order of accuracy, but provide a "Milne error estimate" of the form $\text{error}(y_p^{n_1}) \approx C (y_q^{n+1} - y_p^{n+1})$ where $C$ is known. Should I take some combination of $y_q^{n+1}$ and $y_p^{n+1}$ to obtain a higher-order estimate?
  • $y_q^{n+1}$ could come from an implicit "corrector" step (perhaps of higher-order) and $y_p^{n+1}$ could come from an explicit "predictor" step in a predictor-corrector scheme.

Q: In all these cases, what would be a smart choice for $y^{n+1}$? Should I set $y^{n+1} := y_p^{n+1} + \text{error}(y_p^{n+1})$? Should I set $y^{n+1}:= y_q^{n+1}$?. Should I do something else entirely?

I have also seen that one can choose $y_q^{n+1}$ and $y_p^{n+1}$ to have the same leading error terms but with opposite signs, so that their average is higher order. I am not sure how applicable this technique is in general.


What I mean by "sequence acceleration" is: obtaining local solutions (at the next time step) $y^{n+1}$ which are more accurate than the solution $y^{n+1}_p$ for which we have an error estimate of.


2 Answers 2


You should state clearly what you mean by sequence acceleration. But if I understand you correctly, what you're asking about is exactly what extrapolation codes do. A sequence of low-order approximations of the new solution are computed and then combined (via extrapolation) to produce high-order approximations. The highest-order and second-highest-order approximations are also used to produce an error estimate. The canonical reference for this is the 2-volume series by Hairer, Norsett, & Wanner. If you don't have the books and are looking for a straightforward explanation, you might look at this paper of mine. A canonical implementation is Hairer's ODEX code.

In the early days of error estimation, people used the lower-order estimate as starting point for the next step, but since at least the 1980's everyone uses the higher-order estimate; this choice is sometimes referred to as local extrapolation.


What you're talking about is local extrapolation. Local extrapolation is the idea of getting an error estimate between two methods and then continuing (accepting the new value) from the higher order method. Dormand-Prince is an influential paper which created a locally extrapolating 5th order method where the 5th order method minimizes its truncation error coefficients (A family of embedded Runge-Kutta methods). This Dormand-Prince method was found to be modestly more efficient (for free) than the original adaptive Fehlberg algorithm which uses a 4/5 pair but optimizes the coefficients on the 4th order method and then uses that to advance. This then became the dopri5 method of Hairer, and it went on to be the standard method in most places for a long time (for example, it's ode45).

The reason why local extrapolation was used here is because it's stable. The Runge-Kutta methods were built in such a way that they had similar stability regions, and so the error estimator was stable when continuing with either one. I think stability of adaptive algorithms is covered in Hairer II, where there's also details about using PI-controlled stepsizes to further increase the stability with adaptivity. This isn't true for every algorithm though. But Shampine studied the effects on the most popular algorithms and local extrapolation was shown to either not be harmful or in some cases it increases the stability of the adaptive algorithm, so local extrapolation is used in most modern codes for this reason.

  • $\begingroup$ When using any sort of implicit method, stability is not really a concern. In that case, one should always do local extrapolation by continuing from the higher order method, right? $\endgroup$ Jul 1, 2018 at 20:34
  • $\begingroup$ I believe the paper by Shampine, which I cannot access at this moment, gives some rule of thumb for when one should perform local extrapolation. If I remember correctly, it involved some coefficient $\alpha \approx 1/15$. Is this "usually" still considered? Also, I believe the type of local extrapolation he considered added the estimated error to the lower-order solution, rather than continuing from the higher order method. Why do this instead of just continuing from the higher-order solution? $\endgroup$ Jul 1, 2018 at 21:40
  • 1
    $\begingroup$ "When using any sort of implicit method, stability is not really a concern". No, that's not true. For example, with the standard Rosenbrock W 2/3 method, the third order method is not L-stable with the 2nd order one is, so the thrid order is only used for error control. In many cases you can have a smaller stability in the error control. Another case of this is the canonical TRBDF2. $\endgroup$ Jul 1, 2018 at 23:08
  • 1
    $\begingroup$ Stability in for implicit methods isn't that easy anyways. It can also be do to the quality of step predictors. There's A-stability which is the weakest stability, but then L-stability and B-stability. Error is not that easy either. Order only tells you what happens asymptotically, while the size of truncation error coefficients matters in real integrations and can cause the lower order method to have less error. $\endgroup$ Jul 1, 2018 at 23:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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