Explanation and proof of the 4th order Runge-Kutta method

The 4th order Runge-Kutta (RK4) method is a numerical technique used to solve ordinary differential equations (ODEs) of the following form

$$\frac{dy}{dx} = f(x,y), \qquad y(0)=y_0$$

It gives $$y_{i+1}$$ in the form

$$y_{i+1} = y_i+(a_1k_1+a_2k_2+a_3k_3+a_4k_4)h$$

where the constants are found to be:

$$y_{i+1} = y_i+\frac{1}{6}(k_1+2k_2+2k_3+k_4)h$$

and:

\begin{aligned} k_1 &= f(x_i, y_i)\\ k_2 &= f \left(x_i+\frac{1}{2}, y_i+\frac{1}{2}k_1h \right)\\ k_3 &= f \left(x_i+\frac{1}{2}, y_i+\frac{1}{2}k_2h \right)\\ k_4 &= f(x_i+h, y_i+k_3h) \end{aligned}

I also learnt that these are derived from the first four terms Taylor series:

$$y_{i+1}=y_i+f(x_i,y_i)h+\frac{1}{2!}f'(x_i,y_i)h^2+\frac{1}{3!}f''(x_i,y_i)h^3+\frac{1}{4!}f'''(x_i,y_i)h^4$$

These are the things that I understood. However, here are what I cannot:

• I cannot see why this method "works"?
• How are the $$a_ik_i$$ terms derived and computed from the Taylor series?
• What justifies or proves this result?

Please note that I have seen many articles on these things things, but none of them contained a full proof whatsoever... I would also much like to get some examples where the f function is actually a function of both x, y not only x. Could you please send me some links with further advanced and concise links of this question?

There are many Runge–Kutta methods. The one you have described is (probably) the most popular and widely used one. I am not going to show you how to derive this particular method – instead I will derive the general formula for the explicit second-order Runge–Kutta methods and you can generalise the ideas.

In what follows we will need the following Taylor series expansions.

\begin{align} \tag1\label{eq1} & x(t+\Delta t) = x(t) + (\Delta t)x'(t) + \frac{(\Delta t)^{2}}{2!}x''(t) + \text{higher order terms}. \\ \tag2\label{eq2} & f(t + \Delta t, x + \Delta x) = f(t,x) + (\Delta t)f_{t}(t,x) + (\Delta x)f_{x}(t,x) + \text{higher order terms}. \end{align}

We are interested in the following ODE:

$$x'(t) = f(t,x(t)).$$

The value $x(t)$ is known and $x(t+h)$ is desired. We can solve this ODE by integrating: $$x(t+h) = x(t) + \int_{t}^{t+h}f(\tau,x(\tau))\text{d}\tau.$$ Unfortunately, actually doing that integration exactly is often quite hard (or impossible), so we approximate it using quadrature: $$x(t+h) \approx x(t) + h\sum_{i=1}^{N}\omega_{i}f(t + \nu_{i}h,x(t + \nu_{i}h)).$$ The accuracy of the quadrature depends on the number of terms in the summation (the order of the Runge–Kutta method), the weights, $\omega_{i}$, and the position of the nodes, $\nu_{i}$.

Even this quadrature can be quite difficult to calculate, since on the right-hand side we need $x(t + \nu_{i}h)$, which we don’t know. We get around this problem in the following manner:

Let $\nu_1=0$, which makes the first term in the quadrature $K_{1} := hf(t,x(t)).$ This we do know and we can also use it to approximate $x(t + \nu_{2}h)$ by writing the second term in the quadrature as $K_{2} := hf(t+\alpha h,x(t) + \beta K_{1})$.

With this, the quadrature formula is: $$\tag3\label{eq3} x(t+h) = x(t) + \omega_{1}K_{1} + \omega_{2}K_{2}.$$

Some notes:

