Double Hopf Bifurcation with Strong Resonances in Delayed Systems with Nonlinerities

An efficient method is proposed to study delay-induced strong resonant double Hopf bifurcation for nonlinear systems with time delay. As an illustration, the proposed method is employed to investigate the 1 : 2 double Hopf bifurcation in the van der Pol system with time delay. Dynamics arising from the bifurcation are classified qualitatively and expressed approximately in a closed form for either square or cubic nonlinearity. The results show that 1 : 2 resonance can lead to codimension-three and codimension-two bifurcations. The validity of analytical predictions is shown by their consistency with numerical simulations.


Introduction
To realistically model physical and biological systems with interactions, we generally need to take into account the influence of the past states.Such a system can be modeled by a set of delay differential equations DDEs .Unlike a differential equation, the characteristic equation of a DDE has an infinite number of eigenvalues due to the time delay.This implies that an infinite number of modes may occur in the system governed by a DDE.It is possible that a double Hopf bifurcation is induced by the time delay if the characteristic equation at an equilibrium has eigenvalues with strict negative real parts except a pair of purely imaginary eigenvalues.The study of double Hopf bifurcation has recently attracted particular attention as the bifurcation exhibits rich dynamical behavior and creates complicated motions which reflect the intrinsic properties of the system.In our previous research 1-3 , it was found that even a second-order controlled system with delayed feedback may undergo a double Hopf bifurcation when the time delay and feedback gain are varied.Such phenomenon was also observed by Reddy et al. 4 ,Balanov et al. 5 , and Ma et al. 6 .These researches showed that the time delay in various systems has a qualitative effect on dynamics.Physically, energy between two modes is exchanged when a double Hopf bifurcation occurs in the system under consideration.Therefore, a study of the effect of time delay on the dynamics of mathematical models gives insight into possible mechanism behind the observed behavior.
There are three cases: nonresonances, weak high-order , and strong low-order resonances of the critical frequencies for a double Hopf bifurcation 7 .Let ±iω − and ±iω be the two pairs of purely imaginary eigenvalues of the characteristic equation for a DDE and all other eigenvalues have negative real parts such that ω ≥ ω − > 0, k 1 : k 2 ω − : ω at the critical value of parameters.The double Hopf bifurcation is called nonresonance if k 1 : k 2 is irrational.If k 1 : k 2 is rational, the bifurcation is a strong low-order resonance when k 1 1 and k 2 1, 2, 3, or 4, and a weak high-order resonance otherwise.Many methods and techniques from the geometric theory of dynamical systems on ordinary differential equations may be extended to study the former two cases.However, the problem for the last case becomes complicated since a strong resonance may be a codimension-3 phenomenon which requires three parameters for its unfolding.In fact, one of two single modes induced by the 1 : 2 resonance may undergo saddle-node bifurcation, which leads to so-called mixed modes.Thus, the 1 : 2 resonance is treated as a degenerate double bifurcation or codimension-3 case when the saddle-node bifurcation point and double Hopf bifurcation one coincide in three-parameters space.It is well known that dynamics derived from strong resonance cannot be predicted from the theory of nonresonant double Hopf interactions.The investigation of delay-induced double Hopf bifurcation with strong resonances is both challenging and important.The first successful research for studying 1 : 2 strong resonant double Hopf bifurcation of a DDE was performed by Campbell and LeBlanc 8 in terms of the center manifold reduction CMR .Using CMR for the study of resonant double Hopf bifurcation is tedious.Moreover, the unfolding parameters obtained usually do not carry any physical meaning.However, for the case of DDEs, the unfolding parameters are related to the quantities of a physical system such as time delay and feedback gains in controlled systems.For these issues, we propose an efficient but simple method in 9 to avoid the tedious computation encountered in CMR.The result will also bring out the physical significance of unfolding parameters in higher codimensional bifurcations.
In 9 , we studied weak high-order resonant double Hopf bifurcations which arise from time delay and other physical parameters in a type of two first-order delay differential equations given by

