Frank Vega

ON UP VERSUS NP AND INFINITE TURING MACHINES

Keywords: 2000 Mathematics Subject Classification. Primary 68Q15, Secondary 11A07 P, NP, UP, NP-complete, Quadratic Congruences

We define a problem that we call General Quadratic Congruences. We show General Quadratic Congruences is an NP-complete problem. Moreover, we prove General Quadratic Congruences is also in UP. In this way, we demonstrate that UP = NP.

Introduction

P versus N P is a major unsolved problem in computer science [START_REF] Fortnow | The Golden Ticket: P, NP, and the Search for the Impossible[END_REF]. This problem was introduced in 1971 by Stephen Cook [START_REF] Arora | Computational complexity: A modern approach[END_REF]. It is considered by many to be the most important open problem in the field [START_REF] Fortnow | The Golden Ticket: P, NP, and the Search for the Impossible[END_REF]. It is one of the seven Millennium Prize Problems selected by the Clay Mathematics Institute to carry a US$1,000,000 prize for the first correct solution [START_REF] Fortnow | The Golden Ticket: P, NP, and the Search for the Impossible[END_REF].

In 1936, Turing developed his theoretical computational model [START_REF] Arora | Computational complexity: A modern approach[END_REF]. The deterministic and nondeterministic Turing machines have become in two of the most important definitions related to this theoretical model for computation. A deterministic Turing machine has only one next action for each step defined in its program or transition function [START_REF] Papadimitriou | Computational Complexity[END_REF]. A nondeterministic Turing machine could contain more than one action defined for each step of its program, where this one is no longer a function, but a relation [START_REF] Papadimitriou | Computational Complexity[END_REF].

Another huge advance in the last century has been the definition of a complexity class. A language over an alphabet is any set of strings made up of symbols from that alphabet [START_REF] Thomas | Introduction to Algorithms[END_REF]. A complexity class is a set of problems, which are represented as a language, grouped by measures such as the running time, memory, etc [START_REF] Thomas | Introduction to Algorithms[END_REF].

In the computational complexity theory, the class P contains those languages that can be decided in polynomial time by a deterministic Turing machine [START_REF] Goldreich | Np-Completeness[END_REF]. The class N P consists in those languages that can be decided in polynomial time by a nondeterministic Turing machine [START_REF] Goldreich | Np-Completeness[END_REF].

The biggest open question in theoretical computer science concerns the relationship between these classes: Is P equal to N P ? In 2002, a poll of 100 researchers showed that 61 believed that the answer was not, 9 believed that the answer was yes, and 22 were unsure; 8 believed the question may be independent of the currently accepted axioms and so impossible to prove or disprove [START_REF] William | The P=?NP poll[END_REF].

Another major complexity class is U P . The class U P has all the languages that are decided in polynomial time by a nondeterministic Turing machines with at most one accepting computation for each input [START_REF] Valiant | Relative Complexity of Checking and Evaluating[END_REF]. It is obvious that P ⊆ U P ⊆ N P [START_REF] Papadimitriou | Computational Complexity[END_REF]. Whether P = U P is another fundamental question that it is as important as it is unresolved [START_REF] Papadimitriou | Computational Complexity[END_REF]. All efforts to solve the P versus U P problem have failed [START_REF] Papadimitriou | Computational Complexity[END_REF]. Nevertheless, we prove U P = N P .

Theoretical notions

Let Σ be a finite alphabet with at least two elements, and let Σ * be the set of finite strings over Σ [START_REF] Arora | Computational complexity: A modern approach[END_REF]. A Turing machine M has an associated input alphabet Σ [START_REF] Arora | Computational complexity: A modern approach[END_REF]. For each string w in Σ * there is a computation associated with M on input w [START_REF] Arora | Computational complexity: A modern approach[END_REF]. We say that M accepts w if this computation terminates in the accepting state, that is, M (w) = "yes" [START_REF] Arora | Computational complexity: A modern approach[END_REF]. Note that M fails to accept w either if this computation ends in the rejecting state, or if the computation fails to terminate [START_REF] Arora | Computational complexity: A modern approach[END_REF].

The language accepted by a Turing machine M , denoted L(M ), has an associated alphabet Σ and is defined by

L(M ) = {w ∈ Σ * : M (w) = "yes"}.
We denote by t M (w) the number of steps in the computation of M on input w [START_REF] Arora | Computational complexity: A modern approach[END_REF].

For n ∈ N we denote by T M (n) the worst case run time of M ; that is

T M (n) = max{t M (w) : w ∈ Σ n }
where Σ n is the set of all strings over Σ of length n [START_REF] Arora | Computational complexity: A modern approach[END_REF]. We say that M runs in polynomial time if there exists k such that for all n, T M (n

) ≤ n k + k [1]. Definition 1.1. A language L is in class P if L = L(M )
for some deterministic Turing machine M which runs in polynomial time [START_REF] Arora | Computational complexity: A modern approach[END_REF].

We state the complexity class N P using the following definition. Definition 1.2. A verifier for a language L is a deterministic Turing machine M , where L = {w : M (w, c) = "yes" for some string c}. We measure the time of a verifier only in terms of the length of w, so a polynomial time verifier runs in polynomial time in the length of w [START_REF] Sipser | Introduction to the Theory of Computation[END_REF]. A verifier uses additional information, represented by the symbol c, to verify that a string w is a member of L. This information is called certificate.

Observe that, for polynomial time verifiers, the certificate is polynomially bounded by the length of w, because that is all the verifier can access in its time bound [START_REF] Sipser | Introduction to the Theory of Computation[END_REF].

Definition 1.3. N P is the class of languages that have polynomial time verifiers [START_REF] Sipser | Introduction to the Theory of Computation[END_REF].

In addition, we can define another complexity class called U P . Definition 1.4. A language L is in U P if every instance of L with a given certificate can be verified by a polynomial time verifier, and this verifier machine only accepts at most one certificate for each problem instance [START_REF] Hemaspaandra | Unambiguous Computation: Boolean Hierarchies and Sparse Turing-Complete Sets[END_REF]. More formally, a language L belongs to U P if there exists a polynomial time verifier M and a constant c such that if x ∈ L, then there exists a unique certificate

y with |y| = O(|x| c ) such that M (x, y) = "yes", if x / ∈ L, there is no certificate y with |y| = O(|x| c ) such that M (x, y) = "yes" [8].
A function f : Σ * → Σ * is a polynomial time computable function if some deterministic Turing machine M , on every input w, halts in polynomial time with just f (w) on its tape [START_REF] Sipser | Introduction to the Theory of Computation[END_REF]. Let {0, 1} * be the infinite set of binary strings, we say that a language L 1 ⊆ {0, 1} * is polynomial time reducible to a language L 2 ⊆ {0, 1} * , written L 1 ≤ p L 2 , if there exists a polynomial time computable function f : {0, 1} * → {0, 1} * such that for all x ∈ {0, 1} * ,

x ∈ L 1 iff f (x) ∈ L 2
where iff means "if and only if". An important complexity class is NP-complete [START_REF] Goldreich | Np-Completeness[END_REF]. A language L ⊆ {0, 1} * is NP-complete if

(1) L ∈ N P , and

(2) L ≤ p L for every L ∈ N P . Furthermore, if L is a language such that L ≤ p L for some L ∈ NP-complete, then L is in NP-hard [START_REF] Thomas | Introduction to Algorithms[END_REF]. Moreover, if L ∈ N P , then L ∈ NP-complete [START_REF] Thomas | Introduction to Algorithms[END_REF]. If any single NP-complete problem can be solved in polynomial time, then every N P problem has a polynomial time algorithm [START_REF] Thomas | Introduction to Algorithms[END_REF]. No polynomial time algorithm has yet been discovered for any NP-complete problem [START_REF] Fortnow | The Golden Ticket: P, NP, and the Search for the Impossible[END_REF]. The distinct prime factors of a positive integer n >= 2 are defined as the ω(n) numbers p 1 , . . . , p ω(n) in the prime factorization

Results

n = p a1 1 × p a2 2 × . . . × p a ω(n)
ω(n) . Lemma 2.5. There will exist a constant α, such that there are infinite positive integers n which complies with ω(n) ≤ α × ln ln n.

Proof. The average order of ω(n) is ω(n) ∼ ln ln n [START_REF] Hardy | Ramanujan: Twelve Lectures on Subjects Suggested by His Life and Work[END_REF]. Consequently, it will exist the constant α. If the congruence x 2 ≡ i 2 ( mod b ) has a solution, that solution is necessarily a solution to each of the prime power congruences x 2 ≡ i 2 ( mod p ei i ) when p ei i divides b [START_REF] Nagell | Introduction to Number Theory[END_REF]. For any prime p r , a necessary condition for x 2 ≡ i 2 ( mod p er r ) to have a solution is for x 2 ≡ i 2 ( mod p r ) to have a solution (to see this, note that if x 2 -i 2 is divisible by p er r then it is certainly divisible by p r ). Now, suppose x 2 ≡ i 2 ( mod p er r ) where p er r is a prime power which divides b . Then x 2 -i 2 ≡ (x -i) × (x + i) ≡ 0( mod p er r ). Thus p er r divides the product (x -i) × (x + i) and so p r divides the product as well. If p r = 2 and p r divides (x -i) × (x + i), then this is because x ≡ i( mod p r ) since the sum and the subtraction of two integers is even when both are even or odd at the same time. If p r is an odd prime and divides both (x -i) and (x + i), then p r would divide both their sum and their difference, 2 × x and -2 × i. Since p r is an odd prime, p r does not divide 2 and so p r would divide both x and i which can be translated to x ≡ i( mod p r ). It follows that p r either divides (x -i) or (x + i) but not both. Since p r divides (x -i) × (x + i), it only divides one of (x -i) and (x + i). Therefore, either x ≡ i( mod p r ) or x ≡ -i( mod p r ).

In this way, we prove for every prime p r that divides b we will have either x ≡ i( mod p r ) or x ≡ -i( mod p r ). Conversely, if we find all the possible solutions to each of the prime congruences, then we can use the Chinese Remainder Theorem to produce a solution to the original problem, that is to find the value of i [START_REF] Thomas | Introduction to Algorithms[END_REF]. Since the Chinese Remainder Theorem can be solved in polynomial time O(ln β b), then the remaining order will depend on the computation of all possible solutions. Since we only have two possible choices for each prime factor, then the order will depend on O(2 ω(b ) ). Since ω(b ) ≤ ω(b) ≤ α × ln ln b, then the final order will be of O(ln ) such that M (a, b, c, d, x) = "yes" [START_REF] Hemaspaandra | Unambiguous Computation: Boolean Hierarchies and Sparse Turing-Complete Sets[END_REF].

β b × 2 α×ln ln b ) = O(ln β b × ln α b) = O(ln k b) for a constant k = β + α.
The constant e exists because SQC ∈ N P . We denote this problem as CQC i .

Theorem 2.10. For every positive integer i we have that CQC i ∈ U P .

