Integration methods

Runge-Kutta-Munthe-Kaas methods

Lie group integrators solve differential equations whose solution evolve on a manifold \(\mathcal{M}\), i.e. the solution is a curve \(y(t)\in\mathcal{M}\) whose tangent at any point coincides with a vector field \(F\in\mathcal{X}(\mathcal{M})\) and passes through a designated initial value \(y_0\) at \(t=t_0\):

(1)\[\begin{align} \dot{y}(t) = F\left(y(t)\right),\qquad y(t_0)=y_0. \end{align}\]

Notice that we restrict the discussion to the case of autonomous vector field (explicit time dependence could easily be included). Let \(G\) be a Lie group acting transitively on \(\mathcal{M}\) via the group action \(\psi:G \times \mathcal{M} \rightarrow \mathcal{M}\), so that \(\mathcal{M}\) is a homogeneous manifold. The underlying idea of Runge-Kutta-Munthe-Kaas (RKMK) methods is to express a vector field \(F\in\mathfrak{X}(\mathcal{M})\) as \(F\vert_m = \psi_*(f(m))\vert_m\) , where \(\psi_*\) is the infinitesimal generator of \(\psi\) and \(f:\mathcal{M}\rightarrow\mathfrak{g}\). This allows us to transform the problem from the manifold \(\mathcal{M}\) to the Lie algebra \(\mathfrak{g}\) of \(G\), on which we can perform a time step integration with a Runge-Kutta method. We then map the result back to \(\mathcal{M}\), and repeat this up to the final integration time. More explicitly, let \(h_n\) be the size of the \(n-th\) time step, we then update \(y_n\in\mathcal{M}\) to \(y_{n+1}\) by

(2)\[\begin{split}\begin{align} \begin{cases} \sigma(0) = 0\in\mathfrak{g},\\ \dot{\sigma}(t) = \textrm{dexp}_{\sigma(t)}^{-1}\circ f\circ \psi \left(\exp(\sigma(t)),y_n\right)\in T_{\sigma(t)}\mathfrak{g}, \\ y_{n+1} = \psi(\exp(\sigma_1),y_n)\in \mathcal{M}, \end{cases} \end{align}\end{split}\]

where \(\textrm{exp}:\mathfrak{g}\rightarrow G\) is the exponential map, and \(\sigma_1\approx \sigma(h_n)\in\mathfrak{g}\) is computed with a Runge-Kutta method.

The transformed differential equation for \(\sigma(t)\) makes use of the derivative of the exponential mapping. The map \(v\mapsto\textrm{dexp}_u(v)\) is linear and invertible when \(u\) belongs to some sufficiently small neighborhood of \(0\in\mathfrak{g}\). It has an expansion in nested Lie brackets and, using the operator \(\textrm{ad}_u(v)=[u,v]\) and its powers \(\textrm{ad}_u^2 v=[u,[u,v]]\) etc, one can write

(3)\[\begin{align} \textrm{dexp}_u(v) = \left.\frac{e^z-1}{z}\right|_{z=\textrm{ad}_u}(v) = v + \frac12[u,v] + \frac16[u,[u,v]] + \cdots. \end{align}\]

The inverse is

(4)\[\begin{align} \textrm{dexp}_u^{-1}(v) =\left.\frac{z}{e^z-1}\right|_{z=\textrm{ad}_u}(v)= v -\frac12[u,v] + \frac1{12}[u,[u,v]]+\cdots. \end{align}\]

To evaluate \(\textrm{dexp}_u^{-1}(v)\) one can either truncate the series (3), or compute its exact expression for the particular Lie algebra under consideration. The exponential maps on the Lie groups SO(3) and SE(3) are implemented in expSO3 and expSE3. The exact expressions for the inverse of the derivative of the exponential map on SO(3) and SE(3) are implemented in dexpinvSO3 and dexpinvSE3.

Examples

Let us consider an s-stage Runge-Kutta (RK) method for the initial value problem (1):

(5)\[\begin{align} y_{n+1}=y_n+h \sum_{i=1}^s b_i k_i, \quad k_i=F\left(y_n+h \sum_{j=1}^s a_{i j} k_j \right), \quad i=1, \ldots, s, \end{align}\]

where \(b_i,\,a_{ij}\, (i,\,j=1,\dots\,s)\) are real numbers called, respectively, the weights and coefficients of the method, and \(c_i=\sum_{j=1}^s a_{ij}\) are called the nodes or abscissae. These constants define a specific RK method and can be collected in the following table, known as Butcher’s tableau:

(6)\[\begin{split}\begin{align} \begin{array}{c|cccc} c_1 & a_{11} & a_{12} & \ldots & a_{1 s} \\ c_2 & a_{21} & a_{22} & \ldots & a_{2 s} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_s & a_{s 1} & a_{s 2} & \ldots & a_{s s} \\ \hline & b_1 & b_2 & \ldots & b_s \end{array} \end{align}\end{split}\]

