# 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?

• I would suggest you check the links below, you can find there documents, videos, matlab (mathematica, maple) codes and examples about rk method 2nd order 4th order
– Ömer
Oct 16, 2013 at 19:50
• @Kyle I am sorry for not being clear, so here is a link: damtp.cam.ac.uk/lab/people/sd/lectures/nummeth98/…, I am interested why these formulas approximate...
– gen
Oct 20, 2013 at 8:41
• It seems to me that group theory is used to derive them: youtube.com/watch?v=fLxSCQ4ZEKc Nov 30, 2020 at 22:40

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, 2013 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. Oct 30, 2013 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, 2015 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? Nov 16, 2015 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? Oct 17, 2017 at 11:14