Proof. For i = 1, then CQC 1 = SQC and thus CQC 1 ∈ U P . Suppose for some i = k, then CQC k ∈ U P . Let's prove CQC k+1 ∈ U P . We will take an arbitrary instance (a, b, c, d) and some prime number p > 2 which does not divide b. The prime p can be taken in polynomial time in relation to log 2 b. Certainly, this can be done choosing a candidate from 3 to log 2 2 b because ω(b) ≤ log 2 b and the nth prime number is approximately equal to n × ln n < n 2 [START_REF] Nagell | Introduction to Number Theory[END_REF]. Let's take the number q = p ln 2 b . Since the congruence property

d × x 2 ≡ a( mod b) complies with q × d × x 2 ≡ q × a( mod q × b)
then Q(a, b, c, d, x) = true if and only if Q(q × a, q × b, c, q × d, x) = true. However, if the instance (a, b, c, d) ∈ CQC k+1 , then the instance (q×a, q×b, c, q×

d) ∈ CQC k because ω(q × b) = ω(b) + 1 ≤ (k + 1) × α × ln ln b + 1. Since p > 2 then q = p ln 2 b > b ln b . Therefore k × α × ln ln(q × b) > k × α × ln ln b ln b and this complies for k > 1 with k × α × ln ln b ln b = k × α × ln(ln b × ln b) = k × α × ln ln 2 b = 2 × k × α × ln ln b > (k + 1) × α × ln ln b + 1 ≥ ω(q × b).
In this way, we can reduce in polynomial time CQC k+1 to CQC k , since the calculation of q will be polynomial in relation to ln b if we use the exponentiating by squaring [START_REF] Thomas | Introduction to Algorithms[END_REF]. Since U P is closed under reductions and CQC k ∈ U P , it follows that CQC k+1 ∈ U P . Hence, by mathematical induction we have proved CQC i ∈ U P for every positive integer i [9]. Definition 2.11. A Turing machine is a quadruple M = (K, Σ, δ, s). K is a finite set of states; s ∈ K is the initial state. Σ is a finite set of symbols (we say Σ is the alphabet of M ). We assume K and Σ are disjoint sets. Σ always contains the special symbols and : The blank and first symbol. Finally, δ is a transition function, which maps K ×Σ to (K ∪{h, "yes", "no"})×Σ×{←, →, -}. We assume that h (the halting state), "yes" (the accepting state), "no" (the rejecting state), and the cursors directions ← for "lef t", → for "right" and -for "stay", are not in K ∪ Σ.

We can define the infinite Turing machine formally using the above notions. Definition 2.12. An infinite Turing machine is a quadruple M = (K, Σ, δ, s), where Σ and s are exactly as an ordinary Turing machine, but the set K must contain infinite states. δ is a transition function that must reflect the complexities of infinite states. Intuitively, δ decides the next state as before, but also contains infinite actions (steps). Nevertheless, the final states are unique, that is there should be only one h (the halting state), "yes" (the accepting state) and "no" (the rejecting state) in M . Definition 2.13. We define the bit length of a Turing machine M as the string length of the encoding of M by the binary alphabet of the Universal Turing machine. A Turing machine M has an infinite bit length if M is an infinite Turing machine. We say that a Turing machine has a finite bit length when M is an ordinary Turing machine.

Theorem 2.14. If a language L is not in U P , then it would have a polynomial time verifier M j of infinite bit length such that this verifier machine only accepts at most one certificate for each problem instance of L.

Proof. In this case, the certificate is not longer useful and thus we can always pass a single constant string as certificate such that we will have M j (x, ρ) = "yes" if and only if x ∈ L and ρ = or M j (x, ρ) = "no" if and only if x / ∈ L or ρ = . Indeed, M j (x, ρ) can be decided in polynomial time because M j could decide every x ∈ L in polynomial time since it may store an infinite amount of program actions for the decision of elements in L and this will be a valid verification too. In addition, for every string ρ, M j can compute in polynomial time whether ρ = . Furthermore, since is a constant string, then there exists a constant e such that | | = O(|x| e ) for every x ∈ L. Therefore, this kind of Turing machine M j of infinite bit length only accepts at most one certificate for each problem instance of L. Lemma 2.15. A language L is in U P if and only if it would have a polynomial time verifier M j of finite bit length such that this verifier machine only accepts at most one certificate for each problem instance of L.

Proof. This is a consequence of the self definition of the class U P . Theorem 2.16. GQC ∈ U P .

Proof. For every positive integer i the set CQC i contains infinite elements. Indeed, for some positive integer b there are infinite numbers n such that ω(b) = ω(n), because there are infinite prime numbers [START_REF] Nagell | Introduction to Number Theory[END_REF]. Since for every positive integer i we have that CQC i ∈ U P , then every language CQC i will have a polynomial time verifier M i , and this verifier machine only accepts at most one certificate for each problem instance of CQC i [START_REF] Hemaspaandra | Unambiguous Computation: Boolean Hierarchies and Sparse Turing-Complete Sets[END_REF]. We denote l(M i ) as the bit length of M i . We also denote |CQC i | as the cardinality of CQC i . Certainly, lim i→∞ Theorem 2.17. U P = N P .

Proof. Since GQC will be complete for N P , thus all language in N P will reduce to U P . Since U P is closed under reductions, it follows that U P = N P .

Conclusions

There is a previous known result which states that P = U P if and only if there are no one-way functions [START_REF] Papadimitriou | Computational Complexity[END_REF]. Indeed, for many years it has been accepted the P versus U P question as the correct complexity context for the discussion of the cryptography and one-way functions [START_REF] Papadimitriou | Computational Complexity[END_REF]. For that reason, the proof of Theorem 2.17 negates this current idea and also the belief that U P = N P is a very unlikely event. In addition, this demonstration might be a shortcut to prove P = N P , because if anybody proves that P = U P , then he will be proving the outstanding and difficult P versus N P problem at the same time [START_REF] Fortnow | The Golden Ticket: P, NP, and the Search for the Impossible[END_REF]. Furthermore, if we obtain a possible proof of P = N P , then this work would also contribute to show P = U P .

Definition 2 . 1 .

 21 Given five positive integers a, b ,c, d and x, the boolean function Q(a, b, c, d, x) is true if and only if x < c and d × x 2 ≡ a( mod b) [9]. Definition 2.2. QUADRATIC CONGRUENCES INSTANCE: Positive integers a, b and c, such that we have the prime factorization of b. QUESTION: Is there a positive integer x such that Q(a, b, c, 1, x) = true? We denote this problem as QC. QC ∈ NP-complete [4]. Let's define another problem. Definition 2.3. GENERAL QUADRATIC CONGRUENCES INSTANCE: Positive integers a, b, c and d, such that we have the prime factorization of b. QUESTION: Is there a positive integer x such that Q(a, b, c, d, x) = true? We denote this problem as GQC. Theorem 2.4. GQC ∈ NP-complete. Proof. Since we can check Q(a, b, c, d, x) = true in polynomial time, then GQC ∈ N P . Indeed, the certificate x will be polynomially bounded by any instance (a, b, c, d) when Q(a, b, c, d, x) = true because x < c. In addition, we can reduce every instance (a, b, c) of QC into an instance (a, b, c, 1) of GQC in polynomial time where (a, b, c) ∈ QC iff (a, b, c, 1) ∈ GQC. Since QC ∈ NP-complete then GQC ∈ NP-complete.

Theorem 2 . 6 .

 26 Given four positive integers a, b, c and d, such that we have the prime factorization of b and ω(b) ≤ α × ln ln b, then we can check whether a positive integer x is the minimum that complies Q(a, b, c, d, x) = true in order O(ln k b) for a constant k. Proof. Suppose we have a positive integer i such that 0 < i < x and Q(a, b, c, d, i) = true. Hence, we will obtain d × x 2 ≡ d × i 2 ( mod b). Moreover, by a property of congruences we have x 2 ≡ i 2 ( mod b ) where b = b (d,b) and (d, b) is the greatest common divisor of d and b [9]. We can find (d, b) in polynomial time in relation to ln b just multiplying into a single number each maximum prime power p ei i that divides b when also p ei i divides d. This is possible because we have the prime factorization of b. We are going to assume b = 1, because in case of (d, b) = b then x should be necessarily equal to 1.

Definition 2 . 7 .Theorem 2 . 8 .

 2728 SIMPLE QUADRATIC CONGRUENCES INSTANCE: Positive integers a, b, c and d, such that we have the prime factorization of b and ω(b) ≤ α × ln ln b. QUESTION: Is there a positive integer x such that Q(a, b, c, d, x) = true? We denote this problem as SQC. SQC ∈ U P . Proof. We show a polynomial time verifier, and this verifier machine only accepts at most one certificate for each problem instance of SQC [8]. Given five positive integers a, b, c, d and x, we define the verifier machine M for SQC as follows: M (a, b, c, d, x) = "yes"iff x is the minimum such that Q(a, b, c, d, x) = true. SQC belongs to U P because the verifier M can run in polynomial time as we proved in Theorem 2.6 and there will be a constant e such that if (a, b, c, d) ∈ SQC, then there is a unique certificate x with |x| = O(|(a, b, c, d)| e ) such that M (a, b, c, d, x) = "yes", if (a, b, c, d) / ∈ SQC, there is no certificate x with |x| = O(|(a, b, c, d)| e

Definition 2 . 9 .

 29 COMPLEX QUADRATIC CONGRUENCES ON I INSTANCE: Positive integers a, b, c and d, such that we have the prime factorization of b and ω(b) ≤ i × α × ln ln b for a positive integer i. QUESTION: Is there a positive integer x such that Q(a, b, c, d, x) = true?

  every CQC i has infinite elements and the bit length of M i is finite becauseCQC i ∈ U P . Thus lim i→∞ l(Mi) |CQCi| = limi→∞ l(Mi)limi→∞ |CQCi| = 0. Moreover, we can assure that lim i→∞ |CQC i | = |GQC|, because lim i→∞ CQC i = GQC. Therefore, we obtain limi→∞ l(Mi) |GQC| = 0. By the definition of M i , we get lim i→∞ M i = M GQC where M GQC is a polynomial time verifier such that this verifier machine only accepts at most one certificate for each problem instance of GQC. However, M GQC might have an infinite bit length. Indeed, M GQC has a finite bit length if and only if GQC ∈ U P . Consequently, limi→∞ l(Mi) |GQC| = l(M GQC ) |GQC| = 0. Since GQC contains infinite elements and l(M GQC ) |GQC| = 0, then M GQC has a finite bit length and thus GQC ∈ U P .


A time-stepping scheme for inelastic collisions Numerical handling of the nonoverlapping constraint

Introduction

The present work was initially motivated by the necessity to handle particle collisions in the direct simulation of fluid-rigid particle mixtures. By direct simulation we mean that each particle is treated individually, and the Navier-Stokes (or Stokes) equations are solved in the moving domain occupied by the fluid. Many methods have been proposed recently to compute such flows. Some are based on an embedding of the solid phase in a global domain which is covered by a cartesian mesh (see Glowinski [START_REF] Glowinski | Finite element methods for incompressible viscous flow[END_REF]Chapter 8] or Bertoluzza et al. [START_REF] Bertoluzza | The fbm method: Semi-discrete scheme and some numerical experiments[END_REF]); the other class of methods consists in using a conforming mesh of the fluid domain (see Hu [START_REF] Hu | Direct simulation of flows of solid-liquid mixtures[END_REF], Johnson et al. [START_REF] Johnson | Simulation of multiple spheres falling in a liquid-filled tube[END_REF], Maury [START_REF] Maury | Direct simulations of 2D fluid-particle flows in biperiodic domains[END_REF]). All these approaches face the problem of body overlapping, which makes it difficult to guarantee the robustness of the computation. Different strategies have been used. In Hu [START_REF] Hu | Direct simulation of flows of solid-liquid mixtures[END_REF], the mesh is refined in the neighborhood of the interparticle gap, so that lubrication forces can be approximated with accuracy and consequently can be expected to prevent overlapping. Indeed, as the interparticle force between two smooth bodies separated by a viscous fluid behaves like -ε/ε (see Kim [START_REF] Kim | Microhydrodynamics: principles and selected applications[END_REF]), where ε is the distance, there cannot be any contact in finite time for reasonably regular external forces. Another commonly used strategy consists in adding short range repulsive forces between particles, which tend to prevent particles to overlap (see Glowinski et al. [START_REF] Glowinski | Direct simulation of the motion of neutrally buoyant circular cylinders in plane Poiseuille flow[END_REF]). Those methods have proved to behave quite satisfactorily in many situations, but they necessitate a fine tuning of some numerical parameters, and the actual minimal distance between bodies cannot be controlled a priori. In [START_REF] Maury | Direct simulations of 2D fluid-particle flows in biperiodic domains[END_REF], we introduced a heuristic method to control this minimal distance, by running at each time step a minimization procedure on a global functional of the particle positions, such that the minimal distance corresponding to the resulting configuration is greater than a prescribed "safety distance" ε > 0. This method, although robust in practice, even in the case of multiple contacts, is not consistent from the energy point of view, and its long-term effects on the computed flow are difficult to estimate. The approach we proposed in [START_REF] Maury | A many-body lubrication model[END_REF] is based on a first-order approximation of the lubrication forces exerted by the fluid in the interparticle gap. It is more respectful of the underlying physics, but the singular character of those forces (proportional to the reciprocal of the distance) introduces new constraints on the time step, and we must add that this approach, despite its very respect of local phenomena, does not seem to improve the overall modeling of fluid-particle mixtures significantly. It made us look for more basic tools, yet respectful of some basic physical properties like momentum balance and decreasing character of energy.

