Mechanical Systems of Rigid Bodies Subject to Unilateral Constraints

. The properties of mechanical systems of rigid bodies subject to unilateral constraints are investigated. In particular, properties of interest for the digital simulation of the motion of such systems are studied. The constraints give rise to discontinuities in the solution. Under general assumptions on the system a unique solution is constructed using the linear complementarity theory of mathematical programming. A numerical method for solution of these problems and generalizations of the constraints studied in this paper are briefly discussed.

1. Introduction. The motion of a system of rigid bodies subject to constraints is of interest in many applications, e.g., biomechanics, King and Chou (1976), the performance of vehicles, Magnus (1978), analysis of mechanisms and machines, Paul (1975). It is of practical importance to know how to simulate such systems on digital computers. Sometimes the number of kinematical constraints satisfied as equalities varies in the time-interval of interest. A simple example will illustrate this. Let a rigid bar have one end hinged to the ceiling and the other end free. The hinge constrains the motion of the bar in that the end is not allowed to leave the ceiling, a bilateral constraint. The distance between the end and the ceiling is zero. The motion of the other end may be constrained by the floor. Either the distance between the floor and that end is positive or it is zero when the bar is in contact with the floor, a unilateral constraint. Another simple example is a ball rotating about a point by means of a string. If the string is slack, then no force from the string acts on the ball. The distance between the point and the ball is less than the length of the string. Otherwise, there is a force in the string preventing the distance between the ball and the point from being greater than the length of the string. More complicated systems of this kind are collapsing buildings, Lotstedt and Dahlquist (1978); and sliding rocks in a tunnel or on a slope, Cundall (1974).
Some mathematical properties of mechanical systems subject to unilateral constraints will be investigated here. The theory of complementarity in mathematical programming has proven to be useful in this context. This theory and applications to mechanics are surveyed in Cottle (1979). Other related areas where differential equations and a number of constraints are to be satisfied are discretized problems in variational inequalities, Duvaut and Lions (1972), and discretized contact problems, Kalker (1977).
In § 2 the governing equations are stated and the discontinuities in the solution are discussed. The boundedness of the solution is studied in § 3. The equations in timeintervals where the number of constraints satisfied as equalities is constant are established in § 4. In § 5 the solution in the neighborhood of points at which this number is changed is constructed. The final section contains possible generalizations and some comments on the numerical solution of the problem.
2. The governing equations. The system of ordinary differential equations and inequalities describing the motion of a system of rigid bodies in two or three space dimensions subject to scleronomic, holonomic and nonholonomic constraints are, for t>O, (2.1a) (2.1b) (2.1c) M(q )ij = f(q, q, t) + G(q )A (t), 4J(q)>O, Gf(q)q>O.
The equations without constraints are found in Wittenburg (1977). Kilmister and Reeve (1966) and Peres (1953) give conditions to be satisfied when unilateral constraints are included. q E R m is the vector of coordinates, and fER m contains driving and inertia forces. We assume that ME Rmxm is symmetric, positive definite. However, it is possible that for certain values of q, M is only semidefinite; see, e.g., Huston and Passerello (1976) and Wittenburg (1977). This difficulty can be avoided by redefining the coordinate system in the neighborhood of such a singularity. The components of 4J ERN are the holonomic, unilateral constraint functions. G 2 E Rmxp defines the nonholonomic constraints. The transposition of a matrix or vector c is written cr. Let the Jacobian of 4J be denoted by aq and a column of G 1 and G 2 by gu and g 2 ; respectively. Furthermore, let the constraints be numbered consecutively Gt = (gtb g12, · · ·, gtN) and G2 = (g2.N+b g2.N+2, · · · , g2.N+p).
As an example, {jJ; > 0 may be the geometrical condition preventing a corner of one body from penetrating the edge or surface of another body. A; is proportional to the normal force acting on the bodies at the contact point.
A constraint i is called active if {jJ; = 0 or g~q = 0, while for a passive constraint j, 4Ji > 0 or gL4 > 0. Let those constraints which were active and passive, respectively, in the beginning of the interval T, remain in the same state throughout the whole interval. Then is called the active set in T. The Lagrange multiplier Ai corresponding to a constraint i, satisfies the following complementarity conditions The physical interpretation of (2.2a) in the example above is that if the two bodies are in contact then {jJ; = 0 and A; > 0. If the corner is not in contact with the other body then {jJ; > 0 and A; = 0.
Let us investigate what happens when an active constraint i < N becomes passive at t = to. Assume that the first and second derivatives of 4J exist. For i E J (0, t 0 ], Suppose that for t >to we have l/Ji > 0. Then, at least for a short interval, (t 0 , t 0 + ~tb), ~tb > 0, g[q > 0 since l/Ji(t) = g[(t)q(t)(t-to)+ O((tt 0 ) 2 ) > 0. Similarly, it follows that there is a at 0 < atb such that in (t 0 , to+ ~t 0 ), Ai is nonnegative as long as l/Ji = 0, but at t = t 0 , Ai(t 0 ) = 0, which prevails until l/Ji vanishes again. If ag2i/ aqi exist for all j, then the same arguments can be applied to the constraints (2.1c). Hence we have with gi = gii' i E J(O, t 0 ], j = 1 or 2, in (0, t 0 + ~tb ], and with 8i = g Tii + g T q in (0, to+ ~to], If two or more constraints become passive simultaneously each one of them satisfies the relations (2.3) and (2.4). Note that in [0, oo), Ai, g!q and 8i fulfil the relations Ai > 0, A;g T q = 0, A;8i = 0, without the nonnegativity conditions on g T q for i < N and on 8; for all i. When l/Ji > 0, both g T q and 8i may be negative.
Insert q from (2.la) in the definition of 8;. Let G consist of the columns g;, i E J{O, (2.5) is a system of second order ordinary differential equations (2.5a) and a linear complementarity problem (LCP) (2.5b). See Cottle and Dantzig (1968) and Karamardian (1972) for an introduction to LCP. Since A;= 0, i E J(O, t 0 ]\J(t 0 , t 0 + ~t 0 ], in (t 0 , t 0 +~to] the corresponding column gi in G can be removed without affecting the solution. The complementarity problem (2.2a), which is valid for all t, is in general nonlinear in A since l/>i = l/J;(q), where q is determined by (2.5a).
A constraint is added to the set of active constraints when a previously positive constraint function l/J; or gi;q attains zero. In general, some derivative of the solution is discontinuous when a constraint is added. When the corner reaches the surface of the other body in the example above, a collision occurs and the velocity q and the relative velocity g~q have jumps. There is no problem in deciding when and which constraint to add to the active set since there is a constraint function l/J; or gi;q to base the decision on. The problem of dropping constraints is more delicate and will be treated further in § 5.
In order to discuss the discontinuities in the solution due to the changes in the active set let us consider the system (2.6a) Let a superscript 1 or 2 on a variable denote its value at to-0 or to+ 0, respectively, where t 0 is the point at which the active set is altered. Assume that M, f and G and their first derivatives and the second derivatives of G are continuous functions of their arguments.
After inserting ij from (2.6a) into (2.6c), the equation for A and 8 is as in (2.5), (2.7) GTM-1 GA + GTM-1 / + (;Tq = 8. M-1 has the same continuity properties as M by Lemma A (in the Appendix), q and G are continuous and thus, is continuous as well. It follows from the above arguments that A 2 -A 1 -:f:. 0 and A is discontinuous. Furthermore, (2.6c) implies that ij has a jump.
3. Suppose that no constraint becomes active in the neighborhood of to and that 8i is continuous when a previously active constraint j becomes passive at to. Then 8i(t) = 0, t < t 0 , and 8i(t) > 0, t > t 0 • In order to satisfy the complementarity conditions on Ai (2.4), the following holds: Ai > 0, t < t 0 and Ai = 0, t > t 0 • By (2.7) and the continuity of q and q, A 2 and A 1 satisfy Hence, A 2 -A 1 EN( G), the null space of G, and since 0 EN( G), A can be chosen continuous at t 0 • In particular, Ai(t 0 ) = 0. However, A is usually discontinuous; A]< 0, but AT= 0. q is continuous by (2.6a). By the smoothness of M, G and/, b(q, q, t) and d/ dt The time derivative of (2. 7), after introducing (2.8), is In general, B is discontinuous, and from (2.6a), we infer that ci also is discontinuous.
In the first case the jump condition for q is, Kilmister and Reeve (1966), where G consists of those columns g; associated with active constraints at t 0 and FE Rm contains externally generated impulses at t 0 • External impulses are the counterparts of external forces in (2.1a). The Lagrange multipliers A correspond to impulses in the rigid body system. In the inelastic case they are determined such that the constraints, after the discontinuity, are fulfilled; (2.10b) or, equivalently, (2.11) It is shown by Ingleton (1966) and Cottle (1968) that this LCP always possesses a solution A and a unique il 2 (cf. Lemma 5.2). There are no particular jump conditions in the other two cases, when ij and ci are discontinuous. Equation (2.5) is to be satisfied on both sides of the discontinuity.
In addition to these jumps caused by the changes in the active set, G and f may be sources of discontinuities in 4, ij and i:i-A jump in 4 determined by (2.10) will sometimes be necessary to fulfil all the constraints after such a discontinuity. Consider, e.g., the case when g~ ~ gT, i EJ(O, to). (gf4) 1 = 0 and if (gT)Tq 1 >0 then the constraint i is no longer active fort> t 0 • q 1 is a solution of (2.10) with A= 0. Hence, with F = 0 in (2.10a), then 4 2 =4 1 and 4 is continuous. Conversely, if (gT)T4 1 <0, then 4 is discontinuous satisfying a jump condition such as (2.10). Simple examples are particles sliding on plane surfaces which have L-shaped corners. Also externally generated impulses cause discontinuities in 4. It follows from (2.6a) that ij has a jump when, for some i, A~ > 0, but A 7 = 0 because of the geometry of the system. The corner of a body may, e.g., slide off its support.

Bounds on the velocity vector.
In this section simple bounds on q are derived which are independent of the changes of the active set. A norm is introduced which is related to the kinetic energy of the system. For vERn, llvll denotes the Euclidean norm and for A E Rnxn, IIAII denotes the spectral norm. Furthermore, assume that there are no external impulses, no discontinuities in G and that no holonomic constraints are added to the active set. Then the discussion in the end of the previous section shows that q is continuous.
Since M is symmetric and positive definite, the Cholesky factorization of M exists; (3.3)

M=CCT.
It follows from the construction of C, Dahlquist and Bjorck (1974), that Cis nonsingular and has the same smoothness properties as M. Let A bound on S will be computed and if t is sufficiently small then q and q will remain within the bounds (3.1).
By (3.2), S satisfies (3.5) all4ll < S < bll4ll, and by (2.1a), The term qTGA vanishes identically since A and GTq are complementary (2.3) (or "constraint forces GA perform no work"). If 11411 satisfies (3.1) then by (3.2) there is a bound d' such that liM II< d'. Therefore, Suppose that S > 0, divide the first and the last part of the relation (3.6) by S and use Gronwall's lemma to establish the estimate In the neighborhood of S = 0, and thus q = 0, S 2 is majorized by W satisfying for some k > 0. Hence, when S is small Combining (3.5) and (3.7) yields an upper bound on lliJ(t)ll: where Ct and c 2 are constants. This bound is independent of the number of active constraints and would be the same if the description of the mechanical system (2.1) contained no algebraic constraints at all. Moreover, a change in the rank of G, which sometimes occurs without q being discontinuous, has no influence on the bound.
When holonomic constraints l/Ji have become active, after discontinuities in G and external impulses, A in (2.10) can, e.g., be determined according to Newton's impulse law, Kilmister and Reeve (1966), such that where ei is the coefficient of restitution and superindices are to be interpreted as in § 2. If ei = 1 then we have the perfectly elastic case, and if ei = 0 then the collision is inelastic.
Let lo = {il the impulse law is (3.9), l/Ji(t 0 ) = 0 and (gii)Tqt < 0 or (g~;)Tqt < 0}, and assume that the impulses at the other active constraints at t 0 and the other constraints i for which (g~i) T q 1 < 0 satisfy the complementarity condition,
After substituting A from (4.5) into (4.3a), the governing system of differential equations in (tt, t2) is Using (3.3) and the notation B = c-1 G, a more convenient form is (4.8) and R (A) and N(A) denote the range and nullspace of A, respectively. It follows that P 2 = P and pT = P. P is an orthogonal projection matrix, Ben-Israel and Greville (1974), where

R(P) =N(BT) =N(GTC-T),
Hence, the first term on the right-hand side of (4.7) lies in N(GTC-T) and the second term belongs to the orthogonal space R(C-1 G).
The classical uniqueness theorem for ordinary differential equations can be applied to the particular case (4.6) to obtain: THEOREM 4.1. Suppose that M(q), G(q), aG(q)jaq and f(q, q, t) are bounded and Lipschitz continuous in q and q, that f(q, q, t) is continuous in t and that G has full column rank in the neighborhood of toE (tt, t2), qo and iJo. Then there is a unique solution q(t), q(t) of (4.6) in a neighborhood of to satisfying q(to) = qo and q(to) = iJo.
Proof. By Lemma A and Theorem 2.3 in Coddington and Levinson (1955) the statement follows.

0
The elements of M and G are often sums of products of constants, sin qi and cos qi, where qi represents some angle. Then M, G and aGjaq will all be bounded and Lipschitz continuous. Thus, if the columns of G are linearly independent in [tt, t 2 ] and f(q, q, t) is smooth and bounded, ij will remain bounded as well as q and q. A (t) in (4.5) is continuous, and if Ai(t1) > 0 for all active constraints at t1 then there is certainly an interval (t1, t2) where Ai > 0 and the active set is unaltered.
5. The existence of A and a complementary B. In the previous section we assumed that G had linearly independent columns. A simple example where this restriction is not satisfied is a table with four or more legs on a floor. Therefore, that assumption will be removed when the complementarity problem (2.5) is studied further to show that with sufficiently smooth M, G and f a solution to (2.5) can always be constructed.
To begin with, three lemmas are stated and proved. They will be used in the sequel to investigate the properties of the solution of (2.5).
LEMMA 5 .1. Let x > 0 be the solution of (5.1) G E Rmxn, rank (G)= r <min (m, n) and ME Rmxm is symmetric, positive definite. There are k < r columns of G which are linearly independent and an x' E Rk, x~ > 0, such that for G' E Rmxk, consisting of these k columns, Gx = G'x'.

x-x*EN(G).
There is also a nonsingular G'TM-1 G' with rank ( G'TM-1 G') =rank ( C-1 G') = rank ( G') = k, the number of x *i > 0. Gx = Gx * = G' x' and the lemma is proved. 0 In some mechanical applications the following problem is of interest (see Cottle (1979) LEMMA 5 .2. Assume that rank ( Gt) = q <min (m, p) and rank ( G 2 ) = s < min (m, r). Then a solution to (5.2) always exists with the following properties: 1. 8 and GA are unique. If G has full column rank, then A is also unique. 2. Let A; =A~ +A7 and (Ai)T = ((A~)r, (A~)T), i = 1, 2, where A 1 ER(GT), A 2 E N( G). Then At is unique and A 2 can be chosen arbitrarily but such that A 2 =A~+ A~ > 0.

4.
There is a unique solution Ao such that I lA oii <I lA I I for every solution A.
Proof. 1. To simplify the notation, let Bt = c-tGt and B 2 = c-tG 2 . Any particular A2, At in (5.2a) can be regarded as the solution of a linear least squares problem (5.3) min liB tAt+ B2A2 + CTbll.

A1
The solution of (5.3) satisfying min IIA 1 I I can be written, Ben-Israel and Greville (1974), where A+ denotes the unique Moore-Penrose inverse of A. Furthermore, after eliminating At from (5.2b), the LCP to be solved for A 2 and 8 2 is GA is also unique. If G has full column rank, then GTM-tG is nonsingular, and since G'[b, Gib and 8 2 are unique in (5.2), At and A2 are also unique.
2. GA is unique. Hence, z in (5.6) is also unique.
From linear algebra we know that A 1 is unique. 8 E R ( G T) and for any A 2 , 0 = 8TA 2 = 8JA~.
Since 0 = 8fA2 = 8fA~ +8JA~ = 8fA~, (5.2c) is satisfied for any A~ such that A2 =A~ +A~ >0. 3. It is easily shown that for any solution A to (5.6) there is a solution A' and a matrix G~ with linearly independent columns and rank ( G') = k < s such that GtA t = G~A~.
From (5.5) and part 1 of the proof it follows that Z3 in (5.7) is unique. Apply Lemma 5.1 with M =I, the identity matrix, to (5. 7). Let B ~ consist of those columns in B2 such that (PB2)' = PB~. rank (PB~) satisfies According to the definition of P, R(P) = N(B[) and N(P) = R(Bt). A column bi of B~ can be split into two parts; Since the columns of PB ~ are linearly independent, b T, i = 1, 2, · · · , l, are linearly independent. Therefore, the columns of B~ and G~ are linearly independent. Furthermore, since the columns of C-1 G~ belong to R(B 1 ), G' has linearly independent columns and rank ( G') = k + l. 4. Choose R > 0 such that IIA2II < R, where A2 is the lower part of a solution A of (5.6). The lower parts of the solutions of (5.6) satisfying IIA2II < R and A2 > 0 form a closed, convex set. The convex function IIA2II 2 has a unique minimum A2o on that set. B2A2 in (5.4) is unique. Hence, A to from (5.4) is unique and Ao =(A io, AJo).

0
The next lemma is used to show that the LCP in (2.5) has a solution at a point t = to. LEMMA 5.3. Assume that the rank of G(t) is constant in the neighborhood oft= to and that G has a sufficient number of derivatives. Then the following holds at t = to: If dp T dtv (G v)=O, p =0, 1, · · ·, k-1, k>1, then 1. If X(t) has at least one derivative, then (dv I dtv)(Xr G r v) = xr (dv I dtv)( G r v) = 0. p=O, 1, 2, · · ·, k-1.
2. Any G E Rmxn with rank (G)= r > 0 in the neighborhood of t 0 can be written after permutation of the columns G = (G,, G,X), G, E Rmxr, X E Rrx(n-r\ where G, has full column rank. The kth derivative of XrG;v is according to (5.8) with p = k -1. Hence, dk T dk r dtk (G,v) ) (5.9) dtk (G v) = r dk r = (iry . X dtk (G,v) Since the columns of G, are not linearly dependent, there is a z such that y = o;z. The conclusion from (5.9) is that The complementarity problem to be solved at a certain time point t 0 is (2.5b). We can assume that G T q = 0 at to, since if g i q > 0, then the ith constraint is no longer active for t > to and if g i q < 0, then either cPi > 0 or an impulse must be computed, see § 2 and § 3. Applying Lemma 5.3 with k = 1 to GTq = 0, we find that GTij + GTq E R(GT) and consequently GTq E R(Gr). Then (2.5b) has the form of (5.2), with p = 0 and G 1 =nil. Its solution has the properties listed in Lemma 5.2. In particular, GA is unique, and therefore, ij in (2.5) is unique.
A natural question to pose is: If G T q = 0 at t 0 , is it possible to construct a solution to (2.5) which is valid in an interval [t 0 , t 0 + ~t), ~t > 0? The following observation is important. It is not sufficient to solve (2.5b) at t 0 and let the active set for t > t 0 be l1 = {iiAi(to) > 0} or J2 = {il8i(to) = 0}. If the first alternative is chosen, there may be an index j for which Ai(t 0 ) = 0, but Ai(t 0 ) > 0 and 8i(t) = 0 which should belong to the active set. If j E l2, but Bi(to) > 0 and Ai(t) = 0, then j should not be in the active set of J(t 0 , The theorem in this section proves that under certain smoothness assumptions, there is a solution to (2.5) in [t 0 , t 0 + ~t) for ~t > 0 sufficiently small. The initial conditions at t 0 are We assume that qo and 4o are compatible with the constraints such that cPi(qo) > 0 and g~(qo)4o > 0, and if cPi(qo) = 0 then g~(qo)4o > 0 (possibly after determining an appropriate A in (2.10)).
THEOREM 5.4. Let J = {ilc/Ji = 0 and g[q = 0 or gi4 = 0 at t = t 0 }, and let G consist of those n columns gifor which i E J. Assume that M(q ), G(q) and f(q, q, t) are analytic in t and all qi, qi separately, and that M-1 , G, f and aGj aq are bounded in the neighborhood of q 0 , 4o and to. Furthermore, assume that R =rank (G) is constant in the vicinity of qo.
Then the active set is constant in (to, to+ ~t) for a sufficiently small ~t. There is a unique solution q(t) of (2.5) which is analytic in [to, t 0 + ~t). If G(q) has full column rank, A (t) is also unique and analytic. There is a subset J' c J of cardinality r < R such that Proof. All variables are evaluated here at t 0 unless otherwise stated. Let x<n denote dix/ dti and define Solve the LCP (2.5b) at t = t 0 to determine lAo and 1 80 • According to Lemma 5.2, the solution A can be chosen so that the rank of G' = (gip gi 2 , • • • , gis ), ii E lAo, is s, the number of positive Ai. Moreover, GA and ij are unique. Let us make the following induction hypothesis: For i=O, 1, 2, · · ·, k-1, k>1, (i+ 2 ) . k d .
q ts nown an untque; Gk = (gip gi 2 , • • · , giq), ii Elk= 1\la,k-b and (Gk(q))U) are known and unique; Gkt = (gip gi 2 , • • • , giJ, ii E lA,k-b has full column rank; AU) is known, but not necessarily unique. Let Ak and 8k consist of those components Ai and 8i of A and 8, respectively, for which i Elk. For i elk, A ~n = 0, j = 0, 1, 2, · · · . The equation for the (k + 2) ith time derivative of q is by (2.5a) and the definition of Gk (5.10) G is analytic and all derivatives qU) in (5.11) are known except for q(k+ 2 ). Split Gk into two parts Gkt and Gk2 such that after permutation of the columns and partition Ak and 8k accordingly. It follows from (5 .10) and (5 .11) that there exists a bk such that the system of equations to be solved for A kk) and 8kk) can be written in the partitioned form (5.12) 8~'1 must be equal to zero since Ak 1 (t) > 0 and 8kt(t) = 0 at least for a short time interval T = (t 0 , t 0 +at'). The sign of A ~k!i is of no importance here, since there is a p < k -1 such that A k~)i > 0 and hence in T, Furthermore, A kki and 8k'1 must be complementary since A ~1 = 8~1 = 0, i = 0, 1, ... ,k-1. By Lemma 5.2, (5.12) has a solution where 8~k) is unique. There is a matrix G~2 such that (5.13) has full column rank.
Determine J>..k, lsk and Jk+1· By the definitions of these sets we have J>..i c J>..,i+b lsi c ls,i+1 and Ji+1 c Ji. q<k+ 2 ) is known and unique, and A (k) is known. Gk+1,1 is specified in (5.13). The inductive step is completed.
Since the number of indices in J is finite, eventually in a step the iterative procedure for determination of the active set can be terminated for three reasons: 1. Jk\JA.k = 0. A '(t) consists of those Ai for which A ~o) > 0 or A ~n > 0, but A ~k) = 0, k = 0, 1, · · · , j -1. A'(t) satisfies (cf., (4.5)) A'(t) = -(G'TM-1 G')-1 (G'TM-1 f + G'Tq), and is analytic. After permutation of the columns in G, is a feasible solution. The computed A ~0 and 8 ~0 determine A (i) and 8 U). If the iteration terminated due to the third criterion, 8i(t) may become negative for some i when using (5.16) for q(t), viz., if 8~n = O,j = 0, 1, · · ·, p > k but 8~p+ 1 ) <0. The iterative procedure must be continued. In that situation compute A lk), A lk+ 1 ), · · · with A&k) = A&k+ 1 ) = · · · = 0 until p + 1 is reached. Determine the unique q<k+ 2 ) and the sets J>..k and lsk and proceed as before. The iteration is eventually terminated, since the number of 8i is finite, 8(t) is analytic and if 8i(t) > 0 in T there must be an i(j) such that 8 } 0 > 0 for each j.
When the iteration has terminated due to criterion 1, 2 or 3, then The derivatives q<n(t 0 ), i = 0, 1, 2, · · · are calculated uniquely using the iterative process. They will coincide with the derivatives at to for the solution of (5.16). If no previously passive constraint becomes active in T, then the solution of (5 .16) is also the unique solution of (2.5) in f. However, since the number of passive constraints is finite, there is a tit 0 <tit< tit' such that the active set in (to, to+ tit) is {il5~n = 0, j = 0, 1, 2, · · · }. A possible but not always unique A (t) is given by (5.17). If G has full column rank, then G~2 in (5.13) is unique in each step and G' as well as A are unique.

0
The proof can easily be modified to allow some constraints i to be bilateral, i.e., there is no sign restriction on Ai and 5i = 0, t > 0.
The active set has been determined at a point t1 E (to, t 0 +tit) and in its neighborhood. If G has not full column rank then we can find a G', where G'TM-1 0' is nonsingular, yielding the same solution. As long as the active set and rank (G) are constant, the system of differential equations to be solved is (4.6) with G = G' commencing from t1 with q(t1) and q(t1) as initial values. The uniqueness of the solution is then also given by Theorem 4.1.
It follows from the remark after Theorem 4.1 that M and G are often analytic in all their arguments qi separately. In many applications f is at least piecewise analytic.
The rank of G is not necessarily constant between the updates of the active set. If G is slightly perturbed so that the rank is changed at t 0 , the solution before to may be completely different from that obtained after t 0 although the active set is constant. It is easy to construct simple two-dimensional examples which illustrate this fact. Fortunately, these problems do not seem to arise very frequently in practice. For most systems G has the representation G(q)=A(q)B(q), AERmxR, BERRxn, where the columns of A are linearly independent and rank(B) = R. Such a G satisfies the rank condition in the theorem.
6. Concluding remarks. A solution to the linear complementarity problem (2.5b) is also a solution to the quadratic programming problem (6.1) and vice versa, see Cottle and Dantzig (1968). The dual problem of (6.1) is, Dorn (1960), (6.2) After inserting the expression for GA from (2.5a) into (6.2), the new minimization problem is (6.3) G T··+a·T·>o q q_ .
In (6.1-3) q and q are fixed variables. Equation (6.3) is a generalization of the principle of least constraint by Gauss, Whittaker (1944), to rigid body problems with unilateral constraints. The formulation best suited for numerical computation seems to be (6.1) or (6.3). In order to compute approximations to q(t) and q(t) at discrete time points, the time-derivatives are replaced by difference approximations. Then there are reliable algorithms for solving the optimization problem. For further details see Lotstedt (1979b).
In some applications there are upper and lower bounds on the multipliers (6.4) The corresponding constraint function g T q satisfies if A L < Ai <Au then g T q = 0, (6.5) if Ai =A L then gFq > 0, if A;= Au then g T q < 0.
The problem treated in § § 2-5 is a special case of (6.4) and (6.5), where A L = 0 and Au= +oo. The Coulomb friction forces, Kilmister and Reeve (1966), fulfil the relations (6.4) and (6.5). ALand Au depend on the normal force AN> 0 associated to the friction force in the following way where JL > 0 is the friction coefficient. The results concerning the discontinuities in § 2 are applicable also when (2.5b) is replaced by conditions like (6.4) and (6.5). When A L and A u are constants, theorems similar to those in § 4 and § 5 can be proved to cover also systems governed by (2.5a), (6.4) and (6.5). In the Coulomb friction case, however, there are systems where no solution A exists and if it exists q is not necessarily unique.
This fact was first pointed out by Painleve (1895).
Appendix. The lemma is stated for reference purposes only. It can be proved by using Cramer's rule for inversion of matrices. LEMMA A. Let A(q)ERnxn be nonsingular for llq-q*ll<b and have one of the following properties: I. A hasp> 0 continuous derivatives in all arguments qi, i = 1, 2, · · · , m.

II.
A is bounded and Lipschitz continuous in llqq *II< b. III. A is bounded in llqq *II< b and analytic in all the arguments separately in llq-q*ll<b.
Then A -l also has the same property.