Integration methods

Runge-Kutta-Munthe-Kaas (RKMK) methods with variable step size

The underlying idea of RKMK methods is to express a vector field \(F\in\mathfrak{X}(\mathcal{M})\) as \(F\vert_m = \Psi_*(f(m))\vert_m\) , where \(m \in \mathcal{M}\) , \(\Psi_*\) is the infinitesimal generator of \(\Psi\), a transitive action on \(\mathcal{M}\), 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}\), on which we can perform a time step integration. 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

(1)\[\begin{split}\begin{align} \begin{cases} \sigma(0) = 0\in\mathfrak{g},\\ \dot{\sigma}(t) = \textrm{dexp}_{\sigma(t)}^{-1}\circ f\circ \Psi (\exp(\sigma(t)),y_n)\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 \(\sigma_1\approx \sigma(h_n)\in\mathfrak{g}\) is computed with a Runge-Kutta method.

One approach for varying the step size is based on embedded Runge-Kutta pairs for vector spaces. This approach consists of a principal method of order \(p\), used to propagate the numerical solution, together with some auxiliary method, of order \(\tilde{p}<p\), that is only used to obtain an estimate of the local error. This local error estimate is in turn used to derive a step size adjustment formula that attempts to keep the local error estimate approximately equal to some user-defined tolerance \(\textrm{tol}\) in every step. Both methods are applied to solve the ODE for \(\sigma(t)\) in (1), yielding two approximations \(\sigma_1\) and \(\tilde{\sigma}_1\) respectively, using the same step size \(h_n\). Now, some distance measure between \(\sigma_1\) and \(\tilde{\sigma}_1\) provides an estimate \(e_{n+1}\) for the size of the local truncation error. Thus, \(e_{n+1}=C h_{n+1}^{\tilde{p}+1}+\mathcal{O}(h^{\tilde{p}+2})\). Aiming at \(e_{n+1}\approx\textrm{tol}\) in every step, one may use a formula of the type

(2)\[\begin{align} h_{n+1} = \theta\left(\frac{\textrm{tol}}{e_{n+1}}\right)^{\tfrac{1}{\tilde{p}+1}}\, h_n \end{align}\]

where \(\theta\) is typically chosen between \(0.8\) and \(0.9\). If \(e_n>\textrm{tol}\), the step is rejected. Hence, we can redo the step with the step size obtained by the same formula.

In our code, we perform the numerical experiments on the N-fold 3D pendulum, in which we compare the performance of constant and variable step size methods. We consider the RKMK pair coming from Dormand–Prince method (DOPRI 5(4), which we denote by RKMK(5,4)). We set a tolerance of \(10^{-6}\) and solve the system with the RKMK(5,4) scheme. Fixing the number of time steps required by RKMK(5,4), we repeat the experiment with RKMK of order 5 (denoted by RKMK5). The comparison occurs at the final time \(T=3\) using the Euclidean norm of the ambient space \(\mathbb{R}^{6N}\), where \(N\) is the number of the connected pendula. The quality of the approximation is measured against a reference solution obtained with ODE45 from MATLAB with a strict tolerance.