## 3 Solution Procedure
### The GDQ Method  
The GDQ method approximates the fundamental system of differential equations, by discretizing the derivatives of a function \(J(x)\) respect to a spatial variable at a given discrete grid distribution, by means of the weighting coefficients. For a one-dimensional problem where the whole domain \(0\leq x\leq L\) is discretized in \(N\) grids points, the approximation of the \(n\)th-order derivatives of \(J\) function respect to \(x\) variable can be expressed as [49]  
\[\frac{d^{n}J\big{(}x_{p}\big{)}}{dx^{n}}=\sum_{r=1}^{N}\chi_{pr}^{(n)}J(x_{r}) \quad n=1,2,\ldots,N-1, \tag{29}\]  
\(\chi_{pr}^{(n)}\) being the weighting coefficients, defined as follows [49]  
\[\chi_{pq}^{(1)}=\frac{\mathrm{Y}\big{(}x_{p}\big{)}}{\big{(}x_{p}-x_{q}\big{)} \mathrm{Y}\big{(}x_{q}\big{)}},\;\;p\
eq q;\;p,q=1,2,\ldots,N \tag{30}\]  
and \(\mathrm{Y}\big{(}x_{p}\big{)}\) is the Lagrangian operator expressed as [49]  
\[\mathrm{Y}\big{(}x_{p}\big{)}=\prod_{p\
eq q,\;q=1}^{N}\big{(}x_{p}-x_{q}\big{)}. \tag{31}\]  
For higher order derivatives of the weighting coefficients it is [49]  
\[\chi_{pq}^{(n)}=\eta\Bigg{(}\chi_{pp}^{(n-1)}\chi_{pq}^{(1)}-\frac{\chi_{pq}^{ (n-1)}}{\big{(}x_{p}-x_{q}\big{)}}\Bigg{)} \tag{32}\]  
It is well known in literature that the type of grid distribution within the domain can affect significantly the accuracy of the proposed method [50]. In what follows we apply a Chebyshev-Gauss-Lobatto non-uniform pattern, due to its great performances, as verified by [PERSON] [49]  
\[x_{p}=\frac{L}{2}\Bigg{(}1-\cos\Bigg{(}\frac{x_{p}-1}{N-1}\Bigg{)}\pi\Bigg{)} \quad p=1,2,\ldots,N. \tag{33}\]  
Thus, the governing differential equations of motion and boundary conditions are discretized according to the GDQ approach, as detailed in Appendix B. Let us denote the periodic axial compressive load as  
\[F_{a}(t)=\alpha F_{cr}+\beta F_{cr}\cos(\omega t) \tag{34}\]  
where \(\alpha\) and \(\beta\) refer to the static and dynamic load factors, respectively. Furthermore, \(F_{cr}\) denotes the critical static load and \(\omega\) stands for the excitation frequency. By substituting Equation (34) into the third Equation (A3) from the Appendix A, and by combining the discretized equations of motion along with the associated boundary conditions, the problem can be redefined in the following matrix form  
\[\left(\begin{array}{cc}M_{bb}&M_{bd}\\ M_{db}&M_{dd}\end{array}\right)\left(\begin{array}{c}\ddot{\Gamma}_{b}\\ \ddot{\Gamma}_{d}\end{array}\right)+\left(\left(\begin{array}{cc}K_{bb}&K_{ bd}\\ K_{db}&K_{dd}\end{array}\right)+F_{cr}(\alpha+\beta\cos(\omega t))\left(\begin{array}{ cc}K_{bb}^{G}&K_{bd}^{G}\\ K_{db}^{G}&K_{dd}^{G}\end{array}\right)\right)\left(\begin{array}{c}\Gamma_{b }\\ \Gamma_{d}\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right) \tag{35}\]  
where \(M\), \(K\), and \(K^{G}\) are the mass, stiffness, and geometric stiffness matrixes, respectively, and \(\Gamma=\left[\overline{U},\overline{V},\overline{W},\overline{\Phi}_{x}, \overline{\Phi}_{y}\right]^{T}\) denotes the unknown dynamic displacement vector. In addition, indexes \(b\) and \(d\) indicate the boundary points and domain points, respectively.