Symbolic computation for center manifolds and normal forms of Bogdanov bifurcation in retarded functional differential equations

Nonlinear Analysis 66 (2007) 2833–2851 www.elsevier.com/locate/na

Symbolic computation for center manifolds and normal forms of Bogdanov bifurcation in retarded functional differential equations R. Qesmi a,∗ , M. Ait Babram b , M.L. Hbid a a D´epartement de Math´ematiques, Facult´e des Sciences Semlalia, Universit´e Cadi Ayyad, B.P. S15, Marrakech, Morocco b D´epartement de Math´ematiques, Facult´e des Sciences et t´echniques Gu´eliz, Universit´e Cadi Ayyad,

Marrakech, Morocco Received 28 February 2006; accepted 13 April 2006

Abstract In this paper, we present explicit formulas for computing the coefficients of a center manifold for the Bogdanov singularity in autonomous retarded functional differential equations with parameters. As a consequence, normal forms associated with the flow on a center manifold up to an arbitrary order are derived. The explicit formulas have been implemented using the computer algebra system Maple. The paper ends with some examples given in order to show the applicability of the methodology and the convenience of the computer software. c 2006 Elsevier Ltd. All rights reserved.  Keywords: Retarded differential equation; Delay; Center manifold; Bogdanov

1. Introduction Center manifolds theory is of fundamental importance in the study of nonlinear dynamical systems when analyzing bifurcations of a given type. In fact, this theory allows us to reduce the study of a differential equation with delay near a non-hyperbolic equilibrium point to that of an ordinary differential equation on a finite-dimensional invariant manifold. This approach is particularly interesting when the starting point is an infinite-dimensional problem, such as a functional differential, partial differential or integro-differential equation, in fact, the reduction ∗ Corresponding author. Tel.: +212 6 19 37 259.

E-mail address: [email protected] (R. Qesmi). c 2006 Elsevier Ltd. All rights reserved. 0362-546X/$ - see front matter  doi:10.1016/j.na.2006.04.010

2834

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

also forms a qualitative simplification. Although center manifolds exist, they need not be unique, but 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 regularity; such an expansion may give an idea of the local dynamics near the steady state. Center manifolds for FDE have been considered in various papers (see Henry [13], Vanderbauwhede and Ioos [16], Diekmann and Van Gils [8] and the references therein). In [1,2] Ait Babram et al. have considered retarded differential functional equations without parameters (RFDE), and they derive an initial value problem (IVP) in finite dimensions to compute the terms of center manifolds. However, some of the ingredients of the IVP obtained do not have an explicit form, and the techniques used to compute the solution of the IVP are theoretical, and cannot be programmed on a machine. In a recent paper [15], Qesmi et al. have developed a method of computation of center manifolds for RFDE associated with Hopf singularity without parameters. In this paper we consider RFDE with several delays, depending on parameters, which have as the origin an equilibrium and such that for special (critical) values of the parameters, the linearization has double eigenvalues at zero and no spectrum in the right half-plane. If the geometric multiplicity of this eigenvalue is 1, then this is called a Bogdanov singularity. We give an algorithmic scheme for computing center manifolds for this class of FDE. The aim is to prove that the terms of a center manifold associated with a Bogdanov singularity up to the order of regularity for the RFDE satisfy efficient and explicit recursive formulas, and to give a method for obtaining these terms explicitly by using the bordering technique (see [14] for more details) which, at the same time, proves the uniqueness of the solutions of the recursive formulas and gives an efficient algorithm which can be easily implemented using the computer algebra system Maple. Having computed the equation restricted to a center manifold, the normal form theory, which plays an important role in the study of differential equations related to complex behavior such as bifurcation, comes to mind. It provides a convenient tool for computing a simple form of the original differential equations without changing the fundamental dynamical behavior of the system in the vicinity of the equilibrium. The normal form theory is usually applied together with the center manifold theory: Given a nonlinear system, the center manifold theory is applied before using the normal form theory. Many results related to this area may be found in, for example, Refs. [5–7]. Normal forms have been considered for RFDE by Faria and Magalha´es [9, 10]. The authors gave a method for obtaining normal forms of RFDE by considering an RFDE as an ODE in an appropriate infinite-dimensional phase space. However, the proposed method is difficult to apply in order to compute normal forms of higher orders because it involves the necessity of solving functional equations, which is not easy in most cases. For isolating the difficulties, the evaluation of normal forms after computing a center manifold of a given system is a most welcome process. This motivates the introduction of normal forms as a second stage of our work. It also allows us to show the applicability of the method as regards the computation of center manifolds. The paper is organized as follows. The preliminary Section 2 provides, among other things, elementary notions and a background for the theory of retarded functional differential equations. A theorem for the computational scheme of a center manifold and its proof are given in Section 3, and results of computing normal forms are given in the same section. Section 4 outlines the symbolic computation procedure using Maple, and examples are given in Section 5 to demonstrate the applicability of the results. Conclusions are drawn in Section 6. The results of the examples including the center manifolds and normal forms obtained by executing the Maple programs are listed in the Appendices.

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

2835

2. Notation and background We follow the notation of [11], where the general facts of the FDE theory used here can be found. Let r ≥ 0 be a given real number and C be the space of continuous functions mapping [−r, 0] into Rn with the uniform norm. If x is a continuous function taking [σ −r, σ + α], α ≥ 0, into Rn , then, for each t ∈ [σ, σ +α], we let x t ∈ C be defined by x t (θ ) = x(t +θ ), −r ≤ θ ≤ 0. We will need to consider autonomous RFDE and separate linear from nonlinear terms in the form d x(t) = Lx t + g(x t ), (2.1) dt where L is a bounded linear operator from C to Rn and g a sufficiently smooth function mapping C into Rn such that g(0) = 0, Dg(0) = 0. The linearized system of Eq. (2.1) at x = 0 has the form d x(t) = Lx t . (2.2) dt 0 Such an operator L : C → Rn has the form L(φ) = −r dη(θ )φ(θ ), where η is an n × n matrix valued function of bounded variation defined on [−r, 0]. Define T (t)φ := x t where x is the solution of Eq. (2.2) with initial data φ ∈ C. Then (T (t))t ≥0 is an eventually compact semi-group on C. Let A be its infinitesimal generator, then its spectrum coincides with the point spectrum σ p (A), and λ ∈ σ (A) if and only if it satisfies the characteristic equation  0 n dη(θ )eλθ . (2.3) det Δ(λ) = 0 where Δ(λ) = λIR − −r

