Table of Contents Author Guidelines Submit a Manuscript
Journal of Applied Mathematics
Volume 2012, Article ID 638546, 12 pages
Research Article

A Direct Eigenanalysis of Multibody System in Equilibrium

1School of Mathematical Sciences, Beijing Normal University, Beijing 100875, China
2School of Aerospace, Tsinghua University, Beijing 100084, China

Received 4 September 2011; Accepted 2 November 2011

Academic Editor: Massimiliano Ferronato

Copyright © 2012 Cheng Yang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


This paper presents a direct eigenanalysis procedure for multibody system in equilibrium. The first kind Lagrange’s equation of the dynamics of multibody system is a set of differential algebraic equations, and the equations can be used to solve the equilibrium of the system. The vibration of the system about the equilibrium can be described by the linearization of the governing equation with the generalized coordinates and the multipliers as the perturbed variables. But the multiplier variables and the generalize coordinates are not in the same dimension. As a result, the system matrices in the perturbed vibration equations are badly conditioned, and a direct application of the mature eigensolvers does not guarantee a correct solution to the corresponding eigenvalue problem. This paper discusses the condition number of the problem and proposes a method for preconditioning the system matrices, then the corresponding eigenvalue problem of the multibody system about equilibrium can be smoothly solved with standard eigensolver such as ARPACK. In addition, a necessary frequency shift technology is also presented in the paper. The importance of matrix conditioning and the effectiveness of the presented method for preconditioning are demonstrated with numerical examples.

1. Introduction

Modal analysis of multibody system [1, 2] has important application in structure analysis and modal reduction method [36]. The eigenanalysis methods for constrained multibody system can be mainly divided into two categories, namely, the eigenanalysis method by transforming the equations of motion from DAEs to ODEs by selecting independent coordinates [79] and the direct eigenanalysis method [1012]. Compared with structural systems, the governing equation of constrained multibody system is differential algebraic equations, instead of the ordinary differential systems. In the reduction methods, the differential algebraic equations of multibody system are reduced into ordinary differential equations for independent general coordinates, through looking for a set of orthogonal basis of constraint Jacobian matrix to eliminate the constraints from the system governing equations. The reduction method is accurate, efficient, and stable provided that the degree of freedom of the system and the number of constraints are small. However, the reduction method is considered not suitable for large scale problems due to the spoilage of sparse matrices [13].

In the direct methods, the nonlinear differential and algebraic equations are simultaneously linearized about a given state of the system, that is, the algebraic constraint equations are kept during the eigenanalysis and the sparsity of system matrices is preserved. Depending on how the problem is actually formulated, the enlarged mass matrix may become structurally singular (eigenvalues associated to constraints become infinite); other approaches may result in a structurally singular enlarged stiffness matrix (eigenvalues associated to constraints become 0). Usually a frequency shift is applied, and structurally nonsingular enlarged matrices are resulted. But in most cases, matrices’ condition number of the frequency shifted system is very large and numerical singular due to the constraints and a corresponding precondition is needed. The root cause for numerical singularity can be understood by the fact that the dynamic equations and the Lagrange multipliers carry a dimension of force while the constraint equations and the general coordinate a dimension of length or angular. The literature [14, 15] suggested multiplying a parameter to the constraints’ corresponding rows and columns of system matrices as a preconditioning. The time step integration takes a parameter associate with time step and integration format to do the precondition. But the eigenanalysis about a static equilibrium has nothing to do with time integration, and we need to find an appropriate parameter.

For the direct eigenanalysis of constrained multibody system in static equilibrium, this paper proposes a procedure and a necessary preconditioning method of system matrices to guarantee correct eigensolution. This paper is organized as follows. Section 2 schematically describes the formulation of constrained multibody system based on Lagrange equation with multipliers and the direct eigenvalue problem about equilibrium. Section 3 presented a detailed condition number analysis of the system matrices and a preconditioning scaling technique. Section 4 introduces the selection of eigensolver and numerical examples.

2. Governing Equations of Multibody System and the Direct Eigenvalue Problem

