Newmark-beta method

<h2 id="stability-analysis">Stability Analysis</h2>
A time-integration scheme is said to be stable if there exists an integration time-step 
 
 
 
 Δ
 
 t
 
 0
 
 
 >
 0
 
 
 {\displaystyle \Delta t_{0}>0}
 
 so that for any 
 
 
 
 Δ
 t
 ∈
 (
 0
 ,
 Δ
 
 t
 
 0
 
 
 ]
 
 
 {\displaystyle \Delta t\in (0,\Delta t_{0}]}
 
, a finite variation of the state vector 
 
 
 
 
 q
 
 n
 
 
 
 
 {\displaystyle q_{n}}
 
 at time 
 
 
 
 
 t
 
 n
 
 
 
 
 {\displaystyle t_{n}}
 
 induces only a non-increasing variation of the state-vector 
 
 
 
 
 q
 
 n
 +
 1
 
 
 
 
 {\displaystyle q_{n+1}}
 
 calculated at a subsequent time 
 
 
 
 
 t
 
 n
 +
 1
 
 
 
 
 {\displaystyle t_{n+1}}
 
. Assume the time-integration scheme is 

 
 
 
 
 q
 
 n
 +
 1
 
 
 =
 A
 (
 Δ
 t
 )
 
 q
 
 n
 
 
 +
 
 g
 
 n
 +
 1
 
 
 (
 Δ
 t
 )
 
 
 {\displaystyle q_{n+1}=A(\Delta t)q_{n}+g_{n+1}(\Delta t)}

The linear stability is equivalent to 
 
 
 
 ρ
 (
 A
 (
 Δ
 t
 )
 )
 ≤
 1
 
 
 {\displaystyle \rho (A(\Delta t))\leq 1}
 
, here 
 
 
 
 ρ
 (
 A
 (
 Δ
 t
 )
 )
 
 
 {\displaystyle \rho (A(\Delta t))}
 
 is the <a href="/facts/Spectral_radius/j8T2xrLu">spectral radius</a> of the update matrix 
 
 
 
 A
 (
 Δ
 t
 )
 
 
 {\displaystyle A(\Delta t)}
 
.
For the linear structural equation 

 
 
 
 M
 
 
 
 u
 ¨
 
 
 
 +
 C
 
 
 
 u
 ˙
 
 
 
 +
 K
 u
 =
 
 f
 
 
 ext
 
 
 
 
 
 
 {\displaystyle M{\ddot {u}}+C{\dot {u}}+Ku=f^{\textrm {ext}}\,}

here 
 
 
 
 K
 
 
 {\displaystyle K}
 
 is the stiffness matrix. Let 
 
 
 
 
 q
 
 n
 
 
 =
 [
 
 
 
 
 u
 ˙
 
 
 
 
 n
 
 
 ,
 
 u
 
 n
 
 
 ]
 
 
 {\displaystyle q_{n}=[{\dot {u}}_{n},u_{n}]}
 
, the update matrix is 
 
 
 
 A
 =
 
 H
 
 1
 
 
 −
 1
 
 
 
 H
 
 0
 
 
 
 
 {\displaystyle A=H_{1}^{-1}H_{0}}
 
, and 

 
 
 
 
 
 
 
 
 H
 
 1
 
 
 =
 
 
 [
 
 
 
 M
 +
 γ
 Δ
 t
 C
 
 
 γ
 Δ
 t
 K
 
 
 
 
 β
 Δ
 
 t
 
 2
 
 
 C
 
 
 M
 +
 β
 Δ
 
 t
 
 2
 
 
 K
 
 
 
 ]
 
 
 
 
 H
 
 0
 
 
 =
 
 
 [
 
 
 
 M
 −
 (
 1
 −
 γ
 )
 Δ
 t
 C
 
 
 −
 (
 1
 −
 γ
 )
 Δ
 t
 K
 
 
 
 
 −
 (
 
 
 1
 2
 
 
 −
 β
 )
 Δ
 
 t
 
 2
 
 
 C
 +
 Δ
 t
 M
 
 
 M
 −
 (
 
 
 1
 2
 
 
 −
 β
 )
 Δ
 
 t
 
 2
 
 
 K
 
 
 
 ]
 
 
 
 
 
 
 
 
 {\displaystyle {\begin{aligned}H_{1}={\begin{bmatrix}M+\gamma \Delta tC&\gamma \Delta tK\\\beta \Delta t^{2}C&M+\beta \Delta t^{2}K\end{bmatrix}}\qquad H_{0}={\begin{bmatrix}M-(1-\gamma )\Delta tC&-(1-\gamma )\Delta tK\\-({\frac {1}{2}}-\beta )\Delta t^{2}C+\Delta tM&M-({\frac {1}{2}}-\beta )\Delta t^{2}K\end{bmatrix}}\end{aligned}}}

