Newmark Methods#

We are looking to approximate second order differential equations of the form

\[ \mathbf{M} \ddot{\boldsymbol{u}} + \boldsymbol{p}(\boldsymbol{u}, \dot{\boldsymbol{u}}) = \boldsymbol{f} \]

Newmark Equations#

The Newmark-$\beta$ scheme is defined by the following equations:

\[\begin{split} \boxed{ \begin{aligned} &\boldsymbol{v}_{n+1}= \boldsymbol{v}_n + \Delta t~ \left( (1 - \gamma) \boldsymbol{a}_n + \gamma \boldsymbol{a}_{n+1} \right)\\ &\boldsymbol{d}_{n+1}= \boldsymbol{d}_n + \Delta t~\boldsymbol{v}_n + \tfrac{1}{2}\Delta t^2 \left(\left(1 - 2\beta\right) \boldsymbol{a}_n + 2\beta \boldsymbol{a}_{n+1}\right)\\ % &m \boldsymbol{a}_{n+1} + c \boldsymbol{v}_{n+1} + p(\boldsymbol{d}_{n+1}) = \boldsymbol{f}_{n+1}^{\textrm{ext}} \end{aligned} } \end{split}\]

for scalar parameters $\beta$ and $\gamma$ and the following notation is used:

\[\begin{split} \begin{aligned} \boldsymbol{d}_n &\approx \boldsymbol{u}(t_n) \\ \boldsymbol{v}_n &\approx \dot{\boldsymbol{u}}(t_n) \\ \boldsymbol{a}_n &\approx \ddot{\boldsymbol{u}}(t_n) \\ \end{aligned} \end{split}\]

Remarks#

$\forall \gamma$ and $\forall \beta$, the following holds true $$ \begin{gathered} \int_{t_n}^{t_{n+1}} \ddot{\boldsymbol{u}}(\tau) d \tau=(1-\gamma) h \ddot{\boldsymbol{u}}_n+\gamma h \ddot{\boldsymbol{u}}_{n+1}+\mathbf{r}_n \\ \int_{t_n}^{t_{n+1}} \ddot{\boldsymbol{u}}(\tau)\left(t_{n+1}-\tau\right) d \tau=\left(\frac{1}{2}-\beta\right) h^2 \ddot{\boldsymbol{u}}_n+\beta h^2 \ddot{\boldsymbol{u}}_{n+1}+\mathbf{r}_n^{\prime} \end{gathered} $$ where $$ \mathbf{r}_n=\left(\frac{1}{2}-\gamma\right) h^2 \boldsymbol{u}^{(3)}(\tilde{\tau})+\mathcal{O}\left(h^3 \boldsymbol{u}^{(4)}\right) ; \quad \mathbf{r}_n^{\prime}=\left(\frac{1}{6}-\beta\right) h^3 \boldsymbol{u}^{(3)}(\tilde{\tau})+\mathcal{O}\left(h^4 \boldsymbol{u}^{(4)}\right) $$ and $t_n<\tilde{\tau}<t_{n+1}$

The Newmark-$\beta$ follows from Taylor’s theorem with explicit remainder which states that the first time derivative can be expressed as:

\[\begin{split} \begin{aligned} \dot{\boldsymbol{u}}_{n+1} &= \dot{\boldsymbol{u}}_n + \Delta t~ \ddot{\boldsymbol{u}}_\gamma \\ \boldsymbol{u}_{n+1} &= \boldsymbol{u}_n + \Delta t~ \dot{\boldsymbol{u}}_n + \begin{matrix} \frac 1 2 \end{matrix} \Delta t^2~ \ddot{\boldsymbol{u}}_\beta. \end{aligned} \end{split}\]

for some $\ddot{\boldsymbol{u}}_\gamma$ and $\ddot{\boldsymbol{u}}_\beta$ with the form:

