\begin{table}[!ht]
  \begin{tabular}{|l|l|}
    \hline
    author  & V. Acary\\
    \hline
    date    & Sept, 20, 2011 \\ 
    \hline
    version &  \\
    \hline
  \end{tabular}
\end{table}



This section is devoted to the implementation and the study  of the algorithm. The interval of integration is $[0,T]$, $T>0$, and a grid $t_{0}=0$, $t_{k+1}=t_{k}+h$, $k \geq 0$, $t_{N}=T$ is constructed. The approximation of a function $f(\cdot)$ on $[0,T]$ is denoted as $f^{N}(\cdot)$, and is a piecewise constant function, constant on the intervals $[t_{k},t_{k+1})$. We denote $f^{N}(t_{k})$ as $f_{k}$. The time-step is $h>0$. 


\section{Various second  order dynamical systems with input/output relations}



\subsection{Lagrangian dynamical systems}


The class {\tt LagrangianDS}  defines  and computes a generic ndof-dimensional 
Lagrangian Non Linear Dynamical System of the form :

\begin{equation}
  \begin{cases}
    M(q,z) \dot v + N(v, q, z) + F_{Int}(v , q , t, z) = F_{Ext}(t, z) + p \\
    \dot q = v
  \end{cases}
\end{equation}
 where 
 \begin{itemize}
 \item  $q \in R^{ndof} $ is the set of the generalized
   coordinates, 
 \item $ \dot q =v \in R^{ndof} $ the velocity,
   i. e. the time derivative of the generalized coordinates
   (Lagrangian systems).
 \item $ \ddot q =\dot v \in R^{ndof} $ the
   acceleration, i. e. the second time derivative of the generalized
   coordinates.  
 \item $ p \in R^{ndof} $ the reaction forces due to
   the Non Smooth Interaction.  
 \item $ M(q) \in R^{ndof \times ndof}
   $ is the inertia term saved in the SiconosMatrix mass.  
 \item $
   N(\dot q, q) \in R^{ndof}$ is the non linear inertia term saved
   in the {\tt SiconosVector \_NNL}.  
 \item $ F_{Int}(\dot q , q , t) \in
   R^{ndof} $ are the internal forces saved in the SiconosVector
   fInt.  
 \item $ F_{Ext}(t) \in R^{ndof} $ are the external forces
   saved in the SiconosVector fExt.  
 \item $ z \in R^{zSize}$ is a
   vector of arbitrary algebraic variables, some sort of discrete
   state.
 \end{itemize}

 
  The equation of motion is also shortly denoted as:
  \begin{equation}
  M(q,z) \dot v = F(v, q, t, z) + p
\end{equation}
 
  where  $F(v, q, t, z) \in R^{ndof} $ collects the total forces
  acting on the system, that is 
  \begin{equation}
    F(v, q, t, z) =  F_{Ext}(t, z) -  NNL(v, q, z) + F_{Int}(v, q , t, z) 
\end{equation}

 This vector is stored in the  {\tt SiconosVector \_Forces  }  

\subsection{Fully nonlinear case}
Let us introduce the following system,
\begin{equation}
  \label{eq:FullyNonLinear}
  \begin{cases}
    M(q,z) \dot v = F(v, q, t, z) + p  \\
    \dot q = v \\
    y = h(t,q,\lambda) \\
    p = g(t,q,\lambda)
  \end{cases}
\end{equation}
where $\lambda(t) \in \RR^m$  and $y(t) \in \RR^m$ are  complementary variables related through a multi-valued mapping. According to the class of systems, we are studying, the function $F$ , $h$ and $g$ are defined by a fully nonlinear framework or by affine functions. This fully nonlinear case is not  implemented in Siconos yet. This fully general case is not yet implemented in Siconos.



\subsection{Lagrangian Rheonomous relations}

\begin{equation}
  \label{eq:RheonomousNonLinear}
  \begin{cases}
    M(q,z) \dot v = F(v, q, t, z) + p \\
    \dot q = v \\
    y = h(t,q) \\
    p = G(t,q)\lambda)
  \end{cases}
\end{equation}

\subsection{Lagrangian Scleronomous relations}

\begin{equation}
  \label{eq:ScleronomousNonLinear}
  \begin{cases}
    M(q,z) \dot v  = F(v, q, t, z) + p  \\
    \dot q = v \\
    y = h(q) \\
    p = G(q)\lambda
  \end{cases}
\end{equation}


\paragraph{Fully Linear case}

\begin{equation}
  \label{eq:FullyLinear}
  \begin{cases}
    M \dot v   +C v + Kq = F_{Ext}(t, z) + p  \\
    \dot q = v \\
    y = C q + e + D\lambda  + F z \\
    p = C^T\lambda
  \end{cases}
\end{equation}




\section{Moreau--Jean event-capturing scheme} 

