Bernard Rousselet

Periodic solutions of o.d.e. systems with a lipchitz non linearity

In this report, we address differential systems with Lipschitz non linearities; this study is motivated by the subject of vibrations of structures with unilateral springs or non linear stress-strain law close to the linear case. We consider existence and solution with fixed point methods; this method is constructive and provides a numerical algorithm which is under study. We describe the method for a static case example and we address periodic solutions of differential systems arising in the vibration of structures.

Introduction

This work stems from non destructive testing experiments reported in [START_REF] Vanderborck | Application of non-linear ultrasonic spectroscopy to health monitoring and damage detection in structures,. 38p[END_REF] and from the experiments and the computations of a beam model of the dynamic of satellite solar arrays 1 [HFR09], [START_REF] Hazim | Vibrations of a beam with a unilateral spring.Periodic solutions -Nonlinear normal modes[END_REF]. In these two situations the structure is submitted to an harmonic force; the behavior of the structure is locally non linear in particular due to unilateral elastic contacts.

During these studies, it became clear that an in-depth understanding of the forced dynamic behavior relies on the study of periodic solution of the free system; this can be considered as an extension of the use of normal modes of a linear free system in order to study the dynamics of the forced associated linear system as explained in [START_REF] Géradin | Théorie des vibrations. Application la dynamique des structures[END_REF] (in English see [START_REF] Géradin | Mechanical vibrations : theory and application to structural dynamics[END_REF]). Periodic solutions for non linear o.d.e. arising in structural dynamics have been considered in [BB05, BGKM01, JR09, KPGV09, JPS04, ABBC06, Mik10] and in many others.

In [START_REF] Junca | The method of strained coordinates for vibrations with weak unilateral springs[END_REF], the Lindstedt-Poincarré method was used in order to derive approximate non linear normal modes (a periodic solution close to a linear normal mode) for small non linearity.

At this point it is worth to recall a fundamental remark of Henri Poincarré in his report "Sur le problème des trois corps et les équations de la dynamique": "Ce qui rend ces solutions périodiques aussi précieuses, c'est qu'elle sont, pour ainsi dire, la seule brèche par où nous puissions pénétrer dans une place jusqu'ici réputée inabordable" ( [START_REF] Poincarré | Méthodes nouvelles de la mécanique céleste[END_REF]).

In English "What renders these periodic solutions so precious is that they are, so to speak, the only breach through which we may try to penetrate a stronghold previously reputed to be impregnable" Although this remark is mainly devoted to celestial mechanic, it is also meaning-full for vibration of structures.

In this report we study periodic solutions of non linear systems of ODE arising in free vibrations of structures subjected to unilateral springs; the ideas and the methods may be used for systems modeling other situations in particular when the non linearity is only Lipschitz. We recall a simple modeling process in section 2. Then in section 3, we start by presenting some ideas of quasi Newton type algorithm for solving equilibrium of structures in the static case with some Lipschitz non linearities; for example, these non linearities arise in cases where structures are subjected to unilateral springs; these unilateral springs are often simplified models of structures like bumpers usually made of viscoelastic materials like in [START_REF] Hazim | Numerical and experimental study for a beam system with local unilateral contact modeling satellite solar arrays[END_REF]; the same model will be addressed in the dynamic case. We derive a constructive proof of existence of periodic solutions for a one degree of freedom system in section 4; then the case of several degrees of freedom is addressed in section 5.

Vibrations with non linearities have been considered from an experimental point of view in [DGL + 03, VL04, DET02, MCG02, ZFGC05]; asymptotic expansions and numerical methods have been used in [RV04a, RV06, RV04b, RV05, VR05, RV06, Rou08, Rou10, JR09, BR09, Bra, JR10, Bra10].

The work of Lyapunov [START_REF] Lyapunov | The general problem of the stability of motion[END_REF] is often cited as a basis for the existence of periodic solutions which tends towards linear normal modes as amplitudes tend to zero; this result uses the hypothesis of analycity of the function involved in the differential system. As we address the case where the non linearity is only Lipschitz, this hypothesis of analyticity is obviously violated. Some general existence results, based on calculus of variations, for convex Hamiltonian systems are presented in [START_REF] Ekeland | Convexity methods in Hamiltonian mechanics[END_REF]; the case of analytic functions in connexion with the use of normal forms is considered in [START_REF] Jezequel | Analysis of non linear dynamic systems by the normal form theory[END_REF][START_REF] Touzé | Nonlinear normal modes for damped geometrically nonlinear systems: Application to reduced-order modelling of harmonically forced structures[END_REF][START_REF] Iooss | Approximate invariant manifolds up to exponentially small terms[END_REF]. Energy pumping has been addressed in several papers, for example in [START_REF] Schmidt | Energy pumping for mechanical systems involving non-smooth saint-venant terms[END_REF].

The case of vibration of structures with unilateral springs or more generally systems modeled with Lipschitz non linearities are considered in some recent studies (experiments, asymptotic expansions, numerical computations): [Haz, HR09a, HR09b, HR08, HFR09, Haz10, JR09, JR10, JPS04, VLP08]. A study of dynamics of elastic shocks is in [START_REF] Attouch | The dynamics of elastic shocks via epigraphical regularization of a differential inclusion. barrier and penalty approximations[END_REF]. The case of distributed systems modeled by partial differential equations has been addressed for example in [START_REF] Shaw | Normal modes of vibration for non-linear continuous systems[END_REF]; general references for perturbation methods are among many others [START_REF] Sanders | Averaging methods in nonlinear dynamical systems[END_REF][START_REF] Nayfeh | Introduction to perturbation techniques[END_REF] for differential equations and [START_REF] Kevorkian | Perturbation methods in applied mathematics[END_REF] for partial differential equations.

A review of periodic solutions of non autonomous ordinary differential equations may be found in [START_REF] Mawhin | Periodic solutions of non-autonomous ordinary differential equations[END_REF] for details, see also [START_REF] Rouche | Equations diffrentielles ordinaires. Tome II : Stabilité et solutions priodiques[END_REF].

In the static case, for example, a stability result concerning the obstacle problem for a plate is considered in [START_REF] Pozzolini | A stability result concerning the obstacle problem for a plate[END_REF].

Non smooth optimization is an active field of research; we only cite some references connected to this research: [START_REF] Asen | Implicit functions and solution mappings: : a view from variational analysis[END_REF], [START_REF] Kunisch | Semi-smooth newton methods for nondierentiable optimization problems[END_REF], [lei04] and the references therein.

Here we intend to use methods from non smooth optimization to prove existence of periodic solutions ans to derive a numerical algorithm to find them; the precise description and implementation of the algorithm will be addressed in a forthcoming paper.

Equations from structural mechanics

Consider a spring-mass system with some weak broken springs which are acting only in compression. Let u denote, the displacement of the masses, and γ j the strain of the springs; it is related to the displacement by γ j = u j+1 -u j or in vector form with an incident matrix B γ = Bu if we assume that the material is elastic linear, the stress-strain law is

σ = Eγ + ǫE ′ (γ + d) -
where ǫ is a small parameter; d denotes some backlash; in other words for some springs in small compression or in traction, there is no induced stress. In the static case, the force applied to each mass is denoted by Y and we have the equilibrium equation

