Skip to content

Embedded Runge-Kutta (RKF45)

integrator.method: rkf45 selects an embedded Runge-Kutta pair with adaptive step-size control: one set of stages yields two solutions of different order, their difference estimates the local error, and the step size is adjusted to hold that error at a tolerance. This page derives the scheme; How HINEC uses it records what the tracker actually implements, and RKF45 Usage covers configuration.

Embedded Runge–Kutta (RKF / Dormand–Prince–type) Method

1. Problem Statement

We seek a numerical solution to the initial value problem (IVP) for a (possibly vector-valued) ODE: \( \frac{dy}{dx} = f(x, y), \quad y(x_0) = y_0, \) where \(y \in \mathbb{R}^m\) and \(f: \mathbb{R} \times \mathbb{R}^m \to \mathbb{R}^m\).

2. Explicit s-Stage Runge–Kutta Method

An explicit \(s\)-stage Runge–Kutta (RK) method advances the solution from \((x_n, y_n)\) to \((x_{n+1}, y_{n+1})\) with step size \(h\) by:

  1. Compute stages \(k_i\): \( \begin{aligned} k_1 &= f(x_n, y_n), \\ k_2 &= f\!\bigl(x_n + c_2 h,\; y_n + h a_{21} k_1 \bigr), \\ k_3 &= f\!\bigl(x_n + c_3 h,\; y_n + h (a_{31}k_1 + a_{32}k_2) \bigr), \\ &\quad \vdots \\ k_s &= f\!\bigl(x_n + c_s h,\; y_n + h (a_{s1}k_1 + \cdots + a_{s,s-1}k_{s-1}) \bigr). \end{aligned} \)

  2. Update solution: \( y_{n+1} = y_n + h \sum_{i=1}^s b_i k_i. \)

The coefficients \(a_{ij}, b_i, c_i\) define the scheme and are usually presented in a Butcher tableau:

\[ \begin{array}{c|cccc} c_1 & 0 & & & \\ c_2 & a_{21} & 0 & & \\ c_3 & a_{31} & a_{32} & 0 & \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_s & a_{s1} & a_{s2} & \cdots & a_{s,s-1} \\ \hline & b_1 & b_2 & \cdots & b_s \end{array} \]

For explicit methods we have \(a_{ij} = 0\) for \(j \ge i\).

3. Embedded RK Pair (RKF / RK5(4) Type)

An embedded Runge–Kutta pair uses a single set of stages \(k_i\) but two different sets of weights:

  • One gives an approximation \(\tilde{y}_{n+1}\) of order \(p\) (lower order),
  • The other gives an approximation \(\hat{y}_{n+1}\) of order \(q = p+1\) (higher order).

Given stages \(k_i\) from above, we form: \( \begin{aligned} \hat{y}_{n+1} &= y_n + h \sum_{i=1}^s \hat{b}_i k_i \quad \text{(higher order, e.g. 5th)}, \\ \tilde{y}_{n+1} &= y_n + h \sum_{i=1}^s b_i k_i \quad \text{(lower order, e.g. 4th)}. \end{aligned} \)

Here \(\hat{b}_i\) and \(b_i\) are the two sets of weights (often called “embedded weights”). Both approximations reuse the same \(k_i\).

In the Dormand–Prince RK5(4)7M pair (Dormand & Prince, 1980), we have:

  • \(s = 7\) stages,
  • a 5th-order formula (weights \(\hat{b}_i\)),
  • an embedded 4th-order formula (weights \(b_i\)),
  • the coefficients are chosen to:
  • satisfy the order conditions up to order 5,
  • minimize the principal truncation term for the 5th-order solution,
  • and provide a reasonably large region of absolute stability.

4. Local Error Estimate

The difference between the two embedded solutions is used as an estimate of the local truncation error:

\[ e_{n+1} = \hat{y}_{n+1} - \tilde{y}_{n+1}. \]

For an RK5(4) pair, this difference is \(O(h^5)\) and is therefore a practical estimator of the local error incurred by the 4th-order solution.

In practice, for vector-valued \(y\), a norm is used, e.g.:

\[ \text{err} = \| e_{n+1} \|. \]

