Hamilton–Pontryagin spectral-collocation methods for the orbit propagation

Abstract

According to the discrete Hamilton–Pontryagin variational principle, we construct a class of variational integrators in the real vector spaces and extend to the Lie groups for the left-trivialized Lagrangian mechanical systems by employing the spectral-collocation method to discretize the corresponding Lagrangian and kinematic constraints. The constructed framework can be transformed easily to the well-known symplectic partitioned Runge–Kutta methods and the higher order symplectic partitioned Lie Group methods by choosing same interpolation nodes and quadrature points. Two numerical experiments about the orbit propagation of Kepler two-body system and the rigid-body flow propagation of a free rigid body are conducted, respectively. The simulating results reveal that the constructed update schemes can possess simultaneously the excellent exponent convergence rates of spectral methods and the attractive long-term structure-preserving properties of geometric numerical algorithms.

Graphic abstract

Appendices

Appendix-A

We know that Sect. 3.2 takes a “reverse order” interpolation style to construct the HP spectral-collocation variational integrators in Eq. (14). We would give the results of taking a “positive order” style in the following.

Firstly, adopt the Lagrange interpolation polynomial to approximate the variable $q(t)$ on the $k{\text{th}}$ time subinterval $[kh,(k + 1)h]$, i.e.,

$$q_{\text{d}} (t_{\sigma } ) = \sum\limits_{j = 0}^{s} {l_{j,s} (\sigma )q_{k}^{j} } ,$$

(A1)

where ${{t_{\sigma } = (\sigma + 1)h} \mathord{\left/ {\vphantom {{t_{\sigma } = (\sigma + 1)h} 2}} \right. \kern-\nulldelimiterspace} 2}{\text{ belongs to }}[0,h]$, $l_{j,s} (\sigma )$ is the interpolation basis function, and $q_{k}^{j}$ represent the control nodes in interval $[0,h]$.

Then we can obtain the approximation curve of velocity term $\dot{q}_{\text{d}} (t)$ by using the following differential formula

$$\dot{q}_{\text{d}} (t_{\sigma } ) = \frac{2}{h}\sum\limits_{j = 0}^{s} {\dot{l}_{j,s} (\sigma )q_{k}^{j} } ,$$

(A2)

where $\dot{l}_{j,s} (\sigma )$ is the differential of classic Lagrange interpolation polynomial. One can choose the Barycentric Lagrange interpolation, which possesses a better numerical stability [25].

Given a Lagrangian $L:TQ \to {\mathbb{R}}$ with second order continuous partial derivatives, we then define a discrete path space ${\mathcal{C}}_{\text{d}} = \left\{ {\left[ {q,p,\left( {Q^{i} ,V^{i} ,p^{i} } \right)_{i = 1}^{m} } \right]_{\text{d}} :\left( {t_{k} } \right)_{k = 0}^{N} \to T^{ * } Q \times } \right.$ $\left. {\left( {TQ \times T^{ * } Q} \right)^{m} \left| {q(t_{0} ) = q_{0} ,q(t_{N} ) = q_{N} } \right.} \right\}$. Then the discrete HP action sum ${\mathfrak{S}}_{{\text{d}}} :{\mathcal{C}}_{\text{d}} (q_{0} ,q_{N} ) \to {\mathbb{R}}$ takes the form as following

$$\begin{aligned} {\mathfrak{S}}_{{\text{d}}} & = \sum\limits_{k = 0}^{N - 1} \left\{ \frac{h}{2}\sum\limits_{i = 1}^m {{\omega _i}L(Q_k^i,V_k^i)} + \frac{h}{2}\sum\limits_{i = 1}^m {\left\langle {p_k^i,V_k^i - \frac{2}{h}\sum\limits_{j = 0}^s {{{\dot l}_{j,s}}({\sigma _i})q_k^j} } \right\rangle } \right. \\ & \quad + {\left\langle {{p_k},q_k^ + - {q_k}} \right\rangle + \left\langle {{p_{k + 1}},{q_{k + 1}} - q_{k + 1}^ - } \right\rangle } \Bigg\}, \end{aligned}$$

