Adaptive step size

Example

For simplicity, the following example uses the simplest integration method, the Euler method; in practice, higher-order methods such as Runge–Kutta methods are preferred due to their superior convergence and stability properties.

Consider the initial value problem $${\displaystyle y^{\prime }(t)=f(t,y(t)),y(a)=y_a}$$ where y and f may denote vectors (in which case this equation represents a system of coupled ODEs in several variables).

We are given the function f(t,y) and the initial conditions (a, y_a), and we are interested in finding the solution at t = b. Let y(b) denote the exact solution at b, and let y_b denote the solution that we compute. We write ${\displaystyle y_b+\epsilon =y(b)}$, where ${\displaystyle \epsilon }$ is the error in the numerical solution.

For a sequence (t_n) of values of t, with t_n = a + nh, the Euler method gives approximations to the corresponding values of y(t_n) as $${\displaystyle y_{n+1}^{(0)}=y_n+hf(t_n,y_n)}$$ The local truncation error of this approximation is defined by $${\displaystyle {\tau }_{n+1}^{(0)}=y(t_{n+1})-y_{n+1}^{(0)}}$$ and by Taylor's theorem, it can be shown that (provided f is sufficiently smooth) the local truncation error is proportional to the square of the step size: $${\displaystyle {\tau }_{n+1}^{(0)}=c{h}^2}$$ where c is some constant of proportionality.

We have marked this solution and its error with a ${\displaystyle (0)}$.

The value of c is not known to us. Let us now apply Euler's method again with a different step size to generate a second approximation to y(t_n+1). We get a second solution, which we label with a ${\displaystyle (1)}$. Take the new step size to be one half of the original step size, and apply two steps of Euler's method. This second solution is presumably more accurate. Since we have to apply Euler's method twice, the local error is (in the worst case) twice the original error. $${\displaystyle \begin{array}{cc}{y}_{n+\frac{1}{2}}& =y_n+\frac{h}{2}f(t_n,y_n)\\ y_{n+1}^{(1)}& =y_{n+\frac{1}{2}}+\frac{h}{2}f(t_{n+\frac{1}{2}},y_{n+\frac{1}{2}})\\ {\tau }_{n+1}^{(1)}& =c{\left(\frac{h}{2}\right)}^2+c{\left(\frac{h}{2}\right)}^2=2c{\left(\frac{h}{2}\right)}^2=\frac{1}{2}c{h}^2=\frac{1}{2}{\tau }_{n+1}^{(0)}\\ y_{n+1}^{(1)}+{\tau }_{n+1}^{(1)}& =y(t+h)\end{array}}$$

Here, we assume error factor ${\displaystyle c}$ is constant over the interval ${\displaystyle [t,t+h]}$. In reality its rate of change is proportional to ${\displaystyle y^{(3)}(t)}$. Subtracting solutions gives the error estimate:

$${\displaystyle y_{n+1}^{(1)}-y_{n+1}^{(0)}={\tau }_{n+1}^{(1)}}$$

This local error estimate is third order accurate.

The local error estimate can be used to decide how stepsize ${\displaystyle h}$ should be modified to achieve the desired accuracy. For example, if a local tolerance of ${\displaystyle \text{tol}}$ is allowed, we could let h evolve like:

$${\displaystyle h\to 0.9\times h\times min\left(max\left({\left(\frac{\text{tol}}{2\left|{\tau }_{n+1}^{(1)}\right|}\right)}^{1/2},0.3\right),2\right)}$$

The ${\displaystyle 0.9}$ is a safety factor to ensure success on the next try. The minimum and maximum are to prevent extreme changes from the previous stepsize. This should, in principle give an error of about ${\displaystyle 0.9\times \text{tol}}$ in the next try. If ${\displaystyle |{\tau }_{n+1}^{(1)}|<\text{tol}}$, we consider the step successful, and the error estimate is used to improve the solution:

$${\displaystyle y_{n+1}^{(2)}=y_{n+1}^{(1)}+{\tau }_{n+1}^{(1)}}$$

This solution is actually third order accurate in the local scope (second order in the global scope), but since there is no error estimate for it, this doesn't help in reducing the number of steps. This technique is called Richardson extrapolation.

Beginning with an initial stepsize of ${\displaystyle h=b-a}$, this theory facilitates our controllable integration of the ODE from point ${\displaystyle a}$ to ${\displaystyle b}$, using an optimal number of steps given a local error tolerance. A drawback is that the step size may become prohibitively small, especially when using the low-order Euler method.

Similar methods can be developed for higher order methods, such as the 4th-order Runge–Kutta method. Also, a global error tolerance can be achieved by scaling the local error to global scope.