# Event-Capturing schemes¶

## General Principle¶

Roughtly speaking, the event-capturing, a.k.a. time-stepping, method consists in the time-discretisation of the whole system (dynamics + relations + non-smooth laws), leading to a so-called one-step non smooth problem (OSNSP) solved at each time step.

Indeed, the main stages of the process are:

• integrate the dynamics without constraints, to get some “free” solutions
• formalize and solve a OSNSP (a LCP for example)
• update the dynamics with the OSNSP solutions to get the full state.

The discretization process for different dynamical systems, relations and laws is described thereafter. A summary of all the results can be found in section Summary of the time discretized equations

Notations:

In the following sections, the systems are integrated over a time step $$[t_i,t_{i+1}]$$ of constant size $$h$$. The approximation of any function $$F(t,...)$$ at the time $$t_i$$ is denoted $$F_i$$. Note that in the relations writings, upper case letters are used for all variables related to DynamicalSystem objects: $$X , Q, \ldots$$ are concatenation of $$x, q,\ldots$$ of the dynamical systems variables concerned by the relation.

## First order systems¶

### Time Discretisation of the Dynamics¶

#### First Order Non Linear Systems¶

$\begin{split}M\dot x(t) &= f(x,t,z) + r \\ x(t_0) &= x_0\end{split}$

with $$r = r^d = \sum_{\alpha} r^{\alpha}, \alpha \in I_d$$, $$I_d$$ being the set of all relations in which the current dynamical system, number $$d$$, is involved. In the following, the index “d” will be omitted to lighten notations.

The integration of the ODE over a time step $$[t_i,t_{i+1}]$$ of length $$h$$ is :

$M\int_{t_i}^{t_{i+1}}\dot x\,dt = \int_{t_i}^{t_{i+1}} f(t,x,z)dt + \int_{t_i}^{t_{i+1}}r\,dt$

The left-hand term is $$M(x(t_{i+1})-x(t_i)) \approx M(x_{i+1} - x_i)$$ .

Right-hand terms are approximated with a $$\theta$$-method:

$\begin{split}\int_{t_i}^{t_{i+1}} f(t,x,z)dt &\approx h \theta f(t_{i+1},x_{i+1},z) + h (1-\theta) f(t_i,x_i,z) \\ &\approx h \theta f_{i+1} + h (1-\theta) f_i\end{split}$

and the third integral is approximated with:

$\int_{t_i}^{t_{i+1}}r\,dt \approx h r(t_{i+1}) \approx hr_{i+1}$

Then, we get the following “residu”

$\begin{split}\mathcal R(x_{i+1}) &= M(x_{i+1}-x_i) - h \theta f_{i+1} - h (1-\theta) f_{i} - hr_{i+1} = 0 \\ &= \mathcal R^{free}(x_{i+1}) - hr_{i+1}\end{split}$

Note: We introduce the “free” notation for terms related to the smooth part of the system.

A Newton method is used to solve $$\mathcal R(x_{k+1}) = 0$$. The gradient of the residu according to $$x$$ is:

$\nabla_{x}\mathcal R(x) = M - h \theta\cdot\nabla_{x}f(t,x)$

And we get (index k corresponds to the Newton iteration number):

$W_{i+1}^k\cdot (x_{i+1}^{k+1} - x_{i+1}^k) = - \mathcal R(x_{i+1}^k)$

with

$W_{i+1}^k = M - h \theta\left[\nabla_{x}f\right](t_{i+1},x_{i+1}^k)$

If we assume that $$W_{i+1}^k$$ is invertible, we get the solution at Newton iteration k+1:

$\begin{split}x_{i+1}^{k+1} &= x_{i+1}^k - (W_{i+1}^k)^{-1}\mathcal R^{free}(x_{i+1}^{k}) + h(W_{i+1}^k)^{-1}r_{i+1}^{k+1} \\ &= x^{free,k}_{i+1} + h(W_{i+1}^k)^{-1}r_{i+1}^{k+1}\end{split}$

#### First Order Linear Systems¶

$\begin{split}M\dot x(t) &= A(t,z)x(t) + b(t) + r \\ x(t_0) &= x_0\end{split}$