In this section, a time-discretization method of the Lagrange dynamical equation (\ref{eq:11}), consistent with the nonsmooth character of the solution, is presented. It is assumed in this section, as in the other sections, that $v^+(\cdot)=\dot{q}^{+}(\cdot)$  is a locally bounded variation function. The equation of motion reads as,
\begin{equation}
  \label{eq:11-b}
  \begin{cases}
    M(q(t)) {dv} +N(q(t),v^{+}(t)) dt+  F_{\mathrm{int}}(t, q(t), v^+(t))\,dt = F_{\mathrm{ext}}(t)\,dt + dr \\ \\
   v^+(t)=\dot{q}^+(t) \\ \\
  q(0)=q_{0} \in {\mathcal C}(0),\;\dot{q}(0^{-})=\dot{q}_{0}
  \end{cases}    
\end{equation}
We also assume that $F_{\mathrm{int}}(\cdot)$ and $F_{\mathrm{ext}}(\cdot)$ are continuous with respect to time. This assumption is made for the sake of simplicity to avoid the notation $F^+_{\mathrm{int}}(\cdot)$ and $F^+_{\mathrm{ext}}(\cdot)$. Finally, we will condense the nonlinear inertia terms and the internal forces to lighten the notation. We obtain 
\begin{equation}
  \label{eq:11-c}
  \begin{cases}
    M(q(t)) {dv} + F(t, q(t), v^+(t))\,dt = F_{\mathrm{ext}}(t)\,dt + dr \\ \\
   v^+(t)=\dot{q}^+(t) \\ \\
  q(0)=q_{0} \in {\mathcal C}(0),\;\dot{q}(0^{-})=\dot{q}_{0}
  \end{cases}    
\end{equation}

The NSCD method, also known as the Contact Dynamics (CD) is due to the seminal works of J.J.~\cite{Moreau1983,Moreau1985,Moreau1988,Moreau1994,Moreau1999} and M.~\cite{Jean88,Jean1999}  (See also \citep{Jean.Pratt85,Jean.Moreau91,Jean.Moreau92}). A lot of improvements and variants have been proposed over the years. In this Section,  we take  liberties with these original works, but we  choose to present a version of the NSCD method which preserves the essential of the original work. Some  extra developments and interpretations are added which are only under our responsibility. To come back to the source of the NSCD method, we encourage to read the above references.

\subsection{The Linear Time-invariant NonSmooth Lagrangian Dynamics}
\label{section11.1.1}



 For the sake of simplicity of the presentation, the linear time-invariant case is considered first. The nonlinear case will be examined later in this chapter. 
\begin{equation}
  \label{eq:11-a}
  \begin{cases}
    M dv + (K q(t) + C v^+(t))\,dt = F_{\mathrm{ext}}(t)\,dt + dr  \\ \\
    v^+(t)=\dot{q}^+(t)
  \end{cases}
\end{equation}


\subsubsection{Time--discretization of the Dynamics}


Integrating both sides of this equation over a time step $(t_k,t_{k+1}]$ of length $h>0$, one obtains
\begin{eqnarray}
  \begin{cases}
    \displaystyle \int_{(t_k,t_{k+1}]} M dv + \int_{t_k}^{t_{k+1}} (C v^+(t)
      + K q(t)) \,dt = \displaystyle \int_{t_k}^{t_{k+1}} F_{\mathrm{ext}}\,dt +
        \displaystyle \int_{(t_k,t_{k+1}]} dr \:, \\ \\
     q(t_{k+1}) = q(t_{k}) + \displaystyle \int_{t_k}^{t_{k+1}} v^+(t)\,dt 
   \end{cases}
\end{eqnarray}

By definition of the differential measure $dv$, we obtain
\begin{eqnarray}
\label{eq:19}
&  \displaystyle \int_{(t_k,t_{k+1}]} M \,dv = M \int_{(t_k,t_{k+1}]}\,dv = M\,(v^+(t_{k+1})-v^+(t_{k}))  &
\end{eqnarray}
Note that the right velocities are involved in this formulation.  The impulsion $\displaystyle \int_{(t_k,t_{k+1}]} dr$ of the reaction on the time interval $(t_k,t_{k+1}]$ emerges  as a natural unknown. The equation of the nonsmooth motion can be written under an integral form as:
\begin{eqnarray}
  \begin{cases}
     M\,(v(t_{k+1})-v(t_{k})) =   \displaystyle   \int_{t_k}^{t_{k+1}} (- C v^+(t)
      - K q(t) +  F_{\mathrm{ext}}(t))\,dt +
        \displaystyle \int_{(t_k,t_{k+1}]} dr \:, \\ \\
     q(t_{k+1}) = q(t_{k}) + \displaystyle \int_{t_k}^{t_{k+1}} v^+(t)\,dt 
   \end{cases}