An Efficient Method for Studying Strong Resonant Double Hopf Bifurcations
We follow the idea of 9 in this section.It can be seen from 1.1 that Z 0 is always an equilibrium point or trivial solution of the system.The characteristic equation at Z 0 is given by where I is the identity matrix.The characteristic equation 2.1 can be rewritten as

2.3
The roots of the characteristic equation 2.2 are usually called the eigenvalues of the equilibrium point of system 1.1 .The stability of the trivial equilibrium point will change when the system under consideration has zero or a pair of purely imaginary eigenvalues.The former occurs when λ 0 in 2.2 or c d c 2 det D 0, which leads to a static bifurcation of the equilibrium point.The latter deals with the Hopf bifurcation such that the dynamical behavior of the system changes from a static stable state to a periodic motion or vice versa.The dynamics becomes complicated when the system has two pairs of purely imaginary eigenvalues at a critical value of the time delay.We will concentrate on such cases.For this, we let c d c 2 det D / 0. Thus, λ 0 is not a root of the characteristic equation 2.2 in the present paper.Such assumption can be realized in engineering as long as one chooses a suitable feedback controller.
For case a with det D 0 but c d c 2 / 0, substituting λ a iω into 2.2 , and equating the real and imaginary parts to zero yield One can obtain the explicit expressions for the critical stability boundaries by setting a 0 in 2.4 , given by Eliminating τ from 2.5 , we have where ω − < ω and when the following conditions hold:

2.7
Then, two families of surfaces, denoted by τ − and τ which correspond to ω − and ω , respectively, can be derived from 2.4 and are given by

2.8
It follows from 2.6 and 2.8 that ω ± and τ ± are functions of c 1 , d 1 , c 2 , and c d given in 2.3 .
If there exist values of c 1 , d 1 , c 2 , and c d such that τ − τ , 2.9 where k 1 1, k 2 1, 2, 3, or 4, then such point is called the k 1 : k 2 strong or low-order resonant double Hopf bifurcation point in parameter space.Equations 2.9 and 2.10 not only determine the linearized system around the trivial equilibrium which has two pairs of purely imaginary eigenvalues ±iω − and ±iω , but also give a relation between ω − and ω .Equations 2.9 and 2.10 form the necessary conditions for the occurrence of a strong resonant double Hopf bifurcation point.Equation 2.10 yields when condition 2.7 is satisfied.Substituting 2.11 into 2.6 , one can obtain the frequencies in the simple expressions given by

2.12
The other parameters can be determined by 2.9 or the following equation:

2.13
where d 1 is given in 2.11 .The corresponding value of the time delay at the strong resonant double Hopf bifurcation point, noted as τ c , is given by

2.18
where j 1, 2, . . . .The necessary conditions for the k 1 : k 2 occurrence of a strong resonant double Hopf bifurcation can be obtained by setting τ − j τ j and ω − /ω k 1 /k 2 for k 1 1 and k 2 1, 2, 3 or 4.This yields As mentioned above, a strong resonance may be a codimension-3 phenomenon which requires three parameters for its unfolding.For our purpose, the time delay is always chosen as a bifurcation parameter.The other two variable parameters come from the element of C and D, respectively, while the remaining parameters are kept fixed.Consequently, at a double Hopf point, these three variable parameters satisfy 2.11 , 2.12 and 2.14 for case a , and 2.17 and 2.19 for case b .Correspondingly, denote C, D and τ by C c , D c and τ c at a double Hopf point, respectively.To investigate the dynamical behavior in a neighborhood of the double Hopf point, we consider C , D , and τ for a small , given by where Z Z t , Z τ Z t − τ , and

