33. Initial Value Problems for Ordinary Differential Equations, Part 5: Error Control and Variable Step Sizes — Preview#

References:

33.1. The Basic ODE Initial Value Problem#

We consider again the initial value problem

\[ \frac{d u}{d t} = f(t, u) \quad a \leq t \leq b, \quad u(a) = u_0 \]

We now allow the possibility that \(u\) and \(f\) are vector-valued as in the section on Systems of ODEs and Higher Order ODEs, but omitting the tilde notation \(\tilde u\), \(\tilde f\).

33.2. Error Control by Varying the Time Step Size \(h_i\)#

Recall the variable step-size version of Euler’s method:

Input: \(f\), \(a\), \(b\), \(n\)

\(t_0 = a\)
\(U_0 = u_0\)
\(h = (b-a)/n\)

for i in \([0, n)\):

\(\qquad\) Choose step size \(h_i\) somehow!

\(\qquad\) \(t_{i+1} = t_i + h_i\)

\(\qquad\) \(U_{i+1} = U_i + h_i f(t_i, U_i)\)

end for

We now consider how to choose each step size, by estimating the error in each step, and aiming to have error per unit time below some limit like \(\epsilon/(b-a)\), so that the global error is no more than about \(\epsilon\).

As usual, the theoretical error bounds like \(O(h_i^2)\) for a single step of Euler’s method are not enough for quantitative tasks like choosing \(h_i\), but they do motivate more practical estimates.

33.3. A crude error estimator for Euler’s Method: Richardson Extrapolation#

Starting at a point (t, u(t)), we can estimate the error in Euler’s method approximato at a slightly later time \(t_i + h\) by using two approximations of \(U(t + h)\):

  • The value given by a step of Euler’s method with step size \(h\): call this \(U^{h}\)

  • The value given by taking two steps of Euler’s method each with step size \(h/2\): call this \(U_2^{h/2}\), because it involves 2 steps of size \(h/2\).

The first order accuracy of Euler’s method gives \(e_h := u(t+h) - U^{h} \approx 2(u(t+h) - U_2^{h/2})\), so that

\[e_h \approx \frac{U_2^{h/2} - U^{h}}{2}\]

33.3.1. Step size choice#

What do we do with this error information?

The first obvious ideas are:

  • Accept this step if \(e_h\) is small enough, taking \(h_i = h\), \(t_{i+1} = t_i + h_i\), and \(U_{i+1} = U^h\), but

  • reject it and try again with a smaller \(h\) value otherwise.

33.3.1.1. Exercise 1#

Write a formula for \(U_h\) and \(e_h\) if one starts from the point \((t_i, U_i)\), so that \((t_i + h, U^h)\) is the proposed value for the next point \((t_{i+1}, U_{i+1})\) in the approximate solution — but only if \(e_h\) is small enough!

33.3.2. Error tolerance#

One simple criterion for accuracy is that the estimated error per unit time in this step be no more than some overall bound on error per unit time bound, \(\epsilon\). That is, accept the step size \(h\) if

\[\frac{|e_h|}{h} \leq \epsilon.\]

33.3.3. A crude approach to reducing the step size when needed#

If this error tolerance is not met, we must choose a new step size \(h'\), and we can predict roughly its error behavior using the known order natue of the error in Euler’s method: scaling dowen to \(h' = s h\), the error in a step scales with \(h^2\), so the error per unit time scales with \(h\) (in general it scales with \(h^p\) for a method of order \(p\)), and so to reduce the error by the needed factor \(\displaystyle \frac{e_h/h}{\epsilon}\) one needs approximately

\[ s = \frac{\epsilon}{|e_h|/h} \approx \frac{2 h \epsilon}{|U^{h/2} - U^{h}|} \]

and so

\[ h' = \frac{2 h^2 \epsilon}{|U^{h/2} - U^{h}|}. \]

However this new step size might have error that is still slightly too large, leading to a second failure. Another is that one might get into an infinite loop of step size reduction.

So refinement of this choice must be considered.

33.3.4. Increasing the step size when desirable#

If we simply follow the above aproach, the step size, once reduced, will never be increased. This could lead to great inefficiency, through using an unecessarily small step size just because at an earlier part of the time domain, accuracy required very small steps.

Thus, after a successful time step, one might consider increasing \(h\) for the next step. This could be done using exactly the above formula, but again there are risks, so again refinement of this choice must be considered.

One problem is that if the step size gets too large, the error estimate can become unreliable; another is that one might need some minimum “temporal resolution”, for nice graphs and such.

Both suggest imposing an upper limit on the step size \(h\).

This work is licensed under Creative Commons Attribution-ShareAlike 4.0 International