## 2 Basic Theory and Calculation Formula
In order to study the stability of the moving beam, it is usually hypothesized that the form of free vibration response is:  
\[D_{1}\frac{\partial^{3}w}{\partial x^{3}}-I_{1}\frac{\partial^{3}w}{\partial t ^{2}\partial x}-D_{0}\frac{\partial^{2}u}{\partial x^{2}}=0 \tag{12}\]1  
where, \(\Omega\) is the normalized complex frequency and \(\widetilde{W}(\xi)\) is the modal function. By substituting Equation (20) into Equation (18), we have:  
\[D_{1}\frac{\partial^{3}w}{\partial x^{3}}-I_{1}\frac{\partial^{3}w}{\partial t ^{2}\partial x}-D_{0}\frac{\partial^{2}u}{\partial x^{2}}=0 \tag{12}\]2  
The natural frequency \(\Omega\) of the complex form can be expressed in the following form:  
\[D_{1}\frac{\partial^{3}w}{\partial x^{3}}-I_{1}\frac{\partial^{3}w}{\partial t ^{2}\partial x}-D_{0}\frac{\partial^{2}u}{\partial x^{2}}=0 \tag{12}\]3  
The instability type of an axially moving FGM beam can be determined by the signs of the real part (Re) and imaginary part (Im) of all eigenvalues in Equation (21). Specifically, when \(\mathrm{Re}(\Omega)=0\), \(\mathrm{Im}(\Omega)\
eq 0\), it is static stability; when \(\mathrm{Re}(\Omega)\
eq 0\), \(\mathrm{Im}(\Omega)=0\), it is a static instability (divergence instability); and when \(\mathrm{Re}(\Omega)\
eq 0\), \(\mathrm{Im}(\Omega)\
eq 0\), it is a dynamic instability (coupled flutter). In order to facilitate the derivation and simplification of the following equations, let  
\[D_{1}\frac{\partial^{3}w}{\partial x^{3}}-I_{1}\frac{\partial^{3}w}{\partial t ^{2}\partial x}-D_{0}\frac{\partial^{2}u}{\partial x^{2}}=0 \tag{12}\]4  
By substituting Equation (22) into Equation (21), it can be simplified as,  
\[D_{1}\frac{\partial^{3}w}{\partial x^{3}}-I_{1}\frac{\partial^{3}w}{\partial t ^{2}\partial x}-D_{0}\frac{\partial^{2}u}{\partial x^{2}}=0 \tag{12}\]5  
The boundary conditions at both ends of the axially moving FGM beam after normalization are,  
\[D_{1}\frac{\partial^{3}w}{\partial x^{3}}-I_{1}\frac{\partial^{3}w}{\partial t ^{2}\partial x}-D_{0}\frac{\partial^{2}u}{\partial x^{2}}=0 \tag{12}\]6 \[D_{1}\frac{\partial^{3}w}{\partial x^{3}}-I_{1}\frac{\partial^{3}w}{\partial t ^{2}\partial x}-D_{0}\frac{\partial^{2}u}{\partial x^{2}}=0 \tag{12}\]7 (24b) \[D_{1}\frac{\partial^{3}w}{\partial x^{3}}-I_{1}\frac{\partial^{3}w}{\partial t ^{2}\partial x}-D_{0}\frac{\partial^{2}u}{\partial x^{2}}=0 \tag{12}\]8  
To sum up, the natural frequency calculation of the axially moving FGM beam is transformed into the solution of the eigenvalue of the ordinary differential Equation (23)under the boundary condition in Equation (24). Here, the interpolating matrix method is used for solution [22,23]. From a static state, with the increase of beam axial velocity, each order natural frequency of FGM beam decreases gradually. When the axial velocity reaches a certain value so that the first-order natural frequency disappears, the velocity is the critical velocity. If the FGM beam moves in the axial direction at a faster rate than this critical velocity, the beam will lose static stability. According to this law, in Equation (21), if we let \(\Omega=0\), then the governing equation for solving the critical velocity of the axially moving FGM beam is,  
\[D_{1}\frac{\partial^{3}w}{\partial x^{3}}-I_{1}\frac{\partial^{3}w}{\partial t ^{2}\partial x}-D_{0}\frac{\partial^{2}u}{\partial x^{2}}=0 \tag{12}\]9  
Thus, the interpolating matrix method (IMM) is used to solve the ordinary differential Equation (25) under the boundary condition in Equation (24), and the critical velocity of the axially moving FGM beam can be calculated. Calculation of the Natural Frequency of an Axially Moving FGM Beam Using the Interpolating Matrix Method  
The interpolating matrix method (IMM) [22,23] is an effective numerical calculation method for solving eigenvalues and two-point boundary values of ordinary differential equations. The method takes the highest derivative function in the ordinary differential equation as the unknown parameter of the discrete system.