Weakly resonant double Hopf bifurcation in coupled nonlinear systems with delayed freedback and application of homotopy analysis method

Abstract

In this paper, we study the dynamical behaviors of coupled nonlinear systems with delay coupling by the multiple scales method and the homotopy analysis method. Firstly, we analyze the distribution of the eigenvalues of its linearized characteristic equations, and obtain the critical value for the occurrence of double Hopf bifurcation, which is caused by time delay and strength of coupling. Second, we obtain the normal form equations by the multiple scales method, and study the dynamical behaviors around the 3:5 weakly resonant double Hopf bifurcation point by analyzing the normal form equations. Finally, using the homotopy analysis method, we obtain analytical approximate solutions of the system with parameter values located in different regions. The periodic solution obtained by the homotopy analysis method is compared with the periodic solution obtained by the Runge–Kutta method, we found that the Runge–Kutta method does not get unstable periodic solutions, but the homotopy analysis method can be. So the homotopy analysis method is a powerful tool for studying coupled nonlinear systems with delay coupling.

Introduction

Due to the practical needs of engineering applications, more and more researchers are devoted to the study of nonlinear dynamics systems, especially coupled nonlinear systems^1,2 and time-delay nonlinear systems.^3,4 In this paper, we study the coupled nonlinear systems with delay coupling, which is a relatively complicated issue.

In 1994, Bélair and Campbell ⁵ studied that a first-order differential system with two time delays can present double Hopf bifurcation. Then, Xu and Chung ⁶ declared that a time delay and feedback gain also can make double Hopf bifurcation to occur. In 2003, normal form theory was used to study non-resonance double Hopf bifurcation for a kind of differential equation with delay by Buono and Bélair. ⁷ Szalai and Stépán ⁸ presented a closed-form calculation for the analysis of the period-doubling bifurcation in the time-periodic delay-differential equation model of interrupted machining processes such as milling where the nonlinearity is essentially nonsymmetric. He et al. ⁹ investigated the issue of impulsive stabilization and Hopf bifurcation of a new three-dimensional chaotic system. Ma et al. ¹⁰ applied geometric method to study double Hopf bifurcation for van der Pol-Duffing oscillator with parametric delay feedback control. In recent years, more and more scholars have begun to pay close attention to the complex dynamic behaviors of coupling systems with time delay. For example, Yan and Chu ¹¹ studied a delayed Lotka–Volterra predator–prey system, and got the linear stability of equilibrium point and the character of the Hopf bifurcation solutions. Barrón and Sen ² discussed the synchronization of four coupled van der Pol oscillators. Song et al. ¹² and Zhang et al. ¹³ respectively studied multiple Hopf bifurcations of three coupled van der Pol oscillators with delay. Wang and Xu ¹⁴ successfully applied MMS to study double Hopf bifurcation with 1:3 resonance of coupled van der Pol oscillators with delay. That means, multiple scales method (MSM) also applies to nonlinear systems with time-delay coupling. Based on the normal form method and center manifold theorem, Li et al. ¹⁵ studied the double Hopf bifurcation of a kind of coupled systems with time delay, and introduced the bifurcation diagram by numerical simulation. Ge and Xu ¹⁶ studied the weakly resonance double Hopf bifurcation for a class of delay differential equations with four dimensions. Wang and Chen ¹⁷ applied MSM to study the weakly resonance and non-resonance double Hopf bifurcation for a kind of seven coupled van der Pol systems with time delay. Amer et al. ¹⁸ study the nonlinear vibrations of a parametric excited Duffing oscillator with time-delay feedback by the MSM.

At present, the methods of studying the nonlinear system include harmonic balance method,¹⁹ homotopy analysis method (HAM), ²⁰ perturbation incremental method, ²¹ MSM, ²² and so on. In 1992, the HAM was first put forward by Liao. ²³ Different from other analytical approximation method, the HAM not only overcomes the limitation of traditional methods relying on small parameter perturbation but also provides us with an additional way to conveniently adjust and control the convergence region and rate of solution series. Recently, the HAM has been widely applied for coupled nonlinear systems and time-delay nonlinear systems. For instance, You and Xu ²⁴ applied HAM to study periodic solutions of delayed differential equations that describe time-delayed position feedback on the Duffing system. It is found that the current technique leads to higher accurate prediction on the local dynamics of time-delayed systems near a Hopf bifurcation than the energy analysis method. Li ²⁵ applied HAM to solve the coupled KdV equations. Kheiri and Jabbari ²⁶ combined HAM with Padé approximation method to study a kind of two-dimensional coupled Burgers' equations. Tasbozan et al. ²⁷ obtained approximate analytical solutions of fractional coupled modifed Korteweg de Vries (mKdV) equation by HAM. Jafarian et al. ²⁸ adopted HAM to solve coupled Ramani equations, and got satisfying approximate analytical solutions. Eigoli and Khodabakhsh²⁹ studied the limit cycle of the van der Pol oscillator with delayed amplitude limiting by HAM. Bel and Reartes ³⁰ considered a special type of second-order delay differential equations. Cobiaga and Reartes^31,32 applied HAM to search periodic orbits in delay differential equations. Nave et al. ³³ improved HAM to study the pressure driven flame in porous media. Rajeev ³⁴ discussed the mathematical model of Stefan problem with fractional order derivatives and obtained an approximate solution of the problem by the HAM. Mahendra and Patel ³⁵ discussed one-dimensional nonlinear partial differential equation for finger-imbibition phenomenon arising in a double-phase flow through a heterogeneous porous medium during secondary oil recovery process, and solved the equation with appropriate initial and boundary conditions by the HAM. Singh et al. ³⁶ applied the HAM to obtain approximate analytical solutions of the (1 + 1) dimensional nonlinear Boussinesq equation having a fractional time derivative. By taking proper values of the auxiliary and homotopy parameters, its numerical values of the state variables are computed and presented graphically for different particular cases. Fei et al. ³⁷ analyzed a sinusoidal excited piecewise linear–nonlinear oscillator, and proposed an approximate solution for the oscillator by using the HAM and matching method. The validity criteria for the application of the semi-analytic asymptotic methods are exploited, ³⁸ comparison between the solutions obtained by the two asymptotic techniques, that is, the fractional homotopy analysis transform method and the optimal HAM is performed to select the most accurate technique for the stated problem.

In this paper, attention is paid to a system of two coupled nonlinear systems with delay coupling, as follows:

${\ddot{x}}_1(t)+\left(\gamma x_1^2(t)-\beta \right){\dot{x}}_1(t)+{\overline{\omega }}_0^2{x}_1(t)=\alpha {\dot{x}}_2(t-\tau )$

(1a)

${\ddot{x}}_2(t)+\left(\gamma x_2^2(t)-\beta \right){\dot{x}}_2(t)+{\overline{\omega }}_0^2{x}_2(t)=\alpha {\dot{x}}_1(t-\tau )$

(1b)where, ${\overline{\omega }}_0>0$ is the natural frequency, $\alpha <0$ is coupling intensities, $\tau$ is time delay, $\beta$ is the damping coefficient, and $\gamma$ is nonlinear coefficient.

The constitute of the article is as follows: In The critical value for the occurrence of double Hopf bifurcation section, we analyze the distribution of the eigenvalues of its linearized characteristic equations, and get the critical value for the occurrence of double Hopf bifurcation by the Hopf bifurcation theorem, which is caused by time delay and strength of coupling. In section 3, we obtain the normal form equations by the MSM and analyzed it. At the same time, the dynamic behavior around 3:5 weak resonance double Hopf bifurcation point is studied. In application of HAM section, we obtain analytical approximate solutions of the system with parameter values located in different regions using the HAM. In numerical simulation section, via numerical simulation, the periodic solutions obtained by HAM are compared with those of Runge–Kutta method, we find that the HAM is more suitable for this system.

The critical value for the occurrence of double Hopf bifurcation

In this section, the critical value $({\alpha }_c,{\tau }_c)$ where double Hopf bifurcation occurs is obtained by studying the distribution of the eigenvalues of the associated linearizing characteristic equation.

According to the definition of double Hopf bifurcation and the bifurcation analysis of a four-neuron delayed system, the following conclusions hold.

Conclusion 1: If there occurs double Hopf bifurcation of system (1), it must meet ${\alpha }^2>{\beta }^2$ .

Proof: Let $x_1=y_1,{\dot{x}}_1=y_2,x_2=y_3,{\dot{x}}_2=y_4,$ then the system (1) can be represented as:

${\dot{y}}_1=y_2$

(2a)

${\dot{y}}_2=-{{\overline{\omega }}_0}^2{y}_1+\beta y_2-\gamma y_1^2{y}_2+\text{α}{y}_{4\tau }$

(2b)

${\dot{y}}_3=y_4$

(2c)

${\dot{y}}_4=-{{\overline{\omega }}_0}^2{y}_3+\beta y_4-\gamma y_3^2{y}_4+\text{α}{y}_{2\tau }$

(2d)

Let

$Y(t)={(y_1,y_2,y_3,y_4)}^{\text{T}}$

(3a)

$C=\left(\begin{array}{cccc}0& 1& 0& 0\\ -{{\overline{\omega }}_0}^2& \beta & 0& 0\\ 0& 0& 0& 1\\ 0& 0& -{{\overline{\omega }}_0}^2& \beta \end{array}\right)$

(3b)