(A3)

where $q_{k}^{ + }$ is the approximating point of $q_{k}$ from the “right” side, and $q_{k + 1}^{ - }$ is the approximating point of $q_{k + 1}$ from the “left” side in the $k{\text{th}}$ time interval $[t_{k} ,t_{k + 1} ]$.

According to the discrete HP variational principle, we should find a discrete curve $c_{\text{d}} \in {\mathcal{C}}_{\text{d}} (q_{0} ,q_{N} )$, which is a critical point of ${\mathfrak{S}}_{{\text{d}}} :{\mathcal{C}}_{\text{d}} (q_{0} ,q_{N} ) \to {\mathbb{R}}$, i.e.,${\mathbf{d}}{\mathfrak{S}}_{\text{d}} (c_{\text{d}} ) = 0$. Taking the variations of ${\mathfrak{S}}_{\text{d}}$ yields

$$\begin{gathered} {\updelta }{\mathfrak{S}}_{{\text{d}}} = \sum\limits_{k = 0}^{N - 1} {\left\{ {\frac{h}{2}\sum\limits_{i = 1}^{m} {\omega_{i} \left[ {\frac{\partial L}{{\partial q}}(Q_{k}^{i} ,V_{k}^{i} ) \cdot {\updelta }Q_{k}^{i} + \frac{\partial L}{{\partial \dot{q}}}(Q_{k}^{i} ,V_{k}^{i} ) \cdot {\updelta }V_{k}^{i} } \right]} } \right.} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \;\frac{h}{2}\;\sum\limits_{i = 1}^{m} {\left\langle {{\updelta }p_{k}^{i} ,V_{k}^{i} - \frac{2}{h}\sum\limits_{j = 0}^{s} {\dot{l}_{j,s} (\sigma_{i} )q_{k}^{j} } } \right\rangle } \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \;\frac{h}{2}\sum\limits_{i = 1}^{m} {\left\langle {p_{k}^{i} ,{\updelta }V_{k}^{i} - \frac{2}{h}\sum\limits_{j = 0}^{s} {\dot{l}_{j,s} (\sigma_{i} ) \cdot {\updelta }q_{k}^{j} } } \right\rangle } \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left\langle {{\updelta }p_{k} ,q_{k}^{ + } - q_{k} } \right\rangle + \left\langle {p_{k} ,{\updelta }q_{k}^{ + } } \right\rangle \; \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + {\left\langle {{\updelta }p_{k + 1} ,q_{k + 1} - q_{k + 1}^{ - } } \right\rangle + \left\langle { - p_{k + 1} ,{\updelta }q_{k + 1}^{ - } } \right\rangle } \Bigg\}. \hfill \\ \end{gathered}$$

(A4)

Substituting expression $Q_{k}^{i} = \sum\nolimits_{j = 0}^{s} {l_{j,s} (\sigma_{i} )q_{k}^{j} }$, $q_{k}^{ + } = q_{k}^{0}$, and $q_{k + 1}^{ - } = q_{k}^{s}$ into above equation yields

