RK4, and where the magic numbers come from
Four slope samples per step, weighted 1-2-2-1, and suddenly the error falls as the fourth power. The coefficients are not arbitrary — they are what you get when you demand that error terms cancel.
Forward Euler's problem is that it commits to the slope measured at the starting point and uses it for the whole interval. If the slope changes across the step — and it always does — that is wrong from the moment you leave.
The Runge–Kutta idea is to look before you leap: take trial steps inside the interval, sample the slope there, and combine the samples so that the low-order error terms cancel.
Two samples: Heun
Step out with Euler to get a provisional endpoint, sample the slope there, then use the average of the two slopes:
This is the trapezoid rule applied to time-stepping, and it is second order. The reason it works is the same symmetry argument that made the central difference second order: averaging a slope at the start with a slope at the end cancels the term that a one-sided estimate leaves behind.
Heun costs two f-evaluations per step; forward Euler costs one. At the same total number of f-evaluations, which is more accurate on a smooth problem?
Switch the plot below to vs work — that divides out the per-step cost and compares methods at equal f-evaluations, which is the only fair axis.
Four samples: the classic
RK4 samples four times — once at the start, twice at the midpoint, once at the end — and weights them :
Those weights look like folklore. They are not. Expand the numerical step as a Taylor series in , expand the true solution the same way, and demand that the two agree through . That gives a system of order conditions on the coefficients, and is a solution. The pattern is Simpson's rule — which is not a coincidence, since for depending on alone, RK4 is Simpson's rule.
If you have met Simpson's rule for integration, you have already met the weights.
That plot is worth reading carefully, because it contains the previous chapter's lesson in a new guise. The dimmed points at the left are pre-asymptotic: at large the neglected and higher terms are not yet negligible, so the slope is not 4 there. The dimmed points at the right have hit the roundoff floor from the last chapter — around , where the truncation error is finally smaller than the noise in the arithmetic.
Fitting a line through all the points would give an answer near 4.6, which is wrong. The order is a property of the asymptotic regime, and measuring it means finding that regime first. The platform does this by computing the order between each consecutive pair and fitting only over the longest run where those local orders agree.
1/6, 1/3, 1/3, 1/6 weights come from, and why is RK4 fourth order rather than fifth?They come from order conditions. Expand the numerical step in powers of , expand the exact solution in powers of , and require the two to agree through . That imposes a system of polynomial equations on the stage coefficients and weights; the classical scheme is one solution.
It stops at fourth order because the conditions cannot also be satisfied with only four stages. In fact the number of stages needed grows faster than the order beyond this point — four stages buy order 4, but order 5 needs six stages, and order 8 needs eleven. Four stages for fourth order is the last point where the exchange rate is one-for-one, which is a large part of why RK4 became the default.
Build it
Implement one RK4 step. You have f(t, y) → array, t, y (array) and h. Compute the four stages and return the weighted combination.
Every explicit Runge–Kutta method has the same shape:
so a method is its coefficients . The Butcher tableau writes them down:
Left column is (when to sample), the block is (how to get there), the bottom row is (how to combine). Strictly lower-triangular means explicit — each stage uses only earlier ones. Fill the upper triangle and the method becomes implicit, which is the next lesson.
Two conditions are worth internalising:
- — each stage is evaluated at a consistent time. Nearly every method satisfies this.
- — the weights average the slopes rather than scaling them. Any method violating this is not even first order, and it is the first thing to check when a hand-written integrator measures order 0.
Embedded pairs are the practical payoff. Choose two weight rows and over
the same stages, giving methods of order and . Their difference estimates the
local error for almost no extra cost, and that estimate drives adaptive step size
control. Dormand–Prince 5(4) — the ode45 of MATLAB and the default in SciPy — is
exactly this: seven stages, two weight rows, order 5 with a free order-4 error estimate.
RK4 is an excellent general-purpose method and a poor special-purpose one.
- Stiff problems. RK4's stability region is bounded, so a stiff system still forces a tiny . No explicit method escapes this. → implicit methods, BDF.
- Long-time Hamiltonian integration. RK4 conserves nothing. Over a long orbit its energy drifts secularly, and a second-order symplectic method beats it comprehensively. This is the next chapter, and it is the most counterintuitive result in the path.
- Fixed step in general. Real solvers adapt. RK4 with fixed wastes effort where the solution is smooth and under-resolves where it is not.
- Highly oscillatory problems, where exponential integrators or Magnus methods exploit structure RK4 cannot see.
The habit worth building is not "use RK4" but: what do I know about this problem — stiffness, conserved quantities, smoothness — and which method exploits that?
What to carry forward
- Order comes from cancellation, engineered by choosing coefficients to kill Taylor terms.
- Higher order beats smaller steps at equal cost, on smooth problems, by a widening margin.
- Always measure the order you actually got. The measurement is the test suite for a numerical method, and it is how you catch a subtly wrong stage.
- Order is measured in the asymptotic regime — not too coarse, not below the roundoff floor.