text:
\section{#1}\setcounter{equation}{0}

\title{   The Method of Strained Coordinates
          for Vibrations
         with Weak  Unilateral Springs

  }

  \author{St\'ephane Junca\thanks{ Universit\'{e} de Nice,  JAD laboratory, Parc Valrose, 06108 Nice, France,  junca unice.fr}
          \;  \&  \;
           Bernard Rousselet\thanks{ Universit\'{e} de Nice,  JAD laboratory, Parc Valrose, 06108 Nice, France, br unice.fr}
\\ \\
{\it \small  Laboratoire Dieudonn\'e, UMR CNRS 6621, Universit\'{e} de Nice Sophia-Antipolis}
           }

\date{\today}

   \begin{document}

 \maketitle

  \abstract{
 We study some spring mass  models for a structure
   having  a unilateral spring of small rigidity

    $\varepsilon$.
 We obtain and justify
  an asymptotic expansion with the method of strained coordinates
 with new tools to handle such defects,
 including a non negligible cumulative effect over a long time:
$T_\varepsilon  \sim \varepsilon^{-1}$ as usual;  or, for a new critical case,
 we can only expect:  $T_\varepsilon \sim \varepsilon^{-1/2}$.
We check numerically these results and present a purely numerical algorithm to compute ``Non linear Normal Modes'' (NNM); this algorithm provides results close to the asymptotic expansions but enables to compute NNM even when $\epsilon$ becomes larger.
  }  \bigskip \\

    %
   {\bf Keywords}:
   nonlinear vibrations,
  method of
 strained coordinates,
 piecewise linear,
 unilateral spring,
 approximate nonlinear normal mode.
             \\

 $\mbox{ }$\hspace{-6mm}
{\it Mathematics Subject Classification.}
  Primary:  34E15; \\
$\mbox{ }$ \hspace{53mm} Secondary: 26A16, 26A45, 41A80.

    \section{Introduction \label{SI}}

 For spring mass  models, the presence of
 a small piecewise linear rigidity can model a small defect
 which implies unilateral reactions on the structure.
 So, the  nonlinear and
 piecewise linear function $u_+=\max(0,u)$ plays a key role in this paper.
   For nondestructive testing we  study a  singular nonlinear effect
 for large time by asymptotic expansion of the vibrations.
 New features and comparisons with classical  cases
of smooth perturbations are given, for instance, with  the classical
 Duffing equation:
 $ \ddot{u}+u +\varepsilon u^3=0$ and the non classical case:
  $ \ddot{u}+u +\varepsilon u_+=0$.
Indeed, piecewise  linearity is singular:  nonlinear
 and Lipschitz but not differentiable.
 We give some new results to validate such asymptotic expansions.
 Furthermore, these tools are also valid
for a more general   non linearity.
 A nonlinear crack approach for elastic waves
 can be found in \cite{JuLoA}.
  Another approach   in the framework of  non-smooth analysis
 can be found in \cite{Att,Cla,PaoSch}.
\\
    For short time, a linearization procedure is enough
  to compute a good approximation.
   But for large time, nonlinear cumulative effects
  drastically alter the nature of the solution.
    We will consider the classical
  method of strained coordinates  to compute asymptotic expansions.
 The idea goes further back to Stokes, who in 1847 calculated
 periodic solutions for a weakly nonlinear wave propagation problem,
 see \cite{KC68,KC96,Mi,Nayfeh} for more details and references therein.
 Subsequent authors have generally referred to this as the
  Poincar\'e method  or the Lindstedt method.
 It is a simple and efficient method which gives us
  approximate nonlinear normal modes
 with $1$ or more degrees of freedom.
 \\
   In section \ref{s1dof} we present the method on
  an explicit case with  an internal Lipschitz force.
We  focus on an  equation with   one degree of freedom with expansions
 valid for time   of order $\varepsilon^{-1}$
 or, more surprisingly, $ \varepsilon^{-1/2}$ for a degenerate contact.
  \\
   Section \ref{sexp+} contains a tool to expand $(u+\varepsilon v)_+$
 and some accurate estimate for the remainder.
 This is a new key point  to validate
 the  method of strained coordinates with unilateral contact.
 \\
  In Section \ref{sNdof},
  we extend previous results to  systems
 with $N$ degrees of freedom,
first,
with the same accuracy  for approximate nonlinear normal modes,
 then, with less accuracy with all modes.
We check numerically these results and present a purely numerical algorithm to compute ``Non linear normal Modes'' (NNM) in the sense of Rosenberg \cite{Rosenberg}; see \cite{ABBC06} for two methods for the computation of NNM; see \cite{JPS} for a computation of non linear normal mode with unilateral contact and \cite{nnm-kpgv} for a synthesis on non linear normal modes;
 this algorithm provides results close to the asymptotic expansions but enables to compute NNM even when $\epsilon$ becomes larger.
\\
In Section \ref{scos}, we briefly  explain why we only perform expansions with
 even periodic functions to compute the nonlinear frequency shift.
\\
   Section \ref{sA} is an appendix
 containing some technical proofs and results.

    \section{One degree of freedom \label{s1dof}}

    \subsection{Explicit  angular frequency   \label{ss1dofexp}}

We consider a one degree of freedom spring-mass system (see figure \ref{Ressort}): one spring is  classical linear and attached to the mass and to a rigid wall, the second is still linear attached to a rigid wall  but has a unilateral contact with the mass; this is to be considered as a damaged spring.
The force acting on the mass is  $k_1 u + k_2 u_+ $ where $u$ is the displacement of the mass, $k_1$, the rigidity of the undamaged spring and $k_2$, the rigidity of the damaged unilateral spring. We notice that the
 term
\begin{eqnarray*}
 u_+ &
 =
 &
 \max(0,u).
\end{eqnarray*}
is Lipschitz but not differentiable. Assuming that $k_2=\varepsilon k_1$ and dividing by the value of the mass we can consider the equation:
\begin{eqnarray} \label{ode}
 \ddot{u} + \omega_0^2 u + \varepsilon u_+ =0,
\end{eqnarray}
\begin{figure}[!ht]
\begin{center}
\includegraphics[angle=0, width=11cm, height=3cm]{mass-ress-unilat1_6}
\caption{ Two springs,   on the right
    it has only  a unilateral contact. \label{Ressort}}
\end{center}
\end{figure}
 where  $\omega_0$  a positive constant.
This case has  a constant  energy $E$,
 where
 \begin{eqnarray*}  2 E&=& \dot{u}^2 +\omega_0^2 u^2 + \varepsilon(u_+)^2. \end{eqnarray*}
 Therefore, the level sets of $E(u,\dot{u})$ will
 be made of two half ellipses.
 Indeed, for $u<0$ the level set is an half ellipse,
  and for $u>0$ is another  half ellipse.
  Any solution $u(t)$  is confined to a closed level curve
  of $E(u,\dot{u})$ and is necessarily a periodic function of $t$.
 \\
  More precisely, a non trivial solution ($E>0$)
 is on the half  ellipse: $\dot{u}^2 + \omega_0^2 u =2 E$,
 in the phase plane
 during the time $T_C= \pi/\omega_0$,
 and on the half ellipse
   $\dot{u}^2 + (\omega_0^2+\varepsilon) u = 2E$
during the time
 $T_E= \pi/\sqrt{\omega_0^2+\varepsilon}$.
 The period $P(\varepsilon)$ is then
\\
 $ P(\varepsilon) = \displaystyle (1 + \left(1+\varepsilon/\omega_0^2\right)^{-1/2})
            \pi/\omega_0,$
 and the exact angular frequency  is:
\begin{eqnarray} \label{exactomega1}
 \displaystyle
\begin{array}{ccl}
  \omega(\varepsilon)& =&
         2 \omega_0 (1+ \left(1+\varepsilon/\omega_0^2\right)^{-1/2})^{-1}
                  \\
            &= &
   \displaystyle  \omega_0 + \frac{\varepsilon}{(4\omega_0)} -\frac{\varepsilon^2}{(8\omega_0^3)}
                 + \mathcal{O}(\varepsilon^3).
\end{array}
\end{eqnarray}
Let us compare  with
the angular frequency  for Duffing
 equation
 where  the nonlinear term is  $u^3$:
\begin{eqnarray} \label{ode-duffing}
 \ddot{u} + \omega_0^2 u + \varepsilon u^3 =0,
\end{eqnarray}
 which depends on the amplitude $a_0$
 of the solution ( see for example \cite{KC68,KC96,Mi,Nayfeh}):
\begin{eqnarray*}
  \omega_D(\varepsilon)  & = & \omega_0 + \frac{ 3}{8\omega_0^2} a_0^2\varepsilon
                    -\frac{ 15}{256\omega_0^4} a_0^4\varepsilon^2
                   + \mathcal{O}(\varepsilon^3).
\end{eqnarray*}

    \subsection{The method of strained coordinates
             \label{ss1dofmsc}}

 Now, we compute, with  the method of strained coordinates,
  $\omega_\varepsilon$, {\it an approximation}
  of the exact angular frequency  $\omega(\varepsilon)$ which is smooth
 with respect to $\varepsilon$ by exact formula
  (\ref{exactomega1}):

   \begin{eqnarray*}
       \omega(\varepsilon)&=&\omega_\varepsilon~+~\mbox{\rm O}(\varepsilon^3).
 \end{eqnarray*}
 We expose  this case completely to use the  same method
 of strained coordinates  later
    when we will not have such  an explicit formula.
\\
Let us define the new time and rewrite equation  (\ref{ode})
 with the new time
  \begin{eqnarray} \label{newtime} s  = \omega_\varepsilon t,
       \qquad  u_\varepsilon(t)   =  v_\varepsilon(s), \\
  \label{odes}
\omega_\varepsilon^2 v_\varepsilon''(s) + \omega_0^2 v_\varepsilon(s) + \varepsilon (v_\varepsilon(s))_+ =0,
 \end{eqnarray}
To simplify the exposition  and the computations,
  we take following initial conditions for $u_\varepsilon$
\begin{eqnarray}
  u_\varepsilon(0) = a_0 >0,\; \qquad \dot{u}_\varepsilon(0)&=&0,
\end{eqnarray}
i.e.  $v_\varepsilon(0)=a_0$ and $v_\varepsilon'(0)=0$.
 Similar computations are valid for negative $a_0$, see Proposition
 \ref{propexplicit} below.  With more general data,
i.e. when $\dot{u}_\varepsilon(0) \neq 0$, computations are more complicate
 and give the same angular frequency $\omega_\varepsilon$,
 see section \ref{scos}.
\\
 In new time $s$, we use  the following {\bf ansatz}
 \begin{eqnarray}  \label{anzatz}
   v_\varepsilon(s)
        &=&  v_0(s)+\varepsilon v_1(s ) +\varepsilon^2 r_\varepsilon(s).

 \end{eqnarray}
and the following notations:

\begin{eqnarray}
   \omega_\varepsilon& =& \omega_0+\varepsilon \omega_1+\varepsilon^2 \omega_2,

  \quad  \label{defalpha2}
      \omega_\varepsilon^2  =
        \alpha_0+\varepsilon \alpha_1+\varepsilon^2 \alpha_2 +\mathcal{O}(\varepsilon^3)

  \\
  \alpha_0 & =& \omega_0^2,\; \alpha_1 = 2\omega_0\omega_1 ,
  \quad \alpha_2 = \omega_1^2+2\omega_0\omega_2,
  \end{eqnarray}
 where $\omega_1, \omega_2$ or   $\alpha_1, \alpha_2$
  are unknown.
\\
We will   also use the  natural expansion,
$(u+\varepsilon v)_+ = u_+ + \varepsilon H(u)v +  \cdots$,
where $ H$ is the Heaviside function, equal to $1$ if $u>0$ and else $0$.
This expansion is validated in  Lemma \ref{LDLL} below.
\\
 Now, replacing  this  ansatz in  (\ref{odes})
   we obtain   differential
equations and initial data   for $v_0,v_1,r_\varepsilon$; set
\begin{eqnarray} \label{Ldefinition}
 & L(v) = -\alpha_0( v''  + v) &
\\

\label{eq0onedof}
 L(v_0) = & 0  ,
            &  v_0(0)=a_0,  v'_0(0)=0,\\
  \label{eq1onedof}
L( v_1)  = &  (v_0)_+ +  \alpha_1 v_0'',
               & v_1(0)=0, \; \; v'_1(0)=0, \\
 \label{rode}
L( r_\varepsilon)  = &
  H(v_0)v_1 + \alpha_2 v_0'' + \alpha_1 v_1''
 +  R_\varepsilon(s),
      &      r_\varepsilon(0)=0,\; \; r'_\varepsilon(0)=0.
\end{eqnarray}
  We now  compute,
  $\alpha_1$, $v_1$ and then $\alpha_2$.
  We have $v_0(s)= a_0 \cos(s)$. A key point in the method of
  strained coordinates is to keep bounded $v_1$ and $r_\varepsilon$ for large time
 by a  choice of $\alpha_1$ for $u_1$ and $\alpha_2$ for $r_\varepsilon$.
 For this purpose, we avoid resonant or {\bf secular} term
in the right-hand-side of
 equations (\ref{eq1onedof}), (\ref{rode}).
Let us first focus on $\alpha_1$. Notice that,
 $v_0(s)= a_0 \cos(s)$ and $a_0 > 0$, so
 \begin{eqnarray*}
      \displaystyle (v_0)_+ &  = & \displaystyle
   a_0 \left( \frac{\cos s}{2}+  \frac{|\cos s|}{2}\right).
\end{eqnarray*}
  $|\cos(s)|$ has  no term with frequencies $\pm 1$,
     since there are  only even frequencies.
 Thus $( (v_0)_+ - \alpha_1 v_0)=
  a_0\cos(s)(1/2 -\alpha_1) + a_0|\cos(s)|/2$
 has no secular term if and only if $\alpha_1 = 1/2$,
 so $ \omega_1= 1/(4\omega_0)$.
Now, $v_1$ satisfies:
\begin{eqnarray*}
 -\omega_0^2( v_1''  + v_1)  =   a_0|\cos s|/2 ,
 & \qquad v_1(0)=0, \quad v'_1(0)=0.
\end{eqnarray*}
 %
 To  remove secular term in the equation (\ref{rode}) we have to obtain
 the Fourier expansion for $H(v_0)$ and $v_1$.
 Some computations yield:
 \begin{eqnarray}
  \label{eqabscosF}  |\cos(s)| &= & \frac{2}{\pi} -
      \frac{4}{\pi} \sum_{k=1}^{+\infty} \frac{(-1)^k}{4k^2-1}\cos(2ks),\\
 \label{eqv1F}
    v_1(s) & = & \frac{-a_0}{\pi \omega_0^2}
      \left( 1 -
     2
  \sum_{k=1}^{+\infty}\frac{(-1)^k}{(4k^2-1)^2}\cos(2ks)  \right)
   - A \cos (s)
     ,
      \\
\label{eqHF}
    H(v_0) & = &  \displaystyle \frac{1}{2} +
      \frac{2}{\pi}  \sum_{k=1}^{+\infty} \frac{(-1)^j}{2j+1}\cos((2j+1)s)
      \label{H(v0)S},
  \end{eqnarray}
where
  $A = \displaystyle \frac{-a_0}{\pi \omega_0^2}
      \left( 1 -
     2
  \sum_{k=1}^{+\infty}\frac{(-1)^k}{(4k^2-1)^2}\right). $
\\
 To remove secular term of order one in (\ref{rode}), it suffices to take
 $\alpha_2$ such that:
\begin{eqnarray} \label{evalalpha2}
  0 & = & \displaystyle  \int_0^{2\pi}
   \left[H(v_0(s))v_1(s)+\alpha_{\bf 2}v_0''(s)+\alpha_1v_1''(s)
    \right]\centerdot   v_0(s) ds
 \end{eqnarray}
For Duffing equation, see \cite{KC68,KC96,Mi}, the source term  involves
 only few complex exponentials and the calculus of $\alpha_2$ is explicit.
 For general smooth source term, Fourier coefficients decay very fast.
 Here, we have an infinite set of frequencies for $v_1$ and $H(v_0)$,
  with only a small algebraic  rate of decay for Fourier coefficients.
 So, numerical computations need to compute a large number
 of   Fourier coefficients.
  For our first simple example, we can compute explicitly $\alpha_2$.
 After lengthy and  tedious computations involving numerical series,
 we have from (\ref{evalalpha2}) and
  (\ref{eqabscosF}), (\ref{eqv1F}), (\ref{eqHF})

 to evaluate exactly the following  numerical series:
\begin{eqnarray*}
        \alpha_2     &= &\displaystyle -\frac{2}{(\pi \omega_0)^{2}}
              \left(1 -\frac{1}{9} +\sum_{j=1}^{+\infty}
 {(2j+1)^{-1}}{(4j^2-1)^{-2}} - {(4(j+1)^2-1)^{-2}}  \right)
   \\ &=& -3(4\omega_0)^{-2},
 \end{eqnarray*}
thus $\omega_2 = -(2\omega_0)^{-3}$ as we  have yet obtained
 in (\ref{exactomega1}).
\\
\\
When $a_0=1$   we obtain in figure \ref{figspectra}  first modes
 of the infinite Fourier spectra  for
 \begin{eqnarray*}
  v_0(\omega_\varepsilon t) + \varepsilon v_1(\omega_\varepsilon t)  &\simeq &u_\varepsilon(t):
  \end{eqnarray*}

\begin{figure}[!ht]
\begin{center}
\includegraphics[angle=0, width=12cm, height=7cm]{eps01tmax099nbm819200}
\caption{ \label{figspectra}
 $v_1$,
 $ \displaystyle \log_{10} \left(\frac{|\widehat{u}_\varepsilon|}{\max|\widehat{u}_\varepsilon|} \right)  $
}
\end{center}
\end{figure}
 Indeed  we have for negative or positive $a_0$  the
 following result.
\begin{proposition} \label{propexplicit}

 Let $u_\varepsilon$ be  the solution of
 (\ref{ode}) with the   initial data:
    \begin{eqnarray*}
         u_\varepsilon(0) =a_0 + \varepsilon a_1,
    & \quad

    \dot{u}_\varepsilon(0)= 0,
 \end{eqnarray*}
 then  there exists $\gamma > 0$, such that,
  for all $t <T_\varepsilon=\gamma\varepsilon^{-1}$, we have the  following expansion:
\begin{eqnarray*}
  u_\varepsilon(t)   & = & v_0(\omega_\varepsilon t) + \varepsilon v_1( \omega_\varepsilon t)
                       + \mathcal{O}(\varepsilon^2 ),
\\
 v_0(\omega_\varepsilon t)  & =  & a_0 \cos(\omega_\varepsilon t), \\
v_1(\omega_\varepsilon t) & =&  \frac{-|a_0|}{\pi \omega_0^2}
      \left( 1 -
     2
  \sum_{k=1}^{+\infty}\frac{(-1)^k}{(4k^2-1)^2}\cos(2k\omega_\varepsilon t) \right)
   + (a_1-A) \cos(\omega_\varepsilon t),\\
  \displaystyle \omega_\varepsilon   & = & \omega_0  +\frac{1}{4 \omega_0} \varepsilon
     -\frac{1}{(2 \omega_0)^3}\varepsilon^2 ,
\end{eqnarray*}
where
  $A = \displaystyle \frac{-|a_0|}{\pi \omega_0^2}
      \left( 1 -
     2
  \sum_{k=1}^{+\infty}\frac{(-1)^k}{(4k^2-1)^2}\right). $
\\
\end{proposition}
In the Proposition \ref{propexplicit},
 with the  method of strained coordinates,
we recover  an asymptotic expansion
for the exact angular frequency
$ \omega (\varepsilon) = \omega_0 + \omega_1\varepsilon+\omega_2\varepsilon^2 + \mathcal{O}(\varepsilon^3)$
and for the exact solution $u_\varepsilon(t)$.
The term $ v_1$ is explicitly  given by its Fourier expansion.
Notice also that $\omega_2$ is not so easy to compute.
It needs to compute a numerical series.
\\
The technical proof of the   Proposition \ref{propexplicit}
is postponed to   the appendix.
\\
 Examples from  Proposition  \ref{propexplicit}
  have angular frequency  independent of the amplitude.
  Equation (\ref{ode}) is homogeneous .
  Indeed, it is a special   case, as we can see in the non homogeneous
 following cases. In these cases, we assume that the spring is either not in contact with the mass at rest ($b>0$)  or with a prestress at rest ($b<0$).
 %
 \begin{proposition}[Nonlinear dependence of  angular frequency ] \label{PropNdaf} $\mbox{ }$\\
 Let $b$ be a real number  and let $u_\varepsilon$  be the solution of:
 \begin{eqnarray*}
    \ddot{u}+ \omega_0^2 u +\varepsilon a(u-b)_+ = 0, &
   &  u_\varepsilon(0)= a_0 + \varepsilon a_1, \;  \dot{u}_\varepsilon(0)=0.
 \end{eqnarray*}
 If $|a_0| > |b|$
then there exists $\gamma > 0$, such that, for
 $ t \in [0,T_\varepsilon]$, with $T_\varepsilon = \displaystyle \frac{\gamma}{\varepsilon}$,
we have the following expansion in $C^2([0,T_\varepsilon],\mathbb R)$:
\begin{eqnarray*}
      u_\varepsilon(t) & = & v_0(\omega_\varepsilon t) +
                    \varepsilon v_1(\omega_\varepsilon t) +
                     \mathcal{O}(\varepsilon^2)
                     \qquad  with \\
     v_0(s) & =& a_0  \cos(s),  \\
     v_1(s)  &  = &  \displaystyle \sum_{k=1}^{+\infty} d_k \cos(k s), \\
\end{eqnarray*}
\begin{eqnarray*}
      d_0 &= &  -\frac{a |a_0|}{\pi \omega_0^2}
                \left ( \sin(\beta) - \frac{b}{|a_0|} \beta \right)
                     \quad
 \mbox{ where }
 \beta  =  \displaystyle \arccos\left(  \frac{b}{a_0} \right) \in [0,\pi], \\
       d_k &= &  \frac{- a |a_0|}{\pi \omega_0^2(1-k^2)}
                \left (\frac{ \sin((k+1)\beta)}{k+1}
                      +  \frac{ \sin((k-1)\beta)}{k-1}
               -   \frac{2b \sin(k\beta)}{|a_0| k}
                           \right),

\;  k \geq 2,\\
     d_1 &=&  a_1 -A
 \quad \mbox{ where }  A =\displaystyle  \sum_{k\neq 1} d_k, \\
   \omega_\varepsilon &= & \omega_0 + \varepsilon \omega_1 + \varepsilon^2 \omega_2, \\
    \omega_1 &= & \frac{a }{ 2\pi \omega_0}
           \left (  \frac{\sin (2\beta)}{2} + \beta
  - \frac{2b\sin (\beta)}{|a_0|} \right ),\\
 \omega_2 &= & \displaystyle \frac{ \omega_1 d_1}{a_0} - \frac{ a}{ \omega_0 \pi a_0}
  \int_0^{\pi} H(a_0 \cos(s)-b)v_1(s) \cos(s)ds.
\end{eqnarray*}
 \end{proposition}
Notice that if
 $|a_0| < |b|$, there is no interaction with the weak unilateral spring.
 Thus the linearized solution is the exact solution.
\\
\noindent{\rm\bf Proof : }
   There are two similar cases, $a_0$ positive or negative.
\\
 \underline{First case}: assume $a_0 > 0$.
 With  the previous notations, the method of strained coordinates yields
 the following equations:
\begin{eqnarray*}
  v_0''  + v_0 & =&  0, \quad v_0(0)=a_0,\, \dot{v}_0(0)=0
\quad \mbox{ so }  v_0(s) = a_0 \cos(s),
  \\
-\alpha_0  ( v_1''  + v_1) & = &  a(v_0 - b)_+ -  \alpha_1 v_0
 =  a a_0 (\cos(s) -b/a_0)_+  - \alpha_1 a_0 \cos(s),
 \\
 -\alpha_0  ( r_\varepsilon''  + r_\varepsilon) & = &
 aH(v_0 - b)v_1 -  \alpha_2 v_0 -\alpha_1 v_1
   + R_\varepsilon.
\end{eqnarray*}
After some computations
the Fourier expansion of $(\cos(s)-c)_+$ for  $|c| < 1$ is:
\begin{eqnarray}
(\cos(s)-c)_+ &= & \sum_{k=0}^{+\infty} c_k \cos(k s), \label{eq:seriecosmc+}\\
c_0 &= & \frac{\sin(\beta)  -  c \beta }{\pi }
      \qquad  \mbox{ where }
        \beta  =  \displaystyle \arccos\left( c \right) \in [0,\pi], \\
c_1 &= &  \frac{1}{\pi}
           \left ( \frac{\sin(2\beta)}{2} +\beta  -2 c \sin(\beta) \right),\\
c_k &=&  \frac{1}{\pi}
                \left (\frac{ \sin((k+1)\beta)}{k+1}
                      +  \frac{ \sin((k-1)\beta)}{k-1}
               -   \frac{2c \sin(k\beta)}{ k}
                           \right),
            \; k \geq 2.
\end{eqnarray}
The non secular  condition
 $\displaystyle \int_0^\pi (a(v_0 -b)_+ -\alpha_1 v_0) \cos(s) ds = 0$,
with $c=b/a_0$  yields $ \alpha_1 = a c_1$.
  Now, we can compute $\omega_1=\alpha_1/(2\omega_0)$
 and  $v_1$ with a cosines expansion:
$v_1(s)= \displaystyle \sum_{k} d_k \cos(ks)$
with $d_k=  \displaystyle -\frac{aa_0}{\alpha_0} \frac{c_k}{1-k^2}$  for $k \neq 1$.
 The coefficient $d_1$ is then obtained with the initial data $v_1(0)=a_1,$
  $\dot{v}_1(0) =0.$
\\
$\alpha_2$, is obtained with  the non secular condition for $r_\varepsilon$:
\\
$ 0 = \displaystyle
\frac{1}{\pi}\int_0^{\pi} (aH(v_0 -b)v_1 -\alpha_2 v_0 - \alpha_1 v_1) \cos(s) ds.
 $
This condition is rewritten as follow
\\
$ \alpha_2 = \displaystyle \frac{2 \omega_0 \omega_1 d_1}{a_0} - \frac{2 a}{\pi a_0}
  \int_0^{\pi} H(a_0 \cos(s)-b)v_1(s) \cos(s)ds$,
which gives  $ \omega_2$ since
 $ \omega_2 =  \frac{\alpha_2 - \omega_1^2}{2 \omega_0}$.
\\
\underline{ Second case}:
 when $a_0 = -|a_0| < 0 $,  by a similar way, we obtain a similar expansion,
 except that $(v_0(s) -b)_+ = |a_0|(-\cos(s) - b/a_0)_+$.
 And  we  only need the Fourier expansion of
 $(-\cos(s)-c)_+ = \displaystyle \sum_k \tilde{c}_k \cos(ks)$,
\begin{eqnarray*}
\label{eq:c0c1ck}
\tilde{c}_0 &= & -\frac{\sin(\beta)  +  c \beta }{\pi }

     , \\
\tilde{c}_1 &= &  -\frac{1}{\pi}
           \left ( \frac{\sin(2\beta)}{2} +\beta  +2 c \sin(\beta) \right),\\
\tilde{c}_k &=& - \frac{1}{\pi}
                \left (\frac{ \sin((k+1)\beta)}{k+1}
                      +  \frac{ \sin((k-1)\beta)}{k-1}
               +  \frac{2c \sin(k\beta)}{ k}
                           \right),
            \; k \geq 2.
\end{eqnarray*}
 \epro

  When $|a_0|=|b|$,
 we have  another  asymptotic expansion only valid for time of order
$\frac{1}{\sqrt{\varepsilon}}$ when the unilateral spring interacts with the mass.
 \begin{proposition}[The  critical case ]\label{Propcritic} $\mbox{ }$\\
 Let $b$ be a real number, $ b \neq 0$,
  and  consider, the solution $u_\varepsilon$ of:
 \begin{eqnarray*}
    \ddot{u}+ \omega_0^2u +\varepsilon a(u-b)_+ = 0, &
   &  u_\varepsilon(0)= a_0 + \varepsilon a_1, \;  \dot{u}_\varepsilon(0)=0.
 \end{eqnarray*}
 If $|a_0|=|b|$ then  we have
\begin{eqnarray*}
  u_\varepsilon(t) & = & (a_0+\varepsilon a_1)\cos(\omega_0 t) + \mathcal{O}(\varepsilon^2),
 \end{eqnarray*}
 for  $t \leq \displaystyle T_\varepsilon
                       = \left \{
                     \begin{array}{cccc}
 \displaystyle \frac{\gamma}{\sqrt{\varepsilon}} & if & |a_0+\varepsilon a_1| > |b| &
                                            \mbox{where  $\gamma > 0,$} \\
                    + \infty              & \mbox{ else.} & &
                      \end{array}
                      \right.
\end{proposition}
The method of strained coordinates gives us
 the {\it linear }  approximation for $u_\varepsilon(t)$, with $s=t$, i.e.
 $\omega_\varepsilon =1$.
If $|u_\varepsilon(0)| < |b|$, the exact solution is the solution of the linear
 problem $\ddot{u}+\omega_0^2u=0$.

\\
Otherwise,  if $|u_\varepsilon(0)| > |b|$,
since, $|b|$ is the maximum of $v_0(s)=a_0\cos(s)$, a new phenomenon appears,
 during each period, $|u_\varepsilon(t)| > |b|$ on  interval of time
 of order $\sqrt{\varepsilon}$ instead of $\varepsilon$.
Then   $T_\varepsilon$ is smaller  than in Proposition \ref{propexplicit}.
\\
  To explain this  phenomenon, we give precise estimates of the remainder
 when we expand $(v_0+\varepsilon v_1 + \varepsilon^2 r_\varepsilon)_+$ in the next section,
   see  Lemmas \ref{LDLL} and \ref{LHzi} below.

    \section{Expansion of $(u+\varepsilon v)_+$ \label{sexp+}}

 We give some useful lemmas to perform asymptotic expansions and
to estimate precisely the remainder  for the   piecewise linear map
    $ u \rightarrow u_+= \max(0,u).$

 \begin{lemma}

 {\bf [Asymptotic expansion  for $(u+\varepsilon v)_+$

  ] }
 $\mbox{ }$\\
\label{LDLL}
Let be $T>0$, $M>0$, $u, v$ two real valued functions defined on $I=[0,T]$,
 $$J_\varepsilon =\{t \in I, |u(t)|\leq \varepsilon M \},$$
$\mu_\varepsilon(T)$ the measure of the set $J_\varepsilon$ and
  $H$ be the Heaviside  step function,
then
\begin{eqnarray*}
     (u+\varepsilon v)_+  = (u)_+ + \varepsilon H(u)v
    + \varepsilon \chi_\varepsilon(u,v) ,
 & \mbox{ with } &
  H(u) = \displaystyle \left \{ \begin{array}{cc}
            1   & \mbox{ if }  u > 0 \\
            0   &  \mbox{ elsewhere }
         \end{array}\right. ,
\end{eqnarray*}
 and  $\chi_\varepsilon(u,v)$ is a non negative piecewise linear
    function and  1-Lipschitz with respect to $v$,
   which satisfies for all $\varepsilon$,
\\
  if  $ |v(t)|\leq M $  for any $t \in I$:
  \begin{eqnarray} \label{eq:DLL}
        |\chi_\varepsilon(u,v)| \leq   |v| \leq M, &
           &
        %
       \displaystyle \int_0^T \left |\chi_\varepsilon(u(t),v(t))  \right|dt   \leq
                     M \mu_\varepsilon(T).
  \end{eqnarray}
\end{lemma}
The point in inequality (\ref{eq:DLL}) is  the remainder $\varepsilon \chi_\varepsilon$ is only
 of order $\varepsilon$ in $L^\infty$ but of order $\varepsilon \mu_\varepsilon$ in $L^1$.
 In general, $\mu_\varepsilon$ is not better than a constant,
take for instance $u\equiv 0$. Fortunately,
 it is proved below that  $\mu_\varepsilon $ is often of order
 $\varepsilon$, and for some critical cases of order $\sqrt{\varepsilon}$.
\\
\noindent{\rm\bf Proof : } Equality (\ref{eq:DLL})
defines $\chi_\varepsilon$ and can be rewritten as follow:
  \begin{eqnarray}  \label{eqchi}
       \chi_\varepsilon(u,v)

 &=&
     \displaystyle  \frac{  (u+\varepsilon v)_+  -    u_+ - \varepsilon H(u)v}{\varepsilon}.
\end{eqnarray}

  So,  $\chi_\varepsilon$ is non negative since
  $u \rightarrow u_+$ is a convex function.
We also  easily see that the map
$(u,v) \rightarrow \chi_\varepsilon(u,v) $ is piecewise linear, continuous
 except on the line $u=0$ where $\chi_\varepsilon $ has a jump  $-v$.
 This jump comes from the Heaviside step function.
An explicit computations gives us the simple and useful formula:
   \begin{eqnarray}  \label{formulachi}
      0 \leq \varepsilon  \chi_\varepsilon(u,v)

    &=&
 \displaystyle \left \{ \begin{array}{cc}
    \displaystyle  |u+\varepsilon v|   & \mbox{ if } |u+\varepsilon v|  <   |\varepsilon v| \\
            0   &  \mbox{ elsewhere }
         \end{array}\right.
   .
 \end{eqnarray}
  We then have immediately $ 0 \leq \chi_\varepsilon(u,v)  \leq |v| $.
  Let $u $ be  fixed, then $v \rightarrow \chi_\varepsilon(u,v)$
 is one Lipschitz with respect
  to $v$.
Furthermore,
the support of $\chi_\varepsilon $ is included in $J_\varepsilon$,
    which  concludes the proof.
 \epro
\noindent
We now  investigate the size of  $\mu_\varepsilon(T)$, see \cite{BJ,J4}
for similar results about $\mu_\varepsilon(T)$ and other applications.
\\
With notations
 from  Lemma \ref{LDLL} we have.
\begin{lemma}[Order of $\mu_\varepsilon(T)$]\label{LHzi}

 Let $u$ be  a smooth  periodic function, $M$ be a positive constant   and
  $\mu_\varepsilon(T)$ the measure of the set
 $ \{t \in I, |u(t)|\leq \varepsilon M \}$.
\\
    If $u$   has only  simple roots on $I=[0,T]$

then  for some  positive $C$,
     \begin{eqnarray*}
 \mu_\varepsilon(T)  &\leq& C\varepsilon \times T .
 \end{eqnarray*}
More generally, if $u$
    has also double   roots
then  \begin{eqnarray*}  \mu_\varepsilon(T) &\leq &C\sqrt{\varepsilon} \times T.\end{eqnarray*}
\end{lemma}
The measure of such set $J_\varepsilon$ implies
many applications in averaging lemmas,
for a characterization of $\mu_\varepsilon$ in a multidimensional framework
 see \cite{BJ,J4}.
\\
Notice that
   any non zero  solution $u(t)$
  of any linear homogeneous second order ordinary differential equation
  has  always simple  zeros, thus
 for any constant $c$ the map $t \rightarrow u(t)-c$ has at most double roots.
\medskip \\
\noindent{\rm\bf Proof : }
 First assume  $u$ only has simple roots on a  period $[0,P]$,
 and let $Z = \{t_0 \in [0,P],u(t_0)=~0\}$.

 The set $Z$ is  discrete
  since $u$ has only  simple roots which implies that
     roots of $u$ are isolated.
 Thus $Z$ is a finite subset of $[0,P]$:
   $Z=\{t_1, t_2, \cdots,t_N \}$.
We can choose an open  neighborhood $V_j$ of each  $t_j$
 such that $u$ is a  diffeomorphism on $V_j$
  with derivative $|\dot{u}| > |\dot{u}(t_j)|/2$.
On  the compact set $K=[0,P]-\cup V_j$,
  $u$ never vanishes, then
   $\displaystyle \min_{t \in K} |u(t)|=\varepsilon_0 >0.$
 Thus, we have for all $\varepsilon M < \varepsilon_0$, the length of $J_\varepsilon$ in $V_j$
 is
   $\displaystyle  |V_j \cap J_\varepsilon|  \leq \frac{4 \varepsilon M}{|\dot{u}(t_j)|} $.
  As $\mu_\varepsilon$ is additive ($\mu_\varepsilon(P+t)= \mu_\varepsilon(P)+\mu_\varepsilon(t)$),
  its growth is linear.
 Thus, for the case with simple roots, we get
$\mu_\varepsilon(T) = \mathcal{O}(\varepsilon T)$.
\\
For the general case with double roots,
on each small neighborhood of $t_j$: $V_j$,
we have with a Taylor expansion,
  $|u(t_j+s)| \geq  d_j|s|^l$, with $ 1 \leq l \leq 2$, $d_j > 0$, so,
 $\displaystyle  |V_j \cap J_\varepsilon| \leq 2 (\varepsilon M/d_j)^{1/l}  $, then
$\displaystyle \mu_\varepsilon(P) =  \mathcal{O}(\sqrt{\varepsilon}),
 which  is enough
to  conclude the proof.
\epro

    \section{Several  degrees  of freedom \label{sNdof}}

 Now, we investigate the case with $N$ mass.
 We use, the method of strained coordinates in three cases.
 We present  the formal computations for each expansion. The mathematical proofs
 are postponed in the Appendix.
\\
 In subsection \ref{ssNdof1}, the initial condition is near an eigenvector
 such that
 the approximate solution  stays periodic.
  We gives such  initial condition near an eigenvector
 in subsection \ref{sec:annm} to get an approximate nonlinear normal mode
 up to the order $\varepsilon^2$.
 Finally,  in subsection \ref{ssNdof2}, all modes are excited.
 An extension of the method of strained coordinates is still possible but
 only at the first order with less accuracy.
\\
The system studied is the following:

$M\ddot{U} + K U  +\varepsilon (AU -B)_+=0$,
where, for each component,

   $[(AU-B)_+]_k = \displaystyle\left( \sum_{j=1}^N a_{kj}u_j -b_k\right)_+$,
$M$ is  a  diagonal $N \times N$
mass matrix  with positive terms on the diagonal,
$K$ is  the stiffness  matrix which   is  symmetric definite positive.
 It is also  possible to add many
  terms $\varepsilon (AU -B)_+$
 modeling small defects.
   For a such system, endowed with a natural convex
 energy for the linearized part,
  we can control the $\varepsilon$-Lipschitz last term for $\varepsilon$ small enough up to large time,
 so for $\varepsilon <<1$ the solutions remain bounded for  time of the order $\varepsilon^{-1} $.

 \\
To simplify the exposition,
 we deal with the following   diagonalized system for the linear part,
 keeping the same notations, except for the positive diagonal matrix $\Lambda$:
\begin{eqnarray}\label{sysMK+}
   \ddot{U} + \Lambda^2 U  &=& -\varepsilon (AU -B)_+ .
 \end{eqnarray}
In fact, with a linear change of variables, we have in the right hand side of
the previous equation \eqref{sysMK+} a linear combinations of terms
 $\varepsilon (AU -B)_+$.

\subsection{Initial condition near  an eigenvector,
\\  $\mbox{ }$ \hspace{15mm}second order approximation
     \label{ssNdof1}}

  For the system (\ref{sysMK+}), we take
 an initial condition near  an eigenmode of the linearized system
denoted for instance by  index $1$  .
\begin{eqnarray*}
  \left \{
  \begin{array}{ccccc} u_1^\varepsilon(0) & =& a_0+\varepsilon a_1, &\dot{u}_1^\varepsilon(0)=0, & \\
   u_k^\varepsilon(0)&=&0\;\; + \varepsilon a_k,  &\dot{u}_k^\varepsilon(0)=0,&   for \; k \neq 1.
 \end{array}
   \right.
\end{eqnarray*}
We impose $a_2,\cdots, a_N$ later  to have a periodic approximation, but
 $a_1$ is a free constant as $a_0$.
It is a key point to apply the method of strained coordinates.
\\
 We use  the same time $ s= \omega_\varepsilon t$  for each component
  and  the following  notations.
\begin{eqnarray*}
\begin{array}{cclccl}
   \omega_\varepsilon   & =&  \omega_0 + \varepsilon \omega_1+\varepsilon^2 \omega_2,
  &
   \omega_0 & = & \lambda_1,
  \\
    (\omega_\varepsilon)^2 & = & \alpha_0 + \varepsilon \alpha_1+\varepsilon^2 \alpha_2
 + \mbox{\rm O}(\varepsilon^3),
  &
       \alpha_0 &=& \omega_0^2= \lambda_1^2,
  \\
          \alpha_1&=& 2\omega_0 \omega_1,
    &  \alpha_2&=& \omega_1^2 + 2 \omega_0 \omega_2,
   \\
   u_j^\varepsilon ( t)   & = &   v_j^\varepsilon (s)
   =    v_j^0(s) +\varepsilon v_j^1(s)+ \varepsilon^2  r_j^\varepsilon(s),
 & j&=&1,\cdots,N .
 \end{array}
 \end{eqnarray*}
Replacing, this ansatz in the System (\ref{sysMK+}) we have in variable $s$,
\begin{eqnarray*} \label{dasNLNM}
  ( \omega_\varepsilon)^2 (v_k^\varepsilon)''  + \lambda_k^2 v_k^\varepsilon
& =&
  -\varepsilon \displaystyle\left( \sum_{j=1}^N a_{kj}
v_j^\varepsilon\left(s \right) -b_k\right)_+,
 \end{eqnarray*}
 and then performing the expansion for all $k \in\{1,\cdots,N\}$,
\begin{eqnarray*}  \label{eqNLM0}
  L_k v_k^0  &=& \alpha_0 (v_k^0)''  + \lambda_k^2 v_k^0
 = 0 ,
  \\
\label{eqNLM1}
-L_k v_k^1
 &= &  \displaystyle\left( \sum_{j=1}^N a_{kj}
v_j^0 -b_k\right)_+ + \alpha_1(v_k^0)'', \\
\label{eqNLM2}
 -L_k r_k^\varepsilon
 & =&
     \displaystyle H\left( \sum_{j=1}^N a_{kj}
v_j^0 -b_k\right)
\left( \sum_{j=1}^N a_{kj}
v_j^1\right)
+ \alpha_2 (v_k^0)''+ \alpha_1 (v_k^1)''
+R_k^\varepsilon.
\end{eqnarray*}
  First we have $v_1^0(s)= a_0 \cos(s)$.
\\
  Equations for $v_k^0$, for all $k\neq 1$, with zero initial data
 give us $ v_k^0 = 0$.
\\
 In equation for $v_1^1$,
we remove the secular term in the  right hand side,
\begin{eqnarray*}
-\alpha_0( (v_1^1)''+v_1^1)
 =
  \displaystyle\left(  a_{11}
v_1^0 -b_1\right)_+ + \alpha_1(v_1^0)''=r.h.s.
 &
 \quad
 & v_1^1(0)=a_1, \; (v_1^1)'(0)=0 .
\end{eqnarray*}
The orthogonality of the $r.h.s$ with $\cos(s)$ leads to define $\alpha_1$.
For instance, if $b_1  =0$,
 we have as in Proposition \ref{PropNdaf},
 $\alpha_1= \displaystyle \frac{a_{11}}{2}$
 and $ \omega_1=\displaystyle \frac{a_{11}}{4\lambda_1}$.
\\
Now, $\alpha_1$ is fixed,
so $v_1^1$ is a well defined  even $2\pi$ periodic function.
\\
 Then, for $k \neq 1$, we can compute  $v_k^1$  with the
 simplified equation,
\begin{eqnarray} \label{eqN1} \displaystyle
 \alpha_0 (v_k^1)''  + \lambda_k^2 v_k^1
   &=&   - \displaystyle\left(  a_{k1}
v_1^0 -b_k\right)_+.
\end{eqnarray}
Let $\phi_k^1$ be the unique $2\pi$  periodic  solution of equation (\ref{eqN1}).
Such function exists and is unique  if  $\lambda_k \notin \lambda_1 \mathbb Z$.
Furthermore $\phi_k^1$ is an even function as the right hand side of  equation
(\ref{eqN1}). Let $v_k^1$ be $\phi_k^1$, i.e.
               $v_k^1(0)=a_k=\phi_k^1(0)$ and $(v_k^1)'(0)=0$ for all $k\neq 1$. \\
   The term $r_1^\varepsilon$, with null initial data,
 has  a simplified equation
 since $v_k^0 \equiv 0$ for all $k\neq 1$,
 \begin{eqnarray*}
 -L_1 r_1^\varepsilon
 & =&
     \displaystyle H\left(  a_{11}
v_1^0 -b_1\right)
\left(   \sum_{j=1}^N a_{1j}
v_j^1\right)
+ \alpha_2 (v_1^0)''+ \alpha_1 (v_1^1)''
+R_1^\varepsilon.
 \end{eqnarray*}
 We now  can compute numerically
 $\alpha_2$ to avoid secular term in the right hand side, $R_1^\varepsilon$ excepted,
 with the following condition,
\begin{eqnarray*}
 0 &= &
   \int_0^\pi \left[ \displaystyle H\left(  a_{11}
v_1^0 -b_1\right)
\left(
\sum_{k=1}^N a_{1k}
v_k^1\right)
+ \alpha_2 (v_1^0)''+ \alpha_1 (v_1^1)''
 \right] \cdot \cos(s) ds.
 \end{eqnarray*}

Rewriting this condition, we obtain an equation for $\alpha_2$
in Theorem \ref{Th41} below, which gives
$\omega_2$ as in the Proposition \ref{PropNdaf}.
 \\
  For each $k\neq 1$, $\lambda_k$ and $\lambda_1$ are $\mathbb Z$ independent.
  So  $r_k^\varepsilon$ stays bounded
  for large time since there is no resonance of the order one at the first
 order in equations, $-L_k r_k^\varepsilon = (\cdots) + R_k^\varepsilon$.
 This  is the  technical  part of the proof
   to validate rigorously
  and  to find the time of validity of
  such  asymptotic expansion.
The complete proof to bound $(r_1^\varepsilon, \cdots, r_N^\varepsilon)$ for large time
 is to be found in the Appendix.
\\
 We now state our result with previous notations.

\begin{theorem}\label{Th41}
The following expansion is valid
 on $(0,T_\varepsilon)$, with $T_\varepsilon \rightarrow  + \infty$ when $\varepsilon \rightarrow 0$
 under assumption $ \lambda_k $ and  $ \lambda_1$ are $ \mathbb Z$ independent
 for each $k\neq 1$,
 \begin{eqnarray}
\label{eq:uup=0}
  \left\{
   \begin{array}{cccccccc}
   u_1^\varepsilon(t) &  = &  v_1^0(\omega_\varepsilon t) & +& \varepsilon v_1^1(\omega_\varepsilon t)
                   &+& \mathcal{O}(\varepsilon^{2}),&
                \\
 u_k^\varepsilon(t) &  = & 0 & +& \varepsilon v_k^1(\omega_\varepsilon t)
                   &+& \mathcal{O}(\varepsilon^{2}),&   k \neq 1,
   \end{array}
    \right.
 \end{eqnarray}
where $v_1^0(s)$, $\alpha_1$, $\omega_1$, $v_1^1(s)$,
  $v_k^1(s)$, $\phi_k^1$, $\alpha_2$, $\omega_2$
 are defined successively as follows, with differential operators
 $L_k = \displaystyle \lambda_1^2 \frac{d^2}{ds^2} + \lambda_k^2 $,
\begin{eqnarray}
  v_1^0(s) & = & a_0 \cos(s), \\
  \alpha_1 &= & \displaystyle \frac{2}{a_0\pi}\int_0^\pi
                                 (a_{11}v_1^0(s) - b_1)_+ \cos (s) ds, \\
  \omega_1 &= & \frac{\alpha_1}{2 \omega_0}, \\
-L_1 v_1^1   & = &  \displaystyle\left(  a_{11}
v_1^0 -b_1\right)_+ + \alpha_1(v_1^0)'',
\quad  v_1^1(0)=a_1, \; (v_1^1)'(0)=0,
\\
  \label{eqN1th}
-L_k v_k^1   & = &  \displaystyle\left(  a_{k1}
v_1^0 -b_k\right)_+ ,
\quad  v_k^1(0)=a_k:=\phi_k^1(0), \; (v_k^1)'(0)=0, \mbox{ for } k \neq 1, \label{eq:Lkvk1p=0}
\\
 \phi_k^1 & \mbox{is} & \mbox{the unique $2\pi$ periodic solution
  of  (\ref{eqN1th}) for all $k\neq 1$},
\\
  \label{ortth42}
 \alpha_2 &= & \displaystyle \left( \frac{2}{a_0\pi}\int_0^\pi
   \left\{ H(a_{11}v_1^0(s) - b_1)
     \left(a_{11}v_1^1+ \sum_{k\neq 1}^N a_{1k}\phi_k^1 \right)
  \right\}\cos (s) ds \right) - \alpha_1,
\\
  \omega_2 &= & \frac{\alpha_2 - \omega_1^2}{2 \omega_0}.
\end{eqnarray}
 Furthermore, if $(a_{j1} v_1^0-b_j)$ has got only simple roots
for all $j=1,\cdots,N$,
\\ then   $T_\varepsilon$ is of the order $\varepsilon^{-1}$,
\\ else  $T_\varepsilon$ is of the order $\varepsilon^{-1/2}$.
\end{theorem}
In the theorem $\phi_k^1$ is classically obtained by  a Fourier series.
We give some indication of its initial condition in the
 next subsection \ref{sec:annm}.

\subsection{Approximate non linear normal mode}
\label{sec:annm}

The special initial conditions  of the previous subsection
  can   be explicited
 in order to find a solution where all the components are in phase at the
same frequency.
Indeed we shall obtain an approximate curve of initial conditions for  which the solution is periodic up to the order $\varepsilon$ for a time of the order $\displaystyle \frac{1}{\varepsilon}$: this is up to the approximation a non linear normal mode in the sense of Rosenberg \cite{Rosenberg}; see \cite{JPS} for a computation of non linear normal mode with unilateral contact and \cite{nnm-kpgv} for a synthesis on non linear normal modes.

\begin{theorem}
\label{thm:app-nnm}
  For the system (\ref{sysMK+}), we take
 an initial condition close to an eigenmode of the linearized system
denoted for instance by  index $1$  .
\begin{eqnarray*}
  \left \{
  \begin{array}{ccccc} u_1^\varepsilon(0) & =& a_0+\varepsilon a_1, &\dot{u}_1^\varepsilon(0)=0, & \\
   u_k^\varepsilon(0)&=&0+\varepsilon a_k,  &\dot{u}_k^\varepsilon(0)=0,&   for \; k \neq 1.
 \end{array}
   \right.
\end{eqnarray*}
the results of the previous theorem are still valid for $T_\varepsilon = \displaystyle \frac{\gamma}{\varepsilon}$ but now we are looking for
$a_k=v_k^1(0)=$ in \eqref{eq:Lkvk1p=0}.
\\
Moreover for the cases detailed below, the solution is periodic with angular frequency $\omega_{\varepsilon} $  up to the order $\varepsilon$ and  for $0 \le t \le  T_{\varepsilon }  $
\begin{enumerate}
\item assume $b_k=0$, for $ k\ne 1$ with the initial condition defined with
  \begin{equation}
    \label{eq:akpourbk0}
    a_k= \frac{a_{k1}a_0}{2(\lambda_1^2-\lambda_k^2)} - \frac{|a_{k1}a_0|}{\lambda_k^2 \pi} +\frac{a_{k1}a_0}{2} \sum_{l=1}^{+\infty} \frac{(-1)^l}{(4 l^2\lambda_1^2 -\lambda_k^2)(4l^2 -1)}
  \end{equation}

\item assume $0 < \frac{b_k}{|a_{k1}a_0|}<1$, and $a_0a_{k1}<0$ for $ k\ne 1$ then for the initial condition defined with
  \begin{equation}
    \label{eq:akpourbk0bkak1a0pp1p}
    a_k=a_0a_{k1} \left [\frac{1}{\lambda_1^2-\lambda_k^2} - \frac{b_k}{\lambda_k^2 a_{k1}a_0} -\sum_{l=1}^{+\infty} \frac{c_l}{l^2 \lambda_1^2 -\lambda_k^2}
\right]
\end{equation}
where $c_l$ are defined in \eqref{eq:seriecosmc+} with $c=-\frac{b_k}{a_{k1}a_0}$
\item assume $0 < \frac{|b_k|}{|a_{k1}a_0|}<1$, and $a_0a_{k1}>0$ for $k\ne 1 $ then for the initial condition defined with

  \begin{equation}
    \label{eq:akpourbk0bkak1a0pp1m}
    a_k=a_0a_{k1} \left [ \sum_{l=1}^{+\infty} \frac{c_l}{l^2 \lambda_1^2 -\lambda_k^2}
\right]
  \end{equation}
where $c_l$ are defined in \eqref{eq:seriecosmc+} with  $c=\frac{b_k}{a_{k1}a_0}$
\end{enumerate}

\end{theorem}
\begin{remark}
The other cases are less interesting but may solved similarly.
\end{remark}

\noindent{\rm\bf Proof : }
The principle of the proof is simple; $v_k^1$ is solution of the differential equation  \eqref{eq:Lkvk1p=0}
with  $v_k^1(0)=a_k$ and $a_k$ has  to be determined in order that the function $v_k^1$  has an angular frequency equal to one. It is elementary that the solution of \eqref{eq:Lkvk1p=0} is
\begin{equation}
  \label{eq:v_k^1=b...}
  v_k^1=A \cos\left(\frac{\lambda_k}{\lambda_1} s\right)
  +B \sin\left(\frac{\lambda_k}{\lambda_1} s\right) + \phi_k^1(s)
\end{equation}
where $\phi_k^1$ is a particular solution associated to the right hand side which is of angular frequency equal to $1$; note that $B=0$ as the initial velocity is null; we can get a function of angular frequency equal to $1$ by setting $a_k=\phi_k(0)$; this condition may be written explicitly with formulas \eqref{eq:seriecosmc+} which provides the expansion in Fourier series; formulas
\eqref{eq:akpourbk0}, \eqref{eq:akpourbk0bkak1a0pp1p},     \eqref{eq:akpourbk0bkak1a0pp1m}   are then derived easily.
\begin{enumerate}
\item for $b_k=0$ for $ k\ne 1$, \eqref{eq:Lkvk1p=0} is written:
  \begin{equation}
    -L_k v_k^1= a_{k1} a_0 \frac{cos(s)}{2} +  | a_{k1} a_0 |  \frac{cos(s)}{2}
  \end{equation}
we use formula \eqref{eqabscosF} to get the particular solution
\begin{equation}
  \phi_k^1=  a_{k1} a_0 \frac{cos(s)}{2(\lambda_1^2 -\lambda_k^2) }-
\frac{|a_{k1} a_0|}{\lambda_k^2 \pi} +
\frac{|a_{k1} a_0|}{2} \sum_{l=1}^{+\infty} \frac{(-1)^l}{4l^2 \lambda_1^2-\lambda_k^21} \frac{cos(2l s)}{4l^2 -1}
\end{equation}
from which \eqref{eq:akpourbk0} is deduced.
\item for $0 < \frac{b_k}{|a_{k1}a_0|}<1$, and $a_0a_{k1}<0$ for $ k\ne 1$
\eqref{eq:Lkvk1p=0} is written:
\begin{gather}
      -L_k v_k^1= -a_{k1} a_0 \left (-cos(s) + \frac{b_k}{a_{k1} a_0} \right )_+  \text{ or }\\
 -L_k v_k^1= -a_{k1} a_0  \left [ \left (-cos(s) + \frac{b_k}{a_{k1} a_0} \right ) + \left (cos(s)  - \frac{b_k}{a_{k1} a_0} \right )_+  \right]  \text{ we use \eqref{eq:c0c1ck} to obtain } \\
 \phi_k^1=  -a_{k1} a_0 \left [ ( \frac{-cos(s)}{2(\lambda_1^2 -\lambda_k^2) } - \frac{b_k}{\lambda_k^2 a_{k1} a_0} ) + \sum_{l=1}^{+\infty} \frac{c_l cos(ls)}{l^2 \lambda_1^2 -\lambda_k^2} \right ]
 \text{ where } c_l \text{ is defined with }  \\
\beta =arccos( \frac{b_k}{a_{k1} a_0} ) \text{ which defines } c_l \text{ in  \eqref{eq:c0c1ck}  with  } c=\frac{-b_k}{a_{k1}a_0}
\end{gather}
\item for $0 < \frac{b_k}{|a_{k1}a_0|}<1$, and $a_0a_{k1}>0$ for $ k\ne 1$
\eqref{eq:Lkvk1p=0} is written:
\begin{gather}
       -L_k v_k^1= a_{k1} a_0 \left (cos(s) - \frac{b_k}{a_{k1} a_0} \right )_+  \text{ from which }\\
\phi_k^1= a_0 a_{k1} \sum_{l=1}^{+\infty} \frac{c_l cos(ls)}{l^2 \lambda_1^2 -\lambda_k^2}
\text{ where } c_l \text{ is defined  with }  \\
\beta =arccos( \frac{b_k}{a_{k1} a_0} ) \text{ which defines } c_l
\text{ in  \eqref{eq:c0c1ck}  with  }c=\frac{-b_k}{a_{k1}a_0}
\end{gather}
\end{enumerate}

\epro

\subsection{Numerical results of NNM}

\subsubsection{Using numerically Lindstedt-Poincar\'e  expansions}
Here we use the previous results and  compute numerically a solution of system
\eqref{sysMK+} using the approximation \eqref{eq:uup=0}

$$u^{\varepsilon}(t)=v^0(\omega_{\varepsilon} t) +\varepsilon v^1(\omega_{\varepsilon} t) +O(\varepsilon^2)$$
with the initial conditions of theorem \ref{thm:app-nnm}. The first term $v^0$ is easy to obtain; for the second term $v^1$ an explicit  formula is in principle possible using Fourier series such as for one degree of freedom but it is cumbersome so we choose to compute $v^1$ by solving numerically
 \eqref{eq:Lkvk1p=0}
 with a step by step algorithm; we use as a black-box the routine ODE of SCILAB \cite{scilab} to solve equations of theorem \ref{Th41} after  computing by numerical integration $\alpha_1$. We show numerical results for a system of the type:
\begin{gather} \label{sysMKF}
  M \ddot{X} + K X + \varepsilon F(X)=0
\end{gather}
we still denote $\lambda_j^2$ the eigenvalues and $\phi_j$ the eigenvectors of the usual generalized eigenvalue problem
\begin{gather}
\label{eq:kphi-ommphi}
  K \phi_j -\lambda_j^2  M \phi_j=0 \quad \text{with}\\
^t \phi_k M \phi_j=\delta_{kj}
\end{gather}
We set:
\begin{gather}
  X=\sum_j u_j \phi_j =\underline{ \underline{\phi}} \underline{u}
\end{gather}
In this basis, the system may be written componentwise:
\begin{gather}
\ddot{u_k}  +\lambda_k^2  u_k + ^t \phi_k F(\underline{ \underline{\phi}} \underline{u})=0
\end{gather}
For a local non linearity in the system \eqref{sysMKF}, written in the basis of the eigenvectors, we do not  obtain a system precisely of the form \eqref{sysMK+};
so we illustrate it with the non linearity:
\begin{gather}
  F(X)=  M \phi_1(X1-\beta_1)   =M \phi_1(\sum_j u_j \phi_{j1}  -\beta_1)_+
\end{gather}
where we denote by $\phi_1$ {\it any} normalized eigenvector of \ref{eq:kphi-ommphi};
then the system \eqref{sysMKF} is written:
\begin{gather}
\ddot{u_k}  +\lambda_k^2 u_k + (\sum_k\delta_{k1} \phi_{j1} u_j -\delta_{k1} \beta_1)_+
\end{gather}
so that $A$ and $B$ of \eqref{sysMK+} are:
$$a_{kj}=\delta_{k1} \phi_{j1} \quad b_k=\delta_{k1} \beta_1$$

We find in figure \ref{fig:lp-ener0_03} a numerical example of the Linstedt-Poincar\'e approximation for 5 degrees of freedom with $\varepsilon=0.063$ and with an energy of $0.03002$. The left figure shows the 5 components of the solution with respect to time; the right figure, the solution in the configuration space: absissa component 1 and ordinate components 2 to 5; these lines are rectilinear like in the linear case but the non symmetry may be particularly  noticed on the smallest component which corresponds to the mode where the non linearity is active.
\begin{figure}[!ht]
  \centering
\begin{center}
  \includegraphics[ width=7cm, height=7cm,angle=-90]{ener03perio-lp}
  \includegraphics[ width=7cm, height=7cm, angle=-90]{ener03rosen-lp}
\caption{Lindstedt-Poincar\'e, energy=0.03, 5 dof; left: components with respect to time; right: in configuration space  \label{fig:lp-ener0_03} }
\end{center}
\end{figure}

\subsubsection{Using optimization routines}
We find  also in figure \ref{fig:opt-ener0_03} a numerical example with the same energy of $0.03002$; it is computed with a purely numerical method described below. We notice that the solution is quite similar in both cases.

\label{sec:opti-routi}
\begin{figure} \label{fig:opt-ener0_03}
  \centering
\begin{center}
  \includegraphics[ width=7cm, height=7cm, angle=-90]{ener03fft}
  \includegraphics[ width=7cm, height=7cm,angle=-90]{ener03rosen}
\caption{Continuation and Powell hybrid, energy=0.03, 5 dof; left:fourier transform ; right: in configuration space  \label{}}
\end{center}
\end{figure}

The numerical expansions of the previous subsection gives valid results for $\varepsilon$ small enough; in many practical cases such as  \cite{HR}, $\varepsilon$ may be quite large; in this case, it is natural to try to solve numerically the following equations with respect to the period $T$ and the initial condition $X(0)$.
\begin{gather}
\label{eq:X0=XT}
  X(0)=X(T), \quad \dot{X}(0)=\dot{X}(T) \quad
E(X)=e
\end{gather}
In other words, we look for a periodic solution of prescribed energy; this last condition is to ensure to obtain an isolated local solution: the previous expansions show that in general, the period of the solution depends on its amplitude prescribed here by its energy. To try to solve these equations with a black-box routine for nonlinear equations such as ``fsolve'' routine of SCILAB \cite{scilab}(an implementation of a modification of Powell hybrid method which goes back to \cite{powell64}) in general fails to converge.
Even in case of convergence, we should address the question of link of this solution with normal modes of the linearized system.
\par So we prescribe that $e=c\varepsilon$ and for $\varepsilon\rightarrow 0$, the solution is tangent to a linear eigenmode. In the case where all the eigenvalues of the linear system are simple, we define $N$ (the number of degrees of freedom) non linear normal modes for which, it is reasonable to conjecture that they correspond to isolated solutions of  \eqref{eq:X0=XT} at least for small $\varepsilon$ if we enforce for example $\dot{X}(0)=0$.
\paragraph{Algorithm}
This definition of the solution of \eqref{eq:X0=XT} tangent to a prescribed linear eigenmode provides a simple way of numerical approximation: using a continuation method coupled with a routine for solving a system of non linear equations.
Define:
\begin{gather}
\mathcal{F}(\varepsilon,X_0, T)=[X(T)-X_0, E(X)-c\varepsilon] \\
\text{where X is a numerical solution of the differential system}\\
  M \ddot{X} + K X + \varepsilon F(X)=0 \\
X(0)=X_0, \quad \dot{X}(0)=X_1
\end{gather}
\begin{verse}
choose a small initial value of $ \varepsilon$ and an increment $ \delta$\\
choose an eigenvector $\phi_j$ \\
$X_0(0)= A_{\varepsilon} \phi_j$, $X_1(0)=  B_{\varepsilon} \lambda_j \phi_j$\\
with $E(X_0(0),X_1(0))=c \varepsilon$ \\
\end{verse}
\begin{verse}
{\bf for}  iter=1:itermax
$$\varepsilon=\varepsilon + \delta$$
 with $ X_0(iter-1)$
as a first approximation, solve for $X_0(iter)$,
\begin{gather}
 \mathcal{F}(\varepsilon,X_0,T)=0
\end{gather}
{\bf if } $||\mathcal{F}(\varepsilon,X_0,T)||> tolerance$ {\bf then}
 $\varepsilon=\varepsilon - \delta, \quad\delta=\delta/2 $\\
$\quad$   {\bf endif} \\
{\bf endfor}
\end{verse}

This algorithm may be improved by using  not only the solution associated to the previous value of $\varepsilon$ to solve
\begin{gather}
 \mathcal{F}(\varepsilon,X_0,T)=0
\end{gather}
but also the derivative of the solution with respect to $X_0, T$.
\paragraph{Numerical results}
These results are obtained by solving the differential equation with a step by step numerical approximation of the routine $ode$ of Scilab without prescribing the algorithm. As we are looking for a periodic solution, this numerical approximation may be  certainly improved in precision and computing time by using an harmonic balance algorithm.
In figure \ref{fig:opt-ener0_12}, \ref{fig:opt-ener0_2}, the same example with 5 degrees of freedom and energy equal $0.123$ and $0.201$ are displayed.

 On the left of figure \ref{fig:opt-ener0_12} we find the decimal logarithm of the absolute value of the Fourier transform of the solution; the Fourier transform is computed with the fast fourier transform with the routine $fft$ of Scilab; we notice the frequency zero due to the non symmetry of the solution and multiples of the  basic frequency; no other frequency appears; on the right the five components are plotted with respect to time; we still notice the non symmetry.

 On the left  of figure \ref{fig:opt-ener0_2} we find the solution in the configuration space and
on the right the five components are plotted with respect to time; we still notice the non symmetry.

In figure \ref{fig:opt20d-ener0_12} we find results with 20 degrees of freedom, $\varepsilon=0.272$ and energy of $0.129$; the NNM is computed by starting with an  eigenvector associated to the largest eigenvalue . We see on the left in the configuration space that the components are in phase and on the right, the Fourier transform shows zero frequencies and multiple of the basic frequency.

\begin{figure}[!ht]
\begin{center}
\includegraphics[ width=7cm, height=6cm]{ener012fft}
\includegraphics[ width=7cm, height=7cm]{ener012perio}
\caption{energy=0.123, 5 dof; left:fourier transform ; right:  with respect to time \label{fig:opt-ener0_12}}
\end{center}
\end{figure}

\begin{figure}[!ht]
\centering
\includegraphics[angle=0, width=7cm, height=6cm]{ener02rosen}
\includegraphics[angle=0, width=7cm, height=6cm]{ener02perio}
\caption{ energy=0.2, 5 dof; left: configuration space; right:with respect to time \label{fig:opt-ener0_2}}
\end{figure}
\begin{figure}
  \centering
\begin{center}
  \includegraphics[ width=7cm, height=7cm]{ener12rosen-20d}
  \includegraphics[ width=7cm, height=7cm]{ener12fft-20d}
\caption{energy=0.129, 20 dof; left:in configuration space; right: fourier transform \label{fig:opt20d-ener0_12}}
\end{center}
\end{figure}
In figure \ref{fig:opt20d-ener0_29} the energy is $0.29$ and the NNM is computed by starting with an  eigenvector associated to the smallest eigenvalue; we notice on the left, the solution in the configuration space: at zero each dof has a discontinuity in slope which is clear.

In figure \ref{fig:opt20dmode-ener0_28_29}, the shape of the NNM is displayed on the left for the NNM  starting from the eigenvector associated to the smallest eigenvalue and on the right for the NNM starting from the second smallest eigenvalue. We notice that the shape is quite similar to the shape of the linear  mode.

In figure \ref{fig:opt20d-ener0_28}, the NNM is computed  by starting at an  eigenvector associated to the second  smalest eigenvalue.
We notice that the first NNM is more asymmetric than the second one; on the left, the discontinuity of slope is smaller than for the first NNM.

\begin{figure}[!ht]
\centering
\includegraphics[angle=0, width=7cm, height=6cm]{ener0_29rosenlam1-20d}
\includegraphics[angle=0, width=7cm, height=6cm]{ener0_29fftlam1-20d}
\caption{ energy=0.29, 5 dof; left: configuration space; right:fft \label{fig:opt20d-ener0_29}}
\end{figure}
\begin{figure}
  \centering
\begin{center}
  \includegraphics[ width=7cm, height=7cm]{ener0_29mode-20d}
  \includegraphics[ width=7cm, height=7cm]{ener0_28modelam2_20d}
\caption{ 20 dof; left:energy 0.29 mode 1; right: energy 0.28, mode 2 \label{fig:opt20dmode-ener0_28_29}}
\end{center}
\end{figure}

\begin{figure}[!ht]
\centering
\includegraphics[angle=0, width=7cm, height=6cm]{ener0_28rosenlam2_20d}
\includegraphics[angle=0, width=7cm, height=6cm]{ener0_28fftlam2_20d}
\caption{ energy=0.28, 5 dof; left: configuration space; right:fft \label{fig:opt20d-ener0_28}}
\end{figure}

 \subsection{First order asymptotic expansion \label{ssNdof2}}

 In this subsection, we do not particularize the initial data on one eigenmode.
 We adapt  the method  of strained  coordinates since
  all modes are excited.
 We loose one order of accuracy compared to previous results  since
 each mode  does not stay periodic and becomes almost-periodic.
\\
 More precisely,
 the method of strained coordinates is used for each normal component,
  with the following  initial data
 \begin{eqnarray*}  u_k^\varepsilon(0)=a_k,  &\dot{u}_k^\varepsilon(0)=0, & \qquad  k=1,\cdots,N.
  \end{eqnarray*}
 Let us define  $N $ new times $s_k=\lambda_k^\varepsilon t$ and the   following ansatz,
\begin{eqnarray*}
  \begin{array}{ccccccl}
   \lambda_k^\varepsilon     &=&  \lambda_k^0 + \varepsilon \lambda_k^1,
   &&
  \lambda_k^0&=&\lambda_k,
   \\
   u_k^\varepsilon ( t)   & =&   v_k^\varepsilon (\lambda_k^\varepsilon t) &= & v_k^\varepsilon (s_k)

 & =&   v_k^0(s_k) +\varepsilon  r_k^\varepsilon(s_k) .
 \end{array}
 \end{eqnarray*}
The function $v_k^0$ are easily obtained by the linearized equation.
 Indeed, the only measured nonlinear effect for large time is given
 by  $(\lambda_k^1)_{k=1}^N$.
 To obtain these $N$ unknowns, we
replace  the previous ansatz in the system (\ref{sysMK+}),
 \begin{eqnarray*}
\begin{array}{ccc}
  ( \lambda_k^\varepsilon)^2 (v_k^\varepsilon)''(s_k)  + \lambda_k^2 v_k(s_k)
  & = & -\varepsilon \displaystyle\left( \sum_{j=1}^N a_{kj}
v_j^\varepsilon\left(\frac{\lambda_j^\varepsilon}{\lambda_k^\varepsilon}s_k \right) -b_k\right)_+.
\end{array}
 \end{eqnarray*}
The right hand side is written in variable $s_k$ instead of
$\displaystyle s_j$.
Performing the expansion with respect to  epsilon powers yields
\begin{eqnarray}
L_k v_k^0&=& ( \lambda_k^0)^2 (v_k^0)''(s_k)  + \lambda_k^2 v_k^0(s_k)
 = 0 , \\
-L_k  r_k^\varepsilon(s_k)
 &=&    \displaystyle\left( \sum_{j=1}^N a_{kj}
v_j^0\left(\frac{\lambda_j^0}{\lambda_k^0}s_k \right) -b_k\right)_+
+ 2\lambda_k \lambda_k^1 (v_k^0)''
+R_k^\varepsilon.
\label{eqNdofRemainderH}
 \end{eqnarray}
Noting that replacing $ v_j^\varepsilon(s_j)$ by
$v_j^0 \left(\frac{\lambda_j^0}{\lambda_k^0}s_k \right)$
in (\ref{eqNdofRemainderH})
implies a secular term of the order $\varepsilon t$,
since
 $s_j = \displaystyle\frac{\lambda_j^0}{\lambda_k^0}s_k  +\mathcal{O}(\varepsilon t)  $,
 the functions $v_j^0$ are smooth
and the map $S \rightarrow S_+$ is one-Lipschitz.
 These new kind of  errors $\mathcal{O}(\varepsilon t)$ are contained
in   the remainder of each right hand side:
\begin{eqnarray}
  \label{RFose}
    R_k^\varepsilon(t)  & = & \mathcal{O}(\varepsilon t)+ \mathcal{O}(\varepsilon |r^\varepsilon|),
 \qquad   |r^\varepsilon| = \displaystyle \sqrt{ \sum_{k=1}^N(r^\varepsilon_k)^2}.
\end{eqnarray}
If  $b_k=0$,
we identify the secular term
with the Lemma \ref{LSp} below and the relation   $S_+=\displaystyle S/2+|S|/2 $.
Then, we
remove  the resonant term
in the source term for the remainder $r_k^\varepsilon$, which gives us
 $ \lambda_k^1=\displaystyle \frac{a_{kk}}{4\lambda_k}$.
\\
If $b_k\neq 0$, we compute $\lambda_k^1$ numerically
 with the following orthogonality condition to $\cos(s)$  written in the
 framework of almost periodic functions,
\begin{eqnarray*}
0 & = &
  \lim_{T \rightarrow \infty} \frac{1}{T}\int_0^T
\left[
 \displaystyle\left( \sum_{j=1}^N a_{kj}
v_j^0\left(\frac{\lambda_j^0}{\lambda_k^0}s_k \right) -b_k\right)_+
+ 2\lambda_k \lambda_k^1 (v_k^0)''
\right]\cdot \cos(s)ds .
\end{eqnarray*}
The accuracy of the asymptotic expansion depends on the behavior of the
 solution $\phi=(\phi_1,\cdots,\phi_N)$ of the $N$ following decoupled
linear equations with  right coefficients $\lambda_k^1$ to avoid resonance
\begin{eqnarray} \label{linsysacc}
-L_k  \phi_k(s_k)
 &=&    \displaystyle\left( \sum_{j=1}^N a_{kj}
v_j^0\left(\frac{\lambda_j^0}{\lambda_k^0}s_k \right) -b_k\right)_+
+ 2\lambda_k \lambda_k^1 (v_k^0)''.
\end{eqnarray}
Furthermore  each function $r_k^\varepsilon$ depends  on all times $s_j$, $j=1,\cdots,N$
and becomes almost-periodic,
i.e.  $ r_k^\varepsilon =r_k^\varepsilon( s_1, \cdots,s_N)$.
 Thus the method of strained
 coordinates,
 only working for periodic functions,
fails to be continued.
\\
 Nevertheless, we obtain the following result proved in the Appendix.
\begin{theorem}[All modes]\label{Th43} $\mbox{ }$ \\
 If $\lambda_1, \cdots , \lambda_N$ are $\mathbb Z$ independent,
 then,  for any $T_\varepsilon= o(\varepsilon^{-1})$, i.e.
   such that
 \begin{eqnarray*} \displaystyle \lim_{\varepsilon \rightarrow 0} T_\varepsilon = +\infty,
  & \mbox{ and } &
   \displaystyle \lim_{\varepsilon \rightarrow 0} \varepsilon \times T_\varepsilon = 0,
 \end{eqnarray*}
 we have  for all $k=1,\cdots,N$,

\begin{eqnarray*}\label{aao1}
        \displaystyle \lim_{\varepsilon \rightarrow 0} \|  u_k^\varepsilon(t)-
                               v_k^0\left(  \lambda_k^\varepsilon \;  t\right)
                               \|_{W^{2,\infty}(0,T_\varepsilon)} = 0
 \end{eqnarray*}
where
 $ \lambda_k^\varepsilon=  \lambda_k + \varepsilon \lambda_k^1, $
 $v_k^0(s)= a_k \cos(s)$,  and $\lambda_k^1$ is defined by:
 \begin{eqnarray*}  \displaystyle    \lambda_k^1 &=&
\frac{1}{2 \lambda_k a_0}
\lim_{T \rightarrow + \infty}\frac{1}{T} \int_{0}^{ T}
  \left( \sum_{j=1}^N a_{kj}
v_j^0\left(\frac{\lambda_j^0}{\lambda_k^0}s_k \right) -b_k\right)_+
 \cos(s) ds.
\end{eqnarray*}
Furthermore, if  $b_k=0$, the previous integral yields:
$\displaystyle \lambda_k^1=\frac{a_{kk}}{4\lambda_k}$.
\end{theorem}
Notice that accuracy  and large time  are weaker  than
 these obtained in Theorem \ref{Th41}.
It is due to  the inevitable accumulation of the spectrum near the resonance
and the various times using in the expansion.
On the other side we have the following direct
 improvement from  the Theorem \ref{Th41}:
\begin{remark}[Polarisation] \label{polarisation}
 If only  one mode are excited, for instance the number 1,
 i.e.  $ a_1 \neq 0$,  $a_k=0$ for all $k\neq1$, then we have the estimate
 for all  $ t \in [0,\varepsilon^{-1}]$:
\begin{eqnarray*}
 \displaystyle   u_1^\varepsilon(t)=&    v_1^0\left(  \lambda_1^\varepsilon \;  t\right)
                                &+  \mathcal{O}(\varepsilon).
                               \\
         u_k^\varepsilon(t)=&   0
                &+ \mathcal{O}(\varepsilon) \qquad \text{ for all } k \neq 1.
\end{eqnarray*}
\end{remark}

    \section{Expansions with even periodic functions
               \label{scos}}

   Fourier expansion  involving only cosines are used throughout this paper.
 There is never sinus.
  In this short section we explain  why it is simple to work with even periodic
 functions
 and we give some hints  to work with  more
 general initial data.
\\

 First, we want to work only with co-sinus to avoid two secular terms.
 If we return to equation (\ref{eq1onedof}):
$ -\alpha_0  ( v_1''  + v_1)  =   (v_0)_+ +  \alpha_1 v_0''$.
 A priory, we have two secular terms in the right hand side,
 one with $\cos(s)$ and another
 with $\sin(s)$. Only one  parameter $\alpha_1$ seems not enough to delete
 all secular terms.
\\
 Otherwise, if $ v_0 \in \mathbb R$,
 $u, \, S $ are  $2\pi$  periodic even functions,
$g \in C^0(\mathbb R,\mathbb R)$  such that
\\
$\displaystyle
 0=  \int_0^{2\pi} e^{is} (  S(s)+g(u(s))) ds
$
 then  the solution of
\begin{eqnarray*}
  v''  + v  &= & S(s)+g(u), \quad v(0) =v_0, \, v'(0)=0,
\end{eqnarray*}
 is necessarily a  $2\pi$  periodic even function.
  Since  we only work with  $2\pi$  periodic even functions
  we have always at most one secular term
 proportional to  $\cos(s)$.

We now investigate the case involving not necessarily even periodic functions.
 In general, $\dot{u}_0^\varepsilon \neq 0 $ and
   $u_\varepsilon$ is the solution of
\begin{eqnarray*}
\ddot{u_\varepsilon}  + u_\varepsilon+\varepsilon f(u_\varepsilon)  =  0,
  \quad u_\varepsilon(0) =u_0^\varepsilon, \, \dot{u_\varepsilon}(0)=\dot{u}_0^\varepsilon.
\end{eqnarray*}
By the energy $2E = \dot{u}^2 + u^2 + \varepsilon F(u)$, where $F'=2f$ and $F(0)=0$,
we know
 that $u_\varepsilon$ is periodic for $\varepsilon$ small enough,
 for instance  with an implicit function theorem see \cite{Ver}
 also valid for Lipschitz function \cite{Cla} in our case.
Denote by $\tau_\varepsilon$
 the first time such that $\dot{u}_\varepsilon(t)=0$. Such time exists
 thanks to the periodicity of $u_\varepsilon$.  Now, let $U_\varepsilon$
 defined by $U_\varepsilon(t)= u_\varepsilon(t+\tau_\varepsilon)$. $U_\varepsilon $ is the solution
 of
 \begin{eqnarray*}
\ddot{U_\varepsilon}  + U_\varepsilon+\varepsilon f(U_\varepsilon)  =  0,
  \quad U_\varepsilon(0) =U_0^\varepsilon=u_\varepsilon(\tau_\varepsilon), \,
 \dot{U_\varepsilon}(0)=0.
\end{eqnarray*}
The initial data $U_0^\varepsilon$ depends on the initial position and
 initial velocity of $u_\varepsilon$ through the energy,
$(U^\varepsilon_0)^2 + \varepsilon F(U_0^\varepsilon)= (u_0^\varepsilon)^2 + (\dot{u}_0^\varepsilon)^2
   + \varepsilon F(u_0^\varepsilon) $.  For instance, if $u_0^\varepsilon$
 and $\dot{u}_\varepsilon^0$ are positive then $U_0^\varepsilon $ is positive
 and
\\
$U_0^\varepsilon= \displaystyle \sqrt{(u_0^\varepsilon)^2 + (\dot{u}_0^\varepsilon)^2
  + \varepsilon (F(u_0^\varepsilon)- F\left(\sqrt{(u_0^\varepsilon)^2 + (\dot{u}_0^\varepsilon)^2} \right) }
 + \mathcal{O}(\varepsilon^2)
          $ .
\\
We can apply the method of strained coordinates for $U_\varepsilon$
only with  even periodic functions:
  $U_\varepsilon (t)= v_0(\omega_\varepsilon t)
 + \varepsilon  v_1(\omega_\varepsilon t) +\mathcal{O}(\varepsilon^2)$.
The expansion obtained for $u_\varepsilon$ by $U_\varepsilon$,
with $\phi_\varepsilon =-\omega_\varepsilon \tau_\varepsilon$  is:
 \begin{eqnarray*}
 u_\varepsilon(t)& =&   v_0(\omega_\varepsilon t +\phi_\varepsilon )
 + \varepsilon  v_1(\omega_\varepsilon t +\phi_\varepsilon) +\mathcal{O}(\varepsilon^2),
\end{eqnarray*}
   which is a  good ansatz in general for $u_\varepsilon$,
  where $v_0$ and $v_1$ are even $2\pi-$periodic functions.
 The method of strained coordinates becomes to find the following
 unknowns  $\phi_0$, $\omega_1, \phi_1$, $\omega_2, \phi_2$ such that
 \begin{eqnarray*}
  \omega_\varepsilon &= &\omega_0 +\varepsilon \omega_1 +\varepsilon^2 \omega_2 + \cdots,\\
   \phi_\varepsilon &=& \phi_0 +\varepsilon \phi_1 +\varepsilon^2 \phi_2 + \cdots.
\end{eqnarray*}
 Indeed,  we have two parameters to delete two secular terms at each step.
 If one is  only interested by the nonlinear frequency shift, it is simpler
 to work only with cosines.
 \\

 Otherwise, if $f$ is an odd function, we can work only with
 odd periodic function. It is often the case in literature when occurs a cubic
 non-linearity. See for instance \cite{KC96,Mi,Nayfeh} for the Duffing equation,
 the Rayleigh equation or  the Korteweg-de Vries equation.

    \section{Appendix: technical proofs \label{sA}}

We  give some useful results  about energy estimates and almost periodic functions in subsection
\ref{ssul}.
 Next we  complete the  proofs for each previous asymptotic expansions
 in subsection \ref{ssbr}.
 The   point is to bound the remainder for large time
 in each expansion.
\subsection{Useful lemmas} \label{ssul}
The following Lemma is useful to prove an expansion
for large   time with non smooth non-linearity.
\begin{lemma}{\bf [Bounds for large time

   ]}
\label{lemmaTime}
 $\mbox{ }$
\\
 Let   $w_\varepsilon$ be a    solution of
 \begin{eqnarray} \label{eqGeneralRemainder}
 \displaystyle \left\{
\begin{array}{l}
  w_\varepsilon''+  w_\varepsilon
 =  S(s) +  f_\varepsilon(s) + \varepsilon g_\varepsilon(s,w_\varepsilon),
    \\
   w_\varepsilon(0)=0, \quad  w_\varepsilon'(0)=0.
 \end{array}
\right.
 \end{eqnarray}
  If  source terms satisfy the following conditions where
 $M > 0$, $C >0$ are fixed constants :
\begin{enumerate}
\item
 $S(s)$ is a $2\pi$-periodic function orthogonal to
  $e^{\pm i s}$,

and   $|S(s)|\leq M$ for all $s$,
  \item
   $|f_\varepsilon|\leq M $ and for all $T$,
  $\displaystyle \int_0^T |f_\varepsilon(s) | ds  \leq C\varepsilon T \quad (\mbox{resp. } C\sqrt{\varepsilon}T$),
 \item
 for all $R> 0$:

 $ \displaystyle M_R=\sup_{\varepsilon \in (0,1),s >0, R > |u|} |g_\varepsilon(s,u)| < \infty$,
\\ that is to say that
  $g_\varepsilon(s,u)$ is locally bounded  with respect to $u$
  \\
  independently from $\varepsilon \in (0,1)$ and $s \in (0,+\infty)$,
\end{enumerate}
then,   there exists $\varepsilon_0 > 0$  and $\gamma > 0$
    such that, for $ 0 < \varepsilon < \varepsilon_0$,
   $ w_\varepsilon $ is uniformly bounded in $
         \displaystyle W^{2,\infty}\left(0,T_\varepsilon \right)$,
  where $T_\varepsilon = \displaystyle \frac{\gamma}{\varepsilon}
  \quad (\mbox{resp. } \displaystyle \frac{\gamma}{\sqrt{\varepsilon}}) $.
\end{lemma}
 Notice that $f_\varepsilon$ and $g_\varepsilon$ are not necessarily continuous.
 Indeed this a case for our asymptotic expansion,
 see Lemme \ref{LDLL} and its applications throughout the paper.
But in previous sections the right hand side is globally continuous,
 i.e. $S+ f_\varepsilon+\varepsilon g_\varepsilon(.,w_\varepsilon)$ is continuous, so, in this case,
 $w_\varepsilon$ is $C^2$.
\\
\medskip
\\
{\bf Proof of the  Lemma \ref{lemmaTime}}:
 First we remove the non resonant periodic source term which is independent of $\varepsilon$. Second,
 we get $L^\infty $ bound for $w_\varepsilon$ and $w_\varepsilon'$ with an energy estimate.
 Third, with equation (\ref{eqGeneralRemainder}), we get
 an uniform estimate for
 $w''_\varepsilon$ in $L^\infty(0,T_\varepsilon)$ and the $W^{2,\infty}$ regularity.
\\
\underline{Step 1: remove $S$}
 \\
  It suffices to  write $w_\varepsilon = w_1 + w_2^\varepsilon$ where
  $ w_1$ solves the linear problem:
  \begin{eqnarray}
  \label{eqw1}
 w_1''+  w_1   =  S(s),
    \quad   w_1(0)=0,  w_1'(0)=0.
 \end{eqnarray}
$w_1$ and $w_1'$ are uniformly bounded
 in   $L^\infty(0,+\infty)$
  since there is no resonance.
\\
More precisely,
 $w_1 = F(s)
 + A \cos(s) + B \cos(s)$,
 where $F$ is  $2\pi$ periodic.
$ F$ is obtained by Fourier expansion without harmonic $n=\pm 1$
 since $S$ is never resonant:
\begin{eqnarray*}
 F(s) &= &
 \sum_{ n \neq \pm 1}
  \frac{c_n}{1 - n^2}e^{ins}
\quad
 with
  \quad
  S(s) =
 \sum_{ n \neq \pm 1} c_n e^{ins}.
\end{eqnarray*}
$F$ is uniformly bounded,
 with Cauchy-Schwartz inequality set
 $C_0^2 = \displaystyle \sum_{ n \neq \pm 1}
  |n^2-1  |^{-2}$,
 we obtain:
\begin{eqnarray*}
 \|F \|_{L^\infty} & \leq &
 \sum_{ n \neq \pm 1}
  \frac{|c_n|}{|n^2-1|}
\leq  C_0  \|S \|_{L^2(0,2\pi)}
\leq  C_0  \|S \|_{L^\infty(0,2\pi)}.
\end{eqnarray*}
Similarly, set   $D_0^2 = \displaystyle \sum_{ n \neq \pm 1}
 n^2 |n^2-1 |^{-2}$,
we have
 $  \|F' \|_{L^\infty} \leq  D_0  \|S \|_{L^\infty(0,2\pi)}$.
\\
Furthermore, $0=w_1(0)= F(0)+ A$,
  and   $0=(w_1)'(0)= F'(0)+ B$,
 then,  $A$ and $B$  are well defined.
 $w_1^,$  is also bounded, i.e.
 there exists $M_1> 0 $ such that
 $\|w_1\|_{W^{1,\infty}(0,+\infty)} \leq M_1$.
\\
Notice that from equation \eqref{eqw1}, $w_1$ belongs in $W^{2,\infty}$.
\\
\medskip
\\
Then we  get  an equation similar to (\ref{eqGeneralRemainder}) for $w_2^\varepsilon$
 with   $S \equiv 0$ and the same assumption for
  the same $f_\varepsilon$ and the new $g_\varepsilon$:
 $
    \overline{g}_\varepsilon(s,w) =
     g_\varepsilon(s,w_1 +w)  .$
 \begin{eqnarray} \label{eqGeneralRemainder2}
 \displaystyle \left\{
\begin{array}{l}
  (w^\varepsilon_2)''+ (w^\varepsilon_2)
 =    f_\varepsilon(s) + \varepsilon \overline{g}_\varepsilon(s,w_2^\varepsilon),
    \\
   (w_2^\varepsilon)(0)=0, \quad  (w_2^\varepsilon)'(0)=0.
 \end{array}
\right.
 \end{eqnarray}
\\
\underline{Step 2: energy estimate}
 \\
 Second, we get an energy estimate for $w_2^\varepsilon$.

   We fix $R> 0$ such that $ R $ is greater than the uniform bound
 $M_1$
 obtained for $w_\varepsilon^1$,
 $ R= M_1 + \rho$ with $\rho > 0$.

Let us define
\begin{eqnarray*}
  2E(s)=((w_2^\varepsilon)'(s))^2 +  (w_2^\varepsilon)(s)^2 ,
 & \quad
 \displaystyle \overline{E}(s) =\sup_{0<\tau < s} E(\tau) ,
\end{eqnarray*}
 and $T_\varepsilon$ be the first time   $T>0$ such that
 $ 2\overline{E}(T)\geq \rho^2$, i.e.
  $\rho$ estimates the size of $(w_2^\varepsilon)$ and $(w_2^\varepsilon)'$.
\\
Multiplying the differential equation (\ref{eqGeneralRemainder2})
 by $(w_2^\varepsilon)'$, we have  for  all $s < T < T_\varepsilon(\rho)$
 the following inequalities
since
 $\displaystyle \sup_{0<\tau < s} |(w_2^\varepsilon)'(\tau)| \leq \sqrt{2\overline{E}(s)}$,
and
$\displaystyle \int_0^T |f_\varepsilon(s) | ds  \leq C\varepsilon T
$,
\begin{eqnarray*}
E(s)  & = &  \displaystyle \int_0^s f_\varepsilon(\tau ) (w_2^\varepsilon)'(\tau) d\tau
                  +
 \varepsilon \int_0^s \overline{g}_\varepsilon(\tau,(w_2^\varepsilon)(\tau))(w_2^\varepsilon)'(\tau) d\tau,
         \\
      &  \leq &
  C \varepsilon s  \sqrt{2\overline{E}(s)}
          +  \varepsilon s  M_R  \sqrt{2\overline{E}(s)}
        ,\\
\overline{E}(T)   &\leq &  C \varepsilon T  \sqrt{2\overline{E}(T)}
          +  \varepsilon T  M_R  \sqrt{2\overline{E}(T)}, \\
\varepsilon T & \geq & \displaystyle  \frac{\sqrt{\overline{E}(T)/2}}{M_R + C}.
\end{eqnarray*}
Notice that if $ 2 \overline{E}(T) < \rho^2 $ for all  $T > 0$
 then $T_\varepsilon= + \infty $. The critical case is
when $T_\varepsilon$ is finite and
  ${\overline{E}(T)}$ approaches $\rho^2/2$ when $T$ goes
  to $T_\varepsilon(\rho)$.
Thus we have
 $ \displaystyle
   \displaystyle T_\varepsilon \geq \frac{\rho}{2\varepsilon(M_R+C)}
$ and $E(t) \le \frac{\rho^2}{2}$ for $t \le T_{\varepsilon}=\frac{\gamma}{\varepsilon}$ with
$\gamma=\frac{\rho}{2(M_R+C}$.
\\
The proof is similar when
$\displaystyle \int_0^s |f_\varepsilon(\tau )|  d\tau \leq C \sqrt{\varepsilon} T $
 then
 $ \displaystyle
   \displaystyle T_\varepsilon \geq \frac{\rho}{2\sqrt{\varepsilon}(\sqrt{\varepsilon}M_R+C)}
$.
\epro
 For completeness, we state  a similar and straightforward
version of Lemma \ref{lemmaTime}
 useful for systems.
\begin{lemma}{\bf [Bounds for large time for systems
                ]}
\label{lemmaTimeSystem}
 $\mbox{ }$
\\
 Let   $w_\varepsilon= (w_1^\varepsilon,\cdots, w_N^\varepsilon)$
 be the     solution of  the following system:
 \begin{eqnarray} \label{eqGeneralRemainderS}
 \displaystyle \left\{
\begin{array}{l}
 (\lambda_1)^2 (w_k^\varepsilon)''+ (\lambda_k)^2 w_k^\varepsilon
 =  S_k(s) +  f_k^\varepsilon(s) + \varepsilon g_k^\varepsilon(s;w_\varepsilon),
    \\
   w_k^\varepsilon(0)=0, \quad  (w_k^\varepsilon)'(0)=0,
 \quad k=1,\cdots,N.
 \end{array}
\right.
 \end{eqnarray}
  If  source terms satisfy the following conditions where
 $M > 0$, $C >0$ are fixed constants :
\begin{enumerate}
\item
      \underline{non resonance conditions}
       with  $S_k(s)$ are $2\pi$-periodic functions and   $|S_k(s)|\leq M$,
\begin{enumerate}
  \item
    $S_1(s)$ is orthogonal to
  $e^{\pm i s}$, i.e. $\displaystyle \int_0^{2\pi} S_1(s)e^{\pm i s}ds =0,$

   \item
    $ \lambda_k, \lambda_1$ are $\mathbb Z$ independent for all $k\neq 1$,
 \end{enumerate}
  \item
   $|f_k^\varepsilon|\leq M $ and for all $T$,
  $\displaystyle \int_0^T |f_\varepsilon(s) | ds  \leq C\varepsilon T \mbox{ or } C\sqrt{\varepsilon}T$,
 \item
 for all $R> 0$:

 $ \displaystyle M_R=\max_k\sup_{\varepsilon \in (0,1),s >0, w_1^2+\cdots+w_N^2< R^2}
 |g_k^\varepsilon(s;u)| < \infty$,
\end{enumerate}
then,   there exists $\varepsilon_0 > 0$  and $\gamma > 0$
    such that, for $ 0 < \varepsilon < \varepsilon_0$,
   $ w_\varepsilon $ is uniformly bounded in $
         \displaystyle W^{2,\infty}\left(0,T_\varepsilon \right)$,
  where $T_\varepsilon = \displaystyle \frac{\gamma}{\varepsilon}
  \mbox{ or } \displaystyle \frac{\gamma}{\sqrt{\varepsilon}} $.
\end{lemma}
\noindent{\rm\bf Proof : }
 First we remove source terms $S_k$ independent of $\varepsilon$
setting $w_k^\varepsilon = w_{k,1}+w_{k,2}^\varepsilon$ where $w_{k,1}$ is the solution of
\begin{eqnarray*}
     \lambda_1^2 w_{k,1}'' + \lambda_k^2 w_{k,1}=S_k, &
    w_{k,1}(0)= 0, &  w_{k,1}'(0)=0.
\end{eqnarray*}
As in the proof of  Lemma \ref{lemmaTime},
$w_{1,1}$ belongs in $W^{2,\infty}$ thanks to
 the non-resonance condition {\it 1.(a)}.
For $k \neq 1$, there is no resonance since
  $\frac{\lambda_k}{\lambda_1} \notin \mathbb Z,$
 i.e.  the non-resonance condition {\it 1.(b)},
 thus a similar expansion also yields
$w_{k,1}$ belongs in $W^{2,\infty}(\mathbb R,\mathbb R)$.
\medskip \\
Now  $w_{k,2}^\varepsilon$ are solutions of the following system for $k=1,\cdots,N$
 \begin{eqnarray*}
 \displaystyle \left\{
\begin{array}{l}
 \lambda_1^2 (w_{k,2}^\varepsilon)''+ \lambda_k^2(w_{k,2}^\varepsilon)
 =    f_k^\varepsilon(s) + \varepsilon \overline{g}_k^\varepsilon(s;w_2^\varepsilon),
    \\
   (w_{k,2}^\varepsilon)(0)=0, \quad  (w_{k,2}^\varepsilon)'(0)=0,
 \end{array}
\right.
 \end{eqnarray*}
with $w_\varepsilon=w_1+w_2^\varepsilon$, $w_2^\varepsilon=(\cdots,w_{k,2}^\varepsilon,\cdots)$ and
 $
    \overline{g}_k^\varepsilon(s;\cdots,w_k,\cdots) =
     g_k^\varepsilon(s;\cdots,w_{k,1} +w_k,\cdots)  .$
\\
The end of the proof of  Lemma \ref{lemmaTimeSystem}
 is a straightforward generalization of the
the proof of  Lemma \ref{lemmaTime} with the
energy:
 $ 2 E(w_1,\cdots,w_N) = \displaystyle
\sum_{k=1}^N \left(
 ( \lambda_1)^2 (\dot{w}_k)^2 + (\lambda_k)^2 w_k^2
 \right)
.$
\epro
\medskip
 For systems,  we also have to work
 with linear combination of periodic functions
with different periods and nonlinear function of such sum.
 So we  work with the adherence in $L^\infty(\mathbb R,\cit)$ of
 span$\{e^{i\lambda t},\, \lambda \in \mathbb R\}$, namely
 the set  of almost periodic functions
$ C^0_{ap}(\mathbb R,\cit)$,
and the  Hilbert space of almost-periodic function
 is  $ L^2_{ap}(\mathbb R,\cit)$,
  see \cite{Cor},
with the scalar product
\begin{eqnarray*}
 \left <u,v \right >
 & = &\lim_{T \rightarrow + \infty}
                \frac{1}{T}\int_0^T u(t) \overline{v(t)}dt.
\end{eqnarray*}
 We  give an useful Lemma about the spectrum of $|u|$
     for  $u \in C^0_{ap}(\mathbb R,\mathbb R)$.
   Let us recall  definitions for the Fourier coefficients
  of $u$ associated to frequency $\lambda$: $ c_{\lambda}[u]$
 and its spectrum: $Sp\;[u]$,
 \begin{eqnarray}
 \label{psap}
  \displaystyle
 c_{\lambda}[u] =  \left < u,e^{i\lambda t} \right >
                  = \lim_{T \rightarrow + \infty}
                \frac{1}{T}\int_0^T u(t) e^{-i\lambda t}dt,
 & \quad &
 Sp\;[u] = \{\lambda \in \mathbb R, \, c_\lambda[u] \neq 0\}.
\end{eqnarray}
 \begin{lemma} \label{LSp}

 {\bf [Property of the  spectrum of $|u|$  ] }  $\mbox{ }$\\
If   $u \in C^0_{ap}(\mathbb R,\mathbb R)$,
 and if   $u$ has  got a finite spectrum:
\\
$Sp\;[u] \subset \{\pm \lambda_1,\cdots, \pm \lambda_N\}$,
with
$ (\lambda_1,\cdots,\lambda_N)$ being  $\mathbb Z$-independent,
 such that $ 0 \notin Sp\;[u]$,
\\
then
  $\lambda_k \notin Sp\;[\; |u| \;]$  for all $k$.
\end{lemma}
\noindent{\rm\bf Proof : }
   The result is quite obvious for $u^2$.
We first prove the result for $f(u^2)$ where $f$ is smooth.
 Then, we conclude by approximating $|u|$ by a smooth sequence
$f_n(u^2)=\sqrt{1/n+u^2}$,
  and using the $L^\infty$ stability of the spectrum.
\\
Let $E$ be the set of all $\mathbb Z$ linear combinations of elements of
$ S_2 = \{0, \pm \lambda_{kj}^\pm, \,  k,j=1,\cdots N\}$,
where
 \begin{eqnarray*}
   \lambda_{kj}^\pm &=& \lambda_k\pm \lambda_{j}.
\end{eqnarray*}
 Thus $ Sp\;[f(u^2)]$ is a subset of $E$   since   $Sp\;[u^2]  \subset S_2$.
\\
Notice that
$
 \lambda_{jk}^\pm=\pm  \lambda_{kj}^\pm,
           \quad
             \lambda_{kk}^{-} =0,
 \quad
  \lambda_{kk}^+ =2\lambda_k=\lambda_{kj}^++ \lambda_{kj}^-.
$ \\
 Choosing $k=1$ for instance, so $\lambda_1 \neq 0$, it suffices to prove
 that $\lambda_1 \notin E$.
\\
Assume the converse,
   i.e., $\lambda_1 \in E$.
 %
 Then, for $k<j$,
  there exists some integers
 $\displaystyle (c_{kj}^\pm)_{k < j)}$ such that:
\begin{eqnarray*}
 \lambda_1  & =&  \sum_{k < j}( c_{kj}^+ \lambda_{kj}^+ +
                         c_{kj}^- \lambda_{kj}^-).
\end{eqnarray*}
 Therefore, defining $c_{jk}^\pm$ by  $\pm c_{kj}^\pm$ for $k < j$,
  we have:
\begin{eqnarray*}
   \lambda_1  &= & \lambda_1  \sum_{ j\neq 1}( c_{1j}^+  +
                         c_{1j}^-)
         + \lambda_2  \sum_{ j\neq 2}( c_{2j}^+
                        +  c_{2j}^-)
       + \cdots + \lambda_N\sum_{ j\neq 1}( c_{k=Nj}^+
                         + c_{Nj}^-).
\end{eqnarray*}
Using the $\mathbb Z$-independence, with $ d_{kj}= c_{kj}^+  +
                         c_{kj}^-$ for $k \neq j$ and $d_{kk}=0$, we have
 following system:
$\displaystyle
  1  =  D_1= \sum_j d_{1j}, \quad
  0  =  D_k= \sum_j d_{kj}, \qquad  \mbox{for all  } k >1.
$
\\
Summing up, the $N-1$ last equations  in $\displaystyle \frac{\mathbb Z}{2\mathbb Z}$,
 and using the fact: $ d_{jk} \equiv d_{kj}$ modulo  $2$,
 we have:
$
\displaystyle
  0 \equiv \sum_{k=2}^{N} D_k
   \equiv \sum_{j=2}^{N} d_{1j} + 2\sum_{k<j}^{} d_{kj}
      \equiv  \sum_{j=2}^{N} d_{1j},
$
then
$ \displaystyle D_1  \equiv 0$, i.e.  $D_1$ is even.
 It's impossible since $D_1=1$.
   So $\lambda_1 \notin E $
and the proof is complete.
 \epro

\subsection{Bounds for the remainders}\label{ssbr}

Now, we prove each asymptotic expansion given in previous sections, i.e.
 we bound each remainders with energy estimates up to a large time.
\medskip
\\
{\bf Proof of  Proposition \ref{propexplicit} :}
 First we give the outline of the proof.
\\
 Notice that all these computations only involve
 the function $\cos$. Then, the only way to have a secular term
in equations defining  $v_1$ and $v_2$ is  a $\cos(s)$
 in the right-hand side. So, the good choice of $\alpha_1$ and $\alpha_2$,
 is enough to remove secular term with $\cos(s)$.
Now, it suffices to control $r_\varepsilon$ for large time.
A computation  shows that
the remainder  $R_\varepsilon$ of equation (\ref{rode}) satisfies:
\\
$|R_\varepsilon(s)|\leq  C \varepsilon( 1 + |r_\varepsilon(s)|)
                 +| \chi_\varepsilon|(v_0,v_1+\varepsilon r_\varepsilon)$.
\\
Then, $r_\varepsilon$ is like $w_\varepsilon$
 in Lemma \ref{lemmaTime}, and the term $f_\varepsilon$ comes from
 $\chi_\varepsilon$ which  is estimated by

  Lemmas \ref{LDLL}, \ref{LHzi}.
\\
More precisely, an exact computation of $R_\varepsilon$ in equation (\ref{rode})
leads to
\begin{eqnarray*}
   R_\varepsilon &= &
  \chi_\varepsilon(v_0,v_1+\varepsilon r_\varepsilon) + \varepsilon H(v_0) r_\varepsilon
 + \varepsilon \alpha_3^\varepsilon v''_\varepsilon,
\end{eqnarray*}
where
  $\alpha_3^\varepsilon$ is  a real constant, bounded uniformly for all
$\varepsilon \in [0,1]$  such that
\\
 $
 (\omega_\varepsilon)^2
 =
 \alpha_0 + \varepsilon \alpha_1+\varepsilon^2 \alpha_2+\varepsilon^3 \alpha_3^\varepsilon.
$
From (\ref{eqchi}) we also have
\begin{eqnarray*}
   & \chi_\varepsilon(v_0,v_1+\varepsilon r_\varepsilon) \\
  = &
 \displaystyle  \frac{(v_0 +\varepsilon v_1+\varepsilon^2 r_\varepsilon )_+
           - [(v_0)_+ + H(v_0)(\varepsilon v_1+\varepsilon^2 r_\varepsilon)  ]}{\varepsilon}
 \\
   =&
  \frac{(v_0 +\varepsilon v_1 )_+
           - [(v_0)_+ +\varepsilon  H(v_0) v_1  ]}{\varepsilon}
  +
\frac{(v_0 +\varepsilon v_1+\varepsilon^2 r_\varepsilon )_+
           - (v_0 +\varepsilon v_1)_+   }{\varepsilon}
 - \varepsilon H(v_0) r_\varepsilon
 \\
=&
 \chi_\varepsilon(v_0,v_1) -\varepsilon H(v_0) r_\varepsilon + \varepsilon \tilde{g}_\varepsilon(s,r_\varepsilon),
\end{eqnarray*}
 since $ u \rightarrow (u)_+$ is 1-Lipschitz
 \begin{eqnarray*}
\tilde{g}_\varepsilon(s,r_\varepsilon)
 & = &\frac{(v_0 +\varepsilon v_1+\varepsilon^2 r_\varepsilon )_+
           - (v_0 +\varepsilon v_1)_+   }{\varepsilon} , \\
 |\tilde{g}_\varepsilon(s,r_\varepsilon) |
 & \leq &
 |r_\varepsilon|.
\end{eqnarray*}
So, with  $v_\varepsilon = v_0+\varepsilon v_1+\varepsilon^2 r_\varepsilon$,
 we can rewrite $R_\varepsilon$ as follow
\begin{eqnarray*}
  R_\varepsilon &= &
  \chi_\varepsilon(v_0,v_1) + \varepsilon  \tilde{g}_\varepsilon(s,r_\varepsilon)
 + \varepsilon \alpha_3^\varepsilon v''_\varepsilon.
\end{eqnarray*}
Now, we can rewrite equation (\ref{rode}) in the following way
\begin{eqnarray*}
- \alpha_0(r''_\varepsilon +  r_\varepsilon)
 & =&
 S(s) +f_\varepsilon(s) + \varepsilon g_\varepsilon(s,r_\varepsilon),
\\
 S & = &
        \alpha_2 v_0'' + \alpha_1 v_1''  + H(v_0) v_1,
 \\
 f_\varepsilon & =&  \chi_\varepsilon(v_0,v_1)
                  + \varepsilon\alpha_3^\varepsilon(v_0"+\varepsilon v_1"),
 \\
 g_\varepsilon &=& \tilde{g}_\varepsilon
                   + \varepsilon^2\alpha_3^\varepsilon r_\varepsilon,
\end{eqnarray*}
which allow us to conclude with Lemma \ref{lemmaTime}.
\epro
The proof for other propositions \ref{PropNdaf}, \ref{Propcritic}
 in section \ref{s1dof} are similar.
\\
We now complete the proof for the asymptotic expansions for systems
given in section \ref{sNdof}.
\medskip \\
{\bf Proof of Theorem \ref{Th41} :}
 As in the proof of  Proposition \ref{propexplicit},
 the same technique of proof is used component by component for
 Theorems \ref{Th41},
 with   similar  energy estimates
 we  can conclude
 with the  Lemma \ref{lemmaTimeSystem} for system to control all
 $r_k^\varepsilon$.
\\
More precisely a complete computation of the remainder gives us:
\begin{eqnarray*}
 R_k^\varepsilon &= &
  \chi_\varepsilon(a_{k1}v_1^0 - b_k,\sum_j a_{kj}
[v_{kj}^1+\varepsilon r_j^\varepsilon] ) + \varepsilon H(a_{k1}v_1^0 - b_k)\sum_j a_{kj} r_j^\varepsilon
 + \varepsilon \alpha_3^\varepsilon (v_k^\varepsilon)'',
\end{eqnarray*}
 with notation of the  proof of Proposition \ref{propexplicit} and
 $v_k^\varepsilon = v_k^0 +\varepsilon v_k^1+\varepsilon^2 r_k^\varepsilon$.
Let $u,v,w$ be three functions, as previously, we have:
\begin{eqnarray*}
\chi_\varepsilon( u,v+\varepsilon w)+ \varepsilon H(u)w &= &
 \chi_\varepsilon( u,v) + \varepsilon^{-1}( (u+\varepsilon v +\varepsilon^2 w)_+ -(u+\varepsilon v )_+ ) ,
\end{eqnarray*}
 and, since $ w \rightarrow w_+$ is 1-Lipschitz:
\begin{eqnarray*}
|\varepsilon^{-1}( (u+\varepsilon v +\varepsilon^2 w)_+ -(u+\varepsilon v )_+ )|
 &\leq &  \varepsilon |w|.
\end{eqnarray*}
Now, we can rewrite $R_k^\varepsilon$ as follow:
\begin{eqnarray*}
 R_k^\varepsilon &= &
  \chi_\varepsilon(a_{k1}v_1^0 - b_k,\sum_j a_{kj}
 v_{kj}^1 ) +
\varepsilon g_k^\varepsilon(s,r_1^\varepsilon,\cdots,r_N^\varepsilon)
 + \varepsilon \alpha_3^\varepsilon (v_k^\varepsilon)'',
\end{eqnarray*}
where $ g_k^\varepsilon$ is defined by
\begin{eqnarray*}
  g_k^\varepsilon(s,r_1^\varepsilon,\cdots,r_N^\varepsilon)
 & = &
 \varepsilon^{-1}\{ ( V_k^\varepsilon  + \varepsilon^2 \sum_j a_{kj}  r_j^\varepsilon )_+
 -
 ( V_k^\varepsilon )_+ \},
\\
V_k^\varepsilon
 & = & a_{k1}v_k^0 -b_k
+ \varepsilon \sum_j a_{kj}
v_{kj}^1,
\end{eqnarray*}
Notice that $g_k^\varepsilon$ satisfies
$ \displaystyle
| g_k^\varepsilon(s,r_1^\varepsilon,\cdots,r_N^\varepsilon)  |
\leq
 \sum_j | a_{kj}| | r_j^\varepsilon |.
$
\\
A key  ingredient is the  energy
$ 2 E = \displaystyle  \sum_k ( \alpha_0 (r_k')^2 + \lambda_k^2 r_k^2)$
 for the homogeneous system: $ L_k r_k = 0,\; k=1,\cdots,N$ and the
 for the inhomogeneous system:
 \begin{eqnarray*}
  - L_k r_k &=& S_k^\varepsilon(s)+ f_k^\varepsilon(s) + \varepsilon g_k(s,r_1^\varepsilon,\cdots,r_N^\varepsilon),
  \end{eqnarray*}
 for $ k=1,\cdots,N$,
with
 \begin{eqnarray*}
   S_k  &= &H(a_{k1}v_1^0-b_k) \sum_{j=1}^{N} a_{kj}v^1_j
            + \alpha_2 v_k^0" + \alpha_1 v_k^1",
\end{eqnarray*}
 and $\alpha_1,\, \alpha_2$ are well chosen to avoid secular term when $k=1$.
 Thus, all $S_k$ are $2\pi$ periodic.
 $S_1$ is not resonant with $L_1$.
 The $\lambda_k$ are $\mathbb Z $ independent.
 We can apply  Lemma \ref{lemmaTimeSystem}
 which is enough to conclude the proof.
 \epro
\medskip  $\mbox{ }$ \\
{\bf Proof of Theorem \ref{Th43} :}
 The proof follows two steps. First the solution for linear equations
 \eqref{linsysacc} are bounded by $o(t)$.
Second, energy estimates are used to bound $r^\varepsilon$.
\\
At the end we prove remark
                  \ref{polarisation}.
\\
Notice that we do not use  Lemmas \ref{LDLL}, \ref{LHzi}.
Indeed,  we have no term with $\chi_\varepsilon$.
 We only use that functions  $u_+$ and $v_k^0$ are Lipschitz,
 the Lemma \ref{LSp} to identify resonant terms when $b_k=0$ and an
  energy estimate.
But, since all modes are excited, the accuracy is weaker
than the precision obtained in Theorem \ref{Th41}, as in \cite{SP}.
\medskip \\
\underline{ Step 1}: the $N$ problems \eqref{linsysacc} involves
                decoupled equations rewritten as follow with  $\omega > 0$,
       \begin{eqnarray*}
             \phi''(s)+\omega^2\phi(s) &= & S(s) \in C^0_{ap}(\mathbb R,\mathbb R),
   \qquad  \pm \omega \notin Sp[S].
        \end{eqnarray*}
There is no resonance since $\pm \omega$ are not in the spectrum of $S$.
But, $Sp[S]$ is dense in $\mathbb R$.
Indeed $\lambda_1,\cdots,\lambda_N$ are $\mathbb Z$ independent.
In general, we cannot expect that $\phi$ is bounded on the real line,
 see \cite{Cor}, but $\phi $ is less than $\mathcal{O}(s)$ for large time.
 We can compute explicitly $\phi$
\begin{eqnarray*}
   \phi(s) & = & A \cos(\omega s) + B \sin(\omega s) + \psi(s), \\
   \omega \psi(s) & =& \displaystyle \int_0^s S(\sigma) \sin(\omega(s-\sigma)) d\sigma
           \\ &  = &
 \displaystyle \sin(\omega s) \int_0^s S(\sigma) \cos(\omega\sigma) d\sigma
                  -\cos(\omega s) \int_0^s S(\sigma) \sin(\omega\sigma) d\sigma
.\end{eqnarray*}
 The condition $\pm \omega \notin Sp[S]$ is
 $\displaystyle \lim_{s \rightarrow + \infty}
 s^{-1}\int_0^s S(\sigma) \exp(\pm i\omega\sigma) d\sigma =0$.
 That is to say
\\
 $\displaystyle \int_0^s S(\sigma) \exp(\pm i\omega\sigma) d\sigma =~o(s)$
 when $s\rightarrow + \infty$, thus $\psi$ and $\phi$
are negligible compared to $s$ for large time.
\medskip \\
\underline{ Step 2}: Let us decompose the remainder in the following way
 $r_k^\varepsilon = \phi_k + w_k^\varepsilon$. From equation \eqref{RFose} and the previous
 bound for $\phi_k$ we have in variable $t$ instead of $s_k$ for convenience
 \begin{eqnarray*}
      L_k w_k^\varepsilon(t) & = & \mathcal{O}(\varepsilon t) +
    \left(  \mathcal{O}(\varepsilon \phi_k) + \mathcal{O}(\varepsilon |w^\varepsilon|))
          =  \mathcal{O}(\varepsilon t) + \mathcal{O}(\varepsilon |w^\varepsilon|),
    \right)
 \end{eqnarray*}
 since $\phi_k(t)=o(t).$ Now, we remove the first part of the right hand side
 with $w_k^\varepsilon= \tilde{w}_k^\varepsilon+z_k^\varepsilon$ and $ \tilde{w}_k^\varepsilon$ is
 solution of $L_k  \tilde{w}_k^\varepsilon= \mathcal{O}(\varepsilon t)$.
 Classical energy estimates (or explicit computations as for $\phi$)
 yields to
   $ \tilde{w}_k^\varepsilon(t)= \mathcal{O}(\varepsilon t^2)$.
 Thus there exists a constant $C_1> 0$ such that $z_k^\varepsilon$ satisfies
 \begin{eqnarray*}
  \left | L_k z_k^\varepsilon \right| & \leq &
        C (\varepsilon^2 t^2 + \varepsilon |z^\varepsilon|).
 \end{eqnarray*}
Multiplying each inequality by $|(z_k^\varepsilon)'|$, summing up with respect to $k$,
integrating on $[0,T]$,  by Cauchy-Schwarz inequality,
 with $D= 2C (\min(\lambda_k)+ \min(\lambda_k)^2)$ we get
\begin{eqnarray*}
  E(T) & = & \displaystyle \sum_{k=1}^{N}
                    \left( \lambda_1^2 ((z_k^\varepsilon)')^2
                          + \lambda_k^2 (z_k^\varepsilon)^2
                      \right)\\
     & \leq & 2 C \varepsilon^2 T^2 \displaystyle \int_0^T \sum_{k=1}^{N}|(z_k^\varepsilon)'|(t)dt
               + 2C\varepsilon\displaystyle \int_0^T |(z^\varepsilon)' \cdot z^\varepsilon| dt\\
    & \leq & D \varepsilon^2 T^{2.5} \displaystyle\sqrt{ \int_0^T  E(t)dt}
             + D\varepsilon\displaystyle \int_0^T E(t) dt.
\end{eqnarray*}
Let $Y(T)$ be $\displaystyle\int_0^T E(t) dt $, thus $Y(0)=0$ and for all
 $ t \in [O,T]$,
\begin{eqnarray*}
   E(t)=Y'(t) &\leq& D \varepsilon^2 T^{2.5} \sqrt{Y(t)} + D\varepsilon Y(t).
\end{eqnarray*}
 Since
$\displaystyle \int^Y_0 \frac{dy}{A\sqrt{y}+y} =
 2 \ln \left(1 + \frac{\sqrt{y}}{A} \right) $
we obtain $\sqrt{Y(T)} \leq \varepsilon T^{2.5} \exp(D\varepsilon T)$ and
 then
 \begin{eqnarray*} E(T) &\leq& 2D \varepsilon^3 T^5 \exp(D\varepsilon T).
 \end{eqnarray*}
Finally $r_k^\varepsilon = \phi_k+ \tilde{w}_k^\varepsilon+z_k^\varepsilon
 =  o(T) +\mathcal{O}(\varepsilon T^2)+\mathcal{O}(\varepsilon^{1.5} T^{2.5} \exp(D\varepsilon T)) $,
 so for any $T_\varepsilon = o(\varepsilon^{-1}) $
we have in $W^{1,\infty}(0,T_\varepsilon)$ for all $T \leq T_\varepsilon$
\begin{eqnarray*}
      \varepsilon r_\varepsilon(T)  &= & o(\varepsilon T_\varepsilon) + \mathcal{O}(\varepsilon^2 T_\varepsilon^2)
                    +\mathcal{O}(\varepsilon^{2.5} T_\varepsilon^{2.5}),
\end{eqnarray*}
which is enough  to have the convergence in $W^{1,\infty}(0,T_\varepsilon)$.
Furthermore $r_k^\varepsilon$ satisfies  the second order differential equation
\eqref{eqNdofRemainderH}  which is enough to get the convergence
 in $W^{2,\infty}$.
\medskip \\
\underline{ About remark  \ref{polarisation}}:
 From Theorem \ref{Th43}, this result its obvious.
 Let us explain why we cannot go further up to the order $\varepsilon^2$.

Unfortunately $S_k$ is not periodic since $v_j^1$ is quasi-periodic
 for $j\neq 1$.
 Indeed, the following  initial conditions
\begin{eqnarray*}
  v_k^1(0)=0, &  (v_k^1)'(0)=0,&   k \neq 1,
\end{eqnarray*}
yields to  a quasi-periodic function,
sum of two periodic functions with different periods
$2\pi$ and $2\pi \lambda_1/\lambda_k$, thus a globally bounded function
\begin{eqnarray*}

    v_k^1(s) &=& \phi_k^1(s)
   - \phi_k^1(0)\cos\left(\frac{\lambda_k}{\lambda_1}s\right).
\end{eqnarray*}
 So we cannot apply
Lemma \ref{lemmaTimeSystem}.
 \\
 Let us decompose $S_k=P_k+Q_k$  for $k\neq 2$ where $P_k$ is periodic and
  $Q_k$ is almost-periodic
 \begin{eqnarray*}
   Q_k(s) &=& -H(a_{k1}v_1^0(s)-b_k) \sum_{j=1}^{N} a_{kj}
                \phi^1_j(0)\cos\left(\frac{\lambda_j}{\lambda_1}s\right).
 \end{eqnarray*}
Let $w_k$ be a solution of $-L_k w_k = Q_k$ then
$\displaystyle Sp[w_k] \in \bigcup_{j}
 \left\{\pm \frac{\lambda_j}{\lambda_1} + \mathbb Z  \right\}$,
 so the spectrum of $w_k$ is discrete  and there is resonance
 in the $N-1$ equations,
  $-L_k r^\varepsilon_k=S_k+\cdots $, $k\neq 1$ and the expansion does not still
 valid for time of the order $\varepsilon^{-1}$.
 \epro

   %

{\it Acknowledgments}: we thank Alain L\'eger,  Vincent Pagneux and
 St\'ephane Roux
 for their valuable remarks at the
fifth meeting of the GDR US,  Anglet, 2008.
We also  thank G\'erard Iooss for fruitful discussions.

 \bibliographystyle{plain}

redpajama_set_name: RedPajamaArXiv


text:
class Freshmeat < ApplicationRecord
  validates :name, presence: true

  validates :version, presence: true
end

redpajama_set_name: RedPajamaGithub


UnknownViewType 0 The view type is not known.
FrontView 1 Front side view of the part.
TopView 2 Top side view of the part.
BackView 3 Back side view of the part.
BottomView 4 Bottom side view of the part.
EndView 5 End side view of the part.
SectionView 6 Section view of the part.
ModelView 7 Model view of the part.
DetailView 8 Detail view of the part.
_3DView 9 3d view of the part.