Euler, Improved Euler, and Runge–Kutta

Analytic techniques (integrating factors, separation of variables, series, Laplace transforms) solve a limited catalogue of equations. A generic nonlinear first-order problem

\frac{d y}{d t} = f (t, y), y (t_{0}) = y_{0}

usually has no solution expressible in elementary functions, yet the existence theorem still guarantees a unique solution $y = ϕ (t)$ whenever $f$ and $f_{y}$ are continuous near $(t_{0}, y_{0})$ . A numerical method abandons the search for a formula and instead produces a table of approximate values $y_{1}, y_{2}, \dots, y_{N}$ at prescribed points $t_{1}, t_{2}, \dots, t_{N}$ , each meant to approximate $ϕ (t_{n})$ .¹ We fix a uniform step size $h$ , set $t_{n} = t_{0} + nh$ , and abbreviate $f_{n} = f (t_{n}, y_{n})$ . Each method is a rule for advancing from $y_{n}$ to $y_{n + 1}$ .

Euler's method

The oldest and simplest rule comes from three equivalent readings of the same approximation.

Tangent line. At $(t_{n}, y_{n})$ the differential equation gives the slope $f_{n}$ . Follow that tangent for a horizontal distance $h$ : the height changes by $h f_{n}$ .
Difference quotient. Replace $ϕ^{'} (t_{n})$ by the forward difference $(ϕ (t_{n + 1}) - ϕ (t_{n})) / h$ in the equation and solve for $ϕ (t_{n + 1})$ .
Truncated Taylor series. Expand $ϕ (t_{n} + h) = ϕ (t_{n}) + ϕ^{'} (t_{n}) h + \frac{1}{2} ϕ^{''} (\overset{ˉ}{t}_{n}) h^{2}$ and drop the quadratic term.

All three yield the Euler formula.

Geometrically the exact solution curve is replaced by a polygon of tangent segments. For a convex solution ( $ϕ^{''} > 0$ ) every segment lies below the curve, so the approximation drifts steadily downward; a concave solution drifts the other way.

Euler's method threads a chain of tangent segments through the slope field; for a convex solution each segment undershoots and the error

E_{n} = ϕ (t_{n}) - y_{n}

grows step by step.

Written as a program the method is a single loop; only line 4 encodes the choice of rule, and the methods that follow replace that single line.

Algorithm:

\textsc{Euler}(f, t_0, y_0, h, N)

— tangent-line integration on

N

steps

1
$t \gets t_0$ , $y \gets y_0$
2
output $(t, y)$
3
for $n$ from $1$ to $N$ do
4
$y \gets y + h \cdot f(t, y)$
advance along the tangent
5
$t \gets t + h$
6
output $(t, y)$
7
end

A test problem

The linear equation

y^{'} = 1 - t + 4 y, y (0) = 1

has the exact solution $ϕ (t) = \frac{1}{4} t - \frac{3}{16} + \frac{19}{16} e^{4 t}$ , which lets us monitor the error exactly. Its solutions diverge quickly (the $e^{4 t}$ term), so it is deliberately hard to track over a long interval, and differences between methods show clearly.² Euler's method with several step sizes gives the following values, with the percentage error against $ϕ$ shown in parentheses.

$t$	$h = 0.05$	$h = 0.01$	$h = 0.001$	Exact $ϕ (t)$
$0.5$	$7.2902$	$8.3767 (3.8%)$	$8.6771 (0.40%)$	$8.7120$
$1.0$	$45.588$	$60.037 (7.5%)$	$64.383 (0.79%)$	$64.898$
$2.0$	$1745.7$	$3029.3 (14%)$	$3484.2 (1.6%)$	$3540.2$

Reaching $1.6%$ accuracy at $t = 2$ requires $2000$ steps. The accuracy is unimpressive because the error compounds: halving $h$ only halves the error, and the error grows with $t$ as the solution steepens.

A single step is pure arithmetic. Taking the coarser $h = 0.1$ , so the numbers stay legible, and running five steps by hand:

Truncation error and order

Two errors matter, and they must be kept separate.

Local error

e

is the gap a single step opens from a point on the exact curve; global error

E

is the larger accumulated gap of the full polygon, which also carries forward the error already made.

Subtracting the Euler step from the exact Taylor expansion isolates the local error. Since the two share the terms $ϕ (t_{n}) + h f_{n}$ , only the quadratic remainder remains:

e_{n + 1} = ϕ (t_{n + 1}) - y_{n + 1}^{*} = \frac{1}{2} ϕ^{''} (\overset{ˉ}{t}_{n}) h^{2}, t_{n} < \overset{ˉ}{t}_{n} < t_{n + 1} .

The local truncation error is proportional to $h^{2}$ ; a uniform bound is $∣ e_{n} ∣ \leq \frac{1}{2} M h^{2}$ with $M = max ∣ ϕ^{''} ∣$ on the interval.³ The global error accumulates one local error per step over $\sim (T - t_{0}) / h$ steps, which drops one power of $h$ : it can be shown that $∣ E_{n} ∣ \leq K h$ on any finite interval.⁴ This scaling defines the method's order.

For the test problem $ϕ^{''} (t) = 19 e^{4 t}$ , so the local error near $t = 2$ is about $e^{8} / e^{0} \approx 2500$ times larger than near $t = 0$ . A uniform step sized for the hard end wastes effort at the easy end, which motivates adaptive codes that shrink $h$ only where $ϕ^{''}$ is large.

The backward Euler formula

Approximating the same integral by the value at the right endpoint gives the backward Euler formula

y_{n + 1} = y_{n} + h f (t_{n + 1}, y_{n + 1}) .

The unknown $y_{n + 1}$ now appears on both sides, so the formula is implicit: each step requires solving an equation. For the linear test problem this is routine algebra, but for nonlinear $f$ it needs an iterative solve. Backward Euler is no more accurate than forward Euler (still order 1) and more expensive; its advantage is stability, which matters for stiff problems, where an explicit method forces an impractically small step.

The improved Euler method

Integrating the equation over one step turns the problem into approximating an area:

ϕ (t_{n + 1}) = ϕ (t_{n}) + t_{n} \int t_{n + 1} f (t, ϕ (t)) d t .

Euler replaces the integrand by its left-endpoint value $f_{n}$ : the area of a rectangle. A better estimate averages the endpoint slopes, replacing the rectangle by a trapezoid.

Euler approximates the step integral by a left-endpoint rectangle (left); the improved Euler method uses a trapezoid averaging the two endpoint slopes (right), matching the curve far more closely.

The trapezoid rule uses $\frac{1}{2} (f (t_{n}, y_{n}) + f (t_{n + 1}, y_{n + 1}))$ , but $y_{n + 1}$ is unknown. Replacing it by the Euler prediction removes the implicitness and gives an explicit two-stage rule.

The cost is two evaluations of $f$ per step instead of one.⁵ For a fixed $h$ the improved Euler method does twice the work of Euler, or equivalently matches Euler at half the step.

On the test problem the improved Euler method with $h = 0.025$ ( $160$ evaluations) is already slightly more accurate than Euler with $h = 0.001$ ( $2000$ evaluations): about one-twelfth the work for comparable accuracy.

The Runge–Kutta method

Euler and improved Euler are the first two members of the Runge–Kutta family, which sample $f$ at several points inside the step and take a weighted average of the slopes. The classical fourth-order, four-stage method is the standard general-purpose integrator.⁶

Reading the four slopes: $k_{n 1}$ is the slope at the left end; $k_{n 2}$ is a midpoint slope reached by an Euler half-step; $k_{n 3}$ is a second, corrected midpoint slope; $k_{n 4}$ is the slope at the right end reached using $k_{n 3}$ . The weighting counts the two midpoint slopes double.

RK4 samples four slopes across one step — one at the left, two at the midpoint, one at the right — and advances by their

1 : 2 : 2 : 1

weighted average.

When $f$ depends on $t$ alone, the four slopes collapse to Simpson's rule for $\int f d t$ , whose $h^{5}$ error is consistent with RK4's order. On the test problem RK4 with $h = 0.2$ ( $40$ evaluations) matches the accuracy of improved Euler with $h = 0.025$ ( $160$ evaluations) at a quarter of the work ( $1.4%$ versus $1.2%$ error at $t = 2$ ); with $h = 0.05$ the value at $t = 2$ is correct to four digits ( $0.009%$ error).

Comparing the methods

Each added stage raises the order of accuracy and requires more evaluations of $f$ per step. Since evaluating $f$ dominates the run time, the fair comparison holds the total number of evaluations fixed.

Method	Evaluations of $f$ per step	Local error $e_{n}$	Global error $E_{n}$	Order $p$
Euler	$1$	$\propto h^{2}$	$\propto h$	$1$
Backward Euler	$1$ (implicit)	$\propto h^{2}$	$\propto h$	$1$
Improved Euler (Heun)	$2$	$\propto h^{3}$	$\propto h^{2}$	$2$
Runge–Kutta (RK4)	$4$	$\propto h^{5}$	$\propto h^{4}$	$4$

The test problem, evaluated at $t = 2$ with roughly equal computational effort, shows the ranking.

Method	Step $h$	Evaluations to $t = 2$	Value at $t = 2$	Error
Euler	$0.001$	$2000$	$3484.2$	$1.6%$
Improved Euler	$0.025$	$160$	$3496.7$	$1.2%$
Runge–Kutta	$0.05$	$160$	$3539.9$	$0.009%$

At the same $160$ evaluations RK4 is more than a hundred times more accurate than improved Euler. A higher-order method is more efficient: for the same work it delivers better accuracy, and for the same accuracy it needs far less work. On a log–log plot of error against step size, a method of order $p$ produces a straight line of slope $p$ , so halving $h$ divides the error by $2^{p}$ .

Global error versus step size on log–log axes. A method of order

p

traces a line of slope

p

, so RK4's error falls off far faster than Euler's as the step shrinks.

Step size, accuracy, and adaptive methods

A method's order fixes how fast error falls as $h \to 0$ , but it does not fix the right $h$ for a given tolerance. Since the local error is proportional to $h^{p + 1}$ , the step needed for a target tolerance $ϵ$ scales like $ϵ^{1/ (p + 1)}$ , and the constant depends on derivatives of the unknown solution. Practical codes estimate the local error as they run: compute the step two ways (for instance with a $p$ th- and a $(p + 1)$ th-order formula) and use the difference as an error estimate. When the estimate exceeds the tolerance, shrink $h$ and retry; when it is comfortably under, grow $h$ to save work.

Embedding a higher-order formula that reuses the same slope evaluations adds almost no extra work. The Runge–Kutta–Fehlberg pair (RKF45) computes a fourth- and a fifth-order estimate from one shared set of stages plus one extra, and its 1977 Fortran implementation RKF45 became a template for the adaptive integrators in modern software.⁷

Boyce, §8.1 — a numerical method approximates the solution of $y^{'} = f (t, y)$ , $y (t_{0}) = y_{0}$ on a mesh $t_{n} = t_{0} + nh$ rather than seeking a closed form. ↩
Boyce, §8.1, Example 1 — the linear problem $y^{'} = 1 - t + 4 y$ , $y (0) = 1$ with exact solution $ϕ (t) = \frac{1}{4} t - \frac{3}{16} + \frac{19}{16} e^{4 t}$ , used throughout Ch. 8 to compare methods. ↩
Boyce, §8.1 — the Euler local truncation error $e_{n + 1} = \frac{1}{2} ϕ^{''} (\overset{ˉ}{t}_{n}) h^{2}$ , bounded by $\frac{1}{2} M h^{2}$ . ↩
Boyce, §8.1, Problem 20 — the global truncation error bound $∣ E_{n} ∣ \leq K h$ for the Euler method on a finite interval, making it a first-order method. ↩
Boyce, §8.2 — the improved Euler (Heun) formula as a two-stage predictor–corrector with local error proportional to $h^{3}$ ; also Simmons, §75. ↩
Boyce, §8.3 — the classical fourth-order four-stage Runge–Kutta formula, its four-slope weighted average, and the Simpson's-rule reduction when $f = f (t)$ ; also Simmons, §76. ↩
Boyce, §8.3 — adaptive Runge–Kutta methods and the Runge–Kutta–Fehlberg (RKF45) embedded fourth/fifth-order pair. ↩

$n$	$t_{n}$	$y_{n}$	$f_{n} = 1 - t_{n} + 4 y_{n}$	$y_{n + 1} = y_{n} + 0.1 f_{n}$
$0$	$0.0$	$1.0000$	$5.0000$	$1.5000$
$1$	$0.1$	$1.5000$	$6.9000$	$2.1900$
$2$	$0.2$	$2.1900$	$9.5600$	$3.1460$
$3$	$0.3$	$3.1460$	$13.2840$	$4.4744$
$4$	$0.4$	$4.4744$	$18.4976$	$6.3242$