In mathematics and numerical analysis, an adaptive step size is used in some methods for the numerical solution of ordinary differential equations (including the special case of numerical integration) in order to control the errors of the method and to ensure stability properties such as A-stability. Using an adaptive stepsize is of particular importance when there is a large variation in the size of the derivative. For example, when modeling the motion of a satellite about the earth as a standard Kepler orbit, a fixed time-stepping method such as the Euler method may be sufficient. However things are more difficult if one wishes to model the motion of a spacecraft taking into account both the Earth and the Moon as in the Three-body problem. There, scenarios emerge where one can take large time steps when the spacecraft is far from the Earth and Moon, but if the spacecraft gets close to colliding with one of the planetary bodies, then small time steps are needed. Romberg's method and Runge–Kutta–Fehlberg are examples of a numerical integration methods which use an adaptive stepsize.
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
y'(t)=f(t,y(t)), y(a)=ya
We are given the function f(t,y) and the initial conditions (a, ya), and we are interested in finding the solution at t = b. Let y(b) denote the exact solution at b, and let yb denote the solution that we compute. We write
yb+\varepsilon=y(b)
\varepsilon
For a sequence (tn) of values of t, with tn = a + nh, the Euler method gives approximations to the corresponding values of y(tn) as
(0) | |
y | |
n+1 |
=yn+hf(tn,yn)
(0) | |
\tau | |
n+1 |
=y(tn+1)-
(0) | |
y | |
n+1 |
(0) | |
\tau | |
n+1 |
=ch2
We have marked this solution and its error with a
(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(tn+1). We get a second solution, which we label with a
(1)
y | ||||
|
=y | ||||
|
f(tn,yn)
(1) | |
y | |
n+1 |
=y | + | ||||
|
h | |
2 |
f(t | ||||
|
,y | ||||
|
)
(1) | ||
\tau | =c\left( | |
n+1 |
h | |
2 |
| ||||
\right) |
| ||||
\right) |
| ||||
\right) |
| ||||
ch |
(0) | |
\tau | |
n+1 |
(1) | |
y | |
n+1 |
+
(1) | |
\tau | |
n+1 |
=y(t+h)
Here, we assume error factor
c
[t,t+h]
y(3)(t)
(1) | |
y | |
n+1 |
(0) | |
-y | |
n+1 |
=
(1) | |
\tau | |
n+1 |
This local error estimate is third order accurate.
The local error estimate can be used to decide how stepsize
h
tol
h → 0.9 x h x min\left(max\left(\left(
tol | |||||||||
|
\right)1/2,0.3\right),2\right)
The
0.9
0.9 x tol
(1) | |
|\tau | |
n+1 |
|<tol
(2) | |
y | |
n+1 |
=
(1) | |
y | |
n+1 |
+
(1) | |
\tau | |
n+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
h=b-a
a
b
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.
Adaptive stepsize methods that use a so-called 'embedded' error estimate include the Bogacki–Shampine, Runge–Kutta–Fehlberg, Cash–Karp and Dormand–Prince methods. These methods are considered to be more computationally efficient, but have lower accuracy in their error estimates.
To illustrate the ideas of embedded method, consider the following scheme which update
yn
yn+1=yn+hn\psi(tn,yn,hn)
tn+1=tn+hn
hn
hn=g(tn,yn,hn-1)
For embedded RK method, computation of
\psi
\tilde{\psi}
rm{err}n(h)=\tilde{y}n+1-yn+1=h(\tilde{\psi}(tn,yn,hn)-\psi(tn,yn,hn))
rm{err}n
rm{tol}n=rm{Atol}+rm{Rtol} ⋅ max(|yn|,|yn-1|)
En=rm{norm}(rm{err}n/rm{tol}n)
En
hn
hn=hn-1
1/(q+1) | |
(1/E | |
n) |
\tilde{\psi}
The description given above is a simplified procedures used in the stepsize control for explicit RK solvers. A more detailed treatmentcan be found in Hairer's textbook.[1] The ODE solver in many programming languages uses this procedure as the default strategy for adaptive stepsize control, which adds other engineering parameters to make the system more stable.