The dynamic equations of a general multibody system are usually established with Lagrange’s equations with multipliers [16]. The system equations are 𝑑𝑑𝑡𝜕𝐿𝜕̇𝐱𝜕𝐿=𝜕𝐱𝜕𝚽𝜕𝐱𝑇𝝀+𝐟,𝚽(𝐱)=𝟎,(2.1) where 𝐱𝑛 is the vector of general coordinates, 𝝀𝑚 is the vector of Lagrange multipliers, 𝐿 is the Lagrange function associated with (̇𝐱,𝐱), 𝐟 is the vector of general nonconservative forces which are functions of (̇𝐱,𝐱,𝑡), Φ is the vector of constraint functions associated with (𝐱), and (𝜕Φ/𝜕𝐱)𝑇𝝀 is the vector of generalized constraint forces. Denotë̇𝑑𝐅(𝐱,𝐱,𝐱,𝑡)=𝑑𝑡𝜕𝐿𝜕̇𝐱+𝜕𝐿(̈̇𝜕𝐱+𝐟𝐱,𝐱,𝐱,𝑡),𝐆(𝐱)=𝜕𝚽(𝐱)𝜕𝐱𝑇.(2.2) Then (2.1) becomes ̈̇𝚽𝐅(𝐱,𝐱,𝐱,𝑡)+𝐆(𝐱)𝝀=𝟎,(𝐱)=𝟎.(2.3) The equilibrium state (𝐱𝟎,𝝀𝟎,𝑡0) satisfies𝐅𝟎,𝟎,𝐱𝟎𝐱,𝑡+𝐆𝟎𝝀𝟎𝚽𝐱=𝟎,𝟎=𝟎.(2.4) Assume small vibration about the equilibrium ̈̇(𝛿𝐱,𝛿𝐱,𝐱0+𝛿𝐱,𝝀0+𝛿𝝀,𝑡), then 𝐅𝛿̈̇𝐱,𝛿𝐱,𝐱𝟎𝐱+𝛿𝐱,𝑡+𝐆0𝝀+𝛿𝐱0𝚽𝐱+𝛿𝝀=𝟎,𝟎+𝛿𝐱=𝟎.(2.5) Carry out the Taylor expansion and reserve linear terms, we have 𝐌𝐱0𝛿̈𝐂𝐱𝐲+0𝛿̇𝐊𝐱𝐲+0,𝝀0𝛿𝐲=𝟎,(2.6) where 𝐌(𝐱0)=𝐌(𝐱0)𝟎𝟎𝟎, 𝐂(𝐱0)𝐂(𝐱0)𝟎𝟎𝟎, 𝐊(𝐱0,𝝀0)=𝐊(𝐱0,𝝀0)𝐆(𝐱0)𝐆𝑇(𝐱0)𝟎, 𝐌(𝐱0̈)=(𝜕𝐅/𝜕𝐱)(𝐱0), 𝐂(𝐱0̇)=(𝜕𝐅/𝜕𝐱)(𝐱0), 𝐊(𝐱0,𝝀0)=(𝜕(𝐅+𝐆𝝀)/𝜕𝐱)(𝐱0,𝝀0), 𝛿𝐲=𝛿𝐱𝛿𝝀.

Write the solution as 𝛿𝐲(𝑡)=𝑒𝑟𝑡𝐕, where 𝑟 is one of the eigenvalues, 𝐕 is the corresponding eigenvector, and the eigenvalue problem is 𝑟2𝐊𝐌+𝑟𝐂+𝐕=𝟎.(2.7) A frequency shift is applied to (2.7) by let 𝑟=𝑟𝛼, 𝛼+, and then 𝑟2𝐌+𝑟𝐂+𝐊𝐕=𝟎,(2.8) where 𝐌𝐌=, 𝐊=𝛼2𝐊𝐌+𝛼𝐂+, 𝐂𝐂=2𝛼𝐌+. The frequency shift makes us calculate the eigenvalues about 𝛼 and makes sure 𝐊’s reversibility (the probability of singular is zero when 𝛼0) as the system’s energy must be in the form of kinetic energy or deformation energy. The eigenvalue problem can be transformed into the state space form as the following: 𝑟𝐌𝐘=𝐊𝐘,(2.9) where 𝐘=𝐲𝑟𝐲, 𝐌=𝐂𝐌𝐈𝟎, 𝐊=𝐊𝟎𝟎𝐈.

