#### Abstract

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.

#### 1. 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 and 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 , at the critical value of parameters. The double Hopf bifurcation is called nonresonance if is irrational. If is rational, the bifurcation is a strong (low-order) resonance when and or , 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 resonance may undergo saddle-node bifurcation, which leads to so-called mixed modes. Thus, the 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 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

where , and are real constant matrices such that and , represents a nonlinear function in its variable with , 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]:

We found a double Hopf bifurcation with 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.

#### 2. 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 is always an equilibrium point or trivial solution of the system. The characteristic equation at is given by

where is the identity matrix. The characteristic equation (2.1) can be rewritten as

where

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 in (2.2) or , 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 . Thus, 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 but , substituting 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 in (2.4), given by

Eliminating from (2.5), we have

where and when the following conditions hold:

Then, two families of surfaces, denoted by and which correspond to and , respectively, can be derived from (2.4) and are given by

It follows from (2.6) and (2.8) that and are functions of , , and given in (2.3). If there exist values of , , , and such that

where , or , then such point is called the 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 and , 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

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

where is given in (2.11). The corresponding value of the time delay at the strong resonant double Hopf bifurcation point, noted as , is given by

For case (b) with and , , , , , , , the characteristic equation (2.1) becomes

By substituting into (2.15), we get

It follows from (2.16) that is real only when . To obtain the critical boundary, set . This gives

where The necessary conditions for the occurrence of a strong resonant double Hopf bifurcation can be obtained by setting and for and or . 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 and , 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 , and by , and at a double Hopf point, respectively. To investigate the dynamical behavior in a neighborhood of the double Hopf point, we consider , , and for a small , given by

Substituting (2.20) into (1.1) yields

where , and

When , 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 and is determined by (2.12) for case (a) and (2.17) for case (b). Substituting (2.23) and (2.24) into (2.21) when and using the harmonic balance, one obtains that

where , and is the identity matrix. Equations (2.25) and (2.26) are in fact identical. Let and note that . It follows from (2.25) that the periodic solution of (2.21) for is given by

where , for .

In order to obtain an approximate solution of (2.21) when , we now consider a perturbation to (2.27), given by

where , , and is a detuning parameter. The following theorem provides a method to determine in (2.28).

Theorem 2.1. *If is a periodic solution of the equation
**
with , then
*

*Proof. *Multiplying both sides of (2.21) by and integrating with respect to from zero to , one has
where is the period of . Equation (2.31) yields
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 in (2.32). It follows from (2.27) that the periodic solution of (2.29) can been written as

where and are independent constants. Substituting (2.28), and (2.33) into (2.30) and noting the independence of and yield four algebraic equations in , and , which are called the amplitude modulation equations of (2.21). If there exist square terms in (2.21), then plays the role of master amplitude since and would not exist without [10] for the internal resonance where and . Eliminating from the amplitude modulation equations, one can obtain three algebraic equations in , and .

#### 3. Application

As an application, we investigate the delay-induced 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

where , , , and are constants, , , , , time delay, . System (3.1) can be rewritten in the form of (1.1), where and

Substituting (3.2) into (2.3) yields

Thus, (3.1) belongs to case (a) of Section 1 since and . To investigate the strong resonant double Hopf bifurcation, we choose in , in and the delay as the three variable parameters and keep the other parameters fixed. It follows from the analysis in the above section that , and satisfy (2.11), (2.12), and (2.14) at the resonant double Hopf bifurcation point for and , where , , . Correspondingly, and . Any values of , and in the neighborhood of can be expressed as , and . Equation (3.1) can be expressed in the form of (2.21), where , By simple computation, one has , , , . It follows from (2.33) that the solution of (2.29) is given by

Noting that there is a square term in (3.1) when and from discussion at the end of the last section, one may express an approximate solution of (3.1) at as

where . Substituting (3.4) and (3.5) into (2.30) and noting that , , and are independent variables yield the following four algebraic equations in , , , and with three independent bifurcation parameters , , and , given by

Eliminating from (3.6), one has

For , it follows from (3.7) that , denoted as , is always a root. Besides, up to three other roots (in the positive quadrant), denoted as , and , can be obtained analytically with from (3.7) They are given by

where .

It can be seen from (3.7) and (3.8) that , and determine the feature of motions of the system (3.1) when the strong resonant Hopf point at is perturbed by , and for given values of and . Such motion can be amplitude death when is stable, single-mode periodic when or is a stable solution, or double-periodic when there is a nonzero stable solution in . Therefore, it is necessary to classify the solutions of the algebraic equation (3.7) in the neighborhood of the strong double Hopf point . To this end, we represent the critical surfaces in the -space which partition the space into several regions where the number and the stability of the solutions are distinct. In Figure 1, , and are three different surfaces in the -space, defined by

For given values of 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 . The space is separated into six regions noted as (I)–(VI). There are a stable trivial solution and an unstable periodic solution in region (I) which is an amplitude death region. With changing to region (II), the trivial solution loses its stability. Two unstable solutions at (0,0) and exist in region (III). When enters into region (IV), there are three unstable solutions given by (0,0), and . There are a stable nontrivial solution and three unstable solutions , and in region (V), With varying to region (VI), the nontrivial solution disappears and the solution () becomes stable, which results in a stable single-mode 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 with a frequency for a small *.* Similarly, the mixed-mode solution of (3.1) is given by when the parameters under consideration are located in region (V), where and satisfy (3.8) and is a constant. It is seen that the mixed-mode solution consists of two frequencies, that is, and . This implies that the mixed-mode solution is a frequency-doubling or period-doubling solution. Thus, the parameters cross through 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 period-doubling solution occurring in (3.1). Thus, the parameter has strong influence on the critical boundaries. This confirms that the bifurcation at is codimension three.

To have a clearer picture of the analytical prediction of Figure 1, we consider the cross-section at 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 versus . 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). To study the delay-induced bifurcations analytically, we fix the value of at and vary along this vertical dashing line (see Figure 2). The corresponding bifurcation diagram of the analytical solutions derived from (3.5) and (3.8) for in terms of and is shown in Figure 3.

Figure 3 shows that no vibration occurs in the system when the delay is smaller than the value at , that is, . The system undergoes a Hopf bifurcation at and a stable branch of periodic solution in (3.8) bifurcates from the trivial solution at . With the delay increasing, the periodic solution becomes unstable and a period-doubling bifurcation occurs at with a maximum amplitude .

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 in regions displayed in Figure 1. To compare with our analytical predictions mentioned above, we keep , , and the same as those in Figures 1–3. First, we verify the analytical results shown in Figure 3 by keeping and varying . The values of are chosen to be , and , 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 resonance. The analytical results are also in agreement with those in [8].

**(a)**

**(b)**

**(c)**

#### 4. 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 from (3.6) and (3.7) when and for 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 and and the stabilities of the solutions are shown in Figure 5, which are similar to that in [11]. Thus, the Hopf-Hopf interactions with 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 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 . Thus, the bifurcation diagrams of the corresponding truncated normal form in 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 . The corresponding amplitude equation can be also obtained correctly in (3.7) without the square terms or 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 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 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.

#### Acknowledgments

This work is supported by Hong Kong Research Grants Council under CERG Grant (Cityu 1007/05E) and Key Project of National Nature and Science Foundation of China under Grant no. 10532050. The senior author also thanks the support of the National Outstanding Young Funds of China under Grant no. 10625211 and the Program of Shanghai Subject Chief Scientist under Grant no. 08XD14044.