Produces a discrete-time state-space model of the structural equations

\[\begin{split}\mathbf{\ddot{x}} &= \mathbf{M}^{-1}( -\mathbf{C}\,\mathbf{\dot{x}}-\mathbf{K}\,\mathbf{x}+\mathbf{F} ) \\ \mathbf{y} &= \mathbf{x}\end{split}\]

based on the Newmark-\(\beta\) integration scheme. The output state-space model has form:

\[\begin{split}\mathbf{X}_{n+1} &= \mathbf{A}\,\mathbf{X}_n + \mathbf{B}\,\mathbf{F}_n \\ \mathbf{Y} &= \mathbf{C}\,\mathbf{X} + \mathbf{D}\,\mathbf{F}\end{split}\]

with \(\mathbf{X} = [\mathbf{x}, \mathbf{\dot{x}}]^T\)

Note that as the state-space representation only requires the input force \(\mathbf{F}\) to be evaluated at time-step \(n\),the \(\mathbf{C}\) and \(\mathbf{D}\) matrices are, in general, fully populated.

The Newmark-\(\beta\) integration scheme is carried out following the modifications presented by Geradin [1] that render it unconditionally stable. The displacement and velocities are estimated as:

\[\begin{split}x_{n+1} &= x_n + \Delta t \dot{x}_n + \left(\frac{1}{2}-\theta_2\right)\Delta t^2 \ddot{x}_n + \theta_2\Delta t \ddot{x}_{n+1} \\ \dot{x}_{n+1} &= \dot{x}_n + (1-\theta_1)\Delta t \ddot{x}_n + \theta_1\Delta t \ddot{x}_{n+1}\end{split}\]

The stencil is unconditionally stable if the tuning parameters \(\theta_1\) and \(\theta_2\) are chosen as:

\[\begin{split}\theta_1 &= \frac{1}{2} + \alpha \\ \theta_2 &= \frac{1}{4} \left(\theta_1 + \frac{1}{2}\right)^2 \\ \theta_2 &= \frac{5}{80} + \frac{1}{4} (\theta_1 + \theta_1^2) \text{TBC SOURCE}\end{split}\]

where \(\alpha>0\) accounts for small positive algorithmic damping.

The following steps describe how to apply the Newmark-beta scheme to a state-space formulation. The original idea is based on [1].

The equation of a second order system dynamics reads:

\[M\mathbf{\ddot q} + C\mathbf{\dot q} + K\mathbf{q} = F\]

Applying that equation to the time steps \(n\) and \(n+1\), rearranging terms and multiplying by \(M^{-1}\):

\[\begin{split}\mathbf{\ddot q}_{n} = - M^{-1}C\mathbf{\dot q}_{n} - M^{-1}K\mathbf{q}_{n} + M^{-1}F_{n} \\ \mathbf{\ddot q}_{n+1} = - M^{-1}C\mathbf{\dot q}_{n+1} - M^{-1}K\mathbf{q}_{n+1} + M^{-1}F_{n+1}\end{split}\]

The relations of the Newmark-beta scheme are:

\[\begin{split}\mathbf{q}_{n+1} &= \mathbf{q}_n + \mathbf{\dot q}_n\Delta t + (\frac{1}{2}-\beta)\mathbf{\ddot q}_n \Delta t^2 + \beta \mathbf{\ddot q}_{n+1} \Delta t^2 + O(\Delta t^3) \\ \mathbf{\dot q}_{n+1} &= \mathbf{\dot q}_n + (1-\gamma)\mathbf{\ddot q}_n \Delta t + \gamma \mathbf{\ddot q}_{n+1} \Delta t + O(\Delta t^3)\end{split}\]

Substituting the former relation onto the later ones, rearranging terms, and writing it in state-space form:

\[\begin{split}\begin{bmatrix} I + M^{-1}K \Delta t^2\beta \quad \Delta t^2\beta M^{-1}C \\ (\gamma \Delta t M^{-1}K) \quad (I + \gamma \Delta t M^{-1}C) \end{bmatrix} \begin{Bmatrix} \mathbf{\dot q}_{n+1} \\ \mathbf{\ddot q}_{n+1} \end{Bmatrix} = \begin{bmatrix} (I - \Delta t^2(1/2-\beta)M^{-1}K \quad (\Delta t - \Delta t^2(1/2-\beta)M^{-1}C \\ (-(1-\gamma)\Delta t M^{-1}K \quad (I - (1-\gamma)\Delta tM^{-1}C \end{bmatrix} \begin{Bmatrix} \mathbf{q}_{n} \\ \mathbf{\dot q}_{n} \end{Bmatrix} + \begin{Bmatrix} (\Delta t^2(1/2-\beta) \\ (1-\gamma)\Delta t \end{Bmatrix} M^{-1}F_n+ \begin{Bmatrix} (\Delta t^2\beta) \\ (\gamma \Delta t) \end{Bmatrix}M^{-1}F_{n+1}\end{split}\]

To understand SHARPy code, it is convenient to apply the following change of notation:

\[\begin{split}\textrm{th1} = \gamma \\ \textrm{th2} = \beta \\ \textrm{a0} = \Delta t^2 (1/2 -\beta) \\ \textrm{b0} = \Delta t (1 -\gamma) \\ \textrm{a1} = \Delta t^2 \beta \\ \textrm{b1} = \Delta t \gamma \\\end{split}\]


\[\begin{split}A_{ss1} \begin{Bmatrix} \mathbf{\dot q}_{n+1} \\ \mathbf{\ddot q}_{n+1} \end{Bmatrix} = A_{ss0} \begin{Bmatrix} \mathbf{\dot q}_{n} \\ \mathbf{\ddot q}_{n} \end{Bmatrix} + \begin{Bmatrix} (\Delta t^2(1/2-\beta) \\ (1-\gamma)\Delta t \end{Bmatrix} M^{-1}F_n+ \begin{Bmatrix} (\Delta t^2\beta) \\ (\gamma \Delta t) \end{Bmatrix}M^{-1}F_{n+1}\end{split}\]

To finally isolate the vector at \(n+1\), instead of inverting the \(A_{ss1}\) matrix, several systems are solved. Moreover, the output equation is simply \(y=x\).

param Minv:

Inverse mass matrix \(\mathbf{M^{-1}}\)

type Minv:


param C:

Damping matrix \(\mathbf{C}\)

type C:


param K:

Stiffness matrix \(\mathbf{K}\)

type K:


param dt:

Timestep increment

type dt:


param num_damp:

Numerical damping. Default 1e-4

type num_damp:



the A, B, C, D matrices of the state space packed in a tuple with the predictor and delay term removed.




[1] - Geradin M., Rixen D. - Mechanical Vibrations: Theory and application to structural dynamics