From the spectral theory we have that σ p (T (t)) \ {0} = et σ p ( A) for t ≥ r , then σ p (A) contains information about the behavior of the solutions of Eq. (2.2). Let Λ := σ (A) ∩ i R. The generalized eigenspace associated with Λ is denoted by X c . If Λ is a non-empty finite set, and d is the number of solutions of (2.3) in Λ counting multiplicities, then the dimension of X c is dim X c = d. The adjoint theory is used to find a decomposition C = X c ⊕ X s where both X c and X s are invariant under T (t) and A. With C ∗ = C ([0, r ], Rn∗ ), where Rn∗ is the n-dimensional space of row vectors, we consider the adjoint bilinear form on C ∗ × C  0 θ

ψ, φ = ψ(0)φ(0) − ψ(ξ − θ )dη(θ )φ(ξ )dξ. −r

0

In the special case when the linear operator L has the form m  Lφ = L j φ(−r j ) j =0

for each φ ∈ C, where the L j ’s are square matrices of order n, r0 = 0 and r j > 0 for all 0 < j ≤ m, this adjoint bilinear form reads m  0  ψ(ξ + r j )L j φ(ξ )dξ.

ψ, φ = ψ(0)φ(0) − j =0 −r j

2836

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

Let Φ be a basis of the generalized subspace X c and Ψ denote the basis for the dual subspace X c∗ such that Φ, Ψ  = In , we have X s = {φ ∈ C :< Ψ , φ >= 0}. We denote by B the d × d constant matrix such that T (t)Φ = Φet B ; its spectrum coincides with Λ. Using the variation of constants formula as in [11], it is shown that the solution of (2.1) with initial value φ satisfies  t x t = T (t)φ + T (t − s)X 0 g(x s )ds, t ≥ 0 (2.4) 0

where X 0 is given by  I n, θ = 0 X 0 (θ ) = R 0, −r ≤ θ < 0. Also, from the decomposition of C by Λ as C = X c ⊕ X s , then Eq. (2.4) is equivalent to  t x tc = T (t)φ c + T (t − σ )X 0c g(x σ )dσ, (2.5) 0  t x ts = T (t)φ s + T (t − σ )X 0s g(x σ )dσ (2.6) 0

where superscripts c and s designate the projections of the corresponding functions onto the subspaces X c and X s , respectively. One can show that Ψ , X 0  is well defined and Ψ , X 0  = Ψ (0). Therefore, the quantities X 0c and X 0s are given by X 0c = ΦΨ (0),

X 0s = X 0 − X 0c .