In the present work, we propose a new algorithm to handle the nonoverlapping constraint, which fits into the general class of Contact Dynamics methods, introduced by Moreau [START_REF] Moreau | Some numerical methods in multibody dynamics: Application to granular materials[END_REF]. The action of the surrounding fluid shall be simply described by a given external force field, so that from now on we consider "dry" systems of rigid bodies. Given the nature of the interaction between particles moving in a viscous fluid, which is dissipative, we chose to consider here purely inelastic collisions. The evolution problem on which our approach is based fits into the general framework of differential measure inclusions (see Schatzman [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF], or Moreau [START_REF] Moreau | Some numerical methods in multibody dynamics: Application to granular materials[END_REF]). From the theoretical point of view, the critical point is uniqueness, which necessitates strong regularity assumptions on the data: analyticity is required. Uniqueness is lost as soon as the external force is no longer analytic (see in Schatzman [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF] or Ballard [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF] counter-examples to uniqueness with a C ∞ force). This question of uniqueness in regard with the numerical algorithm will be addressed at the end of the paper. We refer to Stewart [START_REF] Stewart | Rigid-body dynamics with friction and impact[END_REF] for a general presentation of mathematical issues raised by those models, and a wide presentation of numerical methods. The most general result concerning convergence analysis of a numerical scheme for granular flows can be found in [START_REF] Stewart | Convergence of a time-stepping scheme for rigid-body dynamics and resolution of Painlevé's problem[END_REF], where Coulomb friction is taken into account.

The situation we consider here is less general than the aforementioned one, but the scheme we propose presents the following features:

1. Collisions are not treated as events 1 . All contacts which are likely to occur during a time step are indeed handled globally, without any prior prediction of actual violations of the non-overlapping constraint, and reactions are computed as a field of Kuhn-Tucker multipliers associated to first-order approximations of those constraints, expressed in terms of velocities. 2. As a consequence, this method can be proved to be unconditionally stable. For any choice of the time step, the kinetic energy of the approximate solution dissipates (when external forces are reduced to 0). Since there is no parameter to tune up (except for the stopping criterium associated to the Uzawa algorithm), this algorithm can be applied to very different situations with no modification. 3. Provided the quadratic minimization step is solved exactly, the scheme produces feasible configurations only (in the case of spherical bodies). 4. Non-smooth forcing terms can be taken into account, as soon as their mean values on subintervals can be computed. This makes it possible to handle multiple solutions numerically (see section 6.1).

The paper is organized as follows. In section 2 we present the rigid-body problem, and we introduce the natural functional spaces for positions, velocities and reaction forces. We then present the time-stepping scheme (section 3). In the next section we prove some general properties of the scheme, all of which are valid for any time step. Section 5 is devoted to the asymptotic behaviour of the scheme: as the time step h goes to 0, a subsequence of the computed solutions is shown to converge in some sense to a solution of the original problem, in the case of a single contact. Some numerical tests are finally presented, to illustrate the behaviour of the algorithm in different contexts: non-uniqueness, sticky particle model, and many-body motion (section 6).

Model problem

Evolution model

We consider the mechanical system of N rigid spheres (or discs, or line segments) in R d (with d = 3, 2, or 1) of radii (r i ) 1≤i≤N and masses (m i ) 1≤i≤N . We introduce the configuration space,

Q = {q = (q 1 , q 2 , . . . , q N ) ∈ R dN }
and the associated feasible set

Q 0 = {q ∈ Q, D ij (q) ≥ 0 ∀i, j , 1 ≤ i < j ≤ N } where D ij (q) = |q i -q j | -(r i + r j )
is the signed distance between spheres i and j (see Figure 2.1). Note that D ij is negative as soon as spheres i and j overlap.

As we do not consider rotations here, our configuration space is endowed with a natural Euclidean structure, so that the tangent space at any q ∈ Q can be identified with Q itself. We shall nevertheless denote by T Q = R dN this tangent space in order to distinguish velocities from positions. We define G ij = ∇D ij as the gradient of the distance between the two spheres i and j :

r i r j D ij q i q j e ij -e ij
G ij = (• • • , 0, -e ij , 0, . . . , 0, e ij , 0, • • • ) , e ij = q j -q i q j -q i , (1)
and C q as the set of feasible directions at q ∈ Q 0 ,

C q = v ∈ T Q , G ij • v ≥ 0 as soon as D ij (q) = 0 ,
and the outward normal cone to Q 0 at q as its polar cone:

N q = C • q = {h ∈ T Q , h • v ≤ 0 ∀v ∈ C q } = - i<j μ ij G ij (q), μ ij = 0 if D ij (q) > 0, μ ij ∈ R + if D ij (q) = 0 .
Let u = ( q1 , q2 , . . . , qN ) ∈ T Q denote the generalized velocity vector. We finally denote by f the force field acting on the spheres, and by M the mass matrix. The problem we are interested in may be phrased formally:

M q + N q f (2) u + = P C q u -, (3)
where u -(resp. u + ) is the left (resp. right) limit of the velocity vector at time t, and P C q is the Euclidean projection onto the closed convex cone C q . This model will be given a proper mathematical framework in the next section. As regards its mechanical sense, equation [START_REF] Bertoluzza | The fbm method: Semi-discrete scheme and some numerical experiments[END_REF], which can be read formally

∃h ∈ -N q such that M q = f + h, (4)
expresses the fact that overlapping is prevented by repulsive forces acting on each sphere along the normal vector at the contact point. When there is no contact, N q reduces down to {0}, so that (2) reads as the ordinary differential equation M q = f. Equation (3) provides the collision model. This is a particular case of a more general model (which includes possibly an elastic behaviour)

