2.7. Adaptive step size control#
The implementations of the Runge-Kutta methods that we have used previously in this chapter have used a constant value of the step length \(h\). Using a fixed step length does not allow the method to take advantage of increasing \(h\) where the behaviour of the solution allows. Adaptive step size control is a method that attempts to control the value of \(h\) based on accuracy requirements. The accuracy of the approximation to the solution to the ODE is dependent upon three factors:
the order of accuracy of the computational method;
the size of the step length h used to advance the solution;
the behaviour of the solution.
Improving the order of accuracy of the method is often not straightforward and can place restrictions on the applicability of the method. Reducing the step length used will improve the accuracy of the approximation but since more steps are required to advance through the domain, this will increase the computational cost. The behaviour of the solution will have an affect on the accuracy of the approximation because computational methods are more accurate where the solution is slowly varying and less accurate where there are rapid variations in the solution.
With adaptive step size control we compute two solutions of the ODE over a single step using a step length \(h\) using an order \(p\) method and an order \(p+1\) method and calculate the difference between the two solutions. If this difference is less than a desired accuracy we say that the step has succeeded else we say that the step has failed. If the step was successful we increase \(h\) and move to the next step. If the step had failed we decrease \(h\) and repeat the same step.
2.7.1. Calculating the new step length#
At each step in the method we will be adjusting the value of the step length \(h\). If a step has failed we want to reduce the value of \(h\) so that the solution using the \(p\) order method is more accurate and therefore closer to the \(p+1\) order solution. If a step has succeeded this indicates the the step length used could have been larger so we want to increase \(h\) but not too much as to result in a failed step in the next step.
Let \(y^{(p)}\) and \(y^{(p+1)}\) denote the numerical solutions calculated using the order \(p\) and \(p+1\) method respectively and \(y\) denote the exact solution then
Subtracting the second equation from the first gives
If we let \(\Delta = |y^{(p)} - y^{(p+1)}|\) and since \(O(h^{p+1}) + O(h^{p+2}) = O(h^{p+1})\) (the larger error dominates the smaller one) then by the definition of \(O(h^n)\)
We want to calculate a new step length \(h_{new}\) which results in a desired accuracy tolerance \(tol\), i.e.,
and defining a ratio \(r\) between the new step length and current step length, \(r = h_{new} / h\), then
If \(\Delta > tol\) then we need to recalculate the current step which incurs additional computational cost. To guard against this as much as possible a safety factor is usually introduced into the computation. If \(h_{new}\) is the estimated step size to to give \(\Delta = tol\) then some smaller value, such as \(0.9r\) is typically used instead. Another consideration is that for robustness we do not want the step length to change too much from one step to the next so we limit \(r\) so that it is in the range \([\frac{1}{2}, 2]\), i.e.,
2.7.2. Embedded Runge-Kutta methods#
German mathematician Erwin Fehlberg (1911 - 1990) derived a Runge-Kutta method where different weightings applied to the stage values can result in a method of fourth or fifth-order without having to repeat the calculations of the ODE function. Since Fehlberg’s original method, several other Runge-Kutta methods that use the same stage values for producing solutions of different orders have been suggested and these methods are known as embedded Runge-Kutta methods.
The Butcher tableau for an embedded Runge-Kutta method takes the form
Where the weights \(b_i\) and \(b_i^*\) that give an order \(p+1\) and \(p\) method respectively, i.e.,
where \(k_i\) are the same for both equations. Some common embedded methods are defined below.
Definition 2.8 (Fehlberg’s order 4(5) Runge-Kutta method (RKF45))
Fehlberg [1969]
Definition 2.9 (Cash-Karp order 4(5) Runge-Kutta method (RKCK45))
Cash and Karp [1990]
Algorithm 2.2 (Solving an IVP using a single step method with adaptive step size control)
Inputs A system of first-order ODEs of the form \(\vec{y}' = f(t,\vec{y})\), a domain \(t \in [t_{\min}, t_{\max}]\), a vector of initial values \(\vec{y}(t_{\min})\) an initial value of the step length \(h\) and a desired accuracy tolerance \(tol\).
Outputs \(\vec{t}\) and \(\vec{y}\).
\(t_0 \gets t_{\min}\)
\(\vec{y}_0 \gets \vec{y}(t_{\min})\)
\(n \gets 0\)
While \(t_n < t_{\max}\)
\(h \gets \min(h, t_{\max} - t_{n})\)
Calculate \(\vec{y}^{(p+1)}\) and \(\vec{y}^{(p)}\) solutions
\(\Delta \gets \max | \vec{y}^{(p+1)} - \vec{y}^{(p)} |\)
If \(\Delta < tol\)
\(\vec{y}_{n+1} \gets \textbf{y}^{(p+1)}\)
\(t_{n+1} \gets t_n + h\)
\(n \gets n + 1\)
\(h \gets \max\left(0.5, \min \left(2, 0.9 \left( \dfrac{tol}{\Delta} \right) ^{1/(p+1)} \right) \right)\)
Return \((t_0, t_1, \ldots, t_n)\) and \((\vec{y}_0, \vec{y}_1, \ldots, \vec{y}_n)\)
Example 2.8
The Bogacki-Shampine 2(3) embedded Runge-Kutta method has the Butcher tableau
where the first row of \(b\) values gives a third-order method and the second row gives a second-order method. Compute the first step of the Bogacki-Shampine method for the following initial value problem using an initial step length of \(h_0 = 0.1\) and an accuracy tolerance of \(tol = 10^{-4}\)
Solution (click to show)
Calculating the stage values
Calculate \(y^{(p+1)}\) and \(y^{(p)}\)
Calculate \(\Delta\)
\(\Delta = 0.010566 > 10^{-4}\) so the step has failed. Recalculating \(h\)
Calculating the stage values using the new value of \(h = 0.05\)
Calculate \(y^{(p+1)}\) and \(y^{(p)}\)
Calculate \(\Delta\)
\(\Delta = 5.83 \times 10^{-5} < 10^{-4}\) so the step has been successful and \(y_1 = y^{(p+1)} = 0.032140\).
Calculate the new value of \(h\) for the next step
The solution to this initial value problem using the Bogacki-Shampine 2(3) method is shown in the table below along with the step lengths used for each step.
\(n\) |
\(t_n\) |
\(y_n\) |
\(h_n\) |
|---|---|---|---|
0 |
0.000000 |
0.000000 |
- |
1 |
0.050000 |
0.032140 |
0.050000 |
2 |
0.103880 |
0.040939 |
0.053880 |
3 |
0.161862 |
0.041599 |
0.057982 |
4 |
0.239599 |
0.039342 |
0.077737 |
5 |
0.333844 |
0.035754 |
0.094244 |
6 |
0.466041 |
0.031259 |
0.132197 |
7 |
0.598661 |
0.027477 |
0.132620 |
8 |
0.725978 |
0.024064 |
0.127317 |
9 |
0.852679 |
0.021364 |
0.126701 |
10 |
0.962172 |
0.019014 |
0.109494 |
11 |
1.000000 |
0.018354 |
0.037828 |
2.7.3. Code#
The code below defines the function solveIVP_SSC() which computes the solution to an initial value problem using step size control. This uses the same inputs as the solveIVP() with an additional input of tol which specifies the accuracy tolerance. Since we do not know beforehand how many steps will be required we assume the number steps will be large a while loop instead of a for loop for stepping through the solution.
The rkf45() function computes a single step of the Fehlberg 4(5) method. This function returns the values of the fourth and fifth-order solutions y4 and y5 as well as the order of the lower order method \(p\) which in this case is 4. This is used to compute the value of the step length for the next step.
def solveIVP_SSC(f, tspan, y0, h, solver, tol=1e-4):
# Define t and y arrays
t = np.zeros(10000)
y = np.zeros((10000,len(y0)))
t[0] = tspan[0]
y[0,:] = y0
# Solver loop
n = 0
while t[n] < tspan[-1]:
# Ensure t does not exceed tmax
h = min(h, tspan[-1] - t[n])
# Calculate order p and p+1 solutions
yp1, yp, p = solver(f, t[n], y[n,:], h)
# Determine whether step was successful or not
delta = np.max(np.abs(yp1 - yp))
if delta < tol:
y[n+1,:] = yp1
t[n+1] = t[n] + h
n += 1
# Calculate new value of h
h *= max(0.5, min(2, 0.9 * (tol / delta) ** (1 / (p + 1))))
return t[:n+1], y[:n+1,:]
def rkf45(f, t, y, h):
k1 = f(t, y)
k2 = f(t + 1/4 * h, y + 1/4 * h * k1)
k3 = f(t + 3/8 * h, y + h * (3/32 * k1 + 9/32 * k2))
k4 = f(t + 12/13 * h, y + h * (1932/2197 * k1 - 7200/2197 * k2 + 7296/2197 * k3))
k5 = f(t + h, y + h * (439/216 * k1 - 8 * k2 + 3680/513 * k3 - 845/4104 * k4))
k6 = f(t + 1/2 * h, y + h * (-8/27 * k1 + 2 * k2 - 3544/2565 * k3 + 1859/4104 * k4 - 11/40 * k5))
y5 = y + h * (16/135 * k1 + 6656/12825 * k3 + 28561/56430 * k4 - 9/50 * k5 + 2/55 * k6)
y4 = y + h * (25/216 * k1 + 1408/2565 * k3 + 2197/4104 * k4 - 1/5 * k5);
return y5, y4, 4
function [t, y] = solveIVP_SSC(f, tspan, y0, h, solver, tol)
% Define t and y arrays
t = zeros(10000);
y = zeros(10000, length(y0));
t(1) = tspan(1);
y(1,:) = y0;
% Loop through the steps
n = 1;
while t(n) < tspan(2)
% Ensure t does not exceed tmax
h = min(h, tspan(2) - t(n));
% Calculate order p and p+1 solutions
[yp1, yp, order] = solver(f, t(n), y(n,:), h);
% Determine whether the step was successful or not
delta = max(abs(yp1 - yp));
if delta < tol
y(n+1,:) = yp1;
t(n+1) = t(n) + h;
n = n + 1;
end
% Calculate new value of h
h = h * max(0.5, min(2, 0.9 * (tol / delta) ^ (1 / (order + 1))));
end
% Remove unused entries from t and y
t(n+1:end) = [];
y(n+1:end,:) = [];
end
function [y5, y4, order] = rkf45(f, t, y, h)
k1 = f(t, y);
k2 = f(t + 1/4 * h, y + 1/4 * h * k1);
k3 = f(t + 3/8 * h, y + h * (3/32 * k1 + 9/32 * k2));
k4 = f(t + 12/13 * h, y + h * (1932/2197 * k1 - 7200/2197 * k2 + 7296/2197 * k3));
k5 = f(t + h, y + h * (439/216 * k1 - 8 * k2 + 3680/513 * k3 - 845/4104 * k4));
k6 = f(t + 1/2 * h, y + h * (-8/27 * k1 + 2 * k2 - 3544/2565 * k3 + 1859/4104 * k4 - 11/40 * k5));
y5 = y + h * (16/135 * k1 + 6656/12825 * k3 + 28561/56430 * k4 - 9/50 * k5 + 2/55 * k6);
y4 = y + h * (25/216 * k1 + 1408/2565 * k3 + 2197/4104 * k4 - 1/5 * k5);
order = 4;
end
Fehlberg’s method has been applied to solve the initial value problem from Example 2.8 and the solution is plotted in Fig. 2.9. The solution had 11 successful steps, 3 failed steps and a total of 84 function evaluations.
Fig. 2.9 Fehlberg method solution to the initial value problem \(y' = -21y + e^{-t}\), \(t\in [0,1]\), \(y(0)=0\) using Fehlberg’s method with an initial step length \(h=0.1\) and an accuracy tolerance \(tol=10^{-4}\).#