3. Ill-Condition in the Matrices of the Perturbed Vibration Equations and the Preconditioning

Consider the shifted eigenvalue problem (2.8) in the undamped case 𝑟2𝐌+𝐊𝐕=𝟎.(3.1) It is more convenient to consider the standard eigenvalue problem than the original generalized one 𝐊1𝐌+𝑟𝟐𝐈𝐕=𝟎.(3.2) As the inverse of matrix 𝐊 is used, cases in large condition number of 𝐊 could make the direct solution of 𝐊1𝐌 result in erroneous eigensolution. As shown in Section 4, the eigenvalue problems (2.8) obtained by the linearization of the Lagrange equations with multipliers for practical multibody systems usually encounter serious condition number problem. The reason mainly lies in the different dimensions of the multiplier variables and constraint equations with that of generalized coordinates and governing equations.

Inspired by the eigenanalysis method through transforming the equations of motion from DAEs to ODEs by selecting independent coordinates, we divide the matrix into blocks and calculate its inverse by Gaussian elimination. Let 𝛼2𝐌+𝐊=𝑐𝐃 where 𝑐=𝛼2𝐌+𝐊. Without loss of generality, let the first 𝑚 (the number of constraints) rows of 𝐆 be linearly independent and be denoted as 𝐆1𝑚×𝑚. Divide the matrix as 𝐊=𝑐𝐃11𝑐𝐃12𝐆1𝑐𝐃21𝑐𝐃22𝐆2𝐆𝑇1𝐆𝑇2𝟎,(3.3) where 𝐃11𝑚×𝑚, 𝐃22(𝑛𝑚)×(𝑛𝑚), 𝐃12𝑚×(𝑛𝑚), 𝐃21(𝑛𝑚)×𝑚, 𝐆2(𝑛𝑚)×𝑚, 𝑛 is the number of freedom of general coordinates.

Then 𝐊1=𝑐1𝐆1𝑇𝐆𝑇2𝐏1𝐆2𝐆11𝐆1𝑇𝐆𝑇2𝐏1𝑐𝐆1𝑇+𝐆1𝑇𝐆𝑇2𝐏1𝐑𝐏1𝐆2𝐆11𝐏1𝑐𝐏1𝐑𝑐𝐆11+𝐐𝐏1𝐆2𝐆11𝑐𝐐𝐏1𝑐2𝐐𝐏1𝐑𝐍,(3.4) where 𝐍=𝐆11𝐃11𝐆1𝑇, 𝐏=𝐃22+𝐆2𝐍𝐆𝑇2𝐃21𝐆1𝑇𝐆𝑇2𝐆2𝐆11𝐃12, 𝐐=𝐆11𝐃12𝐍𝐆𝑇2, 𝐑=𝐃21𝐆1𝑇𝐆2𝐍. See appendix for detail of Gaussian elimination.

Assume 𝐍𝐏1𝐐𝐑=𝑂(1).(3.5) Then it shows 𝐊1=𝑂(𝑐) and by taking condition number cond()=1cond𝐊𝑐𝑂2.(3.6)

Equation (3.5) can be expected when the elements of eigenvector of the underlying ODE system have about the same magnitude, that is, the elastic parameters of the physical system are about the same magnitude in physical description, and this is often the case. The elastic parameters mean the Young’ modulus of flexible bodies, stiffness of spring forces, and so forth.

The condition number of the original system can be changed by scaling the multiplier variables and constraint equations. Let 𝝀=𝑘1𝝀 and multiply the constraint equations by 𝑘 in (2.1) we have 𝑑𝑑𝑡𝜕𝐿𝜕̇𝐱𝜕𝐿𝜕𝐱=𝑘𝜕𝚽𝜕𝐱𝑇𝝀+𝐟,𝑘𝚽(𝐱)=𝟎.(3.7) The corresponding block matrices are𝐊=𝑐𝐃11𝑐𝐃12𝑘𝐆1𝑐𝐃21𝑐𝐃22𝑘𝐆2𝑘𝐆𝑇1𝑘𝐆𝑇2𝟎,𝐊1=𝑐1𝐆1𝑇𝐆𝑇2𝐏1𝐆2𝐆11𝐆1𝑇𝐆𝑇2𝐏1𝑘1𝑐𝐆1𝑇+𝐆1𝑇𝐆𝑇2𝐏1𝐑𝐏1𝐆2𝐆11𝐏1𝑘1𝑐𝐏1𝐑𝑘1𝑐𝐆11+𝐐𝐏1𝐆2𝐆11𝑘1𝑐𝐐𝐏1𝑘2𝑐2𝐐𝐏1.𝐑𝐍(3.8) Then the condition number analysis shows cond𝐊𝑘𝑂max(𝑘,𝑐)×max1𝑐,𝑘2𝑐2.(3.9) Choosing 𝑘𝑐, the estimation of condition number decreases from about 𝑂(𝑐2) to about 𝑂(1).

