Theta 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, \[ \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.

Wilson-θ method

The Wilson-θ method is defined as \[ \left\{\begin{aligned} \frac{\mathbf{q}^{t+\tau}-\mathbf{q}^t}{\tau} &= (1-\theta)\mathbf{\dot q}^{t}+\theta\mathbf{\dot q}^{t+\tau}&&\theta\in\mathbb [0,1]\\ \frac{\mathbf{\dot q}^{t+\tau}-\mathbf{\dot q}^t}{\tau} &= (1-\theta)\mathbf{\ddot q}^{t}+\theta\mathbf{\ddot q}^{t+\tau}\\ \mathbf{M}\mathbf{\ddot q}^{t+\tau} + \mathbf{f}^{t+\tau}&=0 \end{aligned} \right. \]

One can recognize the explicit Euler ($\theta=0$), the implicit midpoint ($\theta=\tfrac{1}{2}$), and the implicit Euler ($\theta=1$). Since the value of $\mathbf{\ddot q}$ is not generally needed, it can be substituted in the $\mathbf{\dot q}$ term \[ \notag \begin{eqnarray*} \left\{\begin{aligned} \mathbf{q}^{t+\tau}&=\mathbf{q}^t + \tau\left((1-\theta)\mathbf{\dot q}^t+\theta\mathbf{\dot q}^{t+\tau}\right)\\ \mathbf{M}\mathbf{\dot q}^{t+\tau} &= \mathbf{M}\mathbf{\dot q}^t - \tau\left((1-\theta)\mathbf{f}^t+\theta\mathbf{f}^{t+\tau}\right) \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\left((1-\theta)\mathbf{\dot q}^t+\theta\mathbf{\dot q}^{t+\tau}\right) \\ &\mathbf{M}\left(\mathbf{\dot q}^{t+\tau} - \mathbf{\dot q}^t\right) + \tau\left((1-\theta)\mathbf{f}^t+\theta\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 + \tau\left((1-\theta)\mathbf{\dot q}^t+\theta(\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 - \tau\left((1-\theta)\mathbf{f}^t+\theta\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) -\tau\left((1-\theta)\mathbf{f}^t+\theta\mathbf{f}^{t+\tau}_{n}\right) -\theta\tau(\mathbf{K}_n\Delta\mathbf{q}+\mathbf{D}_n\Delta\mathbf{\dot{q}}) \] Substituting $\Delta\mathbf{q}$ yields \[ \begin{aligned} \left(\mathbf{M}+\theta\tau\mathbf{D}_n+\theta^2\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\left((1-\theta)\mathbf{f}^t+\theta\mathbf{f}^{t+\tau}_{n}\right) -\theta\tau\mathbf{K}_n\left(\mathbf{q}^t-\mathbf{q}^{t+\tau}_n+\tau\left((1-\theta)\mathbf{\dot q}^t+\theta\mathbf{\dot q}^{t+\tau}\right) \right) \\ &= -\mathrm{R}_{\dot{q},n} +\theta\tau\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}+\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\left((1-\theta)\mathbf{\dot q}^t+\theta\mathbf{\dot q}^{t+\tau}\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} & -\theta\tau\mathbf{I} \\ \theta\tau\mathbf{K}_n & \mathbf{M}+\theta\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}+\theta\tau\mathbf{D}_n+\theta^2\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 + \Delta\mathbf{\dot q} \\ \mathbf{q}^{t+\tau}_{n+1} &= \mathbf{q}^t + \tau\left((1-\theta)\mathbf{\dot q}^t+\theta\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 \[ \notag \mathbf{q}^{t+\tau}_0 \neq \mathbf{q}^{t}+\tau\left((1-\theta)\mathbf{\dot{q}}^{t}+\theta\mathbf{\dot{q}}^{t+\tau}_0\right) \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. \[ \notag \mathbf{f}^{t+\tau}_0 \neq \mathbf{f}(\mathbf{q}^{t+\tau}_0,\mathbf{\dot{q}}^{t+\tau}_0) \]