2.22
When 0, F 0, and a resonant double Hopf bifurcation point occurs in the system under consideration.Thus, the zero-order approximate solution of 2.21 may be expressed as which results in where φ ωt and ω ω − /k 1 ω /k 2 is determined by 2.12 for case a and 2.17 for case b .Substituting 2.23 and 2.24 into 2.21 when 0 and using the harmonic balance, one obtains that where where In order to obtain an approximate solution of 2.21 when / 0, we now consider a perturbation to 2.27 , given by where and σ is a detuning parameter.The following theorem provides a method to determine

Mathematical Problems in Engineering
Proof.Multiplying both sides of 2.21 by W t T and integrating with respect to t from zero to 2π/ω, one has

2.31
where 2π/ω is the period of W t .Equation 2.31 yields

2.32
Equation 2.30 follows from 2.29 and 2.32 and the theorem is proved.
To apply the above theorem to determine , one must know the expression of W t in 2.32 .It follows from 2.27 that the periodic solution of 2.29 can been written as where p k i and q k i are independent constants.Substituting 2.28 , and 2.33 into 2.30 and noting the independence of p k i and q k i yield four algebraic equations in a k i , b k i and σ, which are called the amplitude modulation equations of 2.21 .If there exist square terms

Application
As an application, we investigate the delay-induced 1 : 2 strong resonant double Hopf bifurcation in the van der Pol oscillator with delayed position feedback and show the efficiency and simplicity of the above method.The oscillator is governed by Thus, 3.1 belongs to case a of Section 1 since det D 0 and c d c 2 1 / 0. To investigate the 1 : 2 strong resonant double Hopf bifurcation, we choose α in C, A in D and the delay τ as the three variable parameters and keep the other parameters fixed.It follows from the analysis in the above section that α, A, and τ satisfy 2.11 , 2.12 , and 2.14 at the 1 : 2 resonant double Hopf bifurcation point for k Noting that there is a square term in 3.1 when β 1 / 0 and from discussion at the end of the last section, one may express an approximate solution of 3.1 at O as where ϕ 1/2 σ.Substituting 3.4 and 3.5 into 2.30 and noting that p 1 , p 2 , q 1 and q 2 are independent variables yield the following four algebraic equations in r 1 , a 2 , b 2 , and σ with three independent bifurcation parameters A , τ , and α , given by

3.7
For β 2 0, it follows from 3.7 that 0, 0, 0 , denoted as r 10 , a 20 , b 20 , is always a root.Besides, up to three other roots in the positive quadrant , denoted as r 11 , r 20 , r 10 , r 21 , and r 12 , r 22 , can be obtained analytically with r 2 a 2 2 b 2 2 1/2 from 3.7 They are given by It can be seen from 3.7 and 3.8 that r 1 • , r 2 • , and σ determine the feature of motions of the system 3.1 when the strong resonant Hopf point at A c , α c , τ c is perturbed by α , A , and τ for given values of β 1 and γ.Such motion can be amplitude death when r 10 , a 20 , b 20 is stable, single-mode periodic when r 11 , r 20 or r 10 , r 21 is a stable solution, or double-periodic when there is a nonzero stable solution in r 12 , r 22 .Therefore, it is necessary to classify the solutions of the algebraic equation 3.7 in the neighborhood of the strong double Hopf point A c , α c , τ c .To this end, we represent the critical surfaces in the A, α, τ -space which partition the space into several regions where the number and the stability of the solutions are distinct.In Figure 1, S 1 , S 2 , and S 3 are three different surfaces in the A, α, τ -space, defined by A, α, τ | r 12 0, r 22 0, σ / 0 .