B T σ = Y or more explicitely B T Eγ + ǫB T E ′ (γ + d) -= Y or B T EBu + ǫB T E ′ (Bu + d) -= Y (2.1) Remark 2.1.
• Obviously this system is quite general; many situations of structural mechanics may be cast in this system by using finite elements. Moreover many other physical systems can be modeled with such non linearities.

• To solve this system, the idea is to use some quasi Newton method taylored to systems involving Lipschitz functions.

• The proof of convergence relies on a fixed point method.

The principle of the method may be found in the book of Dontchev-Rockafellar [START_REF] Asen | Implicit functions and solution mappings: : a view from variational analysis[END_REF] as a way of proving the existence of an inverse of a Lipschitz function f or an implicit theorem for Lipschitz functions.

• Here the function f (X) is the left hand side and we want to solve f (X) = Y for some values of Y ; note that here f is only a Lipschitz function; it is not continuously differentiable!

Example with 5 masses and 4 springs on a straight line The incident matrix is

B =     1 -1 0 0 0 0 1 -1 0 0 0 0 1 -1 0 0 0 0 1 -1    
As we have enforced no boundary condition, there is a rigid body movement: u j = c for j = 1, . . . 5 gives the strain γ = 0; for example, we can enforce u 0 = 0 (the mass 0 is attached) and remove it from the degrees of freedom; so we get

B =     -1 0 0 0 1 -1 0 0 0 1 -1 0 0 0 1 -1     e e£ £ £ £ £ £ £ £ £ £ £ £ g g g g g g g g g e e z e e£ £ £ £ £ £ £ £ £ £ £ £ g g g g g g g g g e e z e e£ £ £ £ £ £ £ £ £ £ £ £ g g g g g g g g g e e z e e£ £ £ £ £ £ £ £ £ £ £ £ g g g g g g g g g e e z e e£ £ £ £ £ £ £ £ £ £ £ £ g g g g g g g g g e e z e e£ £ £ £ £ £ £ £ £ £ £ £ g g g g g g g g g e e
Figure 1: 5 masses and 6 springs broken at supports

In the case where the only broken spring is the second:

B T E ′ (Bu + d) -=     -E ′ 2 (u 3 -u 2 + d) - E ′ 2 (u 3 -u 2 + d) - 0 0    
Example with weak broken springs at supports on a straight line We consider the case where the springs are broken at both ends; then for n masses, we have n + 1 springs;

γ 1 = u 1 , γ j = u j -u j-1 , γ n+1 = -u n (2.2)
For n = 5, we get (see figure 2)

B =         1 0 0 0 0 -1 1 0 0 0 0 -1 1 0 0 0 0 -1 1 0 0 0 0 -1 1 0 0 0 0 -1        
In the case of broken springs at both ends,

σ 1 = E ′ 1 (γ 1 + d 1 ) -, σ j = E j γ j , for j = 2 . . . 5, σ 6 = E ′ 6 (γ 6 + d 6 ) - (2.3)
so that, the equilibrium equation are:

Ku + B T E ′ (Bu + d) -= F (2.4)
with K = B T EB, the rigidity matrix, and (2.5)

[B T E ′ (Bu + d) -] 1 = E ′ 1 (u 1 + d 1 ) - (2.6) [B T E ′ (Bu + d) -] j = 0, for j = 2, . . . 4, [B T E ′ (Bu + d) -] 5 = E ′ 6 (-u 5 + d 5 ) - (2.7) K =       E 2 -E 2 0 0 0 -E 2 E 2 + E 3 -E 2 0 0 0 -E 3 E 3 + E 4 -E 4 0 0 0 -E 4 E 4 + E 5 -E 5 0 0 0 -E 5 E 5       Remark 2.2.
1. We notice that K is not invertible; see Remark 3.2 2. Obviously, we can easily extend this example to a bi or tridimensional array of masses and springs, with unilateral springs or to non linear springs with lipschitz non linearity..

3 From a classical Newton method to quasi Newton method with a Lipschitz function

In this introductory section, we recall a constructive convergence proof of the Newton-Raphson algorithm based on the fixed point method; then we extend the approach to a quasi Newton method for a smooth function; finally, we extend it to a sum of a smooth and a lipschitz function.

3.1 Newton and Quasi Newton methods in the smooth case

Newton method

When f is differentiable, in order to solve an equation f (x) = y, the idea is to approximate f by its differential: we approach f (x k+1 ) by f (x k ) + f ′ (x k )(x k+1 -x k ) so that x k+1 is the solution of:

f (x k ) + f ′ (x k )(x k+1 -x k ) = y
To prove the convergence, we consider the equation for the previous index

f (x k-1 ) + f ′ (x k-1 )(x k -x k-1 ) = y
and manipulating, we obtain

f ′ (x k )(x k+1 -x k ) = -[f (x k ) -f (x k-1 ) -f ′ (x k-1 )(x k -x k-1 )] or (x k+1 -x k ) = -[f ′ (x k )] -1 [f (x k ) -f (x k-1 ) -f ′ (x k-1 )(x k -x k-1 )]
Now, it is clear that to show convergence of this sequence, we need mainly the Assumptions

• f ′ (x) -1 ≤ c in a neighborhood of the solution;

• and an easy consequence of differentiability of f in a neighborhood of x for any ǫ > 0, there exists δ such that for

x k -x ≤ δ, x k+1 -x ≤ δ f (x k ) -f (x k-1 ) -f ′ (x k-1 )(x k -x k-1 ) < ǫ x k -x k-1
we choose ǫ = 1 2c and we obtain

(x k+1 -x k ) ≤ 1 2 (x k -x k-1 )
• which proves convergence (if the initial data is close enough to the solution itself close enough to the given x).

• In this method, it is well known that the main trouble is the inequality

f ′ (x) -1 ≤ c
if c is large, we have to start very close to the solution in order to converge; the "tangent subspace" to the graph of f near the solution should not be "too horizontal".

A way to circumvent this trouble is to embed the Newton method in a continuation process.

Quasi Newton method

The idea of a quasi Newton method is to evaluate the derivative of f at a fixed point x.

Still considering the case where f is differentiable, instead of 3.1.1, we use

f (x k ) + f ′ (x)(x k+1 -x k ) = y
in which the derivative of f is evaluated at the fixed point x. To prove convergence, we consider the same equation for the previous index

f (x k-1 ) + f ′ (x)(x k -x k-1 ) = y
manipulating, we get:

f ′ (x)(x k+1 -x k ) = -[f (x k ) -f (x k-1 ) -f ′ (x)(x k -x k-1 )]
We use the same assumption on the inverse of the derivative f ′ and still with differentiability in the neighborhood of x, we get for any ǫ > 0

f (x k ) -f (x k-1 ) -f ′ (x)(x k -x k-1 ) < ǫ x k -x k-1 then, from (x k+1 -x k ) = -[f ′ (x)] -1 [f (x k ) -f (x k-1 ) -f ′ (x k-1 )(x k -x k-1 )]
with the same assumptions, we still get

(x k+1 -x k ) ≤ 1 2 (x k -x k-1 )
which proves convergence (if the initial data is close enough to the solution itself close enough to the given x).

Quasi Newton method in a Lipschitz case

We use some notions introduced in Dontchev-Rockafellar [START_REF] Asen | Implicit functions and solution mappings: : a view from variational analysis[END_REF].

We no longer assume that f is differentiable but that there exists a strict estimator h of f at x: Definition 3.1.

