Implicit Euler

30/12/2020 updated on 19/11/2021  │ Notes
numerical analysis
timestepping

Notations

Let be \(\mathbf{q}\) be a field function of time \(t\) and \(\mathbf{\dot q}, \mathbf{\ddot q}\) be its first and second derivatives with respect to time, \[ \begin{aligned} \mathbf{q} &= \mathbf{q}(t)\\ \mathbf{\dot q} &= \frac{\mathrm{d} \mathbf{q}}{\mathrm{d} t}\\ \mathbf{\ddot q} &= \frac{\mathrm{d} \mathbf{\dot q}}{\mathrm{d} t}\\ \end{aligned} \] In order to make the description of the time integration algorithms less abstract, one can think of $,, $ as position, velocity and acceleration. Moreover, the balance of momentum equation is considered for a concrete derivation of the scheme. Therefore, \[ \notag \mathbf M \mathbf{\ddot q} +\mathbf{f} = 0 \quad \text{where} \quad \mathbf{f} = \mathbf{f}\left(\mathbf{q}(t),\mathbf{\dot q}(t)\right) \] The mass matrix \(\mathbf{M}\) is assumed constant. In addition, to facilitate the notations, the superscripts involving \(t\) will be used to denote time points such that \(\mathbf{q}^{t+\tau} = \mathbf{q}(t+\tau)\) represents the discrete value of \(\mathbf{q}\) at time \(t+\tau\), while subscript involving \(n\) as in \(\mathbf{q}_n^t\) will represent the value of \(\mathbf{q}^t\) at some iteration \(n\). The symbol \(\tau\) represents the timestep.

Implicit Euler

The implicit Euler method is the first of the Backward Differentiation Formula. The scheme is implicit and therefore unconditionally stable, even more so, it is L-stable, which means that the new computed value approaches toward zero in a single step as the timestep tends to infinity. It is the most common integrator for near real-time applications. While it allows to take relatively large timesteps, the computed solution presents severe numerical damping, which reduces its accuracy. The scheme is given as \[ \left\{\begin{aligned} \frac{\mathbf{q}^{t+\tau}-\mathbf{q}^t}{\tau} &= \mathbf{\dot q}^{t+\tau}\\ \frac{\mathbf{\dot q}^{t+\tau}-\mathbf{\dot q}^t}{\tau} &= \mathbf{\ddot q}^{t+\tau}\\ \mathbf{M}\mathbf{\ddot q}^{t+\tau} + \mathbf{f}^{t+\tau}&=0 \end{aligned} \right. \]

Since the value of \(\mathbf{\ddot q}\) is not generally needed, it can be folded in the \(\mathbf{\dot q}\) term \[ \notag \begin{eqnarray*} \left\{\begin{aligned} \mathbf{q}^{t+\tau}&=\mathbf{q}^t+ \tau\mathbf{\dot q}^{t+\tau}\\ \mathbf{M}\mathbf{\dot q}^{t+\tau} &= \mathbf{M}\mathbf{\dot q}^t - \tau\mathbf{f}^{t+\tau} \end{aligned} \right. \end{eqnarray*} \]

With the residual form \[ \begin{eqnarray*} \begin{aligned} \mathrm{R}=\begin{bmatrix}\mathrm{R}_q\\ \mathrm{R}_\dot q\end{bmatrix}= \begin{bmatrix} &\mathbf{q}^{t+\tau}-\mathbf{q}^t - \tau\mathbf{\dot q}^{t+\tau} \\ &\mathbf{M}\left(\mathbf{\dot q}^{t+\tau} - \mathbf{\dot q}^t\right) + \tau\mathbf{f}^{t+\tau} \\ \end{bmatrix} \end{aligned}=0 \end{eqnarray*} \] Now, if we replace \(\mathbf{q}^{t+\tau}\) by \(\mathbf{q}^{t+\tau}_n\!\!\!+\!\Delta\mathbf{q}\) where \(\mathbf{q}^{t+\tau}_n\) is an approximation of \(\mathbf{q}^{t+\tau}\) which we assume known and \(\Delta\mathbf{q}\) a little nudge needed to reach the solution (and do similarly to its time derivatives), we have \[ \left\{\begin{aligned} \mathbf{q}^{t+\tau}_n+\Delta\mathbf{q} &= \mathbf{q}^t+\tau\left(\mathbf{\dot q}^{t+\tau}_n+\Delta\mathbf{\dot q}\right)\\ \mathbf{M}\left( \mathbf{\dot q}^{t+\tau}_n+ \Delta\mathbf{\dot q}\right) &= \mathbf{M}\mathbf{\dot q}^t -\tau\mathbf{f}^{t+\tau}_{n+1} \end{aligned}\right. \quad\text{where} \quad \mathbf{f}^{t+\tau}_{n+1} = \mathbf{f}^{t+\tau}_n\!(\mathbf{q}^{t+\tau}_n\!\!+\!\Delta\mathbf{q},\mathbf{\dot q}^{t+\tau}_n\!\!+\!\Delta\mathbf{\dot q}) \] Taylor approximation of the force term \(\mathbf{f}^{t+\tau}_{n+1}\) yields \[ \begin{aligned} \mathbf{f}^{t+\tau}_{n+1} &\approx \mathbf{f}^{t+\tau}_{n} +\frac{\partial\mathbf{f}^{t+\tau}_n}{\partial\mathbf{q}}\Delta\mathbf{q}+\frac{\partial \mathbf{f}^{t+\tau}_n}{\partial\mathbf{\dot q}}\Delta\mathbf{\dot q}\\ &\approx \mathbf{f}^{t+\tau}_{n} + \mathbf K_n\Delta\mathbf{q}+\mathbf D_n\Delta\mathbf{\dot q} \end{aligned} \] Where \(\mathbf K, \mathbf D\) are recognized as the tangent stiffness and damping matrix, respectively. \[ \notag \mathbf{M}\Delta\mathbf{\dot{q}} = \mathbf{M}(\mathbf{\dot{q}}^t-\mathbf{\dot{q}}^{t+\tau}_n) -\tau\mathbf{f}^{t+\tau}_n -\tau(\mathbf{K}_n\Delta\mathbf{q}+\mathbf{D}_n\Delta\mathbf{\dot{q}}) \] Substituting \(\Delta\mathbf{q}\) yields \[ \begin{aligned} \left(\mathbf{M}+\tau\mathbf{D}_n+\tau^2\mathbf{K}_n\right)\Delta\mathbf{\dot{q}}&= \mathbf{M}\left(\mathbf{\dot q}^t-\mathbf{\dot q}^{t+\tau}_n\right)-\tau\mathbf{f}^{t+\tau}_n-\tau\mathbf{K}\left(\mathbf{q}^t-\mathbf{q}^{t+\tau}_n+\tau\mathbf{\dot q}^{t+\tau}_n \right) \\ &= -\mathrm{R}_{\dot{q},n} +\tau\mathbf{K}\mathrm{R}_{q,n} \end{aligned} \]

The overall procedure consists in repeating the cycle below, given \(\mathbf{q}^{t+\tau}_0\) and \(\mathbf{\dot{q}}^{t+\tau}_0\), until \(\Vert\Delta\mathbf{\dot{q}}\Vert<\text{threshold}\) \[ \begin{equation} \begin{aligned} \left(\mathbf{M}+\tau\mathbf{D}_n+\tau^2\mathbf{K}_n\right)\Delta\mathbf{\dot{q}}&= -\mathrm{R}_{\dot{q},n} +\tau\mathbf{K}\mathrm{R}_{q,n} \\ \mathbf{\dot{q}}^{t+\tau}_{n+1}&=\mathbf{\dot{q}}^{t+\tau}_n+\Delta\mathbf{\dot{q}} \\ \mathbf{q}^{t+\tau}_{n+1} &= \mathbf{q}^t+ \tau\mathbf{\dot q}^{t+\tau}_{n+1} \\ \mathbf{f}^{t+\tau}_{n+1} &= \mathbf{f}\left(\mathbf{q}^{t+\tau}_{n+1},\mathbf{\dot q}^{t+\tau}_{n+1}\right) \end{aligned} \label{eq:iteration_unrolled} \end{equation} \] It is possible to reach the same the solution by applying Newton’s method \[ \notag \mathrm{R}_n+\frac{\mathrm{d}\mathrm R_n}{\mathrm{d}\{\mathbf q^{t+\tau},\mathbf{\dot q}^{t+\tau}\}}\begin{bmatrix}\mathbf{q}^{t+\tau}_{n+1}-\mathbf{q}^{t+\tau}_{n}\\ \mathbf{\dot q}^{t+\tau}_{n+1}-\mathbf{\dot q}^{t+\tau}_{n}\end{bmatrix} = 0 \]

Which yields \[ \notag \begin{bmatrix} \mathbf{I} & -\tau\mathbf{I} \\ \tau\mathbf{K}_n & \mathbf{M}+\tau\mathbf{D}_n \end{bmatrix} \begin{pmatrix} \Delta\mathbf{q}\\ \Delta\mathbf{\dot{q}}\\ \end{pmatrix} = - \begin{pmatrix} \mathrm{R}_{q,n}\\ \mathrm{R}_{\dot{q},n}\\ \end{pmatrix} \]

Remarks

The unrolled procedure \(\eqref{eq:iteration_unrolled}\) effectively projects \(\mathbf{q}\) to fulfill the residual equation at the end of each iteration. Indeed, \(\mathbf{q}^{t+\tau}_{n+1}\) is explicitly computed to fulfill the residual equation \(\mathrm R_{q} =0\). Therefore, its contribution to the main step is null and can be removed, for \(n\geq1\). \[ \begin{aligned} \left(\mathbf{M}+\tau\mathbf{D}_n+\tau^2\mathbf{K}_n\right)\Delta\mathbf{\dot{q}}&= -\mathrm{R}_{\dot{q},n} \\ \mathbf{\dot{q}}^{t+\tau}_{n+1}&=\mathbf{\dot{q}}^{t+\tau}_n+\Delta\mathbf{\dot{q}} \\ \mathbf{q}^{t+\tau}_{n+1} &= \mathbf{q}^t+ \tau\mathbf{\dot q}^{t+\tau}_{n+1} \end{aligned} \] ⚠️ This is not necessarily true for the first iteration \(n=0\). Indeed, if \(\mathbf{q}^{t+\tau}_0\) is chosen independently from \(\mathbf{\dot{q}}^{t+\tau}_0\) then \[ \mathbf{q}^{t+\tau}_0 \neq \mathbf{q}^{t}+\tau\mathbf{\dot{q}}^{t+\tau}_0 \implies \mathrm R_{q,n} \neq 0 \]

It could still be discarded to simplify the algorithm, at the risk of a slower convergence. This might be particularly interesting in non-linear scenarios, where multiple iterations would anyway be needed to converge. However, in a truncated scheme, halted before convergence after a single or few iterations like in near real-time applications, such approximation might be questionable. Nevertheless, it is commonly omitted.

⚠️ The force term \(\mathbf{f}^{t+\tau}_n\) is computed as \(\mathbf{f}(\mathbf{q}_n^{t+\tau},\mathbf{\dot q}_n^{t+\tau})\). However, at the first iteration, \(\mathbf {f}^{t+\tau}_0\) could also chosen independently from the choices of \(\mathbf{q}^{t+\tau}_0\) and \(\mathbf{\dot{q}}^{t+\tau}_0\) in order to accelerate the convergence of the scheme or simplify its formulation. \[ \mathbf{f}^{t+\tau}_0 \neq \mathbf{f}(\mathbf{q}^{t+\tau}_0,\mathbf{\dot{q}}^{t+\tau}_0) \]

Implicit Euler linearized

A common choice of initial guesses are \(\mathbf{q}^{t+\tau}_0= \mathbf{q}^t\), \(\mathbf{\dot q}^{t+\tau}_0=0\) and \(\mathbf{f}^{t+\tau}_0=\mathbf{f}^t\). Assuming that the residual on position is omitted despite such inconsistent input, it yields for the first Newton iteration \[ \begin{aligned} \left(\mathbf{M}+\tau\mathbf{D}_0+\tau^2\mathbf{K}_0\right)\mathbf{\dot{q}}^{t+\tau}_1&= \mathbf{M}\mathbf{q}^t - \tau\mathbf{f}^t \\ \mathbf{q}^{t+\tau}_{1} &= \mathbf{q}^t+ \tau\mathbf{\dot q}^{t+\tau}_{1} \end{aligned} \] Stopping at the first iteration gives the implicit Euler linearized scheme, an integrator that can handle complex phenomena, such as impacts and discontinuities, while being numerically tractable for soft real-time systems.

Further initial guesses

With \(\mathbf{q}^{t+\tau}_0= \mathbf{q}^t\) and \(\mathbf{\dot q}^{t+\tau}_0=\mathbf{\dot q}^{t}\) \[ \begin{aligned} \left(\mathbf{M}+\tau\mathbf{D}+\tau^2\mathbf{K}\right)\mathbf{\dot{q}}^{t+\tau}_1&= (\mathbf{M}+\tau\mathbf{D})\mathbf{\dot{q}}^t - \tau\mathbf{f}^t \end{aligned} \]

With \(\mathbf{q}^{t+\tau}_0= \mathbf{q}^t + \tau\mathbf{\dot q}^t\) and \(\mathbf{\dot q}^{t+\tau}_0=\mathbf{\dot q}^{t}\) \[ \begin{aligned} \left(\mathbf{M}+\tau\mathbf{D}+\tau^2\mathbf{K}\right)\mathbf{\dot{q}}^{t+\tau}_1&= \left(\mathbf{M}+\tau\mathbf{D}+\tau^2\mathbf{K}\right)\mathbf{\dot{q}}^t - \tau\mathbf{f}^t \end{aligned} \]

With \(\mathbf{q}^{t+\tau}_0= \mathbf{q}^t + \tau\mathbf{\dot q}^t\) and \(\mathbf{\dot q}^{t+\tau}_0=0\) \[ \begin{aligned} \left(\mathbf{M}+\tau\mathbf{D}+\tau^2\mathbf{K}\right)\mathbf{\dot{q}}^{t+\tau}_1&= \left(\mathbf{M}+\tau^2\mathbf{K}\right)\mathbf{\dot{q}}^t - \tau\mathbf{f}^t \end{aligned} \]

With \(\mathbf{q}^{t+\tau}_0= \mathbf{q}^t + \tau\mathbf{\dot q}^t\) and \(\mathbf{\dot q}^{t+\tau}_0=\mathbf{\dot q}^t+\tfrac{\tau}{\tau}(\mathbf{\dot q}^{t}-\mathbf{\dot q}^{t-\tau})=2\mathbf{\dot q}^t-\mathbf{\dot q}^{t-\tau}\) \[ \begin{aligned} \left(\mathbf{M}+\tau\mathbf{D}+\tau^2\mathbf{K}\right)\mathbf{\dot{q}}^{t+\tau}_1&= \left(\mathbf{M}+\tau\mathbf{D}+\tau^2\mathbf{K}\right)\mathbf{\dot{q}}^t +\tau\mathbf{D}(\mathbf{\dot q}^t-\mathbf{\dot q}^{t-\tau}) - \tau\mathbf{f}^t \end{aligned} \]

With \(\mathbf{q}^{t+\tau}_0= \mathbf{q}^t + \tau(2\mathbf{\dot q}^t-\mathbf{\dot q}^{t-\tau})\) and \(\mathbf{\dot q}^{t+\tau}_0=2\mathbf{\dot q}^t-\mathbf{\dot q}^{t-\tau}\)

\[ \begin{aligned} \left(\mathbf{M}+\tau\mathbf{D}+\tau^2\mathbf{K}\right)\mathbf{\dot{q}}^{t+\tau}_1&= \left(\mathbf{M}+2\tau\mathbf{D}+2\tau^2\mathbf{K}\right)\mathbf{\dot{q}}^t -\left(\tau\mathbf{D}-\tau^2\mathbf{K}\right)\mathbf{\dot q}^{t-\tau} - \tau\mathbf{f}^t \end{aligned} \]