$$\begin{gathered} {\updelta }{\mathfrak{S}}_{{\text{d}}} = \sum\limits_{k = 0}^{N - 1} {\left\{ {\frac{h}{2}\sum\limits_{i = 1}^{m} {\omega_{i} \left[ {\frac{\partial L}{{\partial q}}(Q_{k}^{i} ,V_{k}^{i} )\sum\limits_{j = 0}^{s} {l_{j,s} (\sigma_{i} ) \cdot {\updelta }q_{k}^{j} } } \right]} } \right.} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \frac{h}{2}\sum\limits_{i = 1}^{m} {\omega_{i} \frac{\partial L}{{\partial \dot{q}}}(Q_{k}^{i} ,V_{k}^{i} ) \cdot {\updelta }V_{k}^{i} } + \frac{h}{2}\sum\limits_{i = 1}^{m} {\left\langle {p_{k}^{i} ,{\updelta }V_{k}^{i} } \right\rangle } \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \frac{h}{2}\sum\limits_{i = 1}^{m} {\left\langle {{\updelta }p_{k}^{i} ,V_{k}^{i} - \frac{2}{h}\sum\limits_{j = 0}^{s} {\dot{l}_{j,s} (\sigma_{i} )q_{k}^{j} } } \right\rangle } \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \frac{h}{2}\sum\limits_{i = 1}^{m} {\left\langle { - \frac{2}{h}\sum\limits_{j = 0}^{s} {\dot{l}_{j,s} (\sigma_{i} )p_{k}^{i} } ,{\updelta }q_{k}^{j} } \right\rangle } \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left\langle {{\updelta }p_{k} ,q_{k}^{0} - q_{k} } \right\rangle + \left\langle {p_{k} ,{\updelta }q_{k}^{0} } \right\rangle \; \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + {\left\langle {{\updelta }p_{k + 1} ,q_{k + 1} - q_{k}^{s} } \right\rangle + \left\langle { - p_{k + 1} ,{\updelta }q_{k}^{s} } \right\rangle } \Bigg\}.\; \hfill \\ \end{gathered}$$

(A5)

Collecting the terms with same variations and eliminating $p_{k}^{i}$ terms gain

$$\begin{aligned} {\updelta }{\mathfrak{S}}_{{\text{d}}} &= \sum\limits_{k = 0}^{N - 1} {\left\{ {\sum\limits_{i = 1}^{m} {\sum\limits_{j = 1}^{s - 1} {\left\langle {\frac{h}{2}\omega_{i} l_{j,s} (\sigma_{i} )\frac{\partial L}{{\partial q}}(Q_{k}^{i} ,V_{k}^{i} )} \right.} } } \right.} \\ &\qquad \qquad + \left. {\omega_{i} \dot{l}_{j,s} (\sigma_{i} )\frac{\partial L}{{\partial \dot{q}}}(Q_{k}^{i} ,V_{k}^{i} ),{\updelta }q_{k}^{j} } \right\rangle \\ & \qquad + \;\sum\limits_{i = 1}^{m} {\left\langle {\frac{h}{2}\omega_{i} l_{0,s} (\sigma_{i} )\frac{\partial L}{{\partial q}}(Q_{k}^{i} ,V_{k}^{i} )} \right.} \\ &\qquad + \left. {\omega_{i} \dot{l}_{0,s} (\sigma_{i} )\frac{\partial L}{{\partial \dot{q}}}(Q_{k}^{i} ,V_{k}^{i} ) + p_{k} ,{\updelta }q_{k}^{0} } \right\rangle \\ &\qquad + \sum\limits_{i = 1}^{m} {\left\langle {\frac{h}{2}\omega_{i} l_{s,s} (\sigma_{i} )\frac{\partial L}{{\partial q}}(Q_{k}^{i} ,V_{k}^{i} )} \right.} \\ & \qquad \left. { + \omega_{i} \dot{l}_{s,s} (\sigma_{i} )\frac{\partial L}{{\partial \dot{q}}}(Q_{k}^{i} ,V_{k}^{i} ) - p_{k + 1} ,{\updelta }q_{k}^{s} } \right\rangle \\ & \qquad + \;\frac{h}{2}\sum\limits_{i = 1}^{m} {\left\langle {{\updelta }p_{k}^{i} ,V_{k}^{i} - \frac{2}{h}\sum\limits_{j = 0}^{s} {\dot{l}_{j,s} (\sigma_{i} )q_{k}^{j} } } \right\rangle } \\ &\qquad + \;\left\langle {{\updelta }p_{k} ,q_{k}^{0} - q_{k} } \right\rangle \\ &\qquad + {\left\langle {{\updelta }p_{k + 1} ,q_{k + 1} - q_{k}^{s} } \right\rangle } \Bigg\}.\; \\ \end{aligned}$$

