## 4 Dynamic Analysis
### Numerical Simulations and Stability Analysis  
Since multiple values appear in the steady-state solution, the stability of the steady-state solution should be investigated. The general approach of the CVHBM is employed here to judge the stability [29]. The response is rewritten in the form of:  
\[\overline{z}=A_{1}(\tau)\cos(\Omega\tau+\varphi(\tau)) \tag{50}\]  
Substituting Equation (50) into (24), and then by using trigonometric transformation, equating the coefficients of equivalent first-order harmonics results in two coupled differential equations:  
\[A_{1^{{}^{\prime\prime}}}=-\left\{\begin{array}{c}\left[\frac{15}{32}-\frac{ 7}{2}\sigma_{2}(\Omega+\varphi^{\prime})^{2}\right]{A_{1}}^{5}+\\ \left[3(\frac{1}{4}+\sigma_{2})+5\sigma_{2}{({A_{1}}^{\prime})}^{2}-4\sigma_{2 }(\Omega+\varphi^{\prime})^{2}\right]{A_{1}}^{3}\\ +\left[3\sigma_{2}{({A_{1}}^{\prime})}^{2}+(8-4\sigma_{1}(\Omega+\varphi^{ \prime})^{2})\right]{A_{1}}\\ +8\zeta{A_{1}}^{\prime}-4u_{0}\Omega^{2}\cos\varphi\end{array}\right\}/(4 \sigma_{1}+3\sigma_{2}{A_{1}}^{2}+\frac{5}{2}\sigma_{2}{A_{1}}^{4}) \tag{51}\]  
\[\varphi^{\prime\prime}=-\frac{(\Omega+\varphi^{\prime})\big{[}\left(8\sigma_ {1}+6\sigma_{2}{A_{1}}^{2}+5\sigma_{2}{A_{1}}^{4}\right){A_{1}}^{\prime}+16 \zeta{A_{1}}\big{]}+8u_{0}\Omega^{2}\sin\varphi}{8\sigma_{1}{A_{1}}+2\sigma_{2 }{A_{1}}^{3}+\sigma_{2}{A_{1}}^{5}} \tag{52}\]  
Rewriting Equations (51) and (52) in the form of the equation of state, the Jacobian matrix can be calculated to be:  
\[J=\left[\begin{array}{cccc}0&0&1&0\\ 0&0&0&1\\ J_{31}&-\frac{8u_{0}\Omega^{2}\sin z_{2}}{5\sigma_{2}{x_{1}}^{4}+6\sigma_{2}{x _{1}}^{2}+8\sigma_{1}}&-\frac{4\left(5\sigma_{2}{x_{1}}^{3}{x_{3}}+3\sigma_{2 }{x_{1}}{x_{3}}+4\zeta\right)}{5\sigma_{2}{x_{1}}^{4}+6\sigma_{2}{x_{1}}^{2}+8 \sigma_{1}}&\frac{(\Omega+x_{3})\left(14\sigma_{2}{x_{1}}^{5}+16\sigma_{2}{x_{ 1}}^{3}+16\sigma_{1}{x_{1}}\right)}{5\sigma_{2}{x_{1}}^{4}+6\sigma_{2}{x_{1}}^ {2}+8\sigma_{1}}\\ J_{41}&-\frac{8u_{0}\Omega^{2}\cos z_{2}}{8\sigma_{1}{x_{1}}{x_{2}}+2\sigma_{2}{x _{1}}^{3}+\sigma_{2}{x_{1}}^{3}}&-\frac{(\Omega+x_{4})\left(5\sigma_{2}{x_{1}} ^{4}+6\sigma_{2}{x_{1}}^{2}+8\sigma_{1}\right)}{8\sigma_{1}{x_{1}}+2\sigma_{2} {x_{1}}^{3}+\sigma_{2}{x_{1}}^{5}}&-\frac{16\xi_{1}{x_{1}}+3\left(5\sigma_{2}{x _{1}}^{4}+6\sigma_{2}{x_{1}}^{2}+8\sigma_{1}\right)}{8\sigma_{1}{x_{1}}+2\sigma _{2}{x_{1}}^{3}+\sigma_{2}{x_{1}}^{5}}\end{array}\right] \tag{53}\]  
where \(A_{1}\to x_{1},\varphi\to x_{2},A_{1}^{\prime}\to x_{3},\varphi^{\prime}\to x_{4}\)  
the excitation amplitude decreases, and the decline rate of their amplitude displacements basically become consistent. Figure 10 shows the frequency response curves under different excitation amplitudes and the red dotted lines of the FRCs represent the unstable regions. It can be seen that the amplitudes of displacement for the frequency response curves are varied under different excitation amplitudes and the corresponding peak frequency variation basically remains unchanged. When the excitation amplitude is larger than a certain value, there are three steady-state solutions occurring at the same frequency, two of which are stable and one unstable.