Eight-stage pseudo-symplectic Runge-Kutta methods of order (4, 8)

TL;DR

This work constructs explicit pseudo-symplectic Runge–Kutta methods that preserve the symplectic structure up to order $$(4,8)$$ by exploiting time-reversal symmetry and parity-based simplifying assumptions, and it derives a 1-dimensional family of 8-stage schemes achieving this order. It also yields a 7-stage member attaining $$(4,9)$$ and analyzes the associated stability function and algebraic nature of the coefficients, all while enforcing nonnegative weights and increasing nodes for practicality. The authors validate the methods on three Hamiltonian benchmarks, showing superior long-time energy and invariant preservation compared with classical explicit schemes, and demonstrating the numerical viability of the new high-order explicit pseudo-symplectic integrators. The results offer a robust, explicit alternative for accurate long-time simulation of Hamiltonian and near-Hamiltonian systems, with improved conservation properties over existing pseudo-symplectic methods.

Abstract

Using simplifying assumptions that are related to the time reversal symmetry, a 1-dimensional family of 8-stage pseudo-symplectic Runge-Kutta methods of order (4, 8), i.e., methods of order 4 that preserve symplectic structure up to order 8, is derived. An example of 7-stage method of order (4, 9) is given.

Figures (8)

• Figure 1: The curve $\zeta(c_2 ,\mkern2mu c_3) = 0$ (solid curve); the line $c_3 = \frac{1}{6}$ (dashed curve); and the hyperbola $6 c_3 - 12 c_3 (c_3 - c_2) - 1 = 0$ (dotted curve, the coefficient $a_{43}$ is infinite on it) which intersects with the curve $\zeta = 0$ at the points $(\frac{1}{6} ,\mkern2mu \frac{1}{6})$, $(\frac{1}{6} ,\mkern2mu \frac{1}{2})$, $\mathsf{Q} = \bigl( (6 - \sqrt{21}) / 15 ,\mkern2mu (9 - \sqrt{21}) / 12 \bigr) = (0.09449..., 0.36811...)$, and $\mathsf{Q}' = \bigl( (6 + \sqrt{21}) / 15 ,\mkern2mu (9 + \sqrt{21}) / 12 \bigr) = (0.70550..., 1.13188...)$. Some branches of the curve $\zeta = 0$ are not shown as they correspond to too large values of $c_2$ or $c_3$ to be visible on the graph. The points $\mathsf{C}$ and $\mathsf{A}$ are indicated by open dots and correspond to the method in eq. (\ref{['method_C']}) and the $1$-dimensional family of methods that includes the method in eq. (\ref{['the_method']}), respectively.
• Figure 2: Quadrature scheme graphical depiction (see also Suz90, where the vertical coordinate is the stage) for the eight methods listed in Table \ref{['table_methods']}. Here ${\Sigma_{\mkern1mu i} \mkern1.5mu {\pmb{b}}} = \sum_{j = 1}^{\mkern1.5mu i} {b_{\mkern-1.25mu j}}$, where $1 \le i \le s$, are the cumulative weights, with $\Sigma_{\mkern1mu 0} \mkern1.5mu {\pmb{b}} = 0$. Due to the order condition ${\pmb{b}} \mkern0.5mu {\pmb{1}} = 1$ for an $s$-stage method one has $\Sigma_{\mkern1mu s} \mkern1.75mu {\pmb{b}} = 1$. On the diagram the points $\bigl( {c_{\mkern0.5mu i}} ,\mkern2mu {\Sigma_{\mkern1mu i - 1} \mkern1.5mu {\pmb{b}}} \bigr)$ and $\bigl( {c_{\mkern0.5mu i}} ,\mkern2mu {\Sigma_{\mkern1mu i} \mkern1.5mu {\pmb{b}}} \bigr)$, where $1 \le i \le s$, are connected by a thick line (with white filling if ${\Sigma_{\mkern1mu i} \mkern1.5mu {\pmb{b}}} < {\Sigma_{\mkern1mu i - 1} \mkern1.5mu {\pmb{b}}}$). To easier follow the progression of the stages, the points $\bigl( {c_{\mkern0.5mu i}} ,\mkern2mu {\Sigma_{\mkern1mu i} \mkern1.5mu {\pmb{b}}} \bigr)$ and $\bigl( c_{\mkern0.5mu i + 1} ,\mkern2mu {\Sigma_{\mkern1mu i} \mkern1.5mu {\pmb{b}}} \bigr)$, where $1 \le i < s$, are connected by dotted lines.
• Figure 3: A solution of Euler's equations with moments of inertia $I_1 = 1$, $I_2 = 2$, and $I_3 = 3$: $(\omega_{\mkern1mu 1} ,\mkern2mu \omega_{\mkern1.5mu 2} ,\mkern2mu \omega_{\mkern1.5mu 3}) ( \mkern0.25mu t \mkern-0.5mu ) = (12 \mkern1mu \textrm{cn} ,\mkern2mu 12 \mkern1mu \textrm{sn} ,\mkern2mu 7 \mkern1mu \textrm{dn}) \bigl( 7 \mkern0.5mu t ,\mkern2mu \frac{48}{49} \bigr)$, where $\textrm{cn}$, $\textrm{sn}$, and $\textrm{dn}$ are the Jacobi elliptic functions. The motion is periodic with the period $4 K \bigl( \frac{48}{49} \bigr) / 7 = 1.9109...$, where $K(m) = {\int_{\mkern1mu 0}^{\mkern1mu \pi/2}} {\rm d} \theta \mkern1.5mu / \bigl( 1 - m \mkern1mu \sin^2 \theta \bigr){}^{1/2}$ is the complete elliptic integral of the $1^{\textrm{st}}$ kind. The circular and elliptic cylinders correspond to the quadratic invariants $Q_1({\pmb{\omega}}) = \omega_{\mkern1mu 1}^{\mkern1mu 2} + \omega_{\mkern1.5mu 2}^{\mkern1mu 2} = 144$ and $Q_2({\pmb{\omega}}) = \omega_{\mkern1.5mu 2}^{\mkern1mu 2} + 3 \omega_{\mkern1.5mu 3}^{\mkern1mu 2} = 147$, respectively. The thin lines on them divide the period into $72$ equal parts. Points of the trajectory that are connected by vertical lines with $\omega_{\mkern1mu 1}$ and $\omega_{\mkern1.5mu 2}$ axes are $(12 ,\mkern2mu 0 ,\mkern2mu 7)$ and $(0 ,\mkern2mu 12 ,\mkern2mu 1)$. The alternating behavior of $\omega_{\mkern1.5mu 2}( \mkern0.25mu t \mkern-0.5mu )$ is due to the rotation around the intermediate principal axis being unstable Poi34, LaLi76.
• Figure 4: Growth of the error in $Q_1({\pmb{\tilde{\omega}}})$ and $Q_2({\pmb{\tilde{\omega}}})$ with time $t$. Panels from left to right: (1) RK4 (curves from top to near bottom) and GL4 (bottom curves); (2) AC36 (lower curves in groups of three), CLMR47 (middle curves), and CCRL47 (upper curves); (3) CV8 and eq. (\ref{['method_C']}) (lower and upper curves in pairs); and (4) eq. (\ref{['the_method']}). The step size is $h = s \mkern1mu h_1$, where $s$ is the number of stages. Upper solid, dotted, and dashed curves correspond to $h_1 = 2^{-7}$, $h_1 = 2^{-8}$, and $h_1 = 2^{-9}$, respectively. Next to upper solid or dotted curves correspond to $h_1 = 2^{-10}$ or $h_1 = 2^{-11}$, respectively, and so on, i.e., going to the next, with a different line style, curve down divides $h_1$ by $2$. Upper solid curves, upper dotted curves, and so on, on different panels correspond to the same value of $h_1$, and thus to the same number of the r.h.s. function ${\pmb{f}}(t ,\mkern2mu {\pmb{x}})$ evaluations.
• Figure 5: The speed of the systematic drift of the Hamiltonian $\mathcal{H}$ as the function of the step size $h = s \mkern1mu h_1$, where $s$ is the number of stages. The number of the r.h.s. function ${\pmb{f}}(t ,\mkern2mu {\pmb{x}})$ evaluations is equal to the time duration of the simulation (which is $10020000$ here) divided by $h_1$, and with $h_1$ being fixed is independent of the method. The curves correspond to RK4, AC36, and eq. (\ref{['the_method']}) (solid curves from top to bottom), CLMR47 and CV8 (upper and lower dotted curves), and CCRL47, eq. (\ref{['method_C']}), and GL4 (dashed curves from top to bottom). Three thin solid lines on this log-log plot have slopes $5$, $7$, and $9$, which corresponds to the speed of the drift being proportional to $h^5$, $h^7$, and $h^9$. The weighted moving average of the Hamiltonian that was used is $\langle \mathcal{H} \rangle ( \mkern0.25mu t \mkern-0.5mu ) = \sum_{\mkern1.5mu n \mkern1.5mu \in \mathcal{S} ( \mkern0.25mu t \mkern-0.5mu )} \sin^2 \bigl( \pi (n h - t) / 20000 \bigr) \mkern1mu \mathcal{H} \bigl( \tilde{p}(n h) ,\mkern2mu \tilde{x}(n h) \bigr) \mkern1mu / \mkern1mu \sum_{\mkern1.5mu n \mkern1.5mu \in \mathcal{S} ( \mkern0.25mu t \mkern-0.5mu )} \sin^2 \bigl( \pi (n h - t) / 20000 \bigr)$, where $\mathcal{S} ( \mkern0.25mu t \mkern-0.5mu ) = \bigl\{ n \mkern2.5mu \vert \mkern3.5mu t < n h < t + 20000 \bigr\}$. The simulations were run at the University of Arizona High Performance Computing center.