Categories: Mathematics, Numerical methods.

# Runge-Kutta method

A Runge-Kutta method (RKM) is a popular approach to numerically solving systems of ordinary differential equations. Let $$\vb{x}(t)$$ be the vector we want to find, governed by $$\vb{f}(t, \vb{x})$$:

\begin{aligned} \vb{x}'(t) = \vb{f}\big(t, \vb{x}(t)\big) \end{aligned}

Like in all numerical methods, the $$t$$-axis is split into discrete steps. If a step has size $$h$$, then as long as $$h$$ is small enough, we can make the following approximation:

\begin{aligned} \vb{x}'(t) + a h \vb{x}''(t) &\approx \vb{x}'(t \!+\! a h) \\ &\approx \vb{f}\big(t \!+\! a h,\, \vb{x}(t \!+\! a h)\big) \\ &\approx \vb{f}\big(t \!+\! a h,\, \vb{x}(t) \!+\! a h \vb{x}'(t) \big) \end{aligned}

For sufficiently small $$h$$, higher-order derivates can also be included, albeit still at $$t \!+\! a h$$:

\begin{aligned} \vb{x}'(t) + a h \vb{x}''(t) + b h^2 \vb{x}'''(t) &\approx \vb{f}\big(t \!+\! a h,\, \vb{x}(t) \!+\! a h \vb{x}'(t) \!+\! b h^2 \vb{x}''(t) \big) \end{aligned}

Although these approximations might seem innocent, they actually make it quite complicated to determine the error order of a given RKM.

Now, consider a Taylor expansion around the current $$t$$, truncated at a chosen order $$n$$:

\begin{aligned} \vb{x}(t \!+\! h) &= \vb{x}(t) + h \vb{x}'(t) + \frac{h^2}{2} \vb{x}''(t) + \frac{h^3}{6} \vb{x}'''(t) + \:...\, + \frac{h^n}{n!} \vb{x}^{(n)}(t) \\ &= \vb{x}(t) + h \bigg[ \vb{x}'(t) + \frac{h}{2} \vb{x}''(t) + \frac{h^2}{6} \vb{x}'''(t) + \:...\, + \frac{h^{n-1}}{n!} \vb{x}^{(n)}(t) \bigg] \end{aligned}

We are free to split the terms as follows, choosing real factors $$\omega_{mj}$$ subject to $$\sum_{j} \omega_{mj} = 1$$:

\begin{aligned} \vb{x}(t \!+\! h) &= \vb{x} + h \bigg[ \sum_{j = 1}^{N_1} \omega_{1j} \, \vb{x}' + \frac{h}{2} \sum_{j = 1}^{N_2} \omega_{2j} \, \vb{x}'' + \:...\, + \frac{h^{n-1}}{n!} \sum_{j = 1}^{N_n} \omega_{nj} \, \vb{x}^{(n)} \bigg] \end{aligned}

Where the integers $$N_1,...,N_n$$ are also free to choose, but for reasons that will become clear later, the most general choice for an RKM is $$N_1 = n$$, $$N_n = 1$$, and:

\begin{aligned} N_{n-1} = N_n \!+\! 2 ,\quad \cdots ,\quad N_{n-m} = N_{n-m+1} \!+\! m \!+\! 1 ,\quad \cdots ,\quad N_{2} = N_3 \!+\! n \!-\! 1 \end{aligned}

In other words, $$N_{n-m}$$ is the $$m$$th triangular number. This is not so important, since this is not a practical way to describe RKMs, but it is helpful to understand how they work.

## Example derivation

For example, let us truncate at $$n = 3$$, such that $$N_1 = 3$$, $$N_2 = 3$$ and $$N_3 = 1$$. The following derivation is very general, except it requires all $$\alpha_j \neq 0$$. Renaming $$\omega_{mj}$$, we start from:

\begin{aligned} \vb{x}(t \!+\! h) &= \vb{x} + h \bigg[ (\alpha_1 + \alpha_2 + \alpha_3) \, \vb{x}' + \frac{h}{2} (\beta_2 + \beta_{31} + \beta_{32}) \, \vb{x}'' + \frac{h^2}{6} \gamma_3 \, \vb{x}''' \bigg] \\ &= \vb{x} + h \bigg[ \alpha_1 \vb{x}' + \Big( \alpha_2 \vb{x}' + \frac{h}{2} \beta_2 \vb{x}'' \Big) + \Big( \alpha_3 \vb{x}' + \frac{h}{2} (\beta_{31} + \beta_{32}) \vb{x}'' + \frac{h^2}{6} \gamma_3 \vb{x}''' \Big) \bigg] \end{aligned}

As discussed earlier, the parenthesized expressions can be approximately rewritten with $$\vb{f}$$:

\begin{aligned} \vb{x}(t \!+\! h) = \vb{x} + h &\bigg[ \alpha_1 \vb{f}(t, \vb{x}) + \alpha_2 \vb{f}\Big( t \!+\! \frac{h \beta_2}{2 \alpha_2}, \; \vb{x} \!+\! \frac{h \beta_2}{2 \alpha_2} \vb{x}' \Big) \\ & + \alpha_3 \vb{f}\Big( t \!+\! \frac{h (\beta_{31} \!\!+\!\! \beta_{32})}{2 \alpha_3}, \; \vb{x} \!+\! \frac{h \beta_{31}}{2 \alpha_3} \vb{x}' \!+\! \frac{h \beta_{32}}{2 \alpha_3} \vb{x}' \!+\! \frac{h^2 \gamma_3}{6 \alpha_3} \vb{x}'' \Big) \bigg] \\ = \vb{x} + h &\bigg[ \alpha_1 \vb{k}_1 + \alpha_2 \vb{f}\Big( t \!+\! \frac{h \beta_2}{2 \alpha_2}, \; \vb{x} \!+\! \frac{h \beta_2}{2 \alpha_2} \vb{k}_1 \!\Big) \\ & + \alpha_3 \vb{f}\Big( t \!+\! \frac{h (\beta_{31} \!\!+\!\! \beta_{32})}{2 \alpha_3}, \; \vb{x} \!+\! \frac{h \beta_{31}}{2 \alpha_3} \vb{k}_1 \!+\! \frac{h \beta_{32}}{2 \alpha_3} \vb{f}\Big( t \!+\! \frac{h \gamma_3}{3 \beta_{32}}, \; \vb{x} \!+\! \frac{h \gamma_3}{3 \beta_{32}} \vb{k}_1 \!\Big) \!\Big) \bigg] \end{aligned}

Here, we can see an opportunity to save some computational time by reusing an evaluation of $$\vb{f}$$. Technically, this is optional, but it would be madness not to, so we choose:

\begin{aligned} \frac{\beta_2}{2 \alpha_2} = \frac{\gamma_3}{3 \beta_{32}} \end{aligned}

Such that the next step of $$\vb{x}$$’s numerical solution is as follows, recalling that $$\sum_{j} \alpha_j = 1$$:

\begin{aligned} \boxed{ \vb{x}(t \!+\! h) = \vb{x}(t) + h \Big( \alpha_1 \vb{k}_1 + \alpha_2 \vb{k}_2 + \alpha_3 \vb{k}_3 \Big) } \end{aligned}

Where $$\vb{k}_1$$, $$\vb{k}_2$$ and $$\vb{k}_3$$ are different estimates of the average slope $$\vb{x}'$$ between $$t$$ and $$t \!+\! h$$, whose weighted average is used to make the $$t$$-step. They are given by:

\begin{aligned} \boxed{ \begin{aligned} \vb{k}_1 &\equiv \vb{f}(t, \vb{x}) \\ \vb{k}_2 &\equiv \vb{f}\bigg( t + \frac{h \beta_2}{2 \alpha_2}, \; \vb{x} + \frac{h \beta_2}{2 \alpha_2} \vb{k}_1 \bigg) \\ \vb{k}_3 &\equiv \vb{f}\bigg( t + \frac{h (\beta_{31} \!\!+\!\! \beta_{32})}{2 \alpha_3}, \; \vb{x} + \frac{h \beta_{31}}{2 \alpha_3} \vb{k}_1 + \frac{h \beta_{32}}{2 \alpha_3} \vb{k}_2 \bigg) \end{aligned} } \end{aligned}

Despite the contraints on $$\alpha_j$$ and $$\beta_j$$, there is an enormous freedom of choice here, all leading to valid RKMs, although not necessarily good ones.

## General form

A more practical description goes as follows: in an $$s$$-stage RKM, a weighted average is taken of up to $$s$$ slope estimates $$\vb{k}_j$$ with weights $$b_j$$. Let $$\sum_{j} b_j = 1$$, then:

\begin{aligned} \boxed{ \vb{x}(t \!+\! h) = \vb{x}(t) + h \sum_{j = 1}^{s} b_j \vb{k}_j } \end{aligned}

Where the estimates $$\vb{k}_1, ..., \vb{k}_s$$ depend on each other, and are calculated one by one as:

\begin{aligned} \boxed{ \vb{k}_m = \vb{f}\bigg( t + h c_m,\; \vb{x} + h \sum_{j = 1}^{m - 1} a_{mj} \vb{k}_j \bigg) } \end{aligned}

With $$c_1 = 1$$ and $$\sum_{j = 1} a_{mj} = c_m$$. Writing this out for the first few $$m$$, the pattern is clear:

\begin{aligned} \vb{k}_1 &= \vb{f}(t, \vb{x}) \\ \vb{k}_2 &= \vb{f}\big( t + h c_2,\; \vb{x} + h a_{21} \vb{k}_1 \big) \\ \vb{k}_3 &= \vb{f}\big( t + h c_3,\; \vb{x} + h (a_{31} \vb{k}_1 + a_{32} \vb{k}_2) \big) \\ \vb{k}_4 &= \:... \end{aligned}

The coefficients of a given RKM are usually compactly represented in a Butcher tableau:

\begin{aligned} \begin{array}{c|ccc} 0 \\ c_2 & a_{21} \\ c_3 & a_{31} & a_{32} \\ \vdots & \vdots & \vdots & \ddots \\ c_s & a_{s1} & a_{s2} & \cdots & a_{s,s-1} \\ \hline & b_1 & b_2 & \cdots & b_{s-1} & b_s \end{array} \end{aligned}

Each RKM has an order $$p$$, such that the global truncation error is $$\mathcal{O}(h^p)$$, i.e. the accumulated difference between the numerical and the exact solutions is proportional to $$h^p$$.

The surprise is that $$p$$ need not be equal to the Taylor expansion order $$n$$, nor the stage count $$s$$. Typically, $$s = n$$ for computational efficiency, but $$s \ge n$$ is possible in theory.

The order $$p$$ of a given RKM is determined by a complicated set of equations on the coefficients, and the lowest possible $$s$$ for a desired $$p$$ is in fact only partially known. For $$p \le 4$$ the bound is $$s \ge p$$, whereas for $$p \ge 5$$ the only proven bound is $$s \ge p \!+\! 1$$, but for $$p \ge 7$$ no such efficient methods have been found so far.

If you need an RKM with a certain order, look it up. There exist many efficient methods for $$p \le 4$$ where $$s = p$$, and although less popular, higher $$p$$ are also available.

## References

1. J.C. Butcher, Numerical methods for ordinary differential equations, 3rd edition, Wiley.

© Marcus R.A. Newman, a.k.a. "Prefetch". Available under CC BY-SA 4.0.