\end{eqnarray}
%
Choosing a numerical method boils down to choose a method of approximation for the remaining integral terms. Since discontinuities of the derivative  $v(\cdot)$ are to be expected if some shocks are occurring, \ie{}. $dr$  has some  atoms within the interval $(t_k,t_{k+1}]$, it is not relevant to use high order approximations  integration schemes for $dr$ (this was pointed out in remark \ref{remark1023}). It may be shown on some examples that, on the contrary, such high order schemes may generate artefact numerical oscillations (see \citep{Vola.Pratt.ea98}). 


%% \paragraph{Numerical Examples}
%% \begin{ndrva}
%% complete with few numerical illustrations of the oscillations
%% \end{ndrva}


The following notation will be used: 

\begin{itemize}
\item $q_{k}$ is an approximation of $q(t_{k})$ and $q_{k+1}$ is an approximation of $q(t_{k+1})$, 

\item $v_{k}$ is an approximation of $v^+(t_{k})$ and $v_{k+1}$ is an approximation of $v^+(t_{k+1})$, 

\item $p_{k+1}$ is an approximation of $ \displaystyle \int_{(t_k,t_{k+1}]} \,dr$. 
\end{itemize}
%
A popular first order numerical scheme, the so called $\theta$-method, is used for the term supposed to be sufficiently smooth:\index{$\theta$-method}
\begin{eqnarray}
  \displaystyle \int_{t_k}^{t_{k+1}} C v + K q \,dt  &\approx& 
  h \left[ \theta (C v_{k+1}+K q_{k+1}) + (1-\theta) (C v_{k}+K q_{k}) \right]   \nonumber \\
  \displaystyle \int_{t_k}^{t_{k+1}} F_{\mathrm{ext}}(t) \,dt &\approx& 
  h\left[\theta  (F_{\mathrm{ext}})_{k+1}+(1-\theta)  (F_{\mathrm{ext}})_{k}  \right]  \nonumber 
\end{eqnarray}
The displacement, assumed to be absolutely continuous, is approximated by:
\begin{eqnarray}
&  q_{k+1} = q_{k} +  h\,\left[\theta v_{k+1}+(1-\theta) v_{k}  \right] & \nonumber
\end{eqnarray}
Taking into account all these discretizations, the following time-discretized equation of motion is obtained:
\begin{equation}
\label{eq:NSCD-discret}
\begin{cases}
    M (v_{k+1}-v_{k}) + h\left[\theta  (C  v_{k+1}+K q_{k+1}) + (1-\theta) (C v_{k}+K q_{k})  \right] = \\ \\
    \quad\quad\quad\quad\quad = h\left[\theta  (F_{\mathrm{ext}})_{k+1}+(1-\theta)  (F_{\mathrm{ext}})_{k}  \right] + p_{k+1} \\  \\
    q_{k+1} = q_{k} +  h\left[\theta v_{k+1}+(1-\theta) v_{k} \right]
\end{cases}
\end{equation}
Finally, introducing the expression of $q_{k+1}$ in the first equation of~(\ref{eq:NSCD-discret}), one obtains:
\begin{eqnarray}
  \label{eq:23}
&  \left[M+h\theta C + h^2 \theta^2 K\right] (v_{k+1} -v_{k}) = - h  C v_{k} - h K q_{k} - h^2 \theta  K v_{k} & \nonumber \\ \nonumber \\
&+  h\left[\theta  (F_{\mathrm{ext}})_{k+1})+(1-\theta)  (F_{\mathrm{ext}})_{k}  \right]  + p_{k+1}  \:, &
\end{eqnarray}
which can be written as:
\begin{eqnarray}
  \label{eq:24}
   v_{k+1} = v_{\mathrm{free}}  + \widehat{M}^{-1} p_{k+1}
\end{eqnarray}
where,
%% \begin{ndrmj}
%% OK, j'ai remis $\widehat{M} = \left[M+h\theta C + h^2 \theta^2 K \right] $ et j'ai mis $\widehat{M}^{-1}$ a la place
%% de $\widehat{W}$.
%% \end{ndrmj}
\begin{itemize}
\item the matrix   
  \begin{equation}
\widehat{M} = \left[M+h\theta C + h^2 \theta^2 K \right]  \label{eq:2002}
\end{equation}
is usually called the iteration matrix\index{Iteration matrix}.
\item The vector   
\begin{equation}
 \label{eq:2003}
\begin{array}{ll}
v_{\mathrm{free}}  & = v_{k} + \widehat{M}^{-1} \left[   - h  C v_{k} - h K q_{k} - h^2 \theta  K v_{k} \right. \\ \\ 
& \left. +  h\left[ \theta  (F_{\mathrm{ext}})_{k+1})+(1-\theta)  (F_{\mathrm{ext}})_{k} \right] \right] 
\end{array}
\end{equation}
%
is the so-called ``free'' velocity, \ie{}, the velocity of the system when reaction forces are null.     
\end{itemize}


