## Abstract

This work highlights the use of half-implicit numerical integration in the context of the index three differential algebraic equations (DAEs) of multibody dynamics. Although half-implicit numerical integration is well established for ordinary differential equations problems, to the best of our knowledge, no formal discussion covers its use in the context of index three DAEs of multibody dynamics. We wish to address this since when compared to fully implicit methods, half-implicit integration has two attractive features: (i) the solution method does not require the computation of the Jacobian associated with the constraint, friction, contact, or user-defined applied forces; and (ii) the solution is simpler to implement. Moreover, for nonstiff problems, half-implicit numerical integration yields a faster solution. Herein, we outline the numerical method and demonstrate it in conjunction with three mechanisms. We report on convergence order behavior and solution speed. The Python software developed to generate the results reported is available as open in a public repository for reproducibility studies.

## 1 Introduction

*a*as

_{n}This is the half-implicit, symplectic Euler method [6]; it is explicit at the velocity level (level one “L1”), but implicit at the position level (level zero “L0”); and, it is a first-order method, but higher order methods can be obtained by combining a suitable explicit integrator for the L1 update, with a suitably chosen implicit integrator for the L0 update. Higher order integrators fall outside the scope of this technical note. Our contribution is with outlining how integrators like the one in Eq. (1) can be used to directly solve the index three DAEs of multibody dynamics.

Several attractive aspects are noted when using a half-implicit solver to handle the index three DAEs of multibody dynamics. First, when using in conjunction with the Absolute Nodal Coordinate Formulation [7], the Jacobian of the internal forces does not need to be evaluated. Additionally, if no kinematic constraints are present in the system, the mass matrix is constant and only needs to be factored once and used throughout the simulation to compute the analog of *a _{n}* in Eq. (1). Second, the solution does not require the Jacobian of some hard-to-handle forces that might appear in the model, e.g., friction, contact, impact, constraint, user-defined. The latter are particularly difficult to handle, since in many cases user-defined forces do not have closed-form expression and are provided as a black box, in a user subroutine. Finally, the solution process for half-implicit integrators is straightforward. The Jacobian is easy to assemble and has a sparsity pattern that can be leveraged for linear-time solution [8,9].

This Technical Note is organized as follows. Section 2.1 states the index three DAEs problem of interest; Sections 2.2 and 2.3 outline the half- and fully-implicit schemes, respectively; Section 3 reports convergence and timing results. We close with a discussion of the main highlights and directions of future work.

## 2 Solution Approach

Herein, we assume the use of absolute coordinates. The location, ** r**, of each body is with respect to a global reference frame. The orientation, or pose, of each rigid body is also with respect to the global reference frame, captured by Euler angles $\epsilon \u2261[\varphi $ −

*θ*− $\psi ]T$ or Euler parameters (quaternions) $p\u2261[e0,e1,e2,e3]T$. However, herein we use Lie-group integration to produce the orientation matrix $A$ directly, since this approach is more attractive than using Euler parameters or Euler angles [10]. The discussion herein focuses on rigid body dynamics but the conclusions carry over to absolute nodal coordinate formulation (ANCF)-based flex body dynamics. Finally, the local reference frame (LRF) associated with each body is assumed to be central and principal, so that the generalized mass matrix is constant and diagonal.

### 2.1 The Constrained Equations of Motion.

The state of a system of *nb* rigid bodies, each of mass *m _{i}* and inertia matrix $J\xafi$, can be described by translational and angular acceleration, $r\xa8i$ and $\omega \xaf\u02d9i$; velocity, $r\u02d9i$ and $\omega \xafi$; and position $ri$ and orientation matrix $Ai$, where $i=1,2,\u2026,nb$. Herein, a letter with a bar $\xb7\xaf$ represents a geometric quantity expressed in the body's LRF.

*t*is time. By applying D'Alembert's principle, one can derive the constrained equations of motion [11] in matrix form as

where $M\u2261diag{m1I3,\u2026,mnbI3},\u2009J\xaf\u2261diag{J\xaf1,\u2026,J\xafnb},\u2009f\u2261[f1T,\u2026,fnbT]T$, and $\tau \xaf\u2261[\tau \xaf1T,\u2026,\tau \xafnbT]T$. Here, $Mr\xa8$ and $J\xaf\omega \xaf\u02d9$ are the inertia force and torque, respectively; **f*** _{i}* is the external force on body

