10.2 ODE solvers: Euler, midpoint, RK4

An ODE $\dot x = f(t, x)$ with initial condition $x(t_0) = x_0$ has a unique solution (for reasonable $f$ — see Foundations 5.4). Numerically integrating it means walking forward from $t_0$ in small steps of size $h$ , producing a sequence $x_0, x_1, x_2, \ldots$ of approximate values at times $t_0, t_0 + h, t_0 + 2h, \ldots$ . The recipe for going from $x_n$ to $x_{n+1}$ is what distinguishes the methods.

This lesson covers the three solvers that span the practical landscape of explicit ODE integrators: forward Euler (first-order, simplest possible), midpoint (second-order, modestly more work), and Runge–Kutta 4 (fourth-order, the workhorse of practical computation). RK4 is the method behind every ODE-integration interactive on this site.

Forward Euler

The simplest possible recipe: approximate the derivative by its current value, treat it as constant over the step:

x_{n+1} \;=\; x_n \;+\; h\, f(t_n, x_n).

This is exactly the forward-difference approximation from 10.1, rearranged. By the Taylor analysis there, the local truncation error per step is $\mathcal{O}(h^2)$ ; after $N \sim 1/h$ steps the accumulated global error is $\mathcal{O}(h)$ . Euler is a first-order method.

Euler is cheap (one evaluation of $f$ per step), trivially easy to code, and almost never the right choice. The reason is stability rather than accuracy. For the test equation $\dot x = -\alpha x$ with $\alpha > 0$ (exponential decay), the Euler update is $x_{n+1} = (1 - \alpha h)\, x_n$ . The numerical solution decays only if $|1 - \alpha h| < 1$ , i.e. $h < 2 / \alpha$ . Beyond that threshold the numerical solution blows up, even though the true solution decays gracefully. For stiff problems (large $\alpha$ ) the stability bound forces tiny step sizes regardless of how accurate you actually need.

Midpoint (RK2)

Midpoint trades one extra $f$ -evaluation per step for an order of accuracy. The recipe:

\begin{aligned} k_1 &\;=\; f(t_n, x_n), \\ k_2 &\;=\; f(t_n + h/2,\; x_n + (h/2)\, k_1), \\ x_{n+1} &\;=\; x_n + h\, k_2. \end{aligned}

The idea: $k_1$ is the slope at the start of the step; using it, you take a half-step to find a better-informed slope $k_2$ at the midpoint; that slope is the one you use for the full step. The local truncation error is $\mathcal{O}(h^3)$ and the global error is $\mathcal{O}(h^2)$ . Second-order method, costing twice what Euler costs per step.

Runge–Kutta 4

The workhorse of explicit ODE integration. Four function evaluations per step, fourth-order accuracy globally:

\begin{aligned} k_1 &\;=\; f(t_n, x_n), \\ k_2 &\;=\; f(t_n + h/2,\; x_n + (h/2)\, k_1), \\ k_3 &\;=\; f(t_n + h/2,\; x_n + (h/2)\, k_2), \\ k_4 &\;=\; f(t_n + h,\; x_n + h\, k_3), \\ x_{n+1} &\;=\; x_n + \frac{h}{6}\, \bigl(k_1 + 2 k_2 + 2 k_3 + k_4\bigr). \end{aligned}

The slopes are evaluated at the start, twice at the midpoint (each using the previous slope), and once at the end. The four slopes are combined in a weighted average with weights $(1/6, 1/3, 1/3, 1/6)$ . The motivation is to make the Taylor expansion of the update match the true expansion of the solution to fourth order. The resulting method has local error $\mathcal{O}(h^5)$ and global error $\mathcal{O}(h^4)$ .

Going from second-order to fourth-order at a cost of 4× the function evaluations rather than 2× sounds expensive, but the gain is enormous: for the same accuracy, RK4 can use a much larger step size. On reasonable problems RK4 with $h = 0.1$ outperforms Euler with $h = 0.001$ .

Worked example

▶ Worked example: one Euler and one RK4 step, by hand

The problem. Take one step of size $h = 0.5$ for the ODE $\dot x = -x$ from initial condition $x(0) = 1$ . The true solution is $x(t) = e^{-t}$ , so $x(0.5) = e^{-0.5} \approx 0.6065$ .

Euler step.

x_1 \;=\; x_0 + h\, f(t_0, x_0) \;=\; 1 + 0.5 \cdot (-1) \;=\; 0.5.