\subsubsection{Comments} 

Let us make some  comments on the above developments: 

\begin{itemize}

\item The iteration matrix $ \widehat{M} = \left[M+h\theta C + h^2 \theta^2 K \right] $ is supposed to be invertible, 
since the mass matrix $M$ is usually positive definite and $h$ is supposed to be small enough. 
The matrices $C$ and $K$ are  usually semi-definite positive since rigid motions are allowed to bodies.

\item  When $\theta=0$, the $\theta$-scheme is the explicit Euler scheme. When $\theta=1$, the $\theta$-scheme is the fully
implicit Euler scheme. When dealing with a plain ODE
\begin{equation}
    M\ddot{q}(t)  + C \dot{q}(t) + K q(t)  = F(t) 
\end{equation}
the $\theta-$scheme is unconditionally stable for $0.5 < \theta \leq 1$. It is conditionally stable otherwise. 


\item The equation (\ref{eq:24}) is a linear form of the dynamical equation. It appears 
as an affine relation between the two unknowns, $v_{k+1}$ that is an approximation of the right derivative of the Lagrange variable 
at time $t_{k+1}$, and the impulse $p_{k+1}$. Notice that this scheme is fully implicit. Nonsmooth laws have to be treated by implicit methods. 


\item From a numerical point of view, two major features appear. First, the different terms in 
the numerical algorithm will keep finite values. When the time step $h$ vanishes, the scheme copes with finite jumps. 
Secondly, the use of differential measures of the time interval $(t_k,t_{k+1}]$, \ie{}., 
$dv((t_{k},t_{k+1}])=v^+(t_{k+1})-v^+(t_{k})$ and $dr((t_{k},t_{k+1}])$, 
offers a rigorous treatment of the nonsmooth evolutions.  It is to be noticed that approximations of the acceleration are ignored. 

\end{itemize}

These remarks on the contact dynamics method might be viewed only as some numerical tricks. In fact, the
mathematical study of the second order MDI by Moreau provides a sound mathematical ground to this numerical scheme. 
It is noteworthy that convergence results have been proved for such time-stepping schemes \cite{Marques1993,Stewart1998,Mabrouk1998,dzonou2007}, see below.

\subsection{The Nonlinear  NonSmooth Lagrangian Dynamics}
\label{section11.1.2}

\subsubsection{Time--discretization of the Dynamics}
Starting from the nonlinear dynamics~(\ref{eq:11-c}), the integration of  both sides of this equation over a time step $(t_k,t_{k+1}]$ of length $h>0$ yields\begin{eqnarray}
  \begin{cases}
    \displaystyle \int_{(t_k,t_{k+1}]} M(q) dv + \int_{t_k}^{t_{k+1}} F(t, q(t), v^+(t)) \,dt = \displaystyle \int_{t_k}^{t_{k+1}} F_{\mathrm{ext}}(t)\,dt +
        \displaystyle \int_{(t_k,t_{k+1}]} dr \:, \\ \\
     q(t_{k+1}) = q(t_{k}) + \displaystyle \int_{t_k}^{t_{k+1}} v^+(t)\,dt 
   \end{cases}
\end{eqnarray}
The first term is generally approximated by
\begin{equation}
\label{eq:19-NL}
  \displaystyle \int_{(t_k,t_{k+1}]} M(q) \,dv \approx  M(q_{k+\gamma})\,(v_{k+1}-v_{k}) 
\end{equation}
where $q_{k+\gamma}$ generalizes the standard notation for $\gamma \in [0,1]$ such that
\begin{equation}
  \label{eq:NL1}
  q_{k+\gamma} = (1-\gamma) q_{k} + \gamma\,  q_{k+1}
\end{equation}
%\begin{ndrva}
%  Is there a equivalent of the mean-value theorem  for differential measure ? 
%\end{ndrva}
The \textit{a priori} smooth terms are evaluated with a $\theta$-method, chosen in this context for its energy conservation ability,
\begin{eqnarray}
  \displaystyle \int_{t_k}^{t_{k+1}} F(t,q,v) \,dt  &\approx& 
  h  \tilde F_{k+\theta} 