\[\begin{split} \begin{aligned} \ddot{\boldsymbol{u}}_\gamma &= (1 - \gamma) \ddot{\boldsymbol{u}}_n + \gamma \ddot{\boldsymbol{u}}_{n+1}~~~~~~~~0\leq \gamma \leq 1\\ \ddot{\boldsymbol{u}}_\beta &= (1 - 2\beta) \ddot{\boldsymbol{u}}_n + 2\beta \ddot{\boldsymbol{u}}_{n+1}~~~~0\leq 2\beta\leq 1 \end{aligned} \end{split}\]
$\gamma=\frac{1}{2}$ and $\beta=\frac{1}{6}$ corresponds to linearly interpolating $\ddot{\boldsymbol{u}}$ in $\left[t_n, t_{n+1}\right]$ $$ \overline{\ddot{\boldsymbol{u}}}_{\text{lin}}(\tau)=\ddot{\boldsymbol{u}}_n+\left(\tau-t_n\right)\left(\frac{\ddot{\boldsymbol{u}}_{n+1}-\ddot{\boldsymbol{u}}_n}{h}\right) $$
$\gamma=\frac{1}{2}$ and $\beta=\frac{1}{4}$ corresponds to averaging $\ddot{\boldsymbol{u}}$ over $\left[t_n, t_{n+1}\right]$: $$ \ddot{\boldsymbol{u}}_{\text{avg}}=\frac{\ddot{\boldsymbol{u}}_{n+1}+\ddot{\boldsymbol{u}}_n}{2} $$ This is known as the constant average acceleration, or mid-point rule
$\gamma=0.5$ and $\beta=0$ corresponds to an explicit central difference scheme

\[ \frac{1}{\Delta t}\left(v_{n+1} - v_n\right) = (1-\gamma) a_n + \gamma a_{n+1} \]

Generalized $\alpha$ Scheme#

The generalized $\alpha$ scheme imposes equilibrium with the following form:

\[ \begin{aligned} \mathbf{M}\, \boldsymbol{a}_{\alpha_{m}} + \boldsymbol{p}(\boldsymbol{d}_{\alpha_{f}}, \boldsymbol{v}_{\alpha_{f}}) &= \boldsymbol{f}(t_{\alpha_{f}}) \end{aligned} \]

where

\[\begin{split}\begin{aligned} t_{\alpha_{f}} &\triangleq (1-\alpha_{f}) \, t_{n+1} + \alpha_{f} \, t_n \\ \boldsymbol{d}_{\alpha_{f}} &\triangleq (1-\alpha_{f}) \, \boldsymbol{d}_{n+1} + \alpha_{f} \, \boldsymbol{d}_n \\ \boldsymbol{v}_{\alpha_{f}} &\triangleq (1-\alpha_{f}) \, \boldsymbol{v}_{n+1} + \alpha_{f} \, \boldsymbol{v}_n \\ \boldsymbol{a}_{\alpha_{m}} &\triangleq (1-\alpha_{m}) \, \boldsymbol{a}_{n+1} + \alpha_{m} \, \boldsymbol{a}_n \end{aligned}\end{split}\]

However, the following is also sometimes used: $$ \mathbf{M} \ddot{\boldsymbol{u}}_{n+1}+(1-\alpha) \boldsymbol{p}\left(\boldsymbol{u}_{n+1}, \dot{\boldsymbol{u}}_{n+1}\right)+\alpha \boldsymbol{p}\left(\boldsymbol{u}_n, \dot{\boldsymbol{u}}_n\right)=(1-\alpha) \boldsymbol{f}\left(\boldsymbol{u}_{n+1}, t\right)+\alpha \boldsymbol{f}\left(\boldsymbol{u}_n, t\right) $$

Implementation#

The Newmark equations allow one to express all three unknowns (i.e., $a_{n+1}, v_{n+1}$ or $d_{n+1}$) in terms of one of these.

The most straight-forward way to implement the Newmark scheme, is to form a single equation in terms of unknown accelerations, $\boldsymbol{a}_{n+1}$:

\[ \mathbf{A} \boldsymbol{a}_{n+1} = \boldsymbol{b} \]

for some known matrix $\mathbf{A}$ and vector $\boldsymbol{b}$. This is obtained by plugging the Newmark equations into the discrete equilibrium statement at $t_{n+1}$. However before deriving this, it is useful rewrite equations [@eq:newmark] as follows:

\[\begin{split} \begin{aligned} \boldsymbol{d}_{n+1} &= \tilde{\boldsymbol{d}} + c_{da} \, \boldsymbol{a}_{n+1} \\ \boldsymbol{v}_{n+1} &= \tilde{\boldsymbol{v}} + c_{va} \, \boldsymbol{a}_{n+1} \end{aligned} \end{split}\]

where the tilde variables allow us to collect everything that is known at the start of a time step. Expanding and plugging into the discretized equilibrium equation furnishes the Newmark scheme in “a-form”:

\[\begin{split}\begin{aligned} \tilde{\boldsymbol{d}} &= \boldsymbol{d}_n + \Delta t \boldsymbol{v}_n + \tfrac{1}{2} \Delta t^2 (1-2\beta) \boldsymbol{a}_n \\ \tilde{\boldsymbol{v}} &= \boldsymbol{v}_n + \Delta t (1-\gamma) \boldsymbol{a}_n \\ (\mathbf{M} + \gamma\Delta t\, \mathbf{C}+\beta \Delta t^2\, \mathbf{K})\boldsymbol{a}_{n+1} &= \boldsymbol{f}_{n+1} - \mathbf{C}\tilde{\boldsymbol{v}} - \mathbf{K}\tilde{\boldsymbol{d}} \\ \end{aligned}\end{split}\]

\[\begin{split}\begin{aligned} \boldsymbol{d}_{n+1} &= \tilde{\boldsymbol{d}} + \beta \Delta t^2 \, \boldsymbol{a}_{n+1} \\ \boldsymbol{v}_{n+1} &= \tilde{\boldsymbol{v}} + \gamma \Delta t \, \boldsymbol{a}_{n+1} \end{aligned}\end{split}\]

Generalized Unknowns#

Alternatively, the Newmark equations can be manipulated to produce a problem in terms of velocity or displacement.

If $\boldsymbol{x}_{n+1}$ denotes our chosen unknown we can write the equations in the following form:

\[\begin{split} \begin{aligned} \boldsymbol{d}_{\alpha} &= \tilde{\boldsymbol{d}}_x + c_{dx} \, \boldsymbol{x}_{n+1} \\ \boldsymbol{v}_{\alpha} &= \tilde{\boldsymbol{v}}_x + c_{vx} \, \boldsymbol{x}_{n+1} \\ \boldsymbol{a}_{\alpha} &= \tilde{\boldsymbol{a}}_x + c_{ax} \, \boldsymbol{x}_{n+1} \\ \end{aligned} \end{split}\]

This generalizes as follows :

\[\begin{split}\begin{aligned} \tilde{\boldsymbol{d}} &= \sum_i b_{di}\boldsymbol{y}_i \\ \tilde{\boldsymbol{v}} &= \sum_i b_{vi}\boldsymbol{y}_i \\ \tilde{\boldsymbol{a}} &= \sum_i b_{ai}\boldsymbol{y}_i \\ (c_a\,\mathbf{M} + c_v\, \mathbf{C} + c_d\, \mathbf{K})\boldsymbol{x}_{n+1} &= \boldsymbol{f}_{n+1} - \left(\mathbf{M}\tilde{\boldsymbol{a}} + \mathbf{C}\tilde{\boldsymbol{v}} + \mathbf{K}\tilde{\boldsymbol{d}} \right) \\ \end{aligned}\end{split}\]

\[\begin{split}\begin{aligned} \boldsymbol{d}_{n+1} &= \tilde{\boldsymbol{d}} + c_d \, \boldsymbol{x}_{n+1} \\ \boldsymbol{v}_{n+1} &= \tilde{\boldsymbol{v}} + c_v \, \boldsymbol{x}_{n+1} \\ \boldsymbol{a}_{n+1} &= \tilde{\boldsymbol{a}} + c_a \, \boldsymbol{x}_{n+1} \end{aligned}\end{split}\]

where $\boldsymbol{y}_i = (\boldsymbol{d}_n, \boldsymbol{v}_n ,\boldsymbol{a}_n)$ coefficients $c_{yx}$ and $b_{ij}$ given below, and where $\tilde{\left(\cdot\right)}_x$ variables encapsulate the information that is known at the start of the time step.

