MULTIBODY DYNAMICS 2005, ECCOMAS Thematic Conference J.M. Goicolea, J. Cuadrado, J.C. Garcı́a Orden (eds.) Madrid, Spain, 21–24 June 2005 THE UDWADIA-KALABA FORMULATION: A REPORT ON ITS NUMERICAL EFFICIENCY IN MULTIBODY DYNAMICS SIMULATIONS AND ON ITS TEACHING EFFECTIVENESS D. de Falco⋆ , E. Pennestrı̀† , and L. Vita† ⋆ Dipartimento di Ingegneria Spaziale e Aeronautica II Università di Napoli Real Casa dell’Annunziata, Via Roma,29 81031 Aversa - Italy e-mail: domenico.defalco@unina2.it † Dipartimento di Ingegneria Meccanica Università di Roma Tor Vergata Facoltà Ingegneria, via del Politecnico, 1 00133 Roma, Italy e-mails: pennestri@mec.uniroma2.it, web page: www.ingegneriameccanica.org vita@ing.uniroma2.it Keywords: Pseudoinverse matrix computation, Udwadia-Kalaba dynamic formulation, Multibody dynamics, Least squares. Abstract. The formulation of the dynamic equations of motion proposed by Udwadia-Kalaba is discussed from the point of view of teaching effectiveness and numerical efficiency. Since this formulation requires the computation of a pseudoinverse matrix, it is investigated the influence of the method of pseudoinverse computation on the dynamic simulation of an overconstrained linkage. 1 D. de Falco, E. Pennestrı̀ and L. Vita 1 INTRODUCTION In a recent series of papers (e.g. [1]) and in a textbook [2] Udwadia and Kalaba, starting from the Gauss Principle of Least Constraint [3, 5, 6], deduced a new formulation of the dynamics equations of motion for a system of constrained particles or rigid bodies. The main features of this formulation are: - the equations of motion can be reduced to a system of ordinary differential equations (ODE), even when a redundant set of coordinates is used; - variations of degrees-of-freedom due to the change of topology or other causes are allowed and do not require special effort in computer programming; - rheonomic and scleronomic constraints are treated in the same way; - forward and inverse dynamics problems can be solved within the same tool; - easy computer implementation, provided that a subroutine for computing the pseudoinverse matrix is available. The above features make the formulation very attractive. In fact, one of the shortcomings associated with the use of a redundant set of coordinates is the integration of a mixed system of differential-algebraic equations (DAE). These systems are different from ordinary differential equations (ODE) and require specialised numerical methods for their solution. The matrices are of high order, but sparse. The computer programming effort required by redundant coordinates dynamic formulations is relatively low and within the capabilities of the average mechanical engineering student. A fundamental step of the Udwadia-Kalaba (UK) formulation is the computation of the Moore-Penrose generalized inverse or pseudoinverse matrix. Arabyan and Wu [7] investigated numerically the advantages of this formulation within the framework of multibody dynamics and compared the number of floating point operations (flops) required for computing the pseudoinverse with Gram-Schmidt (GS) orthogonalization and Singular Value Decomposition (SVD). However, the authors of the present paper believe that the comparison of the methods for computing the pseudoinverse only on the basis of flops operations required, although important, is somewhat limited. For this reason we decided to compare different methods by observing other parameters such as the reliability, accuracy of results and speed of computation. Thus, one of the purposes of this paper is to report about numerical tests, for accuracy and computational efficiency, of different pseudoinverse matrix calculation algorithms with explicit reference to multibody dynamics simulations. These tests do not seem to be readily available in literature. The paper is organized in four parts. The first and secont parts are mainly tutorial. In the first part the main steps of the following methods for computing the pseudoinverse matrix are summarized: - Singular value decomposition (SVD) [22]; - Varga’s algorithm [19] with QR decomposition computed by means of the Householder factorization [9]; - Greville’s algorithm [20, 21]; 2 D. de Falco, E. Pennestrı̀ and L. Vita - Least squares solution with modified Gram-Schmidt (GS) QR factorization [22]; - Least squares solution with Householder QR factorization. [22] In the second part the formulation of Udwadia-Kalaba is concisely deduced. The third part describes the test problem and reports the numerical results. Finally, the fourth part of the paper discuss some issues related to the teaching of the UK formulation in a multibody dynamics course. 2 REVIEW OF METHODS FOR COMPUTING THE MOORE-PENROSE PSEUDOINVERSE MATRIX This section, for completeness, summarizes the main properties of the Moore-Penrose pseudoinverse matrix and the steps of the different algorithms tested during the dynamic simulations. 2.1 Definitions The main properties of the Moore-Penrose pseudoinverse matrix [A]+ of a matrix [A] are: • [A] [A]+
T = [A] [A]+
T • [A]+ [A] = [A]+ [A] • [A] [A]+ [A] = [A] • [A]+ [A] [A]+ = [A]+ When [A] is a square matrix with full rank, then its pseudoinverse coincide with the inverse. The Moore-Penrose pseudoinverse matrix is associated with the least squares solution of the linear system of equations [A] {x} = {b} , (1) where the number m of equations is not equal to the number n of unknowns and [A] does not necessarily has full rank. In particular, the following cases are distinguished: Overdetermined system of equations (m > n) By requiring that (2) h ≡ k[A] {x} − {b}k22 is a minimum, one obtains [A]T [A] {x} = [A]T {b} . (3) Therefore, the solution of (1) can be stated as where {x} = [A]+ {b} , (4) −1 [A] = [A] [A] [A]T (5) +  T is the right pseudoinverse matrix. Undetermined system of equations (m < n) 3 D. de Falco, E. Pennestrı̀ and L. Vita The solution is obtained imposing the minimum of the Euclidean norm g ≡ kxk22 , (6) with {x} subjected to (1). Thus, introduced the new objective function, g ′ ≡ g + {λ}T ([A] {x} − {b}) the solution is achieved solving the system      0 x I AT , = b λ A 0 (7) (8) or where {x} = [A]+ {b} , (9) −1  T [A] = [A] [A] [A] , (10) + is the left pseudoinverse matrix. T 2.2 The least squares method Since there is abundance of software procedures for computing the least squares solution of a system of algebraic equations [A] {x} = {b} , (11) the computation of the pseudoinverse can be reduced to such solution. These procedures are often based on QR decomposition by means of Householder reflections or GS orthogonalization. Let  T {b1 } = 1 0 0 . . . 0 ,  T {b2 } = 0 1 0 . . . 0 ,  T {b3 } = 0 0 1 . . . 0 , ......  T {bm } = 0 0 0 . . . 1 . The procedure differs according to the dimensions of [A]. Case m > n 1. Solve m times1 the system [A]T [A] {x} = [A]T {b} . (12) 2. From the m solutions {x1 }, {x2 },. . . ,{xm } one can form the pseudoinverse matrix   [A]+ = {x1 } {x2 } {x3 } . . . {xm } (13) 1 Computationally, the system is solved only once. The factor matrices used for the first solution are saved and used also for the remaining ones. 4 D. de Falco, E. Pennestrı̀ and L. Vita Case n ≥ m 1. Solve m times the system  I AT A 0  x λ  =  0 b  (14) Also in this case the pseudoinverse is given by equation (13). The square matrices in (12) and (14) are singular or ill conditioned. Thus their solution requires special care. As mentioned, an appropriate way to solve least squares problems is by means of GS orthogonalization or Householder QR factorization. 2.3 Singular Value Decomposition method 1. Apply the Singular Value Decomposition (SVD) to matrix [A] [A] = [U] [Σ] [V ]T , where [U] [U]T = [I]m×m , [V ] [V ]T = [I]n×n is a diagonal matrix. 0 σ1  ..     [Σ] =     0 0 .. . . σn ··· 0 . .. . .. ··· 0 0          2. Compute [A]+ = [V ] [Σ]+ [U]T , where  σ1+     + [Σ] =     and ( 0 σi+ = σi−1 for i = 1, 2, . . . , n. 5 0 0 .. . 0 0 .. . σn+ ··· 0 . .. . .. ··· 0 (σi = 0) (σi 6= 0) (15)          (16) D. de Falco, E. Pennestrı̀ and L. Vita 2.4 Varga’s method 1. Apply the QR factorization to matrix [A]. Thus, one obtains:     G1 [A] = [U] [G] = U1 U2 0 with [G1 ] matrix with full row rank and [U]T [U] = [Im ]. 2. Apply the QR factorization to [G1 ]T T T [G1 ] = [V ]  GT2 0  , (17) with [G2 ]T invertible matrix and [V ] [V ]T = [In ]. 3. Since [A] = [U]  G2 0 0 0  [V ] , the pseudomatrix is readily obtained from   −1 G2 0 + T [A] = [V ] [U]T . 0 0 (18) 2.5 Greville’s method 1. Decompose the matrix [A]m×n into row vectors {ai }, (i = 1, 2, . . . , m)   [A] = aT1 aT2 · · · aTm 2. Let matrix [Ai ]i×n = with [A1 ]1×n = {a1 }1×n .  Ai−1 ai  3. For i = 2, . . . , m compute the matrices [Ai ]+ as i h + T T = [Ai ]+ , [Ai−1 ] − {bi } {di } {bi } n×i where {di }1×(i−1) = {ai } [Ai−1 ]+ {ci }1×n = {ai } − {di } [Ai−1 ]  {c } i  (kci k = 6 0) {ci }{ci }T {bi }1×n = {di }[A+ ]T i−1  (kci k = 0) 1+{di }{di }T ( {a1 }T (ka1 k 6= 0) T + {a 1 }{a1 } [A1 ] = {a1 }T (ka1 k = 0) 4. After m repetitions [Am ]+ gives the pseudoinverse [A]+ n×m of matrix [A]. 6 (19) D. de Falco, E. Pennestrı̀ and L. Vita 3 THE FORMULATION OF UDWADIA-KALABA In this section the formulation of Udwadia-Kalaba is deduced following the approach of Arabyan and Wu [7]. Let us denote with • {F } the vector of external generalized forces; • [M] the mass matrix; • [Ψq ] the Jacobian of constraints equations; • {q} the vector of generalized coordinates; • {γ} = − ([Ψq ] {q̇})q {q̇} − 2 [Ψqt ] {q̇} − {Ψtt }; • the upperscript + the operation of pseudoinverse of a matrix. When a redundant set of coordinates is used, the following system of differential-algebraic system of equations (DAE) is obtained      F q̈ M ΨTq (20) = Ψq 0 γ λ The square matrix at the left-hand can be inverted by partitioning [23]
−1 #
−1 −1 " −1  M ΨTq Ψq M −1 M −1 ΨTq Ψq M −1 ΨTq M − M −1 ΨTq Ψq M −1 ΨTq
−1
−1 . = Ψq 0 Ψq M −1 − Ψq M −1 ΨTq Ψq M −1 ΨTq (21) Introduced the vector {q̈f } = [M]−1 {F } , (22) representing the acceleration vector of the unconstrained system, the solution of (20) follows −1  ({γ} − [Ψq ] {q̈f }) {q̈} = {q̈f } + [M]−1 [Ψq ]T [Ψq ] [M]−1 [Ψq ]T and  {λ} = [Ψq ] [M] If we let −1 [Ψq ] T −1 ({γ} − [Ψq ] {q̈f }) . 1 (24) 1 [M]−1 = [M]− 2 [M]− 2 , [D] = [Ψq ] [M] (23) − 21 , (25) (26) then equation (23) can be rewritten in the form {q̈} = {q̈f } + [M] − 21 1 = {q̈f } + [M]− 2   −1 − 21 T − 21 − 21 T [M] [Ψq ] [Ψq ] [M] [M] [Ψq ] ({γ} − [Ψq ] {q̈f }) ,  −1 [D]T [D] [D]T ({γ} − [Ψq ] {q̈f }) (27) 7 D. de Falco, E. Pennestrı̀ and L. Vita Taken into account the definition of right pseudoinverse matrix formerly stated, we can let −1  [D]+ = [D]T [D] [D]T , (28) and the previous equation can be concisely expressed as follows [7]: 1 {q̈} = {q̈f } + [M]− 2 [D]+ ({γ} − [Ψq ] {q̈f }) . (29) This is the dynamic formulation originally proposed by Udwadia and Kalaba [2] starting from Gauss Principle of Least Constraint. If the Baumgarte stabilization is introduced, then (29) is modified as follows  n o  + − 21 2 (30) {q̈} = {q̈f } + [M] [D] {γ} − 2α Ψ̇ − β {Ψ} − [Ψq ] {q̈f } , where α and β are the gain parameters usually chosen such that α = β. 4 THE PROBLEM OF RANK COMPUTATION A serious numerical problem that arises in the computation of the pseudoinverse is the accurate computation of the matrix rank. At the practical level, the user must establish a numerical treshold under which the numerical values are all considered to be zero. One may set this value equal to the hardware precision. However, since the inverse of these very small values is required by the algorithms, this inflates a perturbation in the numerical solution. Thus, a truncation or other means of regularization are compulsory to obtain acceptable results. In all the examples discussed in the paper the numerical treshold necessary for the computation of the pseudoinverse is set equal to the error tolerance parameter TOL required by the following IMSL numerical integration subroutines • DIVPAG, which solve an ODE initial-value problem using the Adams-Moulton or Gear’s method. • DIVPRK, which solve an ODE initial-value problem using the Runge-Kutta-Verner fifth and sixth order method. During the integration process, both subroutines try to control the norm of the local error such that the global error is proportional to TOL. It should be acknowledged that the adoption of different criterion may alter the results herein presented. 5 SUMMARY OF FORTRAN SUBROUTINES USED Whenever possible the computation of the pseudoinverse has been executed by means of standard math libraries, such as IMSL. For instance, the Householder QR factorization has been carried out through the combined use of IMSL library subroutines DLQRRR and DLQERR. For other relevant matrix computations, the Fortran code available on the web site www.netlib.org, such as LAPACK, has been used. All the Fortran subroutines described in this section and implemented by the authors are available upon request. • Pseudoinverse by means of the SVD method The IMSL library [9] supplies the readyto-use DLSGRR subroutine. 8 D. de Falco, E. Pennestrı̀ and L. Vita • Pseudoinverse by means of the Varga’s method The two Householder QR factorizations required are executed by means of the combined use of the IMSL library subroutines DLQRRR and DLQERR. For the inversion of the triangular matrix G2 the IMSL routine DLINRT is used. Computation of transpose and matrix multiplications are done with DTRNRR and DMRRRR, respectively. • Pseudoinverse by means of Greville’s method A Fortran subroutine has been implemented with matrix operations executed with LAPACK subroutines [4]. • Pseudoinverse by means of Least squares solution and Householder QR factorization For the least squares solution of a linear system of equations, it has been used the IMSL subroutines DLQRRR and DLQRSL. DLQRRR provides information on the QR factorization, whereas DLQRSL compute the least squares solution. • Pseudoinverse by means of Least squares solution and modified Gram-Schmidt QR factorization For the least squares solution of a linear system of equations, it has been used the code written by R.H. Wapler [10]. 6 NUMERICAL EXAMPLE The example chosen to test the different algorithms of pseudoinverse computation is the overconstrained mechanism shown in Figure 1. The length, mass and moment of inertia of each of the three parallel links are L=1, m=1 and I=0.1, respectively. The length, mass and moment of inertia of the coupler are L2 =2, m2 =2, I2 =0.2. The kinematic constraints were modelled by means of the method described in Nikravesh textbook [18]. The position of each body is described by three coordinates, thus there are n=12 coordinates and k=12 algebraic constraints due to the presence of six revolute pairs. Due to the particular link dimension, Jacobian rank is r = 11. Thus the degree-of-freedom is F = n−r = 1. Gravity force is included and dissipative effects neglected. The results were monitored for tf =20 s of simulation time. Sometimes the execution was halted because of precision loss or other causes. The results of the numerical tests are summarized in Tables 2-7. The simulations were carried out with a PC equipped with AMD Athlon 2600+, 512 MB RAM, running under Microsoft XP Professional. The CPU time reported is the average of the CPU times measured executing four times the same simulation. 7 TEACHING EXPERIENCES In this section are discussed some teaching issues related with the actual implementation at the classroom level of the Udwadia-Kalaba formulation. The current formulations in the syllabus of the course Computational Kinematics and Dynamics (Cinematica e Dinamica Computazionale), offered at University of Roma Tor Vergata, are: • Coordinate partitioning (CP) [11, 17, 18, 12]; • Orthogonalization of constraints by means of QR [13] and SVD [14, 15]; • Udwadia-Kalaba (UK) [2]. At textbook level, a thoughtful discussion of the differences between coordinate partitioning method and formulations based on the orthogonalization of constraints is presented by Bayo and de Jalon [16]. 9 ✂✞ ✞ ✄ ✞ ✄ ✄ ✄ ✆ ✁ ☎ ✂ ✂ ✂ ✆ ☎ ✁✁✂✝ ☎ ✆ ✝ ✄ ✝ D. de Falco, E. Pennestrı̀ and L. Vita Figure 1: Five bar parallelogram linkage Table 1: Comparison of dynamic formulations on the basis of students opinion Linear algebra background required Algebra manipulation effort Computer programming effort Dynamic formulation CP SVD QR Average Extra Extra Medium Medium Medium Medium Medium Medium UK Extra High Low From the didactic point of view, the coordinate partitioning method is very appealing. It is very intuitive and the students grasp it very quickly. In fact, from their very first year of engineering curriculum they became acquainted both with theory and computer programming of the Gaussian elimination method. This method can be applied to partition the set of coordinates. The comprehension of the theoretical bases of the remaining mentioned methods (i.e. SVD, QR and Moore-Penrose inverse) requires an extra knowledge of linear algebra. The amount algebraic manipulation work required to obtain the dynamic formulation based on coordinate partitioning is also lighter when compared to the others based on the orthogonalization of constraints. In particular the algebraic manipulation necessary to deduce the equations of dynamics in the form proposed by Udwadia-Kalaba is not straightforward or obvious. The approach presented by Arabyan and Wu [7] (summarized in a previous section) because of its teaching effectiveness is recommended. On the other side, when it comes to computer programming, the Udwadia-Kalaba formulation appears to be preferred by the students once a reliable subroutine for computing the pseudoinverse matrix is made available. The Table 1 summarizes the opinions expressed by students when comparing the different dynamic formulations. Based on the Udwadia-Kalaba formulation a simple 2D multibody dynamics simulation program has been developed. The core of the software is written in Fortran. Preprocessing module has been developed in GTK+ embedded in the Ch scripting language [24]. Mechanism animation is obtained making use of the useful tools available in the Ch Mechanism Toolkit [25]. 10 D. de Falco, E. Pennestrı̀ and L. Vita 8 CONCLUSIONS Computing pseudoinverse by means of the Gram-Schmidt decomposition is usually less expensive than other methods based on the use of the SVD decompositions, Greville algorithm, or Householder reflections. However, these last methods require about twice as much arithmetic, but are more reliable and accurate near the limits of residual reduction. The loss of accuracy usually connected with the Gram-Schmidt method of computing the pseudoinverse requires the integration subroutine a considerable higher number of computational steps in order to maintain the error lower than a fixed limit. The other two methods discussed are very accurate and require less function evaluation calls during numerical integration. Table 2: TOL=10−10. Integration method: Adams-Moulton. Initial time step, ∆t = 0.66 · 10−2 s. No Baumgarte stabilization applied. Method Result SVD Fails1 Varga Fails2 Greville Correct up to 17.5 s LS-Householder Correct up to 16.5 s LS-Modified GS Correct up to 17.0 s 1 2 P i Ψ2i 3.7 · 10−7 9.9 · 10−2 4.5 · 10−5 2.3 · 10−4 1.1 · 10−4 NFUNC NSTEPS CPU 2,961,955 2,633,811 0.93 6,131,159 4,177,395 11.57 3,749,766 3,015,354 36.85 After 9.1 s of correct simulation, the integration subroutine is halted because of repeated error-test failures (see [9], p. 668). After 8.0 s of correct simulation, the integration subroutine is halted after failing to pass the error-test failures (see [9], p. 668). Table 3: TOL=10−8. Integration method: Adams-Moulton. Initial time step, ∆t = 0.66 · 10−2 s. No Baumgarte stabilization applied. Method P Result i Ψ2i SVD Fails1 3.4 · 10−6 Varga Fails2 9.9 · 10−2 Greville Correct up to 17.5 s 1.86 · 10−4 LS-Householder Correct up to 18.0 s 7.9 · 10−4 LS-Modified GS Correct up to 18.0 s 5.2 · 10−4 1 2 NFUNC NSTEPS CPU 1,648,055 1,422,689 1,641,751 1,422,297 1,644,341 1,424,864 0.53 1.28 6.03 After 12.07 s of correct simulation, the integration subroutine is halted because a fatal error 1 from DIVPAG (see [9], p. 668). After 8.06 s of correct simulation, the integration subroutine is halted because the precision loss of constraints equations. 11 D. de Falco, E. Pennestrı̀ and L. Vita Table 4: TOL=10−6. Integration method: Adams-Moulton. Initial time step, ∆t = 0.66 · 10−2 s. No Baumgarte stabilization applied. Method i SVD Correct up to 17.5 s Varga Fails1 Greville Correct up to 19.0 s LS-Householder Correct LS-Modified GS Correct 1 P Result Ψ2i 4.1 · 10−4 9.9 · 10−2 9.7 · 10−4 7.9 · 10−4 6.9 · 10−4 NFUNC NSTEPS CPU 1,033,381 847,294 1,017,558 1,422,689 1,641,751 1,422,297 1,013,795 839,020 0.28 0.37 0.76 3.40 After 7.3 s of correct simulation, the integration subroutine is halted because the precision loss of constraints equations. Table 5: TOL=10−10 . Integration method: Adams-Moulton. Initial time step, ∆t = 0.66 · 10−2 s. Baumgarte stabilization applied, α = 1000, β = 100. Method Result SVD Fails1 Varga Correct Greville Correct LS-Householder Correct LS-Modified GS Correct s 1 P i Ψ2i 4.5 · 10−5 2.1 · 10−7 3.3 · 10−7 4.1 · 10−7 NFUNC 78,491,398 74,043,990 73,707,588 73,910,268 NSTEPS CPU 47,074,599 11.4 44,053,874 18.3 43,948,609 51.5 44,515,071 241.0 After 1.17 s of correct simulation, the integration subroutine is halted because a fatal error 1 from DIVPAG (see [9], p. 668). Table 6: TOL=10−8 . Integration method: Adams-Moulton. Initial time step, ∆t = 0.66·10−2 s. Baumgarte stabilization applied, α = 1000, β = 100. Method Result SVD Fails1 Varga Fails2 Greville Correct LS-Householder Correct LS-Modified GS Correct 1 2 P i Ψ2i 4.90 · 10−7 3.5 · 10−7 4.9 · 10−7 NFUNC NSTEPS 80,109,155 55,569,995 80,197,458 54,468,468 80,112,706 54,604,452 CPU 19.4 55.0 290 After 0.17 s of correct simulation, the integration subroutine is halted because a fatal error 1 from DIVPAG (see [9], p. 668). After 14.97 s of correct simulation, the integration subroutine is halted because the precision loss of constraints equations. 12 D. de Falco, E. Pennestrı̀ and L. Vita Table 7: TOL=10−6 . Integration method: Adams-Moulton. Initial time step, ∆t = 0.66·10−2 s. Baumgarte stabilization applied, α = 1000, β = 100. Method Result SVD Fails1 Varga Fails2 Greville Correct LS-Householder Correct LS-Modified GS Correct 1 2 P i Ψ2i 1.9 · 10−6 1.14 · 10−6 1.06 · 10−6 NFUNC NSTEPS CPU 80,483,163 57,743,095 80,689,071 58,099,798 80,636,983 58,010,479 19.5 55 291 After 0.46 s of correct simulation, the integration subroutine is halted because a fatal error 1 from DIVPAG (see [9], p. 668). After 9.7 s of correct simulation, the integration subroutine is halted because the precision loss of constraints equations. 1e-06 6e-05 5e-05 Ψ2i 6e-07 P i 4e-05 P i Ψ2i 8e-07 3e-05 4e-07 2e-05 2e-07 1e-05 0 0 0 0.2 0.4 0.6 0.8 1 1.2 0 0.2 0.4 0.6 0.8 1 Time (s) b) Time (s) a) Figure 2: Violation of position constraints: a) Without stabilization b) With Baumgarte stabilization REFERENCES [1] Udwadia, F.E., Kalaba, R.E., A new perspective on constrained motion, Proceedings of the Royal Society of London, A439, 1992, pp.407-410. [2] Udwadia, F. and Kalaba, R. Analytical Dynamics a New Approach. Cambridge University Press, Cambridge, 1996. [3] Gauss, C.F., Über ein neues allgemaines Gründgesetz der Mechanik, Journal ür Reine und Angewandte Mathematik, vol.4, 1829, pp.232-235. [4] Anderson, E. et al. LAPACK - Users’s Guide. Society for Industrial and Applied Mathematics, Philadelphia, 1995. 13 1.2 D. de Falco, E. Pennestrı̀ and L. Vita [5] Lanczos, C., The Variational Principle of Mechanics, Dover Publications, Inc., Fourth Edition, 1970, pp.106-110. [6] Chetaev, N., Mécanique Rationelle, Éditions Mir, Moscou, 1993, pp.251-253. [7] Ara Arabyan, Fei Wu, An Improved Formulation for Constrained Mechanical Systems, Multibody System Dynamics, vol.2, 1998, pp.49-69. [8] Itiki, C., Neto, J.J., Complete automation of the generalized inverse method for constrained mechanical systems of particles, Applied Mathematics and Computation, 2004, pp.561-580. [9] IMSL Math library, vols.I, II, 1997. [10] Wampler, R.H., Modified Gram-Schmidt with iterative refinement, ACM TOMS 5, 1979, pp.494-499. Fortran code available at www.netlib.org [11] Wehage, R. and Haug, E. Generalized Coordinate Partitioning for Dimension Reduction in Analysis of Constrained Dynamic Systems. ASME Journal of Mechanical Design, 134:247–255, 1982. [12] Haug, E. and Yen, J. Generalized Coordinate Partitioning Methods for Numerical Integration of Differential-Algebraic Equations of Dynamics. In Haug, E. and Deyo, R., editors, Real-Time Integration Methods for Mechanical System Simulation, volume 69 of NATO ASI Series, pages 97–114. Springer-Verlag, 1991. [13] Kim, S. and Vanderploeg, M. QR Decomposition for State Space Representation of Constrained Mechanical Dynamic Systems. ASME Journal of Mechanisms, Transmissiona, and Automation in Design, 108:176–182, 1986. [14] Mani, N., Haug, E., and Atkinson, K. Singular Value Decomposition for Analysis of Mechanical System Dynamics. ASME Journal of Mechanisms, Transmissions, and Automation in Design, 107:82–87, 1985. [15] Singh, R. and Likins, P. Singular Value Decomposition for Constrained Mechanical Systems. ASME Journal of Applied Mechanics, 52:943–948, 1985. [16] de Jalon, J. and Bayo, E., Kinematic and Dynamic Simulation of Multibody Systems - The Real-Time Challenge, Springer Verlag, New York (1994). [17] Haug, E., Computer-Aided Kinematics and Dynamics of Mechanical Systems, Allyn and Bacon (1989). [18] Nikravesh, P., Computer-Aided Analysis of Mechanical Systems, Prentice-Hall (1988). [19] Varga, A., On computing generalized inverse systems using matrix pencil methods, Int. J. of Applied Mathematics and Computer Science, vol.11, 2001, n.5. [20] Greville, T.N.E., Some Applications of the pseudoinverse matrix, SIAM Rev., 1960, pp.1522. 14 D. de Falco, E. Pennestrı̀ and L. Vita [21] Udwadia, F.E., Kalaba, R.E., A Unified Approach for the Recursive Determination of Generalized Inverses, Computers and Mathematics with Applications, vol.37, 1999, pp.125130. [22] Golub, G.H., Van Loan, C.F., Matrix Computations, 3rd Edition, The John Hopkins University Press, 1996. [23] Stephenson, G., An Introduction to Matrices Sets and Groups for Science Students, Dover Publications Inc., 1965, pp.48-51. [24] The Ch Language Environment, User’s Guide, v.5.0, www.softintegration.com. [25] Ch Mechanism Toolkit, User’s Guide, v.1.5, www.softintegration.com. 15