*i*; the external torque $n\xafi$ appears in $\tau \xafi\u2261n\xafi\u2212\omega \xaf\u0303iJ\xafi\omega \xafi$; $\Phi ,rT\lambda $ and $\Pi \xafT(\Phi )\lambda $ are the reaction forces and torques from the joints, with Lagrange multipliers $\lambda \u2208\mathbb{R}nc$ associated with the kinematic constraints $\Phi $, and $\Phi ,r$ and $\Pi \xaf(\Phi )$ being the first-order variations of $\Phi $ with respect to

**and**

*r***. Note that in $\Phi $, the orientation matrix**

*A***only shows up in the form of $As\xaf$, with $s\xaf$ being a vector in the local body frame. As shown in Ref. [10], $\Pi \xaf(\Phi )$ has cleaner formulations compared to taking the partial derivatives of $\Phi $ with respect to Euler angles or Euler parameters.**

*A*### 2.2 The Half-Implicit Approach.

*t*=

*t*, the discretized equations of motion are

_{n}*t*

_{n}*k*), and solve for the corrections $\delta q$

Then, $N\delta q=e(k)$ is solved iteratively until the correction $\delta q$ meets the stopping criteria $||\delta q||<\u03f5tol$. Once the solver converges, the solution $\lambda n$ is plugged back in Eq. (3) to compute the accelerations $r\xa8n$ and $\omega \xaf\u02d9n$. Note that the Jacobians $\mathbb{P}n,r$ and $\mathbb{P}n,A$ in Eq. (12) are derived for all lower-pair joints in Ref. [10]. Thus, the solution approach can be effectively used in conjunction with any constrained multibody system.

### 2.3 The Fully-Implicit Approach.

*c*) is for improving the condition number of the Newton iteration matrix. The velocity- and position-level quantities are updated using the Backward Euler scheme (compare Eqs. (4

*a*) and (4

*b*) to (14

*a*) and (14

*b*)),

When assembling the Newton iteration matrix for solving Eq. (13), the Jacobians of the reaction forces and torques, $\mathbb{P}n+1,rT\lambda n+1$ and $\mathbb{P}n+1,AT\lambda n+1$, with respect to ** r**, $r\u02d9$,

**and $\omega \xaf$ are dropped to improve the solver's efficiency.**

*A*## 3 Numerical Experiments

The software implemented to generate the results reported herein is available in a public repository for reproducibility studies and further development [12]. Three systems: a double pendulum, a slider–crank, and a four–link mechanism are investigated to compare the performance of the fully- and half-implicit approaches.

The double pendulum is composed of 4 m and 2 m long rectangular bars, connected via revolute joints. The system is released from the configuration shown in Fig. 1(a) and swings owing to the gravitational pull. Both bars have the same square-shaped cross section and density of 7.8 g cm^{−3}.

The slider–crank and four–link systems are kinematically driven. The mass and inertia properties of the slider–crank mechanism are in Table 1. We prescribe the crank (body 1 in Fig. 1(b)) to rotate about its negative *x*-axis at a rate of $2\pi $ rad s^{− 1}. The four–link setup consists of a rotor disk and two bars of the same square-shaped cross section. The density of each body is 1.5 g cm^{−3}. The rotor (body 1 in Fig. 1(c)) is prescribed to rotate about its *z*-axis at a rate of *π* rad s^{− 1}. More parameters are given in Table 2. All simulations were run for 8 s, which allowed for at least three full rotations for the slider–crank and four–link mechanisms.

Body | Mass | $Ix\xafx\xaf$ | $Iy\xafy\xaf$ | $Iz\xafz\xaf$ |
---|---|---|---|---|

1 | 0.12 | $1\xd710\u22124$ | $1\xd710\u22125$ | $1\xd710\u22124$ |

2 | 0.5 | $4\xd710\u22123$ | $4\xd710\u22124$ | $4\xd710\u22123$ |

3 | 2 | $1\xd710\u22124$ | $1\xd710\u22124$ | $1\xd710\u22124$ |

Body | Mass | $Ix\xafx\xaf$ | $Iy\xafy\xaf$ | $Iz\xafz\xaf$ |
---|---|---|---|---|

1 | 0.12 | $1\xd710\u22124$ | $1\xd710\u22125$ | $1\xd710\u22124$ |

2 | 0.5 | $4\xd710\u22123$ | $4\xd710\u22124$ | $4\xd710\u22123$ |

3 | 2 | $1\xd710\u22124$ | $1\xd710\u22124$ | $1\xd710\u22124$ |

