## 3 Mathematical Investigations  
### Perturbation Analysis  
This section employs the multiple scales perturbation technique to determine an approximate solution for the nonlinear dynamical system controlled by the proposed PPF method outlined in Equations (2) and (3). This approach allows us to derive a first-order approximate solution, as detailed in references [25; 26; 27].  
\[\left.\begin{array}{l}t=\varepsilon^{n}T_{n},\\ x(t;\varepsilon)=x_{0}(T_{0},T_{1})+\varepsilon x_{1}(T_{0},T_{1})+O( \varepsilon^{2}),\\ y(t;\varepsilon)=y_{0}(T_{0},T_{1})+\varepsilon y_{1}(T_{0},T_{1})+O( \varepsilon^{2}).\end{array}\right\} \tag{4}\]  
Within this setting, the \"fast scale\" is represented by \(T_{0}\), and the \"slow scale\" is denoted as \(T_{1}=\varepsilon t\). When employing (MTSM), the derivatives assume the following forms:  
\[\left.\begin{array}{l}\frac{d}{dt_{2}}=D_{0}+\varepsilon D_{1}+\varepsilon^ {2}D_{2}+\ldots,\\ \frac{d^{2}}{dt^{2}}=D_{0}^{2}+2\varepsilon D_{0}D_{1}+\ldots,\end{array} \right\}\quad D_{j}=\frac{\partial}{\partial T_{j}}\ \ (j=0,1). \tag{5}\]  
Announcing equivalences (4) and (5) in Equations (2) and (3) such that  
\[\left(D_{0}{}^{2}+{\omega_{1}}^{2}\right)x_{0}+\varepsilon\left(D_{0}{}^{2}+{ \omega_{1}}^{2}\right)x_{1} =\varepsilon f\sin(2\omega t)+\ \lambda_{1}y_{0}-2D_{0}D_{1}x_{0}-\alpha_{1}D_{0}x_{0}- \eta{x_{0}}^{2} \tag{6}\] \[-\beta{x_{0}}^{3}-\gamma{x_{0}}^{4}-\mu{x_{0}}^{5}+O(\varepsilon^ {2}),\]  
\[\left(D_{0}{}^{2}+{\omega_{2}}^{2}\right)y_{0}+\varepsilon\left(D_{0}{}^{2}+{ \omega_{2}}^{2}\right)y_{1}=\varepsilon(\lambda_{2}x_{0}-2D_{0}D_{1}y_{0}- \alpha_{2}D_{0}y_{0})+O(\varepsilon^{2}). \tag{7}\]  
Comparing the quantities of the same power of \(\varepsilon\):  
\[O(\varepsilon^{0})\]  
\[(D_{0}{}^{2}+{\omega_{1}}^{2})x_{0}=0, \tag{8}\]  
\[\left(D_{0}{}^{2}+{\omega_{2}}^{2}\right)y_{0}=0. \tag{9}\]  
\[O(\varepsilon)\]  
\[\left(D_{0}{}^{2}+{\omega_{1}}^{2}\right)x_{1}=f\sin(2\omega t)+\lambda_{1}y_{ 0}-2D_{0}D_{1}x_{0}-\alpha_{1}D_{0}x_{0}-\eta{x_{0}}^{2}-\beta{x_{0}}^{3}- \gamma{x_{0}}^{4}-\mu{x_{0}}^{5}, \tag{10}\]  
\[\left(D_{0}{}^{2}+{\omega_{2}}^{2}\right)y_{1}=\lambda_{2}x_{0}-2D_{0}D_{1}y_{ 0}-\alpha_{2}D_{0}y_{0}. \tag{11}\]  
From Equations (8) and (9) Resolution the homogenous differential equations we obtain the following:  
\[x_{0}(T_{0},T_{1})=A(T_{1})\ e^{i\omega_{1}(T_{0})}+\overline{A}(T_{1})\ e^{-i \omega_{1}(T_{0})}, \tag{12}\]  
\[y_{0}(T_{0},T_{1})=B(T_{1})\ e^{i\omega_{2}(T_{0})}+\overline{B}(T_{1})\ e^{-i \omega_{2}(T_{0})}. \tag{13}\]  
Differentiate Equations (12) and (13) with respect to \(t\) and submit in Equations (10) and (11).  
\[\left(D_{0}{}^{2}+{\omega_{1}}^{2}\right)x_{1} =\left(-2i\omega_{1}DA-\alpha_{1}i\omega_{1}A-3\beta A^{2} \overline{A}-10\mu A^{3}\overline{A}^{2}\right)e^{i\omega_{1}(T_{0})}-2\eta A \overline{A} \tag{14}\] \[-\left(\eta A^{2}+4\gamma A^{3}\overline{A}\right)e^{2i\omega_{1 }(T_{0})}+f\sin(2\omega t)+\lambda_{1}Be^{i\omega_{2}(T_{0})}-6\gamma A^{2} \overline{A}^{2}\] \[-\left(\beta A^{3}+5\mu A^{4}\overline{A}\right)e^{3i\omega_{1}( T_{0})}-4\gamma A^{4}e^{4i\omega_{1}(T_{0})}-\mu A^{5}e^{5i\omega_{1}(T_{0})}+C.C.,\] \[\left(D_{0}{}^{2}+{\omega_{2}}^{2}\right)y_{1}=\lambda_{2}Ae^{i \omega_{1}(T_{0})}-2i\omega_{2}DBe^{i\omega_{2}(T_{0})}-i\alpha_{2}\omega_{2} Be^{i\omega_{2}(T_{0})}+C.C. \tag{15}\]  
The components that are complex conjugates are gathered in the term denoted as C.C.  
After eliminating, the secular terms take the followings forms:  
\[\left(-2i\omega_{1}DA-\alpha_{1}i\omega_{1}A-3\beta A^{2}\overline{A}-10\mu A^{3 }\overline{A}^{2}\right)e^{i\omega_{1}(T_{0})}+f\sin(2\omega t)+\lambda_{1}Be^{ i\omega_{2}(T_{0})}=0, \tag{16}\]  
\[(-2i\omega_{2}DB-i\alpha_{2}\omega_{2}B)e^{i\omega_{2}(T_{0})}+\lambda_{2}Ae^{ i\omega_{1}(T_{0})}=0. \tag{17}\]  
Based on the primary estimate, we determined the resonance cases that follow include  
* SHR: \(\omega\cong\omega_{1}/2\).
* IR: \(\omega_{1}\cong\omega_{2}\).