The elements of matrix 𝐊 stand for the elastic coefficients of the system and the elements of 𝐆 and 𝐆𝑇work similarly as the stiffness of constraints. Taking 𝐊 as an enlarged stiffness matrix, the scaling is equivalent to changing the stiffness of constraints from 𝑂(1) to the same magnitude of physical stiffness.

Note that the condition cannot be improved by only scaling multipliers or by only scaling constraint equations. It is suggested that any different dimension in the procedure of modeling could lead to bad condition of numerical calculation.

The complete process for direct eigenanalysis of constrained deformable multibody system is thus as follows.

Step 1. Calculate the static equilibrium (𝐱0,𝝀0) from (2.4).

Step 2. Calculate 𝐌, 𝐂, 𝐊 at (𝐱0,𝝀0).

Step 3. Precondition 𝐊 with 𝑐 which is the maximum absolute element of 𝐊.

Step 4. Frequency shift with parameter 𝛼(110) with (2.8).

Step 5. Enlarge the matrix into state space as (2.9).

Step 6. Solve the eigenvalues and eigenvectors form (2.9) with eigensolver such as ARPACK.

Step 7. Add 𝛼to the eigenvalues and multiply the eigenvector components of Lagrange multipliers with 𝑐.

4. Selection of Eigensolver and Numerical Examples

The famous iteration method Implicitly Restarted Arnoldi Method (IRAM) [17, 18] is used in the numerical solution of eigenvalue problem (2.9) in the paper. Arnoldi method is one of the Krylov subspace methods [19] which calculate the largest few eigenvalues of a matrix. The presented numerical examples have verified the efficiency of IRAM in practical calculation. The famous parallel sparse linear system solver PARDISO [20, 21] is used to provide 𝐊1 of (2.9) in the paper. A cantilever and a four-bar linkage are calculated to illustrate the effect of the condition number of the system matrix and the effectiveness of the presented preconditioning method. It is shown that the numerical results are physically meaningless without scaling and the precondition presented in the paper is very effective.

4.1. A Cantilever

In this example, beam element and 3-dimensional solid element are separately used to model the cantilever in the multibody frame, as shown in Figure 1. The parameters of the cantilever are shown in Table 1. The models are computed with the presented direct eigenanalysis method in two cases, that is, with and without axial tension, and the results are compared with analytical results based on the mathematical physics, see Tables 2 and 3.

Table 1: Parameters for the cantilever.
Table 2: Eigenvalues comparison for cantilever (𝐹=0 Newton).
Table 3: Eigenvalues comparison for cantilever (𝐹=1𝑒10 Newton).
Figure 1: Two models of cantilever.

The governing equation and boundary conditions for the single-mode-free vibration of the cantilever under axial tension is as the following: 𝑑𝐸𝐼4𝑦𝑑𝑥4𝑑𝐹2𝑦𝑑𝑥2𝜌𝜔2𝑦||𝑦=0,0=0𝑑𝑦|||𝑑𝑥0𝑑=02𝑦𝑑𝑥2||||𝑙𝑑=03𝑦𝑑𝑥3||||𝑙=𝑎𝑑𝑦|||𝑑𝑥𝑙,(4.1) where 𝑥 is the material coordinate, 𝑦 is the lateral vibration displacement mode, 𝐸 is the Young’s modulus, 𝐼 is the sectional moment of inertia, 𝐹 is the axial tension, 𝜌 is the linear density, and 𝜔 is the circular frequency of vibration. 𝑙 is the length of the cantilever and 𝑎=𝐹/𝐸𝐼. With the conventional mathematical procedure, characteristic equation for frequency 𝜔 is 𝑟21𝑟cosh1𝑙+𝑟22𝑟cos2𝑙𝑟31𝑟sinh1𝑙𝑟32𝑟sin2𝑙+𝑎𝑟1𝑟sinh1𝑙𝑟2𝑟sin2𝑙=𝑟21𝑟sinh1𝑙+𝑟1𝑟2𝑟sin2𝑙𝑟31𝑟cosh1𝑙+𝑟1𝑟22𝑟cos2𝑙+𝑎𝑟1𝑟cosh1𝑙𝑟1𝑟cos2𝑙,(4.2) where 𝑏=𝜌𝜔2/𝐸𝐼, 𝑟1=(𝑎2+4𝑏+𝑎)/2, and 𝑟2=(𝑎2+4𝑏𝑎)/2.