\end{eqnarray}
where $\tilde F_{k+\theta}$ is an approximation with the following dependencies
$$ \tilde F(t_k,q_k,v_k,t_{k+1},q_{k+1},v_{k+1},t_{k+\theta},q_{k+\theta},v_{k+\theta}) $$
The mid-values $t_{k+\theta},q_{k+\theta},v_{k+\theta}$ are defined by
\begin{equation}
  \label{eq:NSCD-discret-b}
  \left\{\begin{array}{l}
  t_{k+\theta} = \theta t_{k+1}+(1-\theta) t_{k}\\
  q_{k+\theta} = \theta q_{k+1}+(1-\theta) q_{k}\\
  v_{k+\theta} = \theta v_{k+1}+(1-\theta) v_{k}
  \end{array}\right.,\quad  \theta \in [0,1]
\end{equation}


\begin{remark} \label{eq:Simo}
  The choice of the approximated function $\tilde F(\cdot)$ strongly depends
  on the nature of the internal forces that are modeled. For the
  linear elastic behavior of homogeneous continuum media, this
  approximation can be made by:
\begin{equation}
\tilde F_{k+\theta} = \frac 1 2 K\contract\left[E(q_{k})+E(q_{k+1})\right] \contract F(q_{k+1/2})
\end{equation}
where $E(:cdot)$ is the Green-Lagrange strain tensor, which leads to an energy conserving algorithm as in
\citep{Simo.Tarnow92}. For nonlinear elastic other smooth nonlinear
behaviors, we refer to the work of
\citep{Gonzalez2000,Laursen.Meng2001} and references therein for the choice of the
discretization and the value of $\theta$.
\end{remark} 

The displacement, assumed to be absolutely continuous is approximated by:
\begin{eqnarray}
&  q_{k+1} = q_{k} +  h\,v_{k+\theta}  & \nonumber
\end{eqnarray}


The following nonlinear time--discretized equation of motion is obtained:
\begin{equation}
\label{eq:NSCD-discret-nl}
\begin{cases}
    M(q_{k+\gamma}) (v_{k+1}-v_{k}) + h \tilde F_{k+\theta} = p_{k+1} \\  \\
    q_{k+1} = q_{k} +  h v_{k+\theta}
\end{cases}
\end{equation}
In its full generality and at least formally, substituting the expression of $q_{k+\gamma},q_{k+1}$ and $q_{k+\theta}$,  the first line of the  problem can be written under the form of a residue $\mathcal R$ depending only on $v_{k+1}$ such that 
\begin{equation}
  \label{eq:NL3}
  \mathcal R (v_{k+1}) = p_{k+1}
\end{equation}
In the last expression, we have omitted the dependence to the known values at the beginning the time--step, \ie{} $q_k$ and $v_k$.

\subsubsection{Linearizing the Dynamics}

The system of equations~(\ref{eq:NL3}) for $v_{k+1}$ and $p_{k+1}$ can be linearized yielding a Newton's procedure  for solving it. This linearization needs the knowledge of the Jacobian matrix  $\nabla \mathcal R (\cdot)$ with respect to its argument to construct the tangent linear model.

 Let us consider that the we have to solve the following equations,
\begin{equation}
  \label{eq:NL4}
  \mathcal R (u) = 0 
\end{equation}
by a Newton's method where
\begin{equation}
  \label{eq:NL6}
    \mathcal R (u) =   M(q_{k+\gamma} ) (v_{k+1}-v_{k}) + h \tilde F_{k+\theta}
\end{equation}
 The solution of this system of nonlinear equations is sought as a limit of the sequence $\{ u^{\tau}_{k+1}\}_{\tau \in \nbN}$ such that
 \begin{equation}
   \label{eq:NL7}
   \begin{cases}
     u^{0}_{k+1} = v_k \\ \\
     \mathcal R_L( u^{\tau+1}_{k+1}) =  \mathcal R (u^{\tau}_{k+1}) + \nabla \mathcal R (u^{\tau}_{k+1} )(u^{\tau+1}_{k+1}-u^{\tau}_{k+1} ) =0
 \end{cases}
\end{equation}
 In practice, all the nonlinearities are not treated in the same manner and the Jacobian matrices for the nonlinear terms involved in the Newton's algorithm are only computed in their natural variables. In the following, we consider some of the most widely used approaches.



\paragraph{The Nonlinear Mass Matrix}
The derivation of the Jacobian of the first term of  $\mathcal R (\cdot)$ implies to compute
\begin{equation}
  \label{eq:NL2000}
   \nabla_u  \left(M(q_{k+\gamma}(u) ) (u-v_{k})\right) \text{ with } q_{k+\gamma}(u) = q_k + \gamma h[(1-\theta) v_k+ \theta u].
\end{equation}
One gets
\begin{equation}
  \label{eq:NL8}
  \begin{array}{ll}
    \nabla_u  \left(M(q_{k+\gamma}(u) ) (u-v_{k})\right) &=   M(q_{k+\gamma}(u))  + \left[ \nabla_u M(q_{k+\gamma}(u) ) \right] (u-v_{k}) \\ \\
                                         &=    M(q_{k+\gamma}(u)) + \left[h \gamma\theta \nabla_{q} M(q_{k+\gamma}(u))\right]  (u-v_{k}) 
\end{array}
\end{equation}

\begin{remark}
The notation $\nabla_{u}M(q_{k+\gamma}(u))(u-v_{k})$ is to be understood as follows: 

$$\nabla_{u}M(q_{k+\gamma}(u))(u-v_{k})=\frac{\partial}{\partial u}[M(q_{k+\gamma}(u))(u-v_{k})]$$

which is denoted as $\frac{\partial M_{ij}}{\partial q^{l}}(q_{k+\gamma}(u))(u^{l}-v_{k}^{l})$ in tensorial notation.
\label{remarkBABAS}
\end{remark}



A very common approximation consists in considering that the mass matrix evolves slowly with the configuration in a single time--step, that is, the term $\nabla_{q} M(q_{k+\gamma})$ is neglected and one gets,
\begin{equation}
  \label{eq:NL9}
    \nabla_u  (M(q_{k+\gamma}(u) ) (u-v_{k})) \approx  M(q_{k+\gamma}(u) )
\end{equation}
The Jacobian matrix $\nabla \mathcal R (\cdot)$ is evaluated in $u^{\tau}_{k+1}$ which yields for the  equation~(\ref{eq:NL9})
\begin{equation}
  \label{eq:NL10}
    \nabla_u  (M(q_{k+\gamma} ) (u^{\tau}_{k+1}-v_{k})) \approx  M(q_k + \gamma h [(1-\theta)v_k+\theta u^{\tau}_{k+1}] ) )