1. We recall that the Lipschitz modulus of e at x is lip(e, x) = lim sup

x,x ′ →x,x =x ′ e(x) -e(x ′ ) x -x ′ 2. h is a strict estimator of f at x, when f (x) = h(x) + e(x)
with the Lipschitz modulus of e satisfying lip(e, x) ≤ µ < +∞

We assume that ȳ = f (x)

Moreover, we assume that h is invertible around ȳ for x and lip(h

-1 , ȳ) ≤ κ (3.1)
Remark 3.1. It should be pointed out that it stems directly from the previous definition that the assumption, the Lipschitz modulus at x of a function be finite, implies that this function is Lipschitz in a neighborhood of x.

Examples

• f (x) = (x) -, lip(f, x) = 1 for x ≤ 0 but lip(f, x) = 0 for x > 0 • e(x) = x --αx with 0 ≤ α ≤ 1 lip(e, x) ≤ max(α, 1 -α) • e = (Bu + d) --Λ(Bu + d) where Λ is a diagonal matrix with 0 ≤ λ j ≤ 1. set d j (u) = B jk u k + d j as e j (u) -e j (u ′ ) = d j--λ j d j , we get |e j (u) -e j (u ′ )| ≤ max(Λ j , 1 -λ j )|d j (u) -d j (u ′ )| ≤ max(λ j , 1 -λ j ) B u -u ′ lip(e, x) ≤ max j max(λ j , 1 -λ j ) B
consider now the case with

f (u) = B T EBu + ǫB T E ′ (Bu + d) -. (3.2)
It is the left hand side of our example from structural mechanics; as the first term is linear, the function e is e = ǫ[(Bu

+ d) --Λ(Bu + d)] and lip(e, x) ≤ ǫ max j max(λ j , 1 -λ j ) B h(u) = B T EBu + ǫB T E ′ Λ(Bu + d) (3.3) is invertible when the rigidity matrix K ǫ = B T EB + ǫB T E ′ B is invertible and its inverse is: h -1 (y) = K -1 ǫ [y -B T E ′ Λd] and lip(h -1 , y) ≤ K -1 ǫ So when K ǫ is invertible, for ǫ small enough, lip(e, x)lip(h -1 , y) < 1
this property is crucial for convergence but it should be pointed out that K 0 = B T EB is not required to be invertible.

Algorithm

Now the idea is to build a sequence to solve f (x) = y by using the strict estimator h:

h(x k+1 ) + e(x k+1 ) = y
to compute x k+1 solution of this equation is as difficult as the initial problem but it may be approximated by replacing e(x k+1 ) by e(x k )

h(x k+1 ) + e(x k ) = y
to prove convergence, we use the formula at the previous iteration

h(x k ) + e(x k-1 ) = y (3.4) we get x k+1 = h -1 (y -e(x k )) x k = h -1 (y -e(x k-1 ))
and

x k+1 = x k + h -1 (y -e(x k )) -h -1 (y -e(x k-1 ))
• with the assumption (3.1), we get that for e(x k ) and e(x k-1 ) small enough and y close enough to ȳ,

h -1 [y -e(x k )] -h -1 [y -e(x k-1 )] ≤ λ -e(x k ) + e(x k-1 )
• on the other hand, h being a strict estimator:

e(x k ) -e(x k-1 ) ≤ ν x k -x k-1
• moreover, we assume that µ κ < 1 so that for x k and x k-1 close enough to x, we may assume that λ ν < 1

• From 3.2.2, we get

x k+1 -x k ≤ λ ν x k -x k-1
Usually, we have to show that if x 0 is close enough to x, x 1 lies in a ball in which h -1 is Lipschitz; here it is obvious; it proves the geometric convergence when

λ ν < 1
We have obtained Proposition 3.1. Consider the equation

f (x) = y (3.5) Assume that f (x) = h(x) + e(x)
with the Lipschitz modulus satisfying

lip(e, x) ≤ µ < +∞ lip(h -1 , ȳ) ≤ κ < +∞ with µκ < 1, f (x) = ȳ and e(x) = 0 then the sequence h(x k+1 ) + e(x k ) = y
with x 0 is given close to the solution, converges geometrically to the solution of (3.5) close to x 0 . This is the case for our example from structural mechanics when ǫ is small enough.

Continuation process This method may be embedded in a continuation process:

• Set a step size δǫ and an initial value of ǫ 1 and a starting point x1 :

• For l from 1 to L

• Use the non smooth quasi Newton method to find the solution x * l of the equation for ǫ l with the starting point xl ;

• increase ǫ l+1 = ǫ l + δǫ and set xl+1 = x * l • end For
Extensions for differential equations are considered below; control problems are possible extensions.

An example Consider the iterations (3.4) for the example of broken springs at both ends (2), we get with (3.2) and (3.3):

Ku k+1 + ǫB T E ′ Λ(Bu k+1 + d) = -ǫ B T E ′ (Bu k + d) --B T E ′ Λ(Bu k + d)) + F (3.6)
We notice that the matrix of the left hand side is the sum of K and ǫB T E ′ ΛB; it is usually invertible.

Remark 3.2. Had we use the "natural" iterative method

Ku k+1 = -ǫB T E ′ Λ(Bu k + d) + F (3.7)
we notice that the matrix K is not invertible so that this algorithm is not working at all; in cases where the smallest eigenvalue of K is not zero but small as in problem arising in large aerospace structures, the algorithm will converge for ǫ quite smaller than with the previous algorithm (3.6)

It could be proved that the matrix is better conditioned with the "quasi Newton" algorithm than with the natural algorithm.

4 Periodic solution of a one degree of freedom spring-mass system with Lipschitz non linearity

The differential equation, periodic solution

We set ẋ = dx dt (4.1)

and we consider a one d.o.f. spring-mass system with a linear and a weak linear unilateral spring:

m ẍ + kx + ǫk(x) -= 0 with x -= x -|x| 2 (4.2)
or more generally a linear and a weak non linear spring

m ẍ + kx + ǫkg(x) = 0 (4.3)
We set ω 2 = k m , so this equation may be written ẍ + ω 2 x + ǫω 2 (x) -= 0 (4.4)

As we suspect that the frequency of a periodic solution depends on ǫ, we perform the change of variable:

θ = ω ǫ t with 1 ω 2 ǫ = 1 -ǫη(ǫ) ω 2
Remark 4.1. This is a classical change of variable, e.g. see [START_REF] Verhulst | Nonlinear differential equations and dynamical systems[END_REF]; it is used in the method of strained coordinates also called Linsted-Poincarré; in particular, it has been used in [START_REF] Junca | The method of strained coordinates for vibrations with weak unilateral springs[END_REF] to derive an asymptotic expansion for a quite similar system of differential equations. Instead, here we use it to derive a constructive proof of existence of a periodic solution; we are only looking for periodic solutions whereas in [START_REF] Junca | The method of strained coordinates for vibrations with weak unilateral springs[END_REF], we can obtain a quasi periodic expansion.

and we set x ′ = dx dθ (4.5)

If we set x(θ) = x(t), equation (4.4) may be written:

x" + (1 -ǫη)x + ǫ(1 -ǫη)x -= 0 and in the general case (4.6)

