Numerical Methods: Euler to Runge–Kutta

Study Sheet

Numerical Methods: Euler to Runge–Kutta

Heun, RK4, error orders, and stability ceilings

Computing Solutions Honestly

Tip
The ladder of methods

EULER: follow the tangent, yn+1=yn+hf(tn,yn)y_{n+1} = y_n + hf(t_n, y_n) — global error O(h)O(h). IMPROVED EULER (Heun): predict with Euler, evaluate the slope at the prediction, use the AVERAGE of the two slopes — O(h2)O(h^2): on y′=yy' = y, y(0)=1y(0)=1, h=1h = 1 it gives 2.52.5 where Euler gives 22 and the truth is e≈2.718e \approx 2.718.

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

Tip
Stability: the other failure mode

Accuracy is not the only concern. On decay y′=−λyy' = -\lambda y, Euler iterates yn+1=(1−λh)yny_{n+1} = (1 - \lambda h)y_n: if h>2λh > \tfrac{2}{\lambda} the multiplier exceeds 11 in size and the computed solution OSCILLATES AND GROWS while the true one dies quietly.

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

Example
Worked: watching the error order with your own eyes

Solve y′=yy' = y, y(0)=1y(0) = 1 to t=1t = 1 (truth: e≈2.71828e \approx 2.71828).

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

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

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

Step 4 — Heun with h=0.25h = 0.25: error ≈0.022\approx 0.022, close to a QUARTER — O(h2)O(h^2) verified numerically. Running two step sizes and comparing is the standard self-check every real computation should include.

Side note
Significance: most equations are solved this way

Weather forecasts, spacecraft trajectories, protein folding, and car-crash simulations never see a closed-form solution — they are RK-style steppers running billions of steps. The stability theory of this topic is not academic: an unstable integrator once meant exploding simulations and, historically, real engineering failures.

Try it
Try it: one honest Euler step

For y′=t+yy' = t + y, y(0)=1y(0) = 1, take one Euler step with h=0.1h = 0.1, then one Heun step. Work: Euler: y1=1+0.1(0+1)=1.1y_1 = 1 + 0.1(0 + 1) = 1.1. Heun: predictor slope 11; corrector slope f(0.1,1.1)=1.2f(0.1, 1.1) = 1.2; y1=1+0.05(1+1.2)=1.11y_1 = 1 + 0.05(1 + 1.2) = 1.11.

True value: 2e0.1−0.1−1≈1.11032e^{0.1} - 0.1 - 1 \approx 1.1103 — Heun lands within 0.00030.0003 at a tenth the step cost of high-accuracy Euler.

Proofs & Why It Matters

Tip
Proof: Euler's global error is O(h)

One step from the true solution commits local error y(t+h)−y(t)−hy′(t)=h22y′′(ξ)y(t + h) - y(t) - hy'(t) = \tfrac{h^2}{2}y''(\xi) — the first Taylor term the tangent line drops.

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

Tip
Proof: Heun's extra order

Heun's update is y+h2[f(y)+f(y+hf(y))]y + \tfrac h2\left[f(y) + f(y + hf(y))\right]. Expand the second slope: f(y+hf(y))=f+hff′+O(h2)f(y + hf(y)) = f + hff' + O(h^2), so the update is y+hf+h22ff′+O(h3)y + hf + \tfrac{h^2}{2}ff' + O(h^3).

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

Going Deeper: Worked Problems

Example
Worked: RK4 by hand, one step

Apply one RK4 step with h=0.2h = 0.2 to y′=y−ty' = y - t, y(0)=2y(0) = 2 (exact solution y=t+1+ety = t + 1 + e^t).

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

Step 2 — weighted average: y1=2+0.26(2+2⋅2.1+2⋅2.11+2.222)=2+0.26(12.642)=2.4214y_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+e0.2=2.42140…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.42.4 (error 0.020.02).

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

Example
Worked: catching instability before it bites

For y′=−50yy' = -50y (a stiff decay), find the largest Euler step that stays stable, and compare with backward Euler.

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

Step 2 — backward Euler: yn+1=yn−50hyn+1y_{n+1} = y_n - 50hy_{n+1}, so yn+1=yn1+50hy_{n+1} = \dfrac{y_n}{1 + 50h}; the multiplier is below 11 for EVERY h>0h > 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.