# Efficient application of Gauss' principle to generic mechanical systems

A Fumagalli, P Masarati

► To cite this version:

A Fumagalli, P Masarati. Efficient application of Gauss' principle to generic mechanical systems. Journal of Multi-body Dynamics, 2009, 223 (2), pp.121 - 131. 10.1243/14644193JMBD182. hal-01375988

HAL Id: hal-01375988

https://hal.science/hal-01375988

Submitted on 5 Oct 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Efficient application of Gauss' principle to generic mechanical systems

A Fumagalli and P Masarati*

Politecnico di Milano, Dipartimento di Ingegneria Aerospaziale, Milano, Italy

**Abstract:** This article discusses one of the solutions proposed in the literature to Gauss' principle known as principle of 'least constraint', proposing a clarifying interpretation that allows significant computational improvements when considering its application to problems with non-diagonal mass matrix. The case of non-symmetric mass matrix is briefly discussed as well.

**Keywords:** analytical dynamics, Gauss' principle, least constraint principle, Cholesky decomposition, non-symmetric inertia matrix

## 1 INTRODUCTION

In computational mechanics, the importance of constrained systems modelling is highlighted by its use in many fields of application. Although the problem of formulating the laws that systems of bodies must obey during their motion has been solved many years ago by great scientists and mathematicians like Newton, d'Alembert, Lagrange, and Euler, the very same principles can be viewed from different perspectives and formulated in different, yet equivalent, manners, as pointed out by Gauss himself in reference [1]. This is the case, e.g. of Hamilton's, Gauss', Jourdain's, and Gibbs-Appell's principles, and of Maggi-Kane's equations, as illustrated in many textbooks and articles on analytical dynamics [2, 3].

One of the reasons those principles have been formulated for is the capability to describe the motion of constrained systems without the need to explicitly take into account the details of the kinematics of the constraints.

This article discusses the principle of 'least constraint', as formulated by Gauss [1], and its solution proposed in recent years by Udwadia and Kalaba [4, 5]. The latter formulation is presented in the literature as an elegant means of simulating mechanical systems with redundant or varying constraints [6]. In reference [7], Arabyan and Wu state that the greatest merit

of the Udwadia-Kalaba formulation is that it 'is more advantageous for systematic generation of equations of motion because the user of a program based on this formulation does not have to worry about the possible redundancy of a pair of constraints'.

Although redundant constraints do not need special care, the solution proposed in references [4] and [5] needs some care in the choice of the coordinates and in formulating the inertia properties of the model. In fact, the Udwadia-Kalaba formulation makes intense use of the square root of the mass matrix. Its computation is straightforward when the matrix is diagonal. However, when dealing with non-diagonal matrices, the computation of the square root of the matrix can be very expensive, possibly reducing the appeal of the formulation. This may happen, e.g. when modelling rigid and/or flexible multibody systems, using a consistent inertia formulation, or accounting for the rigid-body motion of bodies whose inertia tensor is not isotropic, or whose motion is not referred to the centre of mass.

This work illustrates a modification to the original solution that reformulates it in a computationally efficient manner in those cases where a non-diagonal, possibly non-symmetric mass matrix must be dealt with, while maintaining all the qualities of the original one.

In section 2, Gauss' principle of least constraint is presented, while section 3 illustrates the expression of the corresponding accelerations as given by Udwadia and Kalaba in references [4] and [5]. A computationally more efficient formulation when dealing with non-diagonal mass matrices is presented in section 4,

<footnote>* Corresponding author: Dipartimento di Ingegneria Aerospaziale, Politecnico di Milano, via La Masa 34, Milano 20156, Italy.
email: masarati@aero.polimi.it</footnote>
which helps giving a deeper insight into the meaning of Udwadia-Kalaba's solution. Section 5 illustrates how the solution of section 3 can be obtained by a constrained minimization of the Gaussian functional, augmented by the constraint equations in a Lagrangian multipliers' sense. Appendix 2 shows how the principle and the proposed solution can be applied to those cases characterized by a non-symmetric inertia matrix, illustrated by practical examples.

## 2 GAUSS' PRINCIPLE

Gauss (Carl Friedrich Gauss (1777–1855), the great German mathematician), in his 1829 paper [1], formulated the principle of least constraint in mechanics, a generalization of an intuition of Maupertuis (Pierre-Louis Moreau de Maupertuis (1698–1759) was a French philosopher and mathematician who, in a dissertation presented in 1741 to the Academy of Sciences of Paris, ‘Loi du repos des corps’ (Law of bodies at rest), formulated an essentially philosophical principle about a quantity, generically called ‘action’, that tends to be minimized in all natural phenomena, thus the name of ‘least action’).

The principle states that the functional

$$G(\ddot{x}) = \frac{1}{2}(\ddot{x} - a)^{\mathrm{T}}\mathbf{M}(\ddot{x} - a) \qquad (1)$$

is minimal, where

$$a = \mathbf{M}^{-1}f \qquad (2)$$

are the accelerations caused by the external forces f applied to the corresponding unconstrained system of inertia matrix M. The functional $G(\ddot{x})$ is called the Gaussian function.

A strict minimum of $G(\ddot{x})$ must exist, since matrix M is positive-definite (matrix M is positive definite iff $v^T M v > 0 \forall v \neq 0$; if the matrix is semi-definite, namely $v^T M v \geq 0 \forall v \neq 0$, then the order of the system can be reduced by eliminating the degrees of freedom with no associated inertia, resulting in a positive-definite reduced matrix). In the case of an unconstrained system, the accelerations $\ddot{x}$ are coincident with those of equation (2), and the functional is exactly zero.

If the system is constrained by ideal bilateral constraints in the form

$$A(x, \dot{x}, t)\ddot{x} = b(x, \dot{x}, t) \qquad (3)$$

the minimization of $G(\ddot{x})$ can be used to find the right value of $\ddot{x}$ among those that comply with equation (3). This constraint equation can either be the second derivative (in this case, matrix A would possibly

depend on $x$ and $t$, but not on $\dot{x}$) of an equation

$$\phi(x, t) = 0 \qquad (4a)$$

$$(\phi_{/x})\dot{x} = -\phi_{/t} \qquad (4b)$$

$$\underbrace{\phi_{/x}}_{A} \ddot{x} = \underbrace{-(\phi_{/x}\dot{x} + \phi_{/t})\dot{x} - (\phi_{/x}\dot{x} + \phi_{/t})}_{b} \qquad (4c)$$

describing a set of holonomic constraints, or the first derivative (In this case, matrix A could also depend on $\dot{x}$ if the constraint expressed by $\psi$ is non-linear in $\dot{x}$; this extension is attributed to Udwadia and Kalaba [8]) of an equation

$$\psi(x, \dot{x}, t) = 0 \qquad (5a)$$

$$\underbrace{\psi_{/\dot{x}}}_{A} \ddot{x} = \underbrace{-\psi_{/x}\dot{x} - \psi_{/t}}_{b} \qquad (5b)$$

describing a set of non-holonomic constraints.

The principle is equivalent to d'Alembert-Lagrange equations, but presents them in a different, interesting variational light. A formal verification is given in section 5.

## 3 UDWADIA-KALABA'S FORMULATION

According to the formulation proposed by Udwadia and Kalaba [4, 5], the only acceleration $\ddot{x}$ that simultaneously complies with the constraints of equation (3), and minimizes Gauss' functional expressed by equation (1), is

$$\ddot{x} = a + \mathbf{M}^{-1/2}(\mathbf{A}\mathbf{M}^{-1/2})^{+}(b - Aa) \qquad (6)$$

where $(\cdot)^{+}$ indicates the Moore-Penrose pseudoinverse of a matrix [9, 10], whose properties are well described in reference [11].

This formulation requires the computation of the square root of the inertia matrix M, and of its inverse, an expensive operation when the matrix is not diagonal, although symmetric and positive definite.