(A6)

Finally, one can arrive at a set of discrete update schemes, which can be written in a compact form as

$$\left\{ {\begin{array}{*{20}l} {q_{k}^{0} = q_{k} ,} \hfill \\ {Q_{k}^{i} = \sum\limits_{j = 0}^{s} {l_{j,s} (\sigma_{i} )q_{k}^{j} } ,} \hfill \\ {V_{k}^{i} = \frac{2}{h}\sum\limits_{j = 0}^{s} {\dot{l}_{j,s} (\sigma_{i} )q_{k}^{j} } ,} \hfill \\ {p_{k} = - \frac{h}{2}\sum\limits_{i = 1}^{m} {\omega_{i} \left[ {l_{0,s} (\sigma_{i} )\frac{\partial L}{{\partial q}}(Q_{k}^{i} ,V_{k}^{i} ) + \frac{2}{h}\dot{l}_{0,s} (\sigma_{i} )\frac{\partial L}{{\partial \dot{q}}}(Q_{k}^{i} ,V_{k}^{i} )} \right]} ,} \hfill \\ {0 = \frac{h}{2}\sum\limits_{i = 1}^{m} {\omega_{i} \left[ {l_{j,s} (\sigma_{i} )\frac{\partial L}{{\partial q}}(Q_{k}^{i} ,V_{k}^{i} ) + \frac{2}{h}\dot{l}_{j,s} (\sigma_{i} )\frac{\partial L}{{\partial \dot{q}}}(Q_{k}^{i} ,V_{k}^{i} )} \right],} } \hfill \\ {p_{k + 1} = \frac{h}{2}\sum\limits_{i = 1}^{m} {\omega_{i} \left[ {l_{s,s} (\sigma_{i} )\frac{\partial L}{{\partial q}}(Q_{k}^{i} ,V_{k}^{i} ) + \frac{2}{h}\dot{l}_{s,s} (\sigma_{i} )\frac{\partial L}{{\partial \dot{q}}}(Q_{k}^{i} ,V_{k}^{i} )} \right]} ,} \hfill \\ {q_{k + 1} = q_{k}^{s} ,} \hfill \\ \end{array} } \right.$$

(A7)

for $i = 1,{ 2, } \cdots , \, m$; $j = 1,{ 2, } \cdots , \, s - 1$; and $k = 0,{ 1, } \cdots , \, N - 1$.

One can find that the discrete update framework in Eq. (A-7) is just the well-known spectral variational integrators presented in Refs. [3, 25] (i.e., Eq. (8–10) of Ref. [3] and Eq. (11a-11c) of Ref. [25], respectively), where the result is derived from the discrete variational mechanics, but here we achieve it from the HP perspective.

Appendix-B

When the Lagrangian is left-invariant, i.e., $L(g,v) = L(hg,hv)$ holds for all $h \in G$, then the reduced Lagrangian ${\ell }(\xi ) = L(e,\xi ) = L(g,g\xi )$ (with $h = g^{ - 1}$) can be obtained. One can derive the discrete spectral-collocation update schemes as described in Eq. (39).

The continuous reduced HP action integral can be written as

$${\mathcal{S}} = \int_{{t_{0} }}^{{t_{\text{f}} }} {{\ell }(\xi ){\text{d}}t} { + }\int_{{t_{0} }}^{{t_{\text{f}} }} {\left\langle {\mu ,g^{ - 1} \dot{g} - \xi } \right\rangle {\text{d}}t} .$$

(B1)

The approximation strategy is similar to the one in Sect. 3.3. We define the discrete reduced path space as ${\mathcal{C}}_{\text{d}} = \left\{ {\left[ {g,\Lambda ,\left( {\Theta^{i} ,\Xi^{i} } \right)_{i = 1}^{s} ,\left( {\lambda^{j} } \right)_{j = 1}^{m} } \right]_{\text{d}} :\left( {t_{k} } \right)_{k = 0}^{N} \to \left( {G \times {\mathfrak{g}}^{ * } } \right) \times } \right.$$\left. {\left( {{\mathfrak{g}} \times {\mathfrak{g}}^{ * } } \right)^{s} \times \left( {{\mathfrak{g}}^{ * } } \right)^{m} \left| {g(t_{0} ) = g_{0} ,g(t_{N} ) = g_{N} } \right.} \right\}.$

Then the discrete reduced HP action sum ${\mathcal{S}}_{{\text{d}}} :{\mathcal{C}}_{\text{d}} \to {\mathbb{R}}$ takes the following form

$$\begin{gathered} {\mathcal{S}}_{{\text{d}}} = \sum\limits_{{k = 0}}^{{N - 1}} {\frac{h}{2}\left\{ {\sum\limits_{{j = 1}}^{m} {\omega _{j} {\ell }(\Xi _{k}^{j} )} + \sum\limits_{{j = 1}}^{m} {\left\langle {\lambda _{k}^{j} ,{{2\Theta _{k}^{j} } \mathord{\left/ {\vphantom {{2\Theta _{k}^{j} } h}} \right. \kern-\nulldelimiterspace} h} - \sum\limits_{{i = 1}}^{s} {{\mathcal{A}}_{i}^{j} V_{k}^{i} } } \right\rangle } + } \right.} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\left\langle {\Lambda _{{k + 1}} ,{{2\tau ^{{ - 1}} \left( {g_{k}^{{ - 1}} g_{{k + 1}} } \right)} \mathord{\left/ {\vphantom {{2\tau ^{{ - 1}} \left( {g_{k}^{{ - 1}} g_{{k + 1}} } \right)} h}} \right. \kern-\nulldelimiterspace} h} - \sum\limits_{{i = 1}}^{s} {\varpi _{i} V_{k}^{i} } } \right\rangle } \right\}. \hfill \\ \end{gathered}$$

