Numerical Methods: Euler to Runge–Kutta

RK4 blends four slope samples per step for $O(h^4)$ — the workhorse default. Order buys compounding accuracy: halving $h$ improves Euler $2\times$ but RK4 $16\times$.

Stiff problems (mixing fast and slow decay) force explicit methods to tiny steps for stability alone — the reason implicit methods (backward Euler: $y_{n+1} = y_n + hf(t_{n+1}, y_{n+1})$, stable for every $h$) exist.

Step 1 — Euler with $h = 0.5$: $(1.5)^2 = 2.25$, error $0.468$.

Step 2 — Euler with $h = 0.25$: $(1.25)^4 \approx 2.441$, error $0.277$ — roughly halved, as $O(h)$ promises.

Step 3 — Heun with $h = 0.5$: multiplier $1 + h + \tfrac{h^2}{2} = 1.625$ per step, $(1.625)^2 \approx 2.641$, error $0.078$.

Step 4 — Heun with $h = 0.25$: multiplier $1.28125$ per step, $(1.28125)^4 \approx 2.695$, error $\approx 0.023$ — about $0.3$ of the previous error, approaching a quarter as $h$ shrinks — $O(h^2)$ verified numerically. Running two step sizes and comparing is the standard self-check every real computation should include.

Why it approximates: Taylor's theorem gives $y(t + h) = y(t) + hy'(t) + \tfrac{h^2}{2}y''(\xi)$, and $y'(t) = f(t, y(t))$, so an Euler step is the Taylor expansion with the $h^2$ term dropped. The dropped term bends toward the concave side: where solutions curve upward, Euler lands below them, and where they curve downward, above. Laying the computation out as a table — $t_n$, $y_n$, slope, next $y$ — keeps every step checkable.

Step 1 — the first step heads down with slope $-1$, to $0.5$.

Step 2 — at $(0.5, 0.5)$ the slope is $0$, so the second step is flat: Euler cannot see that the true solution, which turns around near $t = 0.69$, is already rising again.

Step 3 — from there the slopes are positive and each step climbs.

Step 4 — read the errors: they grow while the solution bends sharply (the first second), then shrink, because this equation pulls every solution toward the line $y = t - 1$ (the difference decays like $e^{-t}$), which also pulls the Euler values toward it. Whether errors accumulate or heal depends on the equation, not only on the method.

Step 1 — compute the radius: $x_{n+1}^2 + y_{n+1}^2 = (x_n + hy_n)^2 + (y_n - hx_n)^2 = (1 + h^2)(x_n^2 + y_n^2)$, because the cross terms $2hx_ny_n$ cancel.

Step 2 — so every step multiplies the radius by $\sqrt{1 + h^2} = \sqrt{1.25} \approx 1.118$.

Step 3 — after $12$ steps ($t = 6$, just short of one period $2\pi$) the radius is $1.25^6 \approx 3.81$ instead of $1$.

Step 4 — the reason is geometric: each step moves along the tangent to the circle, and a tangent segment always ends outside the circle. For a model of a frictionless system this is a qualitative failure — Euler adds energy every step — and a smaller $h$ only slows the spiral, it never removes it.

Why these weights: when $f$ depends on $t$ alone, the step is $\int_{t_n}^{t_n + h}f\,dt$ approximated by Simpson's rule, $\tfrac h6\big(f(t_n) + 4f(t_n + \tfrac h2) + f(t_n + h)\big)$, with the midpoint weight $4$ split between the two midpoint estimates $k_2$ and $k_3$. In general, the stages are arranged so that the step agrees with the Taylor expansion of the true solution through the $h^4$ term; the local error is $O(h^5)$, and the global error, after $\tfrac Th$ steps, is $O(h^4)$.

Step 1 — RK4's second step, from $(0.5, \tfrac{137}{192})$, gives $y(1) \approx 0.736342$ (exactly $\tfrac{54289}{73728}$): error $0.00058$.

Step 2 — Heun with $h = \tfrac14$ gives $\tfrac{390625}{524288} \approx 0.745058$: error $0.0093$.

