Sequential Regularization Methods for Simulating Mechanical Systems with Many Closed Loops.
The numerical simulation problem of large multibody systems has often been treated in two separate stages: (i) the forward dynamics problem for computing system accelerations from given force functions and constraints and (ii) the numerical integration problem for advancing the state in time. For the forward dynamics problem, algorithms have been given with optimal, linear complexity in the number of bodies, in case the system topology does not contain many closed loops.But the interaction between these two stages can be important. Using explicit time integration schemes, we propose a sequential regularization method (SRM) that has a linear complexity in the number of bodies per time step, even in the presence of many closed loops. The method also handles certain types of constraint singularity.
Introduction
There has been a growing interest in the development of more efficient algorithms for
multibody dynamics simulations. The increase in size and complexity of spacecraft
and robotic systems is one motivation for this development; another is physically-based
modeling in computer graphics. The numerical simulation process has been
Institute of Applied Mathematics and Department of Computer Science, University of British
(ascher@cs.ubc.ca). The work of this author was
partially supported under NSERC Canada Grant OGP0004306.
y Institute of Applied Mathematics, University of British Columbia, Vancouver, B.C., Canada
V6T 1Z2. New address: Program in Scientific Computing and Computational Math, Stanford
University, Durand 262, Stanford, CA 94305-4040, USA (plin@sccm.stanford.edu).
typically treated as two separate stages. The first stage consists of the forward dynamics
problem for computing system accelerations, given the various constraints,
torque and force functions. For tree-structured multibody systems, algorithms have
been proposed with optimal O(n) complexity, where n is the total number of rigid
bodies in the system (see e.g. [10, 13, 18, 6, 21]). They have also been extended
to cope with systems with a small number of closed loops compared with the total
number n of links [10, 18]. For a system with m closed loops (m ! n), a typical
complexity of O(n obtained using a cut-loop technique. But it appears to
be hard to find an O(n) algorithm for chains with a large number (e.g.
closed loops.
The second stage of the simulation algorithm design addresses the numerical integration
problem for advancing the state in time, obtaining generalized body positions
and velocities from the computed accelerations. Explicit or implicit time discretization
schemes can generally be used.
While these two stages are usually treated separately, there are situations in which
the specific treatment of one affects the other, so a global, unified view is beneficial
(e.g. [6, 1]). In this paper we use such a global view of the simulation process and
devise a method which requires O(n) operations per time step even in the presence of
many closed loops. Specifically, we propose using a sequential regularization method
(SRM) [4, 5] for this purpose, combined with an explicit time integration scheme. The
method produces iterates which get arbitrarily close to the solution of the discretized
differential system, and it also handles certain types of constraint singularity.
The mathematical modeling of constrained multibody systems yields differential-algebraic
equations (DAEs) of index 2 or 3 [9, 11]. For tree-structured systems (i.e.
no closed loops), one can formulate the model in terms of a minimal set of relative
solution coordinates, obtaining a system of ordinary differential equations (ODEs),
see e.g. [14, 10]. Some existing commercial software packages (e.g. SD/FAST 1 ) utilize
this approach. Forward dynamics algorithms of complexity O(n) can be interpreted
then as imbedding the ODE in a DAE, at a given time, and eliminating some of the
unknowns locally in the larger but sparser system [6].
In this work, however, we consider the equations of motion in descriptor form (see,
e.g., [9, 14]). These are formulations in non-minimal (redundant) sets of coordinates
which yield an often simpler, albeit larger, DAE, even in the tree-structured case.
This DAE is typically treated by differentiating the constraints to the acceleration
level and using one of the well-known stabilization techniques (see, e.g. [2, 3]). An
O(n) forward dynamics algorithm for this formulation is recalled in x2, following [18].
A system with closed loops may now be considered as being composed of a tree-structured
system plus a set of loop-closing constraints. The latter are treated using
an SRM technique. The method is described in x3, where we prove that the number
of operations needed per time step when using an explicit time discretization scheme
1 SD/FAST is a trademark of Symbolic Dynamics, Inc., 561 Bush Street, Mountain View, CA
USA.
remains O(n), for any m - n. The method also handles certain types of constraint
singularities.
Finally, in x4 we demonstrate our algorithm on two closed-loop chain examples.
Algorithms with optimal computational complexity
for tree-structured problems
2.1 General multibody systems
Consider an idealized multibody system consisting of rigid bodies and point masses
with a kinematic tree-structure. We use redundant, world coordinates p for the positions
of the system, i.e. for describing the position and orientation of each individual
body. The set of feasible positions, which correspond to physically possible geometric
configurations, is given by holonomic constraints,
Differentiations of the constraints with respect to time yield corresponding constraints
for the velocities -
and the accelerations -
Here G(p) is the constraint Jacobian matrix g 0 (p). The Euler-Lagrange equations are
In combination with the constraint equations at the acceleration level (2.3) the Euler-Lagrange
equations form a DAE of index 1 for the differential variables p, v (= -
and the algebraic variables -:
The constraints on the position level (2.1) and on the velocity level
define an invariant manifold for (2.5).
It is well-known that simply simulating (2.5) numerically may cause severe drift
off the constraints manifold, manifesting a mild instability in this formulation. This
Figure

2.1: An example of tree-structured chains
phenomenon can be prevented by using stabilization techniques (e.g. invariant stabi-
lization, Baumgarte's stabilization and projection methods [2, 3, 7]) or regularization
techniques (e.g. penalty and sequential regularization methods [17, 19, 12, 15, 16,
4, 5, 20]) during the numerical integration. In this paper, we assume that such
a stabilization technique has already been applied and is included in (2.5), and we
concentrate on aspects of solving this system.
2.2 Algorithms for tree-structured problems
Consider the case of a multibody system with a kinematic tree structure. As in [18],
consider a graph whose nodes correspond to the joints and whose edges represent the
bodies in the system (see, e.g. Figure 2.1). When there are no closed loops, the graph
consists of trees. In each tree one joint is singled out as the root. Every other joint
then has a unique father in the tree, which is its neighboring body on the path to
the root. We introduce a node "0" as the father of the roots (if there is a fixed joint,
we usually number it as "0"). This gives one tree, and we label its nodes as follows:
joints are numbered from 1 to n such that the label of a joint is always greater than
that of its father. Bodies (links) are numbered such that a body connecting a father
joint and a son joint has the number of the son. For example, in Figure 2.1, joint
is the father of joint 7, and these two joints are connected by body 7.
For this kind of tree-structured systems, the mass matrix M is symmetric positive
definite and block-diagnal with blocks which are symmetric
positive definite. The constraint equations in (2.3) of a body connecting two joints k
father of k are of the form (see [18]):
where
is assumed to have full row rank
for each k. Then, as analyzed by many articles (see, e.g. [10, 13, 18, 21]), a recursive
algorithm can be derived to obtain an O(n) algorithm to invert the left-hand matrix
in (2.5b).
To explain the algorithm, we start at a terminal joint i.e. one which has no
Setting -
, we can obtain
Next, consider the equations for a joint j which has terminal sons k:
k=sons of j
father of j. For each son, (2.9) holds. Substitution into (2.10) yields:
where
k=sons of j
k=sons of j
The system (2.11) has the same form as (2.8), so a recursive algorithm results:
Algorithm
step 1: Climb down the tree (towards the root) recursively, repeatedly forming
f j by (2.12).
step 2: Starting from the root of the tree, solve each local system (2.11) for u j
(with climbing up the tree recursively.
Note that in Step 1, we replace the original problem by n local problems (2.11).
This algorithm has O(n) complexity because all operations are local for one body
and we only climb down and up the tree once. The work for different subtrees can
be done in parallel as well.
The result of the forward dynamics algorithm just described is an expression for
v in terms of
and p. This is an ODE system that must be integrated in
time. Note that the dimension of each of u; v and p in 3D is 6n. This is contrasted
with the treatment using minimal coordinates, where u; v or p may have dimension
as small as n. However, the saving in the latter formulation is chiefly in the time
advancement, and in a faster forward dynamics algorithms for small n, but not in the
O(n) forward dynamics algorithm (see [6]).
3 Algorithms for systems with many closed loops
We consider now multibody systems with loops. The problem can be seen as a tree-structured
system plus some extra loop-closing constraints in the form of
where x k and x j represent the generalized positions of joints k and j, respectively,
and are components of p, and m represents the number of loops. We can write (3.1)
as
where D is an m\Thetan constant matrix which has a full row rank. In the closed-loop case,
we need to invert a matrix -
instead of the matrix
in (2.5b), where -
(D 0). Block Gaussian elimination gives
where the Schur complement -
D T can be obtained using the O(n) algorithm
described in the last section. Using usual techniques for inverting -
, an O(n+
results which is satisfactory when m is small (cf. [18]).
For problems with a large number of loops, it appears difficult to extend this
O(n) algorithm (cf. [10]). In particular, no O(m) algorithm for inverting the Schur
complement is available. Also, we have found no discussion in the literature about
the performance of the method when the constraint Jacobian matrix
is rank
deficient. This situation often happens in closed-loop chains, for example the slider-crank
problem [4, 5] and the first example in x4. The sequential regularization method
proposed in [4, 5] provides a possibility to address these two challenges.
3.1 The algorithm
We now write down the Euler-Lagrange equations for a multibody system with loops
as
(constraints corresponding to a tree structure) ; (3.3c)
where the number of rows of D may be large. We assume that G and D have full
rank. This assumption is generally true for chain problems (see examples in x4). We
already have an O(n) algorithm for the tree-structured problem, i.e. for the system
(3.3a), (3.3b) (without D T -) and (3.3c). So, to develop an O(n) algorithm for the
tree-structured part of the problem we keep the structure of (3.3a), (3.3b) and (3.3c).
For the loop-closing constraints, one possibility is a penalty method:
and
is the penalty parameter. However, it is well known that we have
to use implicit integration schemes for a system like (3.4) since ffl has to be very
small and then (3.4) is a stiff system. Therefore we have to invert a matrix like
for a discretization stepsize h
D T D is
not block diagonal any more, the previous algorithm will be difficult to apply.
In [4, 5] we proposed a new, sequential regularization method (SRM), which is
a modification of the usual penalty method and allows us to use explicit schemes to
solve the regularized problem. Hence, SRM combined with explicit time discretization
makes it possible to obtain an O(n) algorithm for problems with many loops.
Applying the SRM to problem (3.3) to treat the loop-closing constraints, we obtain
a new algorithm: given - 0 (t), for such that
and
If we apply a stabilization method (e.g. Baumgarte's stabilization [7], or the
first stage of the invariant stabilization [2, 3], or projection methods [2]) first for the
constraints (3.3c), then the following system must be considered,
and (3.6c) becomes:
An even simpler algorithm was suggested in [5] for general multibody systems
without singularities. For (3.3), it has the form:
and
This form has O(n) computational complexity and it is simpler than the form (3.6)
since we only need to invert the block diagonal matrix M . However, for problems
with singularities (or ill-conditioning) which often appear in closed-loop chains, the
form (3.10) is not recommended, as we indicated already in [5]. In x4, a numerical
example shows that its performance is worse than that of the form (3.6). Also, it
is indicated in [5] that for the form (3.10) the best choice of ffl is generally ch (for
the reason of stability of the explicit difference schemes), where h is the step size
of the chosen difference scheme and c is a constant dependent on the eigenvalues of
the matrix
G
G
. Hence, the best choice of ffl for this form changes in
time because G is in general dependent on the time t. We will see an example in x4
where ffl has to be chosen to be quite large (hence more SRM iterations are needed)
to make the explicit discretization stable for the whole time interval that we compute
on. Therefore we still recommend the form (3.6) since it performs better as shown
in x4 and since the best choice of ffl for this form depends only on the eigenvalues of
the matrix DM (which is constant if M is a constant matrix, as is the case for
many chain problems, including the examples in x4). That is, the best choice of ffl is
often independent of t for the form (3.6). In comparison with the usual stabilization
methods, an additional advantage of the from (3.6) is that we never invert a singular
matrix even if
is singular, since G has full rank as we assumed.
To summarize, our algorithm consists of applying an explicit time discretization
scheme and the O(n) forward dynamics algorithm for the underlying tree-structured
system to solve for
where
Next we discuss the convergence of the SRM and prove that the iteration number
of the method is independent of the body count n. Hence, per time step our method
is an O(n) algorithm even for chains with many closed loops.
3.2 Convergence of the stabilization-SRM form
Now we want to analyze the convergence of the iterative procedure (3.12). The
method of analysis is based on that in [5]. The differences here are that we only
apply the SRM to part of the constraints (i.e. the loop-closing constraints (3.1)) and
that the problem we consider may have a large number of unknowns (linear in n).
As indicated in [4, 5], the system (3.12) is clearly singularly perturbed for
ffl - 1. Starting from an arbitrary - 0 (t) we may therefore expect an initial layer. For
simplicity of the convergence analysis we assume (as we did in [4, 5])
Assumption 3.1
For initial value problems it is possible to obtain the exact - j (0);
in advance [4, 5]. For a general semi-explicit DAE, under this assumption there are
no initial layers for the solutions of (3.12) up to the lth derivatives.
Now let us make assumptions on the boundedness of the solution. For an n-body
chain, it is reasonable to assume that the solution and its derivatives are bounded
linearly in n. In other cases the solution may be bounded independently of n. We
make corresponding assumptions on p s , v s and - s since our regularization is a nearby
problem to the original system and there will be no initial layers involved in the
solution under the condition (3.13) (letting l - 0).
Assumption 3.2 The solutions p s , v s and - s of (3.12) satisfy
where - may depend on n. From (3.12b) and (3.12c) we can solve for - s in terms of
i.e. we can write
We further assume that the Jacobian matrix of   with respect to p s , v s and - s satisfies
If - s also satisfies (3.14) then (3.14) and (3.16) imply that
represents the maximum norm and K is a generic positive constant
independent of n.
In Remark 3.1 below we give convergence bounds for the case where we only assume
We have noted that the mass matrix M is positive definite and block-diagonal
and that the closing-loop constraint matrix D is block-sparse. From the examples in
x4, we can see that the Jacobian matrix G for the constraints (2.1) is block-sparse
too. We thus assume:
Assumption 3.3 The matrices M , D and G and their finite multiplications (i.e. the
number of multiplications is independent of n) are all block-sparse. More precisely, if
we multiply such a matrix and a vector with norm O(n k ), then the product norm is
still O(n k ), where k is any positive number.
In fact the closing-loop constraints matrix D is not only block-sparse but also block-
row-orthogonal. More concretely, there is at most one nonzero element in each column
(see, e.g. (4.6) in x4). Combining this with the properties of M we thus obtain that
DM positive definite and block-diagonal. Hence the eigenvalues of the matrix
DM are the union of the eigenvalues of diagonal blocks of the matrix. Thus
all these eigenvalues are positive and independent of n since the size and elements of
diagonal blocks of DM are independent of n. We write these as a lemma:
Lemma 3.1 If the mass matrix M and the closing-loop constraint matrix D have
the above properties then the eigenvalues of the matrix DM \Gamma1 D T are all positive and
independent of n.
We finally assume the following perturbation bound (cf. the perturbation index
[11] and corresponding discussion in [5]).
Assumption 3.4 Let -
. Then there exists a
constant K independent of n such that:
where z is the solution of (3.8) and - z satisfies the perturbed (3.8):
Here, although inverting
needs O(n) operations (see the algorithm in
x2), these operations will not involve the perturbations '(t) and ffi(t). So we believe
the bound K in (3.18) to be independent of n.
Now we can consider the convergence of the stabilization-SRM form (3.12). We
choose - 0 such that - 0 and -
- 0 are of O(-). Let u . Then we
have drift equations:
or
with the initial conditions u s Applying (3.21) for
using the Assumptions and Lemma 3.1 we obtain u
(3.21a) further yields
Comparing (3.19) with the
stabilization-SRM form, we need to bound
and their derivatives appearing in (3.18). We already have that
From (3.21a) we obtain -
Using the condition - 0
DM
can be obtained from (3.12b) and (3.12d). So, from (3.21b), -
Differentiating (3.21b) we have
where we use - O(-) obtained from (3.12d) and w
obtained from (3.17). Hence -
Differentiating (3.21a) we next have
This implies
We thus use (3.18) and obtain the desired conclusion for
Subtracting (3.12a), (3.12b) from (3.8a), (3.8b), respectively, and using (3.22) and
(3.16) we have
Also, by differentiating (3.12d), we conclude that -
1 is of O(-).
For first have as for the
and hence also
The estimate (3.23) also holds for 2. This yields that the right hand side of
(3.21b) is O(-ffl 2 ), so
Now we want to get a better estimate for -
Differentiating (3.21b) we have
and -
obtained as for the
We want to show that -
O(-ffl). For this purpose we must estimate -
Using
the condition -
-(0), and the fact that
are also exact at
0, we can obtain
Hence -
Differentiating (3.24) again we now obtain precisely as
when estimating -
Here we need to use -
which can be obtained from a differentiation
of (3.12), applications of (3.14), (3.16) and (3.17) as well as the estimates for -
obtained above. Noting
we thus have
Our previous estimates allow the conclusion that
hence we can conclude that -
Repeating this procedure, we can finally obtain:
Theorem 3.1 Let all the assumptions described at the beginning of this section hold.
Then, for the solution of the stabilization-SRM form (3.12) with - 0 and -
by O(-), we have the following error estimates:
is the solution of the multibody
problem (3.3).
Remark 3.1 If the bounds in (3.14) are only assumed for p s , v s and - s , then we may
generally obtain -
and the j-th derivatives would be O(- j+1 ). This
finally yields the following error estimates:
This result is obviously weaker than (3.27), corresponding to weaker boundedness assumptions

3.3 Computational complexity of the stabilization-SRM form
Consider our algorithm for (3.12). At each regularization iteration we can use the
O(n) algorithm described in x2. So, to consider the computational complexity of our
algorithm we need only study the number of iterations s required.
It is simple to show that s is independent of n. Given the worst error estimate
(3.28), we must choose ffl small enough so that
Then each SRM iteration reduces the error, viz. the difference between the solution
of (3.12) and the solution of (3.3), by at least a factor ff, and so a fixed number of
iterations s, independent of n, is needed to reduce this error below any given tolerance.
Note, though, that - may grow with n, depending on the problem being simulated.
Hence the range of choices for ffl (ffl - ff=-) is restricted depending on n. Since the
time discretization scheme is explicit, the step size h must be restricted by absolute
stability requirements to satisfy
for an appropriate constant fl of moderate size. The number of time steps required
may therefore depend on n, too (cf. [21]) 2 . Still, the number of iterations s required
to obtain a given accuracy obviously remains independent of n, hence operation count
per time step remains O(n).
For the worst error estimate (3.28) suppose that
where fi is a given positive constant, simplicity. Now
let us apply an r-th order explicit difference scheme to the stabilization-SRM form.
At the s-th iteration, the worst combined error for our algorithm is O(- s ffl s
(see Remark 3.1). Trying to roughly equate the two sources of error, we set
hence
Using (3.29) we then obtain a rough upper bound for the number of iterations s:
A requirement such as likely to arise also from accuracy, not only stability
considerations.
Remark 3.2 For the error estimate (3.27) the condition (3.29) can be weakened
significantly to read
where s 0 is an arbitrary integer independent of n. In this case, (3.30) is replaced by
4 Numerical experiments
We now present a couple of examples to demonstrate the algorithm that was proposed
and analyzed in previous sections. At first, we build up the system for a special kind
of n-body chains (see Figure 2.1) which include our two examples. We use the method
described in [21].
Consider a chain consisting of n bodies. Each body is modeled as a line segment
of length l j and mass m j , with uniform mass distribution. We choose Cartesian
coordinates of the joints, x j , and the vectors connecting the joints along the links,
in order to describe the position p of the chain:
n. For 3-D chains, p j has six components and for 2-D chains
four components. The labeling of the chain joints has been shown in Figure 2.1 and
explained in x2. Hence the holonomic constraints include length conditions
and connection conditions
father of j. Therefore
Due to the uniform mass distribution, the center of mass of each body coincides with
its geometric center, x and the moment of inertia about the center of mass
is . The total kinetic energy is given by
I
3 I
where
I

Figure

4.1: A square chain and its labeling
and we let the joint "0" is fixed. The potential energy due to gravity is given by
ge T
Here e v is the unit vector along the vertical axis. For the 2-D case e
for the 3-D case e . Note that e T
is the height of the center
of mass of body j above zero level, and -
is the gravity constant. This
gives the force vector
f =B @
Hence, we know the mass matrix M , the force term f and the constraints g and
can write down the system (2.5) for this kind of tree-structured chain problems. For
chains with loops, we only need to impose some additional geometric constraints onto
their corresponding tree-structured formulation (see the discussion at the beginning
of x3).
Next we consider two specific examples.
Example 4.1 Consider a 2D square chain with unit length and mass for each link
(i.e. 1). The labeling of a corresponding tree structure is shown in Figure
4.1. Under this labeling,
Hence we can write down the Jacobian matrix G(p) for the tree-structure constraints
g(p). For example, when
I \GammaI
and
Here I is a 2 \Theta 2 unit matrix. The extra loop-closing constraints d(p) are:
where m is equal to the number of squares. We thus can form its Jacobian matrix D.
For example, when
The Jacobian matrix (G(p) T ; D T ) T of all constraints has rank deficiency when all four
links of a square are on a line. We let the chain fall freely from the position shown
in

Figure

4.1 where the joint "0" is fixed. As we have mentioned before, for each
SRM iteration our algorithm (3.12) solves a tree-structured problem using the O(n)
algorithm of x2. Here, we only demonstrate that the number of SRM iterations is
independent of n. This means that the computational complexity of our algorithm
is O(n) per time step. We choose step size regularization parameter
apply an explicit second-order Runge-Kutta method to the regularized
problem at each iteration. We do the computations for because to clearly see a
relation between the number of iterations and n we want to avoid the singularity which
happens around and after whose error situation is very complicated. We
count the number of SRM iterations until the errors in the constraints do not exhibit
obvious improvement. Table 4.1 lists iteration counts and constraint errors for various
n.
# of
iterations
# of
iterations 3 3 3 3 3 3 3

Table

4.1: Relation between the number of iterations and n for the square chain
Theoretically we expect that two SRM iterations be sufficient since the difference
scheme is of second order and O(h). From the table we see that at the second
and the third successive iterations the maximum drifts are almost the same, especially
when n - 1=h. Additional experiments with
the need of only two iterations for the algorithm (3.12) to achieve the second order
discrete accuracy.
Next we compare the performance of the algorithms (3.12) (AlgI) and (3.10) (Al-
gII). We set computational results show that the first singularity occurs
after We use the simple forward Euler scheme and take
for both algorithms. We can still take for the algorithm (3.12). But for the
algorithm (3.10) we cannot take :5h. The algorithm is unstable immediately after
the first time step when we take :005. The algorithm blows up around
when we take It becomes stable when we take This agrees with
our expectation about the algorithm in x3.1. That is, for the sake of stability of the
difference scheme the smallest ffl we can choose depends on t and it is often larger
than that needed for the algorithm (3.12). Hence, more iterations are needed for the
simpler algorithm (3.10).
algorithm # of iterations ffl t=.4 t=.5 t=.6 t=1.0
AlgI 1 .005 7.763e-3 1.289e-2 1.744e-2 2.959e-2
AlgI 2 .05 7.763e-3 1.289e-2 1.744e-2 3.107e-2
AlgII 2 .05 1.008e-1 1.640e-1 2.270e-1 4.794e-1
AlgII 4 .05 5.869e-2 8.514e-2 1.129e-1 1.155e-1

Table

4.2: Maximum drifts of two algorithms
We list in Table 4.2 the maximum drifts produced by these two algorithms at
various times. From the table we see that the overall performance of AlgI is better
than that of AlgII. Also it seems that for this example the error improvement of AlgII
by iteration is much slower. The motion of this square chain with
to 2:2 is described in Figure 4.2.
Example 4.2 In this example we consider the motion of a square net (veil) started
by a breeze. The corresponding tree structure of the net and its labeling are shown in

Figure

4.3, where are fixed joints. Under this labeling,
we can determine the father of any given joint and then the Jacobian matrix G(p) for
the constraints The extra loop-closing constraints are:
ae x
For this square net the number of bodies l. Again we set the step size
and the regularization parameter 0:5h, and apply the O(n) algorithm at
each time step of the second-order Runge-Kutta method for the regularized problem at
each SRM iteration (3.12). After a breeze we assume that the net moves to the right
with the largest angle -to the verticle axis. We consider the motion of a net whose
initial position is located at where it has the angle -with respect to the vertical line.
We list the number of iterations for various n at in Table 4.3.
From the table we do not see drift error improvement after the second iteration.
only two iterations are again needed for the algorithm (3.12), independent of n.
The motion of this square chain with
described in Figure 4.4.
Figure

4.2: Motion of the square chain. Time increases in increments \Deltat = :2 from
left to right, top to bottom.
2l
4l
2l +3
2l +4
2(l +l)-3
2(l +l)-2
2(l +l)
2l +l
l =(w-1)l
l =(w-2)l
Figure

4.3: Tree structure and its labeling of the square net
# of
iterations
# of
iterations 3 3 3 3 3 3 3

Table

4.3: Relation between the number of iterations and n for the square net



--R

Recent Advances in the Numerical Integration of Multibody Systems
Stabilization of invariants of discretized differential systems
Stabilization of DAEs and invariant manifolds
Sequential regularization methods for higher index DAEs with constraint singularities: Linear index-2 case
Sequential regularization methods for nonlinear higher index DAEs
formulation stiffness in robot simulation
Stabilization of constraints and integrals of motion in dynamical systems

Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations
Robot Dynamics Algorithms
Solving Ordinary Differential Equations II
Regularization of differential
Unified formulaion of dynamics for serial rigid multibody systems



On a penalty function method for the simulation of mechanical systems subject to constraints

Stabilization of computational procedures for constrained dynamical systems
Numerical solution of differential-algebraic equations with ill-conditioned constraints
Fast recursive SQP methods for large-scale optimal control problems
--TR
robot simulation;stabilization;regularization;constraint singularities;differential-algebraic equations;multibody systems;higher index