Exercises — Runge-Kutta 4th order (RK4) — derivation

$k_3=hf\left(x_n+\frac{h}{2},y_n+\frac{k_2}{2}\right)$. So it is sampled at the midpoint $x_n+\frac{h}{2}$, and it reuses $k_2$ (a refined midpoint slope). Look at figure s01: the two amber arrows ($k_2$ and $k_3$) share the same vertical dashed line at the midpoint; only the height of the estimated $y$ differs.

It should be $\frac{1}{6}$, not $\frac{h}{6}$. In our convention each $k_i$ already contains one factor of $h$ (because $k_i=hf$). Putting $\frac{h}{6}$ would multiply by $h$ twice, giving an $h^2$ step — wrong dimensions. $y_{n+1}=y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4).$ The $\frac{h}{6}$ form is only correct in the other convention where $k_i$ is a pure slope (no $h$).

Here $f(x,y)=x+y$, $x_0=0$, $y_0=1$, $h=0.2$.

• $k_1=0.2f(0,1)=0.2(0+1)=0.2$ — slope at left edge.
• $k_2=0.2f(0.1,1+0.1)=0.2f(0.1,1.1)=0.2(0.1+1.1)=0.24$ — midpoint, using $\frac{k_1}{2}=0.1$.
• $k_3=0.2f(0.1,1+0.12)=0.2f(0.1,1.12)=0.2(0.1+1.12)=0.244$ — refined midpoint, using $\frac{k_2}{2}=0.12$.
• $k_4=0.2f(0.2,1+0.244)=0.2f(0.2,1.244)=0.2(0.2+1.244)=0.2888$ — right edge, using full $k_3$. $y_1=1+\frac{1}{6}(0.2+2(0.24)+2(0.244)+0.2888)=1+\frac{1.4568}{6}=\mathrm{1.2428.}$ Exact: $y=2e^x-x-1\Rightarrow y(0.2)=2e^{0.2}-1.2=1.242806$. Error $\approx 6\times 10^{-6}$.

$f(x,y)=-2y$, $y_0=3$, $h=0.5$.

• $k_1=0.5(-2\cdot 3)=-3$.
• $k_2=0.5(-2(3+\frac{-3}{2}))=0.5(-2\cdot 1.5)=-1.5$.
• $k_3=0.5(-2(3+\frac{-1.5}{2}))=0.5(-2\cdot 2.25)=-2.25$.
• $k_4=0.5(-2(3+(-2.25)))=0.5(-2\cdot 0.75)=-0.75$. $y_1=3+\frac{1}{6}(-3+2(-1.5)+2(-2.25)+(-0.75))=3+\frac{-11.25}{6}=\mathrm{1.125.}$ Let's check the bracket: $-3+2(-1.5)+2(-2.25)+(-0.75)=-3-3-4.5-0.75=-11.25$, and $3-\frac{11.25}{6}=3-1.875=1.125$. Exact: $y=3e^{-2x}\Rightarrow y(0.5)=3e^{-1}=1.103638$. Here the step $h=0.5$ is large for a fast-decaying equation, so the error ($\approx 0.021$) is visible — a reminder that "4th order" is a statement about small $h$.

When $f$ depends only on $x$, the $y$-arguments are irrelevant, so: $k_1=hf(x_n),k_2=k_3=hf\left(x_n+\frac{h}{2}\right),k_4=hf(x_n+h).$ Substituting into the update: $y_{n+1}-y_n=\frac{1}{6}(hf(x_n)+2\cdot hf(x_n+\frac{h}{2})+2\cdot hf(x_n+\frac{h}{2})+hf(x_n+h)).$ $=\frac{h}{6}(f(x_n)+4f\left(x_n+\frac{h}{2}\right)+f(x_n+h)).$ That last line is exactly Simpson's rule $\frac{h}{6}(f_0+4f_{mid}+f_1)$ for the integral over $[x_n,x_n+h]$. So RK4 is Simpson's rule, generalized to allow $y$-dependence. This is why the weights are $1,2,2,1$: the two midpoint samples merge into the single "$4$" of Simpson.

To cross a fixed length $L=b-a$ you take $N=\frac{L}{h}$ steps — that is, $N\propto \frac{1}{h}$.

What "each step adds $\sim C{h}^5$" means. On a single step the method's answer differs from the true curve by an amount bounded by (some constant) $\times h^5$. Call that constant $C$: it packages together how curvy the solution is (roughly, the size of its 5th derivative over the interval) and the fixed numerical factor $\frac{1}{120}$ that RK4 leaves behind. $C$ does not shrink as $h$ shrinks — it depends on the problem, not the step — so the whole $h$-dependence of one step's error lives in the $h^5$. That is why we write "local error $=C{h}^5$".

Why the errors add up (error propagation). Picture the calculation as a relay: each step starts from the previous step's (already slightly wrong) answer and adds its own fresh $C{h}^5$ mistake. So two things happen at every step — (i) a new $\sim C{h}^5$ error is committed, and (ii) the errors already present get carried forward. For a well-behaved (stable) equation the carried-forward errors don't blow up; they grow at most by a bounded factor. In that case the total is dominated by simply summing the $N$ fresh contributions, one per step: $\text{global error}\sim \underset{\text{how many steps}}{\underbrace{N}}\times \underset{\text{one step’s error}}{\underbrace{C{h}^5}}=\frac{L}{h}\cdot C{h}^5=CL{h}^4=O(h^4).$ One power of $h$ is "eaten" by the growing number of steps: you make more steps when $h$ is small, so more little errors pile up, cancelling one factor of $h$. See Local vs Global Truncation Error. This is why RK4 is called 4th order (the global statement), even though each step is $5$th-order accurate.

$f(x,y)=y$, $h=0.5$.