u + = u --(1 + e)P N q u -, ( 5
)
where e is a restitution coefficient, set to 0 in the case of a non-elastic impact law. Indeed, as N q and C q are mutually polar, it holds I -P N q = P C q , where I is the identity. (see Moreau [START_REF] Moreau | Décomposition orthogonale d'un espace Hilbertien selon deux cônes mutuellement polaires[END_REF]). Note that, in the non-elastic model we consider, some energy is lost during each collision. This kinetic energy is actually transformed into heat within the bodies. The reader may refer to Frémond [START_REF] Frémond | Non-smooth thermomechanics[END_REF], where this increase in temperature is explicitly integrated into a general thermomechanical model for collisions between rigid bodies.

Remark 1 For 1D problems (like the one which is presented in Section 6.2), it is in fact more natural to define Q 0 as a connected component of the set of all configurations q with no overlapping. Indeed, the evolution problem does not allow q to pass from a connected component to another (the bodies cannot leap across each other). With this convention, the feasible set

Q 0 = {q = (q 1 , . . . , q N ) ∈ R N , q i+1 -q i ≥ r i+1 + r i }
is closed and convex, and in this particular situation, N q identifies with the subdifferential of the indicatrix of Q 0 :

N q = ∂I Q 0 (q) with I Q 0 (q) = 0 if q ∈ Q 0 +∞ if q / ∈ Q 0
and hence q -→ N q is a maximal monotone operator. In general Q 0 is not convex, but it is not far from being convex: it can be proved easily to be uniformly η-prox-regular in the sense defined in Colombo [START_REF] Colombo | Sweeping by a continuous prox-regular set[END_REF], as soon as 1/η < 4 min(r i ).

Functional framework

We consider the time interval I = (0, T ) and we introduce the following functional spaces: W 1,1 = set of all R dN -valued functions which are absolutely continuous over the interval I ; BV = set of all dN vector-valued functions having bounded variation over I : BV is the set of all functions t → u(t) ∈ T Q , such that each component u of u verifies sup

S∈ N S n=1 |u(t n ) -u(t n-1 )| < ∞,
where S = (t 0 , t 1 , . . . , t N S ) runs over the set of increasing subdivisions of the time interval I ; M 1 = set of N(N -1)/2 vector2 -valued bounded measures on I : it contains all μ = (μ ij ) 1≤i<j ≤N such that μ ij is a continuous linear functional over the set C 0 (I ) of continuous functions over I , vanishing at 0 and T . The set of component-wise positive measures will be denoted by

M 1 + = {μ = (μ ij ) 1≤i<j ≤N ∈ M 1 , μ ij , ϕ ≥ 0 ∀ϕ ∈ C 0 (I ), ϕ ≥ 0}.
As regards the forcing term, we shall consider functions of t and q, with Carathéodory type regularity conditions (see Coddington et al. [START_REF] Coddington | Theory of ordinary differential equations[END_REF]) : let the function f be defined in Q 0 × I , measurable in t for fixed q and continuous in q for fixed t, such that

∃F ∈ L 1 (I ) s.t. |f(q, t)| ≤ F (t) ∀(q, t) ∈ Q 0 × I, ( 6
)
and f is uniformly Lipschitz with respect to the space variable in the feasible set:

∃k , f(q , t) -f(q, t) ≤ k q -q a.e. in I ∀q, q ∈ Q 0 . ( 7
)
In order to avoid a special treatment of the initial time, we shall pick an initial position q 0 in the interior of Q 0 , so that N q = {0}. It follows that the velocity is continuous at 0 + , and u(0) is defined without ambiguity.

We may now state the problem: Given a time interval I = (0, T ), a force f (which meets regularity assumptions ( 6)-( 7)), and initial conditions (q 0 , u

0 ) ∈ Q 0 × T Q , with D ij (q 0 ) > 0 for all i = j , Find (q, u, μ) ∈ W 1,1 × BV × M 1 + such that u(0) = u 0 , ( 8
)
q(t) = q 0 + t 0 u(s) ds ∈ Q 0 ∀ t ∈ I (9) M u = f(q, t) + i<j μ ij G ij (q(t)), ( 10
)
∀ i < j supp(μ ij ) ⊂ {t, D ij (q(t)) = 0}, (11)
u + (t) = u -(t) -P N q u -(t) ∀t ∈ I, ( 12
)
where P N q is the euclidean projection onto the closed convex cone N q . Equation [START_REF] Glowinski | Direct simulation of the motion of neutrally buoyant circular cylinders in plane Poiseuille flow[END_REF] is to be understood in the sense of distributions (it identifies two 0-th order distributions).

Remark 2 Note that f(•, t) is required to be in W 1,∞ (Q 0 ) and not necessarily in W 1,∞ (Q). This makes it possible to apply this approach to model interaction forces like

f = -K∇V (|q 2 -q 1 |),
where V : (0, +∞) → R + is smooth, decreasing and convex, possibly singular at 0 : since the distance between any q ∈ Q 0 and 0 is kept apart from 0 (it is always greater than twice the smallest radius), such a function is indeed lipschitzian. For instance an electrostatic potential, V (d) = K/d, meets the requirements.

Remark 3 As q has the W 1,1 regularity in time, we can consider that q is continuous. Similarly, any u ∈ BV (or u h in the next section), shall be identified with the right-continuous element of the class. With this convention, all pointwise identities or inequalities will be meant in the classical sense (i.e. everywhere).

Remark 4 Obstacles like walls or fixed bodies may be taken into account. For the sake of simplicity, we shall not introduce here extra notations for that purpose, but the set of Kuhn-Tucker multipliers μ could be extended to include reaction forces against obstacles.

3 Numerical scheme

Approximation spaces

Let h = T /N h be the time step, and let P 0 h denote the set of all those functions which are piecewise constant over I according to the uniform subdivision (0, h, 2h, . .

. , (N h -1)h, T ). Similarly, P 1

h denotes the set of continuous, piecewise affine functions according to the same subdivision. We introduce the following approximation spaces for trajectories, velocities, and reaction forces, respectively:

X h = {q h = (q i ) : I → Q , q i ∈ (P 1 h ) d , 1 ≤ i ≤ N}, V h = {u h = (u i ) : I → T Q , u i ∈ (P 0 h ) d , 1 ≤ i ≤ N}, R h = {μ h = (μ ij ) : I → R r , μ ij ∈ P 0 h , 1 ≤ i < j ≤ N }, where r = N(N -1)/2 is the number of constraints. For any u h ∈ V h (resp. μ h ∈ R h ) we shall denote by u n h (resp. μ n h ) its constant value in the subinterval [(n -1)h, nh).
Similarly, q n h denotes q h (nh), for any q h ∈ X h . As for external forces, we shall use the following mean value approximation

f h in V h of f(q, •) defined for any q ∈ Q by its N h values t ∈ [(n -1)h, nh) -→ f n h (q) = 1 h nh (n-1)h
f(q, s) ds.

Time-stepping scheme

The sequence of approximated fields (q h , u h , μ h ) ∈ X h ×V h ×R h is built according to the following scheme, to which we shall refer as scheme (S):

1. Initialization (q 0 h , u 0 h ) = (q 0 , u 0 ). ( 13
)
2. Compute u n+1 h as the solution to the constrained minimization problem min

u∈K h (q n h ) 1 2 u -u n h -hM -1 f n+1 h (q n h ) 2 M ( 14
)
with

|v| 2 M = Mv • v, and K h (q n h ) = {u ∈ T Q , D ij (q n h ) + hG ij (q n h ) • u ≥ 0}. ( 15
)
The approximate reaction field

μ n+1 h = (μ n+1 ij )
is the dual component of a solution to the associated saddle-point problem

L(u n+1 h , λ) ≤ L(u n+1 h , μ n+1 h ) ≤ L(v, μ n+1 h ), ( 16
)
for all λ ∈ R N(N-1)/2 + , v ∈ T Q , where L is the Lagrangian L(v, λ) = 1 2 v -u n h -hM -1 f n+1 h (q n h ) 2 M - 1≤i<j ≤N λ ij (D ij q n h ) + hG ij (q n h ) • v . ( 17
) Note that u n+1 h and μ n+1 h = (μ n+1 ij ) 1≤i<j ≤N are related by Mu n+1 h = Mu n h + hf n+1 h (q n h ) + h 1≤i<j ≤N μ n+1 ij G ij (q n h ). ( 18
)
3. Update the positions

q n+1 h = q n h + hu n+1 h . ( 19
)
Remark 5 The velocity field u n+1 h is uniquely defined (see [START_REF] Kim | Microhydrodynamics: principles and selected applications[END_REF]: it minimizes a strictly convex fonctional over a closed convex set), whereas μ n+1 h is not. So in a strict sense scheme (S) is an algorithm as far as q and u are concerned only.

The following remark explains why the previous scheme can be expected to approximate the original problem.

Remark 6

Step 2, which is the essential (and time-consuming) part of the scheme, can be written

Mu n+1 h -Mu n h h + ∂I K h (q n h ) u n+1 h f n+1 h (q n h ) ( 20
)
where ∂I K h (q n h ) is the subdifferential of the indicatrix function of K h (q n h ). Now observe that, for any q ∈ Q 0 , N q is the sum of all half lines -R + G ij for indices verifying D ij (q) = 0, and ∂I K h (q) (u) can be written as the same sum for indices i and j such that

D ij (q) + hG ij • u = 0.
As the latter quantity is a Taylor expansion of D ij (q + hu), the set ∂I K h (q n h ) (u n+1 h ) can be seen as a prediction of N q n+1 h . Therefore the scheme we propose can be considered, at least formally, as a semi-implicit time discretization of inclusion (2),

M du dt + N q f.
Note that the collision law (3) does not appear explicitly in the scheme. It is actually implicitly contained in [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF], which selects the velocity corresponding to the collision law we considered here. For this reason, this scheme is intrinsically dedicated to inelastic impact models, and we must say that there is no obvious way to adapt it to other impact laws.

Saddle-point problem

In the present approach we solve the saddle-point problem ( 16) by a Uzawa algorithm. Note that this choice is not mandatory, since the way the minimization procedure is performed is completely independant from the time-stepping scheme itself. Other algorithms could be implemented to perform this task, and the theoretical results presented in the next sections would remain unchanged. We introduce, for this section only, the vectors and matrices which are involved in the Uzawa procedure:

u = u n+1 h ∈ R dN , F = Mu n h + hf n+1 h (q n h ) ∈ R dN , μ = (μ n+1 ij ) 1≤i<j ≤N ∈ R r , D = (D ij (q n h )) 1≤i<j ≤N , C ∈ M r,dN (R) such that -C T μ = h 1≤i<j ≤N μ n+1 ij G ij (q n h ),
where r = N(N -1)/2 is the number of constraints. The problem can be put in the classical saddle-point form

Minimize J (v) = 1 2 Mv • v -F • v, over K = {v , Cv ≤ D}.
The Lagrangian of the problem is

L(v, λ) = 1 2 Mv • v -F • v + λ • (Cv -D),
and any saddle-point of L, i.e. any couple (u, μ) verifying

sup λ∈R r + L(u, λ) = L(u, μ) = inf v∈R dN L(v, μ),
is such that u minimizes J over K.

The Uzawa algorithm (see Ciarlet [START_REF] Ciarlet | Introduction à l'analyse numérique matricielle et à l'optimisation[END_REF]) consists in approximating a saddle-point (u, μ) of L by sequences (u k ) and (μ k ). A step decomposes into two substeps:

1. Solve the primal problem

Mu k+1 = F -C T μ k .
2. Update the Kuhn-Tucker multipliers

μ k+1 = + μ k + ρ(Cu k+1 -D) ,
where + is the orthogonal projection onto R r + :

μ ∈ R r -→ + (μ) = max(μ ij , 0) 1≤i<j ≤N .
Proposition 1 We suppose that

0 < ρ < ρ max = 2α C T C 2 ,
where α is the smallest eigenvalue of the matrix M. Then the sequence (u k ) converges to the solution u to the constrained minimization problem. The sequence

(μ k ) converges toward some μ such that (u, μ) is a saddle-point of L.
Proof The convergence of the primal sequence (u k ) is a classical result (see Ciarlet [START_REF] Ciarlet | Introduction à l'analyse numérique matricielle et à l'optimisation[END_REF]). As for the Kuhn-Tucker multipliers, the sequence (μ k ) can be shown to meet the assumptions of Opial's Lemma (see Haraux [START_REF] Haraux | Nonlinear evolution equations -global behaviour of solutions[END_REF]), and thus it converges (weakly) to some μ in the dual solution set. As the problem is finite-dimensional, the convergence is strong.

We conclude this section by two remarks on the actual implementation of this algorithm.

Remark 7 It is not necessary to assemble the whole matrix C. At a given time step, positions and velocities are known, so that most constraints (which correspond to particles far away from each other) can be eliminated a priori, in the spirit of what is presented in Sigurgeirsson [START_REF] Sigurgeirssons | Algorithms for particle-field simulations with collisions[END_REF]. For monodisperse situations (all radii are identical), the total number of potentially active constraints (i.e. number of rows of C) is κN/2, with κ = 6 in two dimensions (each disc is close to 6 other discs, at most), and κ = 12 in three dimensions.

Remark 8 A sufficient condition on ρ for the algorithm to converge can be obtained explicitly in standard situations. We give here a rough lower bound for ρ max in the monodisperse case. In this situation, the maximal number of likely contacts per body is smaller than the aforementioned parameter κ. The ellipticity constant α is simply min(m i ) 1≤i≤N . Using expression [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF] of G ij , we can write C T C ∈ M dN (R) by blocks:

C T C = h 2 (C ij ) 1≤i,j ≤N ,
where C ij is a d × d matrix defined by

C ii = j ∈J i e ij ⊗ e ij , C ij = -e ij ⊗ e ij for i = j,
with J i denoting the set of indices of the spheres which may get into contact with sphere i. By Gerschgorin circle theorem, the spectral radius of C T C (i.e. ||C T C|| 2 ) is smaller than the maximum among the 1-norms of its rows. Since e ij = 1, it follows that ||C T C|| 2 ≤ 2κh 2 √ d, where κ is the upper bound of (J i ) (cardinal number of J i ). Finally,

ρ max ≥ min(m i ) 2κh 2 √ d ,
where κ is 12 for 3D problems, as mentioned in the previous remark. It is to be noted that this lower bound is independent of N , the number of bodies.

Properties of the scheme

This section is devoted to some general properties of the algorithm in the general case (multiple contacts are allowed). Some of those properties will be used in the next section to prove convergence in the case of a single contact.

Proposition 2 (Feasibility) If the minimization step is solved exactly, then the scheme (S) produces feasible configurations only:

q h (t) ∈ Q 0 ∀h > 0, ∀t.
Proof As Step 2 of Scheme (S) is supposed to be solved exactly, then for any h, n, i = j ,

u n+1 h ∈ K h (q n h ) ⇒ D ij (q n h ) + hG ij (q n h ) • u n+1 h ≥ 0, therefore D ij (q n+1 h ) = D ij (q n h + hu n+1 h ) ≥ D ij (q n h ) + hG ij (q n h ) • u n+1 h (since D ij is convex) ≥ 0. ( 21
)
We established that q n h ∈ Q 0 for any h, n. Now for any t = nh+θh, with θ ∈ [0, 1], one has

D ij (q h (t)) = D ij (q n h + θhu n+1 h ) ≥ D ij (q n h ) + θhG ij (q n h ) • u n+1 h (since D ij is convex) = (1 -θ)D ij (q n h ) + θ D ij (q n h ) + hG ij (q n h ) • u n+1 h (22) ≥ 0. ( 23
)
Note that this property is no longer true if the bodies are not spherical, because the functions q -→ D ij (q) are not convex in general.

Proposition 3 (Stability) Let (q h , u h ) h>0 be a sequence built according to Scheme (S), for a given f. Assuming again that the minimization procedure (Step 2) is performed exactly, then u h verifies the following estimate:

|u h (t)| M ≤ |u 0 | M + 1 √ α t+h 0 F (s) ds, ( 24
)
where |•| M is the euclidean norm associated to the symmetric positive definite matrix M:

|u| M = √ Mu • u,
F is the function which dominates f (see condition ( 6)), and α is the smallest eigenvalue of M (α = min(m i )).

Proof Let us first establish that, for every n, i = j ,

μ n+1 ij G ij (q n h ) • u n+1 h ≤ 0. ( 25
)
It is a direct consequence of the fact that a Kuhn-Tucker multiplier may be different from zero only if the corresponding constraint is active. Indeed, either μ n+1 ij = 0 (the constraint is non-active), or, if μ n+1 ij > 0, then necessarily (the constraint is active)

D ij (q n h ) + hG ij (q n h ) • u n+1 h = 0, so that G ij (q n h ) • u n+1 h = -D ij (q n h )/ h ≤ 0 by Proposition 2.
We perform the scalar product of the relation [START_REF] Moreau | Décomposition orthogonale d'un espace Hilbertien selon deux cônes mutuellement polaires[END_REF] with the velocity field u n+1 h :

u n+1 h 2 M = Mu n h • u n+1 h + hf n+1 h (q n h ) • u n+1 h +h 1≤i<j ≤N μ ij G ij • u n+1 h ≤0 ≤ u n h M u n+1 h M + h f n+1 h (q n h ) • u n+1 h with h f n+1 h (q n h ) • u n+1 h = (n+1)h nh f(q n h , s) • u n+1 h ds = (n+1)h nh M -1 f(q n h , s) • Mu n+1 h ds ≤ 1 √ α u n+1 h M (n+1)h nh F (s) ds.
We finally get, by summing up over steps 0, 1, . . . , n,

u n+1 h M ≤ |u 0 | M + 1 √ α (n+1)h 0 F (s) ds,
which gives inequality (24).

Remark 9 In the case f ≡ 0, the previous proposition expresses the fact that kinetic energy decreases at the discrete level, for any time step.

Proposition 4 (Momentum balance) If the forcing term does not depend on q, then the total momentum balance is exactly verified at the right end of each subinterval by the approximated solution.

Proof Let us denote by P(u) the total momentum associated to the velocity field u = (u 1 , . . . , u N ):

P(u) = N i=1 m i u i ,
and by f i the force acting on sphere i. As the approximated reaction forces verify the law of action and reaction, it follows from [START_REF] Moreau | Décomposition orthogonale d'un espace Hilbertien selon deux cônes mutuellement polaires[END_REF] that

P(u n+1 h ) = P(u n h ) + (n+1)h nh f(s) ds,
where f = f i is the resultant force, so that

lim t→(n+1)h -P(u h (t)) = P(u n+1 h ) = P(u 0 ) + (n+1)h 0 f(s) ds,
which is exactly the momentum balance.

Remark 10 Note that the previous property remains valid even if the saddle-point problem is not solved exactly.

Remark 11 As a direct consequence of the previous property, the computed trajectory q h of the center of mass of the system turns out to be an interpolation of its exact trajectory (which is uniquely defined, even if particle trajectories are not), which can be expressed with obvious notations

q = q 0 + u 0 t + t 0 τ 0 M -1 f(s) ds dτ.
We shall complete this section by a proposition concerning the Kuhn-Tucker multipliers. As u n+1 h is uniquely defined, so is the global reaction term

G q n h (μ n+1 h ) = 1≤i<j ≤N μ n+1 ij G ij (q n h ).
As for the vector μ n+1 h itself, uniqueness does not hold in general, as shown by the following example: considering a two-dimensional collection of N identical discs in a cristal-like configuration such that most discs are in contact with 6 others, the number of active constraints is asymptotically 3N, whereas the number of degrees of freedom is 2N . This non-uniqueness of the vector of Kuhn-Tucker multipliers is known to lead to numerical instabilities in some situations. We show next that the risk of instabilities is actually limited, because the solution set is bounded, as will be deduced from the following lemma.

Lemma 1 We consider

q ∈ Q 0 , μ = (μ ij ) ∈ R r
+ , and we define

F = G q (μ) = 1≤i<j ≤N μ ij G ij (q) ∈ R dN .
The set

F = λ = (λ ij ) ∈ R r + , G q (λ) = F is bounded.
Proof Let us first establish the uniqueness for the homogeneous problem. For

q ∈ Q 0 , let λ = (λ ij ) ∈ R r + be such that G q (λ) = 1≤i<j ≤N λ ij G ij (q) = 0.
Let i 0 denote the index of an extremal vertex of the convex hull conv{q i , 1 ≤ i ≤ N}. By Hahn-Banach's theorem, the compact {q i 0 } and the closed convex set conv{q i , 1 ≤ i ≤ N, i = i 0 } can be separated in a strict sense by a plane in R d . We denote by x an element of this plane, and by v a normal vector to it. One has

q i 0 -x • v > 0 , q j -x • v < 0 ∀j = 1, . . . , N , j = i 0 ,
so that (q i 0q j ) • v > 0 for j = i 0 . Now the balance of contact forces exerted upon sphere i 0 in the direction v reads

j =i 0 λ ji 0 e ji 0 (q) • v = j =i 0 λ ji 0 q i 0 -q j q i 0 -q j • v
which is positive unless λ ji 0 = 0 for all j = i 0 . Therefore all multipliers associated to a contact with sphere i 0 are equal to 0, and this approach can be iterated for the reduced family (q j , j = i 0 ). By downward induction on the number of active spheres, we establish that 0 is actually reduced to {0}.

As for the non-homogeneous problem, F is obviously convex. The asymptotic cone of F is defined as (see e.g. Bourbaki [START_REF] Bourbaki | Sticky particles and scalar conservation laws[END_REF])

C F = s>0 s( F -λ),
where λ is any element of F (for example μ). It can be checked that (the proof, which is elementary, can be found in Maury [START_REF] Maury | Analyse fonctionnelle, exercices et problèmes corrigés[END_REF]), in a finite-dimensional space, the asymptotic cone of a convex set is reduced to {0} if and only if the set is bounded. If F is not bounded, then its asymptotic cone C F contains a half line R + ξ , with ξ ∈ R r + , which yields μ + R + ξ ⊂ F , and consequently ξ is a solution to the homogeneous problem, so that ξ = 0. Therefore, F is bounded.

Proposition 5 At each time step of Scheme (S), the solution set for

μ n+1 h is bounded.
Proof This is a direct consequence of lemma 1, with q = q n h and

hF = M(u n+1 h -u n h ) -hf n+1 h (q n h ).

Convergence theorem (case of a single contact)

We establish in this section a convergence result for the proposed time-stepping scheme, in the case of a single constraint (two spheres) in R 3 . As uniqueness does not hold, neither does continuity with respect to initial conditions. Thus we cannot expect any error estimate for the scheme (S). Moreover, the necessity to extract a subsequence to establish convergence is not purely technical: different subsequences may actually converge to distinct solutions to the same problem (see Section 6.1).

In order to alleviate notations, we consider the system of two spheres with unit mass and same radius r, but this assumption is not needed in the proof. The result remains valid for the case of two different spheres, or a sphere and an obstacle. The constrained minimization problem is again supposed to be solved exactly at each time step.

Two-body system

We consider here the mechanical sytem of two identical, hard spheres, described by

Q = {q = (q 1 , q 2 ) ∈ R 3 × R 3 },
and the associated feasible set

Q 0 = {q ∈ Q , D(q) ≥ 0},
where D(q) = |q 2 -q 1 |-2r is the distance between the two bodies. As previously, the tangent space at any q ∈ Q is denoted by T Q = R 6 . The convex cone N q is now

N q = {0} if D(q) > 0, -R + G = {-λG , λ ∈ R + } if D(q) = 0
where G = ∇D is the gradient of D. In order to avoid special care of the initial instant, we shall pick an initial configuration q 0 such that D(q 0 ) > 0.

We may now state the single-contact problem: given a time interval I = (0, T ), a force field (q, t) → f(q, t) ∈ R 6 verifying the regularity conditions ( 6)-( 7) and the initial conditions

(q 0 , u 0 ) ∈ Q 0 × T Q with D(q 0 ) > 0, find (q, u, μ) ∈ W 1,1 × BV × M 1 + such that u(0) = u 0 , (26) q(t) = q 0 + t 0 u(s) ds ∈ Q 0 ∀t ∈ I, (27) u(t) = f(q(t), t) + μ G(q(t)), (28) supp(μ) ⊂ {t , D(q(t)) = 0}, ( 29
) u + = u --P N q u -∀t ∈ I, ( 30
)
where ( 28) is understood in the sense of distributions, and P N q is the Euclidean projection onto the closed convex cone N q .

Let us now rewrite the overall time stepping-scheme in the present case of a single contact. The discretization spaces for pathlines (X h ), velocities (V h ), and reactions (R h ), are

X h = {q h = (q i ) : I → Q , q i ∈ (P 1 h ) 3 , i = 1, 2}, V h = {u h = (u i ) : I → T Q , u i ∈ (P 0 h ) 3 , i = 1, 2}, R h = μ h : I → R , μ h ∈ P 0 h
The approximated force field f h (q) ∈ V h is defined for any q ∈ Q by

t ∈ [(n -1)h, nh) -→ f n h (q) = 1 h nh (n-1)h f(q, s) ds.
For a given h > 0, we shall consider

(q h , u h , μ h ) ∈ X h × V h × R h
built according to Scheme (S):

1. Initialization

(q 0 h , u 0 h ) = (q 0 , u 0 ). ( 31
)
2. Time step

u n+1 h = u n h + hf n+1 h (q n h ) + hμ n+1 h G(q n h ) (32)
where u n+1 h is the solution to the constrained minimization problem min

u∈K h (q n h ) |u -u n h -hf n+1 h (q n h )| 2 (33) with K h (q n h ) = {u ∈ T Q , D(q n h ) + hG(q n h ) • u ≥ 0}, ( 34
)
and μ n+1 h ∈ R + is the associated Kuhn-Tucker multiplier.

Update the positions

q n+1 h = q n h + hu n+1 h . ( 35
)

Convergence

Theorem 1 Let (q h , u h , μ h ) h be a sequence of solutions obtained by Scheme (S) (equations ( 31) to ( 35)), with h → 0. Then there exists a subsequence of time steps (still denoted by h), and

(q, u, μ) ∈ W 1,1 × BV × M 1 + such that q h -→ q in W 1,1 , μ h μ in M 1 ,
and (q, u, μ) is a solution to problem ( 26)-( 30).

Proof The proof will be decomposed into 9 steps, which we shall briefly describe below. The critical ones are steps 3 and 9.

1. The scheme produces feasible configurations only, i.e.,

q h (t) ∈ Q 0 ∀h, t.

The family

(u h ) is uniformly bounded, i.e., ∃C ∞ , |u h (t)| ≤ C ∞ ∀t ∈ [0, T ] , ∀h > 0.
3. The fields u h have uniform bounded variation, i.e.,

∃C var , var(u h ) = N n=1 |u n h -u n-1 h | ≤ C var ∀h.

The familly (u h

) is relatively compact in L 1 (I ). One can extract a subsequence (still denoted q h ) such that q h -→ q in W 1,1 . The limit velocity u = q = lim qh is in BV, and the limit motion q is feasible. 5. The sequence (μ h ) h is bounded in L 1 , so that one can extract a subsequence which converges weak-to a vector-valued bounded measure μ ∈ M 1 . 6. The pair (u, μ) verifies (28). 7. The complementarity slackness condition (29) holds true. 8. The initial condition (26) is verified. 9. The jump equation (30) is verified.

Step 1 (Feasibility) This is exactly Proposition 2, which was established for any number of contacts:

D(q h (t)) ≥ 0 ∀h, t ∈ [0, T ].
Step 2 (u h is uniformly bounded) By Proposition 3,

|u h (t)| M ≤ |u 0 | M + 1 √ α T 0 F (t) dt ≤ C ∞ .
Step 3 (u h has uniformly bounded variation)

We shall first establish a lemma, which can be seen as a mono-dimensional version of the property.

Lemma 2 We consider sequences (u h ) h , (ρ h ) h , and

(g h ) h in P 0 h such that u n+1 h = hg n h + ρ n h u n h , ρ n h ∈ [0, 1] ∀n , ∀h, with u 0 h = u 0 .
We suppose furthermore that there exists a constant C ≥ 0 such that

h N-1 n=0 g n h ≤ C ∀h. ( 36
)
Then the familly (u h ) has uniform bounded variation:

∃C var , var(u h ) = N-1 n=0 u n+1 h -u n h ≤ C var ∀h.
Proof A straightforward calculation leads to

u n+1 h = ρ 0 h ρ 1 h . . . ρ n h u 0 +hg n h + hρ n h g n-1 h + hρ n h ρ n-1 h g n-2 h + • • • + hρ n h ρ n-1 h . . . ρ 1 h g 0 h = ρ 0 h ρ 1 h . . . ρ n h u 0 + n k=0 hλ k h g k h ,
with λ k h ∈ [0, 1] for any h, k. Therefore function u h can be written as a sum of the homogeneous solution v h defined by

v n+1 h = ρ n h v n h = ρ n h ρ n-1 h . . . ρ 0 h u 0

and a function w h determined by

w n+1 h = n k=0 hλ k h g k h .
The first one v h is monotonous and keeps a constant sign, so that var(v h ) ≤ |u 0 |.

The second one verifies var(w

h ) = N -1 n=0 w n+1 h -w n h = N-1 n=0 hλ n h g n h ≤ C, so that finally var(u h ) ≤ |u 0 | + C,
which ends the proof of Lemma 2.

Remark now that, as u n+1 h is defined as the Euclidean projection of

u n h + hf n+1 h (q n h ) onto the half space of T Q K h (q n h ) = {u ∈ T Q , D(q n h ) + h G(q n h ) • u ≥ 0}, whose boundary is a hyperplane of T Q with normal G = G/ |G|, one has u n+1 h -(u n h + hf n+1 h (q n h )) = α G
where α ≥ 0 is some real parameter. If the constraint is non-active (i.e. u n h + hf n+1 h (q n h ) ∈ K h (q n h )), then α = 0. Otherwise, as D(q n h ) ≥ 0,

u n h + hf n+1 h (q n h ) • G ≤ 0, so that α G can be written α G = -λ n h ((u n h + hf n+1 h (q n h )) • G) G,
where λ n h is a non-negative real number. Let us show that λ n h ≤ 1. Again, as D(q n h ) is non-negative, 0 lies in K h (q n h ), therefore

u n+1 h -(u n h + hf n+1 h (q n h )), u n+1 h -0 ≤ 0, which yields λ n h (λ n h -1) (u n h + hf n+1 h (q n h )) • G 2 ≤ 0,
and therefore λ n h ≤ 1. We established that

u n+1 h -(u n h + hf n+1 h (q n h )) = -λ n h ((u n h + hf n+1 h (q n h )) • G) G, with λ n h ∈ [0, 1].

Let us now introduce the Euclidean decomposition

u n h = U n h G(q n h ) + v n h ,
where v n h is the projection of u n h onto G⊥ . We verify now that U n h meets the assumptions of Lemma 2 :

U n+1 h = u n+1 h • G q n+1 h = u n+1 h • G q n h + u n+1 h • G q n+1 h -G q n h = (1 -λ n h ) u n h • G q n h + (1 -λ n h ) h f n+1 h (q n h ) • G q n h +u n+1 h • G q n+1 h -G q n h = ρ n h U n h + ρ n h h f n+1 h q n h • G q n h +u n+1 h • G q n+1 h -G q n h , with ρ n h = 1 -λ n h . As q n+1 h = q n h + hu n+1 h , with u n+1 h ≤ C ∞ , and as G is k-Lipschitz in Q 0 , one has u n+1 h • G(q n+1 h ) -G(q n h ) ≤ u n+1 h k C ∞ h ≤ k C 2 ∞ h. Finally, hf n+1 h (q n h ) • G(q n h ) + u n+1 h
• G(q n+1 h ) -G(q n h ) can be written as hg n h , and the corresponding family of piecewise constant functions (g h ) verifies (36). As ρ n h ∈ [0, 1], all assumptions of Lemma 2 are met, which implies that var(U h ) is bounded uniformly in h.

We denote by P n h the projection onto G(q n h ) ⊥ , so that v n h = P n h u n h . We conclude this step by writing

var(u h ) = N -1 n=0 u n+1 h -u n h ≤ N -1 n=0 (u n+1 h -u n h ) • G(q n h ) + N-1 n=0 P n h (u n+1 h -u n h ) ≤ N -1 n=0 U n+1 h -U n h + N-1 n=0 u n+1 h • ( G(q n+1 h ) -G(q n h )) + N -1 n=0 hP n h f n+1 h (q n h ) because P n h G(q n h ) = 0 ≤ var(U h ) + T k C 2 ∞ + ||F|| 1 .
Step 4 (q h converges to a feasible trajectory)

As (u h ), which is bounded in L ∞ , is now proved to have uniform bounded variations over I , it is relatively compact in L 1 , so that one can extract a subsequence, still denoted by (u h ), which converges in L 1 to u. As u h is the Sobolev derivative of q h (in W 1,1 ), then q h converges to a q in W 1,1 , and q = u. The limit velocity u is in BV. Indeed, for any ϕ ∈ (C 1 c (I )) 6 ,

I u • ϕ = lim h→0 I u h • ϕ ≤ sup h var(u h )||ϕ|| ∞ < +∞.
Finally, as the convergence q h -→ q is uniform (W 1,1 → L ∞ ), and D(q h (t)) ≥ 0 for every h, t, with D continuous with respect to q, then D(q(t)) ≥ 0 for every t ∈ I .

Step 5 (μ h converges weak-to μ)

For any h,

|G| I |μ h | ≤ N-1 n=0 u n+1 h -u n h + N-1 n=0 h f n+1 h (q n h ) ≤ var(u h ) + ||F|| 1 ,
so that (μ h ) is bounded in L 1 (|G| is the modulus of G(q h ), which is a positive constant). Therefore, by the De La Vallée Poussin compactness criterion, one can extract a subsequence which converges weakly-to μ in M 1 .

Step 6 (The balance of momentum is verified by the limits u and μ)

Let us first remark that the derivative of u h can be defined3 in the sense of distributions as

uh = N-1 n=1 (u n+1 h -u n h ) • δ n ,
where δ n is the Dirac mass at t n = nh. For any ϕ ∈ (D(I )) 6 , it comes

I uh • ϕ = N-1 n=1 (u n+1 h -u n h ) • ϕ(t n ) = - N n=1 u n h • (ϕ(t n ) -ϕ(t n-1 )) = - N n=1 hu n h • (ϕ (t n ) + o(1)) -→ - I u(t) • ϕ (t) dt (37)
when h goes to 0. The integral I uh • ϕ may also be expressed

I uh • ϕ = N -1 n=1 (hf n+1 h (q n h ) + hμ n+1 h G(q n h )) • ϕ(t n ) = I (f h (q h ) + μ h G h ) • ϕ h dt
where ϕ h ∈ V h is defined by taking the left value of ϕ on each subinterval [t n , t n+1 ), and similarly q h stands for the function of V h which is equal to q n h on [t n , t n+1 ). As f is uniformly Lipschitz with respect to q, and as q h converges to q in L ∞ , f h (q) converges to t → f(t, q(t)) in L 1 . Similarly ϕ h converges to ϕ in L ∞ , so that

I f h (q h ) • ϕ h h→0 -→ I f(t, q(t)) • ϕ(t) dt.
For the second term, we remark that (μ h ) can be seen as a bounded sequence of linear functionals on L ∞ . According to this remark, one may write (G • q designs t -→ G(q(t)), which is the uniform limit of G

• q h ) I μ h G h • ϕ h dt = < μ h , G h • ϕ h > = < μ h , G • q • ϕ > M 1 ,C 0 (I ) + < μ h , G • q • ϕ -G h • ϕ h > (L ∞ ) ,L ∞ -→ < μ, G • q • ϕ > .

Finally we have

I uh • ϕ -→ I f(t, q(t)) • ϕ(t) dt+ < μ, G • q • ϕ > . ( 38
)
From ( 37) and (38), it follows

- I u(t) • ϕ (t) dt = I f(t, q(t)) • ϕ(t) dt+ < μ, G • q • ϕ >
for all ϕ ∈ D(I ): this is equation (28) in the sense of distributions.

Step 7 (The limit reaction field μ is non-active when there is no contact)

For any t 0 ∈ I such that D(q(t 0 )) > 0, there exists η such that

D(q(t)) ≥ a > 0 ∀t ∈ (t 0 -η, t 0 + η),
so that, by uniform convergence of q h to q, there exists h 1 such that

D(q h (t)) ≥ a/2 > 0 ∀t ∈]t 0 -η, t 0 + η[ ∀h ≤ h 1 .
Since the velocities u h are uniformly bounded in L ∞ , there exists h 2 ≤ h 1 such that, for any h ≤ h 2 , D(q h ) + hG(q h ) • u h ≥ a/4 > 0, so that the numerical constraint D + hG • u ≥ 0 is not activated in (t 0η, t 0 + η).

Consequently, for h ≤ h 2 , μ h ≡ 0 in (t 0η, t 0 + η), so that supp(μ) ⊂ (t 0η, t 0 + η) C , which yields the property.

Step 8 (The initial condition is verified)

As q 0 is chosen apart from the boundary of Q 0 and thanks again to the uniform boundedness of u h in L ∞ , there exists η such that μ h vanishes in (0, η) for all values of h. It holds then Step 9 (The collision law is verified)

For t 0 such that D(q(t 0 )) > 0, we already established (see Step 7) that μ ≡ 0 in a neighbourhood of t 0 , so that u is absolutely continuous in this neighbourhood (as it solves (28)), hence continuous at t 0 :

u + = u -= u --P {0} u -.
The interesting case is of course D(q(t 0 )) = 0. We first establish from (28) that the discontinuity occurs along direction G : Lemma 3 At any contact time t 0 there exists λ ∈ R such that

u + (t 0 ) -u -(t 0 ) = λ G(q(t 0 )),
where G is the normalized gradient G/ |G| introduced in step 3.

Proof For η > 0, and H ∈ T Q such that H • G = 0, we consider the test function ψ η H where ψ η is continuous, piecewise affine, zero outside [t 0η, t 0 + η], taking values 0, 1, 0 at t 0η, t 0 , t 0 + η, respectively4 . Now writing (28) against ψ η H yields

- I u • Hψ η = I ψ η H • f + μ, G • q • Hψ η .
The left-hand side can be expressed

- I u • Hψ η = 1 η t 0 +η t 0 u • H - t 0 t 0 -η u • H = u + (t 0 ) -u -(t 0 ) • H + o(η).
The force integral I ψ η H • f goes to zero with η. The last term can be bounded

μ, G • q • Hψ η ≤ sup t∈[t 0 -η,t 0 +η] |G • q(t) • H| μ, ψ η ,
where the sup goes to zero because G • q(t) tends to G(q(t 0 )) as η tends to 0, with H • G(q(t 0 )) = 0. So finally the jump (u +u -) • H is zero for any H orthogonal to G, which ends the proof of Lemma 3.

Unless explicitly mentioned, fields are now taken at t 0 , where t 0 is a time at which contact occurs, and G designs G(q(t 0 )). We first check that u -• G ≤ 0. If otherwise, then D(q) = 0 implies that

D(q(t 0 -ε)) = -εG • u -+ o(ε) is negative for some ε > 0, meaning that q(t 0 -ε) / ∈ Q 0 , which is impossible. So u -• G ≤ 0. Let us show that λ ≥ -u -• G. If not, then u + • G = u -• G + λ G • G < 0,
which implies similarly that q enters the forbidden domain Q C 0 right after t 0 . Therefore

λ ≥ -u -• G. (39)
Let us now notice that the condition (30) characterizes the velocity after the collision as the solution to a constrained minimization problem : it minimizes uu - over C q , the polar cone of N q , which is the closed half space

{u ∈ T Q , u • G ≥ 0}.
As u -• G ≤ 0, the solution to that problem is clearly

ũ = u --( G • u -) G,
which we have to identify with u + to prove that the collision law (30) holds at time t 0 . As we know that u + is u -+ λ G with condition (39) on λ, we just have to check that

u + • G ≤ ũ • G = 0. (40)
We extend notations of Step 3 by introducing U ∈ L 1 (which can be seen as the velocity in the moving G direction) as the limit of u h • G(q h ). Now we recall the property we want to establish : at any t 0 such that D(q(t 0 )) = 0, equation (40) is verified. As t → G(q(t)) is continuous, it can be formulated in terms of U : at any contact time, U + ≤ 0. We propose actually to establish the property (which is stronger as U -≤ 0 at any contact time) :

U -≤ 0 ⇒ U + ≤ 0. ( 41
)
Proof of (41) : We first recall an important property of U h = u h • G(q h ) which was established in Step 3 : there exists sequences g n h and ρ n h such that

U n+1 h = hg n h + ρ n h U n h ,
with ρ n h ∈ [0, 1] and |g h | ≤ g for some g ∈ L 1 . By summing up over time steps between t and t + τ , it comes

U h (t + τ ) = ε(h, t, τ ) + ρ(h, t, τ )U h (t), ( 42
) with |ε(h, t, τ )| ≤ t+τ t-h g(s) ds (43)
and ρ(h, t, τ ) is in [0, 1]. Let t 0 be such that U -(t 0 ) ≤ 0. Then, for any ε > 0, there exists η 1 such that

U (t) ≤ ε ∀t ∈ (t 0 -η 1 , t 0 ).
Now, as g ∈ L 1 (I ), from (43) there exists η 2 ≤ η 1 such that, for any t ∈ (t 0 -η 2 , t 0 ),

|ε(h, t, τ )| ≤ ε ∀τ ∈ (0, η 2 ] ∀h ∈ (0, η 2 ). As U h converges in L 1 (t 0 -η 2 , t 0 ) to U , U h (t) ≤ U (t) + ε ≤ 2ε ∀t ∈ (t 0 -η 2 , t 0 ) \ A h , (44)
for some measurable set A h ⊂ (t 0 -η 2 , t 0 ) with |A h | → 0. Using (42) to "translate" approximately inequality (44), it comes

U h (t) ≤ 2ρε + ε ∀t ∈ (t 0 , t 0 + η 2 ) \ A h ∀h ∈ (0, η 2 ), with A h = |A h |.
The last inequality and L 1 convergence of U h to U (now considered in the translated interval (t 0 , t 0 + η 2 )) imply

U (t) ≤ 2ρε + ε a.e. in (t 0 , t 0 + η 2 ),
so that U + (t 0 ) ≤ 3ε, for any ε > 0. Therefore we have

U + (t 0 ) ≤ 0.
This completes the proof.

6 Numerical experiments

Numerical non-uniqueness

This first set of results is somewhat distant from the original purpose of this work, but it illustrates an interesting property of Scheme (S), which is the capability to compute multiple solutions. We consider the situation of a single material point moving on the real line, subject to the constraint q(t) ≥ 0, with an inelastic impact law at 0 + . The system reads :

u(0) = 0, (45)
q(t) = t 0 u(s) ds ∀t ∈ I, (46)
u(t) = f (t) + μ G(q(t)), (47) supp(μ) ⊂ {t , q(t) = 0}, (48) u + = u --P N q u -∀t ∈ I, (49)

where N q is {0} whenever q > 0, and R -as soon as q = 0. We shall consider a force field f which does not depend explicitly on q, and which is defined in the time interval I = (0, 4) by (see Figure 2) We chose α = 1.8. When the time step h is a negative power of 2, the right-hand side involved in the time-stepping process

f (t) = 1 for t ∈ 1 2 k+1 , 1 2 k+1 + 1 2 k+2 -α for t ∈ 1 2 k+1 + 1 2 k+2 , 1 2 k ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ k ∈ Z , k ≥ -4.
f n+1 h = (n+1)h nh f (t) dt
can be computed exactly. We present two solutions which have been computed with h = 2 -8 and h = 2 -9 , respectively. Those approximations correspond to distinct piecewise polynomial functions, which can be shown to be both solutions of the problem. One can observe that for time steps of the type 2 -2n , we have convergence to a solution (Figure 3), and for h = 2 -2(n+1) , we have convergence to the other solution (which corresponds to Figure 4). Other sequences of time steps were not investigated (f -integral can no longer be computed exactly).

Sticky particles

This second set of computations illustrates the good behaviour of the scheme when the time step is large. The bodies we consider are punctual masses moving on a line, with no external force, so that when two particles meet, they collide and behave like a single particle. We consider the following situation: the initial condition is q 0 = (0, a, 2a, . . . , 1a, 1) with a = 1 N -1 , the initial velocity vector u 0 is chosen arbitrarily in [-1, 1] N , and particles move according to system (8)- [START_REF] Hu | Direct simulation of flows of solid-liquid mixtures[END_REF].

Remark 12 This one-dimensional model is now commonly referred to as sticky particle model (although particles do not stick in a strict sense, as they could be pulled apart with infinitesimal forces). It is used in Brenier [4] to establish existence of solutions to the pressureless gas model same time, so that the individual behaviour of particles is poorly described. Nevertheless, all these computations give exactly the same final state (up to errors due to the constrained minimization procedure). This behaviour is actually a direct consequence of Proposition 4 and Remark 11. Indeed, in this one-dimensional model, the final state (when all the masses are stuck together) is completely determined by the position of the aggregate and its velocity. As the center of mass moves with a constant velocity, the scheme is exact for its trajectory so that, no matter how large the time step is (as soon as all collisions have been taken into account), the approximated final state is exact.

∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 ) = 0.