From equation (2) it follows that one step of the resulting Runge–Kutta–Munthe-Kaas method writes

(7)\[\begin{split}\begin{align} &y_1=\exp \left(h \sum_{i=1}^s b_i k_i\right) \cdot y_0,\\ &k_i=\operatorname{dexp}^{-1}_{h \sum_j a_{i j} k_j} f\left(\exp \left(h \sum_j a_{i j} k_j\right) \cdot y_0\right), \quad i=1, \ldots, s, \end{align}\end{split}\]

where we denote the group action by “\(\cdot\)” for ease of notation.

The simplest Lie group integrator is the Lie-Euler method, based on the classical explicit Euler method, a first-order method with Butcher’s tableau given by

(8)\[\begin{split}\begin{align} \begin{array}{c|c} 0 & 0 \\ \hline & 1 \end{array} \end{align}\end{split}\]

The resulting Lie-Euler method is implemenmted in LieEuler.

An improvement to the Lie-Euler method is the second-order RKMK method based on the tableau of the Heun’s method,

(9)\[\begin{split}\begin{align} \begin{array}{c|cc} 0 & 0 & 0 \\ 1 & 1 & 0 \\ \hline & 1/2 & 1/2 \end{array}, \end{align}\end{split}\]

and implemented in RKMK2Heun

The following Butcher’s tables provide the coefficients for two classical methods of order three (on the left) and order four (on the right):

(10)\[\begin{split}\begin{align} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ {1/2} & {1/2} & 0 & 0 \\ 1 & -1 & 2 & 0 \\ \hline & {1/6} & {2/3} & {1/6} \end{array} \qquad \qquad \quad \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0 \\ {1/2} & {1/2} & 0 & 0 & 0 \\ {1/2} & 0 & {1/2} & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ \hline & {1/6} & {1/3} & {1/3} & {1/6} \end{array} \end{align}\end{split}\]

The corresponding RKMK integrators are implemented in RKMK3 and RKMK4.

Commutator-free methods

The second class of Lie group integrators to be considered here are the commutator-free methods, named this way to emphasize the contrast to RKMK schemes which usually include commutators in the method format. These schemes include the Crouch-Grossman methods and have the format

(11)\[\begin{split}\begin{align} Y_{n,r} &= \exp\Big(h\sum_{k}\alpha_{r,J}^k f_{n,k}\Big)\cdots \exp\Big(h\sum_{k}\alpha_{r,1}^k f_{n,k}\Big) \cdot y_n\\ f_{n,r} &= f(Y_{n,r}) \\[1mm] y_{n+1} &= \exp\Big(h\sum_k \beta_J^k f_{n,k}\Big)\cdots \exp\Big(h\sum_k \beta_1^k f_{n,k}\Big) \cdot y_n \end{align}\end{split}\]

where “\(\cdot\)” denotes the group action. Here the Runge-Kutta coefficients \(\alpha_{r,j}^k\), \(\beta_{j}^r\) are related to a classical Runge-Kutta scheme with coefficients \(a_r^k\), \(b_r\) in that \(a_r^k=\sum_j \alpha_{r,j}^k\) and \(b_r=\sum_j \beta_{j}^r\). The \(\alpha_{r,j}^k\), \(\beta_{j}^r\) are usually chosen to obtain computationally inexpensive schemes with the highest possible order of convergence. The computational complexity of the above schemes depends on the cost of computing an exponential as well as of evaluating the vector field. Therefore it makes sense to keep the number of exponentials \(J\) in each stage as low as possible, and possibly also the number of stages \(s\).

The following example is a generalization of the classical fourth-order Runge–Kutta method in (10) and is implemented in CFree4:

(12)\[\begin{split}\begin{aligned} &Y_1=y_0, \\ &Y_2=\exp \left(\frac{1}{2} k_1\right) \cdot y_0, \\ &Y_3=\exp \left(\frac{1}{2} k_2\right) \cdot y_0 \\ &Y_4=\exp \left(k_3-\frac{1}{2} k_1\right) \cdot Y_2, \\ &y_{\frac{1}{2}}=\exp \left(\frac{1}{12}\left(3 k_1+2 k_2+2 k_3-k_4\right)\right) \cdot y_0, \\ &y_1=\exp \left(\frac{1}{12}\left(-k_1+2 k_2+2 k_3+3 k_4\right)\right) \cdot y_{\frac{1}{2}}, \end{aligned}\end{split}\]

with \(k_i=hf(Y_i)\). We notice that one exponential is saved in computing \(Y_4\) by making use of \(Y_2\). This shows that sometimes it is possible to come up with tricks that allow to reuse exponentials from one stage to another, thereby lowering the computational cost of the scheme.

We refer to (Celledoni, Çokaj, Leone, Murari and Owren, 2021) and references cited therein for further details.