ARAB JOURNAL OF Arab J Math Sci 20(2) (2014), 280–301

MATHEMATICAL SCIENCES

Double Hopf bifurcation in delay differential equations

a,* b
REDOUANE QESMI , MOHAMED AIT BABRAM
a
Ecole Supérieure de Technologie, Sidi Mohamed Ben Abdellah University,
Fez 30000, Morocco
b
Faculté des Sciences et Téchniques Guéliz, Cadi Ayyad University,
Marrakech, Morocco

Received 7 October 2012; revised 15 October 2013; accepted 16 October 2013

Available online 1 November 2013
Abstract. The paper addresses the computation of elements of double Hopf bifurca-
tion for retarded functional differential equations (FDEs) with parameters. We present
an efficient method for computing, simultaneously, the coefficients of center manifolds
and normal forms, in terms of the original FDEs, associated with the double Hopf
singularity up to an arbitrary order. Finally, we apply our results to a nonlinear model
with periodic delay. This shows the applicability of the methodology in the study of
delay models arising in either natural or technological problems.

Keywords: Double Hopf; Delay; Bifurcation; Functional differential equation; Center

manifold; Normal forms; Regenerative cutting tool

1. INTRODUCTION

The double Hopf bifurcation, which is readily located via linear stability analysis as the
codimension-two point at which two pairs of complex conjugate eigenvalues have their real
part simultaneously change sign, has associated very rich nonlinear dynamics in its neigh-
borhood. Depending on the particulars of the system under consideration, there are
around thirty different dynamical scenarios, divided into simple and difficult cases. Center
manifold theory is of fundamental importance in the study of nonlinear dynamical systems
when analyzing bifurcations of such a type. In fact, this theory allows us to reduce the study
of a differential equation with delay near to a non-hyperbolic equilibrium point to that of
an ordinary differential equation on a finite-dimensional invariant manifold. This

* Corresponding author. Tel.: +212 416 900 1978.

E-mail address: qesmir@gmail.com (R. Qesmi).
Peer review under responsibility of King Saud University.

Production and hosting by Elsevier

1319-5166 ª 2013 King Saud University. Production and hosting by Elsevier B.V. All rights reserved.
http://dx.doi.org/10.1016/j.ajmsc.2013.10.002
 Double Hopf bifurcation in delay differential equations 281

approach is particularly interesting when the starting point is an infinite-dimensional prob-

lem, such as a functional, partial or an integro differential equation. In fact, the reduction
forms also a qualitative simplification. Although center manifolds exist, they need not be
unique, there are some points which must always be on any center manifold: all the center
manifolds have the same Taylor expansion near the equilibrium up to their order of reg-
ularity; such an expansion may give an idea of the local dynamics near the steady state.
Center manifolds and normal forms of such bifurcations have been discussed in [5,9] for
ordinary differential equations. Center manifolds for Hopf and fold-Hopf bifurcations in
FDEs with discrete delays have been studied in [6,7]. In [1,2], Faria and Magalhaès have
considered the computation of coefficients of normal forms for functional differential
equations (FDEs) of both Hopf and Bogdanov singularities. However, it is difficult to
apply the method to compute the explicit expressions for these coefficients since it de-
mands much more computation efforts for high-order normal forms.
The attention of this paper will be focused on the development of methodology and
software for computing the center manifolds and normal forms of double Hopf bifur-
cation for general FDEs with parameters. According to the structure of linearized
equation of a retarded system evaluated at an equilibrium, the case considered in this
paper corresponds to two different pairs of purely imaginary eigenvalues.
The aim is twofold: the primary objective of this paper is to prove that the terms of a
parameterized center manifold associated with double Hopf singularity up to the order of
regularity for FDEs satisfy an efficient and explicit recursive formulas, and to give an effi-
cient algorithm. The purpose of the second fold is to obtain, for the general situation of
double Hopf bifurcation for FDEs, explicit formulas giving the coefficients of normal
forms in terms of the coefficients of the original equation. This allows a simpler scheme
that greatly saves computational time and computation memory. Furthermore, the ob-
tained method combines center manifolds and normal form schemes into one step. To
show the applicability and the efficiency of the method, we apply our algorithmic scheme
to obtain normal forms of a nonlinear mechanical model with periodic delay. We achieve
this by augmenting the explicit-time dependent delay terms as new state variables to the
original equations of motions with appropriate initial conditions and applying the meth-
od of computing the center manifolds and normal forms obtained in the present paper.
For the notation, background about the theory of FDEs and all needed results in the
remainder of this paper, we follow [9]. However, we use Cn ¼ Cð½r; 0; Rn Þ; r P 0 since
we need to work in realization spaces with different dimensions, depending on whether
the parameters are incorporated or not incorporated in the realization space variables.
The paper is organized as follows: a computational scheme for a center manifold
associated with double Hopf with parameters are given in Section 2. Normal forms
for FDEs are derived, in terms of the original equation and the center manifold coef-
ficients, in the same section. To illustrate our results, we study the bifurcation of a
model arising from the mechanics in Section 3. Conclusions are drawn in Section 4.

2. MAIN RESULTS

In this section, we present our result concerning the computation of terms of center
manifolds and normal forms for FDEs of the form
282 R. Qesmi, M.A. Babram

_ ¼ LðaÞxt þ Fðxt ; aÞ;

xðtÞ ð2:1Þ
p 1
where a 2 R ; a # LðaÞ is a C function with values in the space of bounded linear
operators from Cn :¼ Cð½r; 0; Rn Þ to Rn . F : Cn  Rp ! Rn is assumed, without loss
of generality, to be sufficiently smooth ðF 2 C1 Þ with Fð0; aÞ ¼ 0 and DFð0; aÞ ¼ 0
for all a 2 Rp .
We denote L0 ¼ Lð0Þ. In the sequel of this paper, we assume that the following
hypothesis is satisfied:

(H): The linear equation x_ ðtÞ ¼ L0 xt has two purely imaginary pairs ðix1 Þ and ðix2 Þ
as eigenvalues and no other characteristic values with zero real part.

One way of considering center manifolds for differential equations with parameters
is to reduce the situation to the case of differential equations without parameters by
considering the system
_ ¼ L0 xt þ ½LðaÞ  L0 xt þ Fðxt ; aÞ;
xðtÞ ð2:2Þ
_ ¼ 0:
aðtÞ ð2:3Þ
The solutions of this system are of the form x~ðtÞ :¼ ðxðtÞ; aðtÞÞ 2 Rnþp , the phase space
e :¼ Cnþp , and the system can be written as
is C
x~_ ðtÞ ¼ Lex~t þ Fð~
e xt Þ; ð2:4Þ
where Lðu; e vÞ :¼ ðL0 u; 0Þ and Fðu; e vÞ :¼ ð½Lðvð0ÞÞ  L0 u þ Fðu; vð0ÞÞ; 0Þ with
e
u 2 Cn ; v 2 Cp . Let A and A denote the infinitesimal generators associated with the
_ ¼ L0 xt and xðtÞ
equations xðtÞ _ ¼ Lx e t , respectively. The equation aðtÞ
_ ¼ 0 has a unique
characteristic value k ¼ 0 with multiplicity p. The associated generalized eigenspace
consists of the elements of Cp which are constant functions, and it is denoted here also
by Rp . Let K :¼ rðAÞ \ iR and consider the spectral decomposition Cn ¼ Xc  Xs as in
[9]. In particular, we consider bases for Xc and Xc denoted by
U ¼ ð/1 ;    ; /4 Þand W ¼ colðw1 ;    ; w4 Þ, respectively, and satisfying ðW; UÞ ¼ I. We
define K e ¼ K Wf0g; X ec ¼ Xc  Rp ; and X es ¼ Xs  R where R ¼ fv 2 Cp : vð0Þ ¼ 0g.
e e 
As bases of X c and X c , we consider, respectively, the columns of the matrix U e and
the rows of the matrix W, e
" # " #
U 0 W 0
e ¼
U ; e ¼
W ;
0 Ip 0 Ip