For undamped case (
 
 
 
 C
 =
 0
 
 
 {\displaystyle C=0}
 
), the update matrix can be decoupled by introducing the eigenmodes 
 
 
 
 u
 =
 
 e
 
 i
 
 ω
 
 i
 
 
 t
 
 
 
 x
 
 i
 
 
 
 
 {\displaystyle u=e^{i\omega _{i}t}x_{i}}
 
 of the structural system, which are solved by the generalized eigenvalue problem

 
 
 
 
 ω
 
 2
 
 
 M
 x
 =
 K
 x
 
 
 
 {\displaystyle \omega ^{2}Mx=Kx\,}

For each eigenmode, the update matrix becomes

 
 
 
 
 
 
 
 
 H
 
 1
 
 
 =
 
 
 [
 
 
 
 1
 
 
 γ
 Δ
 t
 
 ω
 
 i
 
 
 2
 
 
 
 
 
 
 0
 
 
 1
 +
 β
 Δ
 
 t
 
 2
 
 
 
 ω
 
 i
 
 
 2
 
 
 
 
 
 ]
 
 
 
 
 H
 
 0
 
 
 =
 
 
 [
 
 
 
 1
 
 
 −
 (
 1
 −
 γ
 )
 Δ
 t
 
 ω
 
 i
 
 
 2
 
 
 
 
 
 
 Δ
 t
 
 
 1
 −
 (
 
 
 1
 2
 
 
 −
 β
 )
 Δ
 
 t
 
 2
 
 
 
 ω
 
 i
 
 
 2
 
 
 
 
 
 ]
 
 
 
 
 
 
 
 
 {\displaystyle {\begin{aligned}H_{1}={\begin{bmatrix}1&\gamma \Delta t\omega _{i}^{2}\\0&1+\beta \Delta t^{2}\omega _{i}^{2}\end{bmatrix}}\qquad H_{0}={\begin{bmatrix}1&-(1-\gamma )\Delta t\omega _{i}^{2}\\\Delta t&1-({\frac {1}{2}}-\beta )\Delta t^{2}\omega _{i}^{2}\end{bmatrix}}\end{aligned}}}

The characteristic equation of the update matrix is 

 
 
 
 
 λ
 
 2
 
 
 −
 
 (
 
 2
 −
 (
 γ
 +
 
 
 1
 2
 
 
 )
 
 η
 
 i
 
 
 2
 
 
 
 )
 
 λ
 +
 1
 −
 (
 γ
 −
 
 
 1
 2
 
 
 )
 
 η
 
 i
 
 
 2
 
 
 =
 0
 
 
 
 η
 
 i
 
 
 2
 
 
 =
 
 
 
 
 ω
 
 i
 
 
 2
 
 
 Δ
 
 t
 
 2
 
 
 
 
 1
 +
 β
 
 ω
 
 i
 
 
 2
 
 
 Δ
 
 t
 
 2
 
 
 
 
 
 
 
 {\displaystyle \lambda ^{2}-\left(2-(\gamma +{\frac {1}{2}})\eta _{i}^{2}\right)\lambda +1-(\gamma -{\frac {1}{2}})\eta _{i}^{2}=0\,\qquad \eta _{i}^{2}={\frac {\omega _{i}^{2}\Delta t^{2}}{1+\beta \omega _{i}^{2}\Delta t^{2}}}}

As for the stability, we have 
Explicit central difference scheme (
 
 
 
 γ
 =
 0.5
 
 
 {\displaystyle \gamma =0.5}
 
 and 
 
 
 
 β
 =
 0
 
 
 {\displaystyle \beta =0}
 
) is stable when 
 
 
 
 ω
 Δ
 t
 ≤
 2
 
 
 {\displaystyle \omega \Delta t\leq 2}
 
.
Average constant acceleration (Middle point rule) (
 
 
 
 γ
 =
 0.5
 
 
 {\displaystyle \gamma =0.5}
 
 and 
 
 
 
 β
 =
 0.25
 
 
 {\displaystyle \beta =0.25}
 
) is unconditionally stable.

<h2 id="references">References</h2>

<ol>
<li id="fn:1">Newmark, Nathan M. (1959), "A method of computation for structural dynamics", Journal of the Engineering Mechanics Division, 85 (EM3) (3): 67–94, doi:10.1061/JMCEA3.0000098 <a href="/wiki/Nathan_M._Newmark" target="_blank">/wiki/Nathan_M._Newmark</a> <a href="#fnref:1" class="footnote-back-ref">↩</a></li>
</ol>

Newmark-beta method open-in-new

Newmark-beta method