What RK4 does and why it is the default
The classical fourth-order Runge-Kutta method advances an initial value problem dy/dx = f(x, y), y(x₀) = y₀ one step at a time. At each step it samples the slope four times - once at the start, twice near the middle, once at the end - and takes a weighted average of those four slopes as the effective slope for the whole step.
That averaging is what buys the accuracy. The Euler method uses only the slope at the start of the step, so it consistently lags a curving solution; the global error falls only in proportion to h. RK4 matches the Taylor expansion of the true solution through the h4 term, so its global error falls as h4. Halve the step and Euler's error halves; RK4's error drops by a factor of about sixteen.
The other reason RK4 became the default is that it needs nothing but f itself. No derivatives of f, no algebraic solve inside the step, no history of past points - so it starts itself, handles a change of step size without ceremony, and works on any equation you can evaluate. It costs four function evaluations per step, which is the price.
This calculator also runs Euler over the same interval with the same step, and runs RK4 again at half the step to estimate the error, so you can see what the extra work is buying on your particular problem.
The four stages, and where the weights come from
One step from (xₙ, yₙ) to (xₙ₊₁, yₙ₊₁):
k₁ = f(xₙ, yₙ)- the slope at the left-hand end, exactly what Euler would use.k₂ = f(xₙ + h/2, yₙ + h k₁/2)- the slope at the midpoint, reached by taking a half Euler step with k1.k₃ = f(xₙ + h/2, yₙ + h k₂/2)- the midpoint slope again, but reached using the better estimate k2.k₄ = f(xₙ + h, yₙ + h k₃)- the slope at the right-hand end, reached with k3.
and then
yₙ₊₁ = yₙ + (h/6)(k₁ + 2k₂ + 2k₃ + k₄)
The weights 1, 2, 2, 1 divided by 6 are Simpson's rule. That is not a coincidence: if f happens to depend on x only, the differential equation is just an integral, and the four stages collapse to Simpson's rule applied on the step. RK4 is Simpson's rule generalised to the case where the integrand depends on the answer.
The weights are not free either - they are forced. Expand y(xₙ + h) as a Taylor series in h, expand the four-stage formula the same way, and demand that the two agree through h4. That produces a system of eleven order conditions in the stage coefficients, and the tableau above is the classical solution Wilhelm Kutta published in 1901. Other solutions exist - the 3/8 rule is the best known - and they differ in their error constants rather than in their order.
Because the expansion is matched through the fourth power, the error left over in a single step is proportional to h5. Over an interval of fixed length you take a number of steps proportional to 1/h, so those local errors accumulate into a global error proportional to h4. That is the distinction between local truncation error and global error, and it is why a fourth-order method is described by its global exponent.
Worked example: y′ = y − x² + 1, y(0) = 0.5, h = 0.2
This is the standard textbook test problem, and it has a closed-form solution y = (x + 1)² − 0.5ex to check against. Take the first step by hand.
- k₁. f(0, 0.5) = 0.5 − 0 + 1 = 1.5.
- k₂. The half step gives y = 0.5 + 0.2(1.5)/2 = 0.65, so f(0.1, 0.65) = 0.65 − 0.01 + 1 = 1.64.
- k₃. Now y = 0.5 + 0.2(1.64)/2 = 0.664, so f(0.1, 0.664) = 0.664 − 0.01 + 1 = 1.654.
- k₄. The full step gives y = 0.5 + 0.2(1.654) = 0.8308, so f(0.2, 0.8308) = 0.8308 − 0.04 + 1 = 1.7908.
- Combine. k₁ + 2k₂ + 2k₃ + k₄ = 1.5 + 3.28 + 3.308 + 1.7908 = 9.8788. Multiply by h/6 = 0.2/6 = 0.0333333: 9.8788 × 0.0333333 = 0.3292933. Add to y₀: y₁ = 0.5 + 0.3292933 = 0.8292933.
The exact solution at x = 0.2 is (1.2)² − 0.5e0.2 = 1.44 − 0.5(1.2214028) = 1.44 − 0.6107014 = 0.8292986. The RK4 step is short by 0.8292986 − 0.8292933 = 0.0000053, about five parts in a million after one step of size 0.2.
Compare Euler on the same step: y₁ = 0.5 + 0.2(1.5) = 0.8, an error of 0.0292986 - roughly 5,500 times larger, for one quarter of the arithmetic. Continue to x = 2 and the exact answer is 9 − 0.5e² = 9 − 3.6945280 = 5.3054720, which RK4 tracks closely all the way while Euler drifts steadily below it.
Choosing h and reading the error estimate
The error estimate reported here is Richardson extrapolation: run the whole interval at h, run it again at h/2, and take the difference divided by 15. The 15 is 24 − 1, and it comes straight from the fourth-order behaviour - if the error at h is Ch⁴, the error at h/2 is Ch⁴/16, so their difference is 15/16 of the larger error.
Use it as a target-driven step selector. If the estimate is ten times bigger than you can tolerate, cut h by a factor of about 1.8, because 1.84 ≈ 10.5. If it is a thousand times smaller than you need, you are paying for accuracy nobody asked for and can raise h by a factor of about 5.6.
Two situations break that reasoning. The first is a stiff equation, where the solution contains a component decaying far faster than the one you care about - y′ = −1000y + cos x is the classic. RK4 is only conditionally stable, so it needs h small enough to resolve the fast component even after that component has died away; exceed the stability limit and the numerical solution oscillates and explodes rather than gently losing accuracy. An implicit method is the answer, not a smaller explicit step.
The second is a solution that is not smooth. RK4's order argument assumes f has four continuous derivatives along the solution. Across a kink, a switch, or a discontinuity in a forcing term, the method silently drops to first order. Stop the integration at the discontinuity and restart it on the other side.
Watch the Euler comparison as a rough diagnostic of curvature. When Euler and RK4 nearly agree, the solution is close to straight over a step and almost any method would do. When they differ sharply, the solution is curving hard relative to your step size, and it is worth checking the Richardson estimate before trusting the result.
How the error falls with step size
| Method | Order | f evaluations per step | Error factor when h is halved | Steps needed for 10× less error |
|---|---|---|---|---|
| Euler | 1 | 1 | ÷ 2 | 10× as many |
| Heun / improved Euler | 2 | 2 | ÷ 4 | 3.2× as many |
| Midpoint (RK2) | 2 | 2 | ÷ 4 | 3.2× as many |
| Classical RK4 | 4 | 4 | ÷ 16 | 1.8× as many |
| Runge-Kutta 3/8 rule | 4 | 4 | ÷ 16 | 1.8× as many |
| Dormand-Prince RK45 | 5 (4 for the error estimate) | 6, one reused | ÷ 32 | 1.6× as many |
The last column is 10^(1/p): the factor by which the number of steps must rise to cut the global error tenfold for a method of order p.
Pitfalls
- Confusing local and global error. One RK4 step is accurate to O(h5); a whole interval is accurate to O(h4), because you take O(1/h) steps. Quoting the fifth power for the finished answer overstates the accuracy by a factor of 1/h.
- Using RK4 on a stiff problem. The failure is not gradual: below the stability limit the answer is fine, above it the numerical solution oscillates with growing amplitude. If you find yourself needing a step far smaller than the features you care about, switch to an implicit solver.
- Integrating through a discontinuity. A switch in a forcing term inside a step destroys the order argument. Break the integration at the switch point.
- Assuming a smaller h is always better. Truncation error falls as h4, but rounding error accumulates roughly as 1/h. Past a certain point - typically when h is around 10−4 of the interval for double precision - the total error starts rising again.
- Forgetting the interval must divide by h. This calculator shortens the last step to land exactly on your final x and tells you when it has done so, but a hand-written loop will usually overshoot or stop short.
- Reading agreement between RK4 and Euler as confirmation. They agree when the solution is nearly straight over a step, which says nothing about whether either has resolved a feature between the steps.
Where RK4 fits among the solvers
RK4 is the middle of the family. Below it, the Euler method is the one-stage member, useful mainly for teaching and for showing what first-order accuracy costs. Above it sit the embedded pairs - Runge-Kutta-Fehlberg 4(5) and Dormand-Prince 5(4), the latter being MATLAB's ode45 and the default in most scientific libraries - which compute two solutions of different order from nearly the same stages and use their difference to steer an adaptive step size. That is the same Richardson idea this calculator uses, done inside each step instead of over the whole interval.
For stiff problems the explicit family is the wrong family entirely. Implicit schemes - backward Euler, the trapezoidal rule, and the backward differentiation formulas behind ode15s - solve an algebraic equation at each step, usually with a Newton or secant iteration inside the loop, and buy unconditional stability with that extra work.
When the equation happens to be linear and first order, you do not need a numerical scheme at all: the integrating factor method gives an exact answer. For constant-coefficient linear equations of any order the Laplace transform does the same. Numerical integration earns its keep on the equations those methods cannot touch - and it is worth solving one linear problem both ways to see how closely RK4 reproduces the exact answer.
A higher-order system is handled by the same code. Write y″ = g(x, y, y′) as the pair u′ = v, v′ = g(x, u, v) and apply RK4 to the vector (u, v), computing all four stages for both components before either is updated. This calculator handles the single-equation case; the arithmetic for a system is identical, component by component.