1. If we wanted to find a third-order method, we would introduce $K_{3} = hf(t+\tilde{\alpha}h,x+\tilde{\beta}K_{1} + \tilde{\gamma}K_{2})$. If we wanted a fourth-order method, we would introduce $K_{4}$ in a similar way.
2. This is an explicit method, i.e. we have chosen $K_{2}$ to depend on $K_{1}$ but $K_{1}$ does not depend on $K_{2}$. Similarly, $K_{3}$ (if we introduce it) depends on $K_{1}$ and $K_{2}$ but they do not depend on $K_{3}$. We could allow the dependence to run both ways but then the method is implicit and becomes much harder to solve.
3. We still need to choose $\alpha$, $\beta$ and the $\omega_{i}$. We will do that now using the Taylor series expansions I introduced at the very beginning.

In equation \eqref{eq3}, we substitute in the Taylor series expansion \eqref{eq1} on the left-hand side:

$$x(t) + h x'(t) + \frac{h^{2}}{2!}x''(t) + \mathcal{O}\left(h^{3}\right) = x(t) + \omega_{1}K_{1} + \omega_{2}K_{2}.$$

Since $x' = f$ and thus $x'' = f_{t} + ff_{x}$ (suppressing arguments for ease of notation), we have:

$$hf + \frac{h^{2}}{2}(f_{t} + ff_{x}) + \mathcal{O}\left(h^{3}\right) = \omega_{1}K_{1} + \omega_{2}K_{2}.$$

Now substitute for $K_{1}$ and $K_{2}$:

$$hf + \frac{h^{2}}{2}(f_{t} + ff_{x}) + \mathcal{O}\left(h^{3}\right) = \omega_{1}hf + \omega_{2}hf(t + \alpha h, x + \beta K_{1}).$$

Taylor-expand on the right-hand side using \eqref{eq2}:

$$hf + \frac{h^{2}}{2}(f_{t} + ff_{x}) + \mathcal{O}\left(h^{3}\right) = \omega_{1}hf + \omega_{2}(hf +\alpha h^{2}f_{t} + \beta h^{2} ff_{x}) + \mathcal{O}\left(h^{3}\right).$$

Thus the Runge–Kutta method will agree with the Taylor series approximation to $\mathcal{O}\left(h^{3}\right)$ if we choose:

$$\omega_{1} + \omega_{2} = 1,$$ $$\alpha \omega_{2} = \frac{1}{2},$$ $$\beta \omega_{2} = \frac{1}{2}.$$

The canonical choice for the second-order Runge–Kutta methods is $\alpha = \beta = 1$ and $\omega_{1} = \omega_{2} = 1/2.$

The same procedure can be used to find constraints on the parameters of the fourth-order Runge–Kutta methods. The canonical choice in that case is the method you described in your question.

• When writing 'We can solve this ODE by integrating.', how does r come, eg. why do you integrate with respect to r instead of x? – gen Oct 30 '13 at 6:57
• The ODE has the derivative $dx/dt$ so we integrate on both sides with respect to time. On the right-hand side I'm using the dummy variable $\tau$ instead of $t$ since I want to use $t$ for the limits of integration. – in_mathematica_we_trust Oct 30 '13 at 7:13
• Why do you writing $\alpha, \beta$ here $K_{2} := hf(t+\alpha h,x(t) + \beta K_{1}).$ If we expand $x(t+v_2h)$ with Taylor series we get that $\alpha = \beta = v_2.$ I mean, why do you use different letters for the same quantity. – Yola Aug 19 '15 at 16:50
• @in_wolframAlpha_we_trust : I am wondering how did you arrive at $$y(t)+ \beta k_{1}=y(t+\nu _{2}h)$$ , which is part of the second term of the quadrature , i.e. $$h\cdot f(t+\alpha h,y(t)+\beta k_{1})=h\cdot f(t+\nu _{2}h,y(t+\nu _{2}h))$$. Are you assuming function f is linear ? What if it is nonlinear? – Heber Nov 16 '15 at 16:34
• @in_wolframAlpha_we_trust : I also believe this needs a little explanation. The thing is that $y(t+\nu_2 h) = y(t) + \nu_2 h y'(\xi)$ with $\xi$ between $t$ and $t + \nu_2 h$ (the remainder of Taylor in first order) now you can simply define $\beta=\nu_2 y'(t)/y'(\xi)$....maybe not since then $\beta$ would be a function of $t$, right? – Herman Jaramillo Oct 17 '17 at 11:14