x" + (1 -ǫη)x + ǫ(1 -ǫη)g(x) = 0 (4.7) or x" + x + ǫf (x, η, ǫ) = 0 (4.8) with f (x, η, ǫ) = -ηx + (1 -ǫη)
x -and in the general case (4.9)

f (x, η, ǫ) = -ηx + (1 -ǫη)g(x) (4.10)
As the solution of (4.8) is:

x = acos(θ) + bsin(θ) -ǫ θ 0 sin(θ -s)f (x(s), η, ǫ)ds (4.11) or x = acos(θ) + bsin(θ) -ǫr(θ, a, ǫ, η, b) (4.12)
we obtain the following lemma which may be considered as an extension to the Lipschitz case of a classical result which, for example, may be found in [START_REF] Verhulst | Nonlinear differential equations and dynamical systems[END_REF]. in both cases, we get from F 2 η(0) = 1 2 by using for example

| cos(s)| = 2 π - 4 π +∞ k=1 (-1) k 4k 2 -1 cos(2ks), (4.20)
and by anticipating that η is a Lipschitz function of ǫ (see proposition 4.2 )

1 ω 2 ǫ = 1 -ǫ/2 + o(ǫ) ω 2 and x = acos(θ) + bsin(θ) -ǫr(θ, a, 0, η(0), 0) + o(ǫ)
This result might be obtained by direct inspection of the level curves of the energy associated to the differential equation, see [START_REF] Junca | The method of strained coordinates for vibrations with weak unilateral springs[END_REF]; but this approach may be extended to systems of differential equations.

General case For simplicity, we consider b = 0, the equation F 1 (p, y) = 0 is identically satisfied. The second equation yields: Lemma 4.2. Assume that x is solution of (4.8) (4.10) with f and g Lipschitz with respect to all variables;