For the integration of the ODE over a time-step, we proceed as in the previous section for non-linear systems to get:

$\begin{split}\mathcal R(x_{i+1}) &= M(x_{i+1}-x_i) - h \theta(A_{i+1}x_{i+1} + b_{i+1})- h (1-\theta)(A_{i}x_i + b_i) - hr_{i+1} = 0 \\ or \\ (M - h\theta A_{i+1}) x_{i+1} &= (M + h (1-\theta)A_{i})\cdot x_i + h\theta(b_{i+1}-b_i) + hb_i + hr_{i+1}\end{split}$

We denote:

$W_{i+1} = (M - h\theta A_{i+1})$

and assuming it is invertible, we get:

$\begin{split}x_{i+1} &= W_{i+1}^{-1}\left[(M + h (1-\theta)A_{i})\cdot x_i + h\theta(b_{i+1}-b_i) + hb_i\right] + hW_{i+1}^{-1}r_{i+1} \\ &= x^{free}_{i+1} + hW_{i+1}^{-1}r_{i+1}\end{split}$

#### First Order Linear Systems with time invariant coefficients¶

$\begin{split}M\dot x(t) &= Ax(t) + b + r \\ x(t_0) &= x_0\end{split}$

Using the results of the previous section, the discretisation is straightforward:

$\begin{split}x_{i+1} &= x_i + h W^{-1}(A x_i + b) + hW^{-1}r_{i+1} \\ &= x^{free}_{i} + hW^{-1}r_{i+1}\end{split}$

with a W that does not depend on time:

$W = (M - h\theta A)$

### Time discretization of the relations¶

In the following, $$R$$ represents the concatenation of all $$r^{\alpha}$$ vectors for the DS involved in the present relation.

#### First Order (non-linear) Relations¶

$\begin{split}y &= h(X,t,\lambda,Z)\\ R &= g(X,t,\lambda,Z)\\\end{split}$

Then, for the iteration $$k+1$$ of the Newton process, we get:

$\begin{split}y_{i+1}^{k+1} &= h(X_{i+1}^{k+1},t_{i+1},\lambda_{i+1}^{k+1})\\ R_{i+1}^{k+1} &= g(X_{i+1}^{k+1},t_{i+1},\lambda_{i+1}^{k+1})\end{split}$

These constraints are linearized around state $$(X_{i+1}^{k+1},\lambda_{i+1}^{k+1})$$:

$\begin{split}y_{i+1}^{k+1} &= y_{i+1}^k - H_0(S_{i+1}^k)X_{i+1}^{k} - H_1(S_{i+1}^k)\lambda_{i+1}^{k} + H_0(S_{i+1}^k)X_{i+1}^{k+1} + H_1(S_{i+1}^k)\lambda_{i+1}^{k+1} \\ \\ R_{i+1}^{k+1} &= R_{i+1}^k - G_0(S_{i+1}^k)X_{i+1}^{k} - G_1(S_{i+1}^k)\lambda_{i+1}^{k} + G_0(S_{i+1}^k)X_{i+1}^{k+1} + G_1(S_{i+1}^k)\lambda_{i+1}^{k+1}\end{split}$

Where $$S_{i+1}^k$$ stands for $$(X_{i+1}^{k},t_{i+1},\lambda_{i+1}^{k})$$ and

$\begin{split}H_0(X,t,\lambda)=\nabla_X h(X,t,\lambda)&, \ \ H_1(X,t,\lambda)=\nabla_{\lambda} h(X,t,\lambda) \\ &\\ G_0(X,t,\lambda)=\nabla_X g(X,t,\lambda)&, \ \ G_1(X,t,\lambda)=\nabla_{\lambda} g(X,t,\lambda) \\\end{split}$

In the case where :

$x_{i+1}^{k+1} = x^{free,k}_{i+1} + (w_{i+1}^k)^{-1}r_{i+1}^{k+1}$

We can write

$X_{i+1}^{k+1} = X^{free,k}_{i+1} + (W_{i+1}^k)^{-1}R_{i+1}^{k+1}$

