path: root/content/know/concept/runge-kutta-method/index.pdc

---
title: "Runge-Kutta method"
firstLetter: "R"
publishDate: 2022-03-10
categories:
- Mathematics
- Numerical methods

date: 2022-03-07T14:10:18+01:00
draft: false
markup: pandoc
---

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