diff options
Diffstat (limited to 'content/know/concept/runge-kutta-method/index.pdc')
| -rw-r--r-- | content/know/concept/runge-kutta-method/index.pdc | 267 |
1 files changed, 267 insertions, 0 deletions
diff --git a/content/know/concept/runge-kutta-method/index.pdc b/content/know/concept/runge-kutta-method/index.pdc new file mode 100644 index 0000000..ac2eabf --- /dev/null +++ b/content/know/concept/runge-kutta-method/index.pdc @@ -0,0 +1,267 @@ +--- +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. |