where $$(W_{i+1}^k)^{-1}$$, is a diagonal block matrix holding the $$(w_{i+1}^k)^{-1}$$, then, if there is one and only one interaction we have:

$(1-(W_{i+1}^k)^{-1}G_{0,i+1}^k) X_{i+1}^{k+1} = X_{i+1}^{free,k} + (W_{i+1}^k)^{-1} (R_{i+1}^k - G_{0,i+1}^k X_{i+1}^k - G_{1,i+1}^k \lambda_{i+1}^k + G_{1,i+1}^k \lambda_{i+1}^{k+1})$

and finally:

$\begin{split}y_{i+1}^{k+1} &= M_{lcp}\lambda_{i+1}^{k+1} + q_{lcp} \\ M_{lcp} &= H_{1,i+1}^k + H_{0,i+1}^k (1-(W_{i+1}^k)^{-1} G_{0,i+1}^k)^{-1} (W_{i+1}^k)^{-1} G_{1,i+1}^k \\ q_{lcp} &= y_{i+1} -H_{0,i+1}^k X_{i+1}^k - H_{1,i+1}^k \lambda_{i+1}^k + H_{0,i+1}^k (1-(W_{i+1}^k)^{-1} G_{0,i+1}^k)^{-1} [X_{i+1}^{free,k} + (W_{i+1}^k)^{-1} (R_{i+1}^k - G_{0,i+1}^k X_{i+1}^k - G_{1,i+1}^k \lambda_{i+1}^k)]\end{split}$

#### First Order Linear Relations¶

$\begin{split}y &= C(t,Z)X(t) + F(t,Z)Z + D(t,Z)\lambda + e(t,Z) \\ R &= B(t,Z) \lambda\end{split}$

Note: for time-invariant relations, B, C, F, D and e are constant vectors and matrices </em>

The Time discretization of the relations is fully implicit and may be written as :

$\begin{split}y_{i+1} &= C(t_{i+1})X_{i+1} + D(t_{i+1})\lambda_{i+1} + e(t_{i+1}) + F(t_{i+1})Z \\ \\ R_{i+1} &= B(t_{i+1})\lambda_{i+1}\end{split}$

### Discretisation of the non-smooth law¶

#### Complementarity Condition¶

The complementarity condition writes:

$0 \leq y \, &\perp \, \lambda \geq 0$

and the discretisation is straightforward:

$0 \leq y_{i+1} \, &\perp \, \lambda_{i+1} \geq 0$

## Lagrangian systems¶

### Time Discretisation of the Dynamics¶

#### Lagrangian (second order) Non Linear Systems¶

We provide in the following sections a time discretization method of the Lagrangian dynamical systems, consistent with the non smooth character of the solution.

$\begin{split}M(q(t),z) dv &= f_L(t,v^+(t), q(t), z)dt + dr \\ v^+(t) &= \dot q^+(t) \\ q(t_0) &= q_0 \\ \dot q(t_0^-) &= v_0\end{split}$

with

$q(t) = q_0 + \int_{t_0}^t v^+(t)dt$

Remark: recall that $$v^+(t)$$ means $$v(t^+)$$ ie right limit of $$v$$ in t.

Left hand side is discretised by assuming that:

$\int_{t_i}^{t_{i+1}} M(q(t),z)dv \approx M(q*,z)(v_{i+1}-v_{i})$

As for first order non-linear systems, we use a $$\theta$$-method to integrate the other terms, and obtain:

$\int_{t_i}^{t_{i+1}} f_L(t, v^+(t), q(t), z) dt \approx h\theta f_L(t_{i+1}, v_{i+1}, q_{i+1}, z) + h(1-\theta) f_L(t_{i}, v_{i}, q_{i}, z)$

and for the last term, we set a new variable $$p_{i+1}$$ such that:

$\int_{t_i}^{t_{i+1}} dr \approx p_{i+1}$

Finally the full system discretisation results in:

$\begin{split}\mathcal R(v_{i+1}, q_{i+1}) &= M(q*,z)(v_{i+1}-v_{i}) - h\theta {f_L}_{i+1} - h(1-\theta) {f_L}_{i} - p_{i+1} = 0 \\ &= \mathcal R^{free}(v_{i+1},q_{i+1}) - p_{i+1}\end{split}$