### 3.1 Convergence Order Analysis.

The convergence order analysis starts with producing a reference solution. The position, velocity, and acceleration of the double pendulum were obtained by integrating a set of ordinary differential equations that describe the system in 2D using two coordinates. To that end, we use a highly accurate, 8th-order Runge–Kutta solver – DoPri853 [13]—and a constant step size of $h=1\xd710\u22126\u2009s$. For the slider–crank and four–link mechanisms, the reference solution was generated by posing the problem as a kinematic one. To this end, only Eq. (2*c*) and its time derivatives were solved using the Newton–Raphson method, which allowed for the reference solution to be obtained within machine precision.

Figure 2 shows on a log–log scale the difference between the fully- and half-implicit integrators and the reference solution plotted over a range of step sizes. Labels “HI” and “FI” stand for “Half Implicit” and “Fully Implicit,” respectively. The tolerance of the fully-implicit approach is set to $\u03f5tol=10\u221210/h2$; for the half-implicit, $10\u221210$. The scaling factor $1/h2$ is from Eq. (13*c*), in comparison to Eq. (2*c*). The error relative to the reference solution was calculated as $\u2211i=1n(ziDAE\u2212ziref)2/n$, where *n* is the total number of time steps in each simulation and $ziDAE$ and $ziref$ are quantities at time-step *t _{i}* produced by the DAE solvers and the reference solution, respectively.

Figure 2(a) plots the *z*-position error for bar 2. The most important feature noted is a slope of 1.0 for both integrators. In other words, if the step size is reduced by a factor of 10, the error will decrease by a factor of 10. Figures 2(b) and 2(c) plot the velocity order analysis for the slider of the slider–crank mechanism and link body 2 of the four–link mechanism, respectively. All three figures show first-order convergence for both methods, thus meeting expectations. Finally, since the half-implicit integrator is symplectic, it conserves a numerical Hamiltonian, which leads to superior performance in terms of long-term energy preservation [6]. This is evident in Fig. 3—the total energy of the system drains over time using a fully-implicit approach, a highly-damped Backward Euler integrator, but holds steady for the symplectic half-implicit approach.

### 3.2 Performance Analysis.

For a fair timing comparison, each integrator had its own *ϵ*_{tol} to ensure that they achieved similar quality of the solution (relative to the reference solution). More plots of how *ϵ*_{tol} affects the solution quality are available in Ref. [14]. Note that both solvers use SciPy's built-in sparse solver. Each timing test was run five times and the average CPU time was recorded for three different step sizes, $h=1\xd710\u22124,1\xd710\u22123$, and $1\xd710\u22122\u2009s$, as shown in Figs. 4(a), 5(a), and 6(a). The average number of iterations for the Newton solver to converge is illustrated in Figs. 4(b), 5(b), and 6(b). A larger step size results in more iterations because the initial starting point for the Newton method is less accurate. The half-implicit approach is faster for all mechanisms and step sizes. This is a direct consequence of the fact that the Newton iteration matrix is easier to evaluate and more accurate. Indeed, the Newton iteration matrix of the fully-implicit approach requires the computation of the partial derivative of the constraint reaction forces and torques, yet in practice these partial derivatives are ignored since they are exceedingly expensive to evaluate. However, for the half-implicit approach, the constraint reaction forces and torques are expressed at time-step *t _{n}*, independent of the position or orientation at $tn+1$, resulting in a simpler and better approximation of the Newton iteration matrix.

### 3.3 Scaling Analysis.

A scaling analysis was carried out to confirm that the speed-up observed in Sec. 3.2 translates to larger problems. As shown in Fig. 7, *N* rigid links were connected by spherical joints and released under gravity ** g** pointing in a random direction (not vertically down). At time

*t*=

*0, the bars were assigned an angular velocity of random magnitude*

*ω*around their longitudinal axis. Tests were conducted for $N=1,\u20092,\u20094,\u20096,\u20098,\u200916,\u200932$. The timing results reported in Fig. 8 indicate that the half-implicit approach is consistently faster than the fully-implicit one.

_{i}### 3.4 An Example With Joint Friction.

To demonstrate the versatility of the half-implicit scheme, the slider-crank model is enhanced to include joint friction; thus, the external force ** f** in Eq. (2

*a*) changes over time. Herein, Coulomb friction is applied in the translational joint between the slider and the ground. The direction of the friction force opposes the relative motion of the slider with respect to the ground; its magnitude is the product of the friction coefficient,

