Runge-Kutta Methods

The initial-value problem

The setting is an ordinary differential equation written as a first-order system:

dy/dx = f(x, y)        with y(x₀) = y₀

You know the slope f at any point (x, y) in the plane and the value of y at the starting x₀. Your job is to march forward in x, computing y at successive points x₁, x₂, …. The challenge is that the slope changes as you move — so the longer your step, the worse it is to assume the slope stays constant.

Most physical laws come in this form once you remove the time derivative: Newton's second law gives second-order ODEs that decompose into first-order systems by introducing velocity as an extra variable; the Schrödinger equation, the heat equation, every chemical-kinetics rate equation, every radioactive-decay law lives here. The cases that admit closed-form solutions (linear with constant coefficients, separable, exact, certain second-order linear) are a thin slice of the universe of equations engineers actually face. The rest demand numerical integrators, and Runge-Kutta methods are the most-used family on Earth.

Why Euler's method is not enough

Euler's method (RK1) is the simplest possible scheme: assume the slope at the start of the interval holds across the whole interval.

y_{n+1} = y_n + h · f(x_n, y_n)

It is first-order accurate: global error scales as h. Halve the step size, halve the error. To get four-digit accuracy you typically need millions of steps, and the work per step (one f-evaluation) is the cheapest possible — but the iteration count outweighs the per-step savings on any real problem. Worse, Euler is unstable on oscillatory problems: simulating a pendulum with Euler causes the swing to grow exponentially even when the equation conserves energy.