The “free” notation still stands for terms related to the smooth part of the system. The displacement is integrated through the velocity with :

$q_{i+1} &\approx q_i + h\theta v_{i+1} + h(1 - \theta)v_{i}$

Substituing this into the residu leads to a function depending only on $$v_{i+1}$$, since state “i” and “k” are supposed to be known.

A Newton method will be applied to solve $$\mathcal R(v_{i+1}) = 0$$.

That requires to compute the gradient of the residu; assuming that the mass matrix evolves slowly with the configuration in a single time step, we get:

$\nabla_{v_{i+1}}\left[M(q*,z)(v_{i+1}-v_{i})\right] \approx M(q^{*},z)$

and denoting:

$\begin{split}C_t(t,v,q)=-\left[\frac{\partial{f_L(t,v,q)}}{\partial{v}}\right] \\ \\ K_t(t,v,q)=-\left[\frac{\partial{f_L(t,v,q)}}{\partial{q}}\right]\end{split}$

we get (index k corresponds to the Newton iteration number):

$W(t_{i+1}^k,v_{i+1}^k,q_{i+1}^k)\cdot (v_{i+1}^{k+1}-v_{i+1}^k) = - \mathcal R(v_{i+1}^k)$

with

$W(t,v,q) = M(q*,z) + h\theta C_t(t,v,q) + h^2\theta^2 K_t(t,v,q)$

As an approximation for $$q^*$$, we choose:

$\begin{split}q^* &\approx (1-\gamma) q_i + \gamma q_{i+1}^k \\ &\approx q_i + h\gamma\left[ (1-\theta) v_i + \theta v_{i+1}^k\right]\end{split}$

with $$\gamma \in \left[0,1\right]$$. Moreover, if $$M$$ is evaluated at the first step of the Newton iteration, with $$v_{i+1}^0 = v_i$$, we get:

$M(q^*) \approx M(q_i + h\gamma v_i)$

Finally, if $$W$$ is invertible, the solution at iteration k+1 is given by,

$\begin{split}v_{i+1}^{k+1} &= v_{i+1}^k - (W_{i+1}^k)^{-1} \mathcal R^{free}(v_{i+1}^k) + (W_{i+1}^k)^{-1} p_{i+1}^{k+1} \\ &= v^{free,k}_{i+1} + (W_{i+1}^k)^{-1} p_{i+1}^{k+1}\end{split}$

#### Lagrangian (second order) Linear Systems with Time Invariant coefficients¶

$\begin{split}M dv + Cv^+(t) + K q(t) &= F_{ext}(t,z) + p \\ q(t_0) &= q0 \\ \dot q(t_0^-) &= v_0\end{split}$

Proceeding in the same way as in the previous section, with $$M$$ constant and $$f_L(t,v^+(t), q(t), z) = F_{ext}(t) - Cv^+(t) - Kq(t)$$, integration is straightforward:

$\mathcal R(v_{i+1}, q_{i+1}) &= M(v_{i+1}-v_{i}) - h\theta\left[ F_{ext}(t_{i+1}) - Cv_{i+1} - K q_{i+1}\right] - h(1-\theta)\left[ F_{ext}(t_{i}) - Cv_{i} - K q_{i}\right] - p_{i+1} = 0$

Using the displacement integration through the velocity,

$\begin{split}q_{i+1} = q_{i} + h\left[\theta v_{i+1}+(1-\theta) v_{i} \right]\\\end{split}$

we get:

$W(v_{i+1}-v_{i}) &= (- hC - h^2\theta K )v_{i} - h K q_{i} + h\left[\theta F_{ext}(t_{i+1})+(1-\theta) F_{ext}(t_{i}) \right] + p_{i+1}$

with $$W$$ a constant matrix:

$W = \left[M + h\theta C + h^2 \theta^2 K \right]$

and if $$W$$ is invertible,