(B2)

The discrete reduced HP principle states that ${\updelta }{\mathcal{S}}_{{\text{d}}} = 0$ with variations $\left( {{\updelta }g_{k} ,{\updelta }\Lambda_{k} ,{\updelta }\Theta_{k}^{i} ,{\updelta }\Xi_{k}^{i} ,{\updelta }\lambda_{k}^{j} } \right)$ and the fixed endpoint variations $\left( {{\updelta }g_{0} = {\updelta }g_{N} = 0} \right)$. Taking the variation of ${\mathcal{S}}_{{\text{d}}}$ yields

$$\begin{gathered} \delta {\mathcal{S}}_{d} = \sum\limits_{{k = 0}}^{{N - 1}} {\left\{ {\frac{h}{2}\sum\limits_{{j = 1}}^{m} {\omega _{j} \ell ^{\prime}(\Xi _{k}^{j} ) \cdot \delta \Xi _{k}^{j} } } \right.} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\; + \;\frac{h}{2}\sum\limits_{{j = 1}}^{m} {\left\langle {\delta \lambda _{k}^{j} ,{{2\Theta _{k}^{j} } \mathord{\left/ {\vphantom {{2\Theta _{k}^{j} } h}} \right. \kern-\nulldelimiterspace} h} - \sum\limits_{{i = 1}}^{s} {{\mathcal{A}}_{i}^{j} V_{k}^{i} } } \right\rangle } \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \sum\limits_{{j = 1}}^{m} {\left\langle {\lambda _{k}^{j} ,\delta \Theta _{k}^{j} } \right\rangle } - \;\frac{h}{2}\sum\limits_{{j = 1}}^{m} {\sum\limits_{{i = 1}}^{s} {\left\langle {{\mathcal{A}}_{i}^{j} \lambda _{k}^{j} ,\delta V_{k}^{i} } \right\rangle } } \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \frac{h}{2}\left\langle {\delta \Lambda _{{k + 1}} ,{{2\tau ^{{ - 1}} \left( {g_{k}^{{ - 1}} g_{{k + 1}} } \right)} \mathord{\left/ {\vphantom {{2\tau ^{{ - 1}} \left( {g_{k}^{{ - 1}} g_{{k + 1}} } \right)} h}} \right. \kern-\nulldelimiterspace} h} - \sum\limits_{{i = 1}}^{s} {\varpi _{i} V_{k}^{i} } } \right\rangle \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left\langle { - \left( {\text{d}\tau _{{{h \mathord{\left/ {\vphantom {h 2}} \right. \kern-\nulldelimiterspace} 2}\xi _{{k + 1}} }}^{{ - 1}} } \right)^{ * } \Lambda _{{k + 1}} ,\eta _{k} } \right\rangle \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left\langle {\left( {\text{d}\tau _{{ - {h \mathord{\left/ {\vphantom {h 2}} \right. \kern-\nulldelimiterspace} 2}\xi _{{k + 1}} }}^{{ - 1}} } \right)^{ * } \Lambda _{{k + 1}} ,\eta _{{k + 1}} } \right\rangle \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \;\left. {\frac{h}{2}\sum\limits_{{j = 1}}^{m} {\sum\limits_{{i = 1}}^{s} {\left\langle { - \varpi _{i} \Lambda _{{k + 1}} ,\delta V_{k}^{i} } \right\rangle } } } \right\}. \hfill \\ \end{gathered}$$