Many-body computation

Numerical tests have been performed satisfactorily for large numbers of spheres (up to 50 000). We present here its application to a polydisperse collection of 2000 spheres. For readability reasons, we chose a situation such that the spheres remain Fig. 8 Step 0 Fig. 9 Step 3 Fig. 10 Step 6

Fig. 11 Step 9

Fig. 1

 1 Fig. 1 Notations

|u h (t) -u 0 | ≤ t+h 0 F 0 F

 00 (s) ds, so that |u(t)u 0 | ≤ t (s) ds a.e. in (0, η), which implies u(0) = u 0 .

Fig. 2 Fig. 3 Fig. 4

 234 Fig. 2 Force field t → f (t)

Figures 5 ,Fig. 5 Fig. 6 Fig. 7

 5567 Figures 5,[START_REF] Coddington | Theory of ordinary differential equations[END_REF] and 7 represent the computed pathlines of the particles in the time interval [0, 3], for a "small" time step h = 0.01, and two larger time steps h = 0.5 and h = 1.5. For h = 0.01, all collisions are distinctly captured by the scheme. For larger time steps, several collisions are taken into account by the scheme at the

Borrowing the terminology of numerical methods to handle discontinuity (in space) of solutions to hyperbolic PDE's, one may say that collisions (i.e. discontinuities in time for the velocities) are captured, whereas in most other methods, they are tracked.