$\begin{split}v_{i+1} &= v_{i} + W^{-1}\left[(- hC - h^2\theta K )v_{i} - h K q_{i}+ h\theta F_{ext}(t_{i+1})+h(1-\theta) F_{ext}(t_{i}) \right] + W^{-1} p_{i+1} \\ &= v^{free}_i + W^{-1} p_{i+1}\end{split}$

The free velocity $$v^{free}$$ correponds to the velocity of the system without any constraints.

### Time discretization of the relations¶

#### Lagrangian Scleronomous Relations¶

$\begin{split}y &= h(Q,Z) \\ \dot y &= G_0(Q,Z)V \\ P &= G_0^t(Q,Z)\lambda\end{split}$

with

$\begin{split}G_0(Q) &= \nabla_Qh(Q) \\\end{split}$

From now on, to lighten the notations, the parameter $$Z$$ will omitted.

Considering the Newton process introduced above for Lagrangian non linear systems, the constraints write:

$\begin{split}\dot y_{i+1}^{k+1} = G_0(Q_{i+1}^{k+1}))V_{i+1}^{k+1} \\ P_{i+1}^{k+1} = G_0^t(Q_{i+1}^{k+1}))\lambda_{i+1}^{k+1}\end{split}$

To evaluate $$G_0$$ we still use the prediction $$Q^*$$ defined in the previous section:

$Q^*( V_{i+1}^{k+1}) = Q_i + h\gamma \left[ (1-\theta) V_i + \theta V_{i+1}^{k+1} \right]$

Then we get:

$\begin{split}\dot y_{i+1}^{k+1} = G_0(Q^*(V_{i+1}^{k+1}))V_{i+1}^{k+1} \\ \\ P_{i+1}^{k+1} = G_0^t(Q^*(V_{i+1}^{k+1}))\lambda_{i+1}^{k+1}\end{split}$

These constraints are linearized around the point $$V_{i+1}^{k}$$ and we neglect the second order terms in the computation of the jacobians. It leads to:

$\begin{split}\dot y_{i+1}^{k+1} = G_0(Q^*(V_{i+1}^k))V_{i+1}^{k+1} \\ \\ P_{i+1}^{k+1} = G_0^t(Q^*(V_{i+1}^k))\lambda_{i+1}^{k+1}\end{split}$

As for the evaluation of the mass, the prediction of the position, $$Q^*$$ can be evaluated at the first iteration of the Newton process,

$Q^*(V_{i+1}^0) = Q_i + h\gamma V_i$

#### Lagrangian Rheonomous Relations¶

$\begin{split}y &= h(Q,t) \\ \dot y &= G_0(Q,t)V + G_1(Q,t) \\ P &= G_0^t(Q,t)\lambda \\ with\\ G_0(Q,t) &= \nabla_Qh(Q,t) \\ G_1(Q,t) &= \frac{\partial{h(Q,t)}}{\partial{t}} \\\end{split}$

As for scleronomous relations, we get:

$\begin{split}\dot y_{i+1}^{k+1} &= G_0(Q^*(V_{i+1}^k),t_{i+1})V_{i+1}^{k+1} + G_1(Q^*(V_{i+1}^k, t_{i+1})) \\ \\ P_{i+1}^{k+1} &= G_0^t(Q^*(V_{i+1}^k),t_{i+1})\lambda_{i+1}^{k+1}\end{split}$

#### Lagrangian Compliant Relations¶

$\begin{split}y &= h(Q,\lambda(t)) \\ \dot y &= G_0(Q,\lambda(t))V + G_1(Q,\lambda(t))\dot\lambda \\ P &= G_0^t(Q,\lambda(t))\lambda with\\ G_0(Q,\lambda(t)) &= \nabla_Qh(Q,\lambda(t)) \\ G_1(Q,\lambda(t)) &= \nabla_\lambda h(Q,\lambda(t)) \\\end{split}$

Following the same process as in the paragraph above, it comes:

$\begin{split}\dot y_{i+1}^{k+1} &= G_0(Q^*(V_{i+1}^k),\lambda_{i+1}^k)V_{i+1}^{k+1} + G_1(Q^*(V_{i+1}^k, \lambda_{i+1}^k))\lambda_{i+1}^{k+1} \\ \\ P_{i+1}^{k+1} &= G_0^t(Q^*(V_{i+1}^k),\lambda_{i+1}^k)\lambda_{i+1}^{k+1}\end{split}$

