Error control and sequence acceleration at the same time

Question

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.

EDIT:

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.

Answer 1

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.

Answer 2

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.