N(N -1)/2 is the number of possible contacts.

This notation means: for anyv ∈ T Q , ϕ ∈ (D(I )) 6 , < v • δ n , ϕ >= v • ϕ(t n ).

This "hat" function is obviously not in D(I ), but it can be approximated by regular functions in norm W 1,1 , which is enough here as u ∈ L ∞ . Therefore it can be used as a test function.

in a fixed plane of R 3 , so that they behave like discs in R 2 , but the algorithm can be applied indifferently to genuine 3D problems. Radii are distributed in the interval [0.8r 0 , 1.2r 0 ], with r 0 = 1.5 × 10 -2 , and the mass m i of sphere i is proportional to its volume.

As initial condition, we consider a cluster of spheres with uniform downward velocity U 0 = -e z . We simulate the impact of this cluster with a rigid obstacle defined as the half space R × R × R -. To illustrate the stability of the scheme, we run a case where the time step h = 0.1 is such that the displacement of a sphere during one time step is typically 3 times its own size. Figures 8,9, 10, and 11 represent the 2000 spheres at steps 0, 3, 6, and 9.

Conclusion

We presented a scheme to compute the motion of rigid bodies with nonelastic impact law. Because of its stability properties, its robustness (it produces feasible configurations only, even for large time steps), this scheme is particularly suitable to control the minimal distance between rigid particles in the context of fluid-particle simulations, with a controlled influence of the perturbation on the energy.