It is shown in Tables 2 and 3 that the results without the scaling precondition are totally erroneous in both cases. In the two models, through computation 𝑐=1.39167𝑒+06 for 𝐹=0 and 𝑐=1.39252𝑒+06 for 𝐹=1𝑒+10. Precondition parameter 𝑘 is chosen the same as 𝑐 and the condition numbers in the two computed cases are decreased about 𝑂(𝑐2) times after the precondition. The results show that the presented precondition procedure could guarantee correct results for the direct eigenanalysis, and the unscaled models give out meaningless results by ARPACK as the condition numbers are too large for the nonphysical contribution of constraints. Note that model 2 is used to illustrate the precondition’s effectiveness for models including other flexible elements such as 3D solid elements, and the error of scaled model 2 to Euler beam analytical results is partly because the solid model does not agree with the plan cross section assumption of Euler-Bernoulli beam.

4.2. Four-Bar Linkage

The four-bar linkage is shown in Figure 2, and the parameters of the bar are shown in Table 4. The model is modeled with Cartesian coordinates, and the constraints in the model are all revolute joints. The beam element model is used in this example, as shown in Figure 3. The equilibrium under gravity is first calculated with dynamic relaxation method, [22] and the eigenvalues are calculated about the deformed configuration. The equilibrium configuration of the four-bar system is also shown in Figure 3; the computed results with and without the preconditioning scaling are shown in Table 5 compared with results of FEA software ABAQUS with the same model. Again the unscaled results in the case are all meaningless as they have nonphysical real parts. In this problem 𝑐=9.11047𝑒+13 and precondition parameter 𝑘 is chosen, the same value as 𝑐. The condition numbers are reduced about 𝑂(𝑐) times after the preconditioning, and correct results are obtained with the direct eigenanalysis method.

Table 4: Parameters for the linkage.
Table 5: Eigenvalues with or without scaling.
Figure 2: Four-bar linkage discrete model.
Figure 3: Modes of the four-bar linkage.

5. Summary and Conclusions

In this paper, the bad condition of the eigenvalue problem formulated from the direct linearization of the differential-algebraic equations of multibody system is reviewed and then discussed. The root cause may be attributed to the different dimensions of constraint multipliers and equations with physical system. Through analysis, a simple scale of the constraint multiplier variables and equations according to the stiffness of the elastic elements in the system can reduce the condition number of the system matrix to be inverted. After the preconditioning scale, eigensolvers such as ARPACK can be directly applied to the preconditioned system and yield correct result. The soundness of proposed approach and applicability are illustrated by comparing the results of a cantilever and a four-bar linkage with those obtained from either analysis or ABAQUS.


A. Gaussian elimination to get (3.4) from (3.3)



