Adaptive Stepsizes

Whatever the underlying method for solving an ODE, a major problem lies in determining a good stepsize. Ideally, we want to choose $h$ as large as possible - but not so large as to give us an unreasonable amount of error, or worse still, to induce instability. If we choose a fixed stepsize, we can only proceed as fast as the "worst" sections of $x(t)$ will allow. What we would like to do is to vary $h$ as we march forward in time. Whenever we can make $h$ large without incurring too much error, we should do so. When $h$ has to be reduced to avoid excessive error, we want to do that also. This is the idea of adaptive stepsizing: varying $h$ over the course of solving the ODE.

Adaptive Stepsizes - Euler Method

Here we'll be present adaptive stepsizing for Euler's method. The basic idea as follows. Lets assume we have a given stepsize $h$, and we want to know how much we can consider changing it.

Suppose we compute two estimates for $x(t_0+h)$. We compute an estimate $x_a$, by taking an Euler step of size of size $h$ from $t_0$ to $t_0+h$. We also compute an estimate $x_b$ by taking two Euler steps of size $h/2$, from $t_0$ to $t_0+h$. Both $x_a$ and $x_b$ differ from the true value of $x(t_0+h)$ by $O(h^2)$. That means that $x_a$ and $x_b$ differ from each other by $O(h^2)$. As a result, we can write that a measure of the current $e$ is

$$e = \Vert x_a - x_b \Vert$$

Therefore, a new time step $h_{new}$ can be guessed via

$$h_{new} = \sqrt{ \frac{ e_{step} }{ e } } h$$

where $e_{step}$ is the desired accuracy per time step.

A partial implementation of this idea can be found in skeleton_adaptive_euler.c

The code illustrates the calculation of a new step size $h$ and simply returns it via the calling sequence to be used in the next time step. The logic of repeating the step for wich the error was judged to be too large is not implemented for simplicity.

Documentation on this approach can be found in Lecture Notes by Professors Witkin and Baraff from Carnegie Mellon University..

Adaptive Stepsizes for Runge-Kutta 4th order method

The same logic, a single step versus two steps to reach the time $t_o+h$ can be used with whatever ODE integrator. But Numerical Recipies argues for another approach in the case of Runge-Kutta 4th order. It amounts to a comparison between the Runge-Kutta 4th order and a solution from a carefully crafted Runge-Kutta 5th order method. Please refer to Numerical Recipies Section 16.2 for a detailed explanation of this method.