Beside this particular purpose, the fact that it can be used to handle situations where uniqueness does not hold, and its special behaviour in the case of simultaneous or quasi-simultaneous collisions even for large time steps, should make it an efficient tool to model more general granular flows.


Editorial (J Clin Neurophysiol. 2005 Oct;22(5):285-7) "Neurocomputational models in the study of epileptic phenomena" Fabrice Wendling INSERM, U642, Rennes, F-35000, France Université de Rennes 1, LTSI, Campus de Beaulieu, Rennes, F-35042, France.

Epilepsy is a complex dynamical disease characterized by recurring seizures. Although we know that seizures are related to abnormal excessive firing and synchronization of neurons in brain networks, the cause or causes of many forms of epilepsy are still unknown.

The need to better understand basic mechanisms that lead to epileptogenesis and ictogenesis has led to the development of a variety of reliable recording techniques of brain activity.

Electrophysiology still occupies an important place since it allows the activity of neurons or networks of neurons to be recorded with high time resolution. Functional imaging is also a rapidly developing field which now provides a way to acquire high-resolution spatial maps of epileptic activity. In addition to new experimental techniques, both in animal (in vivo and in vitro) and in human, in recent years novel theoretical approaches have enriched this research field since it is becoming clear that in order to integrate the plethora of experimental data available it is necessary to develop computer models with the purpose of better understanding the dynamics of neuronal networks leading to epilepsy. Although these models are used to investigate epileptic phenomena at different levels of complexity (from cellular to network levels) within specific situations, they share the common feature of providing insight into mechanisms involved in the generation of epileptic activity either at a microscopic (subcellular to cellular) or a macroscopic (multicellular to system) level.