Cheng Yang and Zhengru Zhang were supported by National NSF of China under Grant 11071124.


  1. A. A. Shabana, “Flexible multibody dynamics: review of past and recent developments,” Multibody System Dynamics, vol. 1, no. 2, pp. 189–222, 1997. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  2. A. A. Shabana, Dynamics of Multibody Systems, Cambridge University Press, Cambridge, UK, 2005.
  3. W. H. A. Schilders, H. A. Van Der Vorst, and J. Rommes, Model Order Reduction: Theory, Research Aspects and Applications, vol. 13 of Mathematics in Industry, Springer-Verlag, Berlin, Germany, 2008. View at Publisher · View at Google Scholar
  4. M. Ripepi and P. Masarati, “Reduced order models using generalized eigenanalysis,” Proceedings of the Institution of Mechanical Engineers, Part K, vol. 225, no. 1, pp. 52–65, 2011. View at Publisher · View at Google Scholar
  5. K. Holm-Joslashrgensen and S. R.K. Nielsen, “A component mode synthesis algorithm for multibody dynamics of wind turbines,” Journal of Sound and Vibration, vol. 326, no. 3–5, pp. 753–767, 2009. View at Publisher · View at Google Scholar
  6. P. Koutsovasilis and M. Beitelschmidt, “Comparison of model reduction techniques for large mechanical systems: a study on an elastic rod,” Multibody System Dynamics, vol. 20, no. 2, pp. 111–128, 2008. View at Publisher · View at Google Scholar
  7. Y. Zhang, B. Wen, and S. Chen, “Eigenvalue problem of constrained flexible multibody systems,” Mechanics Research Communications, vol. 24, no. 1, pp. 11–16, 1997. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  8. M. Da Lio and R. Lot, “Analisi modale di sistemi multibody descritti in coordinate naturali,” pp. 6-C9, 1999.
  9. V. Cossalter, R. Lot, and F. Maggio, “The modal analysis of a motorcycle in straight running and on a curve,” Meccanica, vol. 39, no. 1, pp. 1–16, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  10. P. Masarati, “Direct eigenanalysis of constrained system dynamics,” Proceedings of the Institution of Mechanical Engineers, Part K, vol. 223, no. 4, pp. 335–342, 2010. View at Publisher · View at Google Scholar
  11. V. Boschi, R. De Salvo, P. Masarati, G. Quaranta, and V. Sannibale, “Seismic attenuation system synthesis by reduced order models from multibody analysis,” in Proceedings of the Multibody Dynamics, Milano, Italy, June 2007.
  12. J. L. Escalona and R. Chamorro, “Stability analysis of vehicles on circular motions using multibody dynamics,” Nonlinear Dynamics, vol. 53, no. 3, pp. 237–250, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  13. M. Géradin and A. Cardona, Flexible Multibody Dynamics: A Finite Element Approach, Wiley & Sons, New York, NY, USA, 2001.
  14. O. A. Bauchau, A. Epple, and C. L. Bottasso, “Scaling of constraints and augmented lagrangian formulations in multibody dynamics simulations,” Journal of Computational and Nonlinear Dynamics, vol. 4, no. 2, pp. 1–9, 2009. View at Publisher · View at Google Scholar
  15. C. L. Bottasso, D. Dopico, and L. Trainelli, “On the optimal scaling of index three DAEs in multibody dynamics,” Multibody System Dynamics, vol. 19, no. 1-2, pp. 3–20, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  16. C. Lanczos, The Variational Principles of Mechanics, Dover Publications, New York, NY, USA, 1970.
  17. R. B. Lehoucq and D. C. Sorensen, “Deflation techniques for an implicitly restarted Arnoldi iteration,” SIAM Journal on Matrix Analysis and Applications, vol. 17, no. 4, pp. 789–821, 1996. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  18. H. Saberi Najafi and E. Khaleghi, “A new restarting method in the Arnoldi algorithm for computing the eigenvalues of a nonsymmetric matrix,” Applied Mathematics and Computation, vol. 156, no. 1, pp. 59–71, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  19. Y. Saad, Numerical Methods for Large Eigenvalue Problems, Manchester University Press, Manchester, UK, 1992.
  20. O. Schenk and K. Gartner, “Solving unsymmetric sparse systems of linear equations with PARDISO,” in Computational Science-ICCS 2002, Part II, P. Sloot, C. J. K. Tan, J. J. Dongarra, and A. G. Hoekstra, Eds., vol. 2330, pp. 355–363, Springer, Berlin, Germany, 2002. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  21. O. Schenk and K. Gartner, “On fast factorization pivoting methods for sparse symmetric indefinite systems,” Electronic Transactions on Numerical Analysis, vol. 23, pp. 158–179, 2006. View at Google Scholar · View at Zentralblatt MATH
  22. W. J. Lewis, Tension Structures: Form and Behaviour, Thomas Telford Services, 2003.