# 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

### A specific answer to the posed question

The classical RK4, invented by Kutta following the patterns of the methods by Heun, is based on the Simpson method. For quadrature this would be $$y(x_{n+1})-y(x_n)=\int_{x_n}^{x_{n+1}}y'(x)\,dx\approx\frac{h}{6}\Bigl(y'(x_n)+4y'(x_{n+1/2})+y'(x_{n+1})\Bigr)$$

For a differential equation, here autonomous for simplicity, we have $$y'(x)=f(y(x))$$. Then, seen from $$x_n$$, the points $$y(x_{n+1/2})$$, $$y(x_{n+1})$$ are unknown future values, and their slopes or $$f$$ values need to be approximated.

#### A simplistic approximation

Consider the first step where $$y_0=y(x_0)$$, $$y_n$$ the numerical approximations of the values $$y(x_n)$$ of the exact solution. As the most simple approximation use a linear continuation like in the Euler method, $$y(x_0+s)\approx y_0+f_0s$$, $$f_0=k_1=f(y_0)$$. Then $$f_{0+1/2}=k_2=f(y_0+\frac h2 f_0)$$ and $$f_1=k_3=f(y_0+h f_0)$$ have errors $$O(h^2)$$ against the values for the exact solution. Inserting these into the Simpson formula results in a local error of $$O(h^3)$$, which sums up in the usual manner to a global error $$O(h^2)$$ at any fixed end point $$x_{final}$$ of the integration.

### Raising the error order

To increase the order of the method one has to improve the approximation of the midpoint and the end point of the current interval. The local error order of $$f_{n+1}$$ is $$O(h^3)$$ or better if $$f_{n+1/2}$$ has order $$O(h^2)$$, this follows as an instance of the midpoint method.

To improve the error order of the midpoint one can apply the 2nd order Heun method $$k_1=f_0$$, $$k_2=f(y_0+\frac h2 k_1)$$, $$x_{0+1/2}=\frac h4(k_1+k_2)$$. One could derive a method of order 3 $$\begin{array}{c|ccc} 0&\\ \frac12&\frac12\\ \frac12&\frac14&\frac14\\ 1&0&0&1\\ \hline&\frac16&\frac23&\frac16 \end{array}$$

### A final trick

The increased order at the midpoint $$y_{0+1/2}$$ over the Euler methods is due to the symmetry of the leading error terms in $$k_1$$ and $$k_2$$. This symmetry also occurs in the values of $$f$$ at these points. Thus set $$k_2=f(y_0+\frac h2 k_1)$$ and $$k_3=f(y_0+\frac h2 k_2)$$ so that $$f_{0+1/2}=\frac12(k_2+k_3)$$ is now (at least) an $$O(h^3)$$ approximation for the slope of the exact solution in the midpoint.

Adding $$f_1=k_4=f(x_0+h k_3)$$ as approximation for the slope in the end-point already gives the RK4 scheme, as $$y_1=y_0+\frac h6(f_0+4f_{0+1/2}+f_1)=y_0+\frac h6(k_1+2(k_2+k_3)+k_4)$$. In total this shows that the local error of the RK4 method is at least $$O(h^4)$$, giving a global error order $$3$$ for the method.

The selection of the steps was rather goal-oriented, there are some variants that would achieve the same stated orders, leading to methods of at least order 3.

### For order 4 see ...

For the final step to show that the global error order is indeed 4, a more systematic analysis of the Taylor expansions of the stages of the method and the exact solution is necessary, the simple leverage of the knowledge about lower order methods is no longer sufficient. The expression of the derivatives of $$y$$ in terms of $$f$$ and its derivatives contains increasingly more terms the higher the order. The coefficients of these terms behave independent of each other in the expansions of the stages of the method. The number of order conditions that ensure the cancellation of these terms increases accordingly.

Elements of that analysis can be found at

• Thank you for posting this.
– gen
Mar 26, 2023 at 16:25