(B3)

In terms of expressions $\Xi_{k}^{j} = \sum\nolimits_{i = 1}^{s} {l_{i,s} (\sigma_{j} )\Xi_{k}^{i} }$, $\Theta_{k}^{j} = \sum\nolimits_{i = 1}^{s} {l_{i,s} (\sigma_{j} )\Theta_{k}^{i} }$, and $V_{k}^{i} { = d}\tau_{{ - \Theta_{k}^{i} }}^{ - 1} \Xi_{k}^{i}$, we conduct summation by parts with respect to $\eta_{k}$ and $\eta_{k + 1}$, then the above equation can be rewritten as

$$\begin{aligned} \delta {\mathcal{S}}_{d} & = \sum\limits_{{k = 0}}^{{N - 1}} {\left\{ {\frac{h}{2}\sum\limits_{{j = 1}}^{m} {\sum\limits_{{i = 1}}^{s} {l_{{i,s}} (\sigma _{j} )\omega _{j} \ell ^{\prime}(\Xi _{k}^{j} ) \cdot \delta \Xi _{k}^{i} } } } \right.} \\ &\qquad + \frac{h}{2}\sum\limits_{{j = 1}}^{m} {\left\langle {\delta \lambda _{k}^{j} ,{{2\Theta _{k}^{j} } \mathord{\left/ {\vphantom {{2\Theta _{k}^{j} } h}} \right. \kern-\nulldelimiterspace} h} - \sum\limits_{{i = 1}}^{s} {{\mathcal{A}}_{i}^{j} \text{d}\tau _{{ - \Theta _{k}^{i} }}^{{ - 1}} \Xi _{k}^{i} } } \right\rangle } \\ &\qquad + \sum\limits_{{j = 1}}^{m} {\sum\limits_{{i = 1}}^{s} {\left\langle {l_{{i,s}} (\sigma _{j} )\lambda _{k}^{j} ,\delta \Theta _{k}^{i} } \right\rangle } } \\ &\qquad + \frac{h}{2}\sum\limits_{{j = 1}}^{m} {\sum\limits_{{i = 1}}^{s} {\left\langle { -{\mathcal{A}}_{i}^{j} D\left( {\text{d}\tau _{{ - \Theta _{k}^{i} }}^{{ - 1}} \Xi _{k}^{i} } \right)^{ * } \lambda _{k}^{j} ,\delta \Theta _{k}^{i} } \right\rangle } } \\ &\qquad + \frac{h}{2}\sum\limits_{{j = 1}}^{m} {\sum\limits_{{i = 1}}^{s} {\left\langle { - {\mathcal{A}}_{i}^{j} \left( {\text{d}\tau _{{ - \Theta _{k}^{i} }}^{{ - 1}} } \right)^{ * } \lambda _{k}^{j} ,\delta \Xi _{k}^{i} } \right\rangle } } \\ &\qquad + \frac{h}{2}\left\langle {\delta \Lambda _{{k + 1}} ,{{2\tau ^{{ - 1}} \left( {g_{k}^{{ - 1}} g_{{k + 1}} } \right)} \mathord{\left/ {\vphantom {{2\tau ^{{ - 1}} \left( {g_{k}^{{ - 1}} g_{{k + 1}} } \right)} h}} \right. \kern-\nulldelimiterspace} h} - \sum\limits_{{i = 1}}^{s} {\varpi _{i} \text{d}\tau _{{ - \Theta _{k}^{i} }}^{{ - 1}} \Xi _{k}^{i} } } \right\rangle \\ &\qquad + \;\left\langle {\left( {\text{d}\tau _{{ - h/2\xi _{k} }}^{{ - 1}} } \right)^{ * } \Lambda _{k} - \left( {\text{d}\tau _{{h/2\xi _{{k + 1}} }}^{{ - 1}} } \right)^{ * } \Lambda _{{k + 1}} ,\eta _{k} } \right\rangle \\ & \qquad + \frac{h}{2}\sum\limits_{{j = 1}}^{m} {\sum\limits_{{i = 1}}^{s} {\left\langle { - \varpi _{i} D\left( {\text{d}\tau _{{ - \Theta _{k}^{i} }}^{{ - 1}} \Xi _{k}^{i} } \right)^{ * } \Lambda _{{k + 1}} ,\delta \Theta _{k}^{i} } \right\rangle } } \\ &\qquad \left. { + \frac{h}{2}\sum\limits_{{j = 1}}^{m} {\sum\limits_{{i = 1}}^{s} {\left\langle { - \varpi _{i} \left( {\text{d}\tau _{{ - \Theta _{k}^{i} }}^{{ - 1}} } \right)^{ * } \Lambda _{{k + 1}} ,\delta \Xi _{k}^{i} } \right\rangle } } } \right\}. \\ \end{aligned}$$

(B4)

Introducing the external stage momenta $\mu _{k} = \left( {\text{d}\tau _{{ - {h \mathord{\left/ {\vphantom {h 2}} \right. \kern-\nulldelimiterspace} 2}\xi _{k} }}^{{ - 1}} } \right)^{ * } \Lambda _{k}$ and collecting same variation terms would yield

$$\left\{ {\begin{array}{*{20}l} {\sum\limits_{{j = 1}}^{m} {l_{{i,s}} (\sigma _{j} )\omega _{j} \ell ^{\prime}(\Xi _{k}^{j} )} = \left( {\text{d}\tau _{{ - \Theta _{k}^{i} }}^{{ - 1}} } \right)^{ * } \left( {\sum\limits_{{j = 1}}^{m} {{\mathcal{A}}_{i}^{j} \lambda _{k}^{j} } + \varpi _{i} \Lambda _{{k + 1}} } \right),} \hfill \\ {\sum\limits_{{j = 1}}^{m} {l_{{i,s}} (\sigma _{j} )\lambda _{k}^{j} } = \frac{h}{2}D\left( {\text{d}\tau _{{ - \Theta _{k}^{i} }}^{{ - 1}} \Xi _{k}^{i} } \right)^{ * } \left( {\sum\limits_{{j = 1}}^{m} {{\mathcal{A}}_{i}^{j} \lambda _{k}^{j} } + \varpi _{i} \Lambda _{{k + 1}} } \right),} \hfill \\ {\Theta _{k}^{j} = \frac{h}{2}\sum\limits_{{i = 1}}^{s} {{\mathcal{A}}_{i}^{j} \text{d}\tau _{{ - \Theta _{k}^{i} }}^{{ - 1}} \Xi _{k}^{i} } ,} \hfill \\ {\xi _{{k + 1}} = \sum\limits_{{i = 1}}^{s} {\varpi _{i} \text{d}\tau _{{ - \Theta _{k}^{i} }}^{{ - 1}} \Xi _{k}^{i} } ,} \hfill \\ {g_{{k + 1}} = g_{k} \tau \left( {\frac{h}{2}\xi _{{k + 1}} } \right),} \hfill \\ {\mu _{{k + 1}} = \left( {A\text{d}\tau _{{{h \mathord{\left/ {\vphantom {h 2}} \right. \kern-\nulldelimiterspace} 2}\xi _{{k + 1}} }} } \right)^{ * } \mu _{k} .} \hfill \\ \end{array} } \right.$$

(B5)

Finally, we can obtain a compact form as following

$$\left\{ {\begin{array}{*{20}l} {\sum\limits_{{j = 1}}^{m} {l_{{i,s}} (\sigma _{j} )\omega _{j} \ell ^{\prime}(\Xi _{k}^{j} )} = \left( {\text{d}\tau _{{ - \Theta _{k}^{i} }}^{{ - 1}} } \right)^{ * } \left( {\sum\limits_{{j = 1}}^{m} {{\mathcal{A}}_{i}^{j} \lambda _{k}^{j} } + \varpi _{i} \Lambda _{{k + 1}} } \right),} \hfill \\ {\sum\limits_{{j = 1}}^{m} {l_{{i,s}} (\sigma _{j} )\lambda _{k}^{j} } = \frac{h}{2}D\left( {\text{d}\tau _{{ - \Theta _{k}^{i} }}^{{ - 1}} \Xi _{k}^{i} } \right)^{ * } \left( {\sum\limits_{{j = 1}}^{m} {{\mathcal{A}}_{i}^{j} \lambda _{k}^{j} } + \varpi _{i} \Lambda _{{k + 1}} } \right),} \hfill \\ {\Theta _{k}^{j} = \frac{h}{2}\sum\limits_{{i = 1}}^{s} {{\mathcal{A}}_{i}^{j} \text{d}\tau _{{ - \Theta _{k}^{i} }}^{{ - 1}} \Xi _{k}^{i} } ,} \hfill \\ {\xi _{{k + 1}} = \sum\limits_{{i = 1}}^{s} {\varpi _{i} \text{d}\tau _{{ - \Theta _{k}^{i} }}^{{ - 1}} \Xi _{k}^{i} } ,} \hfill \\ {g_{{k + 1}} = g_{k} \tau \left( {\frac{h}{2}\xi _{{k + 1}} } \right),} \hfill \\ {\mu _{{k + 1}} = \left( {A\text{d}\tau _{{{h \mathord{\left/ {\vphantom {h 2}} \right. \kern-\nulldelimiterspace} 2}\xi _{{k + 1}} }} } \right)^{ * } \mu _{k} ,} \hfill \\ \end{array} } \right.$$

(B6)

for $i = 1,{ 2, } \cdots , \, s;\;j = 1,{ 2, } \cdots , \, m;\;{\text{and}}\;k = 0,{ 1, } \cdots , \, N - 1$. Even though they can be simply viewed as eliminating the terms $\left( {{{\partial {\ell }} \mathord{\left/ {\vphantom {{\partial {\ell }} {\partial g}}} \right. \kern-\nulldelimiterspace} {\partial g}}} \right)(G_{k}^{j} ,\Xi_{k}^{j} )$ in Eq. (39), they correspond to different Lagrangian mechanical systems.

Appendix-C

The first three-stage quadrature weights and integrating matrixes in Eqs. (9) and (10) for scheme 1 are presented in the following Table 1 for the Gauss–Legendre points and the Chebyshev points, respectively. Note that these data are applicable for the time interval $[ - 1,\;1]$.

Table 1 First three-stage Butcher tableaux for HPSCVI in scheme 1