#### Lagrangian Linear Relations¶

$\begin{split}y &= HQ + D\lambda + FZ + b \\ \dot y &= HV + D\lambda \\ P &= H^t\lambda\end{split}$

The discretisation is straightforward:

$\begin{split}\dot y_{i+1} &= HV_{i+1} + D\lambda_{i+1} \\ P_{i+1} &= H^t\lambda_{i+1}\end{split}$

### Time discretization of the Non Smooth laws¶

A natural way of discretizing the unilateral constraint leads to the following implicit discretization :

$0 \leq y_{i+1} \perp \lambda_{i+1} \geq 0$

In the Moreau’s time–stepping, we use a reformulation of the unilateral constraints in terms of velocity:

$If y(t) =0, \ then \ 0 \leq \dot y \perp \lambda \geq 0$

which leads to the following discretisation :

$If \ y^{p} \leq 0, \ then \ 0 \leq \dot y_{i+1} \perp \lambda_{i+1} \geq 0$

where $$y^{p}$$ is a prediction of the position at time $$t_{i+1}$$, for instance, $$y^{p} = y_{i} + \frac{h}{2} \dot y_i$$.

To introduce a Newton impact law, consider an equivalent velocity defined by

$\dot y^{e}_{i+1} = \dot y_{i+1} + e \dot y_{i}$

and apply the constraints directly on this velocity :

$If \ y^{p} \leq 0, \ then \ 0 \leq \dot y^{e}_{i+1} \perp \lambda_{i+1} \geq 0$

## Summary of the time discretized equations¶

### First order systems¶

• Non Linear dynamics:
$\begin{split}x_{i+1}^{k+1} &= x^{free,k}_{i+1} + h(W_{i+1}^k)^{-1}r_{i+1}^{k+1} \\ W_{i+1}^k &= M - h \theta\cdot\nabla_{x}f(x_{i+1}^k,t_{i+1}) \\ x^{free,k}_{i+1} &= x_{i+1}^k - (W_{i+1}^k)^{-1}\mathcal R^{free}(x_{i+1}^{k}) \\ \mathcal R^{free}(x_{i+1}^{k}) &= M(x_{i+1}^k-x_i) - h \theta f(x_{i+1}^k,t_{i+1}) - h (1-\theta) f(x_{i},t_i)\end{split}$
• Linear dynamics:
$\begin{split}x_{i+1} &= x^{free}_{i+1} + hW_{i+1}^{-1}r_{i+1} \\ W_{i+1} &= (M - h\theta A_{i+1}) \\ x^{free}_{i+1} &= W_{i+1}^{-1}\left[(M + h (1-\theta)A_{i})\cdot x_i + h\theta(b_{i+1}-b_i) + hb_i\right]\end{split}$
• Linear dynamics with time-invariant coefficients:
$\begin{split}x_{i+1} &= x^{free}_{i} + hW^{-1}r_{i+1} \\ W &= (M - h\theta A) \\ x^{free}_i &= x_i + h W^{-1}(A x_i + b)\end{split}$
• Non Linear Relations
$\begin{split}y_{i+1}^{k+1} &= y_{i+1}^k - H_0(S_{i+1}^k)X_{i+1}^{k} - H_1(S_{i+1}^k)\lambda_{i+1}^{k} + H_0(S_{i+1}^k)X_{i+1}^{k+1} + H_1(S_{i+1}^k)\lambda_{i+1}^{k+1} \\ \\ R_{i+1}^{k+1} &= R_{i+1}^k - G_0(S_{i+1}^k)X_{i+1}^{k} - G_1(S_{i+1}^k)\lambda_{i+1}^{k} + G_0(S_{i+1}^k)X_{i+1}^{k+1} + G_1(S_{i+1}^k)\lambda_{i+1}^{k+1} \\ \\ S_{i+1}^k &\ for \ (X_{i+1}^{k},t_{i+1},\lambda_{i+1}^{k}) \\ \\ H_0(X,t,\lambda)&=\nabla_X h(X,t,\lambda), \ H_1(X,t,\lambda)=\nabla_{\lambda} h(X,t,\lambda) \\ \\ G_0(X,t,\lambda)&=\nabla_X g(X,t,\lambda), \ G_1(X,t,\lambda)=\nabla_{\lambda} g(X,t,\lambda) \\\end{split}$
• Linear Relations
$\begin{split}y_{i+1} &= C(t_{i+1})X_{i+1} + D(t_{i+1})\lambda_{i+1} + e(t_{i+1}) + F(t_{i+1})Z \\ R_{i+1} &= B(t_{i+1})\lambda_{i+1}\end{split}$

