Implementing ODE solvers from scratch — Euler, RK4

WHAT problem are we solving?

WHY discretize? Most real $f(t,y)$ have no closed-form integral. But the slope $f$ is always computable. So we walk forward in tiny straight-line hops.

Euler's Method — deriving from first principles

HOW do we get from $y_n$ to $y_{n+1}$? Use the definition of a derivative: $y'(t_n) = \lim_{h\to 0}\frac{y(t_n+h)-y(t_n)}{h} \approx \frac{y_{n+1}-y_n}{h}$

Drop the limit (use a finite $h$) and solve for $y_{n+1}$:

Equivalently from a Taylor expansion: $y(t_n+h) = y(t_n) + h\,y'(t_n) + \underbrace{\tfrac{h^2}{2}y''(\xi)}_{\text{error we throw away}}$ The discarded term is $O(h^2)$ per step → local truncation error $O(h^2)$. Over $N = T/h$ steps the errors accumulate to global error $==O(h)==$. Halve $h$ → halve the error. That's slow.

def euler(f, t0, y0, h, n):
    t, y = t0, y0
    ts, ys = [t], [y]
    for _ in range(n):
        y = y + h * f(t, y)   # one slope, whole step
        t = t + h
        ts.append(t); ys.append(y)
    return ts, ys

RK4 — averaging four slopes

WHY do better? Euler uses the slope only at the left of the interval, so it consistently over/undershoots when the curve bends. RK4 samples the slope at the start, the middle (twice), and the end, then takes a weighted average — analogous to Simpson's rule for the integral $\int_{t_n}^{t_n+h} f\,dt$.

Why the weights $1,2,2,1$? Sum $= 6$, so $\frac h6(\dots)$ is a true weighted mean slope. RK4 uses two separate midpoint slopes ($k_2, k_3$), each given weight $2$, because the midpoint best represents the average behaviour over $[t_n, t_n+h]$. This is analogous to Simpson's $1/3$ rule (which weights endpoint–midpoint–endpoint as $1,4,1$): both put extra weight on the middle of the interval. RK4 splits Simpson's single "$4$" into two midpoint evaluations $k_2,k_3$ that together carry weight $2+2=4$. Matching the Taylor series up to $h^4$ forces exactly these coefficients.

• Local truncation error $O(h^5)$, global error $==O(h^4)==$.
• Halve $h$ → error drops by $2^4 = 16\times$. That's why RK4 is the workhorse.

def rk4(f, t0, y0, h, n):
    t, y = t0, y0
    ts, ys = [t], [y]
    for _ in range(n):
        k1 = f(t,         y)
        k2 = f(t + h/2,   y + h/2 * k1)
        k3 = f(t + h/2,   y + h/2 * k2)
        k4 = f(t + h,     y + h   * k3)
        y = y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
        t = t + h
        ts.append(t); ys.append(y)
    return ts, ys

Worked Example 1 — $y' = y$, $y(0)=1$ (true: $e^t$)

Take $h = 1$, one step, find $y_1$.

Euler: $y_1 = y_0 + h f(t_0,y_0) = 1 + 1\cdot 1 = 2$. Why? One left-slope step. True value $e^1 \approx 2.718$ → error $0.718$. Big!

RK4:

• $k_1 = f(0,1)=1$ — why? slope = $y$ = 1.
• $k_2 = f(0.5,\,1+0.5\cdot1)=1.5$ — why? nudge $y$ by half-step using $k_1$.
• $k_3 = f(0.5,\,1+0.5\cdot1.5)=1.75$ — why? refine midpoint using $k_2$.
• $k_4 = f(1,\,1+1\cdot1.75)=2.75$ — why? end slope using $k_3$.
• $y_1 = 1 + \frac16(1 + 3 + 3.5 + 2.75) = 1 + \frac{10.25}{6} \approx 2.7083$.

Error $\approx 0.0099$ — 70× better than Euler with the same step size.

Worked Example 2 — error scaling check ($y'=y$, integrate to $t=1$)

Why this matters: halving $h$ helps Euler twice but helps RK4 sixteen times — the entire reason RK4 dominates practical computing.

Flashcards

Connections

• Taylor Series Expansion — both methods derived by truncating it.
• Simpson's Rule — RK4 weights are analogous to it.
• Finite Difference Approximation of Derivatives — Euler is forward difference.
• Numerical Stability and Stiff ODEs — why explicit methods can blow up.
• Adaptive Step Size (RK45 / Dormand–Prince) — practical extension of RK4.
• scipy.integrate.solve_ivp — library version of what we built.

Concept Map

Hinglish (regional understanding)