Step 3 — Euler with $h = \tfrac18$ gives $\tfrac{5764801}{8388608} \approx 0.687218$: error $0.049$.

Step 4 — at equal cost, each step up the ladder buys roughly one more correct digit here, and the gap widens as the work increases: halving $h$ divides the errors by about $2$, $4$ and $16$. This is why RK4, not a finer Euler grid, is the default.

Step 1 — for the rotation $A$, the powers cycle ($A^2 = -I$), so the step matrix is $c\,I + s\,A$ with $c = 1 - \tfrac{h^2}{2} + \tfrac{h^4}{24}$ and $s = h - \tfrac{h^3}{6}$, and it scales every radius by $\sqrt{c^2 + s^2} = \sqrt{1 - \tfrac{h^6}{72} + \tfrac{h^8}{576}}$.

Step 2 — with $h = 0.5$: $c^2 + s^2 \approx 0.999790$, so the radius shrinks by the factor $0.999895$ per step.

Step 3 — after the same $12$ steps where Euler reached radius $3.81$, RK4's radius is $0.99874$: a drift of about $0.13\%$ instead of $281\%$, for four times the work per step.

Step 4 — note the sign: RK4 loses a trace of energy while Euler gains it. Neither conserves it exactly; long simulations of orbits use methods designed to (symplectic methods).

Reaching a fixed time $T$ takes $\tfrac Th$ steps; local errors accumulate (amplified at most by a bounded factor $e^{\emph{lt}}$ from the Lipschitz constant), totalling $\tfrac Th \cdot O(h^2) = O(h)$. The same bookkeeping with a sharper local error $O(h^{p+1})$ yields global order $p$ — which is exactly how Heun ($p = 2$) and RK4 ($p = 4$) earn their labels. $\blacksquare$

The true solution expands as $y + hy' + \tfrac{h^2}{2}y'' + O(h^3)$, and $y'' = \tfrac{d}{dt}f(y) = f'y' = ff'$ — the $h^2$ terms match. Local error drops to $O(h^3)$, hence global $O(h^2)$: averaging the two slopes secretly reproduces the second Taylor coefficient. $\blacksquare$

Step 1 — the four slopes: $k_1 = f(0, 2) = 2$; $k_2 = f(0.1, 2 + 0.1\cdot2) = f(0.1, 2.2) = 2.1$; $k_3 = f(0.1, 2 + 0.1\cdot2.1) = f(0.1, 2.21) = 2.11$; $k_4 = f(0.2, 2 + 0.2\cdot2.11) = f(0.2, 2.422) = 2.222$.

Step 2 — weighted average: $y_1 = 2 + \dfrac{0.2}{6}(2 + 2\cdot2.1 + 2\cdot2.11 + 2.222) = 2 + \dfrac{0.2}{6}(12.642) = 2.4214$.

Step 3 — exact: $y(0.2) = 1.2 + e^{0.2} = 2.42140\ldots$ — RK4 agrees to about five significant figures in one step; Euler would have given $2.4$ (error $0.02$).

Step 4 — the pattern: slopes at start, two midpoint estimates, and end, weighted $1:2:2:1$ — Simpson's rule in disguise, which is why the order is four.

Step 1 — Euler: $y_{n+1} = (1 - 50h)y_n$; stability needs $|1 - 50h| < 1$, i.e. $h < 0.04$. With $h = 0.05$ the multiplier is $-1.5$: the numbers alternate sign and grow $50\%$ per step — nonsense for a decaying solution.

Step 2 — backward Euler: $y_{n+1} = y_n - 50hy_{n+1}$, so $y_{n+1} = \dfrac{y_n}{1 + 50h}$; the multiplier is below $1$ for every $h > 0$ — unconditionally stable.

Step 3 — the trade: backward Euler costs an algebraic solve per step (trivial here, a Newton iteration in general) but lets the step size follow accuracy instead of stability.

Step 4 — this is why "stiff" is a word every simulation engineer knows: the fast-decaying mode is boring but still dictates explicit step sizes.