The Dynamics of Discrete Mechanical Systems with Perfect Unilateral Constraints

Abstract: The dynamics of discrete mechanical systems with perfect unilateral constraints is formulated in a very general setting. The well-posedness of the resulting evolution problem is studied. It is proved that existence and uniqueness of a maximal solution is ensured provided strong assumptions are made on the regularity of the data: they are supposed to be analytic. Simple examples show that this regularity assumption may not be relaxed. Sufficient conditions to ensure that the maximal solution is defined for all time are supplied. The continuous dependence of the solution on initial conditions is also studied and the numerical computation of the solution is discussed.


Introduction
The aim of the Dynamics of Discrete Mechanical Systems (sometimes called Rational Mechanics or, after Lagrange, Analytical Mechanics) is the prediction of the motion of collections of bodies supposed to be perfectly indeformable. The theory classically distinguishes two types of interactions between the bodies themselves and between the bodies and the rest of the universe: the efforts and the constraints. The constraints are kinematical specifications of the motion with which some efforts are associated. A constraint is said perfect or ideal if the associated efforts do not dissipate energy. A constraint is said bilateral (respectively unilateral) if the kinematical specification gives rise to equalities (respectively inequalities). A typical occurrence of unilateral constraints is the handling of non-penetration conditions.
When all the constraints are bilateral and perfect, the motion is classically governed by a second-order ordinary differential equation on a finite dimensionnal riemannian manifold. When the data are smooth enough, Cauchy-Lipschitz theorem guarantees that a unique motion is associated with any given initial state of the system.
When dealing with the dynamics of discrete mechanical systems with unilateral constraints, there is no such a theorem, although many steps towards this direction have been made during the past twenty years. To my knowledge, the first investigation of this question using modern mathematical tools (i.e. introducing motions whose acceleration is a measure with respect to time) is that of Schatzman [18]. She studied the particular case where the configuration space is R d equipped with its canonical euclidean structure and the admissible configuration set is convex. Her setting was also limited to the elastic impact constitutive equation. Using Yosida type regularization and compactness arguments, she was able to prove the existence of solutions under very weak regularity assumptions. She also discussed uniqueness but proved it only in a very specific case. Further investigation on uniqueness was performed by Percivale in [14] and [15]. He is the first to introduce analyticity hypothesis in this respect. But, his results apply also only to very specific cases. The formulation of the problem with completely inelastic impacts has been extensively studied by Moreau [12]. An existence result was proved by Monteiro-Marques [10] in the particular case in which the configuration space is euclidean R d and the unilateral constraints are described by a single smooth function. Very recently, Schatzman [19] studied the general one degree-offreedom problem with arbitrary impact constitutive law. In this case, she proved uniqueness under analyticity assumption on the data.
None of these results has the generality required by Mechanics. The existence and uniqueness results are proved under assumptions which are obviously not fulfilled in most of discrete mechanical systems which may generally be encountered, except the last result of Schatzman, but it is limited to the one degree-of-freedom problem.
In this paper, the dynamics of discrete mechanical systems with perfect unilateral constraints is formulated in a very general setting. To reach full generality, the configuration space is supposed to be an arbitrary riemannian manifold instead of an euclidean space. However, only the most elementary level of differential geometry is needed. The resulting general evolution problem is studied. The existence and uniqueness of a solution associated with given initial condition is proved provided the data are analytic.
In section 2, we give a precise mathematical definition of what we call discrete mechanical system and system of bilateral constraints. We also recall some basic results connected to these definitions that we shall use subsequently.
In section 3, a formulation of the equations of the dynamics of discrete mechanical systems with perfect unilateral constraints is presented. The content of this section follows very closely the work of Moreau [12]. It is included since Moreau restricts himself to completely inelastic impacts. More generality, including the case of elastic impacts, is obtained here with no supplementary difficulty.
In section 4, we prove a local existence and uniqueness result concerning the general problem of the dynamics of discrete mechanical systems with perfect unilateral constraints, under the single assumption that the data are analytic. Existence and uniqueness of a maximal solution follows immediately. A sufficient condition to ensure that this maximal solution is defined for all time is also presented.
In section 5, three examples are discussed. One is due to Moreau and another one to Schatzman. They are included for sake of completeness. The aim of these examples is to show that the regularity assumptions made in the previous section are, in some sense, minimal.
In section 6, we illustrate the generality of the theorems of section 3 in applying them to simple examples issuing from Mechanics.
In section 7, the continuous dependence of the solution on initial conditions is discussed. Dependence on initial conditions is seen to be not continuous in general. However, a restrictive case where continuity holds is exhibited.
In section 8, the numerical computation of the solution is discussed. Problems arise in connection with non-continuous dependence on initial conditions. However, we recall an algorithm, which was first described by Moreau, and prove its convergence in some restrictive cases.
The main results in this paper were announced in Ballard [3].

Discrete mechanical systems and perfect bilateral constraints
The aim of this section is to give a precise definition of what we call a discrete mechanical system, to introduce notations and to recall some basic results that we shall use later on. For a comprehensive presentation, the reader is referred to Arnold [2] and Abraham & Marsden [1].

Discrete mechanical systems
Definition 1 A discrete mechanical system is: -A Hausdorff, smooth (of class C p with 2 ≤ p ≤ ∞) connected manifold Q of dimension d whose topology has a countable basis. Q is called the configuration space of the discrete mechanical system. d is its number of degrees of freedom. The tangent bundle T Q of Q is called the phase space or the state space. A point q of Q is a configuration of the system and a point of T Q a state of the system. T * Q denotes the cotangent bundle, Π Q : T Q → Q and Π * Q : T * Q → Q the natural projections. T q Q is, as usual, the tangent space at q and, to designate an element v of T Q, we shall often use the redundant notation (q, v) where q = Π Q (v) and v ∈ T q Q. A curve on Q (i.e. a continuous mapping from a real interval I to Q) is also called a motion of the system. If a motion q : I → Q admits a tangent vector at t, it will be denoted by (q(t),q(t)). This notation is an abuse consecrated by tradition. The dot will also be used in general to denote a derivative with respect to time. A local chart on Q is also called a local parametrization of the system. -A riemannian metric on Q denoted by (·, ·) q . The mapping is the kinetic energy of the system.
-A real interval I and a smooth (of class C p with 1 ≤ p ≤ p) mapping f : T Q × I → T * Q such that: The mapping f is called virtual power of internal, external and inertial efforts acting on the system or, in short, the efforts mapping. We will denote by ·, · q the local duality product on T * q Q × T q Q and (and = −1 its inverse) the isomorphism of vector bundles from T Q onto T * Q canonically associated with the riemannian metric on Q.
The Fundamental Principle of Dynamics asserts that any motion of the system is of class C 2 and has to satisfy: where D dt denotes the operator of covariant derivation along q(t) canonically associated with the riemannian metric of Q.
In the sequel, for (U, ψ) a local chart on Q, (e 1 (q), e 2 (q), · · · , e d (q)) and e 1 (q), e 2 (q), · · · , e d (q) will denote the dual basis of T q Q and T * q Q naturally associated with the considered chart. ψ(q), that we shall abusively continue to denote by q, is an element q 1 , q 2 , · · · , q d of R d . If q(t) is a smooth motion on Q, q 1 (t),q 2 (t), · · · ,q d (t) will be the components of its tangent vector (also said velocity) in the local basis:q (t) =q i (t)e i (q(t)), where Einstein's summation convention applies. It will always apply unless explicitly stated. No confusion induced by this notation should be expected since: ∀i ∈ {1, 2, · · · , d} ,q i (t) = d dt q i (t).
In general, we shall use the same notation to denote a function and its representative in a chart. As usual, g ij (q) will denote the covariant components of the metric in the considered chart and g ij (q) its contravariant components. Γ i jk (q) will be the associated Christoffel symbols: Proposition 2 (Lagrange) Let (U, ψ) be a local chart and q(t) a C 2 motion on Q. One has: ,q(t)) e i (q(t)).
Proof. It is straightforward since: Coming back to the equation of motion (2), suppose we are given in supplement an element t 0 of I, said initial instant, and an element (q 0 , v 0 ) of T Q, said initial state. Then, we obtain the following Cauchy problem C on Q: (q(t 0 ),q(t 0 )) = (q 0 , v 0 ) . This result allows us to associate with any discrete mechanical system a dynamical system, that is a two real parameters collection F s,t of mappings from T Q into T Q such that:

