Dormand–Prince Method

Essence. The Dormand–Prince 5(4) method (DP5) is an adaptive Runge–Kutta method of orders 5 and 4, similar in spirit to RKF45 but based on a different choice of coefficients. The Dormand–Prince variant was designed so that the 5th–order solution (the main approximation) has particularly small global error, while the 4th–order solution is used only for estimating the local error. For a system \( \dot{\mathbf{x}} = f(\mathbf{x}, t, \theta) \), the method computes slope vectors \( \mathbf{k}_1, \dots, \mathbf{k}_7 \) according to the Butcher tableau, and then uses them to construct the 5th–order approximation \( \mathbf{y}_5 \), together with a 4th–order embedded estimate. The difference between the two approximations provides the local error estimate for the current step. Based on this estimate the time step \( \Delta t \) is adjusted, just as in RKF45: if the error is below the tolerance, the step is accepted and the step size may increase; if the error exceeds the tolerance, the step is rejected and the solver automatically reduces the step size.

Usage. The Dormand–Prince method is one of the most widely used adaptive schemes in practical ODE solvers (it is the basis of MATLAB's well-known “ode45”). Compared with RKF45, DP5 places stronger emphasis on minimizing the global error of the 5th–order solution, which often yields higher-quality trajectories for a similar number of steps. Thanks to adaptive step-size selection, the method is well suited for visualizing complex attractor trajectories such as Lorenz or Rössler, where the sampling density must automatically adjust to the local dynamics.

Implementation. In the Attractor Builder add-on, each step of the Dormand–Prince 5(4) method is executed by step_dp5(rhs_func, p: Vector, params: dict, current_dt: float, tolerance: float). DP5, like RKF45, uses a Butcher tableau — a structured set of coefficients defining time-shifts (matrix \( A \)) and linear combination weights (\( B \), \( B^* \)) from which two solution approximations are constructed: the 5th-order main solution and the embedded 4th-order estimate. The full Dormand–Prince 5(4) tableau is given in (Hairer, Nørsett & Wanner 1993, Table 5.2, p. 178), while the general theory of Butcher tableaux is summarized in (Butcher 2008, Chapter 5, p. 96). The first stage computes the slopes \( \mathbf{k}_1, \dots, \mathbf{k}_6 \) according to matrix \( A \): k1 = rhs_func(p, params) k2 = rhs_func(p + current_dt * A[1][0] * k1, params) k3 = rhs_func(p + current_dt * (A[2][0]*k1 + A[2][1]*k2), params) k4 = rhs_func(p + current_dt * (A[3][0]*k1 + A[3][1]*k2 + A[3][2]*k3), params) k5 = rhs_func(p + current_dt * (A[4][0]*k1 + A[4][1]*k2 + A[4][2]*k3 + A[4][3]*k4), params) k6 = rhs_func(p + current_dt * (A[5][0]*k1 + A[5][1]*k2 + A[5][2]*k3 + A[5][3]*k4 + A[5][4]*k5), params) Then the 5th–order approximation is computed using the weights \( B \): y5 = p + current_dt * (B[0]*k1 + B[1]*k2 + B[2]*k3 + B[3]*k4 + B[4]*k5 + B[5]*k6) Additionally we compute k7 = rhs_func(y5, params) which is used for the embedded 4th–order estimate via weights \( B^* \). The local error vector is all_k = [k1, k2, k3, k4, k5, k6, k7] error_vec = current_dt * sum([(b - bs) * k for b, bs, k in zip(B, B_star, all_k)], Vector()) error = error_vec.length i.e., the difference between the 5th–order main solution and the embedded 4th–order estimate. From this we compute the optimal next step: optimal_dt = 0.9 * current_dt * (tolerance / error)**0.2 If error <= tolerance, the step is accepted (the next trajectory point is y5). If the error is too large, the step is rejected and retried with optimal_dt. As with RKF45, the function step_dp5 returns the triple: (new_point, proposed_dt, acceptance_flag).

Example. Consider the Lorenz system:

Equations:
dx/dt = a * (y - x)
dy/dt = x * (b - z) - y
dz/dt = x * y - c * z

Parameters:
| a = 10 | b = 28 | c = 8/3 |

Simulation settings:
Initial state: x₀ = 0.01, y₀ = 0.01, z₀ = 0.01
Method: DP5
Tolerance: 0.000001
min_dt: 0.000001
max_dt: 0.05
Recorded steps: 3
Burn-in: 1000
Scale: 1

Using the Raw Data → Export function in the Attractor Builder add-on, we obtain CSV data with columns steps, dt, x, y, z. For the settings above, the first four rows of the CSV file are:

The first row corresponds to the initial state \( (x_0, y_0, z_0) = (0.01, 0.01, 0.01) \), with dt = 0. In the subsequent rows we see how the Dormand–Prince method adapts the step size: after the first step, dt grows from 0.01 to about 0.029 (step 3), indicating that in this region the motion is sufficiently regular to allow larger steps without violating the error tolerance.

Sources:

Butcher, J. C. (2008). Numerical Methods for Ordinary Differential Equations. Wiley.

Hairer, E., Nørsett, S., Wanner, G. (1993). Solving Ordinary Differential Equations I. Springer.