$D=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& \alpha \\ 0& 0& 0& 0\\ 0& \alpha & 0& 0\end{array}\right)$

(3c)

$f\left(Y\left(t\right),Y\left(t-\tau \right)\right)=\left(\begin{array}{c}0\\ -\gamma y_1^2{y}_2\\ 0\\ -\gamma y_3^2{y}_4\end{array}\right)$

(3d)

Then the system (1) can be further represented as

$\dot{Y}(t)=C\cdot Y(t)+D\cdot Y(t-\tau )+f(Y(t),Y(t-\tau ))$

(4)

Hence, the characteristic equation of linear part of the system is as follows:

$\mathit{det}(\lambda \cdot E-C-D\cdot e^{-\lambda \tau })=0$

(5)where $E$ is a unit matrix.

That means

$\begin{array}{cccc}\lambda & -1& 0& 0\\ {\overline{\omega }}_0^2& \lambda -\beta & 0& -\alpha \cdot e^{-\lambda \tau }\\ 0& 0& \lambda & -1\\ 0& -\alpha \cdot e^{-\lambda \tau }& {\overline{\omega }}_0^2& \lambda -\beta \end{array}|=0$

(6)

Simplifying it, we can get

${({\lambda }^2+{\overline{\omega }}_0^2-\beta \lambda )}^2={(\alpha \lambda e^{-\lambda \tau })}^2$

(7)

As we know, if there occurs double Hopf bifurcation of system (1), then it means that equation (7) has two and only two purely imaginary eigenvalues.

Since ${{\overline{\omega }}_0}^2>0$ , $\lambda =0$ is not a solution of equation (7), that is, equation (7) has no zero root.

Assume that the pure imaginary numbers $\lambda =\pm i\overline{\omega }$ , $\overline{\omega }>0$ are pair solutions of equation (7), and substitute them into equation (7). Based on $e^{i\theta }=\text{cos}\theta +i\text{sin}\theta$ , separate real and imaginary parts, then we can obtain

${\beta }^2{\overline{\omega }}^2-{\left({\overline{\omega }}^2-{\overline{\omega }}_0^2\right)}^2={\alpha }^2{\overline{\omega }}^2\text{cos}2\overline{\omega }\tau$

(8a)

$2\beta \overline{\omega }({\overline{\omega }}^2-{\overline{\omega }}_0^2)={\alpha }^2{\overline{\omega }}^2\text{sin}2\overline{\omega }\tau$

(8b)

According to ${\text{sin}}^{2}2\overline{\omega }\tau +{\text{cos}}^{2}2\overline{\omega }\tau =1,$ we can get

${[{\beta }^2{\overline{\omega }}^2-{({\overline{\omega }}^2-{\overline{\omega }}_0^2)}^2]}^2+4{\beta }^2{\overline{\omega }}^2{({\overline{\omega }}^2-{\overline{\omega }}_0^2)}^2={({\alpha }^2{\overline{\omega }}^2)}^2$

(9)

Let $w={\overline{\omega }}^2>0,$ and substitute it into equation (9), we can gain

${[{\beta }^{2}w-{(w-{\overline{\omega }}_0^2)}^2]}^2+4{\beta }^{2}w{(w-{\overline{\omega }}_0^2)}^2={({\alpha }^{2}w)}^2$

(10)

Considering ${\overline{\omega }}_0^2>0,{\alpha }^2\ge 0,{\beta }^2\ge 0,$ we can come to a conclusion: When ${\alpha }^2>{\beta }^2$ , equation (10) has two and only two nonnegative roots, as follows

$w_1=\frac{2{\overline{\omega }}_0^2+{\alpha }^2-{\beta }^2-\sqrt{{\alpha }^2-{\beta }^2}\sqrt{{\alpha }^2-{\beta }^2+4{\overline{\omega }}_0^2}}{2}$

(11)

$w_2=\frac{2{\overline{\omega }}_0^2+{\alpha }^2-{\beta }^2+\sqrt{{\alpha }^2-{\beta }^2}\sqrt{{\alpha }^2-{\beta }^2+4{\overline{\omega }}_0^2}}{2}$

(12)

Accordingly, equation (7) has two and only two purely imaginary eigenvalues

${\lambda }_{1,2}=\pm i{\omega }_{-}$

(13a)

${\lambda }_{3,4}=\pm i{\omega }_{+}$

(13b)where, ${\omega }_{-}=\sqrt{w_1}>0$ , ${\omega }_{+}=\sqrt{w_2}>0$ , ${\omega }_{-}\ne {\omega }_{+}$ .

Hence, a conclusion can be drawn that if there occurs double Hopf bifurcation of system (1), it must meet ${\alpha }^2>{\beta }^2$ . □

Further, let ${\omega }_{-}:{\omega }_{+}=k_1:k_2$ , where $k_1$ and $k_2$ are relatively prime positive integers, and $k_1<k_2$ . Accordingly, $w_1:w_2=k_1^2:k_2^2$ . Substituting them into equations (11) and (12), we can gain

$\frac{k_1^2}{k_2^2}+\frac{k_2^2}{k_1^2}-2=\frac{{({\alpha }^2-{\beta }^2)}^2+4{\overline{\omega }}_0^2({\alpha }^2-{\beta }^2)}{{\overline{\omega }}_0^4}$

(14)

It is clear that when ${\overline{\omega }}_0$ and $\beta$ are selected, based on the above formula, $\alpha$ can be expressed as a function of $k_1:k_2$ , such as $\alpha =\alpha (k_1:k_2)$ .

Substituting $\alpha =\alpha (k_1:k_2)$ into equations (12) and (13), we can accordingly get ${\omega }_{-}={\omega }_{-}(k_1:k_2)$ and ${\omega }_{+}={\omega }_{+}(k_1:k_2)$ . Substituting them into equation (8), we can gain

${\tau }_{-}=\frac{1}{2{\omega }_{-}}(2j\pi -\mathit{arccos}\frac{{\beta }^2{\omega }_{-}^2-{({\omega }_{-}^2-{\overline{\omega }}_0^2)}^2}{{\alpha }^2{\omega }_{-}^2}),j=1,2,\cdots$

(15)

${\tau }_{+}=\frac{1}{2{\omega }_{+}}(2j\pi +\mathit{arccos}\frac{{\beta }^2{\omega }_{+}^2-{({\omega }_{+}^2-{\overline{\omega }}_0^2)}^2}{{\alpha }^2{\omega }_{+}^2}),j=1,2,\cdots$

(16)

Based on the Hopf bifurcation theorem, a sufficient condition of double Hopf bifurcation of system (equation (1)) is ${\tau }_c={\tau }_{-}={\tau }_{+}$ , that is, when equations (15) and (16) are equal. According to this conclusion, we can get the suitable numerical value of $k_1:k_2$ , and then determine whether the system has the double Hopf bifurcation with $k_1:k_2$ resonance.

Conclusion 2: When ${\overline{\omega }}_0=1$ , $\beta =0$ , there exists a double Hopf bifurcation with $k_1:k_2$ resonance of system (1), and the bifurcation point is

$\left({\alpha }_c,{\tau }_c\right)=\left(\frac{k_1-k_2}{k_2}\sqrt{\frac{k_2}{k_1}},\frac{k_1\pi }{k_1-k_2}\sqrt{\frac{k_2}{k_1}}\right)$

(17)where $k_1=2j-1$ , $k_2=2j+1$ , $j=1,2,\cdots$ .

Proof: Simplicity, Let

${\overline{\omega }}_0=1,\beta =0$

(18)then

$\alpha =\frac{k_1-k_2}{k_2}\sqrt{\frac{k_2}{k_1}}$

(19a)

${\omega }_{-}=\sqrt{\frac{k_1}{k_2}}$

(19b)

${\omega }_{+}=\sqrt{\frac{k_2}{k_1}}$

(19c)

And further, we can obtain

${\tau }_{-}=\frac{2j\pi -\pi }{2}\sqrt{\frac{k_2}{k_1}},\hspace{1em}j=1,2,\cdots$

(20)

${\tau }_{+}=\frac{2j\pi +\pi }{2}\sqrt{\frac{k_1}{k_2}},\hspace{1em}j=1,2,\cdots$

(21)

Differentiating $\lambda$ with respect to $\tau$ in equation (7), these can be obtained

${\lambda }^{\prime }(\tau )\left[1-\frac{1}{\lambda {(\tau )}^2}-\frac{\tau \left(\lambda {(\tau )}^2+1\right)}{\lambda (\tau )}\right]=-\lambda {(\tau )}^2-1$

(22)and

${\mathit{Re}\left(\frac{1}{{\lambda }^{\prime }(\tau )}\right)|}_{\tau ={\tau }_{-}}=\frac{k_2(k_1+k_2)}{k_1(k_1-k_2)}<0$

(23a)

${\mathit{Re}\left(\frac{1}{{\lambda }^{\prime }(\tau )}\right)|}_{\tau ={\tau }_{+}}=\frac{k_1(k_1+k_2)}{k_2(k_2-k_1)}>0$

(23b)

Then

${sign\left(\mathit{Re}\left(\frac{\text{d}\lambda (\tau )}{\text{d}\tau }\right)\right)|}_{\tau ={\tau }_{\pm }}={sign\left(\mathit{Re}\left(\frac{1}{{\lambda }^{\prime }(\tau )}\right)\right)|}_{\tau ={\tau }_{\pm }}=\{\begin{array}{c}<0,\tau ={\tau }_{-}\\ >0,\tau ={\tau }_{+}\end{array}$

(24)

Further, when ${\tau }_{-}={\tau }_{+}$ , it yields

$k_1=2j-1,k_2=2j+1j=1,2,\cdots$

(25)accordingly

${\tau }_c=\frac{k_1\pi }{k_1-k_2}\sqrt{\frac{k_2}{k_1}}$

(26)

Therefore, when ${\overline{\omega }}_0=1$ , $\beta =0$ , there exists a double Hopf bifurcation with $k_1:k_2$ resonance of system (1), and the bifurcation point is

$({\alpha }_c,{\tau }_c)=\left(\frac{k_1-k_2}{k_2}\sqrt{\frac{k_2}{k_1}},\frac{k_1\pi }{k_1-k_2}\sqrt{\frac{k_2}{k_1}}\right)$

(27)where $k_1=2j-1$ , $k_2=2j+1$ , $j=1,2,\cdots$ . □

Based on Conclusion 2, we can obtain some double Hopf bifurcation points of system (1), as shown in Figure 1.

Figure 1.

The double Hopf bifurcation points for different $k_1:k_2$ .

MSM and normal form equations

In this section, we obtain the normal form equations by the MSM and analyzed it. At the same time, the dynamic behaviors around 3:5 weak resonance double Hopf bifurcation point is studied.

Based on conclusion 2, it is known that when ${\overline{\omega }}_0=1$ , $\beta =0$ , there exists a double Hopf bifurcation with 3:5 weak resonance of system (1), and the bifurcation point is

$\left({\alpha }_c,{\tau }_c\right)=\left(-\frac{2}{\sqrt{15}},\frac{\pi \sqrt{15}}{2}\right)$

(28)with

${\omega }_{-}=\frac{\sqrt{15}}{5},{\omega }_{+}=\frac{\sqrt{15}}{3}$

(29)

For the convenience of calculation, set $\gamma =1$ .

According to MSM, the solution of equation (2) can be expressed as below:

$y_1(t)=\epsilon y_{11}+{\epsilon }^2{y}_{12}+{\epsilon }^3{y}_{13}+o({\epsilon }^4)$

(30a)

$y_2(t)=\epsilon y_{21}+{\epsilon }^2{y}_{22}+{\epsilon }^3{y}_{23}+o({\epsilon }^4)$

(30b)

$y_3(t)=\epsilon y_{31}+{\epsilon }^2{y}_{32}+{\epsilon }^3{y}_{33}+o({\epsilon }^4)$

(30c)

$y_4(t)=\epsilon y_{41}+{\epsilon }^2{y}_{42}+{\epsilon }^3{y}_{43}+o({\epsilon }^4)$

(30d)where $y_{ij}=y_{ij}(T_0,T_2),$ $T_0=t,$ $T_2={\epsilon }^{2}t,$ $0<\epsilon \ll 1.$

While the bifurcation parameters are ordered as

$\alpha ={\alpha }_c+{\epsilon }^2{\alpha }_{\epsilon },\tau ={\tau }_c+{\epsilon }^2{\tau }_{\epsilon }$

(31)

The delay term in equation (2) can be further Taylor expanded as

$\begin{array}{c}{y}_2(t-\tau )=\epsilon y_{21}(t-\tau ,{\epsilon }^2(t-\tau ))+{\epsilon }^2{y}_{22}(t-\tau ,{\epsilon }^2(t-\tau ))+{\epsilon }^3{y}_{23}(t-\tau ,{\epsilon }^2(t-\tau ))+o({\epsilon }^4)\\ =\epsilon y_{21}(T_0-{\tau }_c,T_2)+{\epsilon }^2{y}_{22}(T_0-{\tau }_c,T_2)+{\epsilon }^3[-{\tau }_{\epsilon }{D}_0{y}_{21\tau }-{\tau }_c{D}_2{y}_{21\tau }+y_{23\tau }]+o({\epsilon }^4)\end{array}$

(32a)

$\begin{array}{c}{y}_4(t-\tau )=\epsilon y_{41}(t-\tau ,{\epsilon }^2(t-\tau ))+{\epsilon }^2{y}_{42}(t-\tau ,{\epsilon }^2(t-\tau ))+{\epsilon }^3{y}_{43}(t-\tau ,{\epsilon }^2(t-\tau ))+o({\epsilon }^4)\\ =\epsilon y_{41}(T_0-{\tau }_c,T_2)+{\epsilon }^2{y}_{42}(T_0-{\tau }_c,T_2)+{\epsilon }^3[-{\tau }_{\epsilon }{D}_0{y}_{41\tau }-{\tau }_c{D}_2{y}_{41\tau }+y_{43\tau }]+o({\epsilon }^4)\end{array}$

(32b)where $y_{ij\tau }=y_{ij}(T_0-{\tau }_c,T_2)$ , $D_0=\partial /\partial T_0,$ $D_2=\partial /\partial T_2.$

By substituting equations (30) –(32) into equation (2), Taylor expanding the new formula as well and equating separately coefficients of the same powers of $\epsilon$ , the following perturbative equations are obtained:

$\hspace{1em}\epsilon \to$

$D_0{y}_{11}=y_{21}$

(33a)

$D_0{y}_{21}=-y_{11}-\frac{2\sqrt{15}}{15}{y}_{41\tau }$

(33b)

$D_0{y}_{31}=y_{41}$

(33c)

$D_0{y}_{41}=-y_{31}-\frac{2\sqrt{15}}{15}{y}_{21\tau }$

(33b)

$\hspace{1em}{\epsilon }^2\to$

$D_0{y}_{12}=y_{22}$

(34a)

$D_0{y}_{22}=-y_{12}-\frac{2\sqrt{15}}{15}{y}_{42\tau }$

(34b)

$D_0{y}_{32}=y_{42}$

(34c)

$D_0{y}_{42}=-y_{32}-\frac{2\sqrt{15}}{15}{y}_{22\tau }$

(34d)

$\hspace{1em}{\epsilon }^3\to$

$D_0{y}_{13}+D_2{y}_{11}=y_{23}$

(35a)

$D_0{y}_{23}+D_2{y}_{21}=-y_{13}-\frac{2\sqrt{15}}{15}\left[-{\tau }_{\epsilon }{D}_0\left(y_{41\tau }\right)-\frac{\sqrt{15}}{2}\pi D_2\left(y_{41\tau }\right)+y_{43\tau }\right]-{y_{11}}^2{y}_{21}+{\alpha }_{\epsilon }{y}_{41\tau }$

(35b)

$D_0{y}_{33}+D_2{y}_{31}=y_{43}$

(35c)

$D_0{y}_{43}+D_2{y}_{41}=-y_{33}-\frac{2\sqrt{15}}{15}\left[-{\tau }_{\epsilon }{D}_0\left(y_{21\tau }\right)-\frac{\sqrt{15}}{2}\pi D_2\left(y_{21\tau }\right)+y_{23\tau }\right]-{y_{31}}^2{y}_{41}+{\alpha }_{\epsilon }{y}_{21\tau }$

(35d)

Equation (33) has the following general solution:

$y_{11}={\text{Α}}_1{e}^{i{\omega }_1{T}_0}{+}_{\text{Α2}}{e}^{i{\omega }_2{T}_0}+{\overline{A}}_1{e}^{-i{\omega }_1{T}_0}+{\overline{\text{A}}}_2{e}^{-i{\omega }_2{T}_0}$

(36a)

$y_{21}={\text{A}}_{1}i\frac{\sqrt{15}}{5}{e}^{i{\omega }_1{T}_0}+{\text{A}}_{2}i\frac{\sqrt{15}}{3}{e}^{i{\omega }_2{T}_0}+{\overline{\text{A}}}_{1}i\frac{\sqrt{15}}{5}{e}^{-i{\omega }_1{T}_0}+{\overline{\text{A}}}_{2}i\frac{\sqrt{15}}{3}{e}^{-i{\omega }_2{T}_0}$

(36b)

$y_{31}=-{\text{A}}_1{e}^{i{\omega }_1{T}_0}-{\text{A}}_2{e}^{i{\omega }_2{T}_0}-{\overline{\text{A}}}_1{e}^{-i{\omega }_1{T}_0}-{\overline{\text{A}}}_2{e}^{-i{\omega }_2{T}_0}$

(36c)

$y_{41}=-{\text{A}}_{1}i\frac{\sqrt{15}}{5}\cdot e^{i{\omega }_1{T}_0}-{\text{A}}_{2}i\frac{\sqrt{15}}{3}\cdot e^{i{\omega }_2{T}_0}-{\overline{\text{A}}}_{1}i\frac{\sqrt{15}}{5}{e}^{-i{\omega }_1{T}_0}-{\overline{\text{A}}}_{2}i\frac{\sqrt{15}}{3}{e}^{-i{\omega }_2{T}_0}$

(36c)where ${\text{A}}_j={\text{A}}_j(T_2),j=1,2$ are complex constants with respect to $T_2$ .

Substituting equations (36) into equations (34) and (35), and eliminating the secular terms $\epsilon$ , we can obtain the equations including ${\text{A}}_1=D_2{\text{A}}_1$ and ${\text{A}}_2=D_2{\text{A}}_2$ :

${\dot{\text{A}}}_1=C_1{\text{A}}_1+C_{11\overline{1}}{\text{A}}_1^2{\overline{\text{A}}}_1+C_{12\overline{2}}{\text{A}}_1{\text{A}}_2{\overline{\text{A}}}_2$

(37a)

${\dot{\text{A}}}_2=C_2{\text{A}}_2+C_{1\overline{1}2}{\text{A}}_1{\overline{\text{A}}}_{\text{1}}{\text{A}}_2+C_{22\overline{2}}{\text{A}}_2^2{\overline{\text{A}}}_2$

(37b)where the complex coefficients $C_l,C_{ijk}$ are expressed as follow:

$C_1=-0.\text{477692}+0.\text{369119}i-\text{(0.185009}+0.\text{157041}i\text{)}\alpha +\text{(0.0628164}-0.0\text{740038}i\text{)}\tau$

(38a)

$C_2=\text{0.39156}+\text{1}.\text{09157}i-\text{(0.252751}-0.\text{1128725}i\text{)}\alpha -\text{(0.0858166}+0.\text{168501}i\text{)}\tau$

(38b)

$C_{11\overline{1}}=-\text{0.157041}+\text{0.185009}i$

(38c)

$C_{12\overline{2}}=-\text{0.314082}+\text{0.370019}i$

(38d)

$C_{1\overline{1}2}=-\text{0.25745}-\text{0.505502}i$

(38e)

$C_{22\overline{2}}=-\text{0.128725}-\text{0.252751}i$

(38f)

For expressing the complex amplitude equations in real form, a mixed form representation of the complex amplitudes is introduced:

${\text{A}}_1=\frac{1}{2}{r}_1{e}^{i{\theta }_1},{\text{A}}_2=\frac{1}{2}{r}_2{e}^{i{\theta }_2}$

(39)

Substituting equation (39) into equation (37) and separating the real and imaginary parts, the normal form equations are drawn:

${\dot{r}}_1=(-0.\text{477692}-0.\text{185009\hspace{0.17em}}\alpha +0.0\text{628164}\tau )r_1-\text{0.0392602}{r}_1^3-\text{0.0785204}{r}_1{r}_2^2$

(40a)

${\dot{r}}_2=\text{(0.39156}-0.\text{252751}\alpha -0.0\text{858166}\tau )r_2-\text{0.0643625}{r}_1^2{r}_2-\text{0.0321812}{r}_2^3$

(40b)

Let ${\dot{r}}_1=0,{\dot{r}}_2=0$ , then we can obtain all possible equilibrium points of equation (40):

$E_0=(0,0),E_1=(r_{\text{11}},0),E_2=(0,r_{\text{21}}),E_3=(r_{\text{12}},r_{\text{22}})$

(41)where

$r_{\text{11}}=\sqrt{-\text{12.1673}-\text{4.71239}\alpha \text{+1.6000}\tau }$

(42a)

$r_{\text{21}}=\sqrt{\text{12.1673}-\text{7.85398}\alpha -\text{2.66667}\tau }$

(42b)

$r_{\text{12}}=\sqrt{\text{12.1673}-\text{3.66519}\alpha -\text{2.31111}\tau }$

(42c)

$r_{\text{22}}=\sqrt{-\text{12.1673}-0.\text{523599}\alpha +\text{1.95556}\tau }$

(42d)

The existence of the above possible equilibrium points is dependent on $\alpha$ and $\tau$ . And according to the stability analysis of existing equilibrium points, the parameters $(\alpha ,\tau )$ plane around the 3:5 weakly resonant double Hopf bifurcation point $({\alpha }_c,{\tau }_c)$ can be divided into six distinct dynamical regions as shown in Figures 2 and 3. The solutions in the same region have the same topological structure.

$\text{I}=\left\{(\alpha ,\tau )|r_{\text{11}}<0,r_{\text{21}}<0,r_{\text{12}}{r}_{\text{22}}<0\right\}$

(43a)

$\text{II}=\left\{(\alpha ,\tau )|r_{\text{11}}>0,r_{\text{21}}<0,r_{\text{12}}>0,r_{\text{22}}<0\right\}$

(43b)

$\text{III}=\left\{(\alpha ,\tau )|r_{\text{11}}>0,r_{\text{21}}>0,r_{\text{12}}>0,r_{\text{22}}<0\right\}$

(43c)

$\text{IV}=\left\{(\alpha ,\tau )|r_{\text{11}}>0,r_{\text{21}}>0,r_{\text{12}}>0,r_{\text{22}}>0\right\}$

(43d)

$\text{V}=\left\{(\alpha ,\tau )|r_{\text{11}}>0,r_{\text{21}}>0,r_{\text{12}}<0,r_{\text{22}}>0\right\}$

(43e)

$\text{VI}=\left\{(\alpha ,\tau )|r_{\text{11}}<0,r_{\text{21}}>0,r_{\text{12}}<0,r_{\text{22}}>0\right\}$

(43f)

Figure 2.

The dynamical behaviors around the 3:5 weakly resonant double Hopf bifurcation point $({\alpha }_c,{\tau }_c)$ .

Figure 3.

Classification of intrinsic dynamic behavior around the 3:5 weakly resonant double Hopf bifurcation point $({\alpha }_c,{\tau }_c)$ .

In region I, there is only a trivial equilibrium point $E_0=(0,0)$ , which is stable. That means the vibration of system (1) is restricted by parameters $\alpha$ and $\tau$ , so region I is usually called as an amplitude death region.

With $(\alpha ,\tau )$ changing into region II, the trivial equilibrium point $E_0=(0,0)$ loses its stability and a stable equilibrium point $E_1=(r_{\text{11}},0)$ appears. Correspondingly, there is a stable periodic solution with a frequency close to ${\omega }_{+}\approx 1.29099$ .

When $(\alpha ,\tau )$ enters into region III, the unstable equilibrium point $E_0=(0,0)$ and the stable equilibrium point $E_1=(r_{\text{11}},0)$ still exist, and an unstable equilibrium point $E_2=(0,r_{\text{21}})$ appears. Correspondingly, there exists a stable periodic solution with a frequency close to ${\omega }_{+}$ and an unstable periodic solution.

In region IV, an unstable nontrivial equilibrium point $E_3=(r_{\text{12}},r_{\text{22}})$ appears together with two stable equilibrium points $E_1=(r_{\text{11}},0)$ and $E_2=(0,r_{\text{21}})$ . Correspondingly, system (1) admits two coexisting periodic oscillations and an unstable quasiperiodic motion.

With $(\alpha ,\tau )$ changing into region V, the equilibrium point $E_3=(r_{\text{12}},r_{\text{22}})$ disappears and there exist three equilibrium points, $E_0=(0,0)$ , $E_1=(r_{\text{11}},0)$ and $E_2=(0,r_{\text{21}})$ . In contrast with the situation in region III, the equilibrium point $E_1=(r_{\text{11}},0)$ is unstable and $E_2=(0,r_{\text{21}})$ is stable. Correspondingly, there exists a stable periodic solution with a frequency close to ${\omega }_{-}$ and an unstable periodic solution.

When $(\alpha ,\tau )$ enters into region VI, the equilibrium point $E_1=(r_{\text{11}},0)$ disappears and there exists an unstable equilibrium point $E_0=(0,0)$ and a stable equilibrium point $E_2=(0,r_{\text{21}})$ . Correspondingly, there is a stable periodic solution with a frequency close to ${\omega }_{-}\approx 0.774597$ .

To demonstrate the above theoretical results, we choose six groups of parameter values $(\alpha ,\tau )$ for system (1) with ${\overline{\omega }}_0=1$ , $\beta =0$ , and the results are given in Figure 4. It is obvious that numerical simulations are well consistent with the theoretical results.

Figure 4.

The phase portraits when $(\alpha ,\tau )$ is located in each region.

$\begin{array}{l}{p}_1:(\alpha ,\tau )=(-0.46,6)\in \text{I},p_2:(\alpha ,\tau )=(-0.5,7)\in \text{II}\\ p_3:(\alpha ,\tau )=(-0.7,6.5)\in \text{III},p_4:(\alpha ,\tau )=(-0.8,6.2)\in \text{IV}\\ p_5:(\alpha ,\tau )=(-0.7,5.8)\in \text{V,\hspace{0.17em}\hspace{0.17em}}{p}_6:(\alpha ,\tau )=(-0.5,5)\in \text{VI}\end{array}$

Application of HAM

Assuming system (1) satisfies the initial condition:

$x_1(t)=A(t),{\dot{x}}_1(t)=B(t),x_2(t)=C,{\dot{x}}_2(t)=0,-\tau \le t\le 0$

(44)and the frequency of initial value solution is $\omega$ .

Let $\theta =\omega t,x_1(t)=u(\theta ),x_2(t)=v(\theta )$ , then system (1) can be expressed as

${\omega }^2{u}^{″}(\theta )+\omega \left(\gamma u^2(\theta )-\beta \right)u^{\prime }(\theta )+{\overline{\omega }}_0^{2}u(\theta )=\alpha \omega v^{\prime }(\theta -\omega \tau )$

(45a)

${\omega }^2{v}^{″}(\theta )+\omega \left(\gamma v^2(\theta )-\beta \right)v^{\prime }(\theta )+{\overline{\omega }}_0^{2}v(\theta )=\alpha \omega u^{\prime }(\theta -\omega \tau )$

(45b)with the initial condition

$u(\theta )=A(\theta ),u^{\prime }(\theta )=B(\theta ),v(\theta )=C,v^{\prime }(\theta )=0,-\omega \tau \le \theta \le 0$

(46)

Suppose that the solutions of equations (45) and (46) can be expressed by the set of base functions

$\{\text{cos}(2k+1)\theta ,\text{sin}(2k+1)\theta |k=0,1,2,\cdots \}$

(47)

The initial guess solution is chosen as

$u_0(\theta )=a_0\text{cos}\theta +b_0\text{sin}\theta ,v_0(\theta )=c_0\text{cos}\theta$

(48)and the linear operator is expressed as

$L\left(\begin{array}{c}U(\theta ,q)\\ V(\theta ,q)\end{array}\right)=\left(\begin{array}{c}\frac{{\partial }^{2}U(\theta ,q)}{\partial {\theta }^2}+U(\theta ,q)\\ \frac{{\partial }^{2}V(\theta ,q)}{\partial {\theta }^2}+V(\theta ,q)\end{array}\right)$

(49)with the property

$L\left(\begin{array}{c}{C}_1\text{cos}\theta +C_2\text{sin}\theta \\ C_3\text{cos}\theta +C_4\text{sin}\theta \end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)$

(50)where $C_1$ , $C_2$ , $C_3$ _, and $C_4$ are constants.

From equation (45), we define a nonlinear operator

$\begin{array}{c}N\left(U(\theta ,q)\right)=W^2(q)\frac{{\partial }^{2}U(\theta ,q)}{\partial {\theta }^2}+W(q)\left(\gamma U{(\theta ,q)}^2-\beta \right)\frac{\partial U(\theta ,q)}{\partial \theta }\\ +{\overline{\omega }}_0^{2}U(\theta ,q)-\alpha W(q)\frac{\partial V\left(\theta -W(q)\tau ,q\right)}{\partial \theta }\end{array}$

(51a)

$\begin{array}{c}N\left(V(\theta ,q)\right)=W^2(q)\frac{{\partial }^{2}V(\theta ,q)}{\partial {\theta }^2}+W(q)\left(\gamma V{(\theta ,q)}^2-\beta \right)\frac{\partial V(\theta ,q)}{\partial \theta }\\ +{\overline{\omega }}_0^{2}V(\theta ,q)-\alpha W(q)\frac{\partial U\left(\theta -W(q)\tau ,q\right)}{\partial \theta }\end{array}$

(51b)

Therefore, the zeroth-order deformation equation can be written as

$(1-q)L\left(\begin{array}{c}U(\theta ,q)-u_0(\theta )\\ V(\theta ,q)-v_0(\theta )\end{array}\right)=q\hslash N\left(\begin{array}{c}U(\theta ,q)\\ V(\theta ,q)\end{array}\right)$

(52)with the initial conditions

$U(\theta ,q)=a(q),\frac{\partial U(\theta ,q)}{\partial \theta }=b(q),V(\theta ,q)=c(q),\frac{\partial V(\theta ,q)}{\partial \theta }=0,-\omega \tau \le \theta \le 0$

(53)

For $q=0$ , the solution of equations (52) and (53) is

$U(\theta ,0)=u_0(\theta ),V(\theta ,0)=v_0(\theta ),W(0)={\omega }_0$

(54)

While for $q=1$ , equations (52) and (53) are equivalent to the original equations, and admit the solution

$U(\theta ,1)=u(\theta ),V(\theta ,1)=v(\theta ),W(1)=\omega$

(55)

Hence, as $q$ changes from 0 to 1, $U(\theta ,q)$ and $V(\theta ,q)$ , respectively vary from the initial guess $u_0(\theta )$ and $v_0(\theta )$ to the unknown $u(\theta )$ and $v(\theta )$ . Likewise, $W(q)$ varies from the initial guess frequency ${\omega }_0$ to the physical frequency $\omega$ .

According to Taylor series expansion, we obtain

$U(\theta ,q)=u_0(\theta )+{\displaystyle \sum _{m=1}^{+\infty }{u}_m(\theta )}{q}^m$

(56a)

$V(\theta ,q)=v_0(\theta )+{\displaystyle \sum _{m=1}^{+\infty }{v}_m(\theta )}{q}^m$

(56b)

$W(q)={\omega }_0+{\displaystyle \sum _{m=1}^{+\infty }{\omega }_m}{q}^m$

(56c)where

${u_m(\theta )=\frac{1}{m!}\frac{{\partial }^{m}U(\theta ,q)}{\partial q^m}|}_{q=0}$

(57a)

${v_m(\theta )=\frac{1}{m!}\frac{{\partial }^{m}V(\theta ,q)}{\partial q^m}|}_{q=0}$

(57b)

${{\omega }_m(\theta )=\frac{1}{m!}\frac{{\partial }^{m}W(q)}{\partial q^m}|}_{q=0}$

(57c)

Note that the zero-order deformation equations have auxiliary parameters. Assuming $\hslash$ is so properly chosen that power series equation (56) converge at $q=1$ , the series solutions of equation (56) can be expressed as

$u(\theta )=u_0(\theta )+{\displaystyle \sum _{m=1}^{+\infty }{u}_m(\theta )}$

(58a)

$v(\theta )=v_0(\theta )+{\displaystyle \sum _{m=1}^{+\infty }{v}_m(\theta )}$

(58b)

$\omega ={\omega }_0+{\displaystyle \sum _{m=1}^{+\infty }{\omega }_m}$

(58c)

For the sake of simplicity, we define the vectors

${\overrightarrow{u}}_m=\left\{u_0(\theta ),u_1(\theta ),\cdots ,u_m(\theta )\right\}$

(59a)

${\overrightarrow{v}}_m=\left\{v_0(\theta ),v_1(\theta ),\cdots ,v_m(\theta )\right\}$

(59b)

${\overrightarrow{\omega }}_m=\left\{{\omega }_0,{\omega }_1,\cdots ,{\omega }_m\right\}$

(59c)

By differentiating the zeroth-order deformation equation (52) $m$ $(m\ge 1,m\in N_{+})$ times with respect to $q$ , then dividing the equation by $m!$ and setting $q=0$ , we arrive at the mth-order deformation equation

$L\left(\begin{array}{c}{u}_m(\theta )-{\chi }_m{u}_{m-1}(\theta )\\ v_m(\theta )-{\chi }_m{v}_{m-1}(\theta )\end{array}\right)=\hslash \left(\begin{array}{c}{R}_{1,m}({\overrightarrow{u}}_{m-1},{\overrightarrow{\omega }}_{m-1})\\ R_{2,m}({\overrightarrow{v}}_{m-1},{\overrightarrow{\omega }}_{m-1})\end{array}\right)$

(60)where

${\chi }_m=\{\begin{array}{cc}0,& m=1\\ 1,& m\ge 2\end{array}$

(61)

$\begin{array}{l}{R}_{1,m}({\overrightarrow{u}}_{m-1},{\overrightarrow{\omega }}_{m-1})={\frac{1}{(m-1)!}\frac{{\partial }^{m-1}N\left(U(\theta ,q)\right)}{\partial q^{m-1}}|}_{q=0}\\ ={\displaystyle \sum _{i=0}^{m-1}{u^{″}}_{m-1-i}(\theta )\left({\displaystyle \sum _{j=0}^i{\omega }_j{\omega }_{i-j}}\right)}+\gamma {\displaystyle \sum _{i=0}^{m-1}{u^{\prime }}_{m-1-i}(\theta )\left({\displaystyle \sum _{j=0}^i{\omega }_{i-j}}\left({\displaystyle \sum _{k=0}^j{u}_k{u}_{j-k}}\right)\right)}\\ -\beta {\displaystyle \sum _{i=0}^{m-1}{\omega }_{m-1-i}{u}^{\prime }(\theta )}+{\overline{\omega }}_0^2{u}_{m-1}(\theta )-\alpha P_{m-1}\end{array}$

(62)

$\begin{array}{l}{R}_{2,m}({\overrightarrow{v}}_{m-1},{\overrightarrow{\omega }}_{m-1})={\frac{1}{(m-1)!}\frac{{\partial }^{m-1}N\left(V(\theta ,q)\right)}{\partial q^{m-1}}|}_{q=0}\\ ={\displaystyle \sum _{i=0}^{m-1}{v^{″}}_{m-1-i}(\theta )\left({\displaystyle \sum _{j=0}^i{\omega }_j{\omega }_{i-j}}\right)}+\gamma {\displaystyle \sum _{i=0}^{m-1}{v^{\prime }}_{m-1-i}(\theta )\left({\displaystyle \sum _{j=0}^i{\omega }_{i-j}}\left({\displaystyle \sum _{k=0}^j{v}_k{v}_{j-k}}\right)\right)}\\ -\beta {\displaystyle \sum _{i=0}^{m-1}{\omega }_{m-1-i}{v^{\prime }}_i}(\theta )+{\overline{\omega }}_0^2{v}_{m-1}(\theta )-\alpha Q_{m-1}\end{array}$

(63)

$\begin{array}{l}{P}_{m-1}={\frac{1}{(m-1)!}\frac{{\partial }^{m-1}\left(W(q)\frac{\partial V\left(\theta -W(q)\tau ,q\right)}{\partial \theta }\right)}{\partial q^{m-1}}|}_{q=0}\\ ={\displaystyle \sum _{i=0}^{m-1}{\omega }_i\left[{v^{\prime }}_{m-1-i}(\theta -{\omega }_0\tau )+{\displaystyle \sum _{j=1}^{m-1-i}\frac{{(-\tau )}^j}{j!}}\left({\displaystyle \sum _{k=0}^{m-1-i-j}{D}_{m-1-i-k,j}}\frac{{\partial }^j{v^{\prime }}_k(\theta -{\omega }_0\tau )}{\partial {\theta }^j}\right)\right]}\end{array}$

(64)

$\begin{array}{l}{Q}_{m-1}={\frac{1}{(m-1)!}\frac{{\partial }^{m-1}\left(W(q)\frac{\partial U\left(\theta -W(q)\tau ,q\right)}{\partial \theta }\right)}{\partial q^{m-1}}|}_{q=0}\\ ={\displaystyle \sum _{i=0}^{m-1}{\omega }_i\left[{u^{\prime }}_{m-1-i}(\theta -{\omega }_0\tau )+{\displaystyle \sum _{j=1}^{m-1-i}\frac{{(-\tau )}^j}{j!}}\left({\displaystyle \sum _{k=0}^{m-1-i-j}{D}_{m-1-i-k,j}}\frac{{\partial }^j{u^{\prime }}_k(\theta -{\omega }_0\tau )}{\partial {\theta }^j}\right)\right]}\end{array}$

(65)

$D_{n,1}={\omega }_n,D_{n,l}={\displaystyle \sum _{j=1}^{n-1}{\omega }_j}{B}_{n-j,l-1},\hspace{1em}l\ge 2$

(66)

and the initial conditions

$u_m(\theta )=a_m,{u^{\prime }}_m(\theta )=b_m,v_m(\theta )=c_m,{v^{\prime }}_m(\theta )=0,-\omega \tau \le \theta \le 0$

(67)

According to the rule of solution expression, $u_m$ and $v_m$ should not contain the so-called secular terms of $\theta \text{sin}\theta$ and $\theta \text{cos}\theta$ . Correspondingly, $\hslash R_{1,m}({\overrightarrow{u}}_{m-1},{\overrightarrow{\omega }}_{m-1})$ and $\hslash R_{2,m}({\overrightarrow{u}}_{m-1},{\overrightarrow{\omega }}_{m-1})$ should not contain the terms of $\text{sin}\theta$ and $\text{cos}\theta$ . Setting their coefficients to zeroes, we obtains

$\frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }\hslash R_{1,m}({\overrightarrow{u}}_{m-1},{\overrightarrow{\omega }}_{m-1})\text{cos}\theta }\text{d}\theta =0$

(68a)

$\frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }\hslash R_{1,m}({\overrightarrow{u}}_{m-1},{\overrightarrow{\omega }}_{m-1})\text{sin}\theta }\text{d}\theta =0$

(68b)

$\frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }\hslash R_{2,m}({\overrightarrow{u}}_{m-1},{\overrightarrow{\omega }}_{m-1})\text{cos}\theta }\text{d}\theta =0$

(68c)

$\frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }\hslash R_{2,m}({\overrightarrow{u}}_{m-1},{\overrightarrow{\omega }}_{m-1})\text{sin}\theta }\text{d}\theta =0$

(68d)

Based on equation (68), $a_{m-1},$ $b_{m-1},$ $c_{m-1}$ and ${\omega }_{m-1}$ can be obtained, for each $m\ge 1.$

Hence, $\hslash R_{1,m}({\overrightarrow{u}}_{m-1},{\overrightarrow{\omega }}_{m-1})$ and $\hslash R_{2,m}({\overrightarrow{u}}_{m-1},{\overrightarrow{\omega }}_{m-1})$ can be expressed as

$\hslash R_{1,m}({\overrightarrow{u}}_{m-1},{\overrightarrow{\omega }}_{m-1})={\displaystyle \sum _{n=1}{a}_{1,n}\text{cos}(2n+1)\theta }+{\displaystyle \sum _{n=1}{a}_{2,n}\text{sin}(2n+1)\theta }$

(69a)

$\hslash R_{2,m}({\overrightarrow{u}}_{m-1},{\overrightarrow{\omega }}_{m-1})={\displaystyle \sum _{n=1}{b}_{1,n}\text{cos}(2n+1)\theta }+{\displaystyle \sum _{n=1}{b}_{2,n}\text{sin}(2n+1)\theta }$

(69b)

Further, by solving equation (60), it yields

$\begin{array}{l}{u}_m(\theta )={\chi }_m{u}_{m-1}(\theta )-{\displaystyle \sum _{n=1}\frac{a_{1,n}}{4n(n+1)}\text{cos}(2n+1)\theta }\\ +{\displaystyle \sum _{n=1}\frac{a_{2,n}}{4n(n+1)}\text{sin}(2n+1)\theta }+C_{1,m}\text{sin}\theta +C_{2,m}\text{cos}\theta \end{array}$

(70a)

$\begin{array}{l}{v}_m(\theta )={\chi }_m{v}_{m-1}(\theta )-{\displaystyle \sum _{n=1}\frac{b_{1,n}}{4n(n+1)}\text{cos}(2n+1)\theta }\\ +{\displaystyle \sum _{n=1}\frac{b_{2,n}}{4n(n+1)}\text{sin}(2n+1)\theta }+C_{3,m}\text{sin}\theta +C_{4,m}\text{cos}\theta \end{array}$

(70b)

According to equation (68), $C_{1,m}$ , $C_{2,m}$ , $C_{3,m}$ , and $C_{4,m}$ can be expressed with respect to $a_m$ , $b_m$ , and $c_m$ . So, as $a_m$ , $b_m$ , and $c_m$ are obtained, $C_{1,m}$ , $C_{2,m}$ , $C_{3,m}$ , and $C_{4,m}$ are obtained successively, and then $u_m(\theta )$ and $v_m(\theta )$ are obtained.

Therefore, through the iteration step by step, the nth-order homotopy analytical approximation can be obtained, as expressed in terms of the summation series

$u_{[n]}(\theta )=u_0(\theta )+{\displaystyle \sum _{m=1}^n{u}_m(\theta )}$

(71a)

$v_{[n]}(\theta )=v_0(\theta )+{\displaystyle \sum _{m=1}^n{v}_m(\theta )}$

(71b)

${\omega }_{[n]}={\omega }_0+{\displaystyle \sum _{m=1}^n{\omega }_m}$

(71c)

Numerical simulation

According to the MSM and normal form equations section, parameters $(\alpha ,\tau )$ plane around the 3:5 weakly resonant double Hopf bifurcation point $({\alpha }_c,{\tau }_c)$ can be divided into six distinct dynamical regions, and when $(\alpha ,\tau )$ is located in different regions, the system (1) presents different dynamic behaviors. Now, HAM is used to obtain analytical approximate solutions of system (1) with $(\alpha ,\tau )$ located in each region. Since region I is an amplitude death region, it is not necessary to use HAM to approximate its exact solution. Consider that when $(\alpha ,\tau )$ is located in region II or VI, system (1), respectively has a stable periodic solution, but the frequency is different, so we put these two regions together to study. So is for regions III and V. For region IV, system (1) admits two coexisting periodic oscillations and an unstable quasiperiodic motion. Because the base functions of HAM in this paper are periodic functions, we just can obtain the two stable periodic analytical approximate solutions. As for the quasiperiodic analytical approximate solution, it will be the follow-up research.

The periodic solution by HAM with (α,τ) located in regions II or VI

We use the Mathematica software to the study of the subject when $(\alpha ,\tau )$ is located in regions IIand VI. For the region II, letting

${\overline{\omega }}_0=1,\beta =0,\gamma =1,\alpha =-0.5,\tau =7$

(72)

We can get the initial value from equation (68) as follows

$u_0(\theta )=\text{1.108515}\text{cos}\theta ,v_0(\theta )=\text{1.108515}\text{cos}\theta ,{\omega }_0=\text{1.2165158}$

(73)

It means that system (1) has one periodic solution, as is in accord with the theoretical analysis in MSM and normal form equations section.

Let $\hslash =-0.5$ , the fifth-order approximation can be obtained:

$\begin{array}{l}{x}_1(t)=1.09691\text{cos}1.214241t-0.00179299\text{cos}3.642723t-0.00153991\text{cos}6.07120t\\ -0.00000390\text{cos}8.49969t\text{+}0.00000394\text{cos}10.92817t\text{+}0.0982843\text{sin}1.214241t\\ -0.0329664\text{sin}3.64272t\text{+}0.000001792\text{sin}6.07120t\text{+}0.0000866711\text{sin}8.49969t\\ \text{\hspace{1em}+}0.0000000218\text{sin}10.9282t-0.0000000948\text{sin}13.3567t\end{array}$

(74a)

$\begin{array}{l}{x}_2(t)=1.09691\text{cos}1.214241t-0.00179299\text{cos}3.642723t-0.00153991\text{cos}6.07120t\\ -0.00000390\text{cos}8.49969t\text{+}0.00000394\text{cos}10.92817t\text{+}0.0982843\text{sin}1.214241t\\ -0.0329664\text{sin}3.64272t\text{+}0.000001792\text{sin}6.07120t\text{+}0.0000866711\text{sin}8.49969t\\ \text{\hspace{1em}+}0.0000000218\text{sin}10.9282t-0.0000000948\text{sin}13.3567t\end{array}$

(74b)

$\omega =1.21424$

(74c)

The phase portrait curves of the fifth-order HAM approximation are shown in Figure 5, compared with that of fourth-order Runge–Kutta method. We found that the approximate solution agrees with the numerical solution.

Figure 5.

Comparison of the phase portrait curves of the fifth-order HAM approximation with that of fourth-order Runge–Kutta method for ${\overline{\omega }}_0=1$ , $\beta =0$ , $\gamma =1$ , $\alpha =-0.5$ , $\tau =7$ .

When $(\alpha ,\tau )$ is located in region VI, let

${\overline{\omega }}_0=1,\beta =0,\gamma =1,\alpha =-0.5,\tau =5$

(75)

According to equation (68) when $m=1$ , a set of initial guess is obtained

$u_0(\theta )=\text{1.0814565}\text{cos}\theta ,v_0(\theta )=\text{1.0814565}\text{cos}\theta ,{\omega }_0=\text{0.8175575}$

(76)

It means that system (1) has one periodic solution, as is in accord with the theoretical analysis in MSM and normal form equations section.

Let $\hslash =-1.5$ , the fifth-order approximation can be obtained:

$\begin{array}{l}{x}_1(t)=1.09211\text{cos}0.814731t-0.00411785\text{cos}2.44419t-0.00385916\text{cos}4.07365t\\ -0.0000252727\text{cos}5.70311t+0.0000299202\text{cos}7.33257t+0.130287\text{sin}0.814731t\\ -0.0446114\text{sin}2.44419t+0.000296495\text{sin}4.07365t+0.000299236\text{sin}5.70311t\\ -0.000000467\text{sin}7.33257t-\text{0.00000}2405\text{sin}8.96204t\end{array}$

(77a)

$\begin{array}{l}{x}_2(t)=1.09211\text{cos}0.814731t-0.00411785\text{cos}2.44419t-0.00385916\text{cos}4.07365t\\ -0.0000252727\text{cos}5.70311t+0.0000299202\text{cos}7.33257t+0.130287\text{sin}0.814731t\\ -0.0446114\text{sin}2.44419t+0.000296495\text{sin}4.07365t+0.000299236\text{sin}5.70311t\\ -0.000000467\text{sin}7.33257t-\text{0.00000}2405\text{sin}8.96204t\end{array}$

(77b)

$\omega =\text{0.814731}$

(77c)

The phase portrait curves of the fifth-order HAM approximation are shown in Figure 6, compared with that of fourth-order Runge–Kutta method. By observing Figure 6, it is clear that the fifth-order HAM approximation of $x_1(t)$ is much matched to that of fourth-order Runge–Kutta method, and the same conclusion is for $x_2(t)$ .

Figure 6.

Comparison of the phase portrait curves of the fifth-order HAM approximation with that of fourth-order Runge–Kutta method for ${\overline{\omega }}_0=1$ , $\beta =0$ , $\gamma =1$ , $\alpha =-0.5$ , $\tau =5$ .

The periodic solution by HAM with (α,τ) located in regions III or V

We use the Mathematica software to the study of the system when $(\alpha ,\tau )$ is located in regions III and V. For the region III, letting

${\overline{\omega }}_0=1,\beta =0,\gamma =1,\alpha =-\text{0.9,\hspace{1em}}\tau =7$

(78)

According to equation (68) when $m=1$ , two sets of initial guess are obtained

$u_0(\theta )=\text{1.74718}\text{cos}\theta ,v_0(\theta )=\text{1.74718}\text{cos}\theta ,{\omega }_0=\text{1.26658535}$

(79a)

$u_0(\theta )=\text{0.761786}\text{cos}\theta ,v_0(\theta )=\text{0.761786}\text{cos}\theta ,{\omega }_0=\text{0.650069}$

(79b)

It means that system (1) has two periodic solutions, as is in accord with the theoretical analysis in MSM and normal form equations section.

Select equation (79a) as the initial guess, and let $\hslash =-0.2$ , then the fifth-order approximation can be obtained:

$\begin{array}{l}{x}_1(t)=1.70033\text{cos}1.25801t+0.012356\text{cos}3.77404t-0.0081641\text{cos}6.29007t\\ -0.0000341501\text{cos}.8061t+8.667\times {10}^{-6}\text{cos}11.3221t+0.421051\text{sin}1.25801t\\ -0.140031\text{sin}3.77404t-0.000569125\text{sin}6.29007+0.000269387\text{sin}8.8061t\\ +3.97\times {10}^{-7}\text{sin}11.3221t-8.18\times {10}^{-8}\text{sin}13.8382t\end{array}$

(80a)

$\begin{array}{l}{x}_2\text{(}t\text{)}=1.70033\text{cos}1.25801t+0.012356\text{cos}3.77404t-0.0081641\text{cos}6.29007t\\ -0.0000341501\text{cos}.8061t+8.667\times {10}^{-6}\text{cos}11.3221t+0.421051\text{sin}1.25801t\\ -0.140031\text{sin}3.77404t-0.000569125\text{sin}6.29007+0.000269387\text{sin}8.8061t\\ +3.97\times {10}^{-7}\text{sin}11.3221t-8.18\times {10}^{-8}\text{sin}13.8382t\end{array}$

(80b)

$\omega =1.25801$

(80c)

Select equation (79b) as the initial guess, and let $\hslash =-1.9$ , then the fifth-order approximation can be obtained:

$\begin{array}{l}{x}_1(t)=0.795013\text{cos}0.648237t-0.0120303\text{cos}1.94471t-0.00185005\text{cos}3.24118t\\ +\text{0.0000139149cos4.53766}t+3.24\times {10}^{-6}\text{cos}5.83413t+\text{0.136882sin0.64827}t\\ -\text{0.0456724sin1.94471}t-\text{0.000114675sin3.24118}t+\text{0.000101239sin4.53766}t\\ +3.4967\times {10}^{-8}\text{sin}5.83413t-5.28\times {10}^{-8}\text{sin}7.1306t\end{array}$

(81a)

$\begin{array}{l}{x}_2(t)=\text{0.795013}\text{cos}\text{0.648237}t-\text{0.0120303cos1.94471}t-0.00185005\text{cos}3.24118t\\ +\text{0.0000139149cos4.53766}t+3.24\times {10}^{-6}\text{cos}5.83413t+\text{0.136882sin0.64827}t\\ -\text{0.0456724sin1.94471}t-\text{0.000114675sin3.24118}t+\text{0.000101239sin4.53766}t\\ +3.4967\times {10}^{-8}\text{sin}5.83413t-5.28\times {10}^{-8}\text{sin}7.1306t\end{array}$

(81b)

$\omega =\text{0.648237}$

(81c)

The phase portrait curves of the fifth-order HAM approximation are shown in Figure 7, compared with that of fourth-order Runge–Kutta method. Figure 7 shows that the Runge–Kutta method cannot get the unstable periodic solution, the HAM can get a stable periodic solution, but also can be unstable periodic solution. As discussed in Ebenezer et al. ³⁹ the multistep HAM allows us with great freedom to choose the convergence–control parameter ħ. Although it is plausible to optimize the HAM in other ways, the use of the convergence–control parameter ħ is the most accessible way to achieve faster convergent homotopy-series solutions. The present manuscript proposes the use of the convergence–control parameter ħ to refine the accuracy of solution series.

Figure 7.

Comparison of the phase portrait curves of the fifth-order HAM approximation with that of fourth-order Runge–Kutta method for ${\overline{\omega }}_0=1$ , $\beta =0$ , $\gamma =1$ , $\alpha =-0.9$ , $\tau =7$ .

Similarly, let’s consider the condition that $(\alpha ,\tau )$ is located in region V. Set

${\overline{\omega }}_0=1,\beta =0,\gamma =1,\alpha =-0.8,\tau =5.5$

(82)

According to equation (68) when $m=1$ , two sets of initial guess are obtained

$u_0(\theta )=\text{1.41094}\text{cos}\theta ,v_0(\theta )=\text{1.41094}\text{cos}\theta ,{\omega }_0=\text{0.7347193}$

(83a)

$u_0(\theta )=\text{0.809123}\text{cos}\theta ,v_0(\theta )=\text{0.809123}\text{cos}\theta ,{\omega }_0=\text{1.4654589}$

(83b)

It means that system (1) has two periodic solutions, as is in accord with the theoretical analysis in MSM and normal form equations section.

Select equation (83a) as the initial guess, and let $\hslash =-0.6$ , then the fifth-order approximation can be obtained:

$\begin{array}{l}{x}_1(t)=\text{1.43556cos0.729193}t-\text{0.00632019cos2.18758}t-\text{0.00950247cos3.64597}t\\ -\text{0.0000183465cos5.10435}t+\text{0.0000183063cos6.56274}t+\text{0.278014sin0.729193}t\\ -\text{0.0942229sin2.18758}t+\text{0.000171301sin3.64597}t+\text{0.000543328sin5.10435}t\\ +\text{8.98836}\times {\text{10}}^{\text{-8}}\text{sin6.56274}t-\text{2.69123}\times {\text{10}}^{\text{-7}}\text{sin8.02113}t\end{array}$

(84a)

$\begin{array}{l}{x}_2(t)=\text{1.43556cos0.729193}t-\text{0.00632019cos2.18758}t-\text{0.00950247cos3.64597}t\\ -\text{0.0000183465cos5.10435}t+\text{0.0000183063cos6.56274}t+\text{0.278014sin0.729193}t\\ -\text{0.0942229sin2.18758}t+\text{0.000171301sin3.64597}t+\text{0.000543328sin5.10435}t\\ +\text{8.98836}\times {\text{10}}^{\text{-8}}\text{sin6.56274}t-\text{2.69123}\times {\text{10}}^{\text{-7}}\text{sin8.02113}t\end{array}$

(84b)

$\omega =\text{0.729193}$

(84c)

Select equation (83b) as the initial guess, and let $\hslash =-0.5$ , then the fifth-order approximation can be obtained:

$\begin{array}{l}{x}_1(t)=\text{0.801463cos1.46473}t-\text{0.000496562cos4.39419}t-\text{0.000205669cos7.32364}t\\ +\text{1.9125}\times {\text{10}}^{-\text{8}}\text{cos10.2531}t+\text{3.27256}\times {\text{10}}^{-\text{8}}\text{cos13.1826}t+\text{0.0262201sin1.46473}t\\ -\text{0.0087362sin4.39419}t-\text{9.75925}\times {\text{10}}^{-\text{6}}\text{sin7.32364}t+\text{5.34464}\times {\text{10}}^{-\text{6}}\text{sin10.2531}t\\ +\text{1.38416}\times {\text{10}}^{-\text{9}}\text{sin13.1826}t-\text{7.53144}\times {\text{10}}^{-\text{9}}\text{sin16.112}t\end{array}$

(85a)

$\begin{array}{l}{x}_2(t)=\text{0.801463cos1.46473}t-\text{0.000496562cos4.39419}t-\text{0.000205669cos7.32364}t\\ +\text{1.9125}\times {\text{10}}^{-\text{8}}\text{cos10.2531}t+\text{3.27256}\times {\text{10}}^{-\text{8}}\text{cos13.1826}t+\text{0.0262201sin1.46473}t\\ -\text{0.0087362sin4.39419}t-\text{9.75925}\times {\text{10}}^{-\text{6}}\text{sin7.32364}t+\text{5.34464}\times {\text{10}}^{-\text{6}}\text{sin10.2531}t\\ +\text{1.38416}\times {\text{10}}^{-\text{9}}\text{sin13.1826}t-\text{7.53144}\times {\text{10}}^{-\text{9}}\text{sin16.112}t\end{array}$

(85b)

$\omega =\text{1.46473}$

(85c)

The phase portrait curves of the fifth-order HAM approximation are shown in Figure 8, compared with that of fourth-order Runge–Kutta method. It is noteworthy that traditional numerical methods such as fourth-order Runge–Kutta method cannot describe the phase portrait curve, but HAM can.

Figure 8.

Comparison of the phase portrait curves of the fifth-order HAM approximation with that of fourth-order Runge–Kutta method for ${\overline{\omega }}_0=1$ , $\beta =0$ , $\gamma =1$ , $\alpha =-0.8$ , $\tau =5.5$ .

The periodic solution by HAM with (α,τ) located in region IV

We use the Mathematica software to the study of the subject when $(\alpha ,\tau )$ is located in region IV, letting

${\overline{\omega }}_0=1,\beta =0,\gamma =1,\alpha =-0.8,\tau =6.2$

(86)

We can get two sets of the initial value from equation (60), as follows

$u_0(\theta )=\text{1.10751}\text{cos}\theta ,v_0(\theta )=\text{1.10751}\text{cos}\theta ,{\omega }_0=\text{0.6966154}$

(87a)

$u_0(\theta )=\text{1.38445}\text{cos}\theta ,v_0(\theta )=\text{1.38445}\text{cos}\theta ,{\omega }_0=\text{1.37035447}$

(87b)

It means that system (1) has two periodic solutions, as is in accord with the theoretical analysis in MSM and normal form equations section.

Choose equation (87a) as the initial guess, and let $\hslash =-0.85$ , then the fifth-order approximation can be obtained:

$\begin{array}{l}{x}_1(t)=\text{1.13937cos0.693119}t-\text{0.0130678cos2.07936}t-\text{0.00486125cos3.46559}t\\ +\text{6.676}\times {\text{10}}^{-\text{6}}\text{cos4.85183}t+\text{6.67815}\times {\text{10}}^{-\text{6}}\text{cos6.23807}t+\text{0.202144sin0.693119}t\\ -\text{0.0676095sin2.07936}t-\text{0.000174728sin3.46559}t+\text{0.000222864sin4.85183}t\\ -\text{1.13555}\times {\text{10}}^{-\text{7}}\text{sin6.23807}t-\text{8.20136}\times {\text{10}}^{-\text{8}}\text{sin7.62431}t\end{array}$

(88a)

$\begin{array}{l}{x}_2(t)=\text{1.13937cos0.693119}t-\text{0.0130678cos2.07936}t-\text{0.00486125cos3.46559}t\\ +\text{6.676}\times {\text{10}}^{-\text{6}}\text{cos4.85183}t+\text{6.67815}\times {\text{10}}^{-\text{6}}\text{cos6.23807}t+\text{0.202144sin0.693119}t\\ -\text{0.0676095sin2.07936}t-\text{0.000174728sin3.46559}t+\text{0.000222864sin4.85183}t\\ -\text{1.13555}\times {\text{10}}^{-\text{7}}\text{sin6.23807}t-\text{8.20136}\times {\text{10}}^{-\text{8}}\text{sin7.62431}t\end{array}$

(88b)

$\omega =\text{0.693119}$

(88c)

Choose equation (87b) as the initial guess, and let $\hslash =-0.5$ , then the fifth-order approximation can be obtained:

$\begin{array}{c}{x}_1(t)=\text{1.36281cos1.36661}t-\text{0.00489138cos4.09983}t-\text{0.00323062cos6.83305}t\\ -\text{0.0000222699cos9.56627}t+\text{0.0000250936cos12.2995}t+\text{0.157406sin1.36661}t\\ -\text{0.0531801sin4.09983}t+\text{0.0000693371sin6.83305}t+\text{0.000258585sin9.56627}t\\ +\text{3.57706}\times {10}^{-9}\text{sin12.2995}t\text{-1.98183}\times {10}^{-6}\text{sin15.0327}t\end{array}$

(89a)

$\begin{array}{l}{x}_2(t)=\text{1.36281cos1.36661}t-\text{0.00489138cos4.09983}t-\text{0.00323062cos6.83305}t\\ -\text{0.0000222699cos9.56627}t+\text{0.0000250936cos12.2995}t+\text{0.157406sin1.36661}t\\ -\text{0.0531801sin4.09983}t+\text{0.0000693371sin6.83305}t+\text{0.000258585sin9.56627}t\\ +\text{3.57706}\times {10}^{-9}\text{sin12.2995}t\text{-1.98183}\times {10}^{-6}\text{sin15.0327}t\end{array}$

(89b)

$\omega =\text{1.36661}$

(89c)

The phase portrait curves of the fifth-order HAM approximation are shown in Figure 9, compared with that of fourth-order Runge–Kutta method. We found that the approximate solution agrees with the numerical solution.

Figure 9.

Comparison of the phase portrait curves of the fifth-order HAM approximation with that of fourth-order Runge–Kutta method for ${\overline{\omega }}_0=1$ , $\beta =0$ , $\gamma =1$ , $\alpha =-0.8$ , $\tau =6.2$ .

We obtain the periodic solution by HAM with (α,τ) located in regions II, III, IV, V, and VI. Figures 5 to 9 show that the HAM can obtain stable periodic solutions and unstable periodic solutions, but the Runge–Kutta method does not get unstable periodic solutions, so the HAM is suitable for coupled nonlinear systems with delay coupling.

Conclusions

Recently, more and more researchers are devoted to the study of nonlinear dynamical systems, especially time-delay nonlinear systems ⁴⁰ and coupled nonlinear systems.^41,42 In this paper, we study the double Hopf bifurcation of coupled nonlinear systems with delayed coupling. Firstly, we obtain the critical value $\left({\alpha }_c,{\tau }_c\right)$ where double Hopf bifurcation occurs. Secondly, we get the normal form equations of 3:5 weakly resonant double Hopf bifurcation by the MSM. By analyzing the bifurcation equations, the dynamical behaviors of the original system are studied. It is shown that the system with ${\overline{\omega }}_0=1,$ $\beta =0$ has a unique zero solution when $\left(\alpha ,\tau \right)$ is located in region I; it has a stable periodic solution when $\left(\alpha ,\tau \right)$ is located in region II or VI; it has a stable periodic solution and an unstable periodic solution when $\left(\alpha ,\tau \right)$ is located in region III or V; it has two stable periodic solutions and an unstable quasiperiodic solution when $\left(\alpha ,\tau \right)$ is located in region IV. Finally, we apply the HAM to study the periodic solutions of region II to VI, and compare them with Runge–Kutta method. The results show that the HAM can obtain stable periodic solutions and unstable periodic solutions, but the Runge–Kutta method does not get unstable periodic solutions. Hence, the HAM is suitable for the system, and is a more powerful tool for the research coupled nonlinear systems with delay coupling. Due to the limit of length and energy, this article does not get the quasiperiodic analytical approximate solution in region IV, and it will be the follow-up research. The method used in this paper is only suitable for weakly resonant double Hopf bifurcation. How to obtain the strong resonant double Hopf bifurcation and use the Akbari–Ganji method ⁴³ is our next research task.