Nonlinear Generalized - $\alpha$#

Applying this to a nonlinear problem yields:

\[ \begin{aligned} r(\boldsymbol{x}_{n+1})&= \mathbf{M}\, \left[\tilde{\boldsymbol{a}} + c_{ax}\,\boldsymbol{x}_{n+1}\right] + \mathbf{C}\, \left[\tilde{\boldsymbol{v}} + c_{vx}\,\boldsymbol{x}_{n+1}\right] + \boldsymbol{p}\left(\tilde{\boldsymbol{d}} + c_{dx}\,\boldsymbol{x}_{n+1}, \tilde{\boldsymbol{v}} + c_{vx}\,\boldsymbol{x}_{n+1}\right) - \boldsymbol{f}(t_{n+1}) \end{aligned} \]

Linearizing for $r_\eta = r(\boldsymbol{x}_{n+1}^i + \eta \, d\boldsymbol{x})$ yields

\[\begin{split} \begin{aligned} \mathcal{L} \, r_{\eta} &= r(\boldsymbol{x}^i_{n+1}) + \mathbf{A} d \boldsymbol{x} \\ \end{aligned} \end{split}\]

where

\[ \mathbf{A} d\boldsymbol{x} = \left.\frac{d}{d \eta}r_\eta\right|_{\eta=0} = \left( c_d \mathbf{K}\left(\tilde{\boldsymbol{d}}+c_{d}\,\boldsymbol{x}_{n+1}\right) + c_v \mathbf{C} + c_a \mathbf{M}\right) d \boldsymbol{x} \]

\[ A(\boldsymbol{x})= \frac{\partial \boldsymbol{f}}{\partial \boldsymbol{u}}\frac{\partial \boldsymbol{u}}{\partial \boldsymbol{x}} +\frac{\partial \boldsymbol{f}}{\partial \dot{\boldsymbol{u}}} \frac{\partial \dot{\boldsymbol{u}}}{\partial \boldsymbol{x}} +M \frac{\partial \ddot{\boldsymbol{u}}}{\partial \boldsymbol{x}} -\frac{\partial \boldsymbol{p}}{\partial \boldsymbol{u}}\frac{\partial \boldsymbol{u}}{\partial \boldsymbol{x}} \]

\[ \frac{\partial \ddot{\boldsymbol{u}}}{\partial \boldsymbol{u}}=\frac{1}{\beta h^2} \boldsymbol{1} \quad \text { and } \quad \frac{\partial \dot{\boldsymbol{u}}}{\partial \boldsymbol{u}}=\frac{\gamma}{\beta h} \boldsymbol{1} \]

For a chosen unknown $\boldsymbol{x}$, we can collect coefficients into a “predictor matrix” $\mathsf{B}_x$, and corrector array $c_{yx}$ so that any of the other unknowns, $\boldsymbol{y} \in (\boldsymbol{d}, \boldsymbol{v} ,\boldsymbol{a})$ can be expressed as:

\[ \boldsymbol{y}_{n+1} = \mathsf{B}_x \, \boldsymbol{y}_n + {c}_{yx}\, \boldsymbol{x}_{n+1} \]

From the Newmark equations, the coefficients $c_{xy}$ follow as:

\[\begin{split} c_{xy} = \begin{pmatrix} 1 & \frac{\gamma}{\beta \Delta t} & \frac{1}{\beta \Delta t^2} \\ \frac{\Delta t \beta}{\gamma} & 1 & \frac{1}{\gamma \Delta t} \\ \beta \Delta t^2 & \gamma \Delta t & 1 \end{pmatrix} \end{split}\]

The matrix $\mathsf{B}$ is given for the acceleration formulation ($\boldsymbol{x} \triangleq \boldsymbol{a}$) by:

\[\begin{split} \begin{pmatrix} \tilde{ \boldsymbol{d}} \\ \tilde{ \boldsymbol{v}} \\ \tilde{ \boldsymbol{a}} \end{pmatrix}_a = \begin{pmatrix} 1 & \Delta t & \Delta t^2 \left(\tfrac{1}{2} - \beta\right) \\ & 1 & \Delta t \left(1 - \gamma\right) \\ && 0 \end{pmatrix} \begin{pmatrix} \boldsymbol{d}_n \\ \boldsymbol{v}_n \\ \boldsymbol{a}_n \end{pmatrix} \end{split}\]

In the velocity formulation ($\boldsymbol{x} \triangleq \boldsymbol{v}$):

\[\begin{split} \begin{pmatrix} \tilde{ \boldsymbol{d}} \\ \tilde{ \boldsymbol{v}} \\ \tilde{ \boldsymbol{a}} \end{pmatrix}_v = \begin{pmatrix} 1 & - \Delta t\frac{\beta}{\gamma}\left(1 - \frac{\gamma}{\beta}\right) & \Delta t^2 \frac{\beta}{\gamma}\left(\frac{\gamma}{2\beta} - 1 \right) \\ & 0 & \\ &\frac{-1}{\gamma \Delta t} & 1 - \frac{1}{\gamma} \end{pmatrix} \begin{pmatrix} \boldsymbol{d}_n \\ \boldsymbol{v}_n \\ \boldsymbol{a}_n \end{pmatrix} \end{split}\]

Finally, for the displacement formulation ($\boldsymbol{x} \triangleq \boldsymbol{d}$):

\[\begin{split} \begin{pmatrix} \tilde{ \boldsymbol{d}} \\ \tilde{ \boldsymbol{v}} \\ \tilde{ \boldsymbol{a}} \end{pmatrix}_d = \begin{pmatrix} 0 \\ -\frac{\gamma}{\beta\Delta t} & 1 - \frac{\gamma}{\beta} & \Delta t \left(1 - \frac{\gamma}{2\beta}\right) \\ \frac{1}{\beta \Delta t^2} & \frac{-1}{\beta\Delta t} & \quad \left(1 + \tfrac{1}{2\beta}\right) \end{pmatrix} \begin{pmatrix} \boldsymbol{d}_n \\ \boldsymbol{v}_n \\ \boldsymbol{a}_n \end{pmatrix} \end{split}\]

\[\begin{split} \begin{array}{lccccc} \hline \text { Algorithm } & \gamma & \beta & \begin{array}{c} \text { Stability } \\ \text { limit } \\ \omega h \end{array} & \begin{array}{c} \text { Numerical } \\ \text { damping } \\ \text { ratio } \bar{\varepsilon} \end{array} & \begin{array}{c} \text { Periodicity } \\ \text { error } \\ \frac{\Delta T}{T} \end{array} \\ \hline \text { Purely explicit } & 0 & 0 & 0 & \frac{-\omega h}{4} & - \\ \text { Central difference } & \frac{1}{2} & 0 & 2 & 0 & \frac{-\omega^2 h^2}{24} \\ \text { Fox \& Goodwin } & \frac{1}{2} & \frac{1}{12} & 2.45 & 0 & O\left(\omega^4 h^4\right) \\ \begin{array}{l} \text { Linear acceleration } \end{array} & \frac{1}{2} & \frac{1}{6} & 3.46 & 0 & \frac{\omega^2 h^2}{24} \\ \begin{array}{c} \text { Average constant } \\ \text { acceleration } \end{array} & \frac{1}{2} & \frac{1}{4} & \infty & 0 & \frac{\omega^2 h^2}{12} \\ \begin{array}{c} \text { Average constant } \\ \text { acceleration } \\ \text { (with damping) } \end{array} & \frac{1}{2}+\alpha & \frac{(1+\alpha)^2}{4} & \infty & \frac{\alpha}{2} \omega h & \left(\frac{1}{12}+\frac{\alpha^2}{4}\right) \omega^2 h^2 \\ \hline \end{array} \end{split}\]

References#

J. Chung, G.M.Hubert. “A Time Integration Algorithm for Structural Dynamics with Improved Numerical Dissipation: The Generalized-$\alpha$ Method” ASME Journal of Applied Mechanics, 60, 371:375, 1993.