The forces applied by the constraints to the otherwise unconstrained system are

$$f_c = \mathbf{M}^{1/2}(\mathbf{A}\mathbf{M}^{-1/2})^{+}(b - Aa) \qquad (7)$$

The solution proposed by Udwadia and Kalaba appears interesting because its use of a pseudoinversion allows to transparently handle overconstrained systems, which are characterized by a rank-deficient matrix A.

This formulation is appealing because it allows to compute the accelerations directly from the external loads as if the problem were not constrained but, unlike other approaches, it does not need to preliminarily reduce the coordinate set. Since the resulting equations are formally independent and differential,
numerical integration techniques for ordinary differential equations (ODE) can be used. In this case, however, the integration may suffer from drift, since only the second derivative of the constraints is enforced, so only the accelerations are guaranteed to be compliant; constraint stabilization is thus required. This can be provided by a variety of techniques, ranging from the early work of Baumgarte [12], extended by other authors (e.g. reference [13]), to augmented Lagrangian [14] and projection methods [3]. An interesting review is presented in reference [15].

## 3.1 The square root of a matrix

Equations (6) and (7) require the square root of the mass matrix. This section briefly discusses the implications of computing the square root of a matrix, in view of eliminating this requirement.

The square root of a matrix **H** is another matrix **H**¹/² that, multiplied by itself, yields the original one, namely

$$ \mathbf{H}^{1/2} \mathbf{H}^{1/2} = \mathbf{H} \qquad (8) $$

The requirements for the square root of a matrix to exist and be real are usually met by inertia matrices, typically symmetric positive-definite.

Several techniques have been developed to compute the square root of a matrix. For example, the one based on matrix diagonalization, and other, more efficient algorithms, like those based on the matrix sign function method [16–18], or on Schur decomposition [19, 20].

The spectral decomposition of the original matrix, when it is diagonalizable (a mass matrix, usually symmetric positive-definite, is) yields

$$ \mathbf{H} = \mathbf{V} \mathbf{D} \mathbf{V}^{-1} \qquad (9) $$

where **V** is a complete basis for **H** and **D** is a diagonal matrix containing the eigenvalues of matrix **H**. This formulation can be extended to generic non-singular (it is assumed in this work that problems with singular inertia matrices can be reduced in order to eliminate the singularity) matrices; in that case matrix **D** will be in Jordan canonical form instead of diagonal [20].

The square root of a diagonal matrix consists of a diagonal matrix whose diagonal coefficients are the square roots of the coefficients of the original matrix. If **H** is positive definite, or at least positive semidefinite, its eigenvalues are real non-negative; in this case, the coefficients of the square root matrix are real as well. So matrix **H**¹/² is

$$ \mathbf{H}^{1/2} = \mathbf{V} \mathbf{D}^{1/2} \mathbf{V}^{-1} \qquad (10) $$

where the positive roots of matrix **D** are arbitrarily chosen for **D**¹/². The computational complexity of this

algorithm goes from 12n³ for a symmetric matrix to 21n³ in the general case.

Similarly, since the inverse of matrix **H** is

$$ \mathbf{H}^{-1} = \mathbf{V} \mathbf{D}^{-1} \mathbf{V}^{-1} \qquad (11) $$

where **D**⁻¹ is a diagonal matrix whose diagonal contains the inverses of the diagonal coefficients of matrix **D**, the inverse of matrix **H**¹/² is

$$ \mathbf{H}^{-1/2} = \mathbf{V} \mathbf{D}^{-1/2} \mathbf{V}^{-1} \qquad (12) $$

The Denman and Beavers iterative algorithm [16] consists of recursively computing matrices

$$ \mathbf{Y}_{k+1} = \frac{1}{2}(\mathbf{Y}_k + \mathbf{Z}_k^{-1}) \qquad (13a) $$

$$ \mathbf{Z}_{k+1} = \frac{1}{2}(\mathbf{Z}_k + \mathbf{Y}_k^{-1}) \qquad (13b) $$

initialized as $\mathbf{Y}_0 = \mathbf{M}$ and $\mathbf{Z}_0 = \mathbf{I}$. Matrices $\mathbf{Y}$ and $\mathbf{Z}$ converge quadratically to $\mathbf{M}^{1/2}$ and $\mathbf{M}^{-1/2}$, respectively. The cost of this algorithm is that of several matrix factorizations, repeated as many times as needed for quadratic convergence.

Another method [19] consists of computing the Schur decomposition of matrix **M** as

$$ \mathbf{M} = \mathbf{Q} \mathbf{S} \mathbf{Q}^{\mathrm{T}} \qquad (14) $$

where **Q** is a unitary matrix and **S** is upper triangular (quasi-upper triangular when matrix **M** is not diagonalizable). Then, a matrix **T**, such that

$$ \mathbf{T}^2 = \mathbf{S} \qquad (15) $$

can be easily and efficiently computed based only on the square root of the diagonal coefficients (the matrix square root of the diagonal blocks when matrix **M** is not diagonalizable) of **S**. Thus the square root of matrix **M** is given by

$$ \mathbf{M}^{1/2} = \mathbf{Q} \mathbf{T} \mathbf{Q}^{\mathrm{T}} \qquad (16) $$

This algorithm has been extended in reference [20] to make use only of real arithmetic when applied to real matrices that have real square roots, also considering numerical stability issues [21]. The computational complexity of the real arithmetic version of this algorithm, according to reference [20], is 15n³ for the Schur decomposition, n³/6 for the computation of **T**, and 3n³/2 for equation (16). The overall complexity is thus 50n³/3.

In conclusion, the formulation proposed by Udwadia and Kalaba could require the solution of a full eigenproblem, or in any case a similarly expensive matrix decomposition, with little chances to be able to exploit any mass matrix sparsity, unless matrix **M** is already diagonal. As a consequence, the direct application of equation (6) could be quite expensive.
## 4 A BETTER COMPUTATIONAL FORM

A slightly different formulation is proposed, which is shown to possess the same properties of that presented by Udwadia and Kalaba, but requires much less computational effort when matrix M is not diagonal.

### 4.1 What property of the matrix square root is actually required?

In Udwadia-Kalaba's formulation, the square root of matrix M is essentially used because of its property expressed by equation (8). Since matrix M is symmetric positive-definite, if the rows of matrix V are normalized with unit norm, then equations (9) and (10) can be, respectively, written as

$$H = VDV^T \qquad (17)$$

and

$$H^{1/2} = VD^{1/2}V^T \qquad (18)$$

which shows that both are symmetric by construction. As a consequence, in this case equation (8) can also be written as

$$(H^{1/2})^T H^{1/2} = H \qquad (19)$$

In fact, what is basically required by the weighted least-squares interpretation of Gauss' principle is to be able to describe equation (1) in the form of the dot product of a vector by itself, namely

$$ 
\begin{aligned}
G(\ddot{x}) &= \frac{1}{2}(\ddot{x} - a)^T M(\ddot{x} - a) \\
&= \frac{1}{2}[M^{1/2}(\ddot{x} - a)]^T [M^{1/2}(\ddot{x} - a)]
\end{aligned}
\qquad (20) $$

This suggests that any other convenient decomposition that allows to write matrix M as the product of a matrix by its transpose could be used as well.

A very convenient decomposition is the one proposed by Cholesky [22], namely

$$M = U^T U \qquad (21)$$

where U is an upper-triangular matrix (this decomposition is occasionally found in the literature in the alternative form $M = LL^T$, where L = U^T is a lower-triangular matrix). This transformation was already suggested in reference [23] as a means to turn a constrained dynamics problem into a minimization problem for the Euclidean norm of a variable w = Ux, corresponding to a non-Euclidean norm for the variable x.

The decomposition of equation (21) is perfectly suitable to matrix M because the matrix is positive-definite

