## 4 Solution Procedure
### Numerical Solution  
The same problem is also solved numerically by means of the GDQ method, due to its capability to yield accurate solutions even with a reduced computational effort [41; 42; 43; 44; 45] while maintaining a certain flexibility when involving any kind of boundary condition along the structural edges. The proposed method allows solving of the problem in a strong form, by discretizing the derivatives of a function in the following form [41]:  
\[\left.\frac{\partial f}{\partial x}\right|_{x=x_{i},y=y_{j}}=\sum\limits_{m=1}^ {N_{x}}\sum\limits_{n=1}^{N_{y}}A_{im}^{x}I_{jn}^{y}f_{mn} \tag{28a}\] \[\left.\frac{\partial f}{\partial y}\right|_{x=x_{i},y=y_{j}}=\sum\limits_{m=1}^ {N_{x}}\sum\limits_{n=1}^{N_{y}}I_{im}^{x}A_{jn}^{y}f_{mn}\] (28b) \[\left.\frac{\partial}{\partial r}\left(\left.\frac{\partial f}{\partial\theta} \right|_{x=x_{i},y=y_{j}}\right)=\sum\limits_{m=1}^{N_{x}}\sum\limits_{n=1}^ {N_{y}}A_{im}^{x}A_{jn}^{y}f_{mn}\] (28c) \[\left.\frac{\partial^{2}f}{\partial x^{2}}\right|_{x=x_{i},y=y_{j}}=\sum\limits_{m=1 }^{N_{x}}\sum\limits_{n=1}^{N_{y}}B_{im}^{x}I_{jn}^{y}f_{mn}\] (28d) \[\left.\frac{\partial^{2}f}{\partial y^{2}}\right|_{x=x_{i},y=y_{j}}=\sum\limits_{m=1 }^{N_{x}}\sum\limits_{n=1}^{N_{y}}I_{im}^{x}B_{jn}^{y}f_{mn} \tag{28e}\]  
where \(I_{im}^{x}\) and \(I_{jn}^{y}\) are equal to one when \(i=m\) and \(j=n\), or equal to zero, otherwise. In addition, \(A_{im}^{x}\), \(A_{jn}^{y}\), \(B_{im}^{x}\) and \(B_{jn}^{y}\) are the weighting coefficients of the first and second-order derivatives along the \(x\) and \(y\)-directions, respectively, defined as:  
\[A_{im}^{(1)}=\left\{\begin{array}{ll}\frac{\xi(y_{i})}{(x_{i}-x_{m})\xi(x_{ m})}&\mbox{when}\ \ \ i\
eq m\\ -\sum\limits_{k=1,k\
eq i}^{N_{x}}A_{ik}^{(1)}&\mbox{when}\ \ \ i=m\end{array}\right.\ \ \ \ \ i,m=1,2,\ldots,N_{x}\] (29a) \[A_{jn}^{(1)}=\left\{\begin{array}{ll}\frac{\xi(y_{j})}{(y_{j}-y_{k})\xi(y_{ n})}&\mbox{when}\ \ \ j\
eq n\\ -\sum\limits_{k=1,k\
eq j}^{N_{y}}A_{jk}^{(1)}&\mbox{when}\ \ \ j=n\end{array}\right.\ \\[B_{jj}^{(2)}=-\sum_{k=1,k\
eq j}^{N_{y}}B_{jk}^{(2)}\ \,\ j=1,2,\ldots,N_{y},\ j=n \tag{31d}\]  
In what follows, we select a Chebyshev distribution of grid points within the domain defined as:  
\[x_{i}=\frac{a}{2}\bigg{(}1-\cos\!\bigg{(}\frac{(i-1)}{(N_{x}-1)}\pi\bigg{)} \bigg{)}\ \ \ \ \ \ i=1,2,3,\ \ldots,N_{x} \tag{32a}\] \[y_{j}=\frac{b}{2}\bigg{(}1-\cos\!\bigg{(}\frac{(j-1)}{(N_{y}-1)}\pi\bigg{)} \bigg{)}\ \ \ \ j=1,2,3,\ \ldots,N_{y} \tag{32b}\]  
Thus, the algebraic eigenvalue problem can be redefined in matrix form as:  
\[\bigg{\{}\left[\begin{array}{cc}\big{[}M_{dd}\big{]}&\big{[}M_{db}\big{]} \\ \big{[}M_{bd}\big{]}&\big{[}M_{bb}\big{]}\end{array}\right]\omega_{mn}^{2}+ \left[\begin{array}{cc}\big{[}K_{dd}\big{]}&\big{[}K_{db}\big{]}\\ \big{[}K_{bd}\big{]}&\big{]}\end{array}\right]\bigg{\}}\bigg{\{}\begin{array}[] {c}\delta_{d}\\ \delta_{b}\end{array}\bigg{\}}=0 \tag{33}\]  
where we distinguish among inner and boundary grid-points, by means of subscripts \(d\) and \(b\), respectively, whereas \(\delta\) stands for the displacement vector. The natural frequencies of the problem are derived as solutions of Equation (33).