### Lagrangian second-order systems¶

• Non Linear Dynamics:
$\begin{split}v_{i+1}^{k+1} &= v^{free,k}_{i+1} + (W_{i+1}^k)^{-1} p_{i+1}^{k+1} \\ q_{i+1}^{k+1} &= q_i + h\theta v_{i+1}^{k+1} + h(1 - \theta)v_{i} \\ v^{free,k}_{i+1} &= v_{i+1}^k - (W_{i+1}^k)^{-1} \mathcal R^{free}(v_{i+1}^k) \\ \mathcal R^{free}(v_{i+1}^k) &= M(q*)(v_{i+1}^k-v_{i}) - h\theta f_L(t_{i+1},v_{i+1}^k,q_{i+1}^k) - h(1-\theta) f_L(t_i,v_i,q_i) \\ W_{i+1}^k &= M(q*) + h\theta C_t(t_{i+1},v_{i+1}^k,q_{i+1}^k) + h^2\theta^2 K_t(t_{i+1},v_{i+1}^k,q_{i+1}^k) \\ q^* &= q_i + h\gamma v_i \\ C_t(t,v,q)&=-\left[\frac{\partial{f_L(t,v,q)}}{\partial{v}}\right] \\ K_t(t,v,q)&=-\left[\frac{\partial{f_L(t,v,q)}}{\partial{q}}\right]\end{split}$
• Linear Dynamics with and Time–Invariant Coefficients
$\begin{split}v_{i+1} &= v^{free}_i + W^{-1} p_{i+1} \\ q_{i+1} &= q_{i} + h\left[\theta v_{i+1}+(1-\theta) v_{i} \right]\\ v^{free}_i &= v_{i} + W^{-1}\left[(- hC - h^2\theta K )v_{i} - h K q_{i}+ h\theta F_{ext}(t_{i+1})+h(1-\theta) F_{ext}(t_{i}) \right] \\ W &= \left[M + h\theta C + h^2 \theta^2 K \right]\end{split}$
• Lagrangian Scleronomous Relations
$\begin{split}\dot y_{i+1}^{k+1} = G_0(Q^*(V_{i+1}^k))V_{i+1}^{k+1} \\ P_{i+1}^{k+1} = G_0^t(Q^*(V_{i+1}^k))\lambda_{i+1}^{k+1}\end{split}$
• Lagrangian Rheonomous Relations
$\begin{split}\dot y_{i+1}^{k+1} &= G_0(Q^*(V_{i+1}^k),t_{i+1})V_{i+1}^{k+1} + G_1(Q^*(V_{i+1}^k, t_{i+1})) \\ P_{i+1}^{k+1} &= G_0^t(Q^*(V_{i+1}^k),t_{i+1})\lambda_{i+1}^{k+1}\end{split}$
• Lagrangian Compliant Relations
$\begin{split}\dot y_{i+1}^{k+1} &= G_0(Q^*(V_{i+1}^k),\lambda_{i+1}^k)V_{i+1}^{k+1} + G_1(Q^*(V_{i+1}^k, \lambda_{i+1}^k))\lambda_{i+1}^{k+1} \\ P_{i+1}^{k+1} &= G_0^t(Q^*(V_{i+1}^k),\lambda_{i+1}^k)\lambda_{i+1}^{k+1}\end{split}$
• Lagrangian Linear Relations
$\begin{split}\dot y_{i+1} &= HV_{i+1} + D\lambda_{i+1} \\ P_{i+1} &= H^t\lambda_{i+1}\end{split}$