The fix is to sample the slope at more than one point within each step. RK2 (Heun's method, the improved Euler) takes two samples — at the start and at the Euler-predicted endpoint — and averages them. The average is a better representative slope across the step, and the global error drops to O(h²). Going further, RK4 takes four samples and weights them carefully to cancel error terms up to and including h⁴.

The classical RK4 formula

One step of RK4 from (x_n, y_n) to (x_n + h, y_{n+1}):

k₁ = f(x_n,         y_n)
k₂ = f(x_n + h/2,   y_n + (h/2) · k₁)
k₃ = f(x_n + h/2,   y_n + (h/2) · k₂)
k₄ = f(x_n + h,     y_n + h · k₃)

y_{n+1} = y_n + (h/6) · (k₁ + 2 k₂ + 2 k₃ + k₄)

Reading the formula: k₁ is the slope at the start of the interval. k₂ is the slope at the midpoint, predicted using k₁. k₃ is the slope at the midpoint again, predicted using the better k₂ estimate. k₄ is the slope at the end, predicted using k₃. The final step is the weighted average (1, 2, 2, 1)/6 — the same Simpson's-rule pattern that integrates a polynomial of degree four exactly.

The local truncation error per step is O(h⁵) and the accumulated global error after N = (b − a)/h steps is O(h⁴). That fourth-power scaling is the entire reason RK4 dominates: at h = 0.01, the error is ~10⁻⁸; at h = 0.001 it's ~10⁻¹², approaching double precision. Halve h and the error drops by 16; cut h by ten and the error drops by ten thousand.

Worked example: dy/dx = y, y(0) = 1

This problem has the closed-form solution y(x) = eˣ, so we can score the numerical method against the truth. Pick h = 0.1 and integrate from x = 0 to x = 1, comparing Euler to RK4.

True value: y(1) = e = 2.71828 18284 59045 ...

Euler (RK1, 10 steps of h=0.1):
  y_0 = 1.0
  y_1 = 1.0 + 0.1 · 1.0           = 1.10000
  y_2 = 1.1 + 0.1 · 1.1           = 1.21000
  ...
  y_10                              = 2.59374     error ≈ 0.124 (4.6 %)

RK4 (10 steps of h=0.1):
  At x = 0:
    k₁ = f(0, 1)         = 1.000000
    k₂ = f(0.05, 1.05)   = 1.050000
    k₃ = f(0.05, 1.0525) = 1.052500
    k₄ = f(0.10, 1.10525) = 1.105250
    y_1 = 1 + (0.1/6)(1 + 2·1.05 + 2·1.0525 + 1.10525)
        = 1 + (0.1/6)(6.31025) = 1.105170833...
  ...
  y_10                              = 2.71828 1235...   error ≈ 5.6 × 10⁻⁶

Euler gets one digit; RK4 gets six. Both algorithms used the same step size and the same problem. The difference is just the four-stage averaging. Halving h to 0.05 (twenty steps) shrinks Euler's error to 0.067 (a factor of 1.85, roughly halving as expected) and RK4's error to about 3.5 × 10⁻⁷ (a factor of 16, the fourth-power scaling).

This is why high-precision applications stretch RK4 over very small h: the cost grows linearly with the step count but the error shrinks like h⁴. There is a sweet spot where roundoff (~10⁻¹⁵ in double precision) starts to dominate truncation error, beyond which decreasing h further only adds noise.

Comparing methods at fixed work

Method	Stages per step	Order	Global error scaling	Steps for 6-digit accuracy on dy/dx=y, y(0)=1, x∈[0,1]
Euler (RK1)	1	1	O(h)	~270,000
Improved Euler (Heun, RK2)	2	2	O(h²)	~1,000
RK3 (Kutta)	3	3	O(h³)	~70
RK4 (classical)	4	4	O(h⁴)	~10
RK6 (Verner)	~7	6	O(h⁶)	~4
Adams-Bashforth 4 (multistep)	1 (after start)	4	O(h⁴)	~10 + start-up
Bulirsch-Stoer (extrapolation)	variable	up to ~10	O(h^k) adaptive	~3–5

The pattern is sharp: each extra stage typically buys an order of accuracy until you hit the Butcher-barrier ceiling (RK4 needs 4 stages, RK5 needs 6, RK6 needs 7, RK8 needs 11). The cost-per-step grows linearly with stage count but the error drops geometrically with order, so increasing order is almost always a win on smooth problems. Most production solvers default to order 4 or 5 because that combination is robust across a wide range of right-hand-sides.

Adaptive step size: RKF45 and DOPRI5

Fixed step size is wasteful: in regions where the solution changes slowly, you can take huge steps; near a sharp transient you need tiny ones. Adaptive Runge-Kutta methods estimate the local error after each step and resize h accordingly.

The trick is to compute two solutions per step at different orders, sharing the same slope evaluations. RKF45 (Fehlberg, 1969) uses six stages to produce both a fourth-order and a fifth-order estimate. The difference between them is a tight estimate of the local truncation error of the fourth-order step. If that error exceeds the user's tolerance, halve h and retry; if it is well under, double h and proceed.

The current default in most production solvers is DOPRI5 (Dormand and Prince, 1980), a seven-stage method whose first stage of the next step equals the last stage of the current step (the FSAL property — First Same As Last). FSAL effectively makes DOPRI5 cost six stages per step on average, identical to RKF45, but with smaller error constants. Matlab's ode45, SciPy's solve_ivp(method='RK45'), Julia's DifferentialEquations.jl default — all DOPRI5.

Where Runge-Kutta shows up

• Spacecraft trajectory propagation. Voyager, Cassini, the James Webb Space Telescope and SpaceX Dragon all integrate orbital state vectors with high-order Runge-Kutta or related methods. NASA's GMAT and JPL's SPICE toolkits use RK7(8) and RK8(9) routines for millimetre-level accuracy over months of flight.
• Climate and weather models. The European Centre for Medium-Range Weather Forecasts (ECMWF) integrates the primitive equations of the atmosphere with semi-implicit Runge-Kutta variants on a 9 km grid. Each forecast advances ~10⁹ grid points 10 days into the future in a few hours of supercomputer time.
• Molecular dynamics. GROMACS, NAMD, and AMBER simulate proteins by integrating Newton's equations for tens of thousands of atoms. They use velocity-Verlet (a symplectic relative of RK methods) for energy conservation, with RK extensions for thermostats and barostats.
• Game-engine physics. Unity, Unreal, Box2D and Bullet integrate rigid-body and soft-body dynamics with semi-implicit Euler or RK2 — fourth-order is overkill at 60 Hz for visual-quality physics, but the framework is identical.
• Pharmacokinetic modelling. NONMEM and Pumas simulate drug concentration over time with RK4 or DOPRI5 inside their compartment-model engines, used at every major pharmaceutical company for FDA submissions.

Variants and extensions

• Runge-Kutta-Fehlberg (RKF45). The classic adaptive variant — fourth- and fifth-order solutions share six stages, the difference giving a step-size error estimate. Still widely taught even though Dormand-Prince has largely replaced it in production.
• Dormand-Prince (DOPRI5, DOPRI8). Modern adaptive workhorses with the FSAL property and tighter error constants than Fehlberg. DOPRI5 is the default in Matlab ode45 and SciPy solve_ivp.
• Implicit Runge-Kutta (Radau, Lobatto, Gauss). The slope evaluations include y_{n+1} on the right-hand side, so each step solves a nonlinear system. Required for stiff ODEs where explicit RK is forced into impractically tiny step sizes.
• Symplectic Runge-Kutta. Variants designed to preserve the Hamiltonian structure of mechanical systems. Used in long-term solar-system simulations where energy must not drift over millions of years.
• Stochastic Runge-Kutta. Extensions to stochastic differential equations dy = f dt + g dW, important in mathematical finance (Black-Scholes path simulation) and statistical physics (Langevin dynamics).

Common pitfalls

• Using RK4 on a stiff problem. Symptoms: the solver takes microscopic steps and crawls. Cause: a fast-decaying mode forces step size below ~1/|λ_max|, where λ_max is the largest eigenvalue magnitude of the Jacobian. Switch to an implicit BDF or Radau method.
• Energy drift in long Hamiltonian simulations. RK4 isn't symplectic. Over millions of orbital periods it slowly heats or cools the system. Use leapfrog or a higher-order symplectic integrator if you need million-year stability.
• Mismatching tolerances and step caps. Modern adaptive solvers expose rtol, atol, hmin, hmax. Setting hmax too small forces tiny steps even where they aren't needed; setting tolerances too loose lets the solver overshoot transients. Tune them to your problem.
• Order vs accuracy confusion. A fourth-order method only achieves O(h⁴) error asymptotically as h → 0. At large h or near discontinuities, RK4 can be less accurate than RK2. Check that you're in the asymptotic regime by halving h and watching the error scale.
• Discontinuities in f. All Runge-Kutta methods assume f is smooth across the step. If f has a kink (e.g. friction switching sign, a control law toggling), the high-order error estimates lose meaning. Detect events explicitly and restart the integrator after each event.

Why RK4 became the default

RK4 is not the most accurate method, the cheapest method, the most stable method, or the most modern method. It is the most balanced: fourth-order accuracy at four function evaluations is the highest order achievable per stage before the Butcher barriers. It generalises trivially to vector-valued ODEs, requires only the right-hand side f (no Jacobian, no linear solver), and works straight out of the box on any non-stiff problem. The same five lines of code that integrate a pendulum integrate the orbit of Mars or the population dynamics of a coral reef.

The hierarchy of go-to ODE solvers in 2026 looks like: DOPRI5 (adaptive RK4/5) for most non-stiff problems; BDF or implicit Runge-Kutta for stiff problems; symplectic integrators for long Hamiltonian runs; and high-order RK or Bulirsch-Stoer for tight-tolerance work. RK4 itself, in its 1901 form, remains the textbook starting point and the right tool for any quick-and-dirty ODE in a script — which is why it gets called "the workhorse" in every numerical-analysis course taught since 1960.