\end{equation}
The prediction of the position which plays an important role will be denoted by
\begin{equation}
  \label{eq:NL555}
  \tilde q^{\tau}_{k+1}= q_k + \gamma h [(1-\theta)v_k+\theta u^{\tau}_{k+1}] 
\end{equation}


Very often, the  matrix $M(q_{k+\gamma})$ is only  evaluated at the first Newton's iteration with $u^{0}_{k+1}= v_k$ leading the approximation for the whole step:
\begin{equation}
M(q_k + \gamma h [(1-\theta)v_k+\theta u^{\tau}_{k+1}] ) )\approx M(q_k + h \gamma v_k)
\label{eq:NL11}
\end{equation}
Another way to interpret the approximation~(\ref{eq:NL11}) is to remark that this evaluation is just an explicit evaluation of the predictive position~(\ref{eq:NL555}) given by $\theta=0$:
\begin{equation}
  \label{eq:NL5}
  \tilde q_{k+1}= q_k + h \gamma v_k
\end{equation}

Using this prediction, the problem~(\ref{eq:NSCD-discret-nl}) is written as follows:
\begin{equation}
\label{eq:NSCD-discret2}
\begin{cases}
    M(\tilde q_{k+1}) (v_{k+1}-v_{k}) + h \tilde F_{k+\theta} = p_{k+1} \\  \\
    q_{k+1} = q_{k} +  h v_{k+\theta} \\ \\
    \tilde q_{k+1}= q_k + h \gamma v_k
\end{cases}
\end{equation}


%% \begin{remark}
%%   Rapid change of the mass matrix == > Solver Hairer.
%% \end{remark}


\paragraph{The Nonlinear Term $F(t,q,v)$}
The remaining nonlinear term is  linearized providing the Jacobian matrices of $F(t,q,v)$ with respect to $q$ and $v$. This expression depends strongly on the choice of the approximation $\tilde F_{k+\theta}$. Let us consider a pedagogical example, which is not necessarily the best as the Remark~\ref{eq:Simo} suggests but which is one of the simplest,
\begin{equation}
  \label{eq:NL13}
  \tilde F_{k+\theta} = (1-\theta) F(t_k,q_k,v_k) + \theta F(t_{k+1},q_{k+1},v_{k+1}) 
\end{equation}
The computation of the Jacobian of $  \tilde F_{k+\theta}(t,q(u),u)$ for $$q(u) = q_k+h[(1-\theta)v_k+\theta u]   $$ is given for this example by
\begin{equation}
  \label{eq:NL12}
  \begin{array}{ll}
    \nabla_u  \tilde F_{k+\theta}(t,q,u) &= \theta \nabla_u  F(t,q(u),u) \\ \\
    &= \theta \nabla_q F(t_{k+1},q(u)   ,u) \nabla_{u} q(u) + \theta \nabla_{u} F(t,q(u),u)    \\ \\
    &= h \theta^2 \nabla_q F(t, q(u)   ,u) + \theta \nabla_{u} F(t,q(u),u) \\   
  \end{array}
\end{equation}
The standard tangent stiffness and damping matrices $K_t$ and $C_t$ are defined by
\begin{equation}
  \label{eq:NL14}
  \begin{array}{ll}
  K_t(t,q,u) &= \nabla_q F(t, q   ,u) \\ \\
  C_t(t,q,u) &= \nabla_u F(t, q   ,u) \\
\end{array}  
\end{equation}
In this case, the  Jacobian of $  \tilde F_{k+\theta}(t,q(u),u)$ may be written as 
\begin{equation}
  \label{eq:NL15}
  \begin{array}{ll}
    \nabla_u  \tilde F_{k+\theta}(t,q,u) &=  h \theta^2  K_t(t,q,u) + \theta C_t(t, q   ,u)  \\   
  \end{array}