Error: $|0.5 - 0.6065| = 0.1065$ .

RK4 step. With $f(t, x) = -x$ (autonomous; doesn’t depend on $t$ ):

\begin{aligned} k_1 &\;=\; f(0, 1) \;=\; -1. \\ k_2 &\;=\; f(0.25,\; 1 + 0.25 \cdot (-1)) \;=\; f(0.25, 0.75) \;=\; -0.75. \\ k_3 &\;=\; f(0.25,\; 1 + 0.25 \cdot (-0.75)) \;=\; f(0.25, 0.8125) \;=\; -0.8125. \\ k_4 &\;=\; f(0.5,\; 1 + 0.5 \cdot (-0.8125)) \;=\; f(0.5, 0.59375) \;=\; -0.59375. \end{aligned}

x_1 \;=\; 1 + \frac{0.5}{6} \bigl(-1 + 2(-0.75) + 2(-0.8125) + (-0.59375)\bigr) \;=\; 1 + \frac{0.5}{6} (-4.71875) \;=\; 1 - 0.3932... \;=\; 0.60677...

Error: $|0.60677 - 0.6065| \approx 3 \times 10^{-4}$ .

Comparison. One RK4 step (4 function evaluations) gives error $\sim 10^{-4}$ . The same accuracy with Euler would need step size $h \sim 10^{-4}/0.1 \sim 10^{-3}$ , so $\sim 500$ Euler steps. RK4 is more than $100 \times$ more efficient on this problem.

Side-by-side comparison

step size h = 0.30 number of steps N = 40

Three solvers integrate the harmonic-oscillator system (ẋ, v̇) = (v, −x) from the initial condition (x, v) = (1, 0). The true trajectory is the unit circle (dashed). Each method takes the same number of steps with the same step size h: Euler spirals visibly outward at every step (it is consistently *adding* energy), midpoint drifts much less because its O(h²) error is smaller, and RK4 has O(h⁴) error and tracks the true circle even at h = 0.6. Crank h up and the disparity becomes dramatic: at h ≈ 0.4 Euler is already noticeably wrong, while RK4 is essentially exact.

Three solvers integrate the harmonic-oscillator system from $(x, v) = (1, 0)$ . The true trajectory is the unit circle (energy conserved). Euler spirals outward at every step — it consistently adds energy to the system. Midpoint drifts much more slowly. RK4 tracks the circle essentially perfectly, even at large step sizes. Crank $h$ up to see the disparity become dramatic.

Stability

Stability is the second axis of solver quality, distinct from accuracy. A solver is stable for a given ODE and step size if numerical solutions stay bounded when the true solution does. For the test equation $\dot x = \lambda x$ (linear, scalar, possibly complex $\lambda$ ), each method has a stability region in the complex plane of $z = \lambda h$ :

Forward Euler is stable for $|1 + z| < 1$ — a disc of radius 1 centred at $-1$ . Practically, $\mathrm{Re}\,\lambda < 0$ requires $|h| < 2 / |\lambda|$ .
RK4 is stable in a larger region that includes more of the negative real axis: stable for $|h \lambda| \lesssim 2.78$ on the real axis.
Implicit methods (backward Euler, trapezoidal rule) have unbounded stability regions — they are stable for any $h > 0$ on $\mathrm{Re}\, \lambda < 0$ . The cost is solving a nonlinear equation at each step rather than just evaluating $f$ .

The implication: for stiff systems — those with widely-separated time scales — explicit methods like RK4 are forced into tiny step sizes by the fastest time scale, even when you only care about the slow dynamics. Stiff problems demand implicit solvers, which is where the heavy machinery of numerical-ODE-software libraries lives (LSODA, SUNDIALS, etc.).

What we use this for

RK4 (sometimes its adaptive-step cousin Cash–Karp or Dormand–Prince) is the integrator behind every ODE simulation on this site:

DirectionField (MATH-30) in Foundations 5.1 — RK4 for the trajectories that thread through the slope field.
FixedPointGallery (MATH-29) in Foundations 5.4 — RK4 for the phase-portrait trajectories around each equilibrium type.
DampedDrivenOscillator (MATH-9) — RK4 for the time-domain response curves.

The next lesson, 10.3, takes the same finite-difference machinery from 10.1 into spatial dimensions and discretises the wave and heat equations.