Many experimental models were developed in the past decades, much before computer models. In vitro models include brain-slice preparations as well as whole-structure (mainly the hippocampus and associated structures) or whole-brain preparations. The experimental situation allows exposition of neuronal tissue to drugs capable of provoking convulsions that act specifically on a precise target, and recording of epileptiform activity such as spikes, transient bursts or seizures. In vitro models have been widely used in cellular neurophysiology and allowed to identify some key features of epileptogenic circuits such as the existence of recurrent excitatory connections between pyramidal cells. In parallel herewith several in vivo models have been developed. A question that is often raised -and not yet solved -is whether animal epilepsies resemble human epilepsies in meaningful ways? Apart from the difficulty related to homology, most people agree on the fact that experimental models have helped a lot in our understanding of the mechanisms leading to epileptic activity. Models of chronic focal epilepsy include electrical kindling or the application of toxic compounds, such as tetanus toxin, kainate, pilocarpine, that may lead to recurrent seizures usually after a latent period of several weeks. In several of these cases an initial status epilepticus is induced. These experimental models allow developmental aspects of epileptogenesis to be studied even if clinical relevance is mainly related to the issue of seizures induced by abnormal already established changes of local brain tissue.

Computational neuroscience is an interdisciplinary field of research connected to neuroscience, applied mathematics, physics and computer sciences. This discipline discusses neurophysiogically or neurobiologically relevant mathematical models and simulation methods that contribute to our understanding of neural mechanisms. It is a rapidly growing field, mainly because of the necessity of integrating structural, functional and dynamic approaches to study the brain. This issue of the Journal of Clinical Neurophysiology focus on computational models developed in the field of epilepsy. These models are complementary to aforementioned experimental models. They can be viewed as a way to integrate new and detailed knowledge coming from neurobiological research to explain experimental findings and to generate experimentally testable hypotheses about possible mechanisms lying at the sub-cellular, neuronal or population level depending on the scale of representation and on the nature of the experimental data under study. By formalizing and relating information across multiple levels of analysis, computational models have this unique integrative property of establishing links between successive levels of reduction. For instance, from the modeling of a network of interconnected detailed neurons and interneurons, computational models can help to simultaneously study mechanisms lying at the cell level and at the network level and to simulate corresponding activities that can be experimentally recorded (using multi-unit or field potential electrodes).

For the study of epileptic phenomena, two complementary approaches were developed over the past decades.

In the first one, referred to as "detailed", single neurons are accurately modeled regarding their structural components (dendrites, soma, axon) and functional properties (voltagedependent channels with kinetics derived from experimental voltage-clamp studies). Neuronal networks are then built from the interconnection of a relatively large number (i.e. several thousands) of principal neurons and different types of interneurons with appropriate synaptic interactions. In these networks, the summated post-synaptic potentials of pyramidal cell HAL author manuscript inserm-00147366, version 1 membranes that correspond to the field activity, can be studied as a function of various parameters such as the types of neurons introduced in the network, network size, connectivity patterns, and conduction delays. According to this approach that was extensively developed by Traub and collaborators since the early 80's, both spatial and temporal properties of the activity in the simulated patch of neuronal tissue are represented. Insights about tissue excitability, role of interneurons, factors leading to hypersynchronization have been put in evidence and confronted to real observations. For instance, combined experimental and theoretical work showed that some features are necessary for an epileptic discharge to occur: the population of neurons must be large enough and, inside this population, excitatory pyramidal neurons must be connected in a synaptic network in which the synapses need to have a sufficiently high probability of driving their targets above threshold. The use of such models has explained some basic mechanisms by which synchronized activity emerges. In particular, these models are able to generate activity patterns, in realistic networks, that closely mimic epileptic activity recorded in vitro as well as different types of electroencephalographic (EEG) activity seen in patients.

According to the second approach, referred to as "macroscopic" or "lumped", a higher level of organization, i.e. the neuronal population level, is considered. It starts from the fact that neurons form populations and that the EEG is a reflection of ensemble dynamics rising from macroscopic statistical interactions between interconnected neuronal sub-populations (pyramidal cells and interneurons). This approach was initially proposed by Wilson and Cowan and by Freeman and co-workers who made substantial progress in the understanding of perceptual processing in the olfactory system. Their studies, spanning the three last decades, are based on experimental data and on computational models in which the dynamics of each neuronal subset are simply represented by a 2 nd order ordinary differential equation As is clear from these brief sections, basic mechanisms of epileptic activity are an active area of research in which computational models play a major integrative role. More recent in the advances also demonstrate that computational neuroscience is becoming an accepted tool in the neurosciences with possible clinical applications. The intent of this special issue on "Neurocomputational models in the study of epileptic phenomena" is precisely to demonstrate that computational model-based methods constitute a way to bridge the gap between basic neuroscience research and clinical research in epilepsy at both microscopic and macroscopic levels. The five studies reported in this issue give a good overview of the use of physiologically-relevant models (following either the first or the second approach described in the above section) with respect to the interpretation of clinical data recorded in patients with epilepsy (partial or absence) or in patients presenting seizures under anesthesia or in the analysis of these data in the context of seizure prediction.

Suffczynski and colleagues (Amsterdam, The Netherlands and Warsaw, Poland) investigate the mechanisms of transition between normal EEG activity (i.e. non-epileptiform) and epileptiform paroxysmal activity using a computational model of thalamocortical circuits based on relevant patho(physiological) data. They compare model predictions with experimental results from different types of epilepsy, including patients with absence seizures. Chernihovskyi and colleagues (Cuernavaca, Mexico and Bonn, Germany) describe a very original approach which complements the topics of the four other contributions. Indeed, while the four other articles mainly present computational models aimed at interpreting mechanisms involved in the generation of real interictal/ictal observations (either at a microscopic or a macroscopic level), this one uses biologically inspired models to process EEG signals.

Authors elaborated a nonlinear, excitable, spatially extended medium composed of coupled model neurons in which EEG recordings are applied as local perturbations that induce transient changes in the dynamics of the perturbed system. Implemented in Cellular Neural Networks (CNN), excitable media are then used to approximate the degree of synchronization in EEG recordings. This approach provides a new direction in signal processing techniques and could help to progress on difficult questions such as the prediction of epileptic seizures.

Traub and colleagues (Brooklyn, Philadelphia, USA and Leeds, U.K.) contributed to this special issue by presenting the benefit of combining experimental with modeling studies (microscopic approach) to progress in the understanding of cellular and network mechanisms of epileptogenesis in vitro and in vivo. Based on their pioneered and recent works, they review some of the means by which cellular intrinsic properties, synaptic interactions, and electrical coupling via gap junctions contribute to the anomalous population activities characteristic of seizures. Their detailed models permit to make precise hypotheses about the neurobiological substrate of certain types of oscillations (observed at the electrophysiological level) that often precede the onset of seizures (very fast oscillations, > 80 Hz). This combined experimental/model approach leads to new clinical perspectives like the development of compounds that would specifically block putative axonal gap junctions which are showed to be critical for epileptogenesis.

Finally, the fifth article by Wendling and collaborators (Rennes, Marseille, France) describes a model-based approach used to study some pathophysiological mechanisms in human partial epilepsy. In this approach, electrophysiological patterns (intracerebral EEG) recorded from the hippocampus in patients with mesial temporal lobe epilepsy are related to mechanisms involved in the transition to seizure through a macroscopic physiologically-relevant model of EEG activity. Model parameters (related to synaptic interactions between subsets of neurons and interneurons represented in the model) are automatically identified using a procedure that minimizes a distance between real and simulated signals. Analysis of identified parameters suggest that the transition from interictal to ictal activity can not be simply explained by an increase in excitation and a decrease in inhibition but rather by time-varying ensemble interactions between pyramidal cells and local interneurons projecting to either their dendritic or perisomatic region (with slow and fast GABA A kinetics). Such approaches open new

  HAL author manuscript inserm-00147366, version 1 having a static nonlinearity identified as a asymmetric sigmoid curve. In this way a model of olfactory area able to produce EEG signals that approximate experimentally recorded EEGs quite accurately was developed. Similar ideas developed at the same time by Lopes da Silva and collaborators led to the development of a lumped-parameter population model able to explain the alpha rhythm of the EEG and more recently to the development of a model of the thalamocortical loop providing useful insights into the interpretation of mechanisms involved in absence epilepsy.

  Their approach allows them to derive hypotheses about physiological mechanisms of ictal transitions. They also reach conclusions about the predictability of absence seizures and open perspectives about the use of electrical stimulations to abort paroxysmal oscillations.Liley and Bojak (Hawthorn, Australia) present recent approaches to theoretically describing the electroencephalographic effects of proconvulsant volatile general anesthetic agents. Using biologically plausible computational and mathematical models, they managed to simulate realistic EEG signals in response to specific and quantifiable physiological changes related to different anesthetic agents. This original study deals with the crucial issue of mechanisms involved in ictogenesis and provides valuable insights into the physiological workings of seizure initiation.

Journal of Clinical Neurophysiology 10/2005; 22(5): 285-7

HAL author manuscript inserm-00147366, version 1