We remark that in view of the fact that X 0 is not continuous, the variation of constants formula (2.4) is not well defined on the space C. However, some care must be exercised in order to achieve a correct use of the formula (2.4). Let   Cˆ = φ : [−r, 0] −→ Rn ; φ is continuous on [−r, 0) and φ(0− ) ∈ Rn exists . ˆ We can show that for each φ ∈ C, ˆ there exists a and denote by . the supremum norm of C. unique continuous solution x := x(φ) of (2.2) on [0, ∞). In this way, we extend the domain of ˆ t ≥ 0. the solution semi-group T (.) of (2.2) to Cˆ by T (t)φ = x t (φ), φ ∈ C, In the sequel, we recall the definition of a local center manifold associated with Eq. (2.1). Definition 2.1 ([11]). Given a C 1 map h from Rd into X s , the graph of h is said to be a local center manifold associated with Eq. (2.1) if and only if h(0) = 0, Dh(0) = 0, and there exists a neighborhood V of zero in Rd such that, for each ξ ∈ V , there exist δ = δ(ξ ) > 0 and a function x defined on ]−δ −r, δ[ such that x 0 = Φξ +h(ξ ) and x verifies Eq. (2.1) on ]−δ, δ[ and satisfies the identity x t = Φz(t) + h(z(t)),

for t ∈ [0, δ[,

where z(t) is the unique solution of the ordinary differential equation ⎧ ⎨d z(t) = Bz(t) + Ψ (0)g(Φz(t) + h(z(t))) dt ⎩z(0) = ξ, ξ ∈ Rd where B is the matrix given by AΦ = Φ B.

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

2837

Geometrically speaking, a center manifold is a graph of a given function mapping a neighborhood at zero in X c into X s which is tangent to X c and locally invariant under the semiflow generated by Eq. (2.1). The behavior of orbits of Eq. (2.1) in C near the singularity at the origin can be completely described by the restriction of the flow to the associated center manifold, which is necessarily finite dimensional. 3. Main results In this section we present our main results concerning the computation of the terms of center manifolds, and normal forms for the following class of RFDE m  d x(t) = L j (α)x(t − r j ) + f (x(t), x(t − r1 ), . . . , x(t − rm ), α), dt j =0

(3.1)

where α ∈ R p , α → L j (α) is a C ∞ function with values in the space of square matrices of order n. f : Rn×(m+1) × R p → Rn is assumed to be sufficiently smooth ( f ∈ C ∞ ) such that: f (0, . . . , 0, α) = 0 and D f (0, . . . , 0, α) = 0 for all α ∈ R p , r0 = 0 and r j > 0 for all 1 < j ≤ m. If we define r := max1≤i≤m {ri }, Cn := C([−r, 0], Rn ) the space of continuous functions from [−r, 0] to Rn , and for φ ∈ Cn L(α)φ := L 0 (α)φ(0) +

m 

L j (α)φ(−r j ),

F(φ, α) := f (φ(0), φ(−r1 ), . . . , φ(−rm ), α)

j =1

then Eq. (3.1) reads as the equation x(t) ˙ = L(α)x t + F(x t , α).

(3.2)

We define L 0 = L(0). Throughout this paper we assume that the following hypothesis is satisfied: (H): The linear equation x(t) ˙ = L 0 x t has λ = 0 as a double characteristic value and no other characteristic values with zero real part. The way in which we will consider center manifolds for the differential equation (3.2) is to include the parameter α as a new dependent variable as follows x(t) ˙ = L 0 x t + [L(α) − L 0 ]x t + F(x t , α) α(t) ˙ = 0. The solutions of this system are of the form (x(t), α(t)) ∈ Rn+ p , the phase space is C˜ := Cn+ p , and the system can be written as ˙˜ = L˜ x˜t + F( ˜ x˜t ), (3.3) x(t)    0 ˜ ˜ v) := L(v(0)) − L 0 u + F(u, v(0)), 0 with u ∈ where L(u, v) := L u, 0 and F(u, Cn , v ∈ C p . It is now possible to apply to (3.3) the center manifold theory introduced in Section 2. Let A and A˜ denote the infinitesimal generators associated with the equations x(t) ˙ = L 0 xt , ˜ t , respectively. The equation in C p , α(t) ˙ = 0, has a unique characteristic value and x(t) ˙ = Lx λ = 0, with multiplicity p, and the associated generalized eigenspace consists of the elements of

2838

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

˜ = C pwhich are constant functions, and it is denoted here also by R p . Let Λ := σ (A) ∩ i R, Λ Λ {0} and consider the decomposition Cn = X c ⊕ X s obtained as in section 2. Then we have the decomposition C˜ = X˜ c ⊕ X˜ s , where X˜ c = X c ×R p , X˜ s = X s × R and R = {v ∈ C p : v(0) = 0}. The bases of X˜ c and X˜ s , respectively, are given by the columns of the matrix Φ˜ and the rows of the matrix Ψ˜ ,     Φ 0 ˜ = Ψ 0 , Φ˜ = , Ψ 0 Ip 0 Ip     0 1 0 ˜ Ψ˜ = In+ p and Φ˙˜ = Φ˜ B˜ with B˜ = 0 0 0 . which satisfy Φ, 0

0

0p

Let h˜ be a graph of a local center manifold of Eq. (3.3). Then from Definition 2.1 there exists a neighborhood V of zero in R2+ p such that, for each ξ˜ ∈ V , there exists δ > 0 such that the ˜ ξ˜ ) exists on the interval ]−δ − r, δ[ and it is given by solution of (3.3) with initial data Φ˜ ξ˜ + h( ˜ z (t)), x˜t = Φ˜ z˜ (t) + h(˜

for t ∈ ]−δ, δ[,

such that z˜ (t) is a solution of the equation ⎧ ⎨ d z˜ (t) = B˜ z˜ (t) + Ψ ˜ (0) F( ˜ Φ˜ z˜ (t) + h(˜ ˜ z (t))) dt ⎩ z˜ (0) = ξ˜ , ξ˜ ∈ R2+ p .

(3.4)

˜ α) = (h(ξ, α), h 0 (ξ, α)) for ξ ∈ R2 and α ∈ R p , then Eq. Remark 3.1. (i) If we write h(ξ, (3.4) is equivalent to d dt



   z(t) Bz(t) = α 0

  0 + Ψ (0)[(L(α(0)+h 0(z(t), α)(0))− L )(Φz(t)+h(z(t), α))+ F(Φz(t)+h(z(t), α), α(0)+ h 0 (z(t), α)(0))] 0

  with B = 00 10 and z(0) = ξ ∈ R2 . Noting that h 0 (z(t), α)(0) = 0, because h 0 (z(t), α) ∈ R, and dropping the auxiliary equations introduced for handling the parameter, we get the reduced equation on the center manifold – represented as a graph over the ξ and α variables of h(ξ, α) for ξ and α sufficiently small – associated with the parameterized equation (3.1) as ⎧ ⎨d z(t) = Bz(t) + Ψ (0)[(L(α) − L 0 )(Φz(t) + h(z(t), α)) + F(Φz(t) + h(z(t), α), α)] dt ⎩z(0) = ξ, ξ ∈ R2 . (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 [17] for more details). 3.1. Computation of terms of a center manifolds for the RFDE In the following we give an analytic characterization of a center manifold associated with Eq. (3.1).

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

2839

Theorem 3.2. Given a C 1 map h from R2+ p into X s with h(0) = 0 and Dh(0) = 0, a necessary condition for the graph of h to be a local center manifold of Eq. (3.1) is that there exists a neighborhood V of zero in R2+ p such that, for each (ξ, α) ∈ V , ∂h(ξ, α) ∂ ∂h(ξ, α) (θ )Ψ (0) [(L(α) − L(0))(Φξ + h(ξ, α))] (θ ) + (h(ξ, α)) (θ ) = ξ2 ∂θ ∂ξ1 ∂ξ ∂h(ξ, α) + (θ )Ψ (0)F(Φξ + h(ξ, α), α) + Φ(θ )Ψ (0)[(L(α) ∂ξ − L(0))(Φξ + h(ξ, α))] + Φ(θ )Ψ (0)F(Φξ + h(ξ, α), α) (3.5) ∂ (h(ξ, α)) (0) = L 0 h(ξ, α) + (L(α) − L 0 )(Φξ + h(ξ, α)) + F(Φξ + h(ξ, α), α). ∂θ (3.6) Proof. Let h be a graph of a local center manifold of Eq. (3.1), then from Remark 3.1 there exists a neighborhood V of zero in R2+ p such that, for each (ξ, α) ∈ V , there exists δ > 0 such that the solution of (3.1) with initial data Φξ + h(ξ, α) exists on the interval ]−δ − r, δ[ and it is given by x t = Φz(t) + h(z(t), α),

for t ∈ ]−δ, δ[,

such that z(t) is a solution of the equation ⎧ ⎨d z(t) = Bz(t) + Ψ (0)[(L(α) − L 0 )(Φz(t) + h(z(t), α)) + F(Φz(t) + h(z(t), α), α)] dt ⎩z(0) = ξ, ξ ∈ R2 ,

  with B = 00 10 . The variation of constants formula of Eq. (3.1) can be written as  t    T (t − s)X 0 L(α) − L 0 x s + F(x s , α) ds, x t = T (t)φ +

t ≥ 0,

0

˙ = L 0 x t . It follows from where (T (t))t ≥0 is the semi-group solution of the linear equation x(t) the decomposition of the phase space by Λ : Cn = X c ⊕ X s that the function h satisfies  t T (t − σ )X 0s F(Φz(σ ) + h(z(σ ), α), α)dσ. h(z(t), α) = T (t)h(z(t), α) + 0

Then 1 1 [T (t)h(ξ, α) − h(ξ, α)] = [h(z(t), α) − h(z(0), α)] t t  1 t − T (t − σ )X 0s [(L(α) − L 0 )(Φz(σ ) + h(z(σ ), α)) t 0 + F(Φz(σ ) + h(z(σ ), α), α)]dσ, which implies, from the fact that h and F are smooth and T (.) is a strongly continuous semigroup on the Banach space Cn , that h(ξ ) is in the domain of A, and   ∂h(ξ, α) {Bξ + Ψ (0)[(L(α) − L 0 )(Φξ + h(ξ, α)) Ah(ξ, α) = ∂ξ + F(Φξ + h(ξ, α), α)]} = −X 0s [(L(α) − L 0 )(Φξ + h(ξ, α)) + F(Φξ + h(ξ, α), α)].

2840

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

Consequently, we have by evaluating the above equation at θ = 0 that   ∂h(ξ, α) ∂ (θ ) {Bξ + Ψ (0)[(L(α) − L 0 )(Φξ + h(ξ, α)) (h(ξ, α)) (θ ) = ∂θ ∂ξ + F(Φξ + h(ξ, α), α)]} + Φ(θ )Ψ (0)[(L(α) − L 0 )(Φξ + h(ξ, α)) + F(Φξ + h(ξ, α), α)], (3.7) which is the formula (3.5) of the theorem. On the other hand, it results from the fact that the semi-flow t −→ x t = Φz(t) + h(z(t), α) exists on the open ]−δ, δ[ that for θ ∈ ]−δ, 0] d d x 0 (θ ) (Φ(θ )ξ + h(ξ, α)(θ )) = dθ dθ d x(θ ) = dθ = L(α)(Φz(θ ) + h(z(θ ), α)) + F(Φz(θ ) + h(z(θ ), α), α). Consequently, by the fact that

d dθ Φ(0)ξ

= Φ(0)Bξ = L 0 (Φξ ), we obtain

∂ (h(ξ, α)) (0) = L 0 h(ξ ) + (L(α) − L 0 )(Φξ + h(ξ, α)) + F(Φξ + h(ξ, α), α), ∂θ which completes the proof of the theorem.

(3.8)



Let us recall that h has the same regularity as F. From this fact and in view of the assumed smoothness on F, for all m ∈ N, we can write h(ξ, α) =

m 

h k (ξ, α) + χ(ξ, α)

for (ξ, α) ∈ V,

k=2

where h k is the homogeneous part of degree k and χ(ξ, α) = o(|(ξ, α)|m ). Let k ∈ N, k ≥ 2. The homogeneous parts of degree k of Eqs. (3.5) and (3.6) are respectively given by ∂h k (ξ, α) ∂h k (ξ, α) = ξ2 + N k−1 (ξ, α), ∂θ ∂ξ1 ∂ (h k (ξ, α)) (0) = L 0 h k (ξ, α) + R k−1 (ξ, α), ∂θ

(3.9) (3.10)

where N k−1 (ξ, α) = ΦΨ (0)R k−1 (ξ, α) +

k−1  ∂h k− j +1 j =2

∂ξ

(ξ, α)Ψ (0)R j −1 ,

and R i−1 is the homogeneous part of degree i of R(ξ, α) = (L(α)− L 0 )(Φξ +h(ξ, α))+ F(Φξ + h(ξ, α), α). In particular, R 1 is the homogeneous part of degree 2 of (L(α) − L 0 )Φξ + F(Φξ, α) which is independent from terms of center manifolds.

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

2841 l

If ξ = (ξ1 , ξ2 ) ∈ R2 , q = (q1 , q2 ), α = (α1 , . . . , α p ), αl = α1l1 α2l2 · · · α pp for l = p (l1 , . . . , l p ) ∈ N p , and |l| = i=1 li , then we can write  q q h k (ξ, α) = h k(q,l) ξ1 1 ξ2 2 αl for some h k(q,l) ∈ X s , (q,l)∈Dk

N

k−1

(ξ, α)(θ ) =



(q,l)∈Dk

and R k−1 (ξ, α) =

 (q,l)∈Dk

q

q

k−1 1 2 l N(q,l) ξ1 ξ2 α

q

q

k−1 1 2 l R(q,l) ξ1 ξ2 α ,

k for some N(q,l) ∈ Xs ,

k−1 for some R(q,l) ∈ Rn ,

(3.11)

(3.12)

where Dk = {(q, l) ∈ N p : |(q, l)| = k}. Theorem 3.3. 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. (3.1) is given in a unique way by the following recursive formulas. For k = 2 and (q, l) ∈ D2 :   θ θ2 2 2 1 , (3.13) h (q1 ,0,l) (θ ) = h (q1 ,0,l) (0) + Φ(0)Ψ (0)R(q 1 ,0,l) 0 θ  θ (q1 + 1)h 2(q1 +1,q2 −1,l) (s)ds h 2(q1 ,q2 ,l) (θ ) = h 2(q1 ,q2 ,l) (0) +  0  θ θ2 1 Φ(0)Ψ (0)R(q for q2 = 0, (3.14) + 1 ,q2 ,l) 0 θ and for k > 2 and (q, l) ∈ Dk : ⎧  θ ⎪ k−1 k k ⎪ h (θ ) = h (q1 ,0,l) (0) + N(q (s)ds, ⎪ ⎪ 1 ,0,l) ⎨ (q1 ,0,l) 0  θ (3.15) k−1 k ⎪h k(q ,q ,l) (θ ) = h k(q ,q ,l) (0) + (q + 1)h (s) + N (s) ds, 1 ⎪ (q1 +1,q2 −1,l) (q1 ,q2 ,l) 1 2 1 2 ⎪ ⎪ 0 ⎩ for q2 = 0. with the vectors h k(q1 ,q2 ,l) (0), (q, l) ∈ Dk are given by solving the following system   k−1   E (q1 ,q2 ,l) h k(q1 ,q2 ,l) (0) = , M k−1 0 v(q 1 ,q2 ,l)

(3.16)

where M is the (n + 1) × (n + 1) matrix defined by   Δ(0) ψ2 (0) M=

ψ1 , IRn  0 and the second member of the system is a vector given by means of the coefficients of the center manifolds already computed. Proof. From relation (3.9), we have ∂h k(q1 ,0,l) (θ ) ∂θ

k−1 = N(q (θ ) 1 ,0,l)

(3.17)

2842

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

and ∂h k(q1 ,q2 ,l) (θ )

k−1 = (q1 + 1)h k(q1 +1,q2 −1,l) (θ ) + N(q (θ ), for q2 = 0, (3.18) 1 ,q2 ,l) ∂θ or equivalently ⎧  θ ⎪ k−1 k k ⎪ N(q (s)ds, ⎪ ⎪h (q1 ,0,l) (θ ) = h (q1 ,0,l) (0) + 1 ,0,l) ⎨ 0  θ (3.19) k−1 k k k ⎪ (q (θ ) = h (0) + + 1)h (s) + N (s) ds, h 1 ⎪ (q1 ,q2 ,l) (q1 ,q2 ,l) (q1 +1,q2 −1,l) (q1 ,q2 ,l) ⎪ ⎪ 0 ⎩ for q2 = 0.

The relation (3.10) is equivalent to ∂h k(q1 ,q2 ,l)

k−1 (0) = L 0 h k(q1 ,q2 ,l) + R(q,l) , ∂θ which gives by using the relations (3.17) and (3.18) that  0 k k−1 k−1 = N(q (0) L h (q1 ,q0,l) + R(q 1 0,,l) 1 ,0,l) k−1 k−1 L 0 h k(q1 ,q2 ,l) + R(q = (q1 + 1)h k(q1 +1,q2 −1,l) (0) + N(q (0) 1 ,q2 ,l) 1 ,q2 ,l)

for q2 = 0.

On the other hand, by taking the values of the linear operator L 0 on both sides of relation (3.19), and from the expression for the characteristic equation we obtain k−1 Δ(0)h k(q1 ,q2 ,l) (0) = E (q 1 ,q2 ,l)

where k−1 E (q 1 ,q2 ,l)

k−1 E (q 1 0,l)

=

k−1 R(q 1 q2 ,,l)

for (q1 , q2 , l) ∈ Dk ,

(3.20)

 .    k−1 k +L (q1 + 1)h (q1 +1,q2 −1,l) (s) + N(q1 ,q2 ,l) (s) ds 0

0

k−1 for q2 = 0, − (q1 + 1)h k(q1 +1,q2 −1,l) (0) − N(q ,q ,l) (0),  1 2  . k−1 k−1 k−1 = R(q + L0 N(q (s)ds − N(q (0). 1 0,l) 1 ,0,l) 1 ,0,l) 0

We will  also use thefact that all center manifolds have rank in the subspace X s , which implies that ψ1 , h k(q1 ,q2 ,l) (.) = 0 for all (q1 , q2 , l) ∈ Dk and k > 1. Consequently, we obtain from relation (3.19) that k−1

ψ1 , IRn  h k(q1 ,q2 ,l) (0) = v(q 1 ,q2 ,l)

for (q1 , q2 , l) ∈ Dk ,

(3.21)

k−1 where the vector v(q is given by 1 ,q2 ,l)   . k−1 k−1 v(q = − ψ1 , N(q (s)ds 1 ,0,l) 1 ,0,l) 0

and

  .  k−1 k−1 k v(q (q = − ψ , + 1)h (s) + N (s) ds 1 1 +1,q −1,l) (q (q ,q ,l) ,q ,l) 1 2 1 2 1 2 0

We introduce the (n + 1) × (n + 1) matrix M defined by   Δ(0) ψ2 (0) , M=

ψ1 , IRn  0

for q2 = 0.

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

2843

then it follows from the systems (3.20) and (3.21) that the vectors h k(q1 ,q2 ,l) (0) satisfy  k−1  E (q1 ,q2 ,l) k M hˆ (q1 ,q2 ,l) (0) = k−1 v(q 1 ,q2 ,l) where hˆ k(q1 ,q2 ,l) (0) =

  h k(q1 ,q2 ,l) (0) 0

.

The uniqueness of h k(q1 ,q2 ,l) (0) is a consequence of the next lemma. Lemma. M is an invertible matrix. Proof. Let x ∈ Rn , η ∈ R be such that   x M = 0, η then



Δ(0)x + ηψ2 (0) = 0

ψ1 , Rn x = 0,

and then from the fact that ψ2 (0)Δ(0) = 0, we have η = 0, and the above system becomes  Δ(0)x = 0

ψ1 , Rn x = 0. From the fact that dim ker Δ(0) = 1 and since we have Δ(0)φ1 (0) = 0, then the first equation implies that x ∈ span{φ1(0)}, it follows that there exists c ∈ R such that x = cφ1 (0), and from the second equation we have ψ1 , cφ1  = 0 which yields c = 0 and finally we have x = 0. Consequently M is invertible. The relations (3.13) and (3.14) are a consequence of the fact that N 1 (ξ, α) is the homogeneous part of second degree of Φ(θ )Ψ (0)R 1 (ξ, α), which completes the proof of the theorem.  3.2. Normal forms To compute the normal form of Eq. (3.1) on the center manifolds already computed, we will use the reduced equation on center manifolds obtained, i.e. d z(t) = Bz(t) + H (z(t), α) dt   where B = 00 10 , z ∈ R2 , α ∈ R p and

(3.22)

H (z, α) = Ψ (0) [(L(α) − L(0))(Φz + h(z, α)) + F(Φz + h(z, α), α)] . The basic idea of normal form theory consists of employing successive, near identity nonlinear transformations to eliminate the so-called non-resonant nonlinear terms, and the terms called resonant which cannot be eliminated are remained in normal forms. Assume that by a nonlinear transformation z = u + T (u, α), z, u ∈ R3 , α ∈ R p

(3.23)

2844

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

the above system is transformed to its normal form d u(t) = Bu(t) + N(u(t), α). dt

(3.24)

Then if we replace (3.23) and (3.24) in Eq. (3.22) we obtain the following formula DT (u, α)Bu − BT (u, α) = H (u + T (u, α), α) − DT (u, α)N(u, α) − N(u, α). (3.25) In view of the regularity of h, we restrict our attention to the terms of degree lower than m for each m ≥ 2: Let T (u, α) =

m  

q

k=2 (q,l)∈Dk

q

k T(q,l) u 11 u 22 α l ,

N(u, α) =

k T(q,l)

=

k,1  T(q,l) , k,2 T(q,l)

 k N(q,l)

=

k,1  N(q,l) k,2 N(q,l)

q

k=2 (q,l)∈Dk

H (u + T (u, α), α) − DT (u, α)N(u, α) = 

m  

m 



k=2 (q,l)∈Dk

and

k C(q,l)

q

k N(q,l) u 11 u 22 α l ,

q

q

k C(q,l) u 11 u 22 α l ,

 =

k,1 C(q,l)



k,2 C(q,l)

.

Then from the fact DT (u, α)Bu =

m 



k=2 (q,l)∈Dk ,q1 =0

=

m 



k=2 (q,l)∈Dk ,q2 =0

and

 BT (u, α) =

q −1 q2 +1 l u2 α

k q1 T(q,l) u 11

q

q

(q1 + 1)T(qk 1+1,q2 −1,l) u 11 u 22 αl ,

 T 2 (u, α) , 0

the coefficients of the nonlinear transformation and the normal form satisfy the following system for (q1 , q2 , l) ∈ Dk : k,1 k,1 (q1 + 1)T(qk,1 − T(qk,2 = C(q − N(q , 1 +1,q2 −1,l) 1 ,q2 ,l) 1 ,q2 ,l) 1 ,q2 ,l) k,2 = C(q (q1 + 1)T(qk,2 1 +1,q2 −1,l) 1 ,q2 ,l) k,2 k,1 k,1 −T(q1 ,0,l) = C(q1 ,0,l) − N(q1 ,0,l) , k,2 k,1 = C(q . N(q 1 ,0,l) 1 ,0,l)

−

k,2 N(q , 1 ,q2 ,l)

for q2 = 0,

for q2 = 0,

k Remark 3.4. Note that for a fixed k ≥ 2, the coefficients C(q , (q1 , q2 , l) ∈ Dk , are given 1 ,q2 ,l) j

j

in terms of T(q1 ,q2 ,l) , N(q1 ,q2 ,l) , j ∈ {2 · · · k − 1}, because of T (0, 0) = DT (0, 0) = N(0, 0) = D N(0, 0) = 0. Finally, the algorithm for computation of normal forms is the following:

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851 j

2845

j

Having computed N(q1 ,q2 ,,l) , T(q1 ,q2 ,l) for all j ∈ {2, . . . , k − 1} and (q1 , q2 , l) ∈ D j then k N(q , T(qk 1 ,q2 ,l) are given as follows: 1 ,q2 ,l) From the first and the third equations, and since there are no interactions for choosing the k,1 k,1 , we can choose N(q = 0 for all (q1 , q2 , l) ∈ Dk . From the second coefficients N(q 1 ,q2 ,l) 1 ,q2 ,l) k,2 and the third equations we can only choose N(q = 0 for q2 ∈ {0, 1}. Since T(qk,2 is well 1 ,q2 ,l) 1 ,0,l)

k,2 k,2 determined from the third equation, it follows that N(q and N(q are given in a unique 1 ,0,l) 1 ,1,l) way by k,2 k,2 = C(q − (q1 + 1)T(qk,2 N(q 1 ,1,l) 1 ,q2 ,l) 1 +1,0,l)

and k,2 k,1 N(q = C(q . 1 ,0,l) 1 ,0,l)

The NLT associated with those choices can be given as follows: k,1 T(qk,2 = −C(q , 1 ,0,l) 1 ,0,l)

T(qk,2 = 1 ,q2 ,l)

k,2 C(q 1 −1,q2 +1,l)

q1

for q1 = 0.

k,2 T(0,q = 0 for q2 = 0. 2 ,l)   k,1 k,2 + T C(q (q1 −1,q2 +1,l) 1 −1,q2 +1,l) T(qk,1 = 1 ,q2 ,l) q1

for q1 = 0.

k,1 T(0,q = 0 for q2 = 0. 2 ,l)

The normal from for parameterized equation (3.2) is given by u˙ 1 = u 2 ! m  u˙ 2 = k=2



(q1 ,0,l)∈Dk

q1 l k,2 N(q u ,0,l) 1 α 1



+ u2

(q1 ,1,l)∈Dk

" q k,2 N(q u 1 u 2 αl 1 ,1,l) 1

+ O(|u 1 | + |u 2 | + |α|)m+1 . If we neglect the third order terms from this equation we obtain u˙ 1 = u 2 u˙ 2 = μ1 u 1 + μ2 u 2 + au 21 + bu 1 u 2 , where μ1 =

 (1,0,l)∈D2

a=



(2,0,0)∈D2

2,2 N(1,0,l) αl ,

2,2 N(2,0,0)

μ2 =

and b =

 (0,1,l)∈D2

2,2 N(0,1,l) αl ,



(1,1,0)∈D2

2,2 N(1,1,0)

This later equation is interesting when ab = 0 since for the Bogdanov singularity the quadratic terms are needed to determine the local phase portrait (see [3,4]).

2846

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

4. Outline of the symbolic computation All the formulas presented in the previous section are given explicitly in terms of the coefficients of the original differential equations, and thus can be easily implemented on a symbolic computation system. The symbolic manipulation language Maple has been used to code these explicit formulas. In this section, we shall outline the computer programs. 4.1. Create the input file (a): Set and define the variables: n: the dimension of the system. m: the number of delays. r [l], l = 1 · · · m: the delays. phi[i, j ] (i = 1, 2; j = 1 · · · n): the j -th component of the i -th element of the basis Φ, of the generalized subspace associated with the critical eigenvalues. m0: the order of the center manifold to be computed. (b): Create the vector field of the original system: In order to create the vector field of the system, we must consider the parameters αl , l = 1 · · · p, as new dependent variables as follows: We put x[1, n + l] := alpha[l], l = 1 · · · p, where x[i, j ] denotes the j -th component of the dependent variables associated with the delay r [i ] (see the examples for more clarification). L P[l, i, j ]: the (i, j )-th coefficient of the linear part of the vector field associated with the delay r [l] of the system. N L[i ]: the i -th component of the nonlinearity of the vector field of the system. 4.2. Compute the center manifold and normal form (a): compute the basis Ψ . (b): compute recursively the coefficients h i(q1 ,q2 ,q3 ,l) , q1 + q2 + q3 + l = i, i = 2, . . . , j of center manifolds:  j −1 • Compute the homogeneous part H j of F(Φξ + i=2 h i (ξ, α), α). • Compute the coefficients of the nonlinear terms N j −1 (ξ, α)(θ ) and R j −1 (ξ, α). j • Compute the initial data h (q1 ,q2 ,q3 ,l) (0) using formulas (3.16). j

• Compute the vector h (q1 ,q2 ,q3 ,l) (θ ) using formulas (3.13)–(3.15). • Construct the j -th homogeneous part h j (ξ, α)(θ ) of center manifolds. (c): compute the normal forms N associated up to the desired order. (d): write the reduced equation on the center manifold and the normal form N in the output file. 5. Examples To show the applicability of the method and the efficiency of the Maple programs, two examples are presented. They are chosen from the problems which have been studied before by many authors since [8,12], in order to verify the results obtained using the programs developed in this paper. All the results including the center manifolds, normal forms obtained by executing the Maple programs on a PC (Pentium 2 - 400 MHz), are given in the Appendices. We describe also how to create the input file of each example (see the Appendices).

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

2847

5.1. Example 1 We consider the scalar singularly perturbed delay equation of the form ε

dx(t) = −x(t) + f (x(t − 1), λ) dt

(5.1)

where ε is a real positive number assumed to be small, λ is a real parameter, and f is an element of C k (R, R), k ≥ 3; more specifically f (x, λ) = −(1 + λ)x + ax 2 + bx 3 + cx 4 + o(x 4 ),

for some a, b, c ∈ R.

With the change of variables w1 (t) = x(−εr t)

and w2 (t) = x(−εr t − 1 − εr ),

Eq. (5.1) reads w˙ 1 (t) = r w1 (t) − r f (w2 (t − 1), λ) w˙ 2 (t) = r w2 (t) − r f (w2 (t − 1), λ). If r = 1 and λ = −2, then the characteristic equation of the linear part is det Δ(μ) = (μ − 1)2 + (μ − 1)e−μ = 0 and it has zero as a double root and  theremaining roots have negative real part. If we let α1 = r − 1, α2 = λ + 2 and w(t) := ww12 (t) (t) then this system can be read as w(t) ˙ = L 1 (α)w(t) + L 2 (α)w(t − 1) + F(wt , α), where

and

  1 0 L 1 (α) = (1 + α1 ) , 0 1

  0 1 L 2 (α) = (1 + α1 )(1 + α2 ) 0 1

  −(1 + α1 )(aφ2 (−1)2 + bφ2 (−1)3 ) F(φ, α) = + O(|φ(−1)|3 ). −(1 + α1 )(aφ2 (−1)2 + bφ2 (−1)3 )

Then executing the Maple program developed in this paper yields the reduced equation on the center manifold and the associated normal form up to fourth order in Appendix A. Note that in the critical case (r = 1, λ = −2), the same results for the center manifold, but only up to second order, have been obtained by Hale and Huang [12]. 5.2. Example 2 Consider the RFDE x(t) ˙ = αx(t) + f (x(t − 1)), where f ∈ C(R, R) and f (0) = 0. Linearization around the zero equilibrium yields x(t) ˙ = αx(t) + f  (0)x(t − 1).

2848

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

This equation has a double eigenvalue zero at (α, f  (0)) = (1, −1). If we introduce the small parameters α1 and α2 through α1 = α − 1

and α2 = f  (0) + 1,

then the RFDE can be read as x(t) ˙ = x(t) − x(t − 1) + α1 x(t) + α2 x(t − 1) + g(x(t − 1)), where 

f (0) 2 x + O(|x|2 ). g(x) = 2 Again we execute the Maple program developed in this paper to find the reduced equation on the center manifolds and the associated normal form, up to third order, given in Appendix B, where  μ1 = f 2(0) . Our result obtained for the center manifold up to second order confirms the one established in [8], p. 281. 6. Conclusion Maple computer programs have been developed for computing the explicit center manifolds of Bogdanov–Takens bifurcations in retarded functional differential equations with parameters. As a consequence, the normal forms on the center manifolds are derived explicitly in terms of the original equation. The formulas are given in an explicit iterative procedure, and thus are implemented using computer algebra. Examples are presented to verify the method and to show that the method is computationally efficient. Appendix A The input file for example 1: n:=2: m:=2: phi[1, 1]:=1: phi[1, 2]:=1: phi[2, 1]:=1 + theta: phi[2, 2]:=1 + theta: r[1]:=0: r[2]:=1: LP[1, 1, 1]:=1: LP[1, 1, 2]:=0: LP[1, 2, 1]:=0: LP[1, 2, 2]:=1: LP[2, 1, 1]:=0: LP[2, 1, 2]:= − 1: LP[2, 2, 1]:=0: LP[2, 2, 2]:= − 1: NL[1]:=x[1, 3]∗x[1, 1] + x[1, 3]∗x[2, 2] + x[1, 4]∗x[2, 2] + (x[1, 3]∗ x[1, 4])∗x[2, 2] − a∗ (x[1, 3] + 1)∗ (x[2, 2]∗∗2) − b∗ (x[1, 3] + 1)∗ (x[2, 2]∗∗3):

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

2849

NL[2]:=x[1, 3]∗x[1, 2] + x[1, 3]∗x[2, 2] + x[1, 4]∗x[2, 2] + (x[1, 3]∗x[1, 4])∗x[2, 2] − a∗ (x[1, 3] + 1)∗ (x[2, 2]∗∗2) − b∗ (x[1, 3] + 1)∗ (x[2, 2]∗∗3): m0:=3: The reduced equation on a center manifold up to the fourth order is  172 2 112 32 2 146 2 4 α1 − α1α2 + α2 − α1 α2 dxi 1/dt = ξ 2 + − α1 − 3 405 405 405 405   32 4 32 14 8 80 28 + α1α22 ξ 2 + − α1 − α2 − α12 − α1α2 − α22 − α12 α2 405 3 3 27 27 27 27   32 5329 6656 14 34 aα1 − α2a + aα12 − α1α2a − α1α22 ξ 1 + − 27 405 405 54675 54675   256 128 4672 224 578 + aα22 ξ 22 + aα1 − α2a + aα12 − α1α2a 54675 405 405 3645 3645   224 574 14 700 4 − aα22 ξ 1ξ 2 + a + 4aα1 + α2a + aα12 + α1α2a 3645 3 9 243 243    49 2336 2 512 2 aα22 ξ 12 + a α1 − a α2 ξ 23 + 243 54675 54675   64 2 20032 2 2336 2 a + a α1 + a α2 ξ 1ξ 22 + 405 54675 54675   32 128 2 38 2 224 2 146 a + a α1 + a α2 + bα1 − bα2 ξ 12 ξ 2 + 405 135 1215 135 135   28 2 4 448 2 14 98 2 40 + − a + b− a α1 − a α2 + bα1 + bα2 ξ 13 27 3 243 243 9 9   256 3 4 1024 3 2336 32 a ξ2 + a ξ 1ξ 23 + − a3 + ba ξ 12 ξ 22 + 54675 54675 54675 135     4 448 3 64 49 14 a + ba ξ 13 ξ 2 + a 3 − ba + c ξ 14 + − 3645 135 243 9 3   16 2 86 2 56 73 2 16 α1 + α1α2 − α2 + α1 α2 − α1α22 ξ 2 dxi 2/dt = 2α1 + 135 135 135 135 135   7 2 14 2 16 2 40 7 + 4α1 + 2α2 + α1 + α1α2 + α2 + α1 α2 + α1α22 ξ 1 9 9 9 9 9   16 5329 128 17 3328 aα1 + α2a − aα12 + α1α2a − aα22 ξ 22 + 135 135 36450 18225 18225   64 2336 112 112 289 aα1 + α2a − aα12 + α1α2a + aα22 ξ 1ξ 2 + − 135 135 1215 1215 1215   350 49 7 287 aα12 − α1α2a − aα22 ξ 12 + −2a − 6aα1 − α2a − 3 81 81 162   1168 2 256 2 a α1 + a α2 ξ 23 + − 18225 18225   32 2 10016 2 1168 2 + − a − a α1 − a α2 ξ 1ξ 22 135 18225 18225

2850

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

  16 64 2 19 2 112 2 73 a − a α1 − a α2 − bα1 + bα2 ξ 12 ξ 2 + − 135 45 405 45 45   7 224 2 49 2 20 128 3 4 14 2 a − 2b + a α1 + a α2 − bα1 − bα2 ξ 13 − a ξ2 + 9 81 81 3 3 18225     512 3 1168 3 16 224 3 32 a ξ 1ξ 23 + a − ba ξ 12 ξ 22 + a − ba ξ 13 ξ 2 − 18225 18225 45 1215 45   49 3 7 a + ba − 2c ξ 14 . + − 162 3 The normal form is: order = 4, dz1/dt = z2  12553 2 1634 1994 2 1304 36224 3 2 dz2/dt = − α1 α2 − α1 − α1α2 − α1α22 − α1 − α1 405 135 135 405 1215 3   86 2 512 3 4 12772 416 1586 α2 − α2 − α2 ξ 2 + − α1α22 − α12 α2 − α12 + 1215 3 135 405 405 27  8 36224 3 14 2 4 512 3 272 α1α2 − α1 − α1 − α2 − α2 + α2 ξ 1 − 27 3 1215 27 3 1215   1225 122 68728 8 484 4328 2 2 a+ α1α2a − aα2 + α2a + aα1 + aα1 ξ 1ξ 2 + 3 27 81 135 405 135   24502 4 148 105052 1547 14 2 2 + α1α2a + a + aα1 + aα1 + α2a − aα2 ξ 12 1215 3 9 1215 9 1215   25163 2 484 2 110 2 226 2927 a + 4b + a α2 + bα2 − a α1 + bα1 ξ 12 ξ 2 + − 135 81 45 405 45   1814 2 4 28 2 200 21752 2 14 + a α2 + bα2 + b − a + bα1 − a α1 ξ 13 1215 9 3 27 9 1215     16 14 220 3 104 907 3 4 3 a − ba + c ξ 1 ξ 2 + − a + c − ba ξ 14 . + − 243 15 3 1215 3 9 Appendix B The input file for example 2: n:=1: m:=2: phi[1, 1]:=1: phi[2, 1]:=theta: r[1]:=0: r[2]:=1: LP[1, 1, 1]:=1: LP[2, 1, 1]:= − 1: NL[1]:=x[1, 2]∗x[1, 1] + x[1, 3]∗ x[2, 1] + mu1∗ x[2, 1]∗∗2 + mu2∗ x[2, 1]∗∗3: m0:=2:

R. Qesmi et al. / Nonlinear Analysis 66 (2007) 2833–2851

2851

The reduced equation on a center manifold up to the third order is   2 2 2 2 4 2 α1 + α2 ξ 1 + μ1ξ 22 − μ1ξ 1ξ 2 + μ1ξ 12 dxi 1/dt = ξ 2 − α2ξ 2 + 3 3 3 3 3 3 2 2 − μ2ξ 23 + 2μ2ξ 1ξ 22 − 2μ2ξ 12 ξ 2 + μ2ξ 13 3 3 2 dxi 2/dt = −2α2ξ 2 + (2α1 + 2α2)ξ 1 + 2μ1ξ 2 − 4μ1ξ 1ξ 2 + 2μ1ξ 12 − 2μ2ξ 23 + 6μ2ξ 1ξ 22 − 6μ2ξ 12 ξ 2 + 2μ2ξ 13 . The normal form is: order = 3, dz1/dt = z2    4 20 20 2 2 2 α1 − α2 − α1α2 − α22 z2 + α1 + α2 dz2/dt = 3 3 9 9 3 3     20 20 44 4 2 8 − α1α2 − α22 z1 + − μ1 − μ1α1 + μ1α2 z1z2 + μ1 9 9 3 9 3 3      28 40 16 2 64 2 μ1 − 4μ2 z12 z2 + − μ12 + μ2 z13 . − μ1α1 − μ1α2 z12 + 9 3 3 9 3 References [1] M. Ait Babram, O. Arino, M.L. Hbid, Approximation scheme of a system manifold for functional differential equations, J. Math. Anal. Appl. 213 (1997) 554–572. [2] M. Ait Babram, O. Arino, M.L. Hbid, Computational scheme of a center manifold for neutral functional differential equations, J. Math. Anal. Appl. 258 (2001) 396–414. [3] R.L. Bogdanov, Versal deformation of a singularity of a vector field on the plane in the case of zero eigenvalues, Funct. Anal. Appl. 9 (1975) 144–145. [4] R.L. Bogdanov, Versal deformation of a singularity of a vector field on the plane in the case of zero eigenvalues, Sel. Math. Sov. 1 (1981) 389–421. [5] H.W. Broer, et al., The Hamiltonian double-zero eigenvalues, Fields Inst. Commun. 4 (1995) 1–20. [6] S.N. Chow, C. Li, D. Wang, Normal Forms and Bifurcation of planar Vector Fields, Cambridge University Press, Cambridge, 1994. [7] L.O. Chua, H. Cocubu, Normal forms for nonlinear vector fields—Part I: Theory and algorithm, IEEE Trans. Circuits Syst. 35 (1988) 836–880. [8] O. Diekmann, S.A. Van Gils, S.M. Verduyn Lunel, H.O. Walther, Delay Equations, Functional-,Complex-, and Nonlinear Analysis, Springer-Verlag, New York, 1995. [9] T. Faria, L.T. Magalh´aes, Normal forms for retarded functional differential equations with parameters and applications to Hopf bifurcation, J. Differential Equations 122 (1995) 181–200. [10] T. Faria, L.T. Magalh´aes, Normal forms for retarded functional differential equations and applications to Bogdanov–Takens bifurcation, J. Differential Equations 122 (1995) 201–224. [11] J.K. Hale, S.M. Verduyn Lunel, Introduction to Functional Differential Equations, Springer-Verlag, New York, 1993. [12] J.K. Hale, W. Huang, Period doubling in singularly perturbed delay equations, J. Differential Equations 114 (1994) 1–23. [13] D. Henry, Geometric theory of semi linear parabolic equations, in: Lect. Notes in Math., vol. 840, Springer-Verlag, NY, 1981. [14] H. Keller, Numerical solution of bifurcation and nonlinear eigenvalue problems, in: P. Rabinowitz (Ed.), Applications of Bifurcation Theory, Academic Press, New York, 1977, pp. 359–384. [15] R. Qesmi, et al., A Maple program for computing a terms of a center manifolds, and elements of bifurcations for a class of retarded functional differential equations with Hopf singularity, J. Appl. Math. Comput. 175 (2006) 932–968. [16] A. Vanderbauwhede, G. Ioos, Center manifold theory in infinite dimension, Dynamics report, vol. 2, 1992. [17] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, in: Texts in Applied Mathematics, vol. 2, Springer-Verlag, New York, 1990.