Exercises — Implementing ODE solvers from scratch — Euler, RK4

Euler evaluates the slope only at the left endpoint of the step, at the point $(t_n,y_n)$, and then assumes that one slope stays constant across the whole interval $[t_n,t_n+h]$. $y_{n+1}=y_n+hf(t_n,y_n)$ Picture a straight line leaving $(t_n,y_n)$ with that slope; you slide along it by a horizontal distance $h$ and read off the new height $y_{n+1}$.

Slope at the start: $f(t_0,y_0)=f(0,0)=2\cdot 0=0$. $y_1=y_0+hf(0,0)=0+0.5\cdot 0=0$ Here $f$ depends only on $t$, and at $t=0$ the slope is genuinely $0$, so Euler predicts no rise on the first step (the true solution $y=t^2$ is also flat at the origin). $y_1=0$.

Step 1 (from $t=0$): slope $=f(0,1)=1$. $y_1=1+0.5\cdot 1=1.5(t_1=0.5)$ Step 2 (from $t=0.5$): slope $=f(0.5,1.5)=1.5$. $y_2=1.5+0.5\cdot 1.5=2.25(t_2=1.0)$ Error $=2.71828-2.25=0.46828$. The curve $e^t$ bends upward, so Euler — using the lower left-endpoint slope each time — undershoots.

Here $f(t,y)=\sin t$, so each $k_i$ is just $\sin$ of a shifted time:

Euler's global error is $O(h)$ — halving $h$ roughly halves the error (the ratio $0.125/0.065\approx 1.9\approx 2$ confirms this). Halving again from $0.05$ to $0.025$: $\text{error}(0.025)\approx \frac{0.065}{2}\approx 0.0325$ Why $O(h)$? Each step throws away a Taylor term of size $O(h^2)$ (the local truncation error). Over $N=T/h$ steps these accumulate: $N\cdot O(h^2)=\frac{T}{h}\cdot O(h^2)=O(h)$. See Taylor Series Expansion.

RK4 has global error $O(h^4)$, so halving $h$ shrinks the error by $2^4=16$: $\text{error}(0.05)\approx \frac{2.1\times 10^{-6}}{16}\approx 1.31\times 10^{-7}$ $\text{error}(0.025)\approx \frac{1.31\times 10^{-7}}{16}\approx 8.2\times 10^{-9}$ This is why RK4 dominates: the same effort per step buys a sixteenfold accuracy gain each halving, versus Euler's mere doubling.

Euler follows the tangent line drawn at the left endpoint. For a curve that is concave up (bending upward, like $e^t$), the tangent lies below the curve, so Euler undershoots (as in L2.1: $2.25<2.718$). For a curve that is concave down, the tangent lies above the curve, so Euler overshoots. The sign of the error is set by the sign of $y^{\prime \prime }$ — precisely the $\frac{h^2}{2}{y}^{\prime \prime }$ Taylor term Euler discards. RK4 largely cancels this because it peeks at the midpoint and end before committing.

Number of steps $N=1/h$, and error $\sim C{h}^p$, so to reach $\epsilon$ we need $h\sim (\epsilon /C{)}^{1/p}$, i.e. $N\sim {\epsilon }^{-1/p}$ (ignoring the constant $C$).

• Euler $p=1$: $N\sim {\epsilon }^{-1}=10^8$ steps — a hundred million.
• RK4 $p=4$: $N\sim {\epsilon }^{-1/4}=10^2=100$ steps. Even though RK4 costs 4 function calls per step, its total work $\approx 4\times 100=400$ evaluations crushes Euler's $10^8$. Verdict: at any demanding accuracy, RK4 wins by orders of magnitude — this is the whole reason it is the default workhorse.

With $f(t,y)=\lambda y$ (write $z=\lambda h$):

Each step multiplies by $R(z)$, so after $n$ steps $y_n=R(z{)}^n{y}_0$. For the numerical solution to decay like the true solution ($e^{\lambda t}\to 0$ when $\lambda <0$), we need $|R(z)|\le 1.$ For real $z=\lambda h<0$ this holds only while $z\gtrsim -2.79$ (the real stability boundary of RK4, the most negative $z$ where $|R(z)|=1$). So for a strongly negative $\lambda$ you must keep $h\lesssim 2.79/|\lambda |$ — a very small step — or the numbers blow up despite the true answer shrinking. This is the door to Numerical Stability and Stiff ODEs: stiff problems have huge $|\lambda |$, forcing tiny $h$ and motivating implicit solvers.

Treat $\mathbf{y}$ and each $k_i$ as vectors; $f$ returns a vector. In Python with NumPy arrays the code is identical — the arithmetic y + h/2 * k1, the weighted sum, all broadcast componentwise:

def rk4_system(f, t0, y0, h, n):   # y0 is a numpy array
    t, y = t0, y0.copy()
    ts, ys = [t], [y.copy()]
    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.copy())
    return ts, ys

The update is linear in the $k_i$, so it applies to each component independently — the same weighted average of slopes, now one slope per dimension. This is how a 2nd-order ODE $y^{\prime \prime }=g$ is solved: rewrite as the system $\mathbf{y}=(y,v)$, ${\mathbf{y}}^{\prime }=(v,g)$.

Let $\mathbf{y}=(y,v)$ with $v=y^{\prime }$. Then ${\mathbf{y}}^{\prime }=(v,-y)$, start ${\mathbf{y}}_0=(1,0)$, $h=\frac{\pi }{2}\approx 1.5708$. Each stage is a vector; we track both components and read off the position at the end.

• (a) $h=0$: every update is $y_{n+1}=y_n+0=y_n$ and $t$ never advances. The loop runs but produces a constant list at $t_0$ — mathematically vacuous (no time passes). Guard against it in code.
• (b) $n=0$: the for loop body never executes; you return just the initial point $([t_0],[y_0])$. Correct — asking for zero steps.
• (c) $f\equiv 0$: all $k_i=0$, so $y_{n+1}=y_n$ forever — exactly the true constant solution. Both Euler and RK4 are exact here (zero error), because a constant is a degree-0 polynomial well within RK4's order-4 exactness.
• (d) discontinuous $f$: RK4's order-4 accuracy assumes $f$ is smooth (it relies on the Taylor expansion of $y$). A jump inside $[t_n,t_n+h]$ breaks that assumption; you silently lose accuracy. Fix: place a grid point exactly at the discontinuity so no step straddles it — the same event-handling that scipy.integrate.solve_ivp provides.