e Ui
which satisfies h W; e ¼ Inþp , where h; i is the bilinear form in C
e  C
e introduced in
[9] and defined by
Z 0 Z h
< w; / >¼ wð0Þ/ð0Þ  wðn  hÞdgðhÞ/ðnÞdn:
r 0

e_ ¼ U
and U e Be with
 Double Hopf bifurcation in delay differential equations 283

2 3
ix1 0 0 0 0
6 7
6 0 ix1 0 0 07
6 7
6 7
Be ¼ 6
6 0 0 ix2 0 07
7:
6 7
6 0 0 0 ix2 07
4 5
0 0 0 0 0p
Then, we obtain the decomposition C e¼X ec  Xes , where X e
ec is the invariant space of A
e
associated with K.
In the sequel, we recall the definition of a local center manifold associated to Eq.
(2.4).

Definition 2.1 [9]. Given a C1 map h~ from R4þp into X es , the graph of h~ is said to be a
local center manifold associated to Eq. (2.4) if and only if hð0Þ ~ ¼ 0; Dhð0Þ ~ ¼ 0, and
4þp
there exists a neighborhood V of zero in R such that, for each n 2 V, there exist
d ¼ dðnÞ > 0 and a function x defined on   d  r; d½ such that x0 ¼ Un e þ hðnÞ~ and x
verify Eq. (2.4) on   d; d½ and satisfy the identity

e zðtÞ þ hð~
xt ¼ U~ ~ zðtÞÞ; for t 2 ½0; d½;

