Fourth-Order Runge-Kutta and Adaptive Steps

Definition The classical Runge-Kutta method, often called RK4, advances from $(t_n, y_n)$ with four stage slopes, computed in this order,

and then sets

The stages $k_2$ and $k_3$ share the midpoint time $t_n + h/2$ but are read at the predicted states $y_n + \tfrac{h}{2}k_1$ and $y_n + \tfrac{h}{2}k_2$, which in general differ, and $k_4$ is read at the end of the step, at a state predicted with $k_3$. Every stage state starts from the same value $y_n$. The predicted states are not consecutive solution points, and only $y_{n+1}$ is kept. Each step costs four evaluations of $f$.

Definition The RK4 update moves from $y_n$ by $h$ times the weighted average

of its stage slopes, with the weights $\tfrac16$, $\tfrac13$, $\tfrac13$ and $\tfrac16$. The weights are positive and add up to $1$, so a constant slope $f = m$ gives the exact change $hm$. Each $k_i$ is a slope, and $h$ multiplies the average once. Some books use the increments $K_i = hk_i$ instead, so that the update adds $(K_1 + 2K_2 + 2K_3 + K_4)/6$ and the second stage is $K_2 = h f(t_n + h/2, y_n + K_1/2)$. Mixing the two conventions applies $h$ twice or not at all. When $f$ depends on $t$ alone, the two middle slopes are equal and the step is Simpson’s rule for the integral of $f$ over $[t_n, t_{n+1}]$, with the weights $\tfrac16$, $\tfrac46$ and $\tfrac16$ at the start, the middle and the end.

Definition Let $f$ have continuous partial derivatives through order four on an open set containing the graph of the exact solution over $[t_0, T]$, a closed interval inside its interval of existence. Then RK4 has local truncation error of order $h^5$, and there are constants $C$ and $h_0 > 0$ such that

for every grid time $t_n \le T$ whenever $0 < h \le h_0$. In a typical problem, once $h$ is small, halving it divides the global error by about $16$. On $y' = y$ one step multiplies by

the Taylor polynomial of $e^h$ of degree four, so $e^h - R(h) = h^5/120 + O(h^6)$. Matching this one equation is a check, not a proof. The general result also needs the order conditions for right sides that depend nonlinearly on $t$ and $y$, and the classical stages and weights satisfy all of them. The order also describes small steps only. On $y' = -ay$ with $a > 0$, one step multiplies by $R(-ah)$, which is positive for every $h$, and the computed values decay exactly when $ah$ is below a boundary value $\beta \approx 2.785$.

Definition A solver estimates the local error of a step by computing a second approximation of different accuracy at the same end time. Step doubling does this with one method. With RK4, compute $y_c$ by one step of size $h$ and $y_f$ by two steps of size $h/2$, both from $(t_n, y_n)$ to $t_{n+1} = t_n + h$. If the coarse local error is about $Ch^5$, the fine result carries two errors of about $C(h/2)^5$ each, so to leading order

where $y$ here means the exact solution through $(t_n, y_n)$. An embedded Runge-Kutta pair gets a similar estimate more cheaply. One set of stage slopes feeds two update formulas of different orders, usually $q$ and $q + 1$, and the difference of the two results estimates the local error of the lower-order one. Well-known pairs are the Fehlberg pair of orders $4$ and $5$ with six stages, and the Dormand-Prince pair of orders $5$ and $4$, which has seven stages but reuses its last stage as the first stage of the next step. Every such estimate rests on the leading term of an error expansion, so it is an approximation and not a guaranteed bound.

Definition An absolute tolerance $\text{atol} \ge 0$ and a relative tolerance $\text{rtol} \ge 0$, not both zero, set the error scale of a step from $y_n$ to $y_{n+1}$,

and the step passes when its estimated local error $\delta$ satisfies $\lvert \delta \rvert \le s$. The absolute part controls the error while the solution is near $0$, and the relative part measures the error against the size of the solution when it is large. Some solvers use $\max(\text{atol}, \text{rtol} \cdot m)$ in place of the sum, which changes the scale by at most a factor of $2$. For a system each component gets its own scale, often with its own absolute tolerance in that component’s units, and the ratios are combined by a maximum or a root mean square. The tolerances control the estimated error of each step. They do not guarantee that the global error at the end is below the tolerance, because local errors add up over many steps and can grow as the solution evolves.

Definition An adaptive solver changes $h$ from step to step to keep its estimated local error near the requested size. After a trial step of size $h$ it forms the error estimate $\delta$ and the scale $s$, and computes the ratio $\rho = \lvert \delta \rvert / s$. If $\rho \le 1$ the step is accepted, and otherwise it is rejected and repeated with a smaller step. In both cases the next trial step is

where $q$ is the lower of the two orders, since that formula’s local error scales like $h^{q+1}$. The factor $0.9$ is a safety margin, and the new step is kept between fixed multiples of the old one, such as $0.2h$ and $5h$. Many solvers keep the higher-order result as the accepted value, so the estimate is cautious for the value actually kept. Known jumps in $f$ or switches in the model belong at the ends of separate integration intervals, or are located as events. Steps that keep shrinking and repeated rejections are information about the problem, for example a nearby singularity or stiffness, and not something to hide by loosening the tolerances.