Cauchy
To illustrate these basic definitions and results, we give a simple example that we shall reuse later on in a slightly different context. Consider a plane system of two homogeneous rigid bars 1 and 2. The bar 1, of length l 1 and mass m 1 is connected to a fixed support by means of a perfect ball-and-socket joint equipped with a spiral spring of stiffness k 1 . The bar 2, of length l 2 and mass m 2 is connected to the free extremity of the bar 1 by means of another ball-and-socket joint also equipped with a spiral spring of stiffness k 2 . A force acts on the free extremity of the bar 2. This force remains parallel to the direction of the bar 2 and is of constant magnitude λ > 0 (see figure 2.1). With this system is associated the following discrete mechanical system:  -The configuration space is R 2 equipped with its canonical structure of C ∞ manifold (it is not the 2-torus since the spiral springs impose to be able to count the 'number of turns'). This manifold may be represented by a single chart; in other terms, there exists a global parametrization of the system. In the sequel, we shall only use the chart (q 1 , q 2 ) defined by the angular measures associated with each of the joints. -The kinetic energy is: This kinetic energy defines a riemannian structure on the configuration space. The expression of the metric tensor in the considered chart is: -The efforts mapping has for expression in the considered chart: Proposition 2 allows us to form easily the equation of motion in the considered chart: (4) The deterministic conclusion of Cauchy-Lipshitz theorem on the dynamic evolution of the system is illusive. Indeed, if we add to the differential system (4) the initial condition: it is easily seen that the maximal solution is the identically vanishing function on the real line. But, Poincaré-Lyapunov theory shows that this solution is unstable for some value of λ and the real motion will differ in this case from this trivial solution.
The correct analysis of the motion should in this case refer to some investigation of topological nature on the dynamical system generated by the equation of motion.
In any case, one has to abandon the objective of predicting exactly the motion of the system. One has to accept to obtain only partial information on this motion: this is a consequence of the over-idealization made during the modelling process. However, Cauchy-Lipschitz theorem is at the basis of any further analysis which has to be performed on the equation of motion. This fact will be discussed with more details in section 7 in the context of the dynamics of discrete mechanical systems with perfect unilateral constraints.