*μ*, and the reaction force in the normal-to-the-motion direction (herein, the

*z*direction). No distinction is made between static and kinetic friction. Keeping with the half-implicit nature of the discretization scheme, the friction force at the current time-step is evaluated explicitly based on the joint normal force from the previous step. Three simulations were carried out using a prescribed motion of the crank, but various friction coefficients –

*μ*= 0, 0.2, and 0.4. In all cases the step size was $1\xd710\u22123$ s. Figure 9 shows the torques required to enforce the prescribed crank motion, in each of the three cases. With increasing

*μ*(and thus joint friction force), the driving torque used to drive the crank experiences larger swings. The sudden jump in the driving torque value results from the change of the slider velocity direction (and that of the friction force).

## 4 Discussion

The half-implicit integrator for index three DAEs advances the simulation from *t _{n}* to $tn+1$ by computing at

*t*a set of Lagrange multipliers $\lambda n$ that enforce the kinematic constraints equations at $tn+1$; i.e., $\Phi (rn+1,An+1,tn+1)=0$. With $\lambda n$ in hand, the acceleration at

_{n}*t*can be obtained by solving a linear system of equations whose matrix is diagonal and constant throughout the simulation. Once the acceleration at

_{n}*t*is known, one uses an explicit integrator to compute the velocities at $tn+1$ (see Eqs. (4

_{n}*a*) and (4

*b*)). The new position and orientation can be computed using an implicit formula with the new velocities (see Eqs. (4

*c*) and (4

*d*)).

The method works for arbitrary multibody systems: with open and closed loops, arbitrary number of bodies, any collection of joints (lower-pair constraints) or bilateral kinematic constraints. Although not discussed herein, the presence of nonholonomic constraints poses no additional challenges, as derived in Sec. 3.4 in Ref. [14]. In terms of the convergence of the Newton solver, the half-implicit method performs slightly better than the fully-implicit one. This can be traced back to the fact that some entries in the Jacobian matrix of the fully-implicit integrator were dropped—they are both expensive to compute and challenging to produce, e.g., Jacobian of constraint forces, friction forces and contact forces. Should one use the exact Jacobian, convergence can improve; however, the fully-implicit solver will be computationally more expensive. Note that since the half-implicit formulation uses reaction forces and torques from the previous time-step, $\mathbb{P}n,rT\lambda n$ and $\mathbb{P}n,AT\lambda n$, it does not require the Jacobian of external or internal generalized forces. Consequently, a simpler yet more accurate Newton matrix is obtained, leading to faster convergence.

Finally, the half-implicit approach is expected to work equally well for the index three DAEs problem formulated using Euler angles or Euler parameters as generalized coordinates. However, as shown in Ref. [10], the fully-implicit ** rA** approach is faster than these alternatives, which justified the choice for Lie-integration in this technical note.

## 5 Conclusions and Future Work

We discuss the use of a half-implicit integration method for the solution of the index three DAEs of multibody dynamics. The approach stands to benefit the simulation of multibody systems that are not exceedingly stiff. In the case of stiff problems, a fully implicit DAE solution is expected to be superior. We demonstrate the half-implicit approach in conjunction with a first-order Lie-group integrator, yet the method can be equally applied to solutions that rely on Euler parameters or Euler angles. The half-implicit scheme treats the reaction forces and torque explicitly in time, leading to a simpler and better approximation of the Jacobian matrix in the Newton step and a more expeditious solution. The scaling analysis demonstrated that the conclusion above carries to larger systems. Finally, the numerical solution is straightforward to implement, see Ref. [12].

There are two things that remain to be investigated. First, provide a formal proof that the half-implicit approach has a guaranteed convergence order when used in the context of index three DAEs of multibody dynamics. Second, we are currently working on a second-order half-implicit integration formula that relies on Lie-group integration in multibody dynamics. Adaptive step size should be considered to further improve the efficiency of the numerical solution.

## Funding Data

Army Research Office (Award No. W911NF1910431; Funder ID: 10.13039/100000183).

Directorate for Computer and Information Science and Engineering (Award No. CISE1835674; Funder ID: 10.13039/100000083).

Division of Civil, Mechanical and Manufacturing Innovation (Award ID: CMMI2153855; Funder ID: 10.13039/100000147).

Office of Cyber infrastructure (Award No. OAC2209791; Funder ID: 10.13039/100000105).