Why is the one-step-error in adaptive ODE methods defined as it is?

Question

If you are using the common adaptive ODE approach to choose time steps then you often enforce that the relative one-step-error

$\frac{|y(t_{i+1})-w_{i+1}|}{h_i} \le \epsilon$

where $y(t_{i+1})$ is the exact solution and $w_{i+1}$ is the approximation from the method and $h_i=t_{i+1}-t_i$. To me this is a little funny, since, if we wanted to control the error at $t_f$, the final time we are approximating the solution out to, then bounding the errors at each step we get

$|y(t_f) - w_{t_f}| \le \sum_i^N |y(t_{i+1})-w_{i+1}| \le \sum_i^N h_i\epsilon = (t_f-t_0)\epsilon$

where N is the number of steps taken to reach $t_f$ and $t_0$ is your starting time. That is, our control on the final error is multiplied by the total time solved over. Why isn't the one step error enforced as

$\frac{y(t_{i+1})-w_{i+1}}{h_i} \le \frac{\epsilon}{(t_f-t_0)}$

so that we can enforce that the error in our approximation at the final time is less than the tolerance?

Answer 1

The description is not complete, you do not control the evolving global error, but the local unit step error, that is, you want $$ \frac{|w_{i+1}-y_i(t_{i+1})|}{h_i}\le ϵ, $$ where $y_i(t)=\phi(t;t_i,w_i)$ is the exact solution with $y_i(t_i)=w_i$.

From local to global error

Then the global error $e_i=|w_i-y_0(t_i)|$ evolves over this step as $$ e_{i+1}\le |w_{i+1}-y_i(t_{i+1})|+|y_i(t_{i+1})-y_0(t_{i+1})| \le ϵh_i+e^{Lh_i}e_i $$ where $L$ is the Lipschitz constant of the rhs function and the Grönwall lemma was applied to the first term.

Apply the factor $e^{-Lt_{i+1}}$ to transform into a telescoping inequality $$ e^{-Lt_{i+1}}e_{i+1}\le e^{-Lt_i}e_i+ϵe^{-Lt_{i+1}}h_i \le…\le e^{-Lt_0}e_0+ϵ\sum_{j=0}^ie^{-Lt_{j+1}}h_j $$ As $e_0=0$ and the second term is a right-sided Riemann sum for a falling function, this continues as $$ e^{-Lt_{i+1}}e_{i+1}\leϵ\int_{t_0}^{t_{i+1}}e^{-Lt}\,dt =ϵ\frac{e^{-Lt_0}-e^{-Lt_{i+1}}}{L}\\~\\ e_{i+1}\leϵ\frac{e^{L(t_{i+1}-t_0)}-1}{L} $$

How to define the target error

In first order this means indeed that $ϵ$ is the desired error tolerance over a unit interval. This has the consequence that when changing the interval length, the sequence of nodes over the common parts of the intervals remains the same. If one wants the other behavior it is easy to divide the error tolerance be the length of the time interval while passing the parameters at the call of the integrator function.

On extrapolation methods

In the classical step size controller, like in the original RKF 4(5) method, one constructs a better result to estimate the error of the base method. This gives a very precise error estimate as long as the step size remains well inside the stability region. (If that becomes a concern, the ODE is called "stiff".) The above estimates are true and followed exactly for some initial segment. The local error can oscillate in its direction so that the global error can also fall back to zero at some points.

It was remarked that with the classical strategy in embedded methods one practically throws away the better, higher order value update. Using that higher order step indeed improves the global error, usually to $ϵ^{1+1/p}$ times the exponential factor.

The extrapolation idea, like used in DoPri (4)5, carries this trend one step further and uses $|w_{i+1}-\tilde w_{i+1}|\le ϵ$ as a proxy for the local unit step error of $\tilde w_{i+1}$, which is the next value of the higher-order method (if the orders are $p,p+1$). With the correct calibration this works out mostly well, especially if the method design is optimized for it.