where z~ðtÞ is the unique solution of the ordinary differential equation (ODE)
(
d
z~ðtÞ ¼ B~ e
ezðtÞ þ Wð0Þ Fð e zðtÞ þ hð~
e U~ ~ zðtÞÞÞ;
dt
ð2:5Þ
z~ð0Þ ¼ n; n 2 R4þp :

~ zÞ ¼ ðhð~
Remark 2.2. (i) If we write hð~ zÞÞ; z~ ¼ ðz; aÞ for z 2 R4 and a 2 Rp , then
zÞ; h0 ð~
Eq. (2.5) is equivalent to
   
d z Bz
¼
dt a 0
!
Wð0Þ½ðLðað0Þ þ h0 ðz;aÞð0ÞÞ  L0 ÞðUz þ hðz;aÞÞ þ FðUz þ hðz;aÞ;að0Þ þ h0 ðz;aÞð0ÞÞ
þ
0

with
2 3
ix1 0 0 0
6 7
6 0 ix1 0 0 7
6 7
B¼6 7
6 0 0 ix2 0 7
4 5
0 0 0 ix2
and zð0Þ ¼ n 2 R4 . Noting that h0 ðzðtÞ; aÞð0Þ ¼ 0 and dropping the auxiliary equations
introduced for handling the parameter, we get the reduced equation on the center man-
ifold – represented as a graph over the n and a variables of hðn; aÞ for sufficiently small
n and a – associated with the parameterized Eq. (2.1) as
284 R. Qesmi, M.A. Babram

(
d
dt
zðtÞ ¼ BzðtÞ þ Wð0Þ½ðLðaÞ  L0 ÞðUzðtÞ þ hðzðtÞ;aÞÞ þ FðUzðtÞ þ hðzðtÞ;aÞ;aÞ;
4
zð0Þ ¼ n; n2R :
ð2:6Þ
(ii) It is noted that the invariance properties of center manifolds guarantee that any
small solutions bifurcating from ð0; 0; 0Þ must lie in any center manifold and thus we
may follow the local evolution of bifurcating families of solutions in this suspended
family of center manifolds (see [18] for more details).

2.1. Calculation of center manifolds

In the following, we give an analytic characterization of a center manifold associated to

Eq. (2.1).

Theorem 2.3. Given a C1 map h from R4þp into Xs with hð0Þ ¼ 0 and Dhð0Þ ¼ 0, a
necessary condition for the graph of h to be a local center manifold of Eq. (2.1) is that
there exists a neighborhood V of zero in R4þp such that, for each ðn; aÞ 2 V,

@ @hðn; aÞ @hðn; aÞ @hðn; aÞ

ðhðn; aÞÞðhÞ ¼ ix1 ðhÞn1 þ ix1 ðhÞn2  ix2 ðhÞn3
@h @n1 @n2 @n3

@hðn; aÞ @hðn; aÞ
þ ix2 ðhÞn4 þ ðhÞWð0Þ½ðLðaÞ  Lð0ÞÞðUn þ hðn; aÞ
@n4 @n

@hðn; aÞ
þ ðhÞWð0ÞFðUn þ hðn; aÞ; aÞ þ UðhÞWð0Þ½ðLðaÞ
@n

 Lð0ÞÞðUn þ hðn; aÞ þ UðhÞWð0ÞFðUn þ hðn; aÞ; aÞ; ð2:7Þ

@
ðhðn; aÞÞð0Þ ¼ L0 hðn; aÞ þ ðLðaÞ  L0 ÞðUn þ hðn; aÞÞ þ FðUn þ hðn; aÞ; aÞ:
@h
ð2:8Þ

Proof. Let h be a graph of a local center manifold of Eq. (2.1). From Definition 2.1,
there exists a neighborhood V of zero in R4þp such that, for each ðn; aÞ 2 V, there exists
d > 0 such that the solution of (2.1) with initial data Un þ hðn; aÞ exists in the interval
  d  r; d½ and is given by

xt ¼ UzðtÞ þ hðzðtÞ; aÞ; for t 2  d; d½;

such that zðtÞ is a solution of the equation

(
d
dt
zðtÞ ¼ BzðtÞ þ Wð0Þ½ðLðaÞ  L0 ÞðUzðtÞ þ hðzðtÞ; aÞÞ þ FðUzðtÞ þ hðzðtÞ; aÞ; aÞ;
zð0Þ ¼ n; n 2 R4 ;
 Double Hopf bifurcation in delay differential equations 285

2 3
ix1 0 0 0
6 0 ix1 0 0 7
with B ¼ 6
4 0
7.
0 ix2 0 5
0 0 0 ix2
The variation of constants formula of Eq. (2.1) can be written as
Z t
xt ¼ TðtÞ/ þ Tðt  sÞX0 ½ðLðaÞ  L0 Þxs þ Fðxs ; aÞds; t P 0;
0

_ ¼ L0 xt . It follows
where ðTðtÞÞtP0 is the semi-group solution of the linear equation xðtÞ
from the decomposition of the phase space K : Cn ¼ Xc  Xs that the function h
satisfies
Z t
hðzðtÞ; aÞ ¼ TðtÞhðzðtÞ; aÞ þ Tðt  rÞXs0 FðUzðrÞ þ hðzðrÞ; aÞ; aÞdr:
0

Then
1 1
½TðtÞhðn; aÞ  hðn;aÞ ¼ ½hðzðtÞ;aÞ  hðzð0Þ;aÞ
t t
Z
1 t
 Tðt  rÞXs0 ½ðLðaÞ  L0 ÞðUzðrÞ þ hðzðrÞ;aÞÞ þ FðUzðrÞ þ hðzðrÞ;aÞ;aÞdr;
t 0
which implies, from the fact that h and F are smooth and TðÞ is a strongly continuous
semi-group on the Banach space Cn , that hðnÞ is in the domain of A, and
 
@hðn; aÞ
Ahðn; aÞ ¼ fBn þ Wð0Þ½ðLðaÞ  L0 ÞðUn þ hðn; aÞÞ þ FðUn
@n
þ hðn; aÞ; aÞg  Xs0 ½ðLðaÞ  L0 ÞðUn þ hðn; aÞÞ þ FðUn þ hðn; aÞ; aÞ:
Consequently, by evaluating the above equation at h – 0, we have
 
@ @hðn; aÞ
ðhðn; aÞÞðhÞ ¼ ðhÞ fBn þ Wð0Þ½ðLðaÞ  L0 ÞðUn þ hðn; aÞÞ
@h @n
þ FðUn þ hðn; aÞ; aÞg þ UðhÞWð0Þ½ðLðaÞ  L0 ÞðUn
þ hðn; aÞÞ þ FðUn þ hðn; aÞ; aÞ; ð2:9Þ
which is the formula (2.7) of theorem.
On the other hand, it results from the fact that the semi-flow
t # xt ¼ UzðtÞ þ hðzðtÞ; aÞ exists on the open   d; d½ such that for h 2  d; 0,

d d
ðUðhÞn þ hðn; aÞðhÞÞ ¼ x0 ðhÞ
dh dh
d
¼ xðhÞ
dh
¼ LðaÞðUzðhÞ þ hðzðhÞ; aÞÞ þ FðUzðhÞ þ hðzðhÞ; aÞ; aÞ:
Consequently, by the fact that d
dh
Uð0Þn ¼ Uð0ÞBn ¼ L0 ðUnÞ, we obtain
286 R. Qesmi, M.A. Babram

@
ðhðn; aÞÞð0Þ ¼ L0 hðnÞ þ ðLðaÞ  L0 ÞðUn þ hðn; aÞÞ þ FðUn þ hðn; aÞ; aÞ; ð2:10Þ
@h
which completes the proof of theorem. h
Let us recall that the function h which represents the center manifold for Eq. (2.1)
has the same regularity as the nonlinearity F. From this fact and in view of the assumed
smoothness on F, for all m 2 N, we can write
X
m
hðn; aÞ ¼ hk ðn; aÞ þ vðn; aÞ for n 2 V; ð2:11Þ
k¼2
m
where hk is the homogeneous part of degree k and vðn; aÞ ¼ oðjðn; aÞj Þ.
Let k 2 N; k P 2. The homogeneous parts of degree k of Eqs. (2.7) and (2.8) are
respectively given by
@ @hk ðn; aÞ @hk ðn; aÞ @hk ðn; aÞ
ðhk ðn; aÞÞðhÞ ¼ ix1 n1 þ ix1 n2  ix2 n3
@h @n1 @n2 @n3
@hk ðn; aÞ k1
þ ix2 n4 þ N ðn; aÞ; ð2:12Þ
@n4

@
ðhk ðn; aÞÞð0Þ ¼ L0 hk ðn; aÞ þ Rk1 ðn; aÞ; ð2:13Þ
@h
where

k1
X
k1
@hkjþ1 ðn; aÞ
N ðn; aÞ ¼ UWð0ÞRk1 ðn; aÞ þ Wð0ÞRj1 ;
j¼2
@n

and R is the homogeneous part of degree i of Rðn; aÞ ¼ ðLðaÞ  L0 Þ ðUn þ hðn; aÞÞ
i1

þFðUn þ hðn; aÞ; aÞ. In particular, R1 is the homogeneous part of degree 2 of

ðLðaÞ  L0 ÞUn þ FðUn; aÞ which is independent from the terms of a center manifold.
l
If n ¼ ðn1 ; n2 ; n3 ; n4 Þ 2 R4P ; q ¼ ðq1 ; q2 ; q3 ; q4 Þ; a ¼ ða1 ;    ; ap Þ; al ¼ al11 al22    app for
p p
l ¼ ðl1 ;    lp Þ 2 N , and jlj ¼ i¼1 li , then we can write
X k q q q q
hk ðn; aÞ ¼ hðq;lÞ n11 n22 n33 n44 al ; for some hkðq;lÞ 2 Xs ;
ðq;lÞ2Dk
k1
X q q q q
ð2:14Þ
N ðn; aÞðhÞ ¼ Nk1 4 l
ðq;lÞ n1 n2 n3 n4 a ;
1 2 3
for some Nkðq;lÞ 2 Xs ;
ðq;lÞ2Dk

and
X q q q q n
Rk1 ðn; aÞ ¼ Rk1 4 l
ðq;lÞ n1 n2 n3 n4 a ;
1 2 3
for some Rk1
ðq;lÞ 2 R ; ð2:15Þ
ðq;lÞ2Dk

where Dk ¼ fðq; lÞ 2 N4  Np : jðq; lÞj ¼ kg.

For simplicity, we put Wq ¼ x1 ðq1  q2 Þ  x2 ðq3  q4 Þ.

Theorem 2.4. Assume that (H) holds. Then the vector of the coefficients of the
homogeneous part of degree k of a local center manifold associated with Eq. (2.1) is given
in a unique way by the following recursive formulas:
 Double Hopf bifurcation in delay differential equations 287

For k ¼ 2; ðq1 ; q2 ; q3 ; q4 ; lÞ 2 D2 :

h2ðq1 ;q2 ;q3 ;q4 ;lÞ ðhÞ


1
¼ eiWq h h2ðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ þ ðeiWq h  eix1 h Þ/1 ð0Þw1 ð0Þ
iðWq þ x1 Þ
1 1
þ ðeiWq h  eix1 h Þ/2 ð0Þw2 ð0Þ þ ðeiWq h  eix2 h Þ/3 ð0Þw3 ð0Þ
iðWq  x1 Þ iðWq þ x2 Þ
1 1
þ ðeiWq h  eix1 h Þ/2 ð0Þw2 ð0Þ þ ðeiWq h  eix2 h Þ/3 ð0Þw3 ð0Þ
iðWq  x1 Þ iðWq þ x2 Þ

1
þ ðeiWq h  eix2 h Þ/4 ð0Þw4 ð0Þ R1ðq1 ;q2 ;q3 ;q4 ;lÞ : ð2:16Þ
iðWq  x2 Þ
For k > 2 and ðq1 ; q2 ; q3 ; q4 ; lÞ 2 Dk :
Z h
hkðq1 ;q2 ;q3 ;q4 ;lÞ ðhÞ ¼ eiWq h hkðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ þ eiWq ðhsÞ Nk1
ðq1 ;q2 ;q3 ;q4 ;lÞ ðsÞds; ð2:17Þ
0

where the vectors hkðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ; ðq1 ; q2 ; q3 ; q4 ; lÞ 2 Dk are given by solving the following
systems:
1
for Wq R fx1 ; x2 g; hkðq1 ;q2 ;q3 ;q4 lÞ ð0Þ ¼ ½DðiWq Þ Ek1
ðq1 ;q2 ;q3 ;q4 ;lÞ ; ð2:18Þ
! !
hkðq1 ;q2 ;q3 ;lÞ ð0Þ Ek1
ðq1 ;q2 ;q3 ;lÞ
for Wq ¼ x1 ; M1 ¼ ; ð2:19Þ
0 vk1
ðq1 ;q2 ;q3 ;lÞ

! !
hkðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ Ek1
ðq1 ;q2 ;q3 ;q4 ;lÞ
for Wq ¼ x1 ; M1 ¼ ; ð2:20Þ
0 vk1
ðq1 ;q2 ;q3 ;q4 ;lÞ

! !
hkðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ Ek1
ðq1 ;q2 ;q3 ;q4 ;lÞ
for Wq ¼ x2 ; M2 ¼ ; ð2:21Þ
0 vk1
ðq1 ;q2 ;q3 ;q4 ;lÞ

and
! !
hkðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ Ek1
ðq1 ;q2 ;q3 ;q4 ;lÞ
for Wq ¼ x2 ; M2 ¼ : ð2:22Þ
0 vk1
ðq1 ;q2 ;q3 ;q4 ;lÞ

M1 and M2 are ðn þ 1Þ  ðn þ 1Þ matrices defined by

!
Dðx1 iÞ w> 1 ð0Þ
M1 ¼
hw1 ; ex1 i: i 0
and
!
Dðx2 iÞ w> 3 ð0Þ
M2 ¼ :
hw3 ; ex2 i: i 0
288 R. Qesmi, M.A. Babram

Ek1 k1
ðq1 ;q2 ;q3 ;lÞ and vðq1 ;q2 ;q3 ;lÞ are vectors given by means of the coefficients of the center

manifolds already computed (see (2.26) and (2.28) in the proof below).

Proof. Let ðq1 ; q2 ; q3 ; q4 ; lÞ 2 Dk , from Eq. (2.12), we have

dhkðq1 ;q2 ;q3 ;q4 ;lÞ ðhÞ

¼ ½ix1 ðq1  q2 Þ  ix2 ðq3  q4 Þhkðq1 ;q2 ;q3 ;q4 ;lÞ ðhÞ þ Nk1
ðq1 ;q2 ;q3 ;q4 ;lÞ ðhÞ; ð2:23Þ
dh

or equivalently
Z h
hkðq1 ;q2 ;q3 ;q4 ;lÞ ðhÞ ¼ eiWq h hkðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ þ eiðhsÞWq Nk1
ðq1 ;q2 ;q3 ;q4 ;lÞ ðsÞds: ð2:24Þ
0

However, for k ¼ 2, we have that

N1 ðn; aÞðhÞ ¼ UðhÞWð0ÞR1 ðn;aÞ
 
¼ eix1 h /1 ð0Þw1 ð0Þ þ eix1 h /2 ð0Þw2 ð0Þ þ eix2 h /3 ð0Þw3 ð0Þ þ eix2 h /4 ð0Þw4 ð0Þ R1 ðn; aÞ:

It follows that
Z h
eiWq ðhsÞ N1ðq1 ;q2 ;q3 ;q4 ;lÞ ðsÞds
0
Z h Z h
¼ eiWq h eiðWq þx1 Þs ds/1 ð0Þw1 ð0Þ þ eiðWq x1 Þs ds/2 ð0Þw2 ð0Þ
0 0
Z h Z h 
iðWq þx2 Þs iðWq x2 Þs
þ e ds/3 ð0Þw3 ð0Þ þ e ds/4 ð0Þw4 ð0Þ R1ðq1 ;q2 ;q3 ;q4 ;lÞ ;
0 0

which implies, by discussing the cases whether Wq  x1;2 ¼ 0 or not, the value of
h2ðq1 ;q2 ;q3 ;q4 ;lÞ ðhÞ is given in theorem.
The Eq. (2.13) is equivalent to

dhkðq1 ;q2 ;q3 ;q4 ;lÞ

ð0Þ ¼ L0 hkðq1 ;q2 ;q3 ;q4 ;lÞ þ Rk1
ðq1 ;q2 ;q3 ;q4 ;lÞ ;
dh

which gives by evaluating the relation (2.23) at h ¼ 0 that

L0 hkðq1 ;q2 ;q3 ;q4 ;lÞ þ Rk1 k k1
ðq1 ;q2 ;q3 ;q4 ;lÞ ¼ iWq hðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ þ Nðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ:

On the other hand, by taking the values of the linear operator L0 on both sides of Eq.
(2.24), we have
Z  
L0 hkðq1 ;q2 ;q3 ;q4 ;lÞ ¼ L0 ðeiWq  Þhkðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ þ L0 eiðsÞWq Nk1
ðq1 ;q2 ;q3 ;q4 ;lÞ ðsÞds :
0

Consequently, we obtain from the expression of the characteristic equation that

DðiWq Þhkðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ ¼ Ek1

ðq1 ;q2 ;q3 ;q4 ;lÞ for ðq1 ; q2 ; q3 ; q4 ; lÞ 2 Dk ; ð2:25Þ
 Double Hopf bifurcation in delay differential equations 289

where
Z  
0
Ek1
ðq1 ;q2 ;q3 ;q4 ;lÞ ¼ Rk1
ðq1 ;q2 ;q3 ;q4 ;lÞ þL e iðsÞWq
Nk1
ðq1 ;q2 ;q3 ;q4 ;lÞ ðsÞds  Nk1
ðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ:
0

ð2:26Þ
We will also use the fact that all center manifolds have rank in the subspace Xs , which
implies in particular that hw1 ; hkðq1 ;q2 ;q3 ;q4 ;lÞ ðÞi ¼ 0 for all ðq1 ; q2 ; q3 ; q4 ; lÞ 2 Dk and k > 1.
Consequently, we obtain from Eq. (2.24) that
hw1 ; eiWq  ihkðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ ¼ vk1
ðq1 ;q2 ;q3 ;q4 ;lÞ ; for ðq1 ; q2 ; q3 ; q4 ; lÞ 2 Dk ; ð2:27Þ
where the vector vk1 ðq1 ;q2 ;q3 ;q4 ;lÞ is given by
 Z 
k1
vðq1 ;q2 ;q3 ;q4 ;lÞ ¼  w1 ; eiðsÞWq Nk1
ðq1 ;q2 ;q3 ;q4 ;lÞ ðsÞds : ð2:28Þ
0

We introduce the ðn þ 1Þ  ðn þ 1Þ matrices M1 and M2 defined by

!
Dðix1 Þ w> 1 ð0Þ
M1 ¼ ;
hw1 ; eix1  i 0
and
!
Dðix2 Þ w>
3 ð0Þ
M2 ¼ :
hw3 ; eix2  i 0
Then it follows from the Eqs. (2.25) and (2.27) that the vectors defined by
!
^k hkðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ
hðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ ¼
0
satisfy
!
Ek1
ðq1 ;q2 ;q3 ;q4 ;lÞ
M1 h^kðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ ¼ for Wq ¼ x1 ;
vk1
ðq1 ;q2 ;q3 ;q4 ;lÞ

!
Ek1
ðq1 ;q2 ;q3 ;q4 ;lÞ
M1 h^kðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ ¼ for Wq ¼ x1 ;
vk1
ðq1 ;q2 ;q3 ;q4 ;lÞ

!
Ek1
ðq1 ;q2 ;q3 ;q4 ;lÞ
M2 h^kðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ ¼ for Wq ¼ x2
vk1
ðq1 ;q2 ;q3 ;q4 ;lÞ

and
!
Ek1
ðq1 ;q2 ;q3 ;q4 ;lÞ
M2 h^kðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ ¼ for Wq ¼ x2 :
vk1
ðq1 ;q2 ;q3 ;q4 ;lÞ

Now to complete the proof of the theorem, we need to show the uniqueness of
hkðq1 ;q2 ;q3 ;q4 ;lÞ ð0Þ, which is a consequence of the next lemma: h
290 R. Qesmi, M.A. Babram

Lemma 1. The ðn þ 1Þ  ðn þ 1Þ matrices M1 and M2 are invertible.

Proof. Let x 2 Rn ; g 2 R such that

 
x
M1 ¼ 0:
g
Then
(
Dðix1 Þx þ gw>
1 ð0Þ ¼ 0;
hw1 ; eix1  ix ¼ 0:
From the fact that w1 ð0ÞDðix1 Þ ¼ 0, we have g ¼ 0, and the above system becomes

Dðix1 Þx ¼ 0
hw1 ; eix1 : ix ¼ 0:
From the fact that dim ker Dðix1 Þ ¼ 1 and since we have Dðix1 Þ/1 ð0Þ ¼ 0, the first
equation implies that x 2 spanf/1 ð0Þg, which follows that there exists c 2 R such that
x ¼ c/1 ð0Þ. From the second equation, we have chw1 ; /1 i ¼ 0 which yields c ¼ 0 and
finally we obtain x ¼ 0. Consequently M1 is invertible. By the same manner as before,
we prove that M2 is invertible. h

2.2. Normal forms calculation

For j 2 N and Y a normed space, Vm j ðYÞ denotes the space of homogeneous polynomi-
als of degree j and m variables z ¼ ðz1 ; z2 ; . . . ; zm Þ with coefficients in Y. In other words
( )
X
Vmj ðYÞ ¼ cq zq : q 2 N m ; cq 2 Y
jqj¼j
Pm
where jqj ¼ i¼1 qi .
Let
e aÞ :¼ Wð0Þ½ðLðaÞ  L0 ÞðUz þ hðz; aÞÞ þ FðUz þ hðz; aÞ; aÞ:
Hðz;
If we suppose that the homogeneous polynomials ðhi ÞiP2 are already computed, the
P
e¼
ODE (2.6) is known at any order. In other words, if H e e
iP2 H i , then H i is known
for all i P 2 once ðhi ÞiP2 are computed.
Consider a nonlinear transformation of the form
z ¼ x þ Vðx; aÞ; ð2:29Þ
m m
where V is mapping from R into R such that Vð0Þ ¼ 0 and DVð0Þ ¼ 0. Assume that
the effects of the above change of variables on the reduced system (2.6) is of the form
x_ ¼ Bx þ Nðx; aÞ: ð2:30Þ
Let
X X
V¼ Vi ; N ¼ Ni
iP2 iP2
 Double Hopf bifurcation in delay differential equations 291

and
X
ðLðaÞ  L0 Þð/Þ þ Fð/; aÞ ¼ fj ð/; aÞ: ð2:31Þ
iP2

Definition 2.5. Let j P 2. We denote by M1j the operator defined by

M1j p ðxÞ ¼ DpðxÞBx  BpðxÞ;

m
for all p 2 Vm
j ðR Þ.

The operators M1j , defined above, are exactly the Lie brackets that appear in com-
puting normal forms for finite dimensional ODEs. It is well known that the space
m
Vm
j ðR Þ may be decomposed as
c
m 1 1
Vm
j ðR Þ ¼ ImðMj Þ  Im Mj ; ð2:32Þ
c
where the complementary space Im M1j is not uniquely determined.
The following theorem summarizes the result for the recursive efficient approach for
computing the kth order normal form and the kth order nonlinear transformation Vk
in terms of center manifold coefficients ðhi Þ26i6k as well as the terms ðNi Þ26i6k1 already
computed.

Theorem 2.6. Fix k P 2. Then, the recursive formula for computing the coefficients of
the normal form associated to the decomposition (2.32) is given by

e k ðx; aÞ  M1 Vk ðx; aÞ
Nk ðx; aÞ ¼ H k

Xk1 n o
þ DHe kþ1i ðx; aÞVi ðx; aÞ  DVi ðx; aÞNkþ1i ðx; aÞ
i¼2

½k2 X ( )
X 1
kj X
j e
þ D H i ðx; aÞ Vl1 Vl2    Vlj ðx; aÞ; ð2:33Þ
j! i¼j
j¼2 l 1 þl2 þþlj ¼kiþj

and for all i P 2,

X
i1
e i ðx; aÞ ¼ Wð0Þfi ðUx; aÞ þ
H Wð0ÞDfj ðUx; aÞhiþ1j ðxÞ
j¼2

½ X
i ( )
X 2
1
ij X
j
þ Wð0Þ D fp ðUx; aÞ hl1 hl2    hlj ðx; aÞ:
j¼2
j! p¼j l1 þl2 þþlj ¼kiþj

e aÞ near h ¼ 0 is given by
Proof. The Taylor expansion of Hðx;
X1
e aÞ ¼ wð0ÞfðUx; aÞ þ wð0ÞDfðUx; aÞhðx; aÞ þ wð0Þ
Hðx;
i
Di fðUx; aÞ½hðx; aÞ :
iP2
i!

So, according to the representations (2.11) and (2.31), the above equation becomes
292 R. Qesmi, M.A. Babram

X X XX
e i ¼ wð0Þ
H fi ðUx; aÞ þ Wð0Þ Dfi ðUx; aÞhj ðx; aÞ
iP2 iP2 jP2 iP2
" #i
X1X X
j
þ Wð0Þ D fi ðUx; aÞ hl ðx; aÞ : ð2:34Þ
jP2
j! iP2 lP2

Comparing the same order terms of (2.34) yields the second formula of the theorem.
On the other hand, substituting the change of variables of (2.29) into (2.6), we have
" #
X X
Bz þ e i ðz; aÞ ¼ ½I þ DVðx; aÞ Bx þ
H Nk ðx; aÞ ;
iP2 kP2

which can be rearranged as

X X X X
Nk ðx; aÞ ¼ e k ðx þ Vðx; aÞ; aÞ þ
H BVk ðx; aÞ  DVk ðx; aÞBx
kP2 kP2 kP2 kP2
" #" #
X X
 DVk ðx; aÞ Nk ðx; aÞ :
kP2 kP2

e þ Vðx; aÞ; aÞ near V ¼ 0 to rewrite the above

Then we use Taylor expansion of Hðx
equation as
X X XX X
Ni ðx;aÞ ¼ e i ðx;aÞ þ
H e i ðx;aÞVj ðx;aÞ þ
DH M1j Vj ðx;aÞ
iP2 iP2 jP2 iP2 jP2
" #" # " #i
X X X1X X
 DVj ðx;aÞ Nj ðx;aÞ þ e i ðx;aÞ
Dj H Vl ðx;aÞ :
jP2 jP2 jP2
j! iP2 lP2

As for (2.34), we compare the same order terms, which gives the formula (2.33) in
theorem. h
The formula (2.33) given in the above theorem can be rewritten in a compact form
Nk ¼ Hk  M1k Vk ; ð2:35Þ
where
k1 n
X o
e k ðx; aÞ þ
Hk ðx; aÞ ¼ H e kþ1i ðx; aÞVi ðx; aÞ  DVi ðx; aÞNkþ1i ðx; aÞ
DH
i¼2

½ X
k ( )
X 2
1
kj X
j e
þ D H i ðx; aÞ Vl1 Vl2    Vlj ðx; aÞ:
j! i¼j
j¼2 l 1 þl2 þþlj ¼kiþj

Now consider the formula (2.33). On one hand, by the aid of (2.32), we can compute an
adequate nonlinear transformation Vk such that
M1k Vk ¼ Pk Hk ; ð2:36Þ
which leads to
Nk ¼ I  pK Hk : ð2:37Þ
 Double Hopf bifurcation in delay differential equations 293

In other words, we can find a change of variables that affects the reduced system by
taking away nonlinear terms (called non-resonant terms) that are in the range
subspaces ImðM1k Þ and conserving only terms (called resonant terms) that are in the
c
complementary subspace ImðM1k Þ .

Remark 1. This theorem is an important approach because it combines center manifolds

and normal form schemes into one step. In fact, the kth order Nk of the normal form
depends upon the known ðNi Þ26i6k1 and the results ðVi Þ26i6k as well as the terms
ðhi Þ26i6k2 of center manifolds. We will see in Section 3 that the above theorem gives a
useful and efficient tool to obtain the normal forms for FDEs in a compact form.

3. APPLICATION

The mathematical model which describes the regenerative cutting tool chatter in
turning was introduced in [12,17] and given by
y00 ðtÞ þ AðaÞy0 ðtÞ þ BðaÞyðtÞ þ CðaÞyðt  r0  qðtÞÞ ¼ fða; yðtÞ; yðt  r0  qðtÞÞ;
ð3:1Þ
where yðtÞ is roughly the displacement normal to the machined surface at time t; a is a
bifurcation parameter,  is a small perturbation parameter and f is a nonlinear function
given by
fða; yðtÞ; yðt  r0  qðtÞÞ
¼ c20 ðaÞy2 ðtÞ þ c11 ðaÞyðtÞyðt  r0  qðtÞÞ þ c02 ðaÞy2 ðt  r0  qðtÞÞ
þ c30 ðaÞy3 ðtÞ þ c21 ðaÞy2 ðtÞyðt  r0  qðtÞÞ
þ c12 ðaÞyðtÞy2 ðt  r0  qðtÞÞ þ c03 ðaÞy3 ðt  r0  qðtÞÞ:
The regenerative effect enters the equation of motion through chip thickness
1
j ¼ yðtÞ  yðt  sÞ; s ¼ r0 þ qðtÞ ¼ ; qðtÞ ¼ l cosðmtÞ
Xðt; Þ
and X is the spindle speed. Here m is the frequency of the periodic fluctuations and l is
the amplitude. The width-of-cut j, and time delay r0 , (which corresponds to the inverse
of the mean value of the spindle speed X) are natural control or bifurcation parameters
in the machine cutting process.
We consider the following hypothesis on the linear part of system (3.1):

(H’): In the absence of periodic perturbation, i.e.,  ¼ 0, FDE exhibits a Hopf

bifurcation at a ¼ 0, with a simple pair of pure imaginary eigenvalues on
the imaginary axis.

3.1. Augmented autonomous system

The explicit time-dependent delay terms are replaced by state-dependent delay terms.
To this end, we define the additional coordinations:
294 R. Qesmi, M.A. Babram

x1 ðtÞ ¼ yðtÞ; x2 ðtÞ ¼ x_1 ðtÞ; x3 ðtÞ ¼ l cosðmtÞ; x4 ðtÞ ¼ l sinðmtÞ:

Then the relationship between x3 ðtÞ and x4 ðtÞ can be given as follows:
x_ 3 ðtÞ ¼ mx4 ðtÞ; x_ 4 ðtÞ ¼ mx3 ðtÞ;
where the parameter l will be identified later as an amplitude. By augmenting x3 ðtÞ and
x4 ðtÞ to the problem, we get the autonomous system
0 1 0 1
0 1 0 0 0 1 0 0
B C B C
B BðaÞ AðaÞ 0 0 C B CðaÞ 0 0 0 C
_ B
XðtÞ ¼ B C XðtÞ þ BB CXðt  r0  x4 ðtÞÞ
@ 0 0 0 m C
A @ 0 0 0 0C A
0 0 m 0 0 0 0 0
0 1
0
B C
BfC
þB
B0C
C ð3:2Þ
@ A
0

where
0 1
x1 ðtÞ
B x ðtÞ C
B 2 C
XðtÞ ¼ B C
@ x3 ðtÞ A
x4 ðtÞ

and
f ¼ c20 ðaÞx21 ðtÞ þ c11 ðaÞx1 ðtÞx1 ðt  r0  x4 ðtÞÞ þ c02 ðaÞx21 ðt  r0  x4 ðtÞÞ
þ c30 ðaÞx31 ðtÞ þ c21 ðaÞx21 ðtÞx1 ðt  r0  x4 ðtÞÞ þ c12 ðaÞx1 ðtÞx21 ðt  r0  x4 ðtÞÞ
þ c03 ðaÞx31 ðt  r0  x4 ðtÞÞ:

Here, it is assumed that jj is small enough. Using Taylor expansion of the term
Xðt  r0  x3 ðtÞÞ and substituting the obtained formula into Eq. (3.2), we obtain the
following system with constant delay
0 1 0 1
0 1 0 0 0 1 0 0
B BðaÞ AðaÞ 0 0 C B CðaÞ 0 0 0 C
_ ¼B
XðtÞ B
C
CXðtÞ þ B
B C
CXðt  r0 Þ
@ 0 0 0 m A @ 0 0 0 0A
0 0 m 0 0 0 0 0
0 1
0
B fðx ; aÞ C
B t C
þB C
@ 0 A
0
ð3:3Þ
 Double Hopf bifurcation in delay differential equations 295

where
X
fðxt ; aÞ ¼ ck1 k2 ðaÞxk11 ðtÞxk12 ðt  r0 Þ þ CðaÞx2 ðt  r0 Þx4 ðtÞ
jkj¼2
X
þ ck1 k2 ðaÞxk11 ðtÞxk12 ðt  r0 Þ
jkj¼3
X
 k2 ck1 k2 ðaÞxk11 ðtÞx1k2 1 ðt  r0 Þx2 ðt  r0 Þx4 ðtÞ
jkj¼2

CðaÞ 4
þ ½BðaÞx1 ðt  r0 Þ þ AðaÞx2 ðt  r0 Þ þ Oðjxj Þ
2
As a FDE in C :¼ Cð½r0 ; 0; R4 Þ, Eq. (3.3) becomes Eq. (2.1), where
0 1 0 1
0 1 0 0 0 1 0 0
B BðaÞ AðaÞ 0 0 C B CðaÞ 0 0 0 C
B C B C
LðaÞ/ ¼ B C/ð0Þ þ B C/ðr0 Þ;
@ 0 0 0 m A @ 0 0 0 0A
0 0 m 0 0 0 0 0
and
0 1
0
B fð/; aÞ C
B C
Fð/; aÞ ¼ B C:
@ 0 A
0
Define the following constants
K1 :¼ ðC1 eir0  iA1  B1 Þ=N; N ¼ A0 þ 2i  r0 C0 eir0 ;
 
C0 ir0 2C0 ð1 þ mÞeið1þmÞr0 2C0 ð1  mÞeið1mÞr0
K2 ðmÞ :¼ e 1 þ
2N TðmÞ TðmÞ

TðmÞ ¼ B0  ð1 þ mÞ2 þ C0 eið1þmÞr0 þ ið1 þ mÞA0

 
c11 ð2c02 þ c11 Þ 2c11 c02 þ c11 c20 þ 2c02 c20 e2ir0
K3 ¼ c12 þ þ
B0 þ C0 Tð1Þ N
  ir0
c11 ð2c20 þ c11 Þ c11 c20 e
þ c12 þ þ
B0 þ C0 Tð1Þ N
 2 2
4c þ c11 þ 4c02 c20 þ 4c11 c20 þ 2c11 c02
þ 3c03 þ 2c21 þ 02
B0 þ C0

2c02 c20 þ 2c11 c20 þ c11 c02 ð3c30 þ 2c12 Þeir0
þ
Tð1Þ N2
4c2 þ c211 þ 4c02 c20 þ 2c11 c20 þ 4c11 c02
þ 20
NðB0 þ C0 Þ
2c þ c11 þ 2c11 c02 e4ir0 þ c211 e3ir0 þ 2c202 e3ir0
2 2
þ 20 :
NTð1Þ
296 R. Qesmi, M.A. Babram

Using the main result of Section 2.1, we obtain the following theorem:

Theorem 3.1. Consider the class of equations with periodic delay (3.1). Then the normal
form up to third order, in polar coordinates, associated to this class of equations is given
by
2
_ ¼ ½KR
rðtÞ R R 2
1 a þ K2 ðmÞjlj þ K3 r ðtÞrðtÞ
_
/ðtÞ
2
¼ 1 þ aKI þ KI ðmÞjlj þ KI r2 ðtÞ;
1 2 3

where KR
i and KIi are respectively, the real and imaginary parts of Ki ; i ¼ 1; 2; 3.

Proof. The associated characteristic equation is given by

(
k2 þ AðaÞk þ BðaÞ þ CðaÞekr0 ¼ 0
k2 þ m2 ¼ 0:
It is obvious that the eigenvalues of the augmented oscillator is k ¼ im. The transcen-
dental equation, at a ¼ 0, has a pair of pure imaginary roots, which we normalize to
one, that is k ¼ i. This normalization of chatter frequency gives the following relation
between the mean delay and the linear critical coefficients
C0 cosðr0 Þ ¼ 1  B0 and C0 sinðr0 Þ ¼ A0 ;

where A0 :¼ Að0Þ; B0 :¼ Bð0Þ and C0 :¼ Cð0Þ. Thus, at a ¼ 0, Eq. (3.3) has two pairs of
eigenvalues on the imaginary axis and all other eigenvalues have negative real parts.
Let
K ¼ fþi; img
and Xc the generalized eigenspace associated with the eigenvalues of K with the basis
U ¼ ½/1 ; /2 ; /3 ; /4  given by
2 ih 3
e eih 0 0
6 ieih ieih 0 0 7
6 7
UðhÞ ¼ 6 7; for h 2 ½r0 ; 0:
4 0 0 eimh eimh 5
0 0 ieimh ieimh

Let W be the basis for the generalized eigenspace (dual space Xc ) of the transposed
equation associated with K such that hW; Ui ¼ I. We have
2 ðA iÞ 3
2 3 0 is
eis eN 0 0
w1 ðsÞ N
6 7
6 w ðsÞ 7 6 ðA0 þiÞ eis eis 0 0 7
6 2 7 6 N N 7
WðsÞ ¼ 6 7¼6 7; for s 2 ½0; r0 :
4 w3 ðsÞ 5 6 0 0 eims eims 7
i 2 5
4 2
w4 ðsÞ 0
ims
0 e 2 i e 2
ims

The long term behavior of solutions of the original FDE is described by the solutions
of the four-dimensional reduced system on a local center manifold
 Double Hopf bifurcation in delay differential equations 297

y0 ðtÞ ¼ By þ R1 ðy; aÞ þ R2 ðy; aÞ þ h:o:t:; ð3:4Þ

1 2
where R ðy; aÞ and R ðy; aÞ are respectively the homogeneous parts of degree two and
three in ðy; aÞ given in Section 2. The coefficients R1ðq1 ;q2 ;q3 ;q4 ;lÞ ; ðq1 ; q2 ; q3 ; q4 ; lÞ 2 D2 of
R1 ðy; aÞ, are given by
R1ð2;0;0;0;0Þ ¼ c20 þ c11 eir0 þ c02 e2ir0 ; R1ð0;2;0;0;0Þ ¼ R1ð2;0;0;0;0Þ ;
R1ð1;1;0;0;0Þ ¼ 2c20 þ 2c11 cosðr0 Þ þ 2c02 ; R1ð0;2;0;0;0Þ ¼ R1ð2;0;0;0;0Þ ;
R1ð1;0;1;0;0Þ ¼ C0 eir0 ; Rð0;1;0;1;0Þ ¼ R1ð1;0;1;0;0Þ ;
R1ð1;0;0;1;0Þ ¼ C0 eir0 ; R1ð0;1;0;0;1Þ ¼ R1ð1;0;0;0;1Þ ;
R1ð1;0;0;0;1Þ ¼ C1 eir0  iA1  B1 ; R1ðq;lÞ ¼ 0; otherwise:
Then, by applying directly the schemes for computing the terms of center manifolds
obtained in Section 2, we obtain all the terms of center manifolds up to second order
as:
2 3
eWq h eih ih
Dðq;lÞ
 NðiþWq Þ
e
 NðiW qÞ
R1ðq;lÞ
6 7
6 Wq e W q h eih
7
6 þi þ eih
R1ðq;lÞ 7
h2ðq;lÞ ðhÞ ¼6 Dðq;lÞ NðiþWq Þ NðiWq Þ 7
6 7
4 0 5
0
 W r
2 
where Dq :¼ ð1 þ ðWq Þ Þ þ C0 e q 0  2i1 ði  Wq Þeir0  2i1 ði þ Wq Þeir0 .
Now, let us compute the normal forms of the reduced system. One can see that

R2 ðy; aÞ ¼ Wð0ÞDF2 ð/y; 0Þh2 ðy; 0Þ þ Wð0ÞF3 ðUy; 0Þ þ Oðajyj2 ; a2 jyjÞ:

2
However, the terms Oðajyj ; a2 jyjÞ are irrelevant to determine the generic double Hopf
bifurcation. Hence, for simplicity, we will omit these terms in the calculations of
normal forms of order 3.
Consider a nonlinear transformation of the form

y ¼ z þ V2 ðzÞ
where V2 ðzÞ is a homogeneous polynomial of degree two. The effects of the above
change of variables on the reduced system is
z0 ¼ Bz þ Nðz; aÞ;
P
with Nðz; aÞ ¼ kP2 Nk ðz; aÞ.
By using the formulas given by Theorem 2.6, we obtain

e j  M1 Vj ;
Nj ¼ H j ¼ 2; 3; ð3:5Þ
j

where
e 2 ðy; aÞ ¼ Wð0ÞðLðaÞ  Lð0ÞÞUy þ Wð0ÞF2 ðUy; 0Þ
H ð3:6Þ
298 R. Qesmi, M.A. Babram

and
e aÞ ¼ Wð0ÞF3 ð/y; 0Þ þ Wð0ÞDF2 ðUy; 0Þ½h2 ðz; aÞ þ UV2 ðz; aÞ:
Hðy;

Note that
80 2 1 0 10 10 10 10 10 10 19
> z1 z2 z1 z3 z4 0 0 0 0 0 0 >
>
>B >
< 0 C B 0 C B z z2C B
z z z C B 0 C B 0 C B 0 C B 0 C>
=
B C B C B 1 2 C B 2 3 4 C B C B C B C B C
kerðM13 Þ ¼ span B CB
; CB
; CB
; CB 2 CB
; ; CB
; CB
; C :
>
>@ 0 A @ 0 A @ 0 A @ 0 A @ z3 z4 A @ z1 z2 z3 A @ 0 A @ 0 A> >
>
: >
;
0 0 0 0 0 0 z3 z24 z1 z2 z4
ð3:7Þ

After simple calculations, one can obtain

21 X 3
N
R1ðq;lÞ zq al
6 ðq;lÞ2D2 7
6 X 7
61 1 q l7
6 R z a 7
H~2 ðz; aÞ ¼ 6 N ðq;lÞ 7:
6 ðq;lÞ2D2 7
6 7
4 0 5
0

So, by virtue of (3.5) and (3.6) it follows that

2 X R2 3
1 ðq;lÞ q l
z a
6 N ðq;lÞ2D ðWq iÞ 7
6 2 7
6 X 1 7
61 Rðq;lÞ
q l7
V2 ðz; aÞ ¼ 6
6N ðWq iÞ
z a 7
7
6 ðq;lÞ2D2 7
6 7
4 0 5
0

and
2 3
R1ð1;0;0;0;1Þ z1 a
6 1 7
6R 7
N2 ðz; aÞ ¼ 6 ð0;1;0;0;1Þ z2 a 7:
6 7
4 0 5
0
On the other hand, we have
21 X 3
N
bðq;lÞ zq al
6 ðq;lÞ2D3 7
6 X 7
61 7
6N bðq;lÞ zq al 7
wð0ÞF3 ð/z; 0Þ ¼ 6 7;
6 ðq;lÞ2D3 7
6 7
4 0 5
0
 Double Hopf bifurcation in delay differential equations 299

where
bð3;0;0;0;0Þ ¼ c30 þ c21 eir0 þ c12 e2ir0 þ c03 e3ir0 ;
bð2;1;0;0;0Þ ¼ 3c30 þ c21 ð2eir0 þ eir0 Þ þ c12 ð2 þ e2ir0 Þ þ 3c03 e3ir0 ;
bð2;0;1;0;0Þ ¼ c11 eir0  2c02 e2ir0 ; bð2;0;0;1;0Þ ¼ c11 eir0 þ 2c02 e2ir0 ;

bð1;2;0;0;0Þ ¼ bð2;1;0;0;0Þ ; bð1;1;1;0;0Þ ¼ bð1;1;0;1;0Þ ; bð0;3;0;0;0Þ ¼ bð3;0;0;0;0Þ ;

bð0;2;1;0;0Þ ¼ bð2;0;0;1;0Þ ; bð0;2;0;1;0Þ ¼ bð2;0;1;0;0Þ ; bð1;1;0;1;0Þ ¼ c11 ðeir0  eir0 Þ;
bðq;lÞ ¼ 0; otherwise;

and
21 X 3
N
cðq;lÞ zq al
6 ðq;lÞ2D3 7
6 7
6 X 7
61 c z q l7
a
6N ðq;lÞ 7
Wð0ÞDF2 ðUy; 0Þ½h2 ðz; 0Þ þ UV2 ðz; 0Þ ¼ 6 ðq;lÞ2D 7
6 3 7
6 7
6 0 7
4 5
0

where cðq;lÞ are complex constants. According to (3.7), the specific coefficients cðq;lÞ that
are needed in the third order normal forms are cð2;1;0;0;0Þ ; cð1;2;0;0;0Þ ; cð1;0;1;1;0Þ and cð0;1;1;1;0Þ .
We have
!
iR1ð2;1;0;0;0Þ R1ð3;0;0;0;0Þ R1ð2;1;0;0;0Þ i2R1ð1;2;0;0;0Þ h 2 i
cð2;1;0;0;0Þ ¼   2
þ hð2;0;0;0;0Þ ð0Þ
N N N jNj 1
h i h i
 ð2c20 þ c11 eir0 Þ þ h2ð2;0;0;0;0Þ ðr0 Þ ðc11 þ 2c02 eir0 Þ þ h2ð1;1;0;0;0Þ ð0Þ
1 1
h i
ir0 2 ir0
 ð2c20 þ c11 e Þ þ hð1;1;0;0;0Þ ðr0 Þ ðc11 þ 2c02 e Þ;
1

iR1ð1;1;1;0;0Þ R1ð2;0;0;1;0Þ iR1ð1;1;0;1;0Þ R1ð2;0;1;0;0Þ h i

cð1;0;1;1;0Þ ¼ 2
 2
þ h2ð0;0;1;1;0Þ ð0Þ ð2c20 þ c11 eir0 Þ
jNj ð2  mÞ jNj ð2  mÞ 1
h i h i h i
þ hð0;0;1;1;0Þ ð0Þ ð2c11 þ c02 eir0 Þ þ iC0 hð1;0;0;1;0Þ ðr0 Þ  h2ð1;0;1;0;0Þ ðr0 Þ :
2 2
1 2 2

and cð1;2;0;0;0Þ ¼ cð2;1;0;0;0Þ ; cð0;1;1;1;0Þ ¼ cð1;0;1;1;0Þ .

This leads to
2 3
ðbð2;1;0;0;0Þ þ cð2;1;0;0;0Þ Þz21 z2 þ bð1;0;1;1;0Þ þ cð1;0;1;1;0Þ Þz1 z3 z4
6 7
6 ðbð2;1;0;0;0Þ þ cð2;1;0;0;0Þ Þz1 z22 þ bð1;0;1;1;0Þ þ cð1;0;1;1;0Þ Þz2 z3 z4 7
N3 ðz; aÞ ¼ 6
6
7:
7
4 0 5
0
300 R. Qesmi, M.A. Babram

Finally, the truncated normal forms of the reduced Eq. (3.4) up to third order reads
8
>
> z01 ðtÞ ¼ ði þ R1ð1;0;0;0;1Þ aÞz1 ðtÞ þ ðbð2;1;0;0;0Þ þ cð2;1;0;0;0Þ Þz21 z2
>
>
>
> þðbð1;0;1;1;0Þ þ cð1;0;1;1;0Þ Þz1 z3 z4
>
>
>
< 0
z2 ðtÞ ¼ ði þ R1ð0;1;0;0;1Þ aÞz2 ðtÞ þ ðbð1;2;0;0;0Þ þ cð1;2;0;0;0Þ Þz22 z1
>
> þðbð0;1;1;1;0Þ þ cð0;1;1;1;0Þ Þz2 z3 z4
>
>
>
> z0 ðtÞ ¼ imz3 ðtÞ
>
> 3
>
: 0
z4 ðtÞ ¼ imz4 ðtÞ

which, by transforming the above equation in polar coordinates, completes the proof
of theorem. h

4. CONCLUSION

Methodology for simultaneous computation of center manifolds and normal forms for
double Hopf bifurcation, in delay differential equations, has been developed. The cal-
culations and formulas are given in an explicit iterative procedure, and thus, are easy to
be implemented on a symbolic computation system. It has been shown by an example
that the method is computationally efficient, particularly, suitable for the computations
of complicated systems and higher-order center manifolds and normal forms.

ACKNOWLEDGMENT

The authors are very grateful to the anonymous referee for his/her valuable comments
and suggestions.

REFERENCES

[1] T. Faria, L.T. Magalhaès, Normal forms for retarded functional differential equations with parameters
and applications to Hopf bifurcation, J. Differential Eqs. 122 (2) (1995) 181–200.
[2] T. Faria, L.T. Magalhàes, Normal forms for retarded functional differential equations and applications to
Bogdanov–Takens singularity, J. Differential Eqs. 122 (2) (1995) 201–224.
[3] J.K. Hale, S.M. Verduyn Lunel, Introduction to Functional Differential Equations, Springer-Verlag, New
York, 1993.
[4] T. Inamura, T. Sata, Stability analysis of cutting under varying spindle speed, J. Faculty Eng. Univ.
Tokyo-B 23 (1) (1975) 13–29.
[5] Yu.A. Kuznetsov, Elements of Applied Bifurcation Theory, Springer-Verlag, New York, 1994.
[6] R. Qesmi, M. Ait Babram, M.L. Hbid, A Maple program for computing a terms of center manifolds, and
element of bifurcations for a class of retarded functional differential equations with Hopf singularity,
Appl. Math. Comput. 175 (2005) 932–968.
 Double Hopf bifurcation in delay differential equations 301

[7] R. Qesmi, M. Ait Babram, M.L. Hbid, Center manifolds and normal forms for a class of retarded
functional differential equations with parameter associated with Fold-Hopf singularity, Appl. Math.
Comput. 181 (2006) 220–246.
[8] S.A. Tobias, Machine-Tool Vibration, Blackie, 1965.
[9] S. Wiggins, Introduction to applied nonlinear dynamical systems and chaosTexts in Applied Mathemat-
ics, vol. 2, Springer-Verlag, New York, 1990.