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.

Keywords

Double Hopf bifurcation nonlinear coupling system method of multiple scales homotopy analysis method periodic solution

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) + (γ x_{1}^{2} (t) - β) {\dot{x}}_{1} (t) + {\bar{ω}}_{0}^{2} x_{1} (t) = α {\dot{x}}_{2} (t - τ)

(1a)

{\ddot{x}}_{2} (t) + (γ x_{2}^{2} (t) - β) {\dot{x}}_{2} (t) + {\bar{ω}}_{0}^{2} x_{2} (t) = α {\dot{x}}_{1} (t - τ)

(1b)

where,

{\bar{ω}}_{0} > 0

is the natural frequency,

α < 0

is coupling intensities,

τ

is time delay,

β

is the damping coefficient, and

γ

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 $(α_{c}, τ_{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 $α^{2} > β^{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} = - {\bar{ω}}_{0}^{2} y_{1} + β y_{2} - γ y_{1}^{2} y_{2} + α y_{4 τ}

(2b)

{\dot{y}}_{3} = y_{4}

(2c)

{\dot{y}}_{4} = - {\bar{ω}}_{0}^{2} y_{3} + β y_{4} - γ y_{3}^{2} y_{4} + α y_{2 τ}

(2d)

Let

Y (t) = {(y_{1}, y_{2}, y_{3}, y_{4})}^{T}

(3a)

C = (\begin{matrix} 0 & 1 & 0 & 0 \\ - {\bar{ω}}_{0}^{2} & β & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & - {\bar{ω}}_{0}^{2} & β \end{matrix})

(3b)

D = (\begin{matrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & α \\ 0 & 0 & 0 & 0 \\ 0 & α & 0 & 0 \end{matrix})

(3c)

f (Y (t), Y (t - τ)) = (\begin{matrix} 0 \\ - γ y_{1}^{2} y_{2} \\ 0 \\ - γ y_{3}^{2} y_{4} \end{matrix})

(3d)

Then the system (1) can be further represented as

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

(4)

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

\det (λ \cdot E - C - D \cdot e^{- λ τ}) = 0

(5)

where

E

is a unit matrix.

That means

\begin{matrix} λ & - 1 & 0 & 0 \\ {\bar{ω}}_{0}^{2} & λ - β & 0 & - α \cdot e^{- λ τ} \\ 0 & 0 & λ & - 1 \\ 0 & - α \cdot e^{- λ τ} & {\bar{ω}}_{0}^{2} & λ - β \end{matrix} | = 0

(6)

Simplifying it, we can get

{(λ^{2} + {\bar{ω}}_{0}^{2} - β λ)}^{2} = {(α λ e^{- λ τ})}^{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 ${\bar{ω}}_{0}^{2} > 0$ , $λ = 0$ is not a solution of equation (7), that is, equation (7) has no zero root.

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

β^{2} {\bar{ω}}^{2} - {({\bar{ω}}^{2} - {\bar{ω}}_{0}^{2})}^{2} = α^{2} {\bar{ω}}^{2} cos 2 \bar{ω} τ

(8a)

2 β \bar{ω} ({\bar{ω}}^{2} - {\bar{ω}}_{0}^{2}) = α^{2} {\bar{ω}}^{2} sin 2 \bar{ω} τ

(8b)

According to ${sin}^{2} 2 \bar{ω} τ + {cos}^{2} 2 \bar{ω} τ = 1,$ we can get

{[β^{2} {\bar{ω}}^{2} - {({\bar{ω}}^{2} - {\bar{ω}}_{0}^{2})}^{2}]}^{2} + 4 β^{2} {\bar{ω}}^{2} {({\bar{ω}}^{2} - {\bar{ω}}_{0}^{2})}^{2} = {(α^{2} {\bar{ω}}^{2})}^{2}

(9)

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

{[β^{2} w - {(w - {\bar{ω}}_{0}^{2})}^{2}]}^{2} + 4 β^{2} w {(w - {\bar{ω}}_{0}^{2})}^{2} = {(α^{2} w)}^{2}

(10)

Considering ${\bar{ω}}_{0}^{2} > 0, α^{2} \geq 0, β^{2} \geq 0,$ we can come to a conclusion: When $α^{2} > β^{2}$ , equation (10) has two and only two nonnegative roots, as follows

w_{1} = \frac{2 {\bar{ω}}_{0}^{2} + α^{2} - β^{2} - \sqrt{α^{2} - β^{2}} \sqrt{α^{2} - β^{2} + 4 {\bar{ω}}_{0}^{2}}}{2}

(11)

w_{2} = \frac{2 {\bar{ω}}_{0}^{2} + α^{2} - β^{2} + \sqrt{α^{2} - β^{2}} \sqrt{α^{2} - β^{2} + 4 {\bar{ω}}_{0}^{2}}}{2}

(12)

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

λ_{1, 2} = \pm i ω_{-}

(13a)

λ_{3, 4} = \pm i ω_{+}

(13b)

where,

ω_{-} = \sqrt{w_{1}} > 0

ω_{+} = \sqrt{w_{2}} > 0

ω_{-} \neq ω_{+}

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

Further, let $ω_{-} : ω_{+} = 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{{(α^{2} - β^{2})}^{2} + 4 {\bar{ω}}_{0}^{2} (α^{2} - β^{2})}{{\bar{ω}}_{0}^{4}}

(14)

It is clear that when ${\bar{ω}}_{0}$ and $β$ are selected, based on the above formula, $α$ can be expressed as a function of $k_{1} : k_{2}$ , such as $α = α (k_{1} : k_{2})$ .

Substituting $α = α (k_{1} : k_{2})$ into equations (12) and (13), we can accordingly get $ω_{-} = ω_{-} (k_{1} : k_{2})$ and $ω_{+} = ω_{+} (k_{1} : k_{2})$ . Substituting them into equation (8), we can gain

τ_{-} = \frac{1}{2 ω_{-}} (2 j π - \arccos \frac{β^{2} ω_{-}^{2} - {(ω_{-}^{2} - {\bar{ω}}_{0}^{2})}^{2}}{α^{2} ω_{-}^{2}}), j = 1, 2, \dots

(15)

τ_{+} = \frac{1}{2 ω_{+}} (2 j π + \arccos \frac{β^{2} ω_{+}^{2} - {(ω_{+}^{2} - {\bar{ω}}_{0}^{2})}^{2}}{α^{2} ω_{+}^{2}}), j = 1, 2, \dots

(16)

Based on the Hopf bifurcation theorem, a sufficient condition of double Hopf bifurcation of system (equation (1)) is $τ_{c} = τ_{-} = τ_{+}$ , 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 ${\bar{ω}}_{0} = 1$ , $β = 0$ , there exists a double Hopf bifurcation with $k_{1} : k_{2}$ resonance of system (1), and the bifurcation point is

(α_{c}, τ_{c}) = (\frac{k_{1} - k_{2}}{k_{2}} \sqrt{\frac{k_{2}}{k_{1}}}, \frac{k_{1} π}{k_{1} - k_{2}} \sqrt{\frac{k_{2}}{k_{1}}})

(17)

where

k_{1} = 2 j - 1

k_{2} = 2 j + 1

j = 1, 2, \dots

Proof: Simplicity, Let

{\bar{ω}}_{0} = 1, β = 0

(18)

then

α = \frac{k_{1} - k_{2}}{k_{2}} \sqrt{\frac{k_{2}}{k_{1}}}

(19a)

ω_{-} = \sqrt{\frac{k_{1}}{k_{2}}}

(19b)

ω_{+} = \sqrt{\frac{k_{2}}{k_{1}}}

(19c)

And further, we can obtain

τ_{-} = \frac{2 j π - π}{2} \sqrt{\frac{k_{2}}{k_{1}}}, j = 1, 2, \dots

(20)

τ_{+} = \frac{2 j π + π}{2} \sqrt{\frac{k_{1}}{k_{2}}}, j = 1, 2, \dots

(21)

Differentiating $λ$ with respect to $τ$ in equation (7), these can be obtained

λ^{'} (τ) [1 - \frac{1}{λ {(τ)}^{2}} - \frac{τ (λ {(τ)}^{2} + 1)}{λ (τ)}] = - λ {(τ)}^{2} - 1

(22)

and

{Re (\frac{1}{λ^{'} (τ)}) |}_{τ = τ_{-}} = \frac{k_{2} (k_{1} + k_{2})}{k_{1} (k_{1} - k_{2})} < 0

(23a)

{Re (\frac{1}{λ^{'} (τ)}) |}_{τ = τ_{+}} = \frac{k_{1} (k_{1} + k_{2})}{k_{2} (k_{2} - k_{1})} > 0

(23b)

Then

{s i g n (Re (\frac{d λ (τ)}{d τ})) |}_{τ = τ_{\pm}} = {s i g n (Re (\frac{1}{λ^{'} (τ)})) |}_{τ = τ_{\pm}} = {\begin{matrix} < 0, τ = τ_{-} \\ > 0, τ = τ_{+} \end{matrix}

(24)

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

k_{1} = 2 j - 1, k_{2} = 2 j + 1 j = 1, 2, \dots

(25)

accordingly

τ_{c} = \frac{k_{1} π}{k_{1} - k_{2}} \sqrt{\frac{k_{2}}{k_{1}}}

(26)

Therefore, when ${\bar{ω}}_{0} = 1$ , $β = 0$ , there exists a double Hopf bifurcation with $k_{1} : k_{2}$ resonance of system (1), and the bifurcation point is

(α_{c}, τ_{c}) = (\frac{k_{1} - k_{2}}{k_{2}} \sqrt{\frac{k_{2}}{k_{1}}}, \frac{k_{1} π}{k_{1} - k_{2}} \sqrt{\frac{k_{2}}{k_{1}}})

(27)

where

k_{1} = 2 j - 1

k_{2} = 2 j + 1

j = 1, 2, \dots

. □

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

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 ${\bar{ω}}_{0} = 1$ , $β = 0$ , there exists a double Hopf bifurcation with 3:5 weak resonance of system (1), and the bifurcation point is

(α_{c}, τ_{c}) = (- \frac{2}{\sqrt{15}}, \frac{π \sqrt{15}}{2})

(28)

with

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

(29)

For the convenience of calculation, set $γ = 1$ .

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

y_{1} (t) = ε y_{11} + ε^{2} y_{12} + ε^{3} y_{13} + ο (ε^{4})

(30a)

y_{2} (t) = ε y_{21} + ε^{2} y_{22} + ε^{3} y_{23} + ο (ε^{4})

(30b)

y_{3} (t) = ε y_{31} + ε^{2} y_{32} + ε^{3} y_{33} + ο (ε^{4})

(30c)

y_{4} (t) = ε y_{41} + ε^{2} y_{42} + ε^{3} y_{43} + ο (ε^{4})

(30d)

where

y_{i j} = y_{i j} (T_{0}, T_{2}),

T_{0} = t,

T_{2} = ε^{2} t,

0 < ε ≪ 1.

While the bifurcation parameters are ordered as

α = α_{c} + ε^{2} α_{ε}, τ = τ_{c} + ε^{2} τ_{ε}

(31)

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

\begin{matrix} y_{2} (t - τ) = ε y_{21} (t - τ, ε^{2} (t - τ)) + ε^{2} y_{22} (t - τ, ε^{2} (t - τ)) + ε^{3} y_{23} (t - τ, ε^{2} (t - τ)) + ο (ε^{4}) \\ = ε y_{21} (T_{0} - τ_{c}, T_{2}) + ε^{2} y_{22} (T_{0} - τ_{c}, T_{2}) + ε^{3} [- τ_{ε} D_{0} y_{21 τ} - τ_{c} D_{2} y_{21 τ} + y_{23 τ}] + ο (ε^{4}) \end{matrix}

(32a)

\begin{matrix} y_{4} (t - τ) = ε y_{41} (t - τ, ε^{2} (t - τ)) + ε^{2} y_{42} (t - τ, ε^{2} (t - τ)) + ε^{3} y_{43} (t - τ, ε^{2} (t - τ)) + ο (ε^{4}) \\ = ε y_{41} (T_{0} - τ_{c}, T_{2}) + ε^{2} y_{42} (T_{0} - τ_{c}, T_{2}) + ε^{3} [- τ_{ε} D_{0} y_{41 τ} - τ_{c} D_{2} y_{41 τ} + y_{43 τ}] + ο (ε^{4}) \end{matrix}

(32b)

where

y_{i j τ} = y_{i j} (T_{0} - τ_{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 $ε$ , the following perturbative equations are obtained:

ε \to

D_{0} y_{11} = y_{21}

(33a)

D_{0} y_{21} = - y_{11} - \frac{2 \sqrt{15}}{15} y_{41 τ}

(33b)

D_{0} y_{31} = y_{41}

(33c)

D_{0} y_{41} = - y_{31} - \frac{2 \sqrt{15}}{15} y_{21 τ}

(33b)

ε^{2} \to

D_{0} y_{12} = y_{22}

(34a)

D_{0} y_{22} = - y_{12} - \frac{2 \sqrt{15}}{15} y_{42 τ}

(34b)

D_{0} y_{32} = y_{42}

(34c)

D_{0} y_{42} = - y_{32} - \frac{2 \sqrt{15}}{15} y_{22 τ}

(34d)

ε^{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} [- τ_{ε} D_{0} (y_{41 τ}) - \frac{\sqrt{15}}{2} π D_{2} (y_{41 τ}) + y_{43 τ}] - {y_{11}}^{2} y_{21} + α_{ε} y_{41 τ}

(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} [- τ_{ε} D_{0} (y_{21 τ}) - \frac{\sqrt{15}}{2} π D_{2} (y_{21 τ}) + y_{23 τ}] - {y_{31}}^{2} y_{41} + α_{ε} y_{21 τ}

(35d)

Equation (33) has the following general solution:

y_{11} = Α_{1} e^{i ω_{1} T_{0}} +_{Α2} e^{i ω_{2} T_{0}} + {\bar{A}}_{1} e^{- i ω_{1} T_{0}} + {\bar{A}}_{2} e^{- i ω_{2} T_{0}}

(36a)

y_{21} = A_{1} i \frac{\sqrt{15}}{5} e^{i ω_{1} T_{0}} + A_{2} i \frac{\sqrt{15}}{3} e^{i ω_{2} T_{0}} + {\bar{A}}_{1} i \frac{\sqrt{15}}{5} e^{- i ω_{1} T_{0}} + {\bar{A}}_{2} i \frac{\sqrt{15}}{3} e^{- i ω_{2} T_{0}}

(36b)

y_{31} = - A_{1} e^{i ω_{1} T_{0}} - A_{2} e^{i ω_{2} T_{0}} - {\bar{A}}_{1} e^{- i ω_{1} T_{0}} - {\bar{A}}_{2} e^{- i ω_{2} T_{0}}

(36c)

y_{41} = - A_{1} i \frac{\sqrt{15}}{5} \cdot e^{i ω_{1} T_{0}} - A_{2} i \frac{\sqrt{15}}{3} \cdot e^{i ω_{2} T_{0}} - {\bar{A}}_{1} i \frac{\sqrt{15}}{5} e^{- i ω_{1} T_{0}} - {\bar{A}}_{2} i \frac{\sqrt{15}}{3} e^{- i ω_{2} T_{0}}

(36c)

where

A_{j} = 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 $ε$ , we can obtain the equations including $A_{1} = D_{2} A_{1}$ and $A_{2} = D_{2} A_{2}$ :

{\dot{A}}_{1} = C_{1} A_{1} + C_{11 \bar{1}} A_{1}^{2} {\bar{A}}_{1} + C_{12 \bar{2}} A_{1} A_{2} {\bar{A}}_{2}

(37a)

{\dot{A}}_{2} = C_{2} A_{2} + C_{1 \bar{1} 2} A_{1} {\bar{A}}_{1} A_{2} + C_{22 \bar{2}} A_{2}^{2} {\bar{A}}_{2}

(37b)

where the complex coefficients

C_{l}, C_{i j k}

are expressed as follow:

C_{1} = - 0. 477692 + 0. 369119 i - (0.185009 + 0. 157041 i) α + (0.0628164 - 0.0 740038 i) τ

(38a)

C_{2} = 0.39156 + 1.09157 i - (0.252751 - 0. 1128725 i) α - (0.0858166 + 0. 168501 i) τ

(38b)

C_{11 \bar{1}} = - 0.157041 + 0.185009 i

(38c)

C_{12 \bar{2}} = - 0.314082 + 0.370019 i

(38d)

C_{1 \bar{1} 2} = - 0.25745 - 0.505502 i

(38e)

C_{22 \bar{2}} = - 0.128725 - 0.252751 i

(38f)

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

A_{1} = \frac{1}{2} r_{1} e^{i θ_{1}}, A_{2} = \frac{1}{2} r_{2} e^{i θ_{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. 477692 - 0. 185009 α + 0.0 628164 τ) r_{1} - 0.0392602 r_{1}^{3} - 0.0785204 r_{1} r_{2}^{2}

(40a)

{\dot{r}}_{2} = (0.39156 - 0. 252751 α - 0.0 858166 τ) r_{2} - 0.0643625 r_{1}^{2} r_{2} - 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_{11}, 0), E_{2} = (0, r_{21}), E_{3} = (r_{12}, r_{22})

(41)

where

r_{11} = \sqrt{- 12.1673 - 4.71239 α +1.6000 τ}

(42a)

r_{21} = \sqrt{12.1673 - 7.85398 α - 2.66667 τ}

(42b)

r_{12} = \sqrt{12.1673 - 3.66519 α - 2.31111 τ}

(42c)

r_{22} = \sqrt{- 12.1673 - 0. 523599 α + 1.95556 τ}

(42d)

The existence of the above possible equilibrium points is dependent on $α$ and $τ$ . And according to the stability analysis of existing equilibrium points, the parameters $(α, τ)$ plane around the 3:5 weakly resonant double Hopf bifurcation point $(α_{c}, τ_{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.

I = {(α, τ) | r_{11} < 0, r_{21} < 0, r_{12} r_{22} < 0}

(43a)

II = {(α, τ) | r_{11} > 0, r_{21} < 0, r_{12} > 0, r_{22} < 0}

(43b)

III = {(α, τ) | r_{11} > 0, r_{21} > 0, r_{12} > 0, r_{22} < 0}

(43c)

IV = {(α, τ) | r_{11} > 0, r_{21} > 0, r_{12} > 0, r_{22} > 0}

(43d)

V = {(α, τ) | r_{11} > 0, r_{21} > 0, r_{12} < 0, r_{22} > 0}

(43e)

VI = {(α, τ) | r_{11} < 0, r_{21} > 0, r_{12} < 0, r_{22} > 0}

(43f)

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 $α$ and $τ$ , so region I is usually called as an amplitude death region.

With $(α, τ)$ changing into region II, the trivial equilibrium point $E_{0} = (0, 0)$ loses its stability and a stable equilibrium point $E_{1} = (r_{11}, 0)$ appears. Correspondingly, there is a stable periodic solution with a frequency close to $ω_{+} \approx 1.29099$ .

When $(α, τ)$ enters into region III, the unstable equilibrium point $E_{0} = (0, 0)$ and the stable equilibrium point $E_{1} = (r_{11}, 0)$ still exist, and an unstable equilibrium point $E_{2} = (0, r_{21})$ appears. Correspondingly, there exists a stable periodic solution with a frequency close to $ω_{+}$ and an unstable periodic solution.

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

With $(α, τ)$ changing into region V, the equilibrium point $E_{3} = (r_{12}, r_{22})$ disappears and there exist three equilibrium points, $E_{0} = (0, 0)$ , $E_{1} = (r_{11}, 0)$ and $E_{2} = (0, r_{21})$ . In contrast with the situation in region III, the equilibrium point $E_{1} = (r_{11}, 0)$ is unstable and $E_{2} = (0, r_{21})$ is stable. Correspondingly, there exists a stable periodic solution with a frequency close to $ω_{-}$ and an unstable periodic solution.

When $(α, τ)$ enters into region VI, the equilibrium point $E_{1} = (r_{11}, 0)$ disappears and there exists an unstable equilibrium point $E_{0} = (0, 0)$ and a stable equilibrium point $E_{2} = (0, r_{21})$ . Correspondingly, there is a stable periodic solution with a frequency close to $ω_{-} \approx 0.774597$ .

To demonstrate the above theoretical results, we choose six groups of parameter values $(α, τ)$ for system (1) with ${\bar{ω}}_{0} = 1$ , $β = 0$ , and the results are given in Figure 4. It is obvious that numerical simulations are well consistent with the theoretical results.

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, - τ \leq t \leq 0

(44)

and the frequency of initial value solution is

ω

Let $θ = ω t, x_{1} (t) = u (θ), x_{2} (t) = v (θ)$ , then system (1) can be expressed as

ω^{2} u^{″} (θ) + ω (γ u^{2} (θ) - β) u^{'} (θ) + {\bar{ω}}_{0}^{2} u (θ) = α ω v^{'} (θ - ω τ)

(45a)

ω^{2} v^{″} (θ) + ω (γ v^{2} (θ) - β) v^{'} (θ) + {\bar{ω}}_{0}^{2} v (θ) = α ω u^{'} (θ - ω τ)

(45b)

with the initial condition

u (θ) = A (θ), u^{'} (θ) = B (θ), v (θ) = C, v^{'} (θ) = 0, - ω τ \leq θ \leq 0

(46)

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

{cos (2 k + 1) θ, sin (2 k + 1) θ | k = 0, 1, 2, \dots}

(47)

The initial guess solution is chosen as

u_{0} (θ) = a_{0} cos θ + b_{0} sin θ, v_{0} (θ) = c_{0} cos θ

(48)

and the linear operator is expressed as

L (\begin{matrix} U (θ, q) \\ V (θ, q) \end{matrix}) = (\begin{matrix} \frac{\partial^{2} U (θ, q)}{\partial θ^{2}} + U (θ, q) \\ \frac{\partial^{2} V (θ, q)}{\partial θ^{2}} + V (θ, q) \end{matrix})

(49)

with the property

L (\begin{matrix} C_{1} cos θ + C_{2} sin θ \\ C_{3} cos θ + C_{4} sin θ \end{matrix}) = (\begin{matrix} 0 \\ 0 \end{matrix})

(50)

where

C_{1}

C_{2}

C_{3}

_, and

C_{4}

are constants.

From equation (45), we define a nonlinear operator

\begin{matrix} N (U (θ, q)) = W^{2} (q) \frac{\partial^{2} U (θ, q)}{\partial θ^{2}} + W (q) (γ U {(θ, q)}^{2} - β) \frac{\partial U (θ, q)}{\partial θ} \\ + {\bar{ω}}_{0}^{2} U (θ, q) - α W (q) \frac{\partial V (θ - W (q) τ, q)}{\partial θ} \end{matrix}

(51a)

\begin{matrix} N (V (θ, q)) = W^{2} (q) \frac{\partial^{2} V (θ, q)}{\partial θ^{2}} + W (q) (γ V {(θ, q)}^{2} - β) \frac{\partial V (θ, q)}{\partial θ} \\ + {\bar{ω}}_{0}^{2} V (θ, q) - α W (q) \frac{\partial U (θ - W (q) τ, q)}{\partial θ} \end{matrix}

(51b)

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

(1 - q) L (\begin{matrix} U (θ, q) - u_{0} (θ) \\ V (θ, q) - v_{0} (θ) \end{matrix}) = q ℏ N (\begin{matrix} U (θ, q) \\ V (θ, q) \end{matrix})

(52)

with the initial conditions

U (θ, q) = a (q), \frac{\partial U (θ, q)}{\partial θ} = b (q), V (θ, q) = c (q), \frac{\partial V (θ, q)}{\partial θ} = 0, - ω τ \leq θ \leq 0

(53)

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

U (θ, 0) = u_{0} (θ), V (θ, 0) = v_{0} (θ), W (0) = ω_{0}

(54)

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

U (θ, 1) = u (θ), V (θ, 1) = v (θ), W (1) = ω

(55)

Hence, as $q$ changes from 0 to 1, $U (θ, q)$ and $V (θ, q)$ , respectively vary from the initial guess $u_{0} (θ)$ and $v_{0} (θ)$ to the unknown $u (θ)$ and $v (θ)$ . Likewise, $W (q)$ varies from the initial guess frequency $ω_{0}$ to the physical frequency $ω$ .

According to Taylor series expansion, we obtain

U (θ, q) = u_{0} (θ) + \sum_{m = 1}^{+ \infty} u_{m} (θ) q^{m}

(56a)

V (θ, q) = v_{0} (θ) + \sum_{m = 1}^{+ \infty} v_{m} (θ) q^{m}

(56b)

W (q) = ω_{0} + \sum_{m = 1}^{+ \infty} ω_{m} q^{m}

(56c)

where

{u_{m} (θ) = \frac{1}{m!} \frac{\partial^{m} U (θ, q)}{\partial q^{m}} |}_{q = 0}

(57a)

{v_{m} (θ) = \frac{1}{m!} \frac{\partial^{m} V (θ, q)}{\partial q^{m}} |}_{q = 0}

(57b)

{ω_{m} (θ) = \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 $ℏ$ is so properly chosen that power series equation (56) converge at $q = 1$ , the series solutions of equation (56) can be expressed as

u (θ) = u_{0} (θ) + \sum_{m = 1}^{+ \infty} u_{m} (θ)

(58a)

v (θ) = v_{0} (θ) + \sum_{m = 1}^{+ \infty} v_{m} (θ)

(58b)

ω = ω_{0} + \sum_{m = 1}^{+ \infty} ω_{m}

(58c)

For the sake of simplicity, we define the vectors

{\vec{u}}_{m} = {u_{0} (θ), u_{1} (θ), \dots, u_{m} (θ)}

(59a)

{\vec{v}}_{m} = {v_{0} (θ), v_{1} (θ), \dots, v_{m} (θ)}

(59b)

{\vec{ω}}_{m} = {ω_{0}, ω_{1}, \dots, ω_{m}}

(59c)

By differentiating the zeroth-order deformation equation (52) $m$ $(m \geq 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 (\begin{matrix} u_{m} (θ) - χ_{m} u_{m - 1} (θ) \\ v_{m} (θ) - χ_{m} v_{m - 1} (θ) \end{matrix}) = ℏ (\begin{matrix} R_{1, m} ({\vec{u}}_{m - 1}, {\vec{ω}}_{m - 1}) \\ R_{2, m} ({\vec{v}}_{m - 1}, {\vec{ω}}_{m - 1}) \end{matrix})

(60)

where

χ_{m} = {\begin{matrix} 0, & m = 1 \\ 1, & m \geq 2 \end{matrix}

(61)

\begin{array}{l} R_{1, m} ({\vec{u}}_{m - 1}, {\vec{ω}}_{m - 1}) = {\frac{1}{(m - 1)!} \frac{\partial^{m - 1} N (U (θ, q))}{\partial q^{m - 1}} |}_{q = 0} \\ = \sum_{i = 0}^{m - 1} {u^{″}}_{m - 1 - i} (θ) (\sum_{j = 0}^{i} ω_{j} ω_{i - j}) + γ \sum_{i = 0}^{m - 1} {u^{'}}_{m - 1 - i} (θ) (\sum_{j = 0}^{i} ω_{i - j} (\sum_{k = 0}^{j} u_{k} u_{j - k})) \\ - β \sum_{i = 0}^{m - 1} ω_{m - 1 - i} u^{'} (θ) + {\bar{ω}}_{0}^{2} u_{m - 1} (θ) - α P_{m - 1} \end{array}

(62)

\begin{array}{l} R_{2, m} ({\vec{v}}_{m - 1}, {\vec{ω}}_{m - 1}) = {\frac{1}{(m - 1)!} \frac{\partial^{m - 1} N (V (θ, q))}{\partial q^{m - 1}} |}_{q = 0} \\ = \sum_{i = 0}^{m - 1} {v^{″}}_{m - 1 - i} (θ) (\sum_{j = 0}^{i} ω_{j} ω_{i - j}) + γ \sum_{i = 0}^{m - 1} {v^{'}}_{m - 1 - i} (θ) (\sum_{j = 0}^{i} ω_{i - j} (\sum_{k = 0}^{j} v_{k} v_{j - k})) \\ - β \sum_{i = 0}^{m - 1} ω_{m - 1 - i} {v^{'}}_{i} (θ) + {\bar{ω}}_{0}^{2} v_{m - 1} (θ) - α Q_{m - 1} \end{array}

(63)

\begin{array}{l} P_{m - 1} = {\frac{1}{(m - 1)!} \frac{\partial^{m - 1} (W (q) \frac{\partial V (θ - W (q) τ, q)}{\partial θ})}{\partial q^{m - 1}} |}_{q = 0} \\ = \sum_{i = 0}^{m - 1} ω_{i} [{v^{'}}_{m - 1 - i} (θ - ω_{0} τ) + \sum_{j = 1}^{m - 1 - i} \frac{{(- τ)}^{j}}{j!} (\sum_{k = 0}^{m - 1 - i - j} D_{m - 1 - i - k, j} \frac{\partial^{j} {v^{'}}_{k} (θ - ω_{0} τ)}{\partial θ^{j}})] \end{array}

(64)

\begin{array}{l} Q_{m - 1} = {\frac{1}{(m - 1)!} \frac{\partial^{m - 1} (W (q) \frac{\partial U (θ - W (q) τ, q)}{\partial θ})}{\partial q^{m - 1}} |}_{q = 0} \\ = \sum_{i = 0}^{m - 1} ω_{i} [{u^{'}}_{m - 1 - i} (θ - ω_{0} τ) + \sum_{j = 1}^{m - 1 - i} \frac{{(- τ)}^{j}}{j!} (\sum_{k = 0}^{m - 1 - i - j} D_{m - 1 - i - k, j} \frac{\partial^{j} {u^{'}}_{k} (θ - ω_{0} τ)}{\partial θ^{j}})] \end{array}

(65)

D_{n, 1} = ω_{n}, D_{n, l} = \sum_{j = 1}^{n - 1} ω_{j} B_{n - j, l - 1}, l \geq 2

(66)

and the initial conditions

u_{m} (θ) = a_{m}, {u^{'}}_{m} (θ) = b_{m}, v_{m} (θ) = c_{m}, {v^{'}}_{m} (θ) = 0, - ω τ \leq θ \leq 0

(67)

According to the rule of solution expression, $u_{m}$ and $v_{m}$ should not contain the so-called secular terms of $θ sin θ$ and $θ cos θ$ . Correspondingly, $ℏ R_{1, m} ({\vec{u}}_{m - 1}, {\vec{ω}}_{m - 1})$ and $ℏ R_{2, m} ({\vec{u}}_{m - 1}, {\vec{ω}}_{m - 1})$ should not contain the terms of $sin θ$ and $cos θ$ . Setting their coefficients to zeroes, we obtains

\frac{1}{π} \int_{0}^{2 π} ℏ R_{1, m} ({\vec{u}}_{m - 1}, {\vec{ω}}_{m - 1}) cos θ d θ = 0

(68a)

\frac{1}{π} \int_{0}^{2 π} ℏ R_{1, m} ({\vec{u}}_{m - 1}, {\vec{ω}}_{m - 1}) sin θ d θ = 0

(68b)

\frac{1}{π} \int_{0}^{2 π} ℏ R_{2, m} ({\vec{u}}_{m - 1}, {\vec{ω}}_{m - 1}) cos θ d θ = 0

(68c)

\frac{1}{π} \int_{0}^{2 π} ℏ R_{2, m} ({\vec{u}}_{m - 1}, {\vec{ω}}_{m - 1}) sin θ d θ = 0

(68d)

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

Hence, $ℏ R_{1, m} ({\vec{u}}_{m - 1}, {\vec{ω}}_{m - 1})$ and $ℏ R_{2, m} ({\vec{u}}_{m - 1}, {\vec{ω}}_{m - 1})$ can be expressed as

ℏ R_{1, m} ({\vec{u}}_{m - 1}, {\vec{ω}}_{m - 1}) = \sum_{n = 1} a_{1, n} cos (2 n + 1) θ + \sum_{n = 1} a_{2, n} sin (2 n + 1) θ

(69a)

ℏ R_{2, m} ({\vec{u}}_{m - 1}, {\vec{ω}}_{m - 1}) = \sum_{n = 1} b_{1, n} cos (2 n + 1) θ + \sum_{n = 1} b_{2, n} sin (2 n + 1) θ

(69b)

Further, by solving equation (60), it yields

\begin{array}{l} u_{m} (θ) = χ_{m} u_{m - 1} (θ) - \sum_{n = 1} \frac{a_{1, n}}{4 n (n + 1)} cos (2 n + 1) θ \\ + \sum_{n = 1} \frac{a_{2, n}}{4 n (n + 1)} sin (2 n + 1) θ + C_{1, m} sin θ + C_{2, m} cos θ \end{array}

(70a)

\begin{array}{l} v_{m} (θ) = χ_{m} v_{m - 1} (θ) - \sum_{n = 1} \frac{b_{1, n}}{4 n (n + 1)} cos (2 n + 1) θ \\ + \sum_{n = 1} \frac{b_{2, n}}{4 n (n + 1)} sin (2 n + 1) θ + C_{3, m} sin θ + C_{4, m} cos θ \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} (θ)$ and $v_{m} (θ)$ 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]} (θ) = u_{0} (θ) + \sum_{m = 1}^{n} u_{m} (θ)

(71a)

v_{[n]} (θ) = v_{0} (θ) + \sum_{m = 1}^{n} v_{m} (θ)

(71b)

ω_{[n]} = ω_{0} + \sum_{m = 1}^{n} ω_{m}

(71c)

Numerical simulation

According to the MSM and normal form equations section, parameters $(α, τ)$ plane around the 3:5 weakly resonant double Hopf bifurcation point $(α_{c}, τ_{c})$ can be divided into six distinct dynamical regions, and when $(α, τ)$ 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 $(α, τ)$ 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 $(α, τ)$ 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 $(α, τ)$ is located in regions IIand VI. For the region II, letting

{\bar{ω}}_{0} = 1, β = 0, γ = 1, α = - 0.5, τ = 7

(72)

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

u_{0} (θ) = 1.108515 cos θ, v_{0} (θ) = 1.108515 cos θ, ω_{0} = 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 $ℏ = - 0.5$ , the fifth-order approximation can be obtained:

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

(74a)

\begin{array}{l} x_{2} (t) = 1.09691 cos 1.214241 t - 0.00179299 cos 3.642723 t - 0.00153991 cos 6.07120 t \\ - 0.00000390 cos 8.49969 t + 0.00000394 cos 10.92817 t + 0.0982843 sin 1.214241 t \\ - 0.0329664 sin 3.64272 t + 0.000001792 sin 6.07120 t + 0.0000866711 sin 8.49969 t \\ + 0.0000000218 sin 10.9282 t - 0.0000000948 sin 13.3567 t \end{array}

(74b)

ω = 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.

When $(α, τ)$ is located in region VI, let

{\bar{ω}}_{0} = 1, β = 0, γ = 1, α = - 0.5, τ = 5

(75)

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

u_{0} (θ) = 1.0814565 cos θ, v_{0} (θ) = 1.0814565 cos θ, ω_{0} = 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 $ℏ = - 1.5$ , the fifth-order approximation can be obtained:

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

(77a)

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

(77b)

ω = 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)$ .

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

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

{\bar{ω}}_{0} = 1, β = 0, γ = 1, α = - 0.9, τ = 7

(78)

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

u_{0} (θ) = 1.74718 cos θ, v_{0} (θ) = 1.74718 cos θ, ω_{0} = 1.26658535

(79a)

u_{0} (θ) = 0.761786 cos θ, v_{0} (θ) = 0.761786 cos θ, ω_{0} = 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 $ℏ = - 0.2$ , then the fifth-order approximation can be obtained:

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

(80a)

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

(80b)

ω = 1.25801

(80c)

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

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

(81a)

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

(81b)

ω = 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.

Similarly, let’s consider the condition that $(α, τ)$ is located in region V. Set

{\bar{ω}}_{0} = 1, β = 0, γ = 1, α = - 0.8, τ = 5.5

(82)

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

u_{0} (θ) = 1.41094 cos θ, v_{0} (θ) = 1.41094 cos θ, ω_{0} = 0.7347193

(83a)

u_{0} (θ) = 0.809123 cos θ, v_{0} (θ) = 0.809123 cos θ, ω_{0} = 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 $ℏ = - 0.6$ , then the fifth-order approximation can be obtained:

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

(84a)

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

(84b)

ω = 0.729193

(84c)

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

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

(85a)

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

(85b)

ω = 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.

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

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

{\bar{ω}}_{0} = 1, β = 0, γ = 1, α = - 0.8, τ = 6.2

(86)

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

u_{0} (θ) = 1.10751 cos θ, v_{0} (θ) = 1.10751 cos θ, ω_{0} = 0.6966154

(87a)

u_{0} (θ) = 1.38445 cos θ, v_{0} (θ) = 1.38445 cos θ, ω_{0} = 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 $ℏ = - 0.85$ , then the fifth-order approximation can be obtained:

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

(88a)

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

(88b)

ω = 0.693119

(88c)

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

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

(89a)

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

(89b)

ω = 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.

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 $(α_{c}, τ_{c})$ 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 ${\bar{ω}}_{0} = 1,$ $β = 0$ has a unique zero solution when $(α, τ)$ is located in region I; it has a stable periodic solution when $(α, τ)$ is located in region II or VI; it has a stable periodic solution and an unstable periodic solution when $(α, τ)$ is located in region III or V; it has two stable periodic solutions and an unstable quasiperiodic solution when $(α, τ)$ 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.

Footnotes

Acknowledgments

We are grateful to the anonymous reviewers for their constructive comments and suggestions.

Authors’ contributions

All the authors contribute equally and significantly in writing this paper. All the authors read and approve the final manuscript.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors gratefully acknowledge the support of the National Natural Science Foundation of China (NNSFC) through grant No.11572288.

References

Kevin

Richard

Howard

Dynamics of three coupled van der Pol oscillators with application to circadian rhythms. Commun Nonlinear Sci Numeric Simulat 2007; 12: 794–803.

Barron

Sen

Synchronization of four coupled van der Pol oscillators. Nonlinear Dyn 2009; 56: 357–367.

Wei

Hopf bifurcation analysis in a delayed Nicholson blowflies equation. Nonlinear Analysis Theory Methods Appl 2005; 60: 1351–1367.

Jiang

Wei

Stability and bifurcation analysis in a delayed SIR model. Chaos Solitons Fractals 2008; 35: 609–619.

Bélair

Campbell

SA.

Stability and bifurcations of equilibria in a multiple-delayed differential equation.

SIAM J Appl Math 1994; 54: 1402–1424.

Chung

KW.

Effects of time delayed position feedback on a van Pol-Duffing oscillator. Physica D 2003; 180: 17–39.

Buono

Bélair

Restrications and unfolding of double Hopf bifurcation in functional differential equations. J Different Equat 2003; 189: 234–266.

Szalai

Stépán

Period doubling bifurcation and center manifold reduction in a time-periodic and time-delayed model of machining. J Vibrat Control 2010; 16: 1169–1187.

, et al. Impulsive control and Hopf bifurcation of a three-dimensional chaotic system. J Vibrat Contr 2013; 20: 1361–1368.

10.

Feng

ZS.

Double Hopf bifurcation for van der Pol Duffing oscillator with parametric delay feedback control. J Mathematica Anal Appl 2008; 338: 993–1007.

11.

Yan

Chu

YD.

Stability and bifurcation analysis for a delayed Lotka-Volterra predator-prey system. J Comput Appl Math 2006; 196: 198–210.

12.

Song

Zhang

TH.

Bifurcation, amplitude death and oscillation patterns in a system of three coupled van der Pol oscillators with diffusively delayed velocity coupling. Chaos 2011; 21: 023111.

13.

Zhang

Zheng

Wang

LC.

Multiple Hopf bifurcation of three coupled van der Pol oscillators with delay.

Appl Math Comput 2011; 217: 7155–7166.

14.

Wang

Multiple scales analysis for double Hopf bifurcation with 1: 3 resonance. Nonlinear Dyn 2011; 66: 39–51.

15.

Jiang

WangH

Double Hopf bifurcation and quasi periodic attractors in delay coupled limit cycle oscillators. J Math Anal Appl 2012; 387: 1114–1126.

16.

Weak resonant double Hopf bifurcations in an inertial four neuron model with time delay. Int J Neural Syst 2012; 22: 63–75.

17.

Wang

Chen

LJ.

Weak and non-resonant double Hopf bifurcations in m coupled van der Pol oscillators with delay coupling. Appl Math Model 2015; 39: 3094–3102.

18.

Amer

El-Sayed

Kotb

AA.

Nonlinear vibration and of the Duffing oscillator to parametric excitation with time delay feedback. Nonlinear Dyn 2016; 85: 2497–2505.

19.

Zhang

LJ.

The principal parametric resonance of coupled van der pol oscillators under feedback control. Shock Vibrat 2012; 19: 365–377.

20.

Liao

Chwang

AT.

Application of homotopy analysis method in nonlinear oscillation. J Appl Mech 1998; 65: 914–922.

21.

Chung

Chan

CL.

An efficient method for studying weak resonant double Hopf bifurcation in nonlinear systems with delayed feedbacks. SIAM J Appl Dyn Syst 2007; 6: 29–60.

22.

Rajasekar

Murali

Resonance behaviour and jump phenomenon in a two coupled Duffing–van der Pol oscillators. Chaos Solitons Fract 2004; 19: 925–934.

23.

Liao

SJ.

The proposed homotopy analysis techniques for the solution of nonlinear problems. PhD dissertation, Shanghai Jiao Tong University, Shanghai, 1992.

24.

You

Analytical approximations for the periodic motion of the Duffing system with delayed feedback. Numer Algor 2011; 56: 561–576.

25.

The application of homotopy analysis method to solve the coupled KdV equations. Int J Nonlinear Sci 2011; 11: 80–85.

26.

Kheiri

Jabbari

Homotopy analysis and Homotopy Padé methods for two-dimensional coupled Burgers’ equations. Iran J Mathemat Sci Informat 2011; 6: 23–31.

27.

Tasbozan

Esen

Yagmurlu

NM.

Approximate analytical solutions of fractional coupled mKdV equation by homotopy analysis method. OJAppS 2012; 2: 193–197.

28.

Jafarian

Ghaderi

Golmankhaneh

, et al. Homotopy analysis method for solving coupled Ramani equations. Roman J Phys 2014; 59: 112–136.

29.

Eigoli

Khodabakhsh

A homotopy analysis method for limit cycle of the van der Pol oscillator with delayed amplitude limiting. Appl Math Computat 2011; 217: 9404–9411.

30.

Bel

Reartes

Isochronous bifurcations in second-order delay differential equations. Electron J Different Equat 2014; 2014: 1–12.

31.

Cobiaga

Reartes

WA.

New approach in the search for periodic orbits. Int J Bifurcation Chaos 2013; 23: 1350186.

32.

Cobiaga

Reartes

Search for periodic orbits in delay differential equations. Int J Bifurcation Chaos 2014; 24: 1450084.

33.

Nave

Gol'dshtein

Ajadi

Singularly perturbed homotopy analysis method applied to the pressure driven flame in porous media. Combust Flame 2015; 162: 864–873.

34.

Rajeev Singh

AK.

Homotopy analysis method for a fractional Stefan problem. Nonlinear Sci Lett A 2017; 8: 50–59.

35.

Patel M

Desai

. Homotopy analysis method for fingero-imbibition phenomenon in heterogeneous porous medium. Nonlinear Sci Lett A 2017; 8: 90–100.

36.

Singh

Prasad

Som

An analytical algorithm for fractional (1 + 1) dimensional nonlinear Boussinesq equation by the homotopy analysis method. Nonlinear Sci Lett A 2017; 8: 276–288.

37.

Arshad

Siddiqui

Sohail

, et al. Comparison of optimal homotopy analysis method and fractional homotopy analysis transform method for the dynamical analysis of fractional order optical solitons. Adv Mech Eng 2017; 9: 168781401769294.

38.

Fei

Lin

Yan

, et al. Approximate solution of a piecewise linear-nonlinear oscillator using the homotopy analysis method. J Vibrat Contr 2017; 22. DOI:10.1177/1077546317729972

39.

Ebenezer

, et al. The multi-step homotopy analysis method for a modified epidemiological model for measles disease. Nonlinear Sci Lett A 2017; 8: 320–332.

40.

Liu

Yue

Zhou

JL.

Piezoelectric optimal delayed feedback control for nonlinear vibration of beams. J Low Freq Noise, Vibrat Active Contr 2016; 35: 25–28.

41.

Shokouhi

Yuan

Zhu

Optimal placement of sensors and piezoelectric friction dampers in the pipeline networks based on nonlinear dynamic analysis. J Low Freq Noise Vibrat Active Contr 2017; 35: 56–82.

42.

Zhao

Zhou

Yue

, et al. Hybrid modelling of driver seat-cushion coupled system for metropolitan bus. J Low Freq Noise Vibrat Active Contr 2017; 36: 214–226.

43.

Nimafar

, et al. Akbari-Ganji method for vibration under external harmonic loads. Nonlinear Sci Lett A 2017; 8: 416–437.