η(0) = 1 aπ 2π 0 cos (
|g(x 2 ) -g(x 1 )| ≤ k|x 2 -x 1 | (4.22) |f (x 2 , η 2 , ǫ 2 ) -f (x 1 , η 1 , ǫ 1 )| ≤ k (|x 2 -x 1 | + |η 2 -η 1 | + |ǫ 2 -ǫ 1 |) (4.23)
assume that the initial data are:

x(0) = a α , x ′ (0) = ǫb α , resp. x ′ (0) = b α with α = 1, 2 (4.24)
then, x and F are Lipschitz with respect to y = (η, b) and p with a modulus of magnitude ǫk resp. k (where k < +∞ ):

∀θ ∈ [0, 2π], x 2 (θ) -x 1 (θ) ≤ ǫk(|b 2 -b 1 | + |η 2 -η 1 |) and (4.25)
resp. when the initial velocity is not of order ǫ (4.26)

∀θ ∈ [0, 2π], x 2 (θ) -x 1 (θ) ≤ k(|b 2 -b 1 | + |η 2 -η 1 |) , and
(4.27) ∀θ ∈ [0, 2π], |g(x 2 (θ)) -g(x 1 (θ))| ≤ ǫk(|b 2 -b 1 | + |η 2 -η 1 |) (4.28)
but we have only (4.29)

|f (x(θ, y 2 ), η 2 , ǫ) -f (x(θ, y 1 ), η 1 , ǫ)| ≤ k (ǫ|b 2 -b 1 | + |η 2 -h 1 |) (4.30) F (p, y 2 ) -F (p, y 1 ) ≤ k(|b 2 -b 1 | + |η 2 -η 1 |) (4.31)
Proof. The proof relies on the formula (4.11), we get:

x 2 (θ) -x 1 (θ) = ǫ(b 2 -b 1 )sin(θ) -ǫ θ 0 sin(θ -s) [f (x 2 , η 2 , ǫ) -f (x 1 , η 1 , ǫ)] ds (4.32) |x 2 (θ) -x 1 (θ)| ≤ ǫ [|b 2 -b 1 | + 2πk|η 2 -η 1 |] + 2πǫk|x 2 -x 1 | (4.33)
from which we get equation of (4.25) and we deduce from (4.23)

|g(x(θ, y 2 ) -g(θ, y 1 )| ≤ ǫk(|b 2 -b 1 | + |η 2 -η 1 |) (4.34)
which is (4.28);but we get only

|f (x(θ, y 2 ), η 2 , ǫ) -f (x(θ, y 1 ), η 1 , ǫ)| ≤ k (ǫ|b 2 -b 1 | + |η 2 -h 1 |) (4.35)
we deduce equations (4.31).

A fixed point approximation method

Consider now the solution of (4.3) , with f given by (4.10); when g(acos(s)) is an even function, a 2π periodic solution with b = 0 satisfies trivially F 1 = 0; with lemma 4.2, we obtain: and for ǫ small enough, this function satisfies

|F(a, ǫ, η 2 ) -F(a, ǫ, η 1 ) ≤ ǫk|η 2 -η 1 | (4.38)
so that the following sequence

η k+1 = (1 -ǫη k )F(a, ǫ, η k ) (4.39)
converges and it proves the existence of a 2π periodic solution of (4.8) with (4.10), hence the existence of a periodic solution of (4.3) with angular frequency (4.4).

Estimator problem

We consider another approach which should be better conditioned (see Remark 3.2 ) and enable to consider arbitrary initial velocity whereas in the previous paragraph, the initial velocity is zero. Equations (4.13) may be considered as defining an implicit function p -→ y (4.40)

As in previous subsection, this function F is not smooth and we are going to show that the implicit function (4.40) may be defined with a fixed point of a contraction maping following general lines of [START_REF] Asen | Implicit functions and solution mappings: : a view from variational analysis[END_REF]. So we introduce a differential equation and a function H that will be proved to be an estimator of F . note that for ǫ = 0, h = f so that H = F .

Method: to solve F (p, y) = 0 following the method introduced in paragraph 3.2, we use the sequence y k , defined by

H(p, y k+1 ) = -E(p, y k ), where E(p, y) = F (p, y) -H(p, y) (4.45)
As H is a smooth map, we prove that H is invertible for ǫ small enough by using the classical implicit function theorem (see e.g. . Remark 4.2. With examples 3.2.1, when f is defined by (4.9) and h by (4.42) then f -h satisfies the general inequality (4.47).

Lemma 4.4. Assume that x is solution of (4.8) (4.10) and ξ solution of (4.41), (4.42) with the same initial conditions and with f and g Lipschitz with respect to all variables and moreover that h is a strict estimator of f , i.e.:

For small ξ, |f (ξ, η, ǫ) -h(ξ, η, ǫ)| ≤ µ|ξ| (4.47)

then |x(θ) -ξ(θ)| ≤ ǫµ c Sup 0≤s≤2π (|ξ|) (4.48) for θ ∈ [0, 2π]
Proof. Indeed, from (4.11) with x and a similar formula for ξ, both having the same initial conditions, we get

x(θ) -ξ(θ) = ǫ θ 0 sin(θ -s)[f (x, η, ǫ) -h(ξ, η, ǫ)]ds (4.49)
and so (4.50) with p 0 = (a, 0) T with a arbitrary, in other words, H is a strict estimator of F , uniformly in p for ǫ ≤ ǭ (following [START_REF] Asen | Implicit functions and solution mappings: : a view from variational analysis[END_REF]) .

x(θ) -ξ(θ) = ǫ θ 0 sin(θ -s)[f (x, η, ǫ) -f (ξ, η, ǫ) + f (ξ, η, ǫ) -h(ξ, η, ǫ)] (4
Proof. To prove it, we use the short hand f = f (x(θ, ỹ), η, ǫ), h = h(ξ(θ, ỹ), η, ǫ) where ỹ = (η, b) is close to y = (η, b). We have to estimate

E(p, ỹ) -E(p, y) = F (p, ỹ) -F (p, y) -(H(p, ỹ) -H(p, y)) (4.55)
where F is defined in (4.15),(4.16) and H in (4.43),(4.44), so the basic point is to consider:

f -f -( h -h) = -ηx + (1 -ǫη)g(x) -(-ηx + (1 -ǫη)g(x)) --η ξ + (1 -ǫη)g( ξ) + (-ηξ + (1 -ǫη)g(ξ)) (4.56)
Then we split the right hand side into 3 expressions that we manipulate separately

|η(x -ξ) -η(x -ξ)| =|η(x -ξ -x + ξ) + (η -η)(x -ξ)| ≤η |x -x| + | ξ -ξ| + |x -ξ||η -η|
with lemma 4.2 and 4.4

≤ǫc(1 + Sup|ξ|)|ỹ -y| (4.57) Similarly |g(x) -g(x) -g( ξ) + g(ξ)| ≤ k |x -x| + | ξ -ξ| ≤ ǫk|ỹ -y| (4.58) and ǫ -ηg(x) + ηg(x) + ηg( ξ) -ηg(ξ) = ǫ (η -η)g(x) + η(g(x) -g(x)) -(η -η)g(ξ) -η(g(ξ) -g( ξ) ] ≤ ǫ (Sup|g(x)| + Sup|g(ξ)|)|η -η| + -g(x)| + |g(ξ) -g( ξ)|) ≤ ǫc|ỹ -y| (4.59)
The last inequality is a consequence of lemma 4.2 and 4.4 So we have obtained

| f -f -( h -h)| ≤ ǫµ|ỹ -y| (4.60)
uniformly with respect to p = (a, ǫ) from which the lemma is obtained using the definitions of E(p, ỹ) and E(p, y) where F is defined in (4.15),(4.16) and H in (4.43),(4.44) .

Proposition 4.2. Equation (4.8), (4.9) (resp. (4.10) ) has a 2π periodic solution for ǫ small enough and it may be computed with the iterative process (4.45); moreover, η, b are Lipschitz function of a, ǫ and so |x -acos(θ) -bsin(θ) + ǫr(θ, a, 0, η(0), 0)| ≤ kǫ 2 (4.61)

Proof. The proof is now simple: write the iterative method (4.45) for k and k -1

H(p, y k+1 ) = -E(p, y k ) (4.62) H(p, y k ) = -E(p, y k-1 ) (4.63)
and by subtraction

H(p, y k+1 ) -H(p, y k ) = E(p, y k-1 ) -E(p, y k ) (4.64)
with previous lemma 4.3, there exists λ

y k+1 -y k ≤ λ E(p, y k-1 ) -E(p, y k ) (4.65)
and as E is a strict estimator (lemma 4.5)

y k+1 -y k ≤ ǫλµ y k -y k-1 (4.66)
As for the static case we should prove that for p close to p 0 and y 0 close to the solution obtained for ǫ = 0, y 1 lies in a ball where y(p, e) defined in lemma 4.3 is Lipschitz. The first part of the proposition is then proved by the fixed point theorem.

The fact that y is a Lipschitz function of p is derived with a somewhat classical argument. With lemma 4.3, we get from H(p 2 , y 2 ) + E(p 2 , y 2 ) = 0, H(p 1 , y 1 ) + E(p 1 , y 1 ) = 0 (4.67)

|y 2 -y 1 | ≤ λ|E(p 2 , y 2 ) -E(p 1 , y 1 )| (4.68)
and with lemma 4.5

|y 2 -y 1 | ≤ λ(k|p 2 -p 1 | + ǫµ|y 2 -y 1 |) (4.69) from which |y 2 -y 1 | ≤ λ 1 -ǫµ (k|p 2 -p 1 |) (4.70)
Then, as r is Lipschitz with respect to p and y, we obtain the inequality (4.61)

5 Vibrations with several degrees of freedom

Differential equation model of a spring-mass system

In the dynamic case, equation (2.1) becomes:

M ü + B T EBu + ǫB T E ′ (Bu + d) -= 0 (5.1) with u = du dt (5.2)
Remark 5.1. Many other mechanical systems with unilateral properties can be cast into this frame in particular when using finite elements see e.g. [START_REF] Jiang | Large-amplitude non-linear normal modes of piecewise linear systems[END_REF]; the case of beams with support on a unilateral spring is considered in H Hazim [START_REF] Hazim | Vibrations of a beam with a unilateral spring.Periodic solutions -Nonlinear normal modes[END_REF]. As for one degree of freedom, we can consider a more general case:

M ü + B T EBu + ǫG(u) = 0 (5.3)
with G(u) a Lipschitz function; as an introduction we shall also consider the case of G(u) a smooth function; it can model a non linear stress-strain law.

Eigenvectors, estimator system

In order to derive approximate periodic solutions, it is convenient to write the differential system in an eigenvector basis; we denote K = B T EB the rigidity matrix and introduce the matrix of generalized eigenvectors Φ = [Φ 1 , . . . , Φ j , . . . , Φ n ] where Φ k are generalized eigenvectors associated to generalized eigenvalues ω 2 k satisfying For unilateral spring, as γ -→ f (γ) = γ -is a concave, Lipschitz function, we can use as strict estimator (in the sense of [START_REF] Asen | Implicit functions and solution mappings: : a view from variational analysis[END_REF]) γ -→ h(γ) = Λγ with Λ = diag(λ j ), 0 < Λ j < 1; if we set e = f -h, we have:

-ω 2 k M Φ k + KΦ k = 0 (5.
lip(e, 0) = max j max(λ j , λ j -1) so a reasonable estimator differential system is ξ + Ω 2 ξ + ǫΦ T B T E ′ Λ(BΦξ + d) = 0
(5.7)

Periodic solutions

In the linear case when ǫ = 0, equation (5.1) has periodic solutions of the form:

u k = (υ k e iω k t + υ k e -iω k t )Φ k (5.8)
where Φ k are generalized eigenvectors defined in (5.4).

Orientation

We intend to prove that for ǫ small enough, there exists periodic solutions to equation (5.1), or in the general case (5.3); the method is constructive and will enable to derive a numerical scheme to approximate it; we use the equivalent system (5.5) or in the general case (5.6) and we introduce the following notations:

p = (a 1 , ǫ) T y = (η, b 1 , a 2 , b 2 , . . . , b n ) T (5.9)
Remark 5.2. Notice that a 1 plays a particular role: we are going to consider a periodic solution close to the harmonic solution of period 2π ω 1 of (5.5) with ǫ = 0; obviously, any other harmonic solution could be chosen by using a permutation of indexes.

For unilateral springs

Following the lines of the 1 d.o.f. case, we introduce a change of variable suited to look for periodic solutions close to the eigenmode of index 1. We set θ = ω ǫ t with

x(θ) = x(t), 1 ω 2 ǫ = 1 -ǫη(ǫ) ω 1 2
and set x ′ = dx dθ (5.10) Equation (5.5) becomes

x" 1 + x 1 -ǫηx 1 + ǫ 1 -ǫη ω 2 1 Φ T 1 B T E ′ (BΦx + d) -= 0
(5.11)

x j " + ω 2 j ω 2 1 x j -ǫη ω 2 j ω 2 1 x j + ǫ 1 -ǫη ω 2 1 Φ T j B T E ′ (BΦx + d) -= 0, j = 2, . . . , n (5.12)
for the estimator differential system,

ξ 1 " + ξ 1 -ǫηξ 1 + ǫ(1 -ǫη) ω 2 1 Φ T 1 B T EΦ T 1 B T Λ(BΦξ + d) = 0 (5.13) ξ j " + ω 2 j ω 2 1 ξ j -ǫη ω 2 j ω 2 1 ξ j + ǫ(1 -ǫη) ω 2 1 Φ T j B T EΛ(BΦξ + d) = 0, j = 2, . . . , n (5.14)
we set:

f 1 (x, η, ǫ) = -ηx 1 + 1 -ǫη ω 2 1 Φ T 1 B T E ′ (BΦx + d) - (5.15) f j (x, η, ǫ) = -η ω 2 j ω 2 1 x j + 1 -ǫη ω 2 1 Φ T j B T E ′ (BΦx + d) -= 0, j = 2, . . . , n (5.16
) so that the system may be written:

x" 1 + x 1 + ǫf 1 (x, η, ǫ) = 0 (5.17)

x" j + ω 2 j ω 2 1 x j + ǫf j (x, η, ǫ) = 0, j = 2, . . . , n (5.18)
similarly, we set:

h 1 (ξ, η, ǫ) = -ηξ 1 + 1 -ǫη ω 2 1 Φ T 1 B T E ′ Λ(BΦξ + d) (5.19) h j (ξ, η, ǫ) = -η ω 2 j ω 2 1 aξ j + 1 -ǫη ω 2 1 Φ T j B T E ′ Λ(BΦξ + d) = 0, j = 2, . . . , n (5.20)
and the "estimator" system is now:

ξ 1 " + ξ 1 + ǫh 1 (ξ, η, ǫ) = 0 (5.21) ξ j " + ω 2 j ω 2 1 ξ j + ǫh j (ξ, η, ǫ) = 0, j = 2, . . . , n (5.22)

General case

In the general case of equation (5.3), in the eigenvector basis (5.6), we still obtain equation (5.17),(5.18), but with

f 1 (x, η, ǫ) = -ηx 1 + 1 -ǫη ω 2 1 Φ T 1 G(Φx) (5.23) f j (x, η, ǫ) = -η ω 2 j ω 2 1 x j + 1 -ǫη ω 2 1 Φ T j G(Φx) = 0, j = 2, . . . , n (5.24)
in the non smooth case, we consider a strict estimator h j (x) of Φ T j G(Φx) and the estimator system is written as previously in (5.21),(5.22)

Condition of periodicity

Lemma 5.1. Consider the solution x of the Cauchy problem of system (5.5) or (5.6) and the solution x of the associated system (5.17), (5.18) after change of variable and with initial conditions: x n (0) -a n = 0, ẋn (0) -b n = 0 (5.27)

x 1 (0) -a 1 = 0, ẋ1 (0)
(5.28)

For ǫ close to zero, x is a solution of (5.5) of period 2π ωǫ with 1

ω 2 ǫ = 1-ǫη(ǫ) ω 1 2
if and only if

F (p, y) = 0 (5.29)
where (p, y) is defined in (5.9) and the function

F : R 2 × R 2n → R 2n (5.30)
is defined by:

F 1 (p, y) = 2π 0 sin(s)f 1 (x (s), η, ǫ)ds (5.31)
F 2 (p, y) = 2π 0 cos(s)f 1 (x(s), η, ǫ)ds (5.32) F 3 (p, y) =x 2 (2π) -a 2
(5.33) . . .

(5.34)

F 2 * n-1 (p, y) =x n (2π) -a n (5.35) F 2 * n (p, y) = ẋn (2π) -b n (5.36)
(5.37)

Proof. The proof is similar to the proof of lemma 4.1

In the non smooth case, it will be proved that the following function H is a strict estimator of F . Definition 5.1. Consider the solution ξ of the Cauchy problem of system (5.7) in the form (5.21), (5.22) with initial conditions:

ξ 1 (0) -a 1 = 0, ξ1 (0)

The smooth case

As indicated in the introduction, the analytic case is well known; here we just assume some differentiability; this type of result may be derived by other methods;we emphasize that this approach is constructive in particular if we approximate x with a Fourier series (harmonic balance or spectral method).

Lemma 5.2. In the case where G is smooth, the Jacobi matrix of y -→ F (p, y) (5.48) is invertible at p 0 = (a, 0) with b 1 = 0; where F is defined by (5.31), (5.37)

Proof. We get for ǫ = 0

x 1 = a 1 cos(θ) + b 1 sin(θ) (5.49)

x j = a j cos( ω j ω 1 θ) + b j ω 1 ω j sin( ω j ω 1 θ) (5.50) ∂x 1 ∂η = 0, ∂x 1 ∂b 1 = sin(θ) (5.51) ∂ ẋ1 ∂η = 0, ∂ ẋ1 ∂b 1 = cos(θ) (5.52) so that for b 1 = 0, ∂f 1 ∂η = -a 1 cos(θ) ∂f 1 ∂b 1 = -ηsin(θ) + 1 ω 2 1 Φ T 1 G ′ (Φx)Φ 1 sin(θ) (5.53) so that ∂F 1 ∂η = 0 ∂F 1 ∂b 1 = -ηπ + 1 ω 2 1 Φ T 1 G ′ (Φx)Φ 1 π (5.54) ∂F 2 ∂η = -a 1 π ∂F 2 ∂b 1 = 0 (5.55) det( ∂(F 1 , F 2 ) ∂(η, b 1 ) = 0 (5.56)
For other indexes

∂x j ∂a j = cos( ω j ω 1 θ), ∂x j ∂b j = ω 1 ω j sin( ω j ω 1 θ) (5.57) ∂ ẋj ∂a j = - ω j ω 1 θsin( ω j ω 1 θ), ∂ ẋj ∂b j = cos( ω j ω 1 θ) (5.58) ∂F 2j-1 ∂a j = ∂x j ∂a j (2π) -1, ∂F 2j-1 ∂b j = ∂x j ∂b j (2π)
(5.59)

∂F 2j ∂a j = ∂ ẋj ∂a j (2π), ∂F 2j ∂b j = ∂ ẋj ∂b j -1 (5.60) ∂F 2j-1 ∂a j = cos(2π ω j ω 1 ) -1 ∂F 2j-1 ∂b j = ω 1 ω j sin(2π ω j ω 1 ) (5.61) ∂F 2j ∂a j = - ω j ω 1 sin(2π ω j ω 1 ) ∂F 2j ∂b j = cos(2π ω j ω 1 ) -1 (5.62) When ω j ω 1 is not an integer for j > 1, the determinant det ∂(F 2 * j-1 , F 2 * j ) ∂(a ′j+1 j , b ′k+1 j ) = 0 (5.63)
ans so the Jacobi of F is invertible.

Proposition 5.1. Equation (5.31) to (5.37) has a solution for ǫ small enough and this solution may be computed with an iterative method.

Proof. We conclude the first part from the classical implicit function theorem; the solution may then be approximated by a Newton method; it may be embedded in a continuation process. at p 0 = (a 1 , 0) T with a 1 arbitrary, (following the general lines of [START_REF] Asen | Implicit functions and solution mappings: : a view from variational analysis[END_REF]) in other words, H is a strict estimator of F uniformly in p close to p 0 .

Strict estimator of F

Proof

The solution of (5.5) after change of variable θ = ω ǫ t, is solution of (5.11) and (5.12) or (5.17), (5.18); it may be written

x 1 =a 1 cos(θ) + b 1 sin(θ) -ǫ θ 0
sin((θ -s))f 1 (x(s), η, ǫ) ds (5.65)

x j =a j cos( ω j ω 1 θ) + b j ω 1 ω j sin( ω j ω 1 θ) -ǫ ω 1 ω j θ 0 sin( ω j ω 1 (θ -s))f j (x(s), η, ǫ) ds (5.66)
and the solution of (5.7) after change of variable is solution of (5.21) and (5.22) with h defined in (5.19), (5.20)

ξ 1 = a j cos(θ) + b j sin(θ) -ǫ θ 0 sin(θ -s)h 1 (ξ(s), η, ǫ) ds
(5.67)

ξ j = a j cos( ω j ω 1 θ) + b j ω 1 ω j sin( ω j ω 1 θ) -ǫ ω 1 ω j θ 0 sin( ω j ω 1 (θ -s))h j (ξ(s), η, ǫ) ds (5.68)
and so

E 1 (p, y) = 2π 0 sin(s) [f 1 (x(s), η, ǫ) -h 1 (ξ(s), η, ǫ)] ds (5.69) E 2 (p, y) = 2π 0 cos(s) [f 1 (x(s), η, ǫ) -h 1 (ξ(s), η, ǫ)] ds
(5.70)

E j (p, y) =x j (T ) -ξ j (T ) = -ǫ ω 1 ω j T 0 sin( ω j ω 1 (T -s)) [f j (x(s), η, ǫ) -h j (ξ(s), η, ǫ)] ds
(5.71)

As the solution of a system of differential equations is Lipschitz with respect to the initial conditions and with respect to time, we obtain the lemma from slight manipulations as in the proof of lemma 4.5; absolute values are replaced by norms and we use the following lemmas.

Some lemmas

Lemma 5.4. Assume that x is solution of (5.17) (5.18) with f and g Lipschitz with respect to all variables;

g(x) -g(x) ≤ k x -x (5.72) f (x, η, ǫ) -f (x, η, ǫ) ≤ k ( x -x + |η -η| + |ǫ -ǫ|) (5.73)
then, x and F are Lipschitz with respect to y = (η, b) and p; more precisely assume that the initial data are:

x(0) = a, x ′ (0) = b, (resp. x(0) = ǎ, x′ (0) = b) (5.74) then: ∀θ ∈ [0, 2π], x(θ) -x(θ) ≤ k( b -b + |η -η|) and (5.75) ∀θ ∈ [0, 2π], g(x(θ)) -g(x(θ)) ≤ k( b -b + |η -η|) (5.76)
and we have (5.77)

f (x(θ, y), η, ǫ) -f (x(θ, y), η, ǫ) ≤ k b -b + |η -η| (5.78) F (p, y) -F (p, y) ≤ k( b -b + |η -η|) (5.79)
Proof. The proof relies on the formula (5.65) , we get:

x1 (θ) -x 1 (θ) = ( b1 -b 1 )sin(θ) -ǫ θ 0 sin(θ -s) [f 1 (x, η, ǫ) -f 1 (x, η, ǫ)] ds (5.80) x1 (θ) -x 1 (θ) ≤ b1 -b 1 + 2πǫk|η -η| + 2πǫk Sup 0≤s≤2π x(s) -x(s) (5.81)
and from (5.66), we get

xj (θ) -x j (θ) = ω 1 ω j ( bj -b j )sin( ω j ω 1 θ) -ǫ θ 0 sin( ω j ω 1 (θ -s)) [f (x, η, ǫ) -f (x, η, ǫ)] ds (5.82) xj (θ) -x j (θ) ≤ ω 1 ω j bj -b j + 2πǫj|η -η| + 2πǫj Sup 0≤s≤2π x(s) -x(s) (5.83)
from which we get equation of (5.75) and we deduce from (5.73)

g(x(θ, y) -g(θ, y) ≤ k( b -b + |η -η|) (5.84)
which is (5.76); we get

f (x(θ, y), η, ǫ) -f (x(θ, y), η, ǫ)| ≤ k b -b + |η -h| (5.85)
we deduce equations (5.79).

Lemma 5.5. Assume that x is solution of (5.17) (5.18) and ξ solution of (5.21), (5.22), with the same initial conditions and with f and g Lipschitz with respect to all variables and moreover that h is a strict estimator of f , i.e.:

For small ξ, f (ξ, η, ǫ) -h(ξ, η, ǫ) ≤ µ ξ (5.86)

then x(θ) -ξ(θ) ≤ ǫµ c Sup 0≤s≤2π ( ξ ) (5.87) for θ ∈ [0, 2π]
Proof. The proof is quite similar to the 1 d.o.f. case of lemma 4.4

Implicit equation with the strict estimator H

In order to compute and to prove the existence of periodic solutions, we consider a constructive approach as suggested in the general lines of [START_REF] Asen | Implicit functions and solution mappings: : a view from variational analysis[END_REF].

Iterative method

Solve

H(p, y k+1 ) = -E(p, y k ) with p = [a 1 , ǫ] T , (5.88)

y k+1 = [η k+1 , b k+1 1 , a k+1 2 , . . . , b k+1 n ] T (5.89)
In other words, set x k (resp. ξ k ) the solution of (5.17) (5.18) (resp. (5.21), (5.22)) for the value y = y k ; similarly set f k = f (x k , η k , ǫ), h k = h(ξ k , η k , ǫ); we have to solve (5.7) for ǫ small enough and with initial conditions (5.38) where the initial value a 1 is prescribed and others are to be found as well as the frequency parameter η such that:

2π 0 sin(s)h k+1 1 = - 2π 0 sin(s)(f k 1 -h k 1 )ds (5.90) 2π 0 cos(s)h k+1 1 = - 2π 0 cos(s)(f k 1 -h k 1 )ds (5.91) ξ k+1 2 (2π) -a k+1 2 = -[x k 2 (2π) -ξ k 2 (2π)]
(5.92) . . .

(5.93)

ξ k+1 n (2π) -a k+1 n = -[x k n (2π) -ξ k n (2π)] (5.94) ξk+1 n (2π) -b k+1 n = -[ ẋk 1 (2π) -ξk n (2π)] (5.95)
As H is a strict estimator of F , we may assume that for ǫ small enough, the right hand side is as small as needed; we denote it α, β: is not an integer for j > 1, the jacobian matrix of y -→ H(p, y) (5.100) defined in (5.42) to (5.47) is invertible at p 0 = (a, 0) with b 1 = 0

α k 1 = 2π 0 sin(s)(f k 1 -h k 1 )ds, β k 1 = 2π 0 cos(s)(f k 1 -h k 1 )ds, (5.96) α k j = x k j (2π) -ξ k j (2π), β k j = ẋk j (2π) -ξk j (
Proof. We get for ǫ = 0 ξ 1 = a 1 cos(θ) + b 1 sin(θ) (5.101)

ξ j =

Computation and existence of periodic solutions

Proposition 5.2. Equation (5.31) to (5.37) has a solution for ǫ small enough and this solution may be computed with the iterative method (5.88).

Remark 5.3. The iterative method may be embeded in a continuation process by increasing the value of ǫ and as E is a strict estimator (lemma 5.64) y k+1 -y k ≤ ǫΛµ y k -y k-1 (5.121) the proposition is proved Proposition 5.3. Equation (5.17),(5.18) with f 1 , . . . f n given by (5.23),(5.24) has a 2π periodic solution for ǫ small enough and it may be computed with the iterative process (4.45); moreover, η, b are lipschitzian functions of p = (a 1 , ǫ)

Proof. The proof is similar to the one of Proposition 4.2

Remark 5.4. This result should be compared with Proposition 2.1 of [START_REF] Junca | The method of strained coordinates for vibrations with weak unilateral springs[END_REF]; in this paper we have an explcit approximate value of the frequency which gives an approximation of the solution in an interval [0, γǫ -1 ]; here we get the existence of a period which gives an approximation of the solution on [0, ∞]; but if this period is computed numerically, the approximation of the solution will only remain valid on a finite interval to be determined. The result of [START_REF] Junca | The method of strained coordinates for vibrations with weak unilateral springs[END_REF] may be viewed as a stability result: if the frequency is computed up to the order ǫ 2 , the approximation of the solution remains valid in an interval [0, γǫ -1 ].

Conclusion

This approach seems to be the first one to provide a rigorous and constructive proof of existence of periodic solutions of non smooth systems of differential equations arising in structural mechanics. Implementation of the algorithm with the harmonic balance principel is in progress. It paves the way to reduced order modelling (see e.g. [START_REF] Touzé | Nonlinear normal modes for damped geometrically nonlinear systems: Application to reduced-order modelling of harmonically forced structures[END_REF]) for vibrating structures with non smooth non linearities such as contact with unilateral springs modeling bumpers.

  Lemma 4.1. Consider the solution x of the Cauchy problem of equation (4.2) (resp. (4.3) ) and the solution x of the associated equation (4.8), (4.9) (resp. (4.10)) after change of variable. For ǫ close to zero, x is a solution of (4.2) (resp. (4.3) ) of period 2π/ω ǫ if and only if x is a solution of (4.8), (4.9) (resp. (4.10)) and r(2π, a, ǫ, η, b) = 0 and ∂r(2π, a, ǫ, η, b) ∂θ = 0 or equivalently F (p, y) = 0 (4.13) with p = (a, ǫ), y = (η, b) (4.14) F 1 (p, y) = 2π 0 sin(s)f (x, η, ǫ)ds (4.15) F 2 (p, y) = 2π 0 cos(s)f (x, η, ǫ)ds (4.16) 4.2 Computation of η(0) Case x -We have f (x, η, 0) = -ηx + x - and for b = 0, f (x, η, 0) = -ηacos(θ) + (acos(θ)) - s) 2 + cos(s)(acos(s))) -ds (4.18)we note that F 1 is identically zero; we remark that if a ≥ 0, (acos(s)) -= a(cos(s)) -and if a ≤ 0, (acos(s)) -= a(cos(s)) + (4.19)

  s)g(acos(s))ds (4.21) 4.3 Perturbation of the solution x with respect to y = (η, b)

Proposition 4. 1 .

 1 For a a 2π periodic solution with b = 0 satisfies η = (1 -ǫη)F(a, ǫ, η)

Definition 4. 1 .

 1 Consider the differential equationξ" + ξ + ǫh(ξ, η, ǫ) = 0 (4.41) with h(x, η, ǫ) = -ηξ + (1 -ǫη)αξ with 0 < α < 1 (4.42)and we introduce the following function:H 1 (p, y)

  Jacobi matrix of y -→ H(p, y) with p = (a, ǫ) and y = (η, b) at ǫ = 0 and b = 0 is 0 (-η(0) + α)π -aπ 0 2. so with α = η(0) = 1 2 and for ǫ small enough, the equation H(p, y) + e = 0 (4.46) has a solution y(p, e) which is differentiable and so locally Lipschitz with respect to its variables p, e for e small and p close to (a, 0). Proof. For ǫ = 0 and b = 0, we notice that ∂ξ ∂η = 0, ∂ξ ∂b = sin(θ), so ∂h ∂η = -acos(θ), ∂h ∂b = (-η + α)sin(θ), from which we deduce the first part of the lemma. The second part is deduced from the classical implicit function theorem.

0 (

 0 .51) with (4.47) and f Lipschitz, |x(θ) -ξ(θ)| ≤ ǫk θ |x -ξ| + µ|ξ|)ds (4.52) with Gronwall lemma (jl2), we get |x(θ) -ξ(θ)| ≤ ǫµSup(|ξ|)e 2πǫk (4.53) for θ ∈ [0, 2π] Lemma 4.5. Assume (4.47), then the function E(p,y)=F(p,y)-H(p,y), satisfies lip(E, y) ≤ ǭµ < +∞ (4.54)

  4) or in matrix form: KΦ = Ω 2 M Φ with Ω 2 the diagonal matrix of eigenvalues ; we assume Φ T M Φ = I consider the change of function u = Φx we obtain: ẍ + Ω 2 x + ǫΦ T B T E ′ (BΦx + d) -= 0 (5.5) and in the general case ẍ + Ω 2 x + ǫΦ T G(Φx) = 0 (5.6)

  ξ n (0) -a n = 0, ξ n (0) -b n = 0 (5.40)(5.41)the function H is defined by:H 1 (p, y) = 2π 0 sin(s)h 1 (ξ(s), η,ǫ)ds (5.42) H 2 (p, y) = 2π 0 cos(s)h 1 (ξ(s), η, ǫ)ds (5.43) H 3 (p, y) =ξ 2 (2π) -a 2 (5.44) . . . (5.45) H 2 * n-1 (p, y) =ξ n (2π) -a n (5.46) H 2 * n (p, y) = ξ n (2π) -b n (5.47)

Lemma 5. 3 .

 3 The function function E(p,y)=F(p,y)-H(p,y), satisfies lip(E, y) ≤ ǫµ < +∞ (5.64)

  2π) (5.97) Equations (5.90) to (5.95) with small right hand side denoted α, β are When the angular frequencies satisfy the property ω j ω 1

  det∂(H 2 * j-1 , H 2 * j ) Jacobian of H is invertible.By using the classical implicit theorem, we get Lemma 5.7. equation (5.88) or (5.90) to (5.95) define a lipschitzian implicit function:p = [a 1 , ǫ] T , [α, β] T → y k+1 = [T k+1 , b k+1 1 , a k+1 2 , . . . , b k+1 n ] T(5.116)

Proof.

  The proof is now simple write the iterative method (5.88) for k and k -1H(p, y k+1 ) = -E(p, y k ) (5.117) H(p, y k ) = -E(p, y k-1 ) (5.118)and by substractionH(p, y k+1 ) -H(p, y k ) = E(p, y k-1 ) -E(p, y k ) (5.119)with previous lemma, there exists Λy k+1 -y k ≤ Λ E(p, y k-1) -E(p, y k ) (5.120)