3.9
For given values of β 1 and γ, the stability of the solutions on these regions may be obtained by observing the eigenvalues of the Jacobian matrix of 3.8 at the corresponding solution.In Figure 1, a three-dimensional view of the asymptotic boundaries is shown for β 1 γ 0.1.The space is separated into six regions noted as I -VI .There are a stable trivial solution r 10 , r 20 0, 0 and an unstable periodic solution r 10 , r 21 in region I which is an amplitude death region.With A, α, τ changing to region II , the trivial solution loses its stability.Two unstable solutions at 0,0 and r 11 , 0 exist in region III .When A, α, τ enters into region IV , there are three unstable solutions given by 0,0 , r 11 , 0 and 0, r 21 .There are a stable nontrivial solution r 12 , r 22 and three unstable solutions 0, 0 , 0, r 21 and r 11 , 0 in region V , With A, α, τ varying to region VI , the nontrivial solution r 12 , r 22 disappears and the solution r 11 , 0 becomes stable, which results in a stable singlemode or periodic solution.It follows form 3.5 and 3.8 that the stable single-mode or periodic solution of 3.1 is expressed as u t 8α − 3τ − 24A π /2γ cos 1/2 σ t with a frequency 1/2 σ for a small .Similarly, the mixed-mode solution of 3.1 is given by u t r 12 cos 1/2 σ t r 22 cos 2 1/2 σ t θ when the parameters under consideration are located in region V , where r 12 and r 22 satisfy 3.8 and θ is a constant.It is seen that the mixed-mode solution consists of two frequencies, that is, 1/2 σ and 2 1/2 σ .This implies that the mixed-mode solution is a frequency-doubling or perioddoubling solution.Thus, the parameters cross through S 3 from region VI to region V , the solution of 3.1 undergoes a period-doubling bifurcation.Different from the nonresonant case 11 , the critical boundaries have strong curvatures when the solution of the single-mode amplitude equation undergoes a saddle-node bifurcation, which corresponds to a perioddoubling solution occurring in 3.1 .Thus, the parameter α has strong influence on the critical boundaries.This confirms that the bifurcation at A c , α c , τ c is codimension three.
To have a clearer picture of the analytical prediction of Figure 1, we consider the cross-section at α 0 which is shown in Figure 2, in which the number and the stability of the analytical solutions in regions I -VI are also displayed by r 1 versus r 2 .It is seen in Figure 2 that there are a stable trivial, periodic and doubling periodic solutions in regions I , VI and V , respectively.No stable solution exists in regions II , III and IV .Figure 3 shows that no vibration occurs in the system when the delay is smaller than the value at P 1 , that is, τ τ c 2π.The system undergoes a Hopf bifurcation at P 1 and a stable branch of periodic solution r 11 in 3.8 bifurcates from the trivial solution at P 1 .With the delay increasing, the periodic solution becomes unstable and a period-doubling bifurcation occurs at P 2 with a maximum amplitude r 2 12 r 2 22 .Now, numerical simulation is employed to qualitatively examine the validity of the present method.The Runge-Kutta scheme is adopted to produce the numerical results of 3.1 for α, A, τ in regions displayed in Figure 1.To compare with our analytical predictions mentioned above, we keep β 2 , β 1 , and γ the same as those in Figures 1-3.First, we verify the analytical results shown in Figure 3 by keeping A −0.36, and varying τ.The values of τ are chosen to be 6.28, 6.285, and 6.295, respectively, as shown in Figures 4 a , 4 b and 4 c . Figure 4 shows that there are a stable trivial, periodic and period-doubling solutions when τ is located in regions I , VI and V , respectively.The numerical simulation also shows no stable solutions in the other regions.Therefore, the analytical prediction is in good agreement with the numerical result, which implies that solution 3.5 from the present method is valid qualitatively.In addition, the period-doubling bifurcation is a direct consequence of the proximity to the 1 : 2 resonance.The analytical results are also in agreement with those in 8 .

