A new mathematical study asks a deceptively simple question about that machinery. If a solver estimates its own error rigorously, using residual-based a posteriori error estimates, can the resulting feedback loop that adjusts the time step remain stable? The answer, published in Computational Science and Engineering by Hendrik Ranocha and Jan Giesselmann, is nuanced. More advanced PI and PID controllers, however, can be designed to stay stable, and the residual-based estimators behave comparably to the embedded error estimates used in production codes while adding a bonus the classical approach cannot offer: a provable upper bound on the global error.
The distinction between the two error-estimation strategies matters. Embedded Runge-Kutta pairs, such as the well-known schemes of Bogacki and Shampine or Dormand and Prince, compute a second, cheaper solution at each step and take the difference between the two as a local error estimate. This is extremely efficient for explicit schemes and is the workhorse of adaptive ODE integration, but it provides no rigorous certificate that the accumulated numerical error actually respects the requested tolerance. Residual-based estimators take a different route. They reconstruct a sufficiently regular function from the discrete solution, measure how badly that reconstruction fails to satisfy the differential equation, and then use a stability argument, often based on a one-sided Lipschitz condition and Gronwall’s lemma, to convert the residual into a genuine bound on the error. The price is computational: residuals are more expensive to evaluate. The reward is reliability, because the resulting estimate is a proven upper bound rather than a heuristic difference.
Ranocha and Giesselmann focused on what happens when these residual-based estimates are plugged into the feedback loop that selects the time step. The authors analyzed this using a scalar linear test problem, deriving the Jacobian of the dynamical system that couples the solution and the step size, and then computing its spectral radius. If that spectral radius exceeds one, the feedback loop amplifies perturbations instead of damping them, and the controller is unstable.
The analysis covered explicit Runge-Kutta methods from first to third order. For the explicit Euler method with a linear reconstruction, and for all explicit second-order two-stage methods with a left-biased quadratic Hermite reconstruction, the spectral radius of the Jacobian under a simple I controller exceeds unity in parts of the complex plane. The same holds for third-order three-stage methods with both left-biased and central cubic Hermite reconstructions. In other words, the classical I controller is not step size control stable for any of these residual-based estimators when the step size sits on the boundary of the stability region.
The picture changes with more sophisticated controllers. For second-order methods, a PI controller with parameters beta equal to 0.6 and minus 0.2 stabilizes the residual-based feedback loop, while the same controller leaves the embedded-method version unstable. For third-order methods, the PI controller with those parameters, recommended in earlier work by Ranocha and colleagues, stabilizes both the residual-based and embedded approaches, whether the residual is measured in the L1 norm with left-biased interpolation or in the L2 norm with central cubic Hermite interpolation. The L2 estimator with central interpolation showed somewhat better damping near the phase angle where the I controller is most badly behaved, though it still did not achieve stability under the I controller. The authors also derived an analytic expression for the L2 residual with central cubic Hermite interpolation, something they could not do for the L1 version, which required numerical quadrature.
Numerical experiments reinforced the theory across a range of problems. On a nonlinear four-dimensional ODE test problem from Krogh, integrated with the third-order Bogacki-Shampine method, the simple I controller rejected a significant number of steps when the problem’s eigenvalues rotated into the unstable region around a particular phase angle, while the PI controller rejected at most one or two steps for all parameter values. On the Euler equations of a rigid body, neither approach had step rejection problems, and both showed the expected convergence behavior: the embedded method achieved tolerance proportionality, meaning the global error scales roughly linearly with the tolerance, while the residual-based approach with error-per-step control produced a global error scaling as tau to the power three quarters, which can be rescaled to recover proportionality. Efficiency, measured as error for a fixed number of steps, was the same for both.
The authors then extended the analysis to partial differential equations discretized by the method of lines. For two-dimensional linear advection discretized with a discontinuous Galerkin spectral element method, both estimators behaved nearly identically, converging to the optimal step size that a manually tuned Courant-Friedrichs-Lewy factor would produce. For a three-dimensional inviscid Taylor-Green vortex simulated with the compressible Euler equations and an entropy-stable flux-differencing discontinuous Galerkin scheme, the step sizes chosen by the two estimators were visually indistinguishable, with the residual approach rejecting three steps and the embedded method rejecting two.
A separate experiment confirmed the reliability claim directly. For a linear equation with a known one-sided Lipschitz constant and a nonlinear equation where that constant is zero, the global error estimate derived from the residual bound proved to be a genuine upper bound on the true error in all tested cases, with the expected scaling of the global error as tau to the power p over p plus one. This is the property that embedded estimators cannot guarantee, since they provide no rigorous link between the local difference of two Runge-Kutta solutions and the accumulated global error.
The residual approach costs more per step but delivers a certified error bound, and the study shows that this added rigor does not destabilize the feedback loop when the controller is designed appropriately. All source code and data are available in a public reproducibility repository, and the implementations use the Julia programming language with established packages for ordinary differential equations, quadrature, and discontinuous Galerkin discretizations.
The work also clarifies a subtlety in how error control and step size selection interact. Error-per-step control, which is standard for embedded methods in local extrapolation mode, does not by itself guarantee that the global error is proportional to the tolerance, but it does guarantee convergence of the error to zero as the tolerance shrinks, provided the corresponding error-per-unit-step control would be proportional. Butcher’s argument that error-per-step control produces near-optimal step size sequences supports this choice in practice.
Subject of Research: Stability of step size control for explicit Runge-Kutta methods using residual-based a posteriori error estimators
Article Title: Stability of step size control based on a posteriori error estimates
Article References: Ranocha, H., & Giesselmann, J. (2024). Stability of step size control based on a posteriori error estimates. Computational Science and Engineering, 1(1), Article 1. https://doi.org/10.1007/s44207-024-00001-0
Image Credits: AI Generated
DOI: 10.1007/s44207-024-00001-0
Keywords: ordinary differential equations, Runge-Kutta methods, step size control, a posteriori error estimates, residuals, PI controller, PID controller, stability-limited regime, embedded methods, computational fluid dynamics, Benjamin-Bona-Mahony equation, discontinuous Galerkin
News Source: Audrey Campbell. (November 6, 2026). Residual-Based Error Estimates Can Stabilize Adaptive Time Steps in Runge-Kutta Solvers. Scienmag.
