Implicit midpoint method

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 \[ \left\{\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}\right. \] In order to make the description of the time integration algorithms less abstract, one can think of \(\mathbf{q},\mathbf{\dot q},\mathbf{\ddot q}\) 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}\left(t\right),\mathbf{\dot q}\left(t\right)\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 midpoint method

The scheme is given as \[ \left\{\begin{aligned} \frac{\mathbf{q}^{t+\tau}-\mathbf{q}^t}{\tau} &= \tfrac{1}{2}\left(\mathbf{\dot q}^{t}+\mathbf{\dot q}^{t+\tau}\right)\\ \frac{\mathbf{\dot q}^{t+\tau}-\mathbf{\dot q}^t}{\tau} &= \tfrac{1}{2}\left(\mathbf{\ddot q}^{t}+\mathbf{\ddot q}^{t+\tau}\right)\\ \mathbf{M}\mathbf{\ddot q}^{t+\tau} + \mathbf{f}^{t+\tau}&=0 \end{aligned} \right. \]

Since the value of \(\mathbf{\ddot q}\) is not often desired, it can be folded in the \(\mathbf{\dot q}\) term \[ \left\{\begin{aligned} \frac{\mathbf{q}^{t+\tau}-\mathbf{q}^t}{\tau} &= \tfrac{1}{2}\left(\mathbf{\dot q}^{t}+\mathbf{\dot q}^{t+\tau}\right)\\ \mathbf{M}\frac{\mathbf{\dot q}^{t+\tau}-\mathbf{\dot q}^t}{\tau} &= -\tfrac{1}{2}\left(\mathbf{f}^{t}+\mathbf{f}^{t+\tau}\right) \end{aligned} \right. \]

One way to solve this problem is to apply Newton’s method. However, I always had issues wrapping my head around the direct application of the Newton’s method, therefore let us go over an unrolled version of the method. First, we reorganize the terms from the scheme for convenience \[ \notag \begin{eqnarray*} \left\{\begin{aligned} \mathbf{q}^{t+\tau}&=\mathbf{q}^t+ \tfrac{\tau}{2}\left(\mathbf{\dot q}^{t}+\mathbf{\dot q}^{t+\tau}\right)\\ \mathbf{M}\mathbf{\dot q}^{t+\tau} &= \mathbf{M}\mathbf{\dot q}^t - \tfrac{\tau}{2}\left(\mathbf{f}^{t}+\mathbf{f}^{t+\tau}\right) \end{aligned} \right. \end{eqnarray*} \] With 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 - \tfrac{\tau}{2}\left(\mathbf{\dot q}^{t}+\mathbf{\dot q}^{t+\tau}\right) \\ &\mathbf{M}\left(\mathbf{\dot q}^{t+\tau} - \mathbf{\dot q}^t\right) - \tfrac{\tau}{2}\left(\mathbf{f}^{t}+\mathbf{f}^{t+\tau}\right) \\ \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+ \tfrac{\tau}{2}\left(\mathbf{\dot q}^{t}+\mathbf{\dot q}^{t+\tau}_n+\Delta\mathbf{\dot q}\right)\\ \mathbf{M}(\mathbf{\dot q}^{t+\tau}_n+\Delta\mathbf{\dot q}) &= \mathbf{M}\mathbf{\dot q}^t - \tfrac{\tau}{2}\left(\mathbf{f}^{t}+\mathbf{f}^{t+\tau}_{n+1}\right) \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) -\tfrac{\tau}{2}(\mathbf{f}^t+\mathbf{f}^{t+\tau}_n) -\tfrac{\tau}{2}(\mathbf{K}_n\Delta\mathbf{q}+\mathbf{D}_n\Delta\mathbf{\dot{q}}) \] Substituting \(\Delta\mathbf{q}\) yields \[ \begin{aligned} \mathbf{M}\Delta\mathbf{\dot{q}} &= \mathbf{M}(\mathbf{\dot{q}}^t-\mathbf{\dot{q}}^{t+\tau}_n) -\tfrac{\tau}{2}(\mathbf{f}^t+\mathbf{f}^{t+\tau}_n) -\tfrac{\tau}{2}\mathbf{K}_n\left(\mathbf{q}^t-\mathbf{q}^{t+\tau}_n+ \tfrac{\tau}{2}\left(\mathbf{\dot q}^{t}+\mathbf{\dot q}^{t+\tau}_n+\Delta\mathbf{\dot q}\right)\right) -\tfrac{\tau}{2}\mathbf{D}_n\Delta\mathbf{\dot{q}} \\ &= \mathbf{M}(\mathbf{\dot{q}}^t-\mathbf{\dot{q}}^{t+\tau}_n) -\tfrac{\tau}{2}(\mathbf{f}^t+\mathbf{f}^{t+\tau}_n) -\tfrac{\tau}{2}\mathbf{K}_n\left(\mathbf{q}^t-\mathbf{q}^{t+\tau}_n+ \tfrac{\tau}{2}\left(\mathbf{\dot q}^{t}+\mathbf{\dot q}^{t+\tau}_n\right)\right) -\tfrac{\tau^2}{4}\mathbf{K}_n\Delta\mathbf{\dot q} -\tfrac{\tau}{2}\mathbf{D}_n\Delta\mathbf{\dot{q}} \\ &= -\mathrm{R}_{\dot{q},n} +\tfrac{\tau}{2}\mathbf{K}_n\mathrm{R}_{q,n} -\tfrac{\tau^2}{4}\mathbf{K}_n\Delta\mathbf{\dot q} -\tfrac{\tau}{2}\mathbf{D}_n\Delta\mathbf{\dot{q}} \\ \left(\mathbf{M}+\tfrac{\tau}{2}\mathbf{D}_n+\tfrac{\tau^2}{4}\mathbf{K}_n\right)\Delta\mathbf{\dot{q}}&= -\mathrm{R}_{\dot{q},n} +\tfrac{\tau}{2}\mathbf{K}_n\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}+\tfrac{\tau}{2}\mathbf{D}_n+\tfrac{\tau^2}{4}\mathbf{K}_n\right)\Delta\mathbf{\dot{q}}&= -\mathrm{R}_{\dot{q},n} +\tfrac{\tau}{2}\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+ \tfrac{\tau}{2}\left(\mathbf{\dot q}^{t}+\mathbf{\dot q}^{t+\tau}_{n+1}\right) \\ \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} & -\tfrac{\tau}{2}\mathbf{I} \\ \tfrac{\tau}{2}\mathbf{K}_n & \mathbf{M}+\tfrac{\tau}{2}\mathbf{D}_n \end{bmatrix} \begin{pmatrix} \vphantom{\tfrac{\tau}{2}}\Delta\mathbf{q}\\ \vphantom{\tfrac{\tau}{2}}\Delta\mathbf{\dot{q}}\\ \end{pmatrix} = - \begin{pmatrix} \vphantom{\tfrac{\tau}{2}}\mathrm{R}_{q,n}\\ \vphantom{\tfrac{\tau}{2}}\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}+\tfrac{\tau}{2}\mathbf{D}_n+\tfrac{\tau^2}{4}\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+ \tfrac{\tau}{2}\left(\mathbf{\dot q}^{t}+\mathbf{\dot q}^{t+\tau}_{n+1}\right) \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}+\tfrac{\tau}{2}\left(\mathbf{\dot{q}}^{t+\tau}_0 + \mathbf{\dot{q}}^{t}\right) \implies \mathrm R_{q,n} \neq 0 \]

It could still be omitted 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 reduce the accuracy of the solution¹.

⚠️ 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) \]

Some choices of initial guesses

\(\mathbf{\dot q}^{t+\tau}_0 = 0,\quad \mathbf{q}^{t+\tau}_0 = \mathbf{q}^{t},\quad \mathbf{f}^{t+\tau}_0 = \mathbf{f}^{t}\)

🚧

\(\mathbf{\dot q}^{t+\tau}_0 = \mathbf{\dot q}^{t}_0, \quad \mathbf{q}^{t+\tau}_0 = \mathbf{q}^{t}, \quad \mathbf{f}^{t+\tau}_0 = \mathbf{f}^{t}\)

🚧