by definition, and usually symmetric, yielding

$$G(\ddot{x}) = \frac{1}{2} [U(\ddot{x} - a)]^T [U(\ddot{x} - a)] \qquad (22)$$

Moreover, this decomposition can also be used in equation (2) to efficiently compute the accelerations of the unconstrained system required in equations (1) and (6), and thus only needs to be computed once for two purposes, allowing one to 'kill two birds with a stone', a highly appreciated art when developing computationally intensive scientific software. The computational complexity of this decomposition is $n^3/3$, thus more than 30 times less than any of the algorithms illustrated in section 3.1 to extract the matrix square root. Moreover, it can be dramatically reduced when matrix sparsity is exploited.

### 4.2 Is the proposed modification to Udwadia-Kalaba's formula fit?

In order to demonstrate that the proposed modification to Udwadia-Kalaba's formulation is correct, one needs to show that the accelerations resulting from the formula

$$\ddot{x} = a + U^{-1}(AU^{-1})^{+}(b - Aa) \qquad (23)$$

simultaneously:

(a) comply with the constraints of equation (3);

(b) minimize the functional of equation (1).

Compliance with requirement (a) is ensured by the fact that, after substituting equation (23) into equation (3), one obtains

$$ 
\begin{aligned}
b &= Aa + AU^{-1}(AU^{-1})^{+}(b - Aa) \\
&= Aa + (AU^{-1})(AU^{-1})^{+}b - \underbrace{(AU^{-1})(AU^{-1})^{+}(AU^{-1})}_{AU^{-1}}(Ua) \\
&= (AU^{-1})(AU^{-1})^{+}b
\end{aligned}
\qquad (24) $$

The requirement is strictly satisfied only when $AU^{-1}$ is full row-rank, implying that A is full row-rank, and thus $(AU^{-1})(AU^{-1})^{+} = I$. However, if equation (24) is not satisfied not only are the constraints redundant, but also the term b is inconsistent. The use of Gauss' principle of least constraint to address inconsistently constrained systems is illustrated in reference [24]. Proof: if matrix A is not full rank, equation (3) can only be satisfied when its right-hand side, vector b, belongs to the rangespace of matrix A. Consider a singular value decomposition (SVD) of matrix A

$$ 
\begin{aligned}
A = [U_{m \times r} | U_{m \times (m-r)}] \begin{bmatrix} \Sigma_{r \times r} & \mathbf{0}_{r \times (n-r)} \\ \mathbf{0}_{(m-r) \times r} & \mathbf{0}_{(m-r) \times (n-r)} \end{bmatrix} \\
\times [V_{n \times r} | V_{n \times (n-r)}]^T
\end{aligned}
\qquad (25) $$

where submatrices $U_{(i) \times (i)}$, $V_{(i) \times (i)}$ are unitary, while submatrix $\Sigma_{r \times r}$ is diagonal and contains the non-null
singular values of **A**. Vector **b** must belong to the space represented by submatrix **U**<sub>m×r</sub>

$$
\mathbf{b} = \mathbf{U}_{m \times r} \mathbf{b}^* \qquad (26)
$$

$\mathbf{b}^* \in \mathbb{R}^r$. The problem can be reformulated by left-
multiplying equation (3) by the transpose of matrix
$\mathbf{U}_{m \times r}$, yielding

$$
\mathbf{A}^* \ddot{\mathbf{x}} = \mathbf{b}^* \quad (27)
$$

with

$$
\begin{align}
\mathbf{A}^* &= \mathbf{U}_{m \times r}^{\mathrm{T}} \mathbf{A} \nonumber \\
&= \boldsymbol{\Sigma}_{r \times r} \mathbf{V}_{n \times r}^{\mathrm{T}} \tag{28}
\end{align}
$$

which is full rank. The problem of equation (27) is for-
mally equivalent to the original one when matrix **A** is
full rank ■.

Compliance with requirement (b) is achieved by
showing that any other acceleration that differs from
the one of equation (23), while complying with
equation (3), implies a larger value of Gauss' func-
tional, equation (1). Proof: consider an acceleration

$$
\ddot{u} = \ddot{x} + \ddot{y} \tag{29}
$$

that complies with the constraints of equation (3),
namely

$$
\begin{align}
\mathbf{b} &= \mathbf{A}\ddot{\mathbf{u}} \nonumber \\
&= \mathbf{A}(\ddot{\mathbf{x}} + \ddot{\mathbf{y}}) \tag{30}
\end{align}
$$

Since $\ddot{x}$ already complies with equation (3), this implies
that $\ddot{y}$ belongs to the nullspace of $\mathbf{A}$; as a consequence

$$
\mathbf{A}\ddot{\mathbf{y}} = \mathbf{0} \tag{31}
$$

This in turn implies that

$$
\mathbf{A}\mathbf{U}^{-1}\mathbf{U}\ddot{\mathbf{y}} = \mathbf{0} \qquad (32)
$$

The pseudoinverse of a null vector is the transpose of
the null vector itself

$$
\mathbf{0}^{+} = \mathbf{0}^{\mathrm{T}} \tag{33}
$$

while the pseudoinverse of a non-null vector is

$$
\mathbf{v}^{+} = \frac{1}{\mathbf{v}^{\mathrm{T}} \mathbf{v}} \mathbf{v}^{\mathrm{T}} \qquad (34)
$$

As a consequence, the pseudoinversion of equation
(31) yields

$$
(\mathbf{A}\ddot{\mathbf{y}})^{+} = \mathbf{0}^{\mathrm{T}} \qquad (35)
$$

while that of equation (32) yields

$$
\mathbf{0}^{\mathrm{T}} = (\mathbf{U}\ddot{\mathbf{y}})^{\mathrm{T}} (\mathbf{A}\mathbf{U}^{-1})^{+} \quad (36)
$$

under the assumption that **U**\ddot{**y**} ≠ **0** for **y** ≠ **0**, guar-
anteed by the positive definiteness of matrix **M**. The

evaluation of Gauss' functional of equation (1) with
the accelerations of equation (29) yields