The norm can be componentwise max-norm, Euclidean norm, or a weighted norm (e.g. combining absolute and relative tolerances per component).

5. Adaptive Step Size Control

Given a user-specified tolerance tol, we compare:

  • If err <= tol:
  • Accept the step: \( y_{n+1} \gets \hat{y}_{n+1} \quad (\text{use higher-order solution}), \)
  • Optionally increase the step size for the next step.
  • If err > tol:
  • Reject the step: \( y_{n+1} \text{ not updated}, \)
  • Decrease the step size and recompute.

A typical step size update rule for a RK5(4) pair is:

\[ h_{\text{new}} = \text{safety} \cdot h \left( \frac{\text{tol}}{\text{err}} \right)^{1/5}, \]

where: - safety is a factor slightly less than 1 (e.g. 0.8–0.9), - the exponent \(1/5\) corresponds to the fact that the local error scales like \(h^5\) for the 4th-order approximation.

Commonly, \(h_{\text{new}}\) is also constrained, e.g.:

\( h_{\text{new}} = \max(h_{\min}, \min(h_{\max}, h_{\text{new}})), \) and sometimes limited in growth (e.g. at most doubled per step).

6. Summary of the RKF / Dormand–Prince Step

Given \((x_n, y_n)\) and current step size \(h\):

  1. Compute stages \(k_1, \dots, k_s\) using coefficients \(a_{ij}, c_i\).
  2. Compute:
    • higher-order solution \(\hat{y}_{n+1}\) using weights \(\hat{b}_i\),
    • lower-order solution \(\tilde{y}_{n+1}\) using weights \(b_i\).
  3. Compute error estimate: \( e_{n+1} = \hat{y}_{n+1} - \tilde{y}_{n+1}, \quad \text{err} = \| e_{n+1} \|. \)
  4. If err <= tol:
    • accept step: \(x_{n+1} = x_n + h\), \(y_{n+1} = \hat{y}_{n+1}\).
  5. Choose new step size: \( h_{\text{new}} = \text{safety} \cdot h \left( \frac{\text{tol}}{\text{err}} \right)^{1/5}, \) with appropriate bounds and growth limits.
  6. If the step was rejected, retry with \(h_{\text{new}}\).

This provides automatic error control and adaptive step size using a single set of stage evaluations per step.

7. How HINEC uses it

nim_tractography_hinec implements the scheme above with the Dormand–Prince tableau, over the streamline ODE \( dx/ds = v(x) \) with \( |v| = 1 \):

  • Seven stages, each an interpolation of the direction field, sign-aligned to the incoming tangent; a stage whose interpolation fails reuses the previous stage's vector.
  • The position is advanced with the higher-order member of the pair, and the difference from the lower-order member is the error estimate, \( \text{err} = \lVert \hat{y} - \tilde{y} \rVert \), in voxels.
  • A step is accepted when err <= integrator.tolerance, and the next step is \( h_{\text{new}} = \text{safety}\cdot h\,(\text{tol}/\text{err})^{1/5} \), clamped to [integrator.step_min, integrator.step_max] and to at most \(2h\) growth. A rejected step is retried from the same position with the smaller step.
  • integrator.step is the initial step. integrator.adaptive: false keeps the pair but fixes the step — a genuine third mode, useful when the comparison of interest is the tableau rather than the step control.

rkf45 is a method name, not an order claim

The schema treats the integrator as a method selector rather than an order, because the two are not the same thing here. This implementation uses Dormand–Prince coefficients and advances on the 5th-order solution, keeping the embedded 4th-order one to size the next step — so the pair is 5(4) with local extrapolation, not a 4th-order method. More decisively, the order any scheme can realise on this problem is capped by the smoothness of the interpolated direction field, not by the tableau. See Solution Verification and High-Order Tractography.

References

  1. Dormand, J. R., & Prince, P. J. (1980). A family of embedded Runge-Kutta formulae. Journal of Computational and Applied Mathematics, 6(1), 19–26.
  2. Hairer, E., Nørsett, S. P., & Wanner, G. (1993). Solving Ordinary Differential Equations I: Nonstiff Problems (2nd ed.). Springer.