Editing Newmark-beta method (section)

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

<math>q_{n+1} = A(\Delta t) q_n + g_{n+1}(\Delta t)</math>

The linear stability is equivalent to <math>\rho(A(\Delta t)) \leq 1</math>, here <math>\rho(A(\Delta t))</math> is the [[spectral radius]] of the update matrix <math>A(\Delta t)</math>.

For the linear structural equation 

<math>M\ddot{u} + C\dot{u} + K u = f^{\textrm{ext}} \,</math>

here <math>K</math> is the stiffness matrix. Let <math>q_n = [\dot{u}_n, u_n]</math>, the update matrix is  <math>A = H_1^{-1}H_0</math>, and 

<math>\begin{aligned}
H_1 = \begin{bmatrix}
M + \gamma\Delta tC & \gamma \Delta t K\\
\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 t K\\
-(\frac{1}{2} - \beta) \Delta t^2 C +\Delta t M & M - (\frac{1}{2} - \beta)\Delta t^2 K
\end{bmatrix}
\end{aligned}</math>

For undamped case (<math>C = 0</math>), the update matrix can be decoupled by introducing the eigenmodes <math>u = e^{i \omega_i t} x_i</math> of the structural system, which are solved by the generalized eigenvalue problem

<math>\omega^2 M x =  K x  \,</math>

For each eigenmode, the update matrix becomes

<math>\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}</math>

The characteristic equation of the update matrix is 

<math>\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}</math>

As for the stability,  we have 

'''Explicit central difference scheme''' (<math>\gamma=0.5 </math> and <math>\beta=0 </math>) is stable when <math>\omega \Delta t \leq 2</math>.

'''Average constant acceleration (Middle point rule)''' (<math>\gamma=0.5 </math> and <math>\beta=0.25 </math>) is unconditionally stable.