Discussion and Conclusion
We have proposed an analytical method to consider various delay-induced solutions derived from double Hopf bifurcation with strong resonance for a type of two first-order delay differential equations.It should be noted that the present method can also be applied to the weak and nonresonant cases.For example, one cannot solve any nontrivial solutions with σ / 0 from 3.6 and 3.7 when β 1 0 and β 2 / 0 for 1 : 2 resonance.This implies that there is no solution in the form of 3.5 .However, the solution can be expressed in the following form: where ϕ 1 1/2 σ 1 and ϕ 2 1 σ 2 and the stabilities of the solutions are shown in Figure 5, which are similar to that in 11 .Thus, the Hopf-Hopf interactions with 1 : 2 resonance become weak.
The results provided in the present paper can be also understood further by comparing those given by Campbell and Leblanc 8 , where they use a center manifold reduction and normal form analysis to deduce a four-dimensional invariant center manifold, the flow on which is a good approximation for the long-term behavior of solutions to the original differential equation.Correspondingly, the truncated amplitude equation on the center manifold is presented in O 2 providing that square terms are not degenerate.Such truncated amplitude equation is also given in the present work, as expressed in 3.7 when the cubic terms are neglected, or β 2 γ 0. Thus, the bifurcation diagrams of the corresponding truncated normal form in O 2 are shwon in Figure 2. When the square terms in the truncated normal form are degenerate, one has to consider the truncated normal form in O 3 .The corresponding amplitude equation can be also obtained correctly in 3.7 without the square terms or β 1 0 in 3.1 .For this case, the dynamical behavior near the double Hopf point is distinct topologically, as shown in Figure 5. Consequently, one can distinguish two cases for strong and weak resonances.
Finally, a brief conclusion is represented to end this paper.An efficient method is proposed to study delay-induced strong resonant double Hopf bifurcation for nonlinear systems with time delay.The method is first described and leads to the main theorem.As an application, the proposed method is employed to investigate the double Hopf bifurcation with 1 : 2 strong resonance in the van der Pol system with time delay.Dynamics arising from the bifurcation are classified qualitatively and expressed approximately in a closed form for either square or cubic nonlinearity.The results show that 1 : 2 resonance can lead to codimension-three and codimension-two bifurcations.The validity of analytical predictions is illustrated by their consistency with numerical simulations.As a result, a simple but efficient method is developed for the classification of unfolding of higher codimensional bifurcations in DDEs.The present research brings out the physical meaning and significance of unfolding parameters in higher codimensional bifurcations.

c 11 c 12 c 21 c 22 and D d 11 d 12 d 21 d 22 ,F
C and D are 2×2 real constant matrices such that C represents a nonlinear function in its variable with F 0, 0 0, is a small parameter representing strength of nonlinearities, and τ the time delay.In addition, we restricted our research to the following two cases in 9 : a det D 0, c 22 d 11 − c 21 d 12 − c 12 d 21 c 11 d 22 c 11 c 22 − c 12 c 21 / 0; b det D / 0, c 11 c 22 , c 12 −c 21 , d 11 d 22 , d 12 d 21 0, c 11 > 0, c 12 > 0, d 11 / 0. We found a double Hopf bifurcation with 1 : 2 strong or low-order resonance occurring in 1.1 .The strong resonant bifurcation was left open since it requires three parameters for the unfolding.In this paper, we extend our investigation to study this open problem.

in 2 . 2 k 1 b 2 k 1 plays
21 , then r 1 a the role of master amplitude since a k 2 , and b k 2 would not exist without r 1 10 for the 1 : 2 internal resonance where k 1 1 and k 2 2. Eliminating σ from the amplitude modulation equations, one can obtain three algebraic equations in r 1 , a 2 , and b 2 .

1 1 and k 2 2 , 1 − 1 0 − 1 1 1 .
where α c 0, A c −3/8, τ c 2π. Correspondingly, ω − 1/2 and ω 1.Any values of A, τ and α in the neighborhood of A c , τ c , α c can be expressed as A be expressed in the form of 2.21 , where C c 0 τc τ −x τc A x τc τ −x α y−x 2 β 1 β 2 x γy .By simple computation, one has N k It follows from 2.33 that the solution of 2.29 is given by
and I is the 2 × 2 identity matrix.