## 3 Stability Analysis  
Stability analysis is significant for a new TFEM, and for this analysis a spectral approach is generally employed because a linear dynamic system can be decoupled into a set of single DOF systems based on the orthogonality of mode functions or mode vectors. Therefore, in the stability analysis of the present method, consider a single DOF system as  
\[\ddot{u}+(\omega h)^{2}u=0, \tag{27}\]  
where \(\omega\) denotes the angular frequency. Here, the two dots over \(\mathbf{u}\) denote the second derivative with respect to time \(\tau\), not to \(t\), where \(\tau=t/h\), and \(h\) and \(t\) stand for time step size and the real time, respectively. For Equation (27), the displacement vector \(\mathbf{u}\) in Equation (6) becomes  
\[\mathbf{u}^{\text{T}}=[ \begin{array}{cccc}u_{1}&u_{2}&u_{3}&\cdots&u_{N}\end{array}  ], \tag{28}\]  
where \(u_{j}\) (\(j=1,\,2,\,\ldots\,\,N\)) denotes the displacements at time sampling points, the subscript '\(i\)' in Equation (6) denotes the spatial DOF number and is omitted here since \(i=1\), and  
\[\mathbf{G}=\mathbf{A}^{\text{T}}\mathbf{C}_{\text{G}}\mathbf{A}-(\omega h)^{2}\mathbf{C}_{\text{G }},\mathbf{U}=\mathbf{u},\mathbf{Q}=\mathbf{0} \tag{29}\]  
The nodal parameter vector in Equation (19) has the form of  
\[\widetilde{\mathbf{U}}^{\text{T}}=\Big{[} \begin{array}{cccccccc}u_{1}&u_{2}&u_ {3}&\cdots&u_{N-2}&u_{N-1}&\dot{u}_{1}\end{array} \Big{]}. \tag{30}\]According to Equation (26), one can obtain  
\[\widetilde{\mathbf{U}}_{\mathbf{d}}=\left[ \begin{array}{c}u_{2}\\ u_{3}\\ \vdots\\ u_{N-1}\end{array} \right]=-\left[ \begin{array}{cc}\mathbf{\alpha}&\mathbf{\beta}\end{array}  \right]\left[ \begin{array}{c}u_{0}\\ \dot{u}_{0}\end{array} \right], \tag{31}\]  
where column vectors \(\mathbf{\alpha}\) and \(\mathbf{\beta}\) are extracted from \(\widetilde{\mathbf{G}}\), and \(u_{0}\) and \(\dot{u}_{0}\) are the initial displacement and velocity, respectively. Moreover, using DQ rule yields  
\[ \begin{array}{l}\dot{u}_{1}=\dot{u}_{0}=A_{11}u_{0}+\left[\begin{array}{ ccc}A_{12}&\cdots&A_{1(N-1)}\end{array}\right]\widetilde{\mathbf{U}}_{\mathbf{d}}+A_{1N}u_{N} \\ \dot{u}_{N}=A_{N1}u_{0}+\left[\begin{array}{cc}A_{N2}&\cdots&A_{N(N-1)}\end{array} \right]\widetilde{\mathbf{U}}_{\mathbf{d}}+A_{NN}u_{N}\end{array} . \tag{32}\]  
From Equations (31) and (32), we have  
\[\left[ \begin{array}{c}u_{N}\\ \dot{u}_{N}\end{array} \right]=J\left[ \begin{array}{c}u_{0}\\ \dot{u}_{0}\end{array} \right], \tag{33}\]  
which is the recurrence relation between \((u_{N},\dot{u}_{N})\) and \((u_{0},\dot{u}_{0})\), and the elements of the Jacobi matrix \(J\) are  
\[\left\{ \begin{array}{l}J_{11}=A_{1N}^{-1}\Big{(}\Big{[}\begin{array}{ccc} A_{12}&\cdots&A_{1(N-1)}\end{array}\Big{]}\mathbf{\alpha}-A_{11}\Big{)}\\ J_{12}=A_{1N}^{-1}\Big{(}\Big{[}\begin{array}{ccc}A_{12}&\cdots&A_{1(N-1)} \end{array}\Big{]}\mathbf{\beta}+1\Big{)}\\ J_{21}=A_{N1}-\Big{[}\begin{array}{ccc}A_{N2}&\cdots&A_{N(N-1)}\end{array} \Big{]}\mathbf{\alpha}+A_{NN}J_{11}\\ J_{22}=A_{NN}J_{12}-\Big{[}\begin{array}{ccc}A_{N2}&\cdots&A_{N(N-1)}\end{array} \Big{]}\mathbf{\beta}\end{array} . \tag{34}\]  
According to the relationship between \(\omega h\) and spectral radius \(\rho(J)\), see Figure 1, the first stable intervals for different \(N\) can be obtained as given in Table 1. From Table 1 and Figure 1 we can make the following observations:  
1. \(\rho(J)\geq 1\) for any \(N\), and when \(N\) becomes larger, the first stable interval also becomes larger. Although there are several stable intervals with \(\rho(J)=1\) except for \(N=3\), the first stable interval is significant and practical; therefore, only the first intervals are listed in Table 1 whose elements are obtained by using the criterion \((\rho-1)<10^{-6}\). 2. \(\rho(J)<1\) in the TFEM by [PERSON] [35], but \(\rho(J)=1\) in DQTFEM for all stable intervals, hence DQTFEM is more accurate than that of [PERSON] [35] since DQTFEM is a high-order method that has a smaller phase error than lower order methods.