\end{equation}

The complete Newton's iteration can then be written as 
\begin{equation}
  \label{eq:NL16}
   \widehat M^{\tau+1}_{k+1} (u^{\tau+1}_{k+1}-u^{\tau}_{k+1})  =  \mathcal R (u^{\tau}_{k+1}) +p^{\tau+1}_{k+1}
\end{equation}
where the iteration matrix is evaluated as
\begin{equation}
 \widehat M^{\tau+1}_{k+1} = (M(\tilde q^{\tau}_{k+1}) +  h^2 \theta^2  K_t(t_{k+1},q^{\tau}_{k+1},u^{\tau}_{k+1}) + \theta h C_t(t, q^{\tau}_{k+1}   ,u^{\tau}_{k+1}))\label{eq:NL17}
\end{equation}

(compare with (\ref{eq:2002})). 

\begin{remark}
  The choice of $\theta=0$ leads to an explicit evaluation of the position and the nonlinear forces terms. This choice can be interesting if the time--step has to be chosen relatively small due  to the presence a very rapid dynamical process. This can be the case in crashes applications or in fracture dynamics~\citep{Acary-Monerie2006}. In this case, the iteration matrix reduces to $\widehat M^{\tau+1}_{k+1} = M(\tilde q^{\tau}_{k+1})$ avoiding the expensive evaluation of the tangent operator at each time--step. 

This choice must not be misunderstood. The treatment of the nonsmooth dynamics continues to be implicit.  
\end{remark}

\section{Schatzman--Paoli 'scheme and its linearizations}


\subsection{The scheme}
\begin{subnumcases}{}
  M(q_{k})(q_{k+1}-2q_{k}+q_{k-1})  - h^2 F(v_{k+\theta}, q_{k+\theta}, t_{k+theta})  =  p_{k+1},\quad\,\\ \notag\\ 
  v_{k+1}=\Frac{q_{k+1}-q_{k-1}}{2h}, \\ \notag \\
  y_{k+1} = h\left(\Frac{q_{k+1}+e q_{k-1}}{1+e}\right) \\
  p_{k+1}= G\left(\Frac{q_{k+1}+e q_{k-1}}{1+e}\right) \lambda_{k+1} \\
  0 \leq y_{k+1}  \perp\lambda_{k+1} \geq 0 .
\end{subnumcases}




\begin{ndrva}
Should we have 
  $$ v_{k+1}=\Frac{q_{k+1}-q_{k-1}}{2h}$$ or  $$ v_{k+1}=\Frac{q_{k+1}-q_{k}}{h}$$ ? This question is particularly important for the initialization and the proposed $\theta$-scheme
\end{ndrva}
\subsection{The Newton linearization}

Let us define the residu on $q$
\begin{equation}
  \label{eq:residu}
  \mathcal R(q) =   M(q)(q-2q_{k}+q_{k-1})  + h^2 F( (\theta v(q)+ (1-\theta) v_k),\theta q+ (1-\theta) q_k),  t_{k+\theta})  -  p_{k+1}
\end{equation}
with 
\begin{equation}
  \label{eq:residu-linq1}
  v(q) = \Frac{q-q_{k-1}}{2h}
\end{equation}
that is
\begin{equation}
  \label{eq:residu-linq2}
  \mathcal R(q) =   M(q)(q-2q_{k}+q_{k-1})  + h^2 F( (\theta \Frac{q-q_{k-1}}{2h} + (1-\theta) v_k),\theta q+ (1-\theta) q_k),  t_{k+\theta})   -  p_{k+1}
\end{equation}

Neglecting $\nabla_q  M(q)$ we get 
\begin{equation}
  \label{eq:iterationmatrix}
 \nabla_q \mathcal R(q^\nu) =   M(q^\nu) + h^2  \theta K(q^\nu,v^\nu) + \Frac 1 2 h  \theta C(q^\nu,v^\nu)
\end{equation}
and we  have to solve
\begin{equation}
  \label{eq:iterationloop}
 \nabla_q \mathcal R(q^\nu)(q^{\nu+1}-q^\nu) = -  \mathcal R(q^\nu) .
\end{equation}



\subsection{Linear version of the scheme}


\begin{subnumcases}{}
  M(q_{k+1}-2q_{k}+q_{k-1})  + h^2 (K q_{k+\theta}+ C v_{k+\theta})  =  p_{k+1},\quad\,\\ \notag\\ 
  v_{k+1}=\Frac{q_{k+1}-q_{k-1}}{2h}, \\ \notag \\
  y_{k+1} = h\left(\Frac{q_{k+1}+e q_{k-1}}{1+e}\right) \\
  p_{k+1}= G\left(\Frac{q_{k+1}+e q_{k-1}}{1+e}\right) \lambda_{k+1} \\
  0 \leq y_{k+1}  \perp\lambda_{k+1} \geq 0 .
