## 3 Analysis of Numerical Stability and Accuracy  
The following single-DOF system is considered to analyze the stability and accuracy of the proposed method:  
\[\ddot{u}+\omega^{2}u=0 \tag{21}\]  
where \(\omega\) is the circular frequency of free vibration. The relationship between the response of the \((n+1)\)-th time step and the \(n\)-th time step can be obtained by referring to Equations (11) and (20), that is  
\[ \begin{array}{l}u_{n+1}-c_{3}\ddot{u}_{n+1}\Delta t=u_{n}+\big{(}c_{1}\dot{u} _{n}+c_{2}\ddot{u}_{n+\beta}\big{)}\Delta t\\ \ddot{u}_{n+1}-c_{3}\ddot{u}_{n+1}\Delta t=\dot{u}_{n}+\big{(}c_{1}\ddot{u}_{n} +c_{2}\ddot{u}_{n+\beta}\big{)}\Delta t\\ \dddot{u}_{n+1}+\omega^{2}u_{n+1}=0\end{array}  \tag{22}\]  
which leads to  
\[\left[ \begin{array}{cc}1&-c_{3}\Delta t\\ c_{3}\omega^{2}\Delta t&1\end{array} \right]\left\{ \begin{array}{cc}u_{n+1}\\ \dot{u}_{n+1}\end{array} \right\}=\left[ \begin{array}{cc}1&c_{1}\Delta t\\ -c_{1}\omega^{2}\Delta t&1\end{array} \right]\left\{ \begin{array}{cc}u_{n} \\ i_{n}\end{array} \right\}+\left[ \begin{array}{cc}0&c_{2}\Delta t\\ -c_{2}\omega^{2}\Delta t&0\end{array} \right]\left\{ \begin{array}{cc}u_{n+ \beta}\\ \dot{u}_{n+\beta}\end{array} \right\} \tag{23}\]  
As we know, the transfer matrix \(\mathbf{A}\) in the following equation is generally used to analyze the stability of the method:  
\[\left\{ \begin{array}{c}u_{n+1}\\ i_{n+1}\end{array} \right\}=\mathbf{A}\left\{ \begin{array}{cc}u_{n}\\ i_{n}\end{array} \right\} \tag{24}\]  
Using \(\rho(\mathbf{A})\) to represent the spectral radius of \(\mathbf{A}\), it is well-known that when \(\rho(\mathbf{A})\leq 1\), the algorithm is stable. Referring to Equation (23), for the proposed two-step algorithm, the transfer matrix \(\mathbf{A}\) can be written into the following formula:  
\[\mathbf{A}=\left[ \begin{array}{cc}1&-c_{3}\Delta t\\ c_{3}\omega^{2}\Delta t&1\end{array} \right]^{-1}\left(\left[ \begin{array}{ cc}1&c_{1}\Delta t\\ -c_{1}\omega^{2}\Delta t&1\end{array} \right]+\left[ \begin{array}{cc}0&c_{2} \Delta t\\ -c_{2}\omega^{2}\Delta t&0\end{array} \right]\mathbf{B}\right) \tag{25}\]  
where \(\mathbf{B}\) is the transfer matrix computed in the first sub-step, i.e., \(\mathbf{B}\) satisfies  
\[\left\{ \begin{array}{c}u_{n+\beta}\\ \dot{u}_{n+\beta}\end{array} \right\}=\mathbf{B}\left\{ \begin{array}{cc}u_{n} \\ \dot{u}_{n}\end{array} \right\} \tag{26}\]  
As the transfer matrix of the GGW method, \(\mathbf{B}\) was derived from Ref. [31]:  
\[\mathbf{B}=\frac{1}{\alpha^{2}\Lambda^{2}+4\alpha\Lambda^{2}+4\Lambda^{2}+16} \left[ \begin{array}{cc}\alpha^{2}\Lambda^{2}-4\Lambda^{2}+16&16\beta\Delta t \\ -16\Lambda^{2}/\beta\Delta t&\alpha^{2}\Lambda^{2}-4\Lambda^{2}+16\end{array} \right] \tag{27}\]  
where \(\Lambda=\omega\beta\Delta t\), \(\alpha\) is the parameter in the GGW method to control the numerical damping. It has already been concluded that the spectral radius \(\rho(\mathbf{B})\leq 1\)[33; 34] when using reduced Gaussian integration in Equation (8) and the GGW method is unconditionally stable [31]. By using Equations (25) and (27), the spectral radius of the present two-step algorithm is derived:  
\[\rho(\mathbf{A})=\left|\frac{(4c_{1}+4c_{2}-\alpha\beta-2\beta)\Omega i+(2c_{ 1}\beta-2c_{2}\beta+\alpha c_{1}\beta+\alpha c_{2}\beta)\Omega^{2}+4}{[(2+ \alpha)\beta+4c_{3}]\Omega i+(2c_{3}\beta+\alpha c_{3}\beta)\Omega^{2}-4}\right| \tag{28}\]  
where \(\Omega=\omega\Delta t\). With \(T\) denoting the period of free vibration, when \(\Delta t/T\) tends to infinity, the spectral radius \(\rho(\mathbf{A})\) turns out to be  
\[ \begin{array}{l}\rho_{\infty}\equiv\rho(\mathbf{A})_{\Delta t\to\infty}= \left|\frac{(2c_{1}-2c_{2}+\alpha c_{1}+\alpha c_{2})}{(2c_{3}+\alpha c_{3})} \right|\\ =\left|\frac{c_{1}}{c_{3}}-\frac{2-\alpha}{2+\alpha}\frac{c_{2}}{c_{3}} \right|=\left|\frac{\beta}{1-\beta}-\frac{2-\alpha}{2+\alpha}\frac{1}{1-\beta}\right| \end{array}  \tag{29}\]  
From Figure 1, we can find that the value of \(\rho_{\infty}\) is controlled by changing \(\alpha\) and \(\beta\).