Bilateral constraints
One may introduce on discrete mechanical system another type of efforts, not taken into account by the efforts mapping f . Indeed, one may specify some efforts by their kinematical effects: one speaks of constraint. A constraint induces a restriction on the admissible motions of the system which is expressed by means of a finite number n of smooth real functions defined on Q: The word constraint in the singular will be used indifferently to speak either of a constraint specifically associated with a single function ϕ i or of the constraint associated with all the functions ϕ i . In this terminology, a set of constraints is still a constraint. In formula (5), the constraint is said holonomic (because it applies on the configuration and not on the state), scleronomic (because it does not depend explicitly on time) and bilateral (because it is expressed only by equalities and not inequalities). We denote by S the following subset of Q: and we add the assumption that the functions ϕ i are functionally independent: for all q in S, the dϕ i (q) (i ∈ {1, 2, · · · , n}) are linearly independent in T * Q. As a consequence, S is a submanifold of Q of dimension d − n. The realization of kinematical specifications (5) necessarily involves a virtual power of reaction efforts mapping R taking values in T * Q. It is a priori unknown. Now, consider an initial instant t 0 in I and an initial state (q 0 , v 0 ) compatible with the constraint (i.e. (q 0 , v 0 ) ∈ T S ⊂ T Q). The evolution problem associated with the discrete mechanical system with bilateral constraint is: find These equations fail to determine the motion of the system: one has to precise the mapping R by means of a phenomenological assumption on the way the constraint acts. A constraint will be said perfect if the associated reaction efforts do not produce work in any virtual velocity compatible with the constraint: As a result: Therefore, if the bilateral constraint is perfect, the evolution problem may be written as: find dt is the operator of covariant derivation on Q. Let q be a point of Q, v a vector of T q Q, and E a subspace of T q Q. The orthogonal projection of v on E for the scalar product of T Q induced by the riemannian structure of Q will be denoted by Proj q [v; E]. Similarly, Proj * q [v * ; E * ] will denote the orthogonal projection of the 1-form v * on the subspace E * of T * q Q. Then, consider the evolution problem E S : find T > t 0 and q ∈ C 2 ([t 0 , T [; S) such that: where T * q S is considered as a subspace of T * q Q and D S dt is the operator of covariant derivation on S equipped with the riemannian structure inherited from that of Q. We have: Proof. First, let us identify T q S and T * q S to subspaces of T q Q and T * q Q. One has T * q S = (T q S). T * q S and n i=1 R dϕ i (q) are complementary orthogonal subspaces of T * q Q and (Chavel [7], p. 54): and, which show that q is a solution of E S . Reciprocally, let q be a solution of E S . From: we deduce the existence of n functions λ i : [t 0 , T [→ R such that: It follows that: where the Gram matrix is invertible because of the assumption on the functions ϕ i . This shows that the function λ i are uniquely determined and that they are continuous. Therefore, q generates a solution of E Q . The second part of Proposition 3 follows from standart results on ordinary differential equations (see for example Coddington & Levinson [8]).
The moral of Proposition 3 is that adding a perfect bilateral constraint to a discrete mechanical system generates another discrete mechanical system with smaller number of degrees of freedom.

Discrete mechanical system with perfect unilateral constraints
This section deals with the formulation of the equation of motion of a discrete mechanical system when some perfect unilateral constraints are added. All the basic ideas of this section are due to Moreau [12]. It is included since Moreau restricts himself to the special case of completely inelastic impacts and also because Moreau does not consider the general case of an arbitrary configuration manifold equipped with an arbitrary riemannian structure.

Kinematical setting
Consider a discrete mechanical system according to section 2.1 and suppose that a finite number n of unilateral constraints are taken into account: where the ϕ i : Q → R are C 1 functions. The closed subset A of Q defined by: is called the admissible configuration set. We define the mapping J by: where P ({1, 2, · · · , n}) denotes the set of all subsets of {1, 2, · · · , n}. The set J(q) is called the set of all active constraints in the configuration q. As in the case of bilateral constraints, a functionally independence assumption is made on the functions ϕ i : As an easy consequence of the regularity assumptions made on the functions ϕ i , the boundary ∂A and the interior Consider a motion in A (i.e. a continuous mapping from a real interval I to A) and assume that a right velocityq + (t) ∈ T q(t) Q exists for all instant t of I. We necessarily have: ∀i ∈ {1, 2, · · · , n} , ∀t ∈ I, ϕ i (q(t)) = 0 =⇒ dϕ i (q(t)),q + (t) q(t) ≤ 0 or, equivalently, where ∇ϕ i (q) is the gradient of ϕ i at q defined by: Thus, if the system has configuration q, then the right velocityq + is necessarily in the closed convex cone V (q) of T q Q defined by: is called the cone of admissible right velocities at the configuration q. In particular, Similarly, if a left velocityq − ∈ T q Q exists, then,

Equation of motion
As for bilateral constraints, the realization of the constraints induces some reaction effort R. The following hypothesis are made: -H1: the unilateral constraints are of type contact without adhesion: -H2: the unilateral constraints are perfect: There results from hypothesis H1 and H2 and Farkas' lemma (see e.g. Rockafellar [16], p. 200) that: Thus, the reaction effort R ∈ T * Q must be such that: (10) N * (q) is a closed convex cone of T * q Q and it is the polar cone of V (q) in the duality T q Q, T * q Q . We will also have to consider the polar cone N (q) of V (q) for the euclidean structure of T q Q: Now, consider a motion q(t) starting at q 0 ∈ • A at time t 0 with velocity v 0 . Assumed to be continuous, q(t) remains in • A on a right neighbourhood of t 0 . By formula (10), the reaction effort R vanishes as long as q(t) is in • A and the motion is governed by the ordinary differential equation: Suppose that the solution of this Cauchy problem meets ∂A at some instant greater than t 0 . Denote by T the smallest of these instants. The motion admits a left velocity vector v − T at time T . Of course, there may happen: v − T ∈ V (q(T )). In this case, no differentiable prolongation of the motion can exist in A for t greater than T . The requirement of differentiability has to be dropped. An instant such T is called an instant of impact. However, we are going to still require the existence of a right velocity vectorq + (t) ∈ V (q(t)) at every instant t. The right velocity need not to be a continuous function of time and the equation of motion should be understood in sense of Schwartz's distribution. Actually, we require R to be a vector valued measure rather than a general distribution. We denote by M M A(I; Q) (motions with measure acceleration) the set of all absolutely continuous motions q(t) from a real interval I to Q admitting a right velocityq + (t) at every instant t of I and such that the functionq + (t) has locally bounded variation over I. Naturally, bounded variation is classically defined only for functions taking values in a normed vector space. However, for any absolutely continuous curve q(t) on a riemannian manifold, parallel translation along q(t) classically provides intrinsic identification of the tangent spaces at different points of the curve and so, the definitions can easily be carried over to this case. The precise mathematical setting is postponed to the appendix. The reader will notice from the appendix that with any motion q ∈ M M A(I; Q) is intrinsically associated the covariant Stieljes measure Dq + of its right velocityq + . The equation of motion takes the form: where dt denotes the Lebesgue measure. We have to give a precise meaning to condition (10) with R being a vector valued measure. By convention, we shall write: to mean: if θ ∈ L 1 loc (I, q, |R| ; T * Q) is the density of measure R with respect to its modulus measure |R| defined by proposition 25 of the appendix, then, one has: for |R| -a.a. t ∈ I.
This requirement is easily seen to be equivalent to the requirement of the existence of n nonpositive real measures λ i such that: Using this convention, the final form of the equation of motion is:

The impact constitutive equation
We begin this section by an example. Consider the one degree-of-freedom mechanical system whose configuration space is R equipped with its canonical euclidean structure. The efforts mapping f vanishes identically and the unilateral constraint is represented by the single function ϕ 1 (q) = q so that the admissible configuration set A is R − . At initial time t 0 = 0, we consider an initial state (q 0 , v 0 ) such that q 0 < 0 and v 0 > 0. It is readily seen from the equation of motion (13) that an impact necessarily occurs at time t = −q 0 /v 0 . At this time, the left velocity is v 0 . But, the right velocity can take any negative value and whatever it is, it is compatible with the equation of motion.
The reason for this indetermination lies in the phenomenological nature of the interaction of the system with the obstacle. Thus, we are led to make the general following hypothesis: -H3: the interaction of the system with the obstacle at time t is completely determined by the present configuration q(t) and the present left velocityq − (t).
In other terms, we postulate the existence of a mapping F : T Q → T Q describing the interaction of the system with the obstacle during an impact: To ensure compatibility with the equation of motion (13), the mapping F should satisfy: First, consider the particular case of a motion with no more than one active constraint at any time (∀t, CardJ(q(t)) ≤ 1). The normal cone N (q(t)) is either {0} or a half-line and hypothesis H3 is equivalent to postulate the existence of an impact function φ : T Q → R such that: . (16) Equation (16) admits the equivalent form: (17) For the general case where more than one constraint may be active at a time, we recall the following (Moreau [11]): Lemma 4 (Moreau) Let V and N be two closed convex polar cones of a real Hilbert space H. Then, As a consequence, the 'impact constitutive equations' (16) and (17) still make sense and are still equivalent when more than one constraint may be active at a time. Therefore, it is natural to retain only the particular forms (16) and (17) of the general impact constitutive equation (14). As a result of this further hypothesis, the phenomenology of the interation of the system with the obstacle during an impact is described by the single impact function φ : T Q → R. The impact function is also often called 'restitution coefficient'. Naturally, the impact function φ cannot be arbitrary and has to satisfy some consistency conditions. For example, the normality condition in (15) requires: But, this is not enough, we have to impose supplementary conditions on φ in order to ensure: With respect to this, we have: Proposition 5 Let V and N be two closed convex polar cones of a real Hilbert Proof. For the 'if' part, suppose φ ≥ 0. By lemma 4, one gets: For the 'only if' part, we have by hypothesis, Evaluating the scalar product with Proj[v − ; N ] and using lemma 4, one gets: H ≤ 0, and therefore the desired conclusion φ ≥ 0.
There results from proposition 5, that we have to require that the impact function φ should be nonnegative. This consistency assumption ensures that conditions (15) and (18) will automatically be fulfilled.
At this stage, it should be underlined that hypothesis H3 implies the general forms (16) or (17) for the impact constitutive equation only in the restrictive case where only at most one constraint is active at a time. In case of multiple impacts, the choice we made is only motivated by aesthetic considerations and also to fix ideas, since the concept of restitution coefficient is so firmly anchored in minds. We shall discuss more completely the relevance of that choice in section 6.4. Now, let us examine another example. Consider the one degree of freedom discrete mechanical system whose configuration space is R equipped with its canonical structure of riemannian manifold. The efforts mapping is supposed to be constant: f (q,q; t) ≡ 2. To this discrete mechanical system, we add the unilateral constraint described by the single function ϕ 1 (q) = q. Thus, A = R − . The impact constitutive equation is given by formula (16) where the impact function is supposed to be the constant 1/2: φ ≡ 1/2. This mechanical system is a formal description of the physical occurence of a single particle subjected to gravity and bouncing on the floor. Consider the initial instant t 0 = 0 and the initial state (q 0 , v 0 ) = (−1, 0). It is readily seen that the function q : , satisfies the equation of motion (13) and also the impact constitutive equation (16). Note, by the way, that this motion exhibits an infinite number of impacts on a compact time subinterval. It could easily be proved that no motion, defined on [0, ∞[, with finite number of impact on every compact interval can exist. Now, we are going to analyse what happens when the flow of time is reversed. Let define q by: Considering the initial state (q 0 , v 0 ) = (0, 0) at t 0 = 0, it is easily seen that q satisfies both the equation of motion and the impact constitutive equation as soon as the impact function is replaced by φ ≡ 2. But, q ≡ 0 is also seen to satisfy the same initial condition, the equation of motion and the impact constitutive equation.
To eliminate this pathology of non-uniqueness, we are led to add the following hypothesis: -H4: the kinetic energy of the system can not increase during an impact: Taking into account the impact constitutive equation (16), condition (19) can be rewritten as: The final form of the impact constitutive equation is therefore: where the impact function φ is an arbitrary function from T Q to [0, 1]. The two extreme cases φ ≡ 0 and φ ≡ 1 are called respectively the completely inelastic and the elastic impact function.

Formulation of the evolution problem
In this subsection, the results of the previous subsections are brought together in order to formulate the resulting evolution problem which will be studied in the subsequent sections. We add an assumption on the regularity of the data: they are supposed to be real-analytic. This assumption will be motivated by the counterexamples of section 5. The precise mathematical setting is: -Q is an analytic riemannian manifold of dimension d.
ϕ i (i = 1, 2, · · · , n) are n real analytic functions defined on Q. One defines: The functions ϕ i are assumed to be functionally independent in the sense that: -The impact function φ is an arbitrary function from T A − into [0, 1]. No regularity assumption is made on φ.

-I is a real interval and O an open neighbourhood of T A + in T Q and the efforts
mapping is supposed to be an analytic mapping from O×I into T * Q such that: -We are given an initial time t 0 in I such that I contains a right neighboorhood of t 0 and an initial state According to the previous subsections, the evolution problem associated with the dynamics of discrete mechanical systems with perfect unilateral constraints can be formulated as: where equation (23) is to be understood in the sense of convention (11). The existence and uniqueness of solutions for problem P will be studied in section 4. Before studying this question, let us state two almost obvious results.
Proposition 6 Any solution (T, q) of problem P satisfies: A , q |J is analytic and: Proof.
A . By equality (9), one has: ∀t ∈ J, -We have: There results from proposition 28 thatq + |J is locally absolutely continuous, and, therefore, by proposition 32. We get again by proposition 28. The conclusion follows by use of classical results on ordinary differential equations.
Proposition 7 (Energy inequality) Any solution (T, q) of problem P satisfies: Proof. We have the following equality between real measures: Integrating over ]t 1 , t 2 ] and using propositions 30 and 32, one gets: D is (at most) countable and therefore Lebesgue-negligible. There results: The dynamics of discrete mechanical systems with perfect unilateral constraints 19 -Let us denote by θ R the density of measure R with respect to its modulus measure |R| provided by proposition 26. Since: we get: and therefore the second integral in equation (26) is nonnegative whereas the third is nonpositive since V (q(t) and N * (q(t)) are polar cones. As a consequence: is nonpositive by virtue of hypothesis H4.
The proposition results from equation (25) and from the estimation of these four integrals.

Existence and uniqueness of solutions for problem P
This section is devoted to prove existence and uniqueness of a maximal solution for problem P. Sufficient conditions to ensure that this maximal solution is defined for all time are also given. More precisely, we are going to prove the following results.
Then, a standart argument yields: denotes an arbitrary solution of problem P, then: We shall say that the maximal solution of problem P is global if it is defined on

Theorem 10
Assume that the configuration space Q is a complete riemannian manifold and that the efforts mapping f admits the estimate: where l(t) is a (necessarily nonnegative) function of L 1 loc (R; R). Then, the maximal solution of problem P is global.
Let us say a word of how the proof of these results is going to be structured. First, we construct T a > t 0 and an analytic function q a : is any other solution, then q and q a coincide identically on a right neighbourhood of t 0 . This is the most difficult part to prove but it is also the crucial one. For the proof of theorem 10, we first notice that for q ∈ M M A([t 0 , T [; Q) (T finite) satisfying the equation of motion (23), boundedness ofq + implies finiteness of Var (q + ; [t 0 , T [): this is the object of proposition 18 of section 4.3. Note that the impact constitutive equation (24) plays no role in this property. Then, theorem 10 is deduced from the energy inequality (proposition 7) and Gronwall-Bellman lemma.
In the proof of these results, we shall use the following notations. If J is any subset of {1, 2, · · · , n}, Gram(J) will be the Gram matrix: If x is an arbitrary element of R J whose components are x i with i ∈ J, then (x i ) i∈J will denote the columm matrix: and T (x i ) i∈J the associated row matrix:

Proof of local existence
Local existence is rather easy to prove in the setting of analytic data. The proof is a little bit lengthy but involves no specific difficulty. We begin by technical lemmas.
Let X(t) be a C ∞ vector field along a C ∞ curve q(t) on Q. The covariant derivative D dt X(t) of X along q defines a C ∞ vector field along q. So, one may consider its covariant derivative along q which will be denoted by D 2 dt 2 X(t). By induction, we get the definition of D i dt i X(t) (i ∈ N * ). One has: Lemma 11 Let X be a C ∞ vector field on Q and q I , q II two C ∞ curves on Q. m being a nonnegative integer, one assumes: Then, Then, which gives the desired conclusion for the case m = 1. For arbitrary m, an easy induction based on the same type of computation in a local chart shows the existence of functions h i : (T Q) i−1 → T Q independent on the considered curve q(t) and such that: Exactly the same technique applies to prove where I denotes a real interval containing t 0 . Let m be an arbitrary nonnegative integer and q I , q II two C ∞ curves on Q such that: and: Then, We denote by q u and q c some local solutions of problems: furnished respectively by Cauchy-Lipschitz theorem and proposition 3. Then: then: Proof. First, from: on one hand, and: on the other hand. Therefore: which is the announced result. Second, assume: An easy induction based on lemmas 11 and 12 gives: Therefore: Proposition 14 Considering the data of problem P, we denote by P the following evolution problem. Problem P : find T ∈ I (T > t 0 ), an analytic curve q : [t 0 , T [→ Q and n analytic functions Then, problem P admits a solution (T, q, λ 1 , · · · , λ n ) unique in the sense that any other solution is either a restriction or an analytic extension of (T, q, λ 1 , · · · , λ n ).
Proof. First, let us precise once for all that the meaning of an analytic function on a non necessarily open set S is that there is an analytic extension to an open set O containing S.
Step 1. Construction of some functions q and λ i . Define: and I 0 = K 0 = ∅. We denote by q (1) a solution of the cauchy problem: Define: (1) be the solution of the variational inequality: furnished by Lions-Stampacchia theorem (see [9]). Let μ and I 1 , J 1 , K 1 by: , I n , J n and K n are constructed. Then, q (n+1) is defined to be a local solution of the Cauchy problem: is defined to be the solution of the variational inequality: and I n+1 , J n+1 , K n+1 by: , I n , J n and K n are defined by induction for n ∈ N * and for all n in N * , I n , J n , K n is a partition of J 0 . Moreover, one has: ∀n ∈ N, Define: It is readily seen that I, J, K form a partition of J 0 . We denote by (q, (λ i ) i∈I ) a local solution of the evolution problem: furnished by proposition 3. The functions q and λ i are analytic. For any i in {1, 2, · · · , n} \ I, the functions λ i are defined to be the identically vanishing function: Step 2. We have: Indeed, applying lemma 13 to Cauchy problems C (1) and C yields, thanks to equation (27), But, by definition of I, and, so, and the conclusion follows for i = 0, since the Gram matrix is positive definite. For i ≥ 1, we only have to apply successively lemma 13 to Cauchy problems C (i+1) and C.
Step 3. The functions q and λ i define a solution of problem P . Indeed, by construction of the real numbers λ (j) i and μ (j) i and by step 2, we have: Each function λ i (t) and ϕ i (q(t)) being real-analytic, there results: Actually, α > 0 is assumed to be sufficiently small to ensure: which is possible simply by continuity. Now, it is easily seen that t 0 + α, q, (λ i ) i∈{1,2,···,n} defines a solution of problem P .
Step 4. Uniqueness part of the proposition. By Cauchy-Lipshitz theorem, q is uniquely determined by the functions λ j (j = 1, 2, · · · , n). Being analytic, these functions λ j are uniquely determined by the collection of real numbers . Therefore, to prove uniqueness, one has only to show that these real numbers are determined by the data of the evolution problem.
Consider an arbitrary analytic solution (T, q, λ 1 , · · · , λ n ) of problem P . A repeated use of lemma 13, similar to the one of step 2 yields: and the conclusion follows.
Proof of the local existence part of theorem 8. Let (T a , q a , λ 1 a , · · · , λ n a ) be an analytic solution of problem P . It is readily seen that (T a , q a ) is a local solution of problem P.

Proof of local uniqueness
Local uniqueness is the most difficult part of theorem 8. First, we recall a standart result: Lemma 15 (Gronwall-Bellman) Consider two functions m 1 ∈ BV ([0, T ]; R) and m 2 ∈ L 1 (0, T ; R) such that: then, Proof. This is almost obvious. Dividing each member of the inequality by t m+2 , we obtain: After integration, Gronwall-Bellman lemma yields: Then, an integration by part gives the desired conclusion.
Proof of the local uniqueness part of theorem 8. Consider, on one hand, the analytic solution (T a , q a , λ 1 a , · · · , λ n a ) of problem P supplied by proposition 14, and on the other hand, an arbitrary solution (T, q) of problem P. We have to prove that q and q a identically coincide on a right neighbourhood of t 0 .
Step 1. Parametrization of the problem and notations.
Consider a local chart ψ : U ⊂ Q → R d on Q centered at q 0 such that the cardJ(q 0 ) first components of ψ(q) are (ϕ i (q)) i∈J(q 0 ) . Recall that such a chart exists since (dϕ i (q 0 )) i∈J(q 0 ) is linearly independent in T * q 0 Q. We choose α > 0, sufficiently small to have: (28) Such a choice for α is possible because: the functions q a (t) and ϕ i (q a (t)) are real analytic, the functions q(t) and ϕ i (q(t)) are continuous.
We denote by f i the components of f in the natural basis (e i ) associated with the chart under consideration. Since q a is an analytic local solution of problem P, we have: after appropriate renumbering of the functions λ i a . In the sequel, d 0 will stand for cardJ(q 0 ). There results from these choices that: λ i a ≡ 0. We denote by |.| the standart euclidean norm on R d . Confusing (abusively) q and ψ(q), we shall write: and: Step 2. There exists some positive real constants C 1 and C 2 such that the following estimate: holds.
To prove this assertion, we first write the equation of motion (23) in the chart under consideration using proposition 29: where the μ j are nonpositive real measures. But, by propositions 29 and 30, Therefore, by formulae (28), and, which is a nonpositive real measure by proposition 7. Therefore, in the sense of ordering of real measures. Integrating over ]t 0 , t] (t ∈ [t 0 , t 0 + α]), we get: The term within the integral sign is an analytic function of the three variables q, q + and s. Therefore, it is also an analytic function of the three variables q − q a , q + −q a and s. It is written under the form: But, each function F i can be decomposed under the form: where the G i are analytic and G i (0, 0; s) ≡ 0. Hence, there exists d positive constants M i such that: Defining M to be the maximum of the constants M i , we have proved: Moreover, by a compactness argument: and therefore: Moreover, by use of Cauchy-Schwartz inequality: We obtain: where formulae (29) have been used. We define: Notice that, actually: and, so, by the analyticity of functions q i a and λ i a : Multiplying both terms of inequality (31) by e −C 2 t and integrating, we get: which is nothing but estimate (30). Step Indeed, by estimate (30): which is, after integration by part: the two members of inequality (32) are nonnegative and, therefore, the inequality is preserved when taking the absolute value of each member. We get: The dynamics of discrete mechanical systems with perfect unilateral constraints 33 We define: . With these notations, we obtain: where the L i are nonnegative real-analytic functions and the Q i are nonnegative continuous functions which all vanish at t = 0 and which are differentiable at the origin. We are going to prove that inequality (33) implies that: The functions L i being nonnegative real-analytic, there exists nonnegative integers n 1 < n 2 < · · · < n m , a partition I 1 , I 2 , · · · , I m of {1, 2, · · · , d 0 }, and nonnegative real-analytic functions G i such that: Inequality (33) may be rewritten as: But, by the analyticity of the functions G i : Therefore, Integrating by part the left member of the inequality, we obtain: Since each function G i (σ)Q i (σ)/σ is bounded over [0, β], there exists a nonnegative real constant H such that: Inequality (34) gives: where H 1 is a non negative real constant. As a consequence of lemma 16, we obtain: Coming back to inequality (34), we get: Applying once more lemma 16, we obtain: Proceeding inductively, we obtain: The dynamics of discrete mechanical systems with perfect unilateral constraints 35 But, by inequality (34), Using a last time lemma 16, we get: which implies: which is nothing but: λ i a (σ)q +i (σ) = 0, and the assertion of step 3 is proved.
Step 4. Conclusion of the proof of local uniqueness.
Bringing together the results of steps 2 and 3, we get: which gives the desired conclusion: ∀t ∈ [t 0 , t 0 + β], q(t) = q a (t).

Global solutions: proof of theorem 10
First, we recall a classical lemma whose proof may be found for example in [5], p. 157.

Proposition 18
The riemannian manifold Q is assumed to be complete. Let (T, q) be a solution of problem P such that: -T ∈ • I (and ,in particular, T = +∞), -q + (t) q(t) is bounded: thenq + has bounded variation over [t 0 , T [: Proof. We denote by d the distance function associated with the metric space Q. Since, Let (U, ψ) a local chart at q T on Q such that the cardJ(q T ) first components of ψ(q) in R d are (ϕ i (q)) i∈J(q T ) . Consider a compact neighbourhood K of q T in Q such that: One defines: Since [t 0 , t 0 ] is compact, one has: therefore, it remains only to prove: λ max (resp. λ min ) will denote the maximum (resp. the minimum) of the greatest (resp. least) eigenvalue of the matrix (g ij (q)) i,j=1,2,···,d when q wanders in K.
With these notations, one obtains immediately: The dynamics of discrete mechanical systems with perfect unilateral constraints 37 We denote by B q (0, V m ) the closed ball of T q Q with radius V m and centered at the origin. Considering the following compact subset K of T Q: we define the following nonnegative real constant: and: Writing inclusion (23) in the local chart (U, ψ), we obtain: where the λ i are d nonpositive real measures on ]t 0 ; T [. Expressing the Christoffel symbols in terms of the metric, one has: One deduces: f (q(s),q + (s); s),q + (s) q(s) ds. Thus, By lemma 17, we obtain: which gives, using the hypothesis of the theorem: (1 + l(s)) d(q(s), q 0 ) + q + (s) q(s) ds.
By Gronwall-Bellman lemma (lemma 15), one gets: (1+l(s)) ds , which shows that the function t → q + (t) q(t) is bounded over [t 0 , T [. By the completeness of Q, we deduce, on one hand that exists in Q and, on the other hand, that thanks to proposition 18. Thus,

Three counterexamples
The existence and uniqueness of solution for problem P has been proved under the assumption of functional independence for the constraint and of analyticity for the data. The three examples which are developed in this section aim at showing that these assumptions cannot be weakened very much. In example 1, we show that, in case that the functional independence of the constraints does not hold, the existence of solution may be lost. For the question of the regularity assumptions on the data, the existence of solution can be proved with much weaker assumptions. However, the uniqueness of solutions is generally lost in such a case as seen in examples 2 and 3. In these examples, the data are supposed to have only regularity C ∞ and two different solutions can be exhibited.
Example 1 is extracted from Moreau [12] and example 2 is due to Schatzman [18], but an earlier counterexample in the same spirit is also to be found in Bressan [4].

Example 1
Consider a discrete mechanical system whose configuration space is euclidean R 3 . The unilateral contraints are kinematically described by the three following functions (n = 3): where q = (q 1 , q 2 , q 3 ) ∈ R 3 . The initial instant is supposed to be t 0 = 0 and the initial state is given by q 0 = (0, 0, 0) and v 0 = (0, 2, −1). It follows that: It is readily seen that v 0 belongs to V (q 0 ).
Let now α > 0 be an arbitrary positive real number. Any motion q(t) in M M A([0, α[; R 3 ) compatible with this initial data may be written as: which can not be compatible with: We deduce that no motion in M M A([0, α[; R 3 ) can be compatible with this initial data whatever α > 0 is. Note that in this particular case: dϕ 1 (q 0 ) = −dϕ 2 (q 0 ) and the unilateral constraints are not functionally independent.

Example 2
Consider a discrete mechanical system whose configuration space is R equipped with its canonical structure of riemannian manifold. This is the configuration space of a particle with unit mass constrained to move along a line. A fixed obstacle at the origin is taken into consideration. It gives rise to a unilateral constraint kinematically described by the single function (n = 1): Therefore, the admissible configuration set is A = R − . It is assumed that the impact constitutive equation is the elastic one: φ (q,q − ) ≡ 1 and that the efforts mapping f does not depend on the state but only on time. It will be denoted by f (t). The initial instant is t 0 = 0 and the initial state is (q 0 , v 0 ) = (0, 0). Denoting by RCLBV (I; R) the space of right continuous functions with locally bounded variation from a real interval I to R, problem P admits here the equivalent formulation: find T > 0 and v ∈ RCLBV ([0, T [; R) such that: dt is a nonpositive real measure such that: We investigate uniqueness under the assumption that f is of class C ∞ . Suppose, in addition, that f is nonnegative: It is readily seen that the null function v ≡ 0 on R + is a solution of problem P whatever is the nonnegative C ∞ function f . Now, we are going to construct an explicit example of such a function f in such a way that the associated problem P admits another solution, different from the identically vanishing one. First, let define a function ρ by: The last assertion comes from the fact that: Consider also the real convergent serie: (n + 5) 2 (n + 1)(n + 2)(n + 3)(n + 4) n∈N We define: (n + 5) 2 (n + 1)(n + 2)(n + 3)(n + 4) > 0 Clearly, a 0 = T and the sequence (a n ) n∈N decreases strictly and converges towards 0 when n tends toward infinity. Actually, a n ∼ 1 n when n → +∞ (40) by a very classical and elementary argument. We denote by (δ n ) n∈N , (f n ) n∈N , (v n ) n∈N the real sequences defined by: δ n = n + 5 (n + 1)(n + 2)(n + 4) (i.e. δ n = n + 3 n + 5 (a n − a n+1 ) < a n − a n+1 ), , and : Proof. The only thing which is not obvious is that f is C ∞ at 0. Since: 1] |ρ(s)| , then, lim t→0 + f (t) = 0 and f is continuous at 0. Now, we are going to prove: which will imply by an easy induction that f ∈ C ∞ ([0, T [; R) and: Let fix an arbitrary r in N, we have: Therefore, to prove (43), it suffices to verify: lim n→∞ f n (n + 5) r a n+1 (a n − a n+1 ) r = 0, but, by estimate (40), one has: Second, we claim that: Proof. It is clear that v is continuous on each interval ]a n+1 , a n [ and right contin- = − 1 (n + 4)! + 1 n! n + 5 (n + 1)(n + 2)(n + 3)(n + 4) Since v is nondecresing on each interval [a n+1 , a n [: Denoting δ t the dirac measure located at t, one has: which is a (bounded) nonpositive measure whose support is {0} ∪ {a n ; n ∈ N * }.
Third, we claim that: if q is defined by: then: Proof. An easy calculation using last assertion of formulae (39) shows: We deduce, that making the choice described by relations (41) for the function f , then the function v defined by relations (42) is a solution of the corresponding problem P whereas the identically vanishing function is also a solution. Therefore, the uniqueness of solution does not hold in general if f and the functions ϕ i are supposed to be of class C ∞ only.

Example 3
In example 2, we considered a particle at rest at the initial instant and in contact with the obstacle. Then, a force acts on the particle, constantly pushing it against the obstacle (f ≥ 0). For the particular choice of the function f we made, immobility is a possible motion whereas a bouncing motion is also possible. It is intuitively clear that the assumed elastic impact constitutive equation plays a central role in such a phenomenon. The question arises to know whether such a pathology is possible with the completely inelastic impact constitutive equation Sticking to notations of example 2, the evolution problem takes in this case the form: find T > 0 and v ∈ RCLBV([0, T [; R) such that: dt is a nonpositive real measure such that: If we still assume in this case that f is nonnegative, then it is easy to see that the only possible motion is immobility. Indeed, if, ∃t 2 , q(t 2 ) < 0, define t 1 = inf {t ∈ R + ; ∀s ∈]t, t 2 ] q(s) < 0}. Then, by continuity of q: t 1 < t 2 and q(t 1 ) = 0. By the completely inelastic impact constitutive equation, one gets: v(t 1 ) = 0, and, so: Nevertheless, we are going to construct an example similar to example 2, which shows that even in case of the completely inelastic impact constitutive equation and f of class C ∞ , we can obtain multiple solutions for the corresponding problem P. Of course, f should not be of constant sign.
The function f is searched under the form: where n ∈ N * . (f 1,n ) n∈N * , (f 2,n ) n∈N * , (δ 1,n ) n∈N * and (δ 2,n ) n∈N * are positive real sequences which are to be determined. We demand: and δ 2,n ≤ 1 2 These sequences are to be determined in such a way that the corresponding problem P admits two distinct solutions v I and v II . We demand that v I , v II and the corresponding functions q I , q II are such that : if n is even, where (q n ) n∈N * and (v n ) n∈N * are positive real sequences which are to be determined.

Illustrative examples and comments
6.1. Punctual particle subjected to gravity and bouncing on the floor. Accumulation of impacts.
Let us come back to the example of section 3.3. The configuration space is R equipped with its canonical structure of riemannian manifold, the unilateral constraint is described by the single function ϕ 1 (q) = q (which gives A = R − ). The efforts mapping is supposed to be constant: f (q,q; t) ≡ 2 and the impact function (restitution coefficient) is the constant 1/2: φ ≡ 1/2. Considering the initial instant t 0 = 0 and the initial state (q 0 , v 0 ) = (−1, 0), we have seen in section 3.3 that the function q : R + → R − defined by: (∀n ∈ N) belongs to M M A(R + ; R − ) and is readily seen to be the maximal solution, according to corollary 9, of the corresponding problem P. The solution q(t) is represented on figure 6.1. It is seen that infinitely many impacts accumulate in any left neighbourhood of instant t = 3. However, as predicted by corollary 9, for each instant t ∈ R + , there exists a right neighbourhood [t, t + η[ of t, such that the restriction of q to [t, t + η[ is analytic. A straightforward and general consequence of this is the following.

Proposition 19
Let q be the maximal solution of problem P furnished by corollary 9. Although infinitely many impacts can accumulate at the left of a given instant, this phenomenon can never occur at the right of any instant. Morover, in the particular case where the impact constitutive equation is the elastic one (φ ≡ 1), the instants of impact are isolated and therefore in finite number in any compact interval of time.
Proof. Since for each instant t ∈ [t 0 , T [, there exists a right neighbourhood [t, t+η[ of t, such that the restriction of q to [t, t + η[ is analytic, we get the first part of the proposition. For the second part, let τ be an arbitrary instant in ]t 0 , T [ and consider the problem P associated with the initial condition (q(τ ), −q − (τ )), the elastic constitutive impact equation and the effort mapping g(q, v; t) defined by: which is analytic. By theorem 8, there exists an analytic function q a : [0, T a [→ Q which is a solution of this problem P. Another solution of problem P coincides with q a on a right neighbourhood of t = 0. Actually, as seen in the proof of local uniqueness (section 4.2), a little bit more is proved: any function q ∈ M M A([0, T [; Q) satisfying the initial condition (21), the unilateral constraint (22), the equation of motion (23) and the energy inequality (proposition 7) has to coincide with q a on a right neighbourhood of t = 0. But, it is readily seen that the function The fact that infinitely many impacts can accumulate at the left of a given instant but not at the right is a specific feature of the analytical setting that is lost in the C ∞ setting as seen in counterexamples 2 and 3. Actually, these counterexamples show that pathologies of nonuniqueness in the C ∞ setting are intimately connected to the possibility of right accumulations of impacts. The fact that the analytical setting prevents from such right accumulations is the thorough reason why we could prove uniqueness in this case.

The double pendulum
In this section, we come back to the double pendulum described in section 2.1 but we add to the system a rigid obstacle on the vertical coordinate axis as represented on figure 6.2. This obstacle may be represented by two analytic functions whose expressions in the global chart of Q described in section 2.1 is: It is readily seen that, except in the particular case where l 1 = l 2 , these constraints are functionally independent: (dϕ i (q)) i∈J(q) is linear independent in T * q Q These unilateral constraints are assumed to be perfect and we consider an impact function φ supposed to be constant on T A − : The constant φ is often called restitution coefficient (of normal velocities). We recall that the particular cases φ = 0 and φ = 1 describe the completely inelastic and the elastic impact constitutive equations. An initial state (q 0 , v 0 ) ∈ T A + is given at time t 0 = 0. This initial state is represented in the considered chart by four real numbers (q 1 0 , q 2 0 ; v 1 0 , v 2 0 ). According to section 3, the motion of the system is governed by the evolution problem: • (q(0),q + (0)) = (q 0 , v 0 ), • ∀t ∈ [0, T [, (q(t),q + (t)) ∈ T A + , where the riemannian structure on Q and the mapping f are those described in section 2.1. Corollary 9 ensures existence and uniqueness of a maximal solution. Now, we are going to check that assuptions of theorem 10 are satisfied so that the maximal solution is defined all over R + . First, Q is a complete riemannian manifold since the quotient topology on the torus T 2 derives from a riemannian structure and T 2 is compact and therefore complete. Second, we have the estimate: where, α = 1 9 m 1 m 2 l 2 1 l 2 2 + 1 12 m 2 2 l 2 1 l 2 2 1 where λ min (q) is the least eigenvalue of the matrix (g ij (q)) i,j=1,2 . But: we get: which achieves the proof of estimate (51). Now, let q I , q II be two points of Q represented by their components in the considered chart (q 1 I , q 2 I ) and (q 1 II , q 2 II ). Q being complete, there is a geodesic g : [s 1 , s 2 ] → Q of minimal length between them. One has: Moreover, recalling: one has: Therefore, By virtue of theorem 10, the motion of the system is defined for all t ∈ R + .

Boltzmann's gas
Consider a collection of N rigid homogeneous balls of mass m and radius R in a rigid rectangular box. The balls cannot interpenetrate. The balls are free of internal or external forces except for the reaction efforts induced by the unilateral constraints. The impact constitutive equation is supposed to be the elastic one. Such a system was introduced by Boltzmann to model the interactions between molecules in a gas in order to perform a statistical analysis to connect the microscopical and macroscopical point of view.
Let us describe the discrete mechanical system associated with this situation. The configuration space is R 3N . Indeed, any configuration is described by the coordinates of the center of the balls in the three-dimensional ambiant space equipped with an origin. Strictly speaking, the configuration space should be R 3N × (SO3) N to incorporate the possible rotations of the balls. But, in this case, it would be readily seen that the rotation velocity of any ball in any motion of the system keeps its value at the initial instant. Therefore, rotations play no role in the motion of the system and we may consider only the restricted configuration space R 3N equipped with its canonical riemannian structure. The forces mapping vanishes identically f (q,q + ; t) ≡ 0. There are N (N + 11)/2 functions ϕ i , since N (N − 1)/2 of them are necessary to express the noninterpenetration constraints: and 6N of them are necessary to express that the balls remains inside the box: where 2a, 2b and 2c are the lengths of the sides of the box. The functions ϕ i are defined by arbitrary numbering. They are easily seen to be analytic and functionally independent. Adding the elastic impact constitutive equation φ(q,q − ) ≡ 1, and an initial condition at time t 0 = 0, the corresponding evolution problem turns out to belong to the class of problem P formulated at the beginning of section 4. Then, corollary 9 and theorem 10 state that, to any initial condition compatible with the constraints, there corresponds a unique maximal motion and it is defined all over R + . By proposition 19, we may also state that there are at most finitely many impacts on any bounded time interval. As a conclusion, the results developed in this paper allow us to associate a dynamical system with Boltzmann's gas. Related to this question, let us mention Boltzmann's famous ergodic hypothesis. Roughly speaking, Boltzmann postulated that in any motion of the system, time averages can be replaced by space averages. The modern mathematical transcript is: for almost every initial condition in an energy level set of the phase space, the associated phase curve spends an amount of time in every measurable piece of the level set proportional to the measure of that piece. The question to know whether Boltzmann's gas is ergodic or not is still an open question. However, a positive answer was given in 1970 by Sinai ([20]) for a two balls gas in a plane rectangular box. Let us underline that this question makes sense only when we are able to associate a dynamical system with Boltzmann's gas.

Newton's balls and the impact constitutive equation
In section 3.3, it has been derived from the two phenomenological assumptions H3 and H4, that the general constitutive impact equatioṅ should satisfy: In the particular case of a motion no more than one active constraint at any time (∀t, CardJ(q(t)) ≤ 1), it has been seen in section 3.3 that the general impact constitutive equation (52) necessarily takes the form:  The principle of Newton's balls experiment is well-known. It is sketched on figure 6.4a. As a result of this multiple impact experiment, we have the familiar picture drawn of figure 6.4b. But, testing the simple impact constitutive equation (54) (with φ ≡ 1) to predict the issue of the experiment, we get the situation drawn on figure 6.4c.
The question arises to known whether the results of section 4 remain true if we abandon the simple impact constitutive equation (54) and adopt the general impact constitutive equation (52) defined by an arbitrary function F fulfilling requirements (53). Actually, a careful examination of the proofs of section 4 shows that the impact constitutive equation is only used through the energy inequality (proposition 7). Moreover, it is readily seen that proposition 7 still holds when the simple impact constitutive equation (24) is replaced by a general one (equation (52)) provided requirements (53) hold true. As a result, all the results of section 4, and in particular, theorem 8, corollary 9 and theorem 10 remain true if we adopt an arbitrary impact constitutive equation instead of equation 24 in the definition of problem P.
A general impact constitutive equation will be said elastic if the last requirement in (53) is replaced by: It is readily seen that proposition 19 still holds with an arbitrary impact constitutive equation. In particular, for a solution of problem P with an arbitrary elastic impact constitutive equation, the impacts are isolated.

Continuous dependence on initial conditions
The theory developed in the previous sections allows us to replace the analysis of the motion of a collection of rigid bodies subjected to perfect constraints either bilateral or unilateral by the analysis of the motion of a point in a piece of a d-dimensional manifold bounded by analytic hypersurfaces which intersect transversally. With appropriate regularity assumptions on the data, the motion is completely determined by the initial condition.
The picture seems to be fairly good and the generalization of the dynamics of discrete systems with bilateral constraints to the case of unilateral constraints seems to be achieved. However, there remains a big difference between unilateral and bilateral dynamics of discrete systems that we want to underline in this section.
A pleasant feature of a dynamical system generated by the flow of an ordinary differential equation is that it is smooth. More precisely, if F t,t 0 is the mapping which associates the state of the system at time t with an arbitrary initial condition at time t 0 , then the mapping F t,t 0 is a local diffeomorphism. In particular, the state of the system at a given instant t depends in a differentiable way of the state at time t 0 . Of course, this smooth dependence may be stiff. In such a case, a small uncertainty on the initial state will produce a big one on the actual state and the motion of the system may turn out to be quantitatively unpredictable on both physical and numerical point of view for large time. In certain circumstances, the theory of smooth dynamical systems allow us to get some qualitative information on the motion for large time.
As we shall see, the picture is strongly different in the case of the dynamics of discrete systems with perfect unilateral constraint. The theorems of section 4 allow us to define a mapping F t,t 0 similar to the flow generated by an ordinary differential equation. But, the mapping F t,t 0 is not smooth any more, it is not even continuous in general. In other terms, the generated dynamical system does not belong to the large class of topological dynamical systems.
Let us check this assertion on a simple example. Consider as a configuration space R 2 supplied with its canonical structure of riemannian manifold. A configuration is denoted by a pair (x, y). No forces act on the system: f ≡ 0. Consider a unilateral constraint associated with the two functions: and the elastic impact constitutive equation φ ≡ 1. At time t 0 = 0, we consider the following set of initial conditions: A straightworward calculation gives the state of the system for all instant in R + . In particular, for t greater than 1, one gets: It is readily seen on this example that, if t is greater than 1, the mapping F t,0 is not continuous at initial condition (−1, −1, 1, 1) (see figure 7). Coming back to x y Fig. 5. The generated dynamical system is not continuous in general.
the two examples of section 6, such a situation occurs if, during the motion of the double pendulum, the two bars hit the obstacle at the same time. In the case of Boltzmann's gas, the pathology occurs when three balls hit at the same time. Let us underline that if we consider an initial condition such as the one in the above example, the solution of the associated problem P has no physical meaning. In such a case, one has to abandon the will of predicting the motion of the system: this is a consequence of the over-idealization made in the indeformability assumption.
However, in the particular case of the one degree-of-freedom problem, where no multiple impacts are possible, Schatzman [19] proved that continuous dependence on initial conditions hold. In the general case, her result admits the following generalization which is proved along the same lines: Theorem 20 Consider the problem P described in section 3.4. Assume furthermore that the impact function φ is constant. Considering the initial condition (q 0 , v 0 ) ∈ T A + at initial instant t 0 , we denote by (T, q) the corresponding maximal solution of problem P. We make the following hypothesis: (with the convention that the empty set is orthogonal). Let us onsider a sequence (q 0n , v 0n ) of elements of T A + converging towards (q 0 , v 0 ). For all n, we denote by (T n , q n ) the maximal solution of the problem P associated with the initial condition (q 0n , v 0n ) at instant t 0 . Then, 2. q n converges towards q uniformly on every compact subset of [t 0 , T [: Proof. The proof of theorem 20 is divided into five steps. Before describing these steps, let us introduce a few notations. We fix, once for all, an arbitrary τ in [t 0 , T [ and a compact neighbourhood K of the compact subset q([t 0 , τ ]) of Q. We define: The subset K of T Q is compact in T Q. We define also: and, Notice that we have δ > 0.
But, this contradicts the definition of t 1 and achieves the proof of step 1.
Step 2. For n large enough, q n is defined on (an interval containing) the interval [t 0 , min(τ, t 0 + δ)]. Moreover, there exists a subsequence of (q n ), also denoted by (q n ), such that: q n converges uniformly on [t 0 , min(τ, t 0 +δ)] towards a function q lim belonging to Proof of step 2. For all q in K ∩ A, there exists a compact neighbourhood K q of q which is included in the domain U q of a local chart (U q , ψ q ) at q such that: -∀q ∈ U q , the cardJ(q) first components of ψ q (q ) are the ϕ i (q ) (i ∈ J(q)).
Being compact, K ∩ A can be covered by a finite number, say L, of • K q l . We denote by λ max (resp. λ min ) the maximum (resp. the minimum) of the greatest (resp. least) eigenvalue of the matrix (g ij (q)) i,j=1,2,···,d when q wanders in K q l and l in {1, 2, · · · , L}. We define: We pick an integer N 0 , large enough to ensure: By step 1, there results: By a compactness argument, we have: As a consequence, for n larger than N 0 , the interval [t 0 , min(τ, t 0 + δ)] is the disjoint union of a finite number, say N n , of intervals I ni such that: ∀i ∈ {1, 2, · · · , N n } , ∃l ∈ {1, 2, · · · , L} , q n (I ni ) ⊂ K q l .
Moreover, the intervals I ni can be constructed in such a way that: Furthermore, recalling: where the λ ni are nonpositive real measures on [t 0 , min(τ, t 0 +δ)], and performing the same job as in the proof of proposition 18 (estimate 38), one obtains: There results: The measures λ ni are uniformly bounded with respect to n. By the equation of motion (55), we obtain that the real numbers Var (q + n ; [t 0 , min(τ, t 0 + δ)]) are uniformly bounded with respect to n, for n larger than N 0 . The assertion of step 2 is now a direct consequence of proposition 34.
Step 3. The function q lim constructed in step 2 satisfies the equation of motion: Moreover, the real measure R lim , q + lim +q − lim q lim is a nonpositive measure on the interval [t 0 , min(τ, t 0 + δ)[. Proof of step 3. We denote by M b ([a, b], R) the Banach space of all bounded real measures on an interval [a, b]. By estimate 56, we can find N bounded real measures λ ilim such that: where another subsequence has been extracted, if necessary. Writing the equation of motion (55) in local charts, we have: Furthermore, lim n→+∞ f q n ,q + n ; t dt = f q lim ,q + lim ; t dt in M b weak*, by Lebesgue's Dominated Convergence Theorem. Therefore, we obtain easily: the weak* topology being Hausdorff. Now, by formulae (12) we have to prove: Consider ]a, b[⊂ [t 0 , min(τ, t 0 + δ)] such that: ∀s ∈]a, b[, ϕ i (q lim (s)) < 0.
The interval ]a, b[ is the union of the compact intervals K j = [a + 1/j, b − 1/j] (j ∈ N * ). Fix j ∈ N * . For n large enough: ∀s ∈ K j , ϕ i (q n (s)) < 0, so λ in|K j = 0. We deduce: Therefore, λ ilim|]a,b[ = 0 and this achieves the proof of inclusion (57) and therefore of the first assertion of step 3.
For the second assertion of step 3, we are going to prove actually: (58) Fix such t 1 , t 2 and arbitrary ε > 0. We have: By the right-continuity of the function t → q + lim (t 2 ) q lim (t) and the results of step 2, we can find t 1 , t 2 ∈ [t 0 , min(τ, t 0 + δ)[ (t 1 < t 2 ) and an integer N 0 such that: Moreover, by Lebesgue's Dominated Convergence Theorem, N 0 may be assumed large enough to ensure: It is easily deduced that: Since ε is arbitrary and ]t 1 ,t 2 ] R n , (q + n +q − n ) q n is nonpositive (proposition 7), the conclusion (assertion (58)) follows.
Step 4. Consider an arbitrary instant t g ∈]t 0 , min(τ, t 0 + δ)[ such that: Then, q lim satisfies the impact constitutive equation at instant t g : Proof of step 4. Consider a local chart (U, ψ) centered at q lim (t g ) such that: the cardJ(q lim (t g )) first components of ψ(q) are α i ϕ i (q) (i ∈ J(q lim (t g ))), where the α i are some fixed positive real constants. -∀q ∈ U, J(q) ⊂ J(q lim (t g )), the matrix (g ij ( lim (t g ))) is the identity matrix: We have to prove: First, we are going to prove: and (61) Fix ε > 0 arbitrary, and pick a positive real number η small enough and an integer N 0 large enough to ensure: ∀t ∈ [t g − η, t g + η], ∀n ≥ N 0 , q lim (t) ∈ U and q n (t) ∈ U.
Let V be a positive real constant, large enough to majorize all the quantities q +i n (t) and Var q +i n ; [t g − η, t g + η] , when i, t and n wander respectively in the sets {1, 2, · · · , d}, [t g − η, t g + η] and {n ∈ N ; n ≥ N 0 }. We may assume that η is small enough and N 0 large enough to ensure: Multiplying the equation of motion (36) by (q +i n +q +i n )/2 and integrating over ]t 1 , t 2 ], we obtain easily: which gives, by lemma 17 and the desired conclusion (60) for sufficiently small η. For the second assertion (61), suppose we have in addition: There results that λ in vanishes on [t 1 , t 2 ] and integration of the equation of motion (37) gives: and, therefore the desired conclusion (61) for sufficiently small η. Now, let us come back to the proof of assertions (59). Fix i ∈ {1, 2, · · · , d}. Only the following four cases are possible: and φ > 0 We examine them successively.
Fix ε > 0 arbitrary. By assertion (61), we can pick a positive real number η small enough and an integer N 0 large enough to ensure: by the choice of the chart we made. Actually, η can be assumed small enough to ensure: q +i lim (t) −q +i lim (t g ) ≤ ε, by proposition 32. By step 2, we can find t 1 ∈ [t g − η, t g [ and t 2 ∈]t g , t g + η] such that: lim n→+∞q +i n (t k ) =q +i lim (t k ) (k = 1, 2) and, therefore, N 0 can be assumed large enough to ensure: 1, 2). Then, we have: Since ε is arbitrary, we get the desired conclusion: . Fix ε > 0 arbitrary. By assertion (60), we can pick a positive real number η small enough and an integer N 0 large enough to ensure: q +i n (t 2 ) ≤ q +i n (t 1 ) +ε. Exactly as in case 1, η is assumed sufficiently small to ensure: q +i lim (t) −q +i lim (t g ) ≤ ε, and N 0 large enough to have: , for some t 1 ∈ [t g − η, t g [ and some t 2 ∈]t g , t g + η]. We get: which gives the desired conclusioṅ q +i lim (t g ) = 0, since ε is arbitrary.
Performing the same job for instant t 0 instead of t 0 , we extend the conclusion to interval [t 0 , min(τ, t 0 + 3δ/2)]. Processing so inductively a large enough number of times, we obtain the desired conclusion.

Remark.
A straightforward modification of the proof of step 4 shows that the conclusions of theorem 20 hold if we only assume that φ is continuous and constant on each fiber: The conclusions of theorem 20 also hold if φ is only assumed continuous and if, moreover, we have: ∀t ∈ [t 0 , T [, CardJ(q(t)) ≤ 1.

Indications on the numerical computation of the solutions
Consider the problem P described in section 3.4. Assume furthermore, for sake of simplicity, that the impact function φ is constant. The maximal solution associated with the initial condition (q 0 , v 0 ) at time t 0 = 0, is denoted by (T m , q). We consider a local chart (U, ψ) at q 0 and a positive real number T such that: ∀t ∈ [0, T ], q(t) ∈ U.
Actually, there may happen that the function q n cannot be defined on [0, T ] if there exists an integer k n such that q n,k n +1 ∈ U . In such a case, the function q n is defined only on [0, t n,k n ].
This type of algorithm was introduced by Moreau and used without further justifications. It should be stressed that one cannot hope that the sequence of approximants (q n ) converges in general towards the solution q, since continuous dependence on initial condition does not hold in general. Actually, it is easy to build an explicit example, in the spirit of the example of section 7, where the sequence (q n ) does not converge pointwisely towards any function at all. However, in the special case where all the multiple impacts are orthogonal, things behave nicely and we have: Theorem 21 Suppose that the solution q is such that all multiple impacts are orthogonal: ∀t ∈ [0, T ], (dϕ i (q(t))) i∈J(q(t)) is orthogonal in T * q(t) Q, (with the convention that the empty set is orthogonal). Then, there exists an integer N 0 such that the function q n is well-defined on [0, T ] for n ≥ N 0 . Moreover, the sequence (q n ) converges uniformly on [0, T ] towards q (or more precisely towards ψ(q)).
Theorem 21 can be proved along the same steps as those of the proof of theorem 20. The necessary adaptation of the details is left to the reader.

Appendix: the class of motion M M A(I; Q)
In this section, we are going to define the concept of vector field with bounded variation along a locally absolutely continuous curve on a riemannian manifold. The definition and basic properties of absolutely continuous functions and functions with bounded variation from a real interval to a finite-dimensional normed vector space are supposed to be known. The reader is refered to Rudin [17] and Moreau [13]. These concepts are intimately connected with measure theory. Two expositions of measure theory compete: the set-theoretic approach (see for example Rudin [17]) and the duality approach (see for example Bourbaki [6]). Both approach are connected by Riesz's representation theorem. In this paper, we stick to the duality approach. If I is a real interval and E a real finite-dimensional normed vector space, C 0 c (I; E) will denote the space of continuous functions from I to E with compact support. A measure on I with values in E (resp. E * ) will be any linear form μ on C 0 c (I; E * ) (resp. C 0 c (I; E)) satisfying the following continuity property: When the constant M a,b can be found independent on a and b, the measure μ is said bounded. For all what concerns measure theory, the reader is refered to Bourbaki [6] where he will note the definition of the support Supp μ of a measure μ (Bourbaki [6], p. 64).
The following list of definitions and propositions aims at carrying these concepts over riemannian manifolds. The cornerstone is, of course, the identification of tangent spaces at different points of a curve by means of parallel translation.
This appendix is also an occasion to state precisely the classical theorems which are used in this paper.
Definition 22 Let I be a real interval and q : I → Q a curve on Q. q is said locally absolutely continuous if, for all t in I, there exists a compact neighbourhood J of t in I and a chart (U, ψ) such that: q(J) ⊂ U , ψ • q : J → R d is absolutely continuous.
Since Q can be covered by a countable collection of chart domains, there results from Lebesgue's theorem that q(t) admits a tangent vectorq(t) ∈ T q(t) Q for dt-almost all t in I where dt denotes the Lebesgue measure on the real line (and also its restriction on I). The riemannian structure on Q and the Cauchy-Lipschitz-Caratheodory theorem allow us to define classically a parallel translation operator along q, τ t,s : T q(s) Q → T q(t) Q (see, for example, Chavel [7], p. 7). τ t,s is defined for all (s, t) ∈ I 2 .
Definition 23 Let q be a locally absolutely continuous curve from I to Q. One calls vector (resp. 1-form) field on q(t) any mapping X (resp. X * ) from I to T Q (resp. T * Q) such that: ∀t ∈ I, Π Q (X(t)) = q(t) (resp. Π * Q (X * (t)) = q(t)).
A vector field X on q(t) (resp. a 1-form field X * on q(t)) will be said locally absolutely continuous (resp. absolutely continuous, resp. locally with bounded variation, resp. with bounded variation) if there exists t 0 in I such that the mapping: is locally absolutely continuous (resp. absolutely continuous, resp. locally with bounded variation, resp. with bounded variation). If X has bounded variation on I, its variation over I is by definition: Var(X(s); I) = Var(τ t 0 ,s (X(s)); I).
Proof. For measure taking values in a finite-dimensional vector space, the above statement is a classical direct consequence of Lebesgue-Radon-Nikodym theorem (see Rudin [17]). It is readily carried over manifolds by means of a locally finite partition of unity modelled on chart domains.
Definition 27 Let X be a vector field with locally bounded variation on an absolutely continuous curve (I, q) and t 0 an arbitrary element of I. We denote by d t 0 X the Stieljes measure (see Moreau [13]) associated to the mapping with locally bounded variation: For Y ∈ C 0 c (I, q; T Q) and Y * ∈ C 0 c (I, q; T * Q), the linear forms: turns out to be independent on a particular choice of t 0 and define measures on q taking respectively values in T * Q and T Q. Proof. This is an immediate consequence of definition 27 and of the similar statement for functions taking values in a finite-dimensional normed vector space.
Proposition 28 ensures the consistency of our notations. Let us now turn to practical calculations in charts.
Proposition 29 Let (U, ψ) be a chart on Q, (I, q) an absolutely continuous curve on Q such that q(I) ⊂ U and X a vector field on (I, q). The components (X i ) (i = 1, 2, · · · d) of X in the natural chart of T Q associated with ψ are real functions defined on I. The vector field X is locally absolutely continuous (resp. absolutely continuous, resp. locally with bounded variation, resp. with bounded variation) if and only if every function X i is locally absolutely continuous (resp. absolutely