\end{subnumcases}

Let us define the residu on $q$
\begin{equation}
  \label{eq:residu-linq}
  \mathcal R(q) =   M(q-2q_{k}+q_{k-1})  + h^2 (K(\theta q+ (1-\theta) q_k))+ C (\theta v(q)+ (1-\theta) v_k))  -  p_{k+1}
\end{equation}
with 
\begin{equation}
  \label{eq:residu-linq1}
  v(q) = \Frac{q-q_{k-1}}{2h}
\end{equation}
that is
\begin{equation}
  \label{eq:residu-linq2}
  \mathcal R(q) =   M(q-2q_{k}+q_{k-1})  + h^2 (K(\theta q+ (1-\theta) q_k)))+  h^2 C (\theta \Frac{q-q_{k-1}}{2h}+ (1-\theta) v_k))  -  p_{k+1}
\end{equation}

In this linear case, assuming that $q^0=q^\nu = q_k$, we get
\begin{equation}
  \label{eq:residu-linq2}
  \mathcal R(q^\nu) =   M(-q_{k}+q_{k-1})  + h^2 (K q_k)+  h^2 C (\theta \Frac{q_k-q_{k-1}}{2h}+ (1-\theta) v_k))  -  p_{k+1}
\end{equation}


\section{What about mixing {\tt OnestepIntegrator} in Simulation?}
\label{Sec:MisingOSI}
Let us consider that we have two simple linear Lagrangian Dynamical systems
\begin{equation}
  \label{eq:FullyLinear1}
  \begin{cases}
    M_1 \dot v_1  = F_{1,Ext}(t) + p_1   \\
    \dot q_1 = v_1 
  \end{cases}
\end{equation}
and
\begin{equation}
  \label{eq:FullyLinear1}
  \begin{cases}
    M_2 \dot v_2   = F_{2,Ext}(t) + p_2  \\
    \dot q_2 = v_2 \\
  \end{cases}
\end{equation}
These Dynamical systems (\ref{eq:FullyLinear1}) and (\ref{eq:FullyLinear1}) might numerically solved by choosing two different time--stepping schemes. Let us choose for instance Moreau's scheme for(\ref{eq:FullyLinear1}) 
\begin{equation}
  \label{eq:FullyLinear1-TS}
  \begin{cases}
    M_1 (v_{1,k+1}-v_{1,k})  = F_{1,Ext}(t_{k+1}) + p_{1,k+1}   \\
    q_{1,k+1} = q_{k}+ h  v_{1,k+\theta} 
  \end{cases}
\end{equation}
and Schatzman--Paoli's sheme for (\ref{eq:FullyLinear1}) 
\begin{equation}
  \label{eq:FullyLinear1-TS}
  \begin{cases}
    M_2(q_{2,k+1}-2q_{2,k}+q_{2,k-1})  = F_{2,Ext}(t_{k+1}) + p_{2,k+1}  \\
    v_{2,k+1} = \Frac{q_{2,k+1}-q_{2,k-1}}{2h} \\
  \end{cases}
\end{equation}


Let us consider known that we have a {\tt LagrangianLinearTIR} between this two DSs such that
\begin{equation}
  \label{eq:LTIR-2DS}
  \begin{array}{l}
  y = q_1-q_2 \geq 0 \\ \\
  p = \left[
  \begin{array}{c}
    1 \\
    -1
  \end{array}\right] \lambda
\end{array}
\end{equation}
and a complementarity condition
\begin{equation}
  \label{eq:CP}
  0\leq y \perp \lambda \geq 0
\end{equation}
Many questions are raised when we want to deal with the discrete systems:
\begin{itemize}
\item Which rules should we use for the discretization of~(\ref{eq:CP}) ?
  \begin{equation}
    \label{eq:CP-TS1}
    \text{ if } \bar y_{k+1}\leq 0, \text{ then }  0\leq \dot y _{k+1} + e \dot y_{k} \perp \hat \lambda_{k+1}\geq 0 
  \end{equation}
  or
  \begin{equation}
    \label{eq:CP-TS2}
    0\leq y _{k+1} + e y_{k-1} \perp \tilde \lambda_{k+1}\geq 0 
  \end{equation}
\item Should we assume that $y_{k+1} = q_{1,k+1}-q_{2,k+1}$ and $\dot y_{k+1} = v_{1,k+1}-v_{2,k+1}$
\item How can we link $\hat \lambda_{k+1}$ and  $\tilde \lambda_{k+1}$ with $p_{1,k+1}$ and $p_{2,k+1}$ ?
\end{itemize}

The third is the more difficult question and is seems that it is not reasonable to deal with two DS related by one interaction with different osi.In practice, this should be avoided in Siconos.




%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "DevNotes"
%%% End: