The Lie group Euler methods of multibody system dynamics with holonomic constraints

Abstract

The Euler methods on Lie group are developed for the differential–algebraic equations of multibody system dynamics with holonomic constraints. The implicit Euler method is used to solve the differential–algebraic equations as Euler–Lagrange equations on Lie group with indices 1, 2, and 3 and the case of overdetermined differential–algebraic equations mixing with configuration space. The symplectic Euler method is used to solve the differential–algebraic equations as constrained Hamilton equations on Lie group. For the discrete mapping between Lie group and Lie algebra, the canonical coordinates of the second kind for implicit first-order Crouch–Grossman Euler methods of differential–algebraic equations are used. A single pendulum and a double pendulum in the space are used to verify the accuracy of the Lie group Euler methods.

Introduction

Geometric integration methods^1–4 have received considerable attention in computational multibody system dynamics due to their numerical stabilities. With the successful developments and applications of symplecticity and energy methods in simulations of dynamic systems, the investigation of Lie group methods has become another research focus.^5–10

The Lie group method^11,12 was originally employed in computational solid dynamics to realize specialized orthogonal group (orientation) preserving for numerical stability.^13,14 It had been studied systematically in 1990s for matrix ordinary differential equations (ODEs) using the canonical coordinates of the first kind (CCFK) and canonical coordinates of the second kind (CCSK), based on which the Munthe-Kaas ¹⁵ (MK) method and the Crouch–Grossman ¹⁶ (CG) method were proposed creatively. Making use of traditional Runge–Kutta algorithm, the MK and CG methods have been applied to computational rigid-body dynamics without constraint-described matrix ODEs^17–20 or discrete variational principles.^21–24

Typical equations of motion of multibody systems include coupled ODEs and algebraic equations of constraints, that is, differential–algebraic equations (DAEs). ²⁵ Various algorithms solving DAEs in vector spaces have been designed successfully in the past three decades,^26,27 but how to deal with DAEs in Lie group space is a challenging problem due to the complex group features. Based on the local mapping between Lie group and its corresponding Lie algebra which is isomorphic with a vector space, Lie group methods are usually designed on Lie algebra spaces using traditional algorithms first, then the Lie groups are reconstructed via local mapping. Arnold and coworkers^5–7 have extended the generalized-alpha method to solve DAEs on Lie group with indices 3 and 2. Müller and coworkers^8–10 have extended the explicit Runge–Kutta method with constraint projection method to solve DAEs with index 1 on Lie group.

In practice, different forms of DAEs of multibody system dynamics demonstrate different numerical properties even if they are solved via the same method.^25–27 In this article, we will use implicit Euler method to solve the DAEs as Euler–Lagrange equations on Lie group with indices 1, 2, and 3 and the case of overdetermined differential–algebraic equations (ODAEs) mixing with configuration space, and also we use symplectic Euler method to solve the DAEs as constrained Hamilton equations on Lie group. For the discrete mapping between Lie group and Lie algebra, we use CCSK coordinates for implicit first-order Crouch–Grossman Euler (CGE) methods of DAEs. The implicit Euler method is a special case of implicit multi-step methods and implicit Runge–Kutta methods with a large stability region, and the implicit Euler methods of different DAEs on Lie group are significant for the design of relevant higher order algorithms and large-scale system simulations.

This article is organized as follows. Section “Motion equations of multibody systems on Lie group” introduces five types of DAEs on Lie group of multibody system dynamics in different indices. In section “Implicit Euler Lie group methods of DAEs,” we design five implicit Euler methods on Lie group for DAEs in section “Motion equations of multibody systems on Lie group.” In section “Numerical examples,” we give some numerical examples to test and compare the algorithms in section “Implicit Euler Lie group methods of DAEs.” Concluding remarks are given finally.

Motion equations of multibody systems on Lie group

For a typical multibody system with $n$ rigid bodies and $l$ constraints, let ${\mathit{q}}_i=(\begin{array}{cc}\hfill {\mathit{x}}_i\hfill & \hfill {\mathit{R}}_i\hfill \end{array}),(i=1,\dots ,n)$ denote the absolute position of mass center and orientation of the ith body, and the bases of each body-fixed frame are parallel to the three principal axes of inertia. Also let the corresponding velocity be ${\mathit{V}}_i=(\begin{array}{cc}\hfill {\mathit{v}}_{\mathit{i}}\hfill & \hfill {\mathbf{\Omega }}_i\hfill \end{array}),(i=1,\dots ,n)$ , where ${\mathit{v}}_i$ is the absolute translation velocity and ${\mathbf{\Omega }}_i$ is the absolute angular velocity expressed in body-fixed frame; ${\stackrel{·}{\mathit{R}}}_i={\mathit{R}}_i{\tilde{\mathbf{\Omega }}}_i$ , where ${\mathit{R}}_i$ is the special orthogonal group SO(3) and ${\tilde{\mathbf{\Omega }}}_i$ is its Lie algebra which is just a skew-symmetric matrix of ${\mathbf{\Omega }}_i$ . In order to express the formulations and the algorithms concisely, we introduce some compact notations as follows.

Let the position vectors of the system $\mathit{x}=[\begin{array}{cccc}\hfill {\mathit{x}}_1^{\mathrm{T}}\hfill & \hfill {\mathit{x}}_2^{\mathrm{T}}\hfill & \hfill \cdots \hfill & \hfill {\mathit{x}}_n^{\mathrm{T}}\hfill \end{array}{]}^{\mathrm{T}}$ , the translation velocity vectors $\mathit{v}=[\begin{array}{cccc}\hfill {\mathit{v}}_1^{\mathrm{T}}\hfill & \hfill {\mathit{v}}_2^{\mathrm{T}}\hfill & \hfill \cdots \hfill & \hfill {\mathit{v}}_n^{\mathrm{T}}\hfill \end{array}{]}^{\mathrm{T}}$ , the orientation velocity vectors $\mathit{R}=\mathrm{diag}(\begin{array}{cccc}\hfill {\mathit{R}}_1\hfill & \hfill {\mathit{R}}_2\hfill & \hfill \cdots \hfill & \hfill {\mathit{R}}_n\hfill \end{array})$ , the angular velocity vectors $\mathbf{\Omega }=[\begin{array}{cccc}\hfill {\mathbf{\Omega }}_1^{\mathrm{T}}\hfill & \hfill {\mathbf{\Omega }}_2^{\mathrm{T}}\hfill & \hfill \cdots \hfill & \hfill {\mathbf{\Omega }}_n^{\mathrm{T}}\hfill \end{array}{]}^{\mathrm{T}}$ , and the related Lie algebras $\tilde{\mathbf{\Omega }}=\mathrm{diag}(\begin{array}{cccc}\hfill {\tilde{\mathbf{\Omega }}}_1\hfill & \hfill {\tilde{\mathbf{\Omega }}}_2\hfill & \hfill \cdots \hfill & \hfill {\tilde{\mathbf{\Omega }}}_n\hfill \end{array})$ . We can introduce the compact forms $\mathit{q}=(\begin{array}{cc}\hfill \mathit{x}\hfill & \hfill \mathit{R}\hfill \end{array})$ , $\stackrel{·}{\mathit{q}}=(\begin{array}{cc}\hfill \stackrel{·}{\mathit{x}}\hfill & \hfill \mathit{R}\tilde{\mathbf{\Omega }}\hfill \end{array})$ , and $\mathit{V}=[\begin{array}{cc}\hfill {\mathit{v}}^{\mathrm{T}}\hfill & \hfill {\mathbf{\Omega }}^{\mathrm{T}}\hfill \end{array}{]}^{\mathrm{T}}$ .

Let ${\mathit{q}}_i=(\begin{array}{cc}\hfill {\mathit{x}}_i\hfill & \hfill {\mathit{R}}_i\hfill \end{array})$ , ${\stackrel{·}{\mathit{q}}}_i=(\begin{array}{cc}\hfill {\stackrel{·}{\mathit{x}}}_i\hfill & \hfill {\stackrel{·}{\mathit{R}}}_i\hfill \end{array})=(\begin{array}{cc}\hfill {\mathit{v}}_i\hfill & \hfill {\mathit{R}}_i{\tilde{\mathbf{\Omega }}}_i\hfill \end{array})$ , ${\mathit{V}}_i=[\begin{array}{cc}\hfill {\mathit{v}}_i^{\mathrm{T}}\hfill & \hfill {\mathbf{\Omega }}_i^{\mathrm{T}}\hfill \end{array}]$ , $(i=1,\dots ,n)$ denote the poses and velocities individually. We can introduce other compact forms $\mathit{q}=(\begin{array}{cccc}\hfill {\mathit{q}}_1\hfill & \hfill {\mathit{q}}_2\hfill & \hfill \cdots \hfill & \hfill {\mathit{q}}_n\hfill \end{array})$ , $\stackrel{·}{\mathit{q}}=(\begin{array}{cccc}\hfill {\stackrel{·}{\mathit{q}}}_1\hfill & \hfill {\stackrel{·}{\mathit{q}}}_2\hfill & \hfill \cdots \hfill & \hfill {\stackrel{·}{\mathit{q}}}_n\hfill \end{array})$ , and $\mathit{V}=[\begin{array}{cccc}\hfill {\mathit{V}}_1^{\mathrm{T}}\hfill & \hfill {\mathit{V}}_2^{\mathrm{T}}\hfill & \hfill \cdots \hfill & \hfill {\mathit{V}}_n^{\mathrm{T}}\hfill \end{array}{]}^{\mathrm{T}}$ .

Similarly, the inertia matrix of the system can be expressed in the following two forms: $\mathit{M}=\mathrm{diag}(\begin{array}{cccc}\hfill {\mathit{M}}_1\hfill & \hfill {\mathit{M}}_2\hfill & \hfill \cdots \hfill & \hfill {\mathit{M}}_n\hfill \end{array})$ and $\mathit{M}=\mathrm{diag}(\begin{array}{cc}\hfill \mathit{m}\hfill & \hfill \mathit{J}\hfill \end{array})$ , where ${\mathit{M}}_i=\left[\begin{array}{cc}\hfill m_i\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill J_i\hfill \end{array}\right]$ , $(i=1,\dots ,n)$ , $\mathit{m}=\mathrm{diag}(\begin{array}{cccc}\hfill m_1\hfill & \hfill m_2\hfill & \hfill \cdots \hfill & \hfill m_n\hfill \end{array})$ , and $\mathit{J}=\mathrm{diag}(\begin{array}{cccc}\hfill J_1\hfill & \hfill J_2\hfill & \hfill \cdots \hfill & \hfill J_n\hfill \end{array})$ . $m_i$ and $J_i$ , $(i=1,\dots ,n)$ are the translation and rotation inertia matrices, respectively, of the body $i$ .

It will not cause different interpretations even if we use the two notations sometimes in the following. Based on these notations, the motion equations of multibody systems can be written as

$$\{\begin{array}{c}\hfill \stackrel{·}{\mathit{q}}=\left(\begin{array}{cc}\hfill \mathit{v}\hfill & \hfill \mathit{R}\tilde{\mathbf{\Omega }}\hfill \end{array}\right)\hfill \\ \multicolumn{1}{c}{\mathit{M}\stackrel{·}{\mathit{V}}+{\mathit{B}}^{\mathrm{T}}(\mathit{q})\mathbf{\lambda }=\mathit{Q}(\mathit{q},\mathit{V})}\\ \multicolumn{1}{c}{\mathbf{\Phi }(\mathit{q})=0}\end{array}$$

(1)

where $\mathit{Q}(\mathit{q},\mathit{V})=[\begin{array}{cccc}\hfill {\mathit{Q}}_1^{\mathrm{T}}(\mathit{q},\mathit{V})\hfill & \hfill {\mathit{Q}}_2^{\mathrm{T}}(\mathit{q},\mathit{V})\hfill & \hfill \cdots \hfill & \hfill {\mathit{Q}}_n^{\mathrm{T}}(\mathit{q},\mathit{V})\hfill \end{array}{]}^{\mathrm{T}}$ is the generalized force vector, ${\mathit{Q}}_i(\mathit{q},\mathit{V})=\left[\begin{array}{c}\hfill 0\hfill \\ \multicolumn{1}{c}{-{\tilde{\mathbf{\Omega }}}_i{\mathit{J}}_i{\mathbf{\Omega }}_i}\end{array}\right]+{\mathit{Q}}_{\mathrm{ie}},(i=1,\dots ,n)$ , $Q_{\mathrm{ie}}$ is the generalized external force vector, $\mathbf{\Phi }(\mathit{q})=[\begin{array}{cccc}\hfill \mathbf{\Phi }(\mathit{q})\hfill & \hfill \mathbf{\Phi }(\mathit{q})\hfill & \hfill \cdots \hfill & \hfill \mathbf{\Phi }(\mathit{q})\hfill \end{array}{]}^{\mathrm{T}}=0$ is the scleronomic constraint equation, $\mathbf{\lambda }$ is the related Lagrangian multiplier vector, and $\mathit{B}(\mathit{q})$ is its Jacobian which can be illustrated as

$$\stackrel{·}{\mathbf{\Phi }}(\mathit{q})={\mathbf{\Phi }}_{\mathit{x}}\stackrel{·}{\mathit{x}}+〈\begin{array}{cc}\hfill {\mathbf{\Phi }}_{\mathit{R}}\hfill & \hfill \stackrel{·}{\mathit{R}}\hfill \end{array}〉\equiv \mathit{B}(\mathit{q})\mathit{V}=0$$

(2)

$$\stackrel{··}{\mathbf{\Phi }}(\mathit{q})\equiv \mathit{B}(\mathit{q})\stackrel{·}{\mathit{V}}+\stackrel{·}{\mathit{B}}(\mathit{q})\mathit{V}=0$$

(3)

According to the classic index theory, the equations given in equation (1) are DAEs of index 3, the stabilized index 2 form, ²⁸ also called as Gear–Gupta–Leimkuhler (GGL) form, of which is

$$\{\begin{array}{c}\stackrel{·}{\mathit{q}}=\left(\begin{array}{cc}\hfill \stackrel{·}{\mathit{x}}\prime \hfill & \hfill \mathit{R}\tilde{\Omega }\prime \hfill \end{array}\right)\hfill \\ \left[\begin{array}{c}\hfill \stackrel{·}{\mathit{x}}\prime \hfill \\ \multicolumn{1}{c}{\mathbf{\Omega }\mathbf{\text{'}}}\end{array}\right]=\mathit{V}-{\mathit{B}}^{\mathrm{T}}(\mathit{q})\mathbf{\mu }\hfill \\ \mathit{M}\stackrel{·}{\mathit{V}}+{\mathit{B}}^{\mathrm{T}}(\mathit{q})\mathbf{\lambda }=\mathit{Q}(\mathit{q},\mathit{V})\hfill \\ \mathbf{\Phi }(\mathit{q})=0\hfill \\ \mathit{B}(\mathit{q})\mathit{V}=0\hfill \end{array}$$

(4)

where $\mathbf{\mu }$ is another new Lagrangian multiplier. The index 1 form of equation (1) is

$$\{\begin{array}{c}\stackrel{·}{\mathit{q}}=\left(\begin{array}{cc}\hfill \mathit{v}\hfill & \hfill \mathit{R}\tilde{\mathbf{\Omega }}\hfill \end{array}\right)\hfill \\ \mathit{M}\stackrel{·}{\mathit{V}}+{\mathit{B}}^{\mathrm{T}}(\mathit{q})\mathbf{\lambda }=\mathit{Q}(\mathit{q},\mathit{V})\hfill \\ \mathit{B}(\mathit{q})\stackrel{·}{\mathit{V}}+\stackrel{·}{\mathit{B}}(\mathit{q})\mathit{V}=0\hfill \end{array}$$

(5)

which considers only the acceleration constraint equations. In order to fulfill the constraints at position and velocity level, equation (1) should augment equations (2) and (3), leading to the following ODAEs

$$\{\begin{array}{c}\stackrel{·}{\mathit{q}}=\left(\begin{array}{cc}\hfill \mathit{v}\hfill & \hfill \mathit{R}\tilde{\mathbf{\Omega }}\hfill \end{array}\right)\hfill \\ \mathit{M}\stackrel{·}{\mathit{V}}+{\mathit{B}}^{\mathrm{T}}(\mathit{q})\mathbf{\lambda }=\mathit{Q}(\mathit{q},\mathit{V})\hfill \\ \mathbf{\Phi }(\mathit{q})=0\hfill \\ \mathit{B}(\mathit{q})\mathit{V}=0\hfill \\ \mathit{B}(\mathit{q})\stackrel{·}{\mathit{V}}+\stackrel{·}{\mathit{B}}(\mathit{q})\mathit{V}=0\hfill \end{array}$$

(6)

Introducing canonical coordinates, the aforementioned DAEs can be transformed into the following Hamilton equations with constraints on Lie group

$$\{\begin{array}{c}\stackrel{·}{\mathit{q}}=\left(\begin{array}{cc}\hfill {\mathit{m}}^{-1}{\mathit{p}}_{\mathit{x}}\hfill & \hfill \mathit{R}{\left({\mathit{J}}^{-1}{\mathit{p}}_{\mathit{R}}\right)}^{~}\hfill \end{array}\right)\hfill \\ \stackrel{·}{\mathit{p}}+{\mathit{B}}^{\mathrm{T}}(\mathit{q})\mathbf{\lambda }=\mathit{Q}(\mathit{q},\mathit{p})\hfill \\ \mathbf{\Phi }(\mathit{q})=0\hfill \\ \mathit{B}(\mathit{q}){\mathit{M}}^{-1}\mathit{p}=0\hfill \end{array}$$

(7)

where the generalized momenta are $\mathit{p}=[\begin{array}{cc}\hfill {\mathit{p}}_{\mathit{x}}^{\mathrm{T}}\hfill & \hfill {\mathit{p}}_{\mathit{R}}^{\mathrm{T}}\hfill \end{array}{]}^{\mathrm{T}}$ , ${\mathit{p}}_{\mathit{x}}=\mathit{mv}$ , and ${\mathit{p}}_{\mathit{R}}=\mathit{J}\mathbf{\Omega }$ .

Implicit Euler Lie group methods of DAEs

In vector spaces, the typical DAEs of a multibody system are as follows

$$\{\begin{array}{c}\hfill \stackrel{·}{\mathit{q}}=\mathit{V}\hfill \\ \multicolumn{1}{c}{\mathit{M}\stackrel{·}{\mathit{V}}+{\mathbf{\Phi }}_{\mathit{q}}^{\mathrm{T}}\mathbf{\lambda }=\mathit{Q}(\mathit{q},\mathit{V})}\\ \multicolumn{1}{c}{\mathbf{\Phi }(\mathit{q})=0}\end{array}$$

(8)

Here, $\mathit{q}$ is a usual vector of generalized coordinates. Using the first-order backward difference scheme, ${\stackrel{·}{\mathit{q}}}^{i+1}$ and ${\stackrel{·}{\mathit{V}}}^{i+1}$ are given as

$${\stackrel{·}{\mathit{q}}}^{i+1}=\frac{{\mathit{q}}^{i+1}-{\mathit{q}}^i}{h}$$

(9)

$${\stackrel{·}{\mathit{V}}}^{i+1}=\frac{{\mathit{V}}^{i+1}-{\mathit{V}}^i}{h}$$

(10)

Based on equations (9) and (10), we can get the following implicit Euler method of equations (equation (8)) easily

$$\{\begin{array}{c}\hfill {\stackrel{·}{\mathit{q}}}^{i+1}={\mathit{V}}^{i+1}\hfill \\ \multicolumn{1}{c}{\mathit{M}{\stackrel{·}{\mathit{V}}}^{i+1}+{\mathbf{\Phi }}_{\mathit{q}}^{\mathrm{T}}({\mathit{q}}^{i+1}){\mathbf{\lambda }}^{i+1}=\mathit{Q}({\mathit{q}}^{i+1},{\mathit{V}}^{i+1})}\\ \multicolumn{1}{c}{\mathbf{\Phi }({\mathit{q}}^{i+1})=0}\end{array}$$

(11)

But, for the DAEs on Lie group in the previous section, $\mathit{R}$ is in SO(3) and $\stackrel{·}{\mathit{R}}$ cannot be discretized as equation (10). Based on the CG method, the first-order discretized version of ${\stackrel{·}{\mathit{R}}}_i={\mathit{R}}_i{\tilde{\mathbf{\Omega }}}_i$ is

$${\mathit{R}}^{i+1}={\mathit{R}}^i\exp \left(h{\tilde{\mathbf{\Omega }}}^{i+1}\right)$$

(12)

So the first-order discretization of the pose is

$${\mathit{q}}^{i+1}=\left(\begin{array}{cc}\hfill {\mathit{x}}^i+h{\mathit{v}}^{i+1}\hfill & \hfill {\mathit{R}}^i\exp \left(h{\tilde{\mathbf{\Omega }}}^{i+1}\right)\hfill \end{array}\right)$$

(13)

But the acceleration is a vector, so we can use the same scheme as equation (10). Substituting equations (13) and (10) into equation (1) at instant $t^{i+1}$

$$\{\begin{array}{c}\hfill \mathit{M}{\stackrel{·}{\mathit{V}}}^{i+1}+{\mathit{B}}^{\mathrm{T}}({\mathit{q}}^{i+1}){\mathbf{\lambda }}^{i+1}=\mathit{Q}({\mathit{q}}^{i+1},{\mathit{V}}^{i+1})\hfill \\ \multicolumn{1}{c}{\mathbf{\Phi }({\mathit{q}}^{i+1})=0}\end{array}$$

(14)

One can get the implicit Euler method of DAEs (equation (1)) with index 3 on Lie group as

$$\mathit{F}({\mathit{V}}^{i+1},{\mathbf{\lambda }}^{i+1})=\left[\begin{array}{c}\hfill {\mathit{F}}_1({\mathit{V}}^{i+1},{\mathbf{\lambda }}^{i+1})\hfill \\ \multicolumn{1}{c}{{\mathit{F}}_2({\mathit{V}}^{i+1})}\end{array}\right]=0$$

(15)

where

$$\{\begin{array}{c}\hfill {\mathit{F}}_1({\mathit{V}}^{i+1},{\mathbf{\lambda }}^{i+1})=\mathit{M}({\mathit{V}}^{i+1}-{\mathit{V}}^i)+h{\mathit{B}}^{\mathrm{T}}({\mathit{q}}^{i+1}){\mathbf{\lambda }}^{i+1}-h\mathit{Q}({\mathit{q}}^{i+1},{\mathit{V}}^{i+1})\hfill \\ \multicolumn{1}{c}{{\mathit{F}}_2({\mathit{V}}^{i+1})=\mathbf{\Phi }({\mathit{q}}^{i+1})}\end{array}$$

(16)

The corresponding Newton–Raphson iterative scheme is

$$\begin{array}{c}\left[\begin{array}{c}\hfill {\mathit{V}}^{i+1,k+1}\hfill \\ \multicolumn{1}{c}{{\mathbf{\lambda }}^{i+1,k+1}}\end{array}\right]=\left[\begin{array}{c}\hfill {\mathit{V}}^{i+1,k}\hfill \\ \multicolumn{1}{c}{{\mathbf{\lambda }}^{i+1,k}}\end{array}\right]\hfill \\ -\nabla {\mathit{F}}^{-1}({\mathit{V}}^{i+1,k},{\mathbf{\lambda }}^{i+1,k})\mathit{F}({\mathit{V}}^{i+1,k},{\mathbf{\lambda }}^{i+1,k})\hfill \end{array}$$

(17)

where

$$\nabla \mathit{F}({\mathit{V}}^{i+1},{\mathbf{\lambda }}^{i+1})=\left[\begin{array}{cc}\hfill {\nabla }_{{\mathit{V}}^{i+1}}{\mathit{F}}_1({\mathit{V}}^{i+1},{\mathbf{\lambda }}^{i+1})\hfill & \hfill h{\mathit{B}}^{\mathrm{T}}({\mathit{q}}^{i+1})\hfill \\ \multicolumn{1}{c}{h\mathit{B}({\mathit{q}}^{i+1})}& \hfill 0\hfill \end{array}\right]$$

(18)

${\nabla }_{{\mathit{V}}^{i+1}}{\mathit{F}}_1({\mathit{V}}^{i+1},{\mathbf{\lambda }}^{i+1})$ can be derived by

$$\begin{array}{c}{\delta }_{{\mathit{V}}^{i+1}}{\mathit{F}}_1({\mathit{V}}^{i+1},{\mathbf{\lambda }}^{i+1})={\partial }_{{\mathit{V}}^{i+1}}{\mathit{F}}_1({\mathit{V}}^{i+1},{\mathbf{\lambda }}^{i+1})·\delta {\mathit{V}}^{i+1}\hfill \\ +〈{\partial }_{{\mathit{q}}^{i+1}}{\mathit{F}}_1({\mathit{V}}^{i+1},{\mathbf{\lambda }}^{i+1}),\delta {\mathit{q}}^{i+1}〉\hfill \\ ={\nabla }_{{\mathit{V}}^{i+1}}{\mathit{F}}_1({\mathit{V}}^{i+1},{\mathbf{\lambda }}^{i+1})·\delta {\mathit{V}}^{i+1}\hfill \end{array}$$

(19)

Using equations (13) and (10), the implicit Euler method of DAEs (equation (4)) with index 2 on Lie group is

$$\{\begin{array}{c}{\mathit{R}}^{i+1}={\mathit{R}}^i\exp ({\tilde{\mathbf{\Theta }}}^{i+1})\hfill \\ \left[\begin{array}{c}\hfill {\mathit{x}}^{i+1}-{\mathit{x}}^i\hfill \\ \multicolumn{1}{c}{{\mathbf{\Theta }}^{i+1}}\end{array}\right]=h{\mathit{V}}^{i+1}-{\mathit{B}}^{\mathrm{T}}({\mathit{x}}^{i+1},{\mathit{R}}^{i+1}){\mathbf{\mu }}^{i+1}\hfill \\ \mathit{M}({\mathit{V}}^{i+1}-{\mathit{V}}^i)+h{\mathit{B}}^{\mathrm{T}}({\mathit{q}}^{i+1}){\mathbf{\lambda }}^{i+1}=h\mathit{Q}({\mathit{q}}^{i+1},{\mathit{V}}^{i+1})\hfill \\ \mathbf{\Phi }({\mathit{q}}^{i+1})=0\hfill \\ \mathit{B}({\mathit{q}}^{i+1}){\mathit{V}}^{i+1}=0\hfill \end{array}$$

(20)

The implicit Euler method of DAEs (5) with index 1 on Lie group is

$$\{\begin{array}{c}{\mathit{R}}^{i+1}={\mathit{R}}^i\exp ({\tilde{\mathbf{\Theta }}}^{i+1})\hfill \\ \left[\begin{array}{c}\hfill {\mathit{x}}^{i+1}-{\mathit{x}}^i\hfill \\ \multicolumn{1}{c}{{\mathbf{\Theta }}^{i+1}}\end{array}\right]=h{\mathit{V}}^{i+1}\hfill \\ \mathit{M}({\mathit{V}}^{i+1}-{\mathit{V}}^i)+h{\mathit{B}}^{\mathrm{T}}({\mathit{q}}^{i+1}){\mathbf{\lambda }}^{i+1}=h\mathit{Q}({\mathit{q}}^{i+1},{\mathit{V}}^{i+1})\hfill \\ \mathit{B}({\mathit{q}}^{i+1})({\mathit{V}}^{i+1}-{\mathit{V}}^i)+h\stackrel{·}{\mathit{B}}({\mathit{q}}^{i+1}){\mathit{V}}^{i+1}=0\hfill \end{array}$$

(21)

The implicit Euler method of ODAEs (6) on Lie group is

$$\{\begin{array}{c}{\mathit{R}}^{i+1}={\mathit{R}}^i\exp ({\tilde{\mathbf{\Theta }}}^{i+1})\hfill \\ \left[\begin{array}{c}\hfill {\mathit{x}}^{i+1}-{\mathit{x}}^i\hfill \\ \multicolumn{1}{c}{{\mathbf{\Theta }}^{i+1}}\end{array}\right]=h{\mathit{V}}^{i+1}\hfill \\ \mathit{M}({\mathit{V}}^{i+1}-{\mathit{V}}^i)+h{\mathit{B}}^{\mathrm{T}}({\mathit{q}}^{i+1}){\mathbf{\lambda }}^{i+1}=h\mathit{Q}({\mathit{q}}^{i+1},{\mathit{V}}^{i+1})\hfill \\ \mathbf{\Phi }({\mathit{q}}^{i+1})=0\hfill \\ \mathit{B}({\mathit{q}}^{i+1}){\mathit{V}}^{i+1}=0\hfill \\ \mathit{B}({\mathit{q}}^{i+1})({\mathit{V}}^{i+1}-{\mathit{V}}^i)+h\stackrel{·}{\mathit{B}}({\mathit{q}}^{i+1}){\mathit{V}}^{i+1}=0\hfill \end{array}$$

(22)

The symplectic Euler method of the constrained Hamilton equations on Lie group is

$$\{\begin{array}{c}{\mathit{R}}^{i+1}={\mathit{R}}^i\exp ({\tilde{\mathbf{\Theta }}}^{i+1})\hfill \\ \left[\begin{array}{c}\hfill {\mathit{x}}^{i+1}\hfill \\ \multicolumn{1}{c}{{\mathbf{\Theta }}^{i+1}}\end{array}\right]=\left[\begin{array}{c}\hfill {\mathit{x}}^i\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]+h{\mathit{M}}^{-1}{\overline{\mathit{p}}}^{i+1}\hfill \\ {\overline{\mathit{p}}}^{i+1}-{\mathit{p}}^i+h{\mathit{B}}^{\mathrm{T}}({\mathit{q}}^{i+1}){\mathbf{\lambda }}^{i+1}=h\mathit{Q}({\mathit{q}}^{i+1},{\overline{\mathit{p}}}^{i+1})\hfill \\ \mathbf{\Phi }({\mathit{q}}^{i+1})=0\hfill \end{array}$$

(23)

$$\{\begin{array}{c}{\mathit{p}}^{i+1}-{\overline{\mathit{p}}}^{i+1}+h{\mathit{B}}^{\mathrm{T}}({\mathit{q}}^{i+1}){\mathbf{\mu }}^{i+1}=0\hfill \\ \mathit{B}({\mathit{q}}^{i+1}){\mathit{M}}^{-1}{\mathit{p}}^{i+1}=0\hfill \end{array}$$

(24)

Numerical examples

Figure 1 shows a three-dimensional (3D) single pendulum, for which the generalized coordinate and the corresponding velocity are $\mathit{q}=(\begin{array}{cc}\hfill \mathit{x}\hfill & \hfill \mathit{R}\hfill \end{array})$ and $\mathit{V}=[\begin{array}{cc}\hfill {\mathit{v}}^{\mathrm{T}}\hfill & \hfill {\mathbf{\Omega }}^{\mathrm{T}}\hfill \end{array}{]}^{\mathrm{T}}$ , respectively, the inertia matrices of the system are $\mathit{M}=\mathrm{diag}(\begin{array}{cc}\hfill \mathit{m}\hfill & \hfill \mathit{J}\hfill \end{array})$ , $\mathit{m}=\mathrm{diag}(\overline{m},\overline{m},\overline{m})$ , and $\mathit{J}=\mathrm{diag}(J_1,J_2,J_3)$ , $l=1\mathrm{m},\overline{m}=5\mathrm{kg},r=0.25\mathrm{m},h=0.05\mathrm{m}$ , the generalized force vector is $\mathit{Q}(\mathit{q},\mathit{V})=\left[\begin{array}{c}\hfill -\mathit{mg}\hfill \\ \multicolumn{1}{c}{-\widetilde{\mathbf{\Omega }\mathit{J}\mathbf{\Omega }}}\end{array}\right],\mathit{g}={\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 9.81\hfill \end{array}\right]}^{\mathrm{T}}$ , and the constraint equations are $\mathbf{\Phi }(\mathit{q})=\mathit{x}+\mathit{RX}=0$ , $\mathit{X}=[\begin{array}{ccc}\hfill 0\hfill & \hfill -l\hfill & \hfill 0\hfill \end{array}{]}^{\mathrm{T}}$ . The initial values are $\mathit{R}(0)=\mathit{E},\mathbf{\Omega }(0)=[\begin{array}{ccc}\hfill 0\hfill & \hfill 150\hfill & \hfill -4.6154\hfill \end{array}{]}^{\mathrm{T}}$ .

Figure 1.

A 3D single pendulum.

The methods presented in this article are named as DAEn-IEL, DAEn-IELS, DAEn-SE LS, ODAE-IEL, and ODAE-SEL, where n (n = 1, 2, 3) denotes the index of DAEs, IE denotes implicit Euler method, L denotes the Lie group, SE denotes the symplectic Euler method, and IELS denotes the stabilized integration with Lie group implicit Euler method. Figures 2 and 3 show the results obtained using the DAE3-IEL method with the time step h = 0.0001 s and the simulation time t = 2 s. Figure 2 shows the space trajectory of the end of the 3D single pendulum, the time history of the absolute position of the mass center, and the time history of the Lagrange multiplier. Figure 3 shows the kinetic, potential, and total energies and the constraints of the system. These results indicate that the simulation trajectory has a higher accuracy, and the total energy and constraints are kept well in the smaller time step case.

Figure 2.

The trajectory of the end of the 3D single pendulum by DAE3-IEL.

Figure 3.

The energies and the constraints of the single pendulum by DAE3-IEL.

Figure 4 shows the errors and the orders of the errors of different time steps using DAE3-IEL. $e{x}_1=x_1-{\tilde{x}}_1$ denotes the error of the coordinate x₁, where $x_1$ is the solution of DAE3-IEL and ${\tilde{x}}_1$ is the reference solution obtained using a very small step size h = 0.000001 s. $e{y}_1,e{z}_1,e{\lambda }_1,e{\lambda }_2,\mathrm{and}e{\lambda }_3$ are defined similarly, where ${\lambda }_1,{\lambda }_2,\mathrm{and}{\lambda }_3$ are the Lagrange multipliers and ||.|| is the norm of the errors. It demonstrates that, as the time step increases, the errors are larger. And the orders which are the slopes of log||ex₁||, log||ey₁||, log||ez₁||, log||eλ₁||, log||eλ₂||, and log||eλ₃|| versus log h remain about 1 when the time step h is smaller. But when the time step h is larger, the orders decreased.

Figure 4.

The errors and the orders of the errors of different time steps using DAE3-IEL.

Figure 5 shows the time history of the nine elements of the rotation matrix R during the simulation of the 3D single pendulum by DAE3-IEL and the orthogonality error of the rotation matrix R , that is, $e\mathit{R}={\mathit{R}}^{\mathrm{T}}\mathit{R}-\mathit{I}$ , where $\mathit{I}$ is the identity matrix, which shows that the orthogonality of the rotation matrix can be kept well.

Figure 5.

The rotation matrix and the orthogonality error of the single pendulum by DAE3-IEL.

Figure 6 shows the total energies and the third constraints of the 3D single pendulum obtained by DAE3-IEL, DAE2-IELS, DAE1-IELS, ODAE-IEL, and DAE2-SELS. Results illustrate that the methods DAE3-IEL and ODAE-IEL have the similar accuracy of obtaining the total energy and constraints. The methods DAE2-IELS and DAE2-SELS have the similar accuracy of obtaining the total energy, while DAE2-IELS has the higher accuracy of obtaining the constraints. The errors of the total energy and constraints by DAE1-IELS are maximum among these methods.

Figure 6.

The total energies and the constraints of the single pendulum by different Lie group Euler methods.

Table 1 lists the errors of the total energies and constraints of the 3D single pendulum obtained by DAE3-IEL, DAE2-IELS, DAE1-IELS, ODAE-IEL, and DAE2-SELS. In order to obtain the results of longer time simulation, the terminal time is set to 10 s. ε(H) means the absolute error of the total energy H, δ(H) means the relative error of the total energy H, and ε(Φ) means the absolute error of the displacement constraints Φ. Results show that the maximum errors of the constraints obtained by the methods are all kept in a very small range. The error of the constraints Φ obtained by DAE3-IEL is smaller than that with the other methods, while the error of the total energy H obtained by DAE2-IELS is the smallest.

Table 1.

The errors of the total energy and constraints of the single pendulum by different Lie group Euler methods (t = 10 s).

	Max|ε(H)|	Max|δ(H)|	Max|ε(Φ)|
DAE3-IEL	47.3571	0.0087	6.3283 × 10–15
DAE2-IELS	15.6541	0.0029	5.6162 × 10–12
DAE1-IELS	270.2068	0.0497	3.5303 × 10–11
ODAE-IEL	47.3571	0.0087	7.2997 × 10–15
DAE2-SELS	17.9126	0.0033	1.4410 × 10–11

Figure 7 shows a 3D double pendulum, for which the inertia matrices of the system are ${\mathit{M}}_i=\mathrm{diag}(\begin{array}{cc}\hfill {\mathit{m}}_i\hfill & \hfill {\mathit{J}}_i\hfill \end{array})$ , ${\mathit{m}}_i=\mathrm{diag}({\overline{m}}_i,{\overline{m}}_i,{\overline{m}}_i)$ , and ${\mathit{J}}_i=\mathrm{diag}(J_{i1},J_{i2},J_{i3})$ , $i=1,2$ , $l_1=l_2=1\mathrm{m}$ , ${\overline{m}}_1={\overline{m}}_2=5\mathrm{kg}$ , the generalized force vector is ${\mathit{Q}}_i(\mathit{q},\mathit{V})=[\begin{array}{cc}\hfill -{\mathit{m}}_i\mathit{g}\hfill & \hfill -{\tilde{\mathbf{\Omega }}}_i{\mathit{J}}_i{\mathbf{\Omega }}_i\hfill \end{array}{]}^{\mathrm{T}}$ , $i=1,2$ , $\mathit{g}=[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 9.81\hfill \end{array}{]}^{\mathrm{T}}$ , and the constraint equations are

$$\begin{array}{c}{\mathbf{\Phi }}_1(\mathit{q})={\mathit{x}}_1+{\mathit{R}}_1{\mathit{X}}_{10},{\mathit{X}}_{10}={\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill -l_1/2\hfill \end{array}\right]}^{\mathrm{T}}\hfill \\ {\mathbf{\Phi }}_2(\mathit{q})={\mathit{x}}_1+{\mathit{R}}_1{\mathit{X}}_{11}-{\mathit{x}}_2-{\mathit{R}}_2{\mathit{X}}_{21},\hfill \\ {\mathit{X}}_{11}={\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill l_1/2\hfill \end{array}\right]}^{\mathrm{T}},{\mathit{X}}_{21}={\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill -l_2/2\hfill \end{array}\right]}^{\mathrm{T}}\hfill \end{array}$$

The initial values are

$$\begin{array}{c}{\mathit{R}}_1(0)=\mathit{E},{\mathit{R}}_2(0)=\mathit{E},{\mathbf{\Omega }}_1(0)={\left[\begin{array}{ccc}\hfill 10\pi \hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]}^T,\hfill \\ {\mathbf{\Omega }}_2(0)={\left[\begin{array}{ccc}\hfill 10\pi \hfill & \hfill 10\pi \hfill & \hfill 20\pi \hfill \end{array}\right]}^T\hfill \end{array}$$

Figure 7.

A 3D double pendulum.

Figure 8 shows the time history of the nine elements of the rotation matrix R during the simulation of the 3D double pendulum by DAE3-IEL and the orthogonality error $e\mathit{R}$ of the rotation matrix, which shows that the orthogonality of the rotation matrix can be kept well.

Figure 8.

The rotation matrix and the orthogonality error of the double pendulum by DAE3-IEL.

Figure 9 shows the total energies and the third constraints of the 3D double pendulum obtained by DAE3-IEL, DAE2-IELS, DAE1-IELS, ODAE-IEL, and DAE2-SELS with the time step h = 0.001 s. The results illustrate that the methods DAE3-IEL and ODAE-IEL have the similar accuracy of obtaining the total energy and constraints. The methods DAE2-IELS and DAE2-SELS have the similar accuracy of obtaining the total energy. The errors of the total energy and constraints obtained by DAE1-IELS lie between those of the methods DAE3 and DAE2. In Figure 10, when the time step h = 0.00001 s, the error of the total energy obtained by DAE1-IELS is the smallest.

Figure 9.

The energies and constraints of the double pendulum by different Lie group Euler methods (h = 0.001 s).

Figure 10.

The energies and constraints of the double pendulum by different Lie group Euler methods (h = 0.00001 s).

Table 2 lists the errors of the total energies and the constraints of the 3D double pendulum obtained by DAE3-IEL, DAE2-IELS, DAE1-IELS, ODAE-IEL, and DAE2-SELS. In order to obtain the results of longer time simulation, the terminal time is set to 1 s. The results show that the maximum errors of the constraints obtained by the methods are all kept in a very small range. The error of the constraints Φ obtained by DAE1-IELS is smaller than those with the other methods; meanwhile, the error of the total energy H obtained by DAE1-IELS is the smallest.

Table 2.

The errors of the total energy and constraints of the double pendulum by different Lie group Euler methods (t = 1 s).

	Max|ε(H)|	Max|δ(H)|	Max|ε(Φ)|
DAE3-IEL	179.2942	0.0237	3.3307 × 10–16
DAE2-IELS	96.0631	0.0127	3.7192 × 10–15
DAE1-IELS	93.0519	0.0123	2.6090 × 10–15
ODAE-IEL	179.2942	0.0237	2.8311 × 10–15
DAE2-SELS	96.2233	0.0127	1.7319 × 10–14

Concluding remarks

Five types of DAEs on Lie group of multibody system dynamics are introduced in this article, which are DAEs of index 3, stabilized index 2, index 1, ODAEs, and the relevant Hamilton equations on Lie group. The corresponding five implicit Euler methods on Lie group for DAEs are designed. Numerical examples are used to test and compare the algorithms. The simulation trajectory has a higher accuracy, and the total energy and the constraints are kept well in the smaller time step case. The methods DAE3-IEL and ODAE-IEL have the similar accuracy of obtaining the total energy and constraints. The methods DAE2-IELS and DAE2-SELS also have the similar accuracy. The errors of total energy and constraints obtained by DAE1-IELS are sensitive to the change of time step. Further studies on designing higher order algorithms are needed to improve the accuracy and efficiency for large-scale system simulations.