Recursive and Residual Algorithms For The Efficient Numerical Integration of Multi-Body Systems
Recursive and Residual Algorithms For The Efficient Numerical Integration of Multi-Body Systems
Recursive and Residual Algorithms For The Efficient Numerical Integration of Multi-Body Systems
Keywords
Semi-recursive and recursive formulations; numerical integration; stiff differential
equations; residual algorithm
Abstract
In this paper a family of algorithms for multi-body dynamic simulation is introduced.
These algorithms use topological formulations to derive the differential equations that
govern the dynamic behavior of the mechanical system. The algorithms used are based
on the concept of velocity transformation, which defines a relationship between the
Cartesian velocities of the bodies and the velocities of the degrees of freedom. These
algorithms handle open-chain and closed-chain systems in a different way. In order to
handle both “non-stiff” and “stiff” problems, two numerical integration schemes have
been considered. “Non-stiff” systems are solved by means of an Adams-Bashforth-
Moulton PECE numerical integration scheme, which requires the computation of the
function derivatives. “Stiff” problems are integrated by using either BDF or NDF
methods, which require the computation of the residual of the equations of motion and,
optionally, the evaluation of the Jacobian matrix. The dynamic formalisms presented in
this paper have been developed so as to fit both numerical integration schemes. The
proposed algorithms have been implemented using an Object-Oriented Programming
approach that makes it possible to re-use the source code, keeping programs smaller,
cleaner and easier to maintain. Practical examples that illustrate the use and
performance of these implementations are included. These examples have also been
solved using a commercial multi-body simulation package and comparative results are
included. In most cases, the algorithms here presented outperform those implemented in
the commercial package leading to important savings in terms of total computation
times.
1 Introduction
In this paper, recursive algorithms that greatly improve the time required to carry out
simulations of multi-body systems (or MBS) are presented. Dynamic simulation of
multi-body systems is an important part of the technique known as Computer-Aided
Engineering (CAE). This is a rather complex problem that requires the solution of a set
of non-linear differential-algebraic equations, which in most practical cases are “stiff”.
According to Schwertassek [1], in order to solve this kind of problems, efficient
computer implementations must be developed, combining knowledge from three
different fields: multi-body dynamics, numerical analysis and computer science. Multi-
2
body dynamics makes it possible to select the most efficient formalisms for the
description of the differential equations of motion; numerical analysis makes for the
selection of the most appropriate numerical integration schemes; and finally, computer
science makes it possible to take full advantage of hardware and software tools
currently available to develop an efficient implementation.
The dynamic formulations used in this work are included within the group of
topological formulations which use the relative coordinates in the joints to model the
mechanical system. These coordinates, together with the Object-Oriented Programming,
present excellent characteristics of simplicity and efficiency for the implementation of
computer applications of dynamic simulation of MBS. The global formulations are still
used in many general-purpose programs (ADAMS, DADS, COMPAMM, …), in spite
of the fact that they are significantly less efficient than the topological formulations.
Probably this is mainly due to its greater simplicity –sometimes apparent– with respect
to the topological methods. In this paper it will be shown that, correctly implemented,
topological methods also compete in this aspect with global methods.
For open-chain systems, degrees of freedom coincide with the relative coordinates
defined at the kinematic joints. Equations of motion can be first written in terms of
Cartesian coordinates and transformed in a state-space representation using the velocity
transformation between the Cartesian velocities and the velocities of the relative
coordinates. For “non-stiff” systems, derivatives requested by the numerical integration
routine are computed following either fully recursive or semi-recursive formulations.
The fully recursive formulation yields to a O(N) implementation while the semi-
recursive formulation is O(N3). In the case of “stiff” systems the proposed algorithm
calculates the residual of the equations of motion following a fully recursive method.
An approximation of the Jacobian matrix is computed when requested by the numerical
integration routine. This approximation can be evaluated explicitly thus leading to
important savings in the number of function evaluations and more efficient computer
codes.
The problem of closed-chain systems is tackled in two steps. In the first step, loops
are opened by removing a kinematic joint or a rigid body. The resulting open-chain
system is solved using the method described in the previous paragraph. In the second
step, closure loop conditions are imposed using a second velocity transformation
between the velocities of the relative coordinates and the velocities of the degrees of
freedom. This velocity transformation (characterized by the matrix Rz) is obtained by
computing numerically at each time step a basis of the null-space of the constraint
equations that impose the closure loop conditions. As for open-chain systems, “non-
stiff” problems require the computation of the function derivatives, which is performed
following a semi-recursive formulation; “stiff” problems require the evaluation of the
residual of the equations of motion, which is carried out by using a fully recursive
approach. Optionally, an approximation of the Jacobian matrix is also evaluated.
Two different implicit integration routines have been used: the well-known DASSL
and the ode15s function. The ode15s integration function of Matlab (Shampine and
Reichelt [2]) has been adapted to perform the simulation based on the residual
algorithm. This modified version has been called ode15sR.
Practical examples that illustrate the proposed formulation (written in C++) are
included, comparing its numerical efficiency with the SIMPACK code (written in
Fortran). It must be noticed that, in all cases, the methods presented are almost twice as
3
fast as the SIMPACK code results. The residual algorithm has proved to be the most
efficient for “stiff” systems.
Another important characteristic of the proposed algorithms is that they do share
most of the functions required for the kinematic and dynamic analyses. This
characteristic permits the implementation of a computer program of reduced size that
incorporates all the algorithms presented in this paper. With the use of the Object-
Oriented Programming paradigm in the implementation, computer codes are smaller,
cleaner and easier to maintain and extend.
The so-called global formulations solve all kinds of multi-body systems in the same
way. This makes easier to develop general-purpose programs, although they are
significantly less efficient than the topological formulations.
On the other hand, topological methods try to take advantage of the system structure.
They emerged in the 1960s and 1970s, especially for space and robotics applications.
These formulations arose from the study of the dynamics of serial robots, and were
extended to other more general open-chain and tree-type multi-body systems, and to
those with a closed-chain configuration. In this paper, special attention is paid to the
Articulated Inertia method (Featherstone [3]), following the notation introduced by Bae
and Haug [4]-[5] and re-formulated by Jiménez [6]. The following sections will briefly
describe several dynamic formalisms, including among them the Articulated Inertia
method. The equations of motion are expressed in the form,
y = f (y, t ) (1)
In this approach, the integration routine requires the derivative of the state vector y.
As will be shown later, this is neither the only nor the most efficient way of giving to
the integration routine the information of the differential equations of the mechanical
system.
The Articulated Inertia method (Featherstone [3]) is a fully recursive formulation that
allows for the computing of the accelerations of a mechanical system without the need
to solve any system of linear equations; consequently this is an O(N) formulation.
In this paper, the Articulated Inertia method is introduced using a set of Cartesian
coordinates composed of the Cartesian velocity of the material point that is located
instantaneously at the origin of the inertial reference frame ( si ) and the angular velocity
of the body ( ωi ). The Cartesian velocities (denoted by Z ) of two consecutives bodies
are given by,
s s
Zi −1 = i −1 Zi = i (2)
ωi −1 ωi
where z i are the velocities of the relative coordinates at the kinematic joint that links
both bodies. The matrix bi is a term that depends on the type of joint considered and is
a function of the joint position. Taking the time derivative of equation (3), a relationship
between Cartesian accelerations of two consecutive bodies is obtained,
=Z
Z + b
i i -1 i z i + di (4)
where di depends again on the joint type and is a function of the joint position and
velocity. Kinematic problems (position problem and velocity and acceleration analyses)
are solved by means of a recursive procedure.
Applying the Principle of Virtual Power, the equations of motion may be written as
follows,
N
∑ Z ( M Z
i=1
*T
i i i - Qi ) = 0 (5)
The mass matrix M i and the forces Qi are defined respectively by,
m I − mi g i
M i = DTi M i Di = i 3 (6)
mi g i J i − mi g i g i
fi − ωi × ( ωi × mi g i )
Qi = DTi ( Qi − M i ei ) = (7)
ni − ωi × J i ωi + g i × ( fi − ωi × ( ωi × mi g i ) )
I −g i
Di = 3 (8)
0 I 3
ω × ( ωi × g i )
ei = i (9)
0
where fi y ni are the forces and torques acting on the center of gravity of the i-th
element, respectively; Di is the matrix that transform the reference point from the centre
of gravity to the origin of the global reference frame; ei is the vector of velocity
dependent accelerations; mi is the mass of the body and J i is the inertia tensor in the
inertial reference frame, which can be expressed in terms of the well-known tensor J i
defined in the body local reference frame as,
J i = A i J i ATi (10)
∑ Z ( M Z
i =1
*T
i i i − Qi ) +Z*NT−1 ( M N −1Z
N −1 − Q N −1 ) + Z N ( M N Z N − Q N ) = 0
*T (11)
5
T for the expressions (3) and (4) particularized for the N-the
Substituting Z*NT y Z N
body, and after some algebraic manipulations, the following result may be easily
obtained,
N −2
∑ Z ( M Z*T
i i i − Qi ) + Z*NT−1 ( M N −1Z
N −1 − Q N −1 ) +
i =1 (12)
+ (Z *T
N −1 + z b ) M N ( Z
*
N
T
N
N −1 + b N z N + d N ) − Q N = 0
z N = ( bTN M N b N ) (b (Q − M N d N ) − bTN M N Z )
-1
T (13)
N N N -1
K N ≡ I 6 − M N b N ( bTN M N b N ) bTN
−1
(14)
ˆ
M (
N −1 ≡ M N −1 + M N − M N b N ( b N M N b N ) b N M N = M N −1 + K N M N
T −1 T
) (15)
N −1 ≡ Q N −1 + ( Q N − M N d N ) − M N b N ( b N M N b N ) b N ( Q N − M N d N ) =
ˆ
Q T −1 T
(16)
= Q N -1 + K N ( Q N − M N d N )
( )
N -2
∑ Z ( MZ
i =1
*T
i i − Qi ) + Z*NT-1 M
ˆ Z ˆ
N −1 N −1 − Q N −1 = 0 (17)
Equation (17) is similar to equation (5). However, two differences can be found: it is
written only for (N-1) bodies and the mass matrix and forces corresponding to the body
(N-1) have been modified by equations (15) and (16). The matrix M ˆ is the
N -1
( )
−1
ˆ b bT M
K i ≡ I6 − M ˆ b bTi (18)
i i i i i
ˆ ≡ M +K M
M ˆ (19)
i -1 i -1 i i
ˆ ≡Q +K Q
Q i −1 i −1 i
ˆ −M
i
ˆ d
i i ( ) (20)
6
( ) (b (Qˆ - Mˆ d ) - b Mˆ Z )
-1
ˆ b
z i = bTi M T T
(21)
i i i i i i i i i -1
This procedure can be followed recursively until the root element (i.e., the ground) is
reached. At this point, the value of the relative acceleration in the first kinematic joint
can be easily computed. A new recursive calculation process is followed to compute the
accelerations of the relative coordinates, using equations (4) and (21). This recursive
procedure starts from the root and follows upwards until the “leaves” of all the
“branches” of the kinematic “tree” are reached. From this description it can be seen that
the Articulated Inertia method solves the dynamic problem following three recursive
steps: solution of the kinematic problems, calculation of the “Articulated terms” (mass
matrices and force vectors) and computation of the accelerations of the relative
coordinates.
2.1.2 Semi-recursive formulations
An alternative to the global formulations and the fully recursive methods is the use of
semi-recursive formulations. For open-chain systems, these methods are based on the
use of two different coordinate systems (Cartesian and relatives) and the analysis is
carried out following a “two stages” procedure. The velocity relationship between the
two coordinate systems is expressed by the matrix R, and is consequently referred to as
velocity transformation method. This matrix constitutes a basis for the null-space of the
Jacobian matrix corresponding to the kinematic constraint equations. This kind of
formulations leads to the equation (see García de Jalón and Bayo [7]):
R T MRz = R T ( Q − MRz
) (22)
where M represents the mass matrix and Q the external forces vector of the MBS.
For open-chain systems, the matrix R that relates Cartesian velocities with the
velocities of the relative coordinates can be computed directly with very few arithmetic
operations avoiding the formation and factorization of the Jacobian matrix. This has
been used by many authors, as Jerkovsky [8], Kim & Vanderploeg [9] and Nikravesh &
Affifi [10]. More information can be found in the works of García de Jalón and Bayo
[7].
As in the Articulated Inertia method (Section 2.1.1), the Cartesian velocity of the
material point that is located instantaneously at the origin of the inertial reference frame
and the angular velocity of the body will be considered as Cartesian coordinates Z for
each body in the MBS, as shown in Figure 1. This set of coordinates has also been used
by Jiménez [6]. The work of Kim [11] showed a similar semi-recursive formulation.
The key point of the proposed semi-recursive formulation can be found on the
special shape of the matrix R, denoted as R when using the Cartesian velocities Z. This
velocity transformation can be written as,
Z = Rz (23)
7
3 5
z 3 z 5
2 4 z 6 6
z 2 z 4
1
z 1
In equation (23), the matrix R can be written in terms of sub-matrices or blocks. The
block R ij represents the Cartesian velocities of j-th body due to a unit velocity in the
relative coordinates defined in the i-th joint, taking the remaining independent velocities
as equal to zero. For the particular multi-body system depicted in Figure 1, this matrix is
expressed as,
R11 0 0 0 0 0
1
R 2 R 22 0 0 0 0
R13 R 2
R 3
0 0 0
R= 1 3 3
(24)
R 4 0 0 R 44 0 0
R1 0 0 R 54 R 55 0
15
R 6 0 0 R 64 0 R 66
It can be easily demonstrated that the selected Cartesian coordinates cause the effect
of a unit velocity in a relative coordinate to produce exactly the same numerical value
for the Cartesian coordinates of any body that is located upwards in the kinematic chain.
Consequently equation (24) leads to,
I 0 0 0 0 0 R1 0 0 0 0 0
I
I 0 0 0 0 0 R2 0 0 0 0
I I I 0 0 0 0 0 R3 0 0 0
R= 4 = TH (25)
I 0 0 I 0 0 0 0 0 R 0 0
I 0 0 I I 0 0 0 0 0 R5 0
I 0 0 I 0 I 0 0 0 0 0 R 6
or in an abbreviate form,
Z = R z = TH z (26)
the transformation matrix equals the number of kinematic joints, thus dramatically
reducing the number of blocks required in equation (24). It can also be demonstrated
that these matrices R i correspond to the bi terms of the Articulated Inertia method.
Using the Principle of Virtual Power, the dynamic equations can be written as,
Z*T ( MZ
− Q) = 0 (27)
where M y Q are diagonal block matrices defined in the equations (6) and (7). Taking
the time derivatives of equation (26), the relation between Cartesian and relative
accelerations can be obtained as,
= R
Z z+R z = TH z
z + TH (28)
Introducing this expression into equation (27), the set of dynamic equations is written
as,
( )
HT TT MT Hz = HT TT ( Q − MTHz(
) ) (29)
In the particular case of the MBS depicted in Figure 1, the transformed mass matrix
T
T MT has the following expression,
M1Σ M Σ2 M 3Σ M Σ4 M 5Σ M Σ6
M Σ2 M 3Σ 0 0 0
M Σ
0 0 0
TT MT = 3
Σ Σ
= MΣ (30)
M 4 M 5Σ M6
sim. M 5Σ 0
M Σ6
where,
M Σ6 = M 6 , M 5Σ = M 5 , M 3Σ = M 3 , M Σ2 = M 2 + M 3Σ
(31)
M Σ4 = M 4 + M 5Σ + M Σ6 , M1Σ = M1 + M Σ2 + M Σ4
The expressions (31) illustrate the concept of “accumulated inertia”: the mass matrix
corresponding to each body is a recursive addition of the mass matrices of the bodies
that are found upwards in the kinematic chain. This result is a consequence of the
selected set of Cartesian coordinates.
A similar result is obtained for the “accumulated forces” Q ∑ ,
9
Q 6∑ = Q 6 − M 6 ( THz
)
6
Q1∑
∑ Q = Q5 − M 5 ( THz
∑
5
)
Q 2
5
Q ∑ Q 4∑ = Q 4 − M 4 ( THz
) + Q∑ + Q∑
5 6
Q ∑ = TT Q = 3∑
4
where, (32)
Q 4 Q3∑ = Q3 − M 3 ( THz
)
3
Q5∑
∑ Q = Q 2 − M 2 ( THz
∑
2
) + Q∑
3
Q 6
2
Q = Q1 − M1 ( THz
∑
1
) + Q∑ + Q∑
2 4
1
method (Section 2.1.1). Consequently the product ( THz ) is the “accumulation” of the
j
terms d i located in the path from body j to the end of the kinematic chain.
It is shown in equation (26) that the position and velocity problems are solved by
applying a forward recursive algorithm which calculates Cartesian positions and
velocities from the known values of the relative coordinates. There is a second recursive
step that computes backwards a sort of accumulated mass matrices and forces. This
second recursive step is represented in equation (29) by the matrix products on which
the “accessibility matrix” T is involved.
Introducing the new notation, the dynamic equations can be re-written as,
H T M ∑ Hz = H T Q ∑ (33)
It must be noticed that the matrix H is a block diagonal matrix, that the matrix M Σ is
symmetric (with a topological structure defined by the accessibility matrix T) and that
therefore this product can be made in a very efficient way.
The velocity transformation presented in the equation (26) can be computed
following a joint-by-joint procedure in an independent way for each kinematic joint.
The final step is made solving a system of equations in a global fashion, which is the
reason why the method is O(N3) order. Jiménez [6] demonstrated that this formulation is
more efficient than other recursive formulations based on body-by-body calculations,
like the one presented by Avello et al. [14].
This is the most difficult and interesting case. In Section 2.1, two of the main families of
topological methods for the dynamic analysis of mechanisms of open-chain systems
have been described. When the system under consideration presents closed kinematic
chains, the relative coordinates in the joints are not independent. Whereas the global
methods deal with open and closed-chain mechanisms in the same way, the topological
formulations need to carry out some additional processes. The topological methods must
eliminate in a first step some kinematic joints in order to transform the mechanism into
an open-chain system. In a second step, the constraint equations corresponding to the
removed kinematic pairs are imposed (closed-chain conditions). This second stage can
be introduced by means of any of the methods used by the global formulations to
10
introduce the constraint equations. The topological methods are still much more
efficient than the global methods, because they impose the restrictions in small size
systems.
2.2.1 Fully recursive formulations
It is possible to extend the full-recursive method described in Section 2.1.1 to the case
of closed-chain mechanisms. In a first step, a joint (sometimes an element) is eliminated
in each loop. The closure condition can be introduced in two different ways: by means
of the Lagrange’s Multipliers (Bae and Haug [5]) or by means of the Penalty
formulation (Jiménez [6]).
Nevertheless, the full-recursive formulations for closed-chain systems lead to very
complicated expressions or with numerical problems, which are very difficult to solve.
These are, among others, the reasons why this family of methods is hardly used.
2.2.2 Semi-recursive formulations: Matrix Rz method
(
R T MRz + ΦTz λ = R T Q − MRz ) (34)
In the system of differential equations (34), the relative coordinates are no longer
independent. The size of the vector of Lagrange’s Multipliers λ is equal to the number
of independent constraint equations which relate the relative coordinates z. This
situation is similar to the one that appears when using global formulations and the same
possibilities can be used to solve the problem: that of solving a system of algebraic-
differential equations of index 3 (or index 1 with some form of stabilization); that of
following a Penalty approach; or that of applying a new velocity transformation to
obtain to a system of ODEs with independent coordinates.
In the equation (34) the Lagrange’s Multipliers can be eliminated with a second
velocity transformation that can be computed in a very efficient way solving a system of
linear equations which is much smaller than the original one. It is possible to select z i
as a subgroup of the relative coordinates z that constitute a system of independent
coordinates in the closed-chain system. These independent coordinates can be chosen
from the factorization of the Jacobian or tangent matrix Φ z according to the
methodology presented by several authors in the bibliography.
The constraint equations corresponding to a generalized joint are expressed in terms
of the Cartesian coordinates of the points and the components of the unit vectors defined
at the joint. The equations can be symbolically written as,
Φ ( r j , u j , rk , u k ) = 0 (35)
The calculation of the Jacobian matrix Φ z is made by deriving in a first step the
constraint equations with respect to the coordinates of the points and the components of
the unit vectors, in which terms the restrictions are expressed; in a second step, the
derivatives of these coordinates with respect to the relative coordinates are evaluated.
Thus, the Jacobian matrix leads to,
11
The partial derivatives of equation (36) represent the velocities induced in the points
and vectors located in the restrictions when giving a unit value to the velocity of a
relative coordinate and zero to the others. Once again, the way of evaluating these
derivatives depends on the type of kinematic joint considered and has an expression
very close to the terms R i or b i of Section 2.1. The constraint equations for velocity
and acceleration can be obtained by differentiating the constraint equations twice with
respect to time,
Φ z z = 0 (37)
Φ z z
z = −Φ (38)
z
Introducing the equation (41) into the equation (34), pre-multiplying by R Tz and
considering the expression (40), the equations of motion are given by,
) − R T M RR
z i = R Tz (R T (Q − M Rz
R Tz R T M RR z z i ) (42)
z
The evaluation of the right hand side of the equation (42) may appear to be costly
and somewhat difficult. From the equation that relates the Cartesian velocities to the
independent coordinates,
Z = Rz = RR z z i (43)
Using the equation (44), it is possible to rewrite (42) in the following way,
d ( RR z ) i
z i = R Tz RT Q − M
R Tz RT M RR z z (45)
dt
d ( RR z )
z i (46)
dt
can be evaluated, from the equation (44), like the vector of Cartesian accelerations
which results from giving a zero value to the independent relative accelerations.
Due to the excellent results obtained by using this method, it is worth at this point to
take the time to describe the calculation algorithm in detail. The chosen variant of the
method solves the stabilization of the constraint equations by means of a system of
algebraic-differential equations of index 2, integrating the dependent velocities and the
{
independent accelerations y T = z T }
z iT . The following algorithm makes it possible to
calculate the derivatives y from the vector state y .
{
1. Start at a time t, when position and velocity are known y T = zT }
z iT .
( )
6. Obtain d ( RR z ) dt z i from equation (44) as the Cartesian open-chain system
accelerations Z evaluated recursively with the relative accelerations evaluated in
the previous step.
7. Form the accumulated inertia mass matrices M ∑j and the accumulated force
vectors Q ∑j for each joint in the mechanical system.
(
10. Form the last matrix products of equation (45), that is, R Tz HT M ∑ H R z and )
R T
z (H Q ) .
T ∑
11. Obtain the independent accelerations zi , which together with the relative
velocities z form the vector y .
13
n M Φ1T
q Φq2T R
m1 Φ1q 0 0
m2 Φq2 0 0
N=n-m1-m2
N
N1 R T MR Sym.
N1 N
Rz
m2 Φq2 R 0 R Tz R T MRR z
This Section follows the developments of Rodríguez [16]. The calculation of the
derivative of the state vector is neither the only way (nor the most efficient) of solving
the equations of motion using a numerical integration routine. Implicit integrators do
not require the computation of the function derivative y . This is because the algorithm
used by these integrators to calculate a new value of the function uses the Newton-
Raphson method to adjust the value of the state vector y , while at the same time
making a numerical estimation of the derivative y . For example, using the implicit
form of the differential equations and defining the corresponding residual rs as,
rs ≡ F(y , y , t ) = 0 (47)
the Newton-Raphson method converts the residual to zero by substituting the derivative
y in terms of y , by means of the following iterative process,
14
∂F (y i , y (y i ), t ) i +1 i
∂y
( )
y − y = −F i (y i , y (y i ), t ) (48)
( M Z − Q ) = 0
N
∑Z
i=1
*T
i i i i (49)
or
(
− Q = 0
Z*T MZ ) (50)
where in order to simplify the expressions, the terms Q incorporate terms c (forces due
to the centrifugal and Coriolis acceleration of the elements) as in the semi-recursive
formulation presented in Section 2.1.2.
In this Section, it will be assumed that the kinematic chains are composed of
kinematic joints with a single degree of freedom of relative motion, although this point
can be easily generalized.
The equation (49) can be re-written in the following way,
( M Z − Q ) + Z ( M )
N −2
*T − Q = 0
∑Z
i =1
*T
i i i i
*T
N −1 N −1 Z N −1 − Q N −1 + Z N
MNZ N N
(51)
By replacing Z*TN with the expression based on the virtual velocity of the previous
( M Z − Q ) + Z ( M )
N −2
∑Z *T
i i i i
*T
N −1 N −1
Z N −1 − Q N −1 + Z N −1 + z N b N
*T
(
− Q = 0 (52)
* T M N Z N N
)
i =1
In this case, unlike the Articulated Inertia formulation, Z is not replaced for its
N
by the numerical integration routine. By grouping the terms pre-multiplied by the virtual
velocity z *N , the following expression is obtained,
( M Z − Q ) + Z ( M Z )
N −2
∑Z
i =1
*T
i i i i
*T
N −1 N −1 N −1 − Q N −1 +
(53)
+Z *T
N −1 ( M Z − Q ) + z b ( M Z
N N N
*
N
T
N N N )
− QN = 0
bTN M N Z N(
− Q = 0
N ) (54)
Because the values of the accelerations given by the integration routine are not exact,
this equation will not be identically null, but that it will have a certain residual value
which will be denominated as rs N and will correspond to the residual of equation N.
For a generalized joint with more than one relative degree of freedom, this residual will
be a vector of size equal to the number of relative degrees of freedom, as follows,
− Q
rs N = bTN M N Z N N ( ) (55)
The residual rs N has a clear physical meaning and represents the projection of the
sum of inertia and external forces acting on the relative degrees of freedom of the input
joint of the element. The sum of these projection forces, which if calculated with the
real accelerations would have to be in balance, is not zero because approximated values
of the accelerations have been used.
* in (53) it is obtained,
Grouping the terms multiplied by Z N −1
( M Z − Q ) + Z ( M )
N −2
∑Z
i =1
*T
i i i i
*T
N −1 N −1 Z N −1 − Q N −1 + M N Z N − Q N = 0 (56)
( M Z − Q ) + Z ( M )
N −2
∑Z
N −1 − Q N −1 + RS N = 0
*T *T
i i i i N −1 N −1 Z (58)
i =1
16
( M Z − Q ) +Z ( M )
N −3
∑Z
i =1
*T
i i i i
*T
N −2 N −2 Z N − 2 − Q N − 2 + M N −1Z N −1 − Q N −1 + RS N +
(59)
+ z b ( M
*T T
N −1 N −1 N −1
Z N −1 )
− Q N −1 + RS N = 0
Once again, since z *NT−1 is independent of the remaining virtual velocities, the bracket
multiplying it would have to be zero, so that the identity for any value of this virtual
velocity must be verified. From this equation, the following residual term is obtained,
(
rs N −1 = bTN −1 M N −1Z N −1 − Q N −1 + RS N ) (60)
where RS N represents the external and inertia forces that acting on the N-th element are
transmitted to the N-1 element.
Similarly, it is possible to define the value of RS N −1 in the form,
RS N −1 ≡ M N −1Z N −1 − Q N −1 + RS N (61)
which represents the accumulated external and inertia forces acting on the two last
elements. The equation (61) indicates that the residual forces acting on different
elements are added directly, with no type of transformation or projection. This is a
consequence of the Cartesian coordinates used by all the bodies in the mechanism (the
material point of each body instantaneously located in the origin of the reference frame).
From the practical point of view equation (61) can be replaced by the corresponding
equation involving the centre of gravity as reference point,
( − Q
RS N −1 = D N −1 M N −1Y N −1 N −1 + RS N ) (62)
( − Q + RS
rsi = bTi M i Z i i i +1 ) (63)
− Q + RS
RSi = M i Z (64)
i i i +1
or, alternatively,
( ( − Q + RS
rsi = bTi Di M i Y i i i +1 ) ) (65)
i (
− Q + RS
RSi = Dt M i Y i i+1 ) (66)
17
The expression (64) represents the accumulation of the external and inertia forces in
a generic body of the MBS. The expression (63) represents the residual motor effort that
should be applied in order to balance the effect of these elements.
This process of condensation or accumulation of errors or residuals is completed
when the base body of the kinematic chain is reached. If the base body is a free solid
with six degrees of freedom, the values of the Cartesian accelerations Z * will be
1
independent and the final residual will correspond to equation,
− Q + RS
rs1 = RS1 = M1Z (67)
1 1 2
If the base body is joined to the ground by means of a kinematic joint which allows
less than six degrees of freedom, the resulting residual equation will be,
(
− Q + RS
rs1 = b1T M1Z 1 1 2 ) (68)
rs = (
R T M ⋅ Z
− Q
) (69)
z − z
(4).
4. The terms M i y Qi are calculated.
5. Calculation of the residual term corresponding to the last body as
( − Q .
rs N = bTN M N Z N N )
− Q , which will be accumulated
6. Calculation of the term RS N as RS N = M N Z N N
The algorithm described considers a MBS with a single kinematic “branch”. In the
case of a MBS with several branches, the formulation can also be applied. It will only
be necessary to consider the sum of the different terms RSi which arrive at the junction
element under consideration. This algorithm is very simple to program, especially when
using Object-Oriented Programming techniques.
The method described in the previous Section for open-chain MBS can be extended to
the case of systems that present closed kinematic chains. The dynamic analysis is made
by transforming the system into an open-chain configuration by opening its loops
(removing certain kinematic joints or eliminating some of the solids of the kinematic
chains, such as rod elements). The resulting system is formed by an open-chain system
plus the constraints associated with the loops. The method used to impose the loop
constraint equations determines the algorithm to be used for the calculation of the
residual. In this paper, the velocity transformation method will be used. In this paper,
this method is called the “Method of the R z matrix”.
By means of two consecutive velocity transformations, it is possible to express the
vector of Cartesian velocities Z in terms of the independent relative velocities z i
(subgroup of the relative velocities) in the form indicated by the equation (see equation
(43)),
Z = RR z z i (70)
In the equation (71) the virtual velocities z i* are independent and can therefore be
eliminated, obtaining the following dynamic equations,
R Tz RT ( M ⋅ Z
− Q) = 0 (72)
8. Once the vector z is known, the Cartesian accelerations Z are evaluated in a
i
( )
d
In this algorithm, steps 3, 4 and 7 make use the same matrix ΦTz Φ z which is
factorized only once in the calculation process. Because a least-squares formulation has
been used in order to consider redundant constraint equations, this matrix is symmetric
and positive definite. Efficient solvers for this kind of matrices are used. In addition, the
matrix will be very small because the number of dependent relative coordinates is small
in most practical cases.
This is therefore quite a efficient process, especially if it is compared to the global
methods where the number of constraint equations is considerably greater.
It is widely accepted that it is necessary to use implicit integration methods with MBS
regarded as “stiff” systems. In previous Sections (3.1 and 3.2) efficient methods to
evaluate the residual of the equations of motion have been introduced. Nevertheless, it
is difficult to find general-purpose integrators that can integrate on the basis of the
calculation of the residual. An exception is found in the MBSSIM library (von Schwerin
and Winckler [17]), to which it has nevertheless been impossible to have access. Of
course, none of the integrators distributed in Matlab 5.3 allows this possibility.
20
Whereas the DASSL code must calculate the Jacobian matrix necessary for the
iterative Newton-Raphson process without knowing the structure of the equations, the
functions that use the expression (73) can pose the calculation of the Jacobian matrix in
a more efficient way.
The function ode15s found in Matlab can use BDF (Backward Differentiation
Formulas) or NDF (Numerical Differentiation Formulas), as described by Shampine
and Reichelt [2]. It is a cuasi-constant step and order up to 5 implementation. The
original version of the function ode15s does not include the possibility of carrying out
the integration process on the basis of the residual evaluation (knowing y e y values).
The core of this function has been modified so that the term to evaluate is r = F(y, y , t ) ,
and therefore, the residual algorithm presented in Sections 3.1 and 3.2 can be used. This
modified version of the integrator has been denominated as ode15sR. The new version
is faster than the original because it is not necessary to form the matrix of the
mechanical system nor to carry out an inversion to evaluate the acceleration.
In this Section an approximated expression for the analytical evaluation of the Jacobian
matrix is presented and used with the implicit DASSL and ode15s routines. Both the
ode15sR function and the ode15s function share the same Jacobian matrix.
3.4.1 The DASSL Function
In order that the user can compute the matrices in equations (74) and (75), DASSL
provides approximate values for the state vector y and the derivative y .
If the MBS has an open-chain topology, the residual can be formulated in terms of
the relative coordinates in the form,
rs1 z − z
rs = = T
( )
(76)
rs 2 R MRz − R Q − MRz
T
21
In order to calculate the Jacobian matrix shown in equation (75) and considering that
the coefficients c j are provided by DASSL, it is necessary to calculate the matrices,
where the derivatives that appear in these expressions can be estimated as,
∂rs1 ∂rs1 ∂rs1 ∂rs1
=0 = −I =I =0
∂z ∂z ∂z z
∂
(78)
∂rs 2 ∂Q ∂rs 2 ∂Q ∂rs 2 ∂rs 2
−RT −RT =0 RT MR
∂z ∂z ∂z ∂z
∂z
∂
z
The values shown correspond to the elimination of all the partial derivatives with
respect to the matrix R T . These terms can be eliminated because its value is normally
small in relation to the terms of forces, especially in the “stiff” cases. For “non-stiff”
cases, this approach is still valid because its influence on the matrix is not important.
The Jacobian matrix is (2Nz×2Nz), where Nz is the number of relative coordinates in
the joints. It is formed by four sub-blocks (Nz×Nz) which, once the equations (78) are
introduced, can be written as,
c j I Nz −I N z
J DASSL = T ∂Q T ∂Q
(79)
−R −R + c j R MR
T
∂z ∂z
The coefficient c j is related to the inversion of the integration step and usually has a
value of between 103 and 105. This value is provided by DASSL and varies throughout
the integration process.
For closed-chain systems, the Jacobian matrix using the method of the R z matrix
can be written as,
c j I Nz −R z
J DASSL = T T ∂Q ∂Q
(80)
−R z R − RTz R T R z + c j R Tz R T MRR z
∂z ∂z
The Jacobian matrix used by the ode15s integration scheme can be expressed as,
∂rs ∂f (y n +1 , tn +1 )
= M (y n +1 , tn +1 ) − hinvk ⋅ (81)
∂y ∂y
where the variation of the matrix M (y n +1 , tn +1 ) with respect to the state vector y n +1 has
been assumed to be of little significance. The ode15s code evaluates the term
M (y n +1 , tn +1 ) and the Jacobian matrix ∂f (y , t ) / ∂y separately. From the equation (76)
22
and using the same approach as in the DASSL case, these two matrices take on the
values,
I 0
M (y, t ) = N z (82)
0 R MR
T
0
I Nz
∂f (y, t )
= T ∂Q ∂Q (83)
∂y R RT
∂z ∂z
By introducing the results (82) and (83) into expression (81), the same result used by
DASSL is obtained (see equation (79)), except that it is multiplied by a different scalar
factor,
I Nz −hinvk ⋅ I N z
J ode15 s =
−hinvk ⋅ RT ∂Q ∂Q (84)
−hinvk ⋅ RT + RT MR
∂z ∂z
It is necessary to point out that ode15s asks for the matrix M (y , t ) and the matrix
∂f (y, t ) / ∂y separately, in order to be able to form the Jacobian matrix internally. In all
the examples studied, the number of evaluations of the Jacobian matrix ∂f (y, t ) / ∂y has
been far below the number of evaluations of DASSL.
As in the open-chain case, the expression of the Jacobian matrix obtained for closed-
chain systems is similar to the matrix used by DASSL.
I Nz − hinvk ⋅ R z
J ode15 sR =
− hinvk ⋅ RTz RT ∂Q T ∂Q
(85)
− hinvk ⋅ R z R
T
R z + R z R MRR z
T T
∂z ∂z
4 Numerical Results
In the following Sections the examples used to validate the proposed formulations and
integration schemes are described. The results try to reach a double objective: to show
the efficiency of these methods and to present numerical results to consider how the
methods described would behave in practical cases of similar characteristics.
DOF: 41
Bodies: 41
Joints: 40R, 1P
DOF: 5
Bodies: 5
Joints: 5R
The five-bar pendulum has five consecutive Revolute joints and five degrees of
freedom. The length of all the links is L=1 m and they have a uniformly distributed
mass M = 1 Kg. The pendulum, initially at rest, is left free and evolves under the effect
of gravity.
The human body shown in Figure 4 has 39 bodies and 45 degrees of freedom. For
convenience, the model used for computation uses only Revolute joints. Universal joints
were represented as two perpendicular Revolute joints and Spherical joints as three
mutually perpendicular Revolute joints. A Prismatic joint has been added between the
backside and the seat and a Revolute joint between the seat and the ground. Therefore,
the global system has 41 bodies, 40 Revolute joints, 1 Prismatic joint and 41 degrees of
freedom. The interaction with the environment (steering wheel, safety belt and seat) has
been modeled by introducing position and velocity dependent forces.
The driver is assumed to be seated in a natural driving position. To model the
behavior of driver’s muscles, torsion springs and dampers have been introduced at the
joints. Before the frontal collision takes place, the driver is traveling at a constant speed
of 5 m/s. The system described has been simulated in two different conditions. First,
damping in the joints has not been introduced, which leads to a “non-stiff” system. In
the second case, torsion damping (1% of its rigidity to torsion) has been introduced in
the joints. This damping, which is considered to be reasonable and not very high,
introduces “stiffness” in the system equations, so the simulation is considered “stiff”.
The simulation time is 250 milliseconds.
4.1.2 Closed-chain examples
Three closed-chain examples have been considered. The first and second examples are
considered as “non-stiff” simulations, whereas the last one is a “stiff” case.
The Bricard mechanism shown in Figure 5 is essentially identical to the five-bar
pendulum seen in the previous section, except that the last body is hinged to the ground.
All the geometrical and inertial values remain unchanged. By cutting the last joint and
imposing the cut-joint constraints using penalty formulation the closed loop has been
eliminated.
24
DOF: 3
Bodies: 8
Joints: 10S, 1R, 1P
DOF: 1
Bodies: 5
Joints: 6R
DOF: 14
Bodies: 11
Joints: 8S, 2R, 3PR, 1P
The tests have been done with a PC Intel Pentium II, 400 MHz and 128 Mbytes of
RAM, running Windows NT Workstation 4.0. This is a low-cost system that is
accessible to any private person or institution. All the results have been obtained with a
local and global error tolerance of 1e-05.
25
In this Section the dynamic results obtained for the three open-chain examples are
shown. For the sake of simplicity, the following is a definition of a specific notation for
the formulations and integrators used:
• DTR: Fully recursive dynamic formulation for the direct evaluation of the derivative
y or computation of accelerations from known positions and velocities (Section
2.1.1).
• DSR: Semi-recursive dynamic formulation for the direct evaluation of the derivative
y (Section 2.1.2).
• RR-JN: Recursive Residual algorithm rs (Section 3.1). Numerical evaluation of the
Jacobian matrix.
• RR-JF: Recursive Residual algorithm rs (Section 3.1). Analytical evaluation of the
Jacobian matrix (Section 3.4).
The results shown in Table I are the processor times to evaluate the derivative y or
the residual rs .
Table I. CPU time in milliseconds to evaluate the derivative ( y ) or the residual (rs).
In the five bar pendulum, the semi-recursive formulation DSR (order O(N3)) is faster
than the fully recursive formulation DTR (order O(N)) to evaluate the derivative y . In
the driver example, where the number of degrees of freedom is much greater, the DTR
formulation improves the time of the DSR. This concurs with the idea that the best order
O(N3) method is faster than the best order O(N) method for N<11 (Featherstone [3]) or
N<8 (Stelzle et al. [20]).
Table II presents the dynamic simulation CPU times of the examples considered. The
simulation time for the five bar pendulum is 5 seconds and 0.25 seconds for the driver.
Figures in bold print represent the best times for each one of the examples.
26
From the results shown in previous tables, it can be concluded for the open-chain
systems:
• The integration function DE (Shampine and Gordon [19]), based on the Adams-
Bashforth-Moulton methods, always gives the best results for “non-stiff” systems.
This concurs with the experience of many years in the group of Computational
Mechanics of the CEIT.
• Using implicit integration functions and large size problems, the analytical
evaluation of the Jacobian or tangent matrix becomes a key point in the simulation
process. This is not so important with a reduced number of equations.
• Due to the process of resolution of the system of non-linear equations by Newton-
Raphson, with implicit integrators it is convenient to evaluate the recursive residual
rs , instead of calculating the derivative y . The residual algorithm (the inverse
dynamic problem) is faster than the evaluation of the derivative vector.
• DTR formulation is not adapted for “stiff” cases because it is not possible (or it is at
least extremely difficult) to calculate a tangent or Jacobian matrix. On the contrary,
from the terms that appear in formulation DSR this approximate matrix can be
easily obtained in an analytical way.
4.2.2 Closed-chain examples
In this Section the dynamic results obtained for the three closed-chain examples are
shown. As in the open-chain case, a specific notation for the formulations and
integrators used has been defined:
• DSR-Rz: Semi-recursive dynamic formulation for the direct evaluation of the
derivative y using the method of the matrix R z to impose the closure loop
conditions (Section 2.2.2).
• RR-Rz-JN: Recursive Residual algorithm rs using the method of the matrix R z
(Section 3.2). Numerical evaluation of the Jacobian matrix.
• RR-Rz-JF: Recursive Residual algorithm rs using the method of the matrix R z
(Section 3.2). Analytical evaluation of the Jacobian matrix (Section 3.4).
27
The simplicity of programming the DSR method, as opposed to DTR, together with
the fact that the differences of CPU times among them are not relevant (see results
obtained by Jiménez [6], comparing the RTdyn2 (DTR) and RTdyn1 (DSR) formulations
for closed-chain systems), has led to the decision not to implement formulation DTR for
closed-chain cases.
The results shown in Table III are the processor times to evaluate the derivative y or
the residual rs .
Table III. CPU time in milliseconds to evaluate the derivative ( y ) or the residual (rs).
Table IV shows the CPU times required to simulate the proposed closed-chain
systems using the matrix R z method. The simulated time was 5 seconds for the Bricard
mechanism and 2 seconds for the 5 point rear suspension (5PS) and the MacPherson
type front suspension (MFS).
From the results shown in the previous table, some conclusions can be extracted
(some of them agree with the conclusions obtained for open-chain kinematic examples):
• Whereas in “non-stiff” cases (Bricard and 5PS), the DE function provides very good
results, for “stiff” problems it is essential to use an implicit integration scheme.
• Dealing with large systems and implicit integrators, it is recommended to calculate
the tangent matrix in an analytical way. This is not important in small sized systems.
• Using implicit integrators, it is advisable to calculate the residual rs instead of the
derivative y , even if it is necessary to evaluate the Jacobian matrix by means of a
numerical process.
The results presented have made it possible to compare the dynamic formulations
developed with conventional algorithms. The importance of considering together the
problem type, dynamic formulation and the integrator used has been demonstrated.
28
These results would not be complete if they did not have some external reference. It
has been decided to include a comparison with the SIMPACK# code, which is at the
moment the commercial tool of analysis of multi-body systems which is a leader in
technology. This program incorporates many of the most recent advances in the area
and constitutes for that reason a very suitable reference. Like the presented programs,
SIMPACK uses the latest generation of topological formulations and numerical
integrators. Nevertheless, unlike the programs presented here, SIMPACK is written in
Fortran, which constitutes an important advantage in terms of numerical efficiency
regarding the C++ language used in this paper (although it will not have the advantages
of the Object-Oriented Programming).
The comparative results in Table V have been executed in a PC (Intel Pentium® III
500 MHz and 128 Mbytes of RAM, running Windows NT Workstation 4.0), different
from the one used in the rest of the examples. The SIMPACK 8.0 version has been
used.
The driver examples presented do not correspond exactly to that described
previously, since some user defined forces have been eliminated because we did not
have access to this SIMPACK capability. It was decided to eliminate these forces in the
testing examples for both programs.
Table V shows the open-chain results. Of all the numerical integrators available in
SIMPACK, the three most efficient have been selected.
As closed-chain examples (Table VI), it have been compared the Bricard mechanism
and both models of the driver but with an additional distance constraint between the
hands.
#
http://www.simpack.de/html/default.html
29
From the results shown in Table V and Table VI, the following considerations can be
made:
• The Adams-Bashford-Moulton explicit integrator family is always faster in “non-
stiff” problems.
• Implicit Runge-Kutta type integrators (RADAU5) used by SIMPACK work
acceptably in “stiff” systems, but they are inferior to the multi-step implicit schemes
in all the examples considered. It is necessary to notice that the Bricard mechanism
using RADAU5 with relative and absolute tolerances of 10-5 could not be integrated.
The result shown in Table VI corresponds to 10-6.
• Using SIMPACK, the introduction of closed loops in “non-stiff” systems (Bricard)
results in a loss of performance when using Adams-Bashforth-Moulton methods.
The same does not happen with the methods presented in this paper.
• The programs described are always faster than SIMPACK. This indicates that the
use of the R z matrix to impose the closed loop constraints is very efficient.
5 Conclusions
Taking into account the results presented in this paper, the following conclusions can be
established:
1. Global formulations are much less efficient than topological formulations, especially
with large size mechanical systems. The Object-Oriented Programming remarkably
simplifies the implementation of topological formulations, situating the difficulty at
a level similar to global formulations.
3. The problems that appear in the different global formulations due to the presence of
algebraic constraint equations, is reduced very remarkably in the topological
formulations. In fact, if the system does not present closed loops, these difficulties
disappear completely. Whether for open or closed-chains the first velocity
transformation between Cartesian and relative coordinates is made without using the
constraint equations directly. In this first stage, the size of state vector y is reduced
considerably. For closed-chain systems, the constraint equations of closure loop
conditions can be imposed in a second phase with any of the methods used in the
global formulations, but with the advantage of working with a reduced number of
coordinates.
6. Using implicit integrators, the residual algorithm is more efficient than the explicit
evaluation of the acceleration. This is because it is not necessary to form the matrix
of the system nor of course to factorize it. The residual algorithm poses the
equations in a similar way to the semi-recursive methods and therefore it is possible
to calculate an expression of the Jacobian matrix used by the implicit integrators.
7. The simulation times obtained compare very favorably with those of the SIMPACK
program.
6 References
[5] Bae, D.-S. and Haug, E. J., ‘A Recursive Formulation for Constrained Mechanical
System Dynamics. Part II: Closed-Loop Systems’, Mechanics of Structures and
Machines, 15, 481-506, (1987-88).
[6] Jiménez, J. M., Formulaciones Cinemáticas y Dinámicas para la Simulación en
Tiempo Real de Sistemas de Sólidos Rígidos, Tesis Doctoral, ETSII de San
Sebastián (Universidad de Navarra), (1993).
[7] García de Jalón, J. and Bayo, E., Kinematic and Dynamic Simulation of Multbody
Systems – the Real-Time Challenge, Springer-Verlag, New York, (1993).
[8] Jerkovsky, W., ‘The Structure of Multibody Dynamic Equations’, Journal of
Guidance and Control, 1, 173-182, (1978).
[9] Kim, S. S. and Vanderploeg, M. J., ‘QR Decomposition for State Space
Representation of Constrained Mechanical Dynamic Systems’, ASME Journal on
Mechanisms, Transmissions and Automation in Design, 108, 183-188, (1986).
[10] Nikravesh, P. E., and Affifi, H. A., ‘Construction of the Equations of Motion for
Multibody Dynamics Using Point and Joint Coordinates’, in Computed Aided
Analysis of Rigid and Flexible Mechanical Systems, ed. by Pereira, M. F. O. S.
and Ambrosio, J. A. C., 31-60, Kluwer, (1994).
[11] Kim, S. S., ‘A Susbsystem Synthesis Method for an Efficient Vehicle Multibody
Dynamics’, in Advances in Computational Multibody Dynamics, ed. By
Ambrosio, J. A. C. and Schiehlen, W. O., Lisbon, Portugal, September, 379-396,
(1999).
[12] Wittenburg, J., Dynamics of Systems of Rigid Bodies, B.G. Teubner, Sttutgart,
(1977).
[13] Wehage, R. A., ‘Application of Matrix Partitioning and Recursive Projection to
O(N) Solution of Constrained Equations of Motion’, ASME Advances in Design
Automation, 221-230, (1988).
[14] Avello, A., Jiménez, J. M., Bayo, E. and García de Jalón, J., ‘A Simple and
Highly Parallelizable Method for Real-time Dynamic Simulation Based on
Velocity Transformations’, Computer Methods in Applied Mechanics and
Engineering, 107, 313-339, (1993).
[15] Strang, G., Linear Algebra an its Applications, Academic Press, New York,
(1976).
[16] Rodríguez, J.I., Análisis Eficiente de Mecanismos 3D con Métodos Topológicos y
Tecnología de Componentes en Internet, Tesis Doctoral, Escuela Superior de
Ingenieros Industriales (Universidad de Navarra), San Sebastián, (2000).
[17] von Schwerin, R. and Winckler, M., A Guide to the Integrator Library MBSSIM -
Version 1.00, IWR-Preprint 94-75, University of Heidelberg, (1994).
[18] Petzold, L. R., ‘A description of DASSL: a Differential/Algebraic System Solver’,
IMACS Transactions on Scientific Computation, 1, (1982).
[19] Shampine, L. F. and Gordon, M., Computer Solution of Ordinary Differential
Equations: The Initial Value Problem, Freeman, San Francisco, (1975).
[20] Stelzle, W., Kecskeméthy, A., and Hiller, M., ‘A comparative Study of Recursive
Methods’, Archive of Applied Mechanics, 66, 9-19, (1995).