Step 1 ($x_0=0,y_0=1$):

• $k_1=0.5(1)=0.5$
• $k_2=0.5(1+0.25)=0.625$
• $k_3=0.5(1+0.3125)=0.65625$
• $k_4=0.5(1+0.65625)=0.828125$ $y_1=1+\frac{1}{6}(0.5+2(0.625)+2(0.65625)+0.828125)=1+\frac{3.890625}{6}=\mathrm{1.6484375.}$

Step 2 ($x_1=0.5,y_1=1.6484375$):

• $k_1=0.5(1.6484375)=0.82421875$
• $k_2=0.5(1.6484375+0.412109\dots )=1.0302734\dots$
• $k_3=0.5(1.6484375+0.515137\dots )=1.0817871\dots$
• $k_4=0.5(1.6484375+1.081787\dots )=1.3651123\dots$ $y_2=1.6484375+\frac{1}{6}(0.82421875+2(1.0302734)+2(1.0817871)+1.3651123).$ The bracket sum $=6.4134277\dots$, so $y_2=1.6484375+1.0689046\dots =2.7173421$. Exact $e=2.7182818$; error $\approx 0.00094$. The single big step gave $2.70833$ (error $\approx 0.010$). Halving the step cut the error by roughly $16=2^4$ — the signature of a 4th-order method.

Global error $\propto h^4$. Going from $h=0.5$ to $h=0.25$ halves $h$, so the error scales by ${\left(\frac{1}{2}\right)}^4=\frac{1}{16}$. $\text{new error}\approx \frac{0.00094}{16}\approx 5.9\times 10^{-5}.$ So we predict $y(1)\approx e-5.9\times 10^{-5}\approx 2.718222$. (Direct computation gives error $\approx 6.5\times 10^{-5}$ — the $2^4$ rule works.)

With $f(x,y)=y$, $y_0=1$:

• $k_1=h\cdot 1=h.$
• $k_2=h\left(1+\frac{k_1}{2}\right)=h\left(1+\frac{h}{2}\right)=h+\frac{h^2}{2}.$
• $k_3=h\left(1+\frac{k_2}{2}\right)=h\left(1+\frac{h}{2}+\frac{h^2}{4}\right)=h+\frac{h^2}{2}+\frac{h^3}{4}.$
• $k_4=h\left(1+k_3\right)=h\left(1+h+\frac{h^2}{2}+\frac{h^3}{4}\right)=h+h^2+\frac{h^3}{2}+\frac{h^4}{4}.$

Now combine with weights $1,2,2,1$: $k_1+2k_2+2k_3+k_4.$ Collect by power of $h$:

• $h^1$: $1+2+2+1=6$ → contributes $6h$.
• $h^2$: $2\cdot \frac{1}{2}+2\cdot \frac{1}{2}+1=1+1+1=3$ → contributes $3h^2$.
• $h^3$: $2\cdot \frac{1}{4}+1\cdot \frac{1}{2}=\frac{1}{2}+\frac{1}{2}=1$ → contributes $h^3$.
• $h^4$: $1\cdot \frac{1}{4}=\frac{1}{4}$ → contributes $\frac{1}{4}{h}^4$.

So $k_1+2k_2+2k_3+k_4=6h+3h^2+h^3+\frac{1}{4}{h}^4$. Divide by $6$: $y_1=1+h+\frac{h^2}{2}+\frac{h^3}{6}+\frac{h^4}{24}.$ These are exactly the first five terms of $e^h=1+h+\frac{h^2}{2}+\frac{h^3}{6}+\frac{h^4}{24}+\frac{h^5}{120}+\cdots$. The first term RK4 misses is $\frac{h^5}{120}$ — confirming local error $O(h^5)$. See Taylor Series Methods.

Vector RK4 is identical to scalar RK4, just applied componentwise. Let ${\mathbf{u}}_0=(1,0)$, $h=0.1$, $\mathbf{f}(x,y)=(y,-x)$.

$k_1=h\mathbf{f}({\mathbf{u}}_0)$: $\mathbf{f}(1,0)=(0,-1)$, so $k_1=(0,-0.1)$.

$k_2=h\mathbf{f}({\mathbf{u}}_0+\frac{k_1}{2})$: argument $=(1+0,0-0.05)=(1,-0.05)$; $\mathbf{f}=(-0.05,-1)$; so $k_2=(-0.005,-0.1)$.

$k_3=h\mathbf{f}({\mathbf{u}}_0+\frac{k_2}{2})$: argument $=(1-0.0025,0-0.05)=(0.9975,-0.05)$; $\mathbf{f}=(-0.05,-0.9975)$; so $k_3=(-0.005,-0.09975)$.

$k_4=h\mathbf{f}({\mathbf{u}}_0+k_3)$: argument $=(1-0.005,0-0.09975)=(0.995,-0.09975)$; $\mathbf{f}=(-0.09975,-0.995)$; so $k_4=(-0.009975,-0.0995)$.

Combine componentwise, ${\mathbf{u}}_1={\mathbf{u}}_0+\frac{1}{6}(k_1+2k_2+2k_3+k_4)$:

• $x$: $1+\frac{1}{6}(0+2(-0.005)+2(-0.005)+(-0.009975))=1+\frac{-0.029975}{6}=\mathrm{0.99500417.}$
• $y$: $0+\frac{1}{6}(-0.1+2(-0.1)+2(-0.09975)+(-0.0995))=\frac{-0.59900}{6}=-\mathrm{0.09983333.}$

So $x(0.1)\approx 0.9950042,y(0.1)\approx -0.0998333$. Exact: $\cos 0.1=0.9950042,-\sin 0.1=-0.0998334$. Both match to about $2\times 10^{-6}$ — RK4 handles systems with the same 4th-order accuracy, no new machinery needed.