$$
\begin{align*}
G(\ddot{\mathbf{u}}) &= \frac{1}{2}(\ddot{\mathbf{u}} - \mathbf{a})^{\mathrm{T}}\mathbf{M}(\ddot{\mathbf{u}} - \mathbf{a}) \\
&= \frac{1}{2}[(\mathbf{A}\mathbf{U}^{-1})^{\mathrm{T}}(\mathbf{b} - \mathbf{A}\mathbf{a}) + \mathbf{U}\ddot{\mathbf{y}}]^{\mathrm{T}} \\
&\quad \times [(\mathbf{A}\mathbf{U}^{-1})^{\mathrm{T}}(\mathbf{b} - \mathbf{A}\mathbf{a}) + \mathbf{U}\ddot{\mathbf{y}}] \\
&= \frac{1}{2}[(\mathbf{A}\mathbf{U}^{-1})^{\mathrm{T}}(\mathbf{b} - \mathbf{A}\mathbf{a})]^{\mathrm{T}}[(\mathbf{A}\mathbf{U}^{-1})^{\mathrm{T}}(\mathbf{b} - \mathbf{A}\mathbf{a})] \\
&\quad + \frac{1}{2} \underbrace{(\mathbf{U}\ddot{\mathbf{y}})^{\mathrm{T}}(\mathbf{A}\mathbf{U}^{-1})^{\mathrm{T}}}_{=0 \text{ (according to equation (36))} (\mathbf{b} - \mathbf{A}\mathbf{a}) \\
&\quad + \frac{1}{2}(\mathbf{b} - \mathbf{A}\mathbf{a})^{\mathrm{T}} \underbrace{[(\mathbf{A}\mathbf{U}^{-1})^{\mathrm{T}}(\mathbf{U}\ddot{\mathbf{y}})]}_{=0 \text{ (according to equation (36))} \\
&\quad + \frac{1}{2}(\mathbf{U}\ddot{\mathbf{y}})^{\mathrm{T}}(\mathbf{U}\ddot{\mathbf{y}}) \\
&= G(\ddot{\mathbf{x}}) + \frac{1}{2}\ddot{\mathbf{y}}^{\mathrm{T}}\mathbf{M}\ddot{\mathbf{y}} \tag{37}
\end{align*}
$$

The value of $G(\ddot{\mathbf{u}})$ in equation (37) can only be larger
than that of $G(\ddot{\mathbf{x}})$, since the value of $\ddot{\mathbf{y}}^{\mathrm{T}}\mathbf{M}\ddot{\mathbf{y}}$ can only be
positive when $\ddot{\mathbf{y}} \neq \mathbf{0}$, as matrix $\mathbf{M}$ is positive-definite ■.

4.3 Is it really the same formula?

The question whether the proposed formula,
equation (23), and the one by Udwadia and Kalaba,
equation (6), are the same arises. In fact, their proof,
which inspired the one presented in this work, clearly
indicates that the formula of equation (6) is unique.
The next paragraph shows that the two formulae are
exactly the same.

The Moore-Penrose pseudoinverse of a matrix $\mathbf{H} \in \mathbb{R}^{m \times n}$, with $n > m$ and rank $r = m$, is

$$
\mathbf{H}^{+} = \mathbf{H}^{\mathrm{T}} (\mathbf{H}\mathbf{H}^{\mathrm{T}})^{-1} \quad (38)
$$

The requirement on *r* = *m* can be relaxed to *r* ≤ *m*;
this means that the constraints expressed by matrix **H**
are not independent; in this case, the pseudoinverse is
defined as the limit

$$
\mathbf{H}^{+} = \lim_{\epsilon \to 0} \mathbf{H}^{\mathrm{T}} (\mathbf{H}\mathbf{H}^{\mathrm{T}} + \epsilon \mathbf{I})^{-1} \quad (39)
$$

and this limit exists for every matrix **H**.

So, in the case of equation (6), the expression
**M**⁻¹/²(**AM**⁻¹/²)⁺ yields

$$
\begin{align*}
\mathbf{M}^{-1/2}(\mathbf{A}\mathbf{M}^{-1/2})^{+} &= \mathbf{M}^{-1/2}(\mathbf{A}\mathbf{M}^{-1/2})^{\mathrm{T}}[(\mathbf{A}\mathbf{M}^{-1/2})(\mathbf{A}\mathbf{M}^{-1/2})^{\mathrm{T}}]^{+} \\
&= \mathbf{M}^{-1/2}(\mathbf{M}^{-1/2})^{\mathrm{T}}\mathbf{A}^{\mathrm{T}}[\mathbf{A}\mathbf{M}^{-1/2}(\mathbf{M}^{-1/2})^{\mathrm{T}}\mathbf{A}^{\mathrm{T}}]^{+} \\
&= \mathbf{M}^{-1}\mathbf{A}^{\mathrm{T}}(\mathbf{A}\mathbf{M}^{-1}\mathbf{A}^{\mathrm{T}})^{+} \tag{40}
\end{align*}
$$
In the case of equation (23), the expression $\mathbf{U}^{-1}(\mathbf{A}\mathbf{U}^{-1})^{+}$ yields

$$ 
\begin{aligned}
\mathbf{U}^{-1}(\mathbf{A}\mathbf{U}^{-1})^{+} &= \mathbf{U}^{-1}(\mathbf{A}\mathbf{U}^{-1})^{\mathrm{T}}[(\mathbf{A}\mathbf{U}^{-1})(\mathbf{A}\mathbf{U}^{-1})^{\mathrm{T}}]^{+} \\
&= \mathbf{U}^{-1}\mathbf{U}^{-\mathrm{T}}\mathbf{A}^{\mathrm{T}}(\mathbf{A}\mathbf{U}^{-1}\mathbf{U}^{-\mathrm{T}}\mathbf{A}^{\mathrm{T}})^{+} \\
&= \mathbf{M}^{-1}\mathbf{A}^{\mathrm{T}}(\mathbf{A}\mathbf{M}^{-1}\mathbf{A}^{\mathrm{T}})^{+}
\end{aligned}
\quad (41) $$

The two expressions in equations (40) and (41) are clearly the same. So the use of Cholesky's decomposition simply represents a much more computationally efficient way to write the same correction to the unconstrained accelerations proposed by Udwadia and Kalaba.

## 5 WHERE DOES IT COME FROM?

The principle of least constraint requires the minimization of the functional of equation (1), subjected to the linear constraint of equation (3). Consider the functional of equation (1) augmented by the linear constraint of equation (3) using the Lagrange multipliers approach

$$ \tilde{G}(\ddot{\mathbf{x}}, \lambda) = \frac{1}{2}(\ddot{\mathbf{x}} - \mathbf{a})^{\mathrm{T}}\mathbf{M}(\ddot{\mathbf{x}} - \mathbf{a}) + \lambda^{\mathrm{T}}(\mathbf{A}\ddot{\mathbf{x}} - \mathbf{b}) \quad (42) $$

and minimize it by imposing the vanishing of its gradient with respect to both the unknown accelerations $\ddot{\mathbf{x}}$ and the Lagrange multipliers $\lambda$

$$ \frac{\partial \tilde{G}}{\partial \ddot{\mathbf{x}}} = \mathbf{M}(\ddot{\mathbf{x}} - \mathbf{a}) + \mathbf{A}^{\mathrm{T}}\lambda = \mathbf{0} \quad (43a) $$

$$ \frac{\partial \tilde{G}}{\partial \lambda} = \mathbf{A}\ddot{\mathbf{x}} - \mathbf{b} = \mathbf{0} \quad (43b) $$

Solve equation (43a) with respect to $\ddot{\mathbf{x}}$

$$ \ddot{\mathbf{x}} = \mathbf{a} - \mathbf{M}^{-1}\mathbf{A}^{\mathrm{T}}\lambda \quad (44) $$

and use it to solve equation (43b) with respect to $\lambda$

$$ \lambda = -(\mathbf{A}\mathbf{M}^{-1}\mathbf{A}^{\mathrm{T}})^{+}(\mathbf{b} - \mathbf{A}\mathbf{a}) \quad (45) $$

where the pseudoinverse of matrix $\mathbf{A}\mathbf{M}^{-1}\mathbf{A}^{\mathrm{T}}$ is used instead of its inverse, to indicate the fact that the inverse either exists as is, or in the limit sense of equation (39). This allows to account for the case of redundant constraints, where the rank of matrix $\mathbf{A} \in \mathbb{R}^{m \times n}$ is $r \le m$. When the multipliers are finally eliminated by exploiting equation (45) in equation (44), one obtains

$$ \ddot{\mathbf{x}} = \mathbf{a} + \mathbf{M}^{-1}\mathbf{A}^{\mathrm{T}}(\mathbf{A}\mathbf{M}^{-1}\mathbf{A}^{\mathrm{T}})^{+}(\mathbf{b} - \mathbf{A}\mathbf{a}) \quad (46) $$

which is Udwadia-Kalaba's formula after removing the need to compute the non-strictly necessary matrix $\mathbf{M}^{1/2}$, and the pseudoinversion is confined to the place

where it could be required, if ever, as indicated in equations (40) and (41).

If $r < m$, the multipliers are not determined; this is not an issue, because the dynamics of the system are not affected. In fact, the constraints of equation (3) are ideal, and thus the constraint reaction forces do not do work for any admissible motion. The pseudoinversion of matrix $\mathbf{A}\mathbf{M}^{-1}\mathbf{A}^{\mathrm{T}}$ computes the minimal norm value of the multipliers, weighted by matrix $\mathbf{M}^{-1}$, among the values that yield the accelerations of equation (46).

Note that equation (46) can be conveniently re-worked as

$$ \ddot{\mathbf{x}} = (\mathbf{I} - \mathbf{P})\mathbf{a} + \mathbf{M}^{-1}\mathbf{A}^{\mathrm{T}}(\mathbf{A}\mathbf{M}^{-1}\mathbf{A}^{\mathrm{T}})^{+}\mathbf{b} \quad (47) $$

where

$$ \mathbf{P} = \mathbf{M}^{-1}\mathbf{A}^{\mathrm{T}}(\mathbf{A}\mathbf{M}^{-1}\mathbf{A}^{\mathrm{T}})^{+}\mathbf{A} \quad (48) $$

is the non-orthogonal projection ($\mathbf{P}^2 = \mathbf{P}$, $\mathbf{P}^{\mathrm{T}} \neq \mathbf{P}$) that maps the accelerations in the space of the physical coordinates onto the rangespace of the constraints $\mathbf{A}$, weighted by the inertia matrix $\mathbf{M}$. An analogous interpretation was given in reference [25], in a different context.

The physics of the projection provided by matrix $\mathbf{P}$, interpreting equation (48) right to left, is as follows:

(a) matrix $\mathbf{A}$ computes the accelerations $\mathbf{a}$ in the rangespace of the constraints;

(b) matrix $(\mathbf{A}\mathbf{M}^{-1}\mathbf{A}^{\mathrm{T}})^{+}$ is the mass matrix projected in the space of the constraints, so it computes the inertia forces in the space of the constraints resulting from the previously described accelerations;

(c) matrix $\mathbf{A}^{\mathrm{T}}$ projects the inertia forces in the space of the constraints back onto the space of the physical coordinates;

(d) finally, matrix $\mathbf{M}^{-1}$ extracts the accelerations from the previously described inertia forces.

Since $\mathbf{P}$ is a projection matrix, the norm of the constrained accelerations, $\|\ddot{\mathbf{x}}\|$, is always less than the norm of the unconstrained ones, $\|\mathbf{a}\|$, when $\mathbf{b} = \mathbf{0}$, as can be shown by considering a spectral decomposition of matrix $\mathbf{P}$.

It is worth noticing that the symmetry of matrix $\mathbf{M}$ is not required by equation (46), nor anywhere in this section. Although unusual, there might be cases where the equations of motion of a system can conveniently be written in a manner that yields a non-symmetric mass matrix, as pointed out in Appendix 2. What is described in this section applies in those cases as well.

## 6 COMPUTATIONAL COMPLEXITY

The claim for more efficiency comes from a series of considerations, reported in the following. They are based on the assumption that full matrices are considered. This is not usually the case for typical mass
matrices related to the physical coordinates of sys-
tems of bodies. Apart from the trivial case of diagonal
mass matrix, typically the mass matrix of a system
of rigid bodies is block diagonal. Coupling might be
more complex when inertia is formulated in a consis-
tent manner, although typically resulting in a banded
matrix. In any case, a discussion addressing issues
related to sparsity handling would probably exceed the
scope of this work.

The key step required to directly compute the con-
strained system accelerations $\ddot{x}$ is represented by
the inversion of the matrix that yields the correc-
tion to the unconstrained system accelerations. The
proposed formulation of equation (23) exploits an
efficient factorization of matrix $\mathbf{M}$ that is required
anyway by equation (2). The original formulation of
equation (6) needs matrix $\mathbf{M}^{-1/2}$; its efficient compu-
tation would provide a relatively efficient factorization
of matrix $\mathbf{M}$ as well (e.g. the spectral decomposition
of $\mathbf{M}$, or its Cholesky's factorization for the first itera-
tion of Denman and Beavers' algorithm, or its Schur's
decomposition).

The matrices on the left-hand sides of equations (40)
and (41) do not actually need to be computed; in
fact, what needs to be computed is their product by
vector (Aa – b). As a consequence, most of the com-
putational complexity is represented by computing
and pseudoinverting the respective matrices (AM⁻¹/²)
and (AU⁻¹), while the remaining operations would
essentially be O(nm) at most.

The pseudoinverse of the previously mentioned
matrices can be computed using the SVD. Its compu-
tational complexity amounts to 4n²m + 8nm² + 9m³.
With respect to pseudoinversion, the two approaches
would have the same complexity; the one proposed
in this work would still be slightly more efficient since
the cost of AU⁻¹ is n²m, while that of AM⁻¹/² is 2n²m
if M⁻¹/² is available as a n × n matrix, or higher if
it is available as any of the factorizations indicated
previously.

It is worth noticing that the ‘regularized’ Lagrange’s
formulation of equation (46) outperforms both the
above mentioned formulations. This is especially true
when $n \gg m$, and thus a problem with many bod-
ies has a limited number of constraints. In fact,
equation (46) requires the pseudoinversion of matrix
($AM^{-1}A^T$). The computation of this matrix can exploit
the Cholesky’s factorization of matrix $\mathbf{M}$, and the sym-
metry of the resulting matrix, for a complexity of
$n^2m + nm^2$. The complexity of the SVD for the pseu-
doinversion reduces to $12m^3$, thanks to the symmetry
of the matrix.

7 CONCLUDING REMARKS

This work illustrates how the solution proposed by
Udwadia and Kalaba to the problem of constrained

system dynamics, based on Gauss' least constraint
principle, does not actually need to deal with the
square root of the mass matrix, nor it requires the mass
matrix to be symmetric.

A computational approach that preserves the orig-
inal formalism of Udwadia and Kalaba is proposed,
based on Cholesky's decomposition. Its advantages in
terms of computational complexity are highlighted.

An alternative, generic projection that directly com-
putes the constrained system accelerations from the
active forces is also presented. The latter formulation
does not require the symmetry of the mass matrix.

Points of strength are:

(a) no coordinate partitioning is required;

(b) the problem is naturally written in the form of ODEs, thus it can be solved by explicit integration schemes;

(c) redundant and inconsistent constraints can be
handled seamlessly.

Points of weakness are:

(a) the solution resulting from the integration of $\ddot{x}$ may suffer from drift, and needs constraint stabilization;

(b) the formulation needs the second derivative (for holonomic constraints) or the first derivative (for non-holonomic ones) of the constraint equation, which might not be trivial to formulate and evaluate;

(c) the matrix pseudo inversion required by the pro-
posed formula can be computationally expensive,
although less than the original one especially for
large bodies to constraint ratios;

(d) the basic formulation does not handle non-ideal
constraints; extensions have been proposed, lim-
ited to the case of independent non-ideal con-
straints, but no definitive solution exists, to the
authors' knowledge.

Examples of non-symmetric mass matrices result-
ing from problem formulations of general usefulness,
including rigid and deformable body dynamics, are
presented and discussed in Appendix 2.

ACKNOWLEDGEMENTS

The authors gratefully acknowledge the fruitful discus-
sions they had with Professors Ettore Pennestrì and
Paolo Mantegazza on the subject of this work.

REFERENCES

1 Gauss, C. F. Ueber ein neues algemeines Grundgesetz der Mechanik. *J. Reine Angew. Math.*, 1829, **4**, 232–235.

2 Finzi, B. *Meccanica razionale, vol. 2: dinamica* (in Italian), 1959 (Zanichelli, Bologna).
3 Schiehlen, W. Multibody system dynamics: roots and perspectives. *Multibody Syst. Dyn.*, 1997, **1**(2), 149–188.
DOI:10.1023/A:1009745432698.

4 Udwadia, F. E. and Kalaba, R. E. A new perspective on constrained motion. *Proc. R. Soc. A*, 1992, **439**(1906), 407–410.

5 Udwadia, F. E. and Kalaba, R. E. Analytical dynamics, 1996 (Cambridge University Press, New York).

6 de Falco, D., Pennestrì, E., and Vita, L. The Udwadia–Kalaba formulation: a report on its numerical efficiency in multibody dynamics simulations and on its teaching effectiveness. In Proceedings of Multibody Dynamics 2005, ECCOMAS Thematic Conference, Madrid, Spain, 21–24 June 2005.

7 Arabyan, A. and Wu, F. An improved formulation for constrained mechanical systems. *Multibody Syst. Dyn.*, 1998, **2**, 49–69. DOI:10.1023/A:1009724704839.

8 Udwadia, F. E. and Kalaba, R. E. An alternate proof for the equation of motion for constrained mechanical systems. *Appl. Math. Comput.*, 1995, **70**, 339–342.

9 Moore, E. H. On the reciprocal of the general algebraic matrix. *Bull. Am. Math. Soc.*, 1920, **26**, 394–395.

10 Penrose, R. A generalized inverse for matrices. *Proc. Camb. Philos. Soc.*, 1955, **51**, 406–413.

11 Ben-Israel, A. and Greville, T. N. E. Generalized inverses: theory and applications, 1974 (Wiley, New York).

12 Baumgarte, J. Stabilization of constraints and integrals of motion in dynamical systems. *Comput. Methods Appl. Mech.*, 1972, **1**, 1–36.

13 Amirouche, F. M. L. and Tung, C.-W. Regularization and stability of the constraints in the dynamics of multibody systems. *Nonlinear Dyn.*, 1990, **1**(6), 459–475.
DOI:10.1007/BF01856949.

14 de Jalon, J. and Bayo, E. Kinematic and dynamic simulation of multibody systems, 1993 (Springer-Verlag, Heidelberg).

15 Neto, M. A. and Ambrósio, J. Stabilization methods for the integration of DAE in the presence of redundant constraints. *Multibody Syst. Dyn.*, 2003, **10**(1), 81–105.
DOI:10.1023/A:1024567523268.

16 Denman, E. D. and Beavers Jr, A. N. The matrix sign function and computations in systems. *Appl. Math. Comput.*, 1976, **2**(1), 63–94. DOI:10.1016/0096-3003(76)00020-5.

17 Benner, P. and Quintana-Ortí, E. S. Solving stable generalized Lyapunov equations with the matrix sign function. *Numer. Algorithms*, 1999, **20**(1), 75–100. DOI: 10.1023/A:1019191431273.

18 Stykel, T. A modified matrix sign function method for projected Lyapunov equations. *Syst. Control Lett.*, 2007, **56**, 695–701. DOI:10.1016/j.sysconle.2007.06.003.

19 Björck, Å. and Hammarling, S. A Schur method for the square root of a matrix. *Linear Algeb. Appl.*, 1983, **52–53**, 127–140.

20 Higham, N. J. Computing real square roots of a real matrix. *Linear Algeb. Appl.*, 1987, **88–89**, 405–430.

21 Higham, N. J. Stable iterations for the matrix square root. *Numer. Algorithms*, 1997, **15**(2), 227–242. DOI:10.1023/A:1019150005407.

22 Benoit, C. Note sur une méthode de résolution des équations normales provenant de l'application de la méthode des moindres carrés à un système d'équations linéaires en nombre inférieur à celui des inconnues (procédé

du commandant Cholesky). *Bull. Géodésique*, 1924, **2**, 67–77.

23 Redon, S., Kheddar, A., and Coquillart, S. Gauss' least constraints principle and rigid body simulations. In Proceedings of the IEEE International Conference on Robotics and Automation, May 2002, vol. 1, pp. 517–522.
DOI: 10.1109/ROBOT.2002.1013411.

24 Franklin, J. Least-squares solution of equations of motion under inconsistent constraints. *Linear Algebr. Appl.*, 1995, **222**, 9–13. DOI:10.1016/0024-3795(93)00267-4.

25 Tseng, F.-C., Ma, Z.-D., and Hulbert, G. M. Efficient numerical solution of constrained multibody dynamics systems. *Comput. Methods Appl. Mech. Eng.*, 2003, **192**, 439–472. DOI:10.1016/S0045-7825(02)00521-2.

26 Ghiringhelli, G. L., Masarati, P., and Mantegazza, P. A multi-body implementation of finite volume beams. *AIAA J.*, 2000, **38**(1), 131–138.

APPENDIX 1

Notation

<table>
  <tr>
    <td>a &isin; R<sup>n</sup></td>
    <td>unconstrained bodies<br>accelerations</td>
  </tr>
  <tr>
    <td>A &isin; R<sup>m&times;n</sup></td>
    <td>constraint matrix</td>
  </tr>
  <tr>
    <td>b &isin; R<sup>m</sup></td>
    <td>constraint equation<br>right-hand side</td>
  </tr>
  <tr>
    <td>D &isin; C<sup>n&times;n</sup><br>(R<sup>n&times;n</sup> &gt; 0)</td>
    <td>diagonal matrix containing<br>the (positive) eigenvalues of<br>matrix H</td>
  </tr>
  <tr>
    <td>f &isin; R<sup>n</sup></td>
    <td>external loads</td>
  </tr>
  <tr>
    <td>G &isin; R &ge; 0</td>
    <td>Gauss' functional (the<br>'Gaussian')</td>
  </tr>
  <tr>
    <td>H &isin; R<sup>n&times;n</sup> (&gt;0)</td>
    <td>generic (positive definite)<br>matrix</td>
  </tr>
  <tr>
    <td>m &isin; N</td>
    <td>number of constraint<br>equations</td>
  </tr>
  <tr>
    <td>M &isin; R<sup>n&times;n</sup> sym &gt; 0</td>
    <td>inertia matrix</td>
  </tr>
  <tr>
    <td>n &isin; N</td>
    <td>order of the problem</td>
  </tr>
  <tr>
    <td>P &isin; R<sup>n&times;n</sup></td>
    <td>constrained accelerations<br>projection matrix</td>
  </tr>
  <tr>
    <td>Q &isin; R<sup>n&times;n</sup></td>
    <td>unitary matrix from Schur<br>decomposition of M</td>
  </tr>
  <tr>
    <td>r &isin; N</td>
    <td>number of independent<br>constraint equations (rank of<br>matrix A; r &le; m)</td>
  </tr>
  <tr>
    <td>S &isin; R<sup>n&times;n</sup></td>
    <td>(quasi-)upper triangular<br>matrix from Schur<br>decomposition of M</td>
  </tr>
  <tr>
    <td>t &isin; R</td>
    <td>the time</td>
  </tr>
  <tr>
    <td>T &isin; R<sup>n&times;n</sup></td>
    <td>square root of matrix S</td>
  </tr>
  <tr>
    <td>ii &isin; R<sup>n</sup></td>
    <td>perturbed accelerations</td>
  </tr>
  <tr>
    <td>U &isin; R<sup>n&times;n</sup></td>
    <td>Cholesky's decomposition of<br>matrix M</td>
  </tr>
  <tr>
    <td>v &isin; R<sup>n</sup></td>
    <td>generic vector</td>
  </tr>
  <tr>
    <td>V &isin; C<sup>n&times;n</sup> (R<sup>n&times;n</sup>)</td>
    <td>matrix of (real) eigenvectors<br>of H</td>
  </tr>
</table>
$\boldsymbol{w} \in \mathbb{R}^n$ normalized accelerations
($\boldsymbol{w} = \mathbf{U}\ddot{\boldsymbol{x}}$)

$\boldsymbol{x} \in \mathbb{R}^n$ coordinates of the bodies
$\mathbf{X} \in \mathbb{R}^{n \times n}$ matrix transformation
between test and trial
variables

$\ddot{\boldsymbol{y}} \in \mathbb{R}^n$ acceleration perturbation
$\mathbf{Y} \in \mathbb{R}^{n \times n}$ estimate of $\mathbf{M}^{1/2}$ in Denman
and Beavers iteration

$\mathbf{Z} \in \mathbb{R}^{n \times n}$ estimate of $\mathbf{M}^{-1/2}$ in Denman
and Beavers iteration

$\varepsilon \in \mathbb{R} > 0$ small positive quantity

$\boldsymbol{\lambda} \in \mathbb{R}^m$ Lagrange's multipliers vector

$\boldsymbol{\phi} \in \mathbb{R}^m$ holonomic constraints vector

$\boldsymbol{\psi} \in \mathbb{R}^m$ non-holonomic constraints
vector

*Note*: symbols used only in the examples are defined in each specific section

# APPENDIX 2

## Non-Symmetric Inertia Matrix

To further generalize the discussion of Gauss' principle of least constraint, consider an inertia matrix $\mathbf{M}^*$ that, although non-singular, is not symmetric. Non-symmetric inertia matrices can result from Petrov-Galerkin approaches where different trial and test functions are used; examples are given at the end of this section.

Since in Lagrangian mechanics an inertia matrix $\mathbf{M}$ results from a consistent definition of the kinetic coenergy of a mechanical system, it must be symmetric. However, there might be reasons to express the equations of motion of the system using a set of trial variables, $\boldsymbol{x}^*$, that differ from the test variables $\boldsymbol{x}$. They are related by the invertible relationship

$$ \boldsymbol{x}^* = \boldsymbol{x}^*(\boldsymbol{x}) \quad (49) $$

whose second-order derivative is linear in the second-order derivative of $\boldsymbol{x}$

$$ \begin{aligned} \ddot{\boldsymbol{x}}^* &= \frac{\partial \boldsymbol{x}^*}{\partial \boldsymbol{x}} \ddot{\boldsymbol{x}} + \frac{\partial}{\partial \boldsymbol{x}} \left( \frac{\partial \boldsymbol{x}^*}{\partial \boldsymbol{x}} \dot{\boldsymbol{x}} \right) \dot{\boldsymbol{x}} \\ &= \mathbf{X}^{-1} \ddot{\boldsymbol{x}} + (\boldsymbol{x}^*_{/x} \dot{\boldsymbol{x}})_{/x} \dot{\boldsymbol{x}} \end{aligned} \quad (50) $$

Since the transformation is invertible, matrix $\mathbf{X}$ is non-singular. The non-symmetric mass matrix results from this transformation, namely

$$ \mathbf{M}^* = \mathbf{M}\mathbf{X} \quad (51) $$

If equations (50) were known, the symmetry of the matrix would be restored by substituting it into the

definition of the fraction of the inertia forces that depend on the second derivative of $\boldsymbol{x}^*$, namely

$$ \mathbf{M}^* \ddot{\boldsymbol{x}}^* - \mathbf{M}\mathbf{X}(\boldsymbol{x}^*_{/x} \dot{\boldsymbol{x}})_{/x} \dot{\boldsymbol{x}} = \mathbf{M}\mathbf{X}\mathbf{X}^{-1} \ddot{\boldsymbol{x}} \quad (52) $$

and thus

$$ \mathbf{M} = \mathbf{M}^* \mathbf{X}^{-1} = \mathbf{X}^{-T} (\mathbf{M}^*)^T \quad (53) $$

Similarly, the same fraction of the inertia forces could be rewritten as generalized forces conjugated to a virtual variation of $\boldsymbol{x}^*$

$$ \delta \boldsymbol{x}^* = \boldsymbol{x}^*_{/x} \delta \boldsymbol{x} \quad (54) $$

by premultiplying them by the transpose of the virtual displacement $\delta \boldsymbol{x}$, namely

$$ \delta \boldsymbol{x}^T \mathbf{M}^* \ddot{\boldsymbol{x}}^* = (\delta \boldsymbol{x}^*)^T \mathbf{X}^T \mathbf{M}^* \ddot{\boldsymbol{x}}^* \quad (55) $$

the corresponding inertia matrix would be

$$ \mathbf{M}^{*,*} = \mathbf{X}^T \mathbf{M}^* = (\mathbf{M}^*)^T \mathbf{X} \quad (56) $$

To be consistent with the relationship between the trial and the test functions of equation (49), the augmented Gaussian of equation (42) needs to be redefined as

$$ G(\ddot{\boldsymbol{x}}^*(\boldsymbol{x}), \boldsymbol{\lambda}) = \frac{1}{2}(\ddot{\boldsymbol{x}} - \boldsymbol{a})^T \mathbf{M}^*(\ddot{\boldsymbol{x}}^* - \boldsymbol{a}^*) + \boldsymbol{\lambda}^T \left( \underbrace{\mathbf{A}\boldsymbol{x}}_{\mathbf{A}^*} \ddot{\boldsymbol{x}}^* - \mathbf{b} \right) \quad (57) $$

with

$$ \begin{aligned} \boldsymbol{a}^* &= (\mathbf{M}^*)^{-1} \boldsymbol{f} + (\boldsymbol{x}^*_{/x} \dot{\boldsymbol{x}})_{/x} \dot{\boldsymbol{x}} \\ &= \mathbf{X}^{-1} \boldsymbol{a} + (\boldsymbol{x}^*_{/x} \dot{\boldsymbol{x}})_{/x} \dot{\boldsymbol{x}} \end{aligned} \quad (58) $$

Recalling equations (56) and (50), equations (43) become

$$ \mathbf{M}^*(\ddot{\boldsymbol{x}}^* - \boldsymbol{a}^*) + \mathbf{A}^T \boldsymbol{\lambda} = \mathbf{0} \quad (59a) $$

$$ \mathbf{A}^* \ddot{\boldsymbol{x}}^* = \mathbf{b} \quad (59b) $$

Equations (59) are identical to equations (43) when equations (50) and (58) are taken into account. After substituting $\ddot{\boldsymbol{x}}^*$ from equation (59a) into equation (59b) to compute the multipliers $\boldsymbol{\lambda}$, one obtains

$$ \ddot{\boldsymbol{x}}^* = \boldsymbol{a}^* - (\mathbf{M}^*)^{-1} \mathbf{A}^T (\mathbf{A}^*(\mathbf{M}^*)^{-1} \mathbf{A}^T)^+ (\mathbf{A}^* \boldsymbol{a}^* - \mathbf{b}) \quad (60) $$

which, premultiplied by $\mathbf{X}$, would yield equation (46). Note that matrix $\mathbf{M}$ does not need to be actually built from equation (53); in general, this is not possible at all, since matrix $\mathbf{X}$ might be unavailable. In this case equation (60) needs to deal with a non-symmetric mass matrix, thus requiring a linear solver for generic
matrices. For typical multibody problems, an efficient sparse LU factorization would fit.

Few examples follow, which illustrate how non-diagonal, non-symmetric inertia matrices can arise from the formulation of non-trivial constrained dynamics problems.

### Example: Newton–Euler equations

Consider the Newton–Euler formulation of the equations of motion of a rigid body

$$m\ddot{x}_{\text{CM}} = f \qquad (61a)$$

$$J_{\text{CM}}\dot{\omega} + \omega \times J_{\text{CM}}\omega = m_{\text{CM}} \qquad (61b)$$

where $m$ is the mass, $x_{\text{CM}}$ is the location of the centre of mass, $f$ are the external forces, $J_{\text{CM}}$ is the inertia tensor with respect to the centre of mass, projected in an inertial reference frame and thus depending on the parameters used to describe the orientation of the body, $\omega$ is the angular velocity of the body, and $m_{\text{CM}}$ are the external moments, including those of the external forces with respect to the centre of mass. Consider now an arbitrary parametrization $g$ of the rotations, such that

$$R = R(g) \qquad (62a)$$

$$\omega = \Gamma(g)\dot{g} \qquad (62b)$$

where $R$ is the orthonormal rotation tensor and $\omega$ is the angular velocity, while matrix $\Gamma(g)$ expresses the relationship between the angular velocity and the derivatives of the rotation parameters, with (The operator $\cdot \times$ is used to indicate the matrix that produces the vector cross-product when multiplies a vector: $a \times b$ when multiplied by vector $b$)

$$\omega \times = \dot{R}R^T \qquad (63)$$

The equations of motion, in the variables $x_{\text{CM}}$ and $g$, become

$$\begin{bmatrix} mI & \mathbf{0} \\ \mathbf{0} & J_{\text{CM}}\Gamma(g) \end{bmatrix} \begin{Bmatrix} \ddot{x}_{\text{CM}} \\ \ddot{g} \end{Bmatrix} = \begin{Bmatrix} f \\ m_{\text{CM}} - J_{\text{CM}}\dot{\Gamma}(g)\dot{g} - \omega \times J_{\text{CM}}\omega \end{Bmatrix} \qquad (64)$$

where the inertia matrix is clearly non-diagonal and non-symmetric because of the $J_{\text{CM}}\Gamma(g)$ block-diagonal term. In this case the transformation matrix $X$ is known; for each body

$$x^* = \begin{Bmatrix} x_{\text{CM}} \\ g \end{Bmatrix} \qquad (65)$$

$$M = \begin{bmatrix} mI & \mathbf{0} \\ \mathbf{0} & J_{\text{CM}} \end{bmatrix} \qquad (66)$$

$$X = \begin{bmatrix} I & \mathbf{0} \\ \mathbf{0} & \Gamma(g) \end{bmatrix} \qquad (67)$$

A symmetric inertia matrix would result when applying the Virtual Work Principle, Lagrange's equations of the first kind or similar formulations. In that case, the moment equations would be premultiplied by $\omega_{\dot{g}}^T = \Gamma(g)^T$, resulting in

$$\begin{bmatrix} mI & \mathbf{0} \\ \mathbf{0} & \Gamma(g)^T J_{\text{CM}} \Gamma(g) \end{bmatrix} \begin{Bmatrix} \ddot{x}_{\text{CM}} \\ \ddot{g} \end{Bmatrix} = \begin{Bmatrix} f \\ \Gamma(g)^T (m_{\text{CM}} - J_{\text{CM}}\dot{\Gamma}(g)\dot{g} - \omega \times J_{\text{CM}}\omega) \end{Bmatrix} \qquad (68)$$

In this latter case, the inertia matrix is symmetric, but the moment equations no longer truly express the equilibrium of the moments applied to the body; they rather express the equilibrium of the generalized forces conjugated to the orientation parameters $g$. Their linearization is further complicated by the need to take into account the $\Gamma(g)^T$ term that premultiplies the moment equations. Also for this reason, implicit formulations of the equations of motion that would in any case yield a non-symmetric Jacobian matrix may find the use of Newton–Euler equations appealing.

### Example: finite volume beam

Consider the finite volume beam formulation proposed in [26]. It basically consists of writing the force and moment equilibrium of finite portions of beam with respect to the node they are related to, including the internal forces and moments at the boundaries of the finite portion.

The equilibrium of an infinitesimal beam slice subjected to distributed loads yields

$$\left\{ \begin{matrix} f \\ m \end{matrix} \right\}_{/\xi} + \left[ \begin{matrix} \mathbf{0} \\ p_{/\xi} \times \mathbf{0} \end{matrix} \right] \left\{ \begin{matrix} f \\ m \end{matrix} \right\} = \left\{ \begin{matrix} \tau \\ \mu \end{matrix} \right\} \qquad (69)$$

where $f$ and $m$ are the internal force and moment vectors, respectively, while $\tau$ and $\mu$ are the distributed forces and moments, and $p(\xi)$ is the position of the generic point at abscissa $\xi$ along the reference line $p$ of the beam. An interesting weak form is obtained by integrating equation (69) weighted by piecewise constant functions. Integration by parts shows that this corresponds to writing the equilibrium of finite portions of beam, subjected to the external loads and to the internal forces and moments evaluated at the span-wise boundaries of each finite portion. The latter, in turn, can be expressed as functions of the beam kinematics by means of appropriate interpolation functions.

The distributed inertia forces and moments $\tau$ and $\mu$ acting on the node $n$, located at $\xi_n$ and comprised between abscissae $a$ and $b$ of each beam portion, yield
the external forces and moments

$$ \left[ \begin{matrix} \mathbf{I} \\ (\mathbf{p}(\xi) - \mathbf{p}(\xi_n)) \times \mathbf{I} \end{matrix} \right] \left\{ \begin{matrix} \mathbf{0} \\ \mathbf{m} \end{matrix} \right\}_a^b = \int_a^b \left[ \begin{matrix} \mathbf{I} \\ (\mathbf{p}(\xi) - \mathbf{p}(\xi_n)) \times \mathbf{I} \end{matrix} \right] \left\{ \begin{matrix} \boldsymbol{\tau} \\ \boldsymbol{\mu} \end{matrix} \right\} \frac{d\mathbf{p}}{d\xi} d\xi \quad (70) $$

and

$$ \frac{d\mathbf{p}}{d\xi} = \sqrt{\frac{d\mathbf{p}^T}{d\xi} \frac{d\mathbf{p}}{d\xi}} \quad (71) $$

Consider the distributed inertia forces and moments

$$ \left\{ \begin{matrix} \boldsymbol{\tau} \\ \boldsymbol{\mu} \end{matrix} \right\} = - \left[ \begin{matrix} \bar{m} \mathbf{l} & \bar{s} \times \mathbf{T} \\ \bar{s} \times & \bar{j} \end{matrix} \right] \left\{ \begin{matrix} \ddot{\boldsymbol{p}} \\ \dot{\omega} \end{matrix} \right\} \quad (72) $$

They result from multiplying the reference linear and angular accelerations of the beam section, respectively, $\ddot{\boldsymbol{p}}$ and $\dot{\omega}$, by the sectional inertia properties (mass $\bar{m}$ and first and second inertia moments $\bar{s}$ and $\bar{j}$ per unit length, the latter evaluated with respect to the reference line, which need not be coincident with the centre of mass). In this case, equation (70) yields the nodal inertia forces and moments.

For example, when the linear and angular accelerations at the generic point of abscissa $\xi$ are linearly

interpolated from nodal values, indicated by the subscript $n$, by means of regular shape functions $N_n(\xi)$, namely

$$ \boldsymbol{p} = N_n(\xi) \boldsymbol{p}_n \quad (73a) $$

$$ \boldsymbol{g} = N_n(\xi) \boldsymbol{g}_n \quad (73b) $$

where $\boldsymbol{p}$ and $\boldsymbol{g}$ are generic position and orientation parameters, and thus

$$ \ddot{\boldsymbol{p}} = N_n(\xi) \ddot{\boldsymbol{p}}_n \quad (74a) $$

$$ \dot{\boldsymbol{\omega}} = \boldsymbol{\Gamma}[\boldsymbol{g}(\xi)] N_n(\xi) \ddot{\boldsymbol{g}}_n + \boldsymbol{\gamma}[\boldsymbol{g}(\xi), \dot{\boldsymbol{g}}(\xi)] \quad (74b) $$

this expression naturally leads to a non-symmetric form for the inertia matrix of the corresponding finite element of beam. In fact, even when $\bar{s} = \mathbf{0}$ equation (70) clearly shows that the nodal force $\mathbf{f}$ only depends on the distributed inertia force, while the nodal moment $\boldsymbol{m}_n$ depends on both the distributed inertia force and moment because of the cross product by the arm $\boldsymbol{p}(\xi) - \boldsymbol{p}(\xi_n)$.

In this latter example, the transformation matrix $\mathbf{X}$ that would allow to symmetrize the inertia matrix cannot be explicitly formulated in a trivial manner. For this reason, the capability to deal with non-symmetric, yet consistently formulated inertia matrices can be important.