Hopf bifurcation and stability of periodic solutions ... - Semantic Scholar

Comment

Report 2 Downloads 74 Views

Nonlinear Analysis 62 (2005) 141 – 165 www.elsevier.com/locate/na

Hopf bifurcation and stability of periodic solutions for van der Pol equation with time delay夡 Wenwu Yu, Jinde Cao∗ Department of Mathematics, Southeast University, Nanjing 210096, China Received 5 January 2005; accepted 9 March 2005

Abstract In this paper, the van der Pol equation with a time delay is considered, where the time delay is regarded as a parameter. It is found that Hopf bifurcation occurs when this delay passes through a sequence of critical value. A formula for determining the direction of the Hopf bifurcation and the stability of bifurcating periodic solutions is given by using the normal form method and center manifold theorem. © 2005 Published by Elsevier Ltd. Keywords: Van der Pol equation; Time delay; Hopf bifurcation; Periodic solutions

1. Introduction The well-known van der Pol equation, which describes the oscillations is the second-order nonlinear damped system governed by x(t) ˙ = y(t) − f (x(t)), (1.1) y(t) ˙ = −x(t). This model is considered as one of the most intensely studied system in nonlinear system [9,12,23] and has served as a basic model in physics, electronics, biology, neurology and 夡 This work was jointly supported by the National Natural Science Foundation of China under Grant 60373067, the 973 Program of China under Grant 2003CB317004, the Natural Science Foundation of Jiangsu Province, China under Grants BK2003053, Qing-Lan Engineering Project of Jiangsu Province, China. ∗ Corresponding author. Tel.: +86 25 83792315; fax: +86 25 83792316. E-mail address: [email protected] (J. Cao).

0362-546X/$ - see front matter © 2005 Published by Elsevier Ltd. doi:10.1016/j.na.2005.03.017

142

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

so on. Many efforts have been made to ﬁnd its approximate solutions [1,10,24]. As we know, in ordinary differential equation, a non-constant periodic solution can arise a Hopf bifurcation. This occurs when a eigenvalue crosses the imaginary axis from left to right as a real parameter in the equation passes through a critical value [15,21]. In [10], the authors proposed a class of relaxation algorithms for ﬁnding the periodic steady-state solution of a van der Pol oscillation. In [1], the periodic solution of van der Pol equation is given in the form of a series converging for all values of the damping parameter. Recently, dynamical systems with time delays have been found in neural networks [2–8,11,18,20]. It is worth noting that the dynamical characteristics (including stable, unstable, oscillatory, and chaotic behavior) of neural networks with time delays have become a subject of intense research activities ([2–8,11,18,20] and the references cited therein), and neural networks involving persistent oscillations such as limit cycle may be applied to pattern recognition and associative memory. Thus it is of great interest to understand the mechanism of neural networks that cause and sustain such periodic activities. However, neural networks are complex and large-scale nonlinear dynamical systems. For simplicity, many researchers have directed their attention to the study of simple systems. This is still useful since the complexity found may be carried over to large networks. In [22], Murakaimi introduced a discrete time delay into the van der Pol equation (1.1) and the following delayed differential equation was obtained: x(t) ˙ = y(t − ) − f (x(t − )), (1.2) y(t) ˙ = −x(t − ). And the author discussed in detail the existence of periodic solution by using the center manifold approaches. However, the stability of bifurcating periodic solution was not discussed. Clearly, if =0, (1.2) can reduce to (1.1). In [17], Liao studied Hopf bifurcation and stability of periodic solutions for van der Pol equation with distributed delay. Therefore, dynamical analysis of time-delay systems is an important topic in many ﬁelds [13,14,16,19,25]. In this paper, we will consider the van der Pol equation with a discrete delay, and study the existence of a Hopf bifurcation and the stability of bifurcating periodic solutions of Eq. (1.2). The obtained results ﬁnd that both of them depend on the parameters a and the delay . The organization of this paper is as follows: In Section 2, we will discuss the stability of the trivial solutions and the existence of Hopf bifurcation. In Section 3, a formula for determining the direction of the Hopf bifurcation and the stability of bifurcating periodic solutions will be given by using the normal form method and center manifold theorem introduced by Hassard et al. [14]. In Section 4, numerical simulations aimed at justifying the theoretical analysis will be reported.

2. Stability of trivial solutions and existence of Hopf bifurcation In this section, we consider the following delayed differential equation: x(t) ˙ = y(t − ) − f (x(t − )), y(t) ˙ = −x(t − ),

(2.1)

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

143

where f (x) = ax + bx 2 , a and b are positive and real. The linear equation of (2.1) are as follows: x(t) ˙ = −ax(t − ) + y(t − ), (2.2) y(t) ˙ = −x(t − ). Clearly, the point (0, 0) is the unique equilibria of (2.2). The characteristic equation of the linearized system (2.2) is

+ ae− det e−

−e−

= 0.

(2.3)

By simple calculation, we obtain the following characteristic equation:

2 + a e− + e−2 = 0.

(2.4)

Lemma 2.1. (i) If a > 2 holds, then (2.4) has a pair of purely imaginary roots ±iw0l when = nl (l = 1, 2; n = 0, 1, 2, . . .), where w0l =

nl =

a±

√

a2 − 4 , 2

2n + /2 (l = 1, 2; n = 0, 1, 2, . . .). w0l

(2.5) (2.6)

(ii) If 0 < a 2 holds, then (2.4) has a pair of purely imaginary roots ±iw0 when = j (j = 0, 1, 2, . . .), where w0 = 1,

j = arcsin

(2.7) a + 2j . 2

(2.8)

Proof. Suppose that = iw, with w >0, is a root of Eq. (2.4), then we obtain −w2 + aiw(cos(w ) − i sin(w )) + (cos(2w ) − i sin(2w )) = 0. Separating the real and imaginary parts, we have

w2 − aw sin(w ) − cos(2w ) = 0, aw cos(w) − sin(2w ) = 0.

(2.9)

Which can be rewritten as the following two equations:

w 2 − aw sin(w ) − cos(2w ) = 0, cos(w ) = 0,

(2.10)

144

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

or

w 2 − aw sin(w ) − cos(2w ) = 0, aw . sin(w ) = 2

(2.11)

Next we discuss it in two cases: Case I: From cos(w)=0, it follows that cos(2w )=−1. If sin(w )=−1, then substituting this into the ﬁrst equation of (2.10), we have w 2 + aw + 1 = 0. Suppose w1 and w2 are the two roots of this equation. Clearly, w1 + w2 = −a < 0,

w1 w2 = 1 > 0.

The equation have two negative roots. Since we suppose w > 0, we choose sin(w ) = 1. Substituting this into the ﬁrst equation of (2.10), we obtain w2 − aw + 1 = 0.

(2.12)

√ For the case a > 2, we have two roots of Eq. (2.12): w0l = a ± a 2 − 4/2, thus nl = 2n + /2/w0l (l = 1, 2; n = 0, 1, 2, . . .). The proof is completed. Case II: From sin(w) = aw/2, we have cos(2w ) = 1 − 2 sin2 (w ) = 1 − a 2 w 2 /2. Substituting these into the ﬁrst equation of (2.11), from 0 < a 2, we obtain w0 = 1 and j = arcsin a/2 + 2j . The proof is completed. Denote

() = () + iw(), is the root of Eq. (2.4) satisfying

(nl ) = 0,

w(nl ) = w0l ,

where nl is deﬁned by (2.6).

Lemma 2.2. If a > 2 holds, then we have d(nl ) > 0. d

(2.13)

Proof. Taking the derivative of with respect to in (2.4), we have 2 + a e− + a (−)e− + a (− )e− − 2( + )e−2 = 0, it follows that: d() a 2 e− + 2e−2 = . d 2 + ae− − a e− − 2e−2

(2.14)

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

145

For the sake of simplicity, we denote w0l and nl by w, , respectively, then −aw 2 e−iw + 2iwe−2iw d = d 2iw + ae−iw − aiw e−iw − 2e−2iw = =

−aw2 [cos(w)−i sin(w)]+2iw[ cos(2w)−i sin(2w)] 2iw+a[cos(w)−i sin(w)]−aiw[cos(w)−i sin(w)]−2[cos(2w)−i sin(2w)] [−aw2 cos(w)+2w sin(2w)]+i[aw2 sin(w)+2w cos(2w)] [a cos(w)−aw sin(w)−2 cos(2w)]+i[2w−a sin(w)−aw cos(w)+2 sin(2w)] .

Let Q = [a cos(w) − aw sin(w ) − 2 cos(2w )]2 + [2w − a sin(w ) − aw cos(w ) + 2 sin(2w )]2 . Q Re

d d

= [−aw 2 cos(w ) + 2w sin(2w )] × [a cos(w) − aw sin(w ) − 2 cos(2w )] + [aw2 sin(w ) + 2w cos(2w )] × [2w − a sin(w ) − aw cos(w ) + 2 sin(2w )]

(2.15)

d = (aw 2 − 2w)(2w − a) Q Re d =nl = w(2aw 2 − a 2 w − 4w + 2a)

= w(2a 2 w − 2a − a 2 w − 4w + 2a) = w 2 (a 2 − 4) > 0, this completes the proof. Denote

() = () + iw(), is the root of Eq. (2.4) satisfying

(j ) = 0,

w(j ) = w0 ,

where j is deﬁned by (2.8).

Lemma 2.3. If 0 < a 2 holds, then d(j ) > 0. d

(2.16)

146

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

Proof. If 0 < a < 2, from (2.15) we have Q Re

d = [−a cos + 2 sin(2)][a cos − a 2 /2 − 2(1 − a 2 )/2] d =j + [a 2 /2 + 2 − a 2 ][2 − a 2 /2 − a cos + 2 sin(2)] = a cos (a cos + a 2 /2 − 2) + (2 − a 2 /2)(2 − a 2 /2 + a cos ) = a 2 cos2 + a cos (a 2 /2 − 2) + (2 − a 2 /2)a cos + (2 − a 2 /2)2 = a 2 (1 − a 2 /4) + (4 + a 4 /4 − 2a 2 ) = a 2 − a 4 /4 + 4 + a 4 /4 − 2a 2 = 4 − a 2 > 0.

If a = 2 holds, then (2.4) is equivalent to the following equation:

+ e− = 0.

(2.17)

Here, w = 1, = /2. Using the same method as Lemma 2.2, we can obtain

d Re = 1 > 0. d =j This completes the proof.

Lemma 2.4. For Eq. (2.4). (I) If a > 2, all the roots of Eq. (2.4) have strictly negative real parts for ∈ [0, 01 ), and Eq. (2.4) has a pair of imaginary roots ±iw01 and all the other roots have strictly negative real parts when = 01 , as well as Eq. (2.4) has at least a pair of roots with positive real parts when > 01 . (II) If 0 < a 2, all the roots of Eq. (2.4) have strictly negative real parts for ∈ [0, 0 ), and Eq. (2.4) has a pair of imaginary roots ±iw0 and all the other roots have strictly negative real parts when = 0 , as well as Eq. (2.4) has at least a pair of roots with positive real parts when > 0 . Proof. (I) Obviously, the roots of Eq. (2.4) continuously depend on the parameter . When = 0, we know that (2.4) has two roots 1 and 2 , the product of the roots satisfy

1 + 2 = −a(a > 0),

1 2 = 1.

We obtain that 1 and 2 both have negative real parts. 01 is the smallest positive value when Eq. (2.4) has a pair of purely imaginary roots. Since the roots of Eq. (2.4) continuously

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

147

depend on the parameter , we know that the roots of Eq. (2.4) have negative real parts when ∈ [0, 01 ). Next we show that Eq. (2.4) has a pair of imaginary roots ±iw01 and all the other roots have strictly negative real parts when = 01 . Suppose, on the contrary, that = u + iv with u >0 is a root of Eq. (2.4) with = 01 . Since the roots of Eq. (2.4) continuously depend on the parameter , there exists ∈ (0, 01 ) such that (2.4) has a purely imaginary root at = , which contradicts with the fact that 01 is the smallest such . Thirdly, we will show Eq. (2.4) has at least a pair of roots with positive real parts when > 01 . From Lemma 2.2, we have d(01 )/d > 0. The proof is completed. (II) The approach is the same as (I), we omit it. Theorem 2.1. For Eq. (2.1). (I) If a > 2 and from Lemma 2.4(I), we obtain its zero solution is asymptotically stable for ∈ [0, 01 ), and unstable for > 01 , and Eq. (2.1) undergoes a Hopf bifurcation at the origin when = 01 . That is, system (2.1) has a branch of periodic solutions bifurcating from the zero solution near = 01 . (II) If 0 < a 2 and from Lemma 2.4(II), we obtain its zero solution is asymptotically stable for ∈ [0, 0 ), and unstable for > 0 , and Eq. (2.1) undergoes a Hopf bifurcation at the origin when = 0 . That is, system (2.1) has a branch of periodic solutions bifurcating from the zero solution near = 0 . Remark. In [17], it is simpler than our models since the characteristic equation of ours is transcendental equation corresponding to polynomial equation of Liao. So he discussed the local stability and existence of Hopf bifurcation using Routh–Hurwitz criteria. The stability and existence of Hopf bifurcation which is studied in our paper are not as simple as his.

3. Stability of bifurcating periodic solutions In this section, formulae for determining the direction of Hopf bifurcation and stability of bifurcating periodic solutions of system (2.1) at 0 shall be presented by employing the normal form method and center manifold theorem introduced by Hassard et al. [14]. For convenience, let t = s , x(s ) = x1 (s), y(s ) = x2 (s) and = 0 + , ∈ R. Denote t = s, then system (2.1) is equivalent to the system: x˙1 (t) = (0 + )(−ax 1 (t − 1) + x2 (t − 1) − bx 21 (t − 1)), (3.1) x˙2 (t) = −(0 + )x1 (t − 1). Its linear part is given by x˙1 (t) = (0 + )(−ax 1 (t − 1) + x2 (t − 1)), x˙2 (t) = −(0 + )x1 (t − 1). The nonlinear part of (3.1) is −bx 21 (t − 1) f (, ut ) = (0 + ) . 0

(3.2)

(3.3)

148

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

Denote C k [−1, 0] = {| : [−1, 0] → R 2 , each component of has k order continuous derivative}. For convenience, denote C[−1, 0] by C 0 [−1, 0]. The solutions map continuous initial data into R 2 . We are interested in periodic solutions. For () = (1 () 2 ())T ∈ C[−1, 0], deﬁne an operator

−a L = (0 + ) −1

1 0

1 (−1) , 2 (−1)

(3.4)

where L is a one-parameter family of bounded linear operators in C[−1, 0] → R 2 . By the Riesz representation theorem, there exists a matrix whose components are bounded variation functions (, ) in [−1, 0] → R 2 , such that L =

0 −1

d (, )().

In fact, we choose

(, ) = (0 + )

−a −1

1

( + 1) 0

(3.5)

(where () is Dirac function), then (3.4) is satisﬁed. For ∈ C 1 [−1, 0], deﬁne d() A() = and

, −1 < 0, 0d −1 d (, )(), = 0

 0  , −1 < 0,  0 R() = 2  −b1 (−1)  (0 + ) , = 0. 0

(3.6)

(3.7)

In order to conveniently study Hopf bifurcation problem, we transform system (3.1) into a operator equation of the form: u˙t = A()ut + Rut ,

(3.8)

where u = (x1 , x2 )T . As in [13], ut = u(t + ), ∈ (−1, 0]. The adjoint operator A∗ of A is deﬁned by ∗

A () =

d(s) , 0 < s 1, − 0 ds T −1 d (s, )(−s), s = 0

where T is the transpose of the matrix .

(3.9)

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

149

The domains of A and A∗ are C 1 [−1, 0] and C 1 [0,1], respectively. For ∈ C[−1, 0] and ∈ C[0, 1]. In order to normalize the eigenvectors of operator A and adjoint operator A∗ , we need to introduce the following bilinear form: 0 T ¯ (0) · (0) − , = ¯ ( − ) d ()( ) d . (3.10) =−1 =0

Here () = (, 0), C 2 is complex plane. And for c and d in C 2 , c · d means where ci and di are components of c and d, respectively. Then, as usual, , A = A∗ ,

2

i=1 ci di ,

(3.11)

for (, ) ∈ D(A) × D(A∗ ). We normalize q and q ∗ by the condition q ∗ , q = 1,

q ∗ , q ¯ = 0.

By discussion in Section 2 and transformation t=s , we know that ±i0 w0 are eigenvalues of A(0) and other eigenvalues have strictly negative real parts. Thus they are also eigenvalues of A∗ . Next we calculate the eigenvector q of A belonging to the eigenvalue i0 w0 and the eigenvector q ∗ of A∗ belonging to the eigenvalue −i0 w0 . Let 1 i0 w0 q() = e , −1 < 0. (3.12) From the above discussion, we know that Aq(0) = i0 w0 q(0), −i w −a 1 e 0 0 1 = i , w 0 0 0 e−i0 w0 −1 0 i.e.,

(−a + )e−i0 w0 = iw0 , −e−i0 w0 = iw0 .

Hence, we obtain

= a + iw0 ei0 w0

or

=

i −i0 w0 e . w0

Suppose that the eigenvector q ∗ of A∗ is 1 1 i0 w0 s e q ∗ (s) = , 0 s < 1.

Then we have the following relationship: A∗ q ∗ (0) = −i0 w0 q ∗ (0). i w −a −1 e 0 0 1 = −i0 w0 , 0 1 0 ei0 w0

(3.13)

(3.14)

150

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

i.e.,

(−a − )ei0 w0 = −iw0 , ei0 w0 = −iw0 .

Hence, we obtain

= −a + iw0 e−i0 w0

=

or

i i0 w0 e . w0

(3.15)

Let q ∗ , q = 1. One can obtain , ∗

q , q = q¯∗ (0) · q(0) −

0

=−1 =0

T q¯∗ ( − ) d ()q( ) d

0 1 1 1 i0 w0 −i0 w0 ( −) ¯ ¯ e (1 )e d () d = (1 + ) − ¯

¯ =−1 =0 0 1 1 −a 1 1 ¯ )− = (1 + 0 (1 ¯ ) −1 0 ¯

¯

=−1 =0 × ei0 w0 ( + 1) d d 1 ¯ ) − 1 0 (a + ¯ − )e−i0 w0 = (1 +

¯

¯ = 1. Hence, we have ¯ ) − 0 (a + ¯ − )e−i0 w0 .

¯ = (1 +

(3.16)

Note that ¯ = −, ¯ = −, using the same method it is easy to proof q ∗ , q ¯ = 0, we omit it. Now we obtain q and q ∗ . Next, we study the stability of bifurcating periodic solutions. As in [14], the bifurcating periodic solutions Z(t, ()) has amplitude O() and non-zero Floquet exponent () with (0) = 0. Under our hypotheses, , are given by = 2 2 + 4 4 + · · · , (3.17) = 2 2 + 4 4 + · · · . The sign of 2 indicates the direction of bifurcation while that 2 determines the stability of Z(t, ()). Z(t, ()) is stable if 2 < 0 and unstable if 2 > 0. In the following, we will show how to derive the coefﬁcients in this expansions, but we compute 2 and 2 only. We ﬁrst construct the coordinates to describe a center manifold 0 near = 0, which is a local invariant, attracting a two-dimensional manifold [14]. Let z(t) = q ∗ , ut

(3.18)

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

151

and W (t, ) = ut − 2 Re[z(t)q()].

(3.19)

Where ut is a solution of (3.8). On the manifold 0 : w(t, ) = w(z(t), z¯ (t), ), where W (z, z¯ , ) = W20 ()

z2 z¯ 2 + W11 ()z¯z + W02 () + · · · . 2 2

(3.20)

In fact, z and z¯ are local coordinates of center manifold 0 in the direction of q and q ∗ , respectively. The existence of center manifold 0 enables us to reduce (3.8) to an ordinary differential equation in a single complex variable on 0 . For the solution ut ∈ 0 of (3.8), since = 0, z˙ (t) = q ∗ , u˙ t = q ∗ , Aut + Rut = q ∗ , Aut + q ∗ , Rut = A∗ q ∗ , ut + q ∗ , Rut = i0 w0 z + q¯ ∗ (0) · f (0, W (t, 0) + 2 Re[z(t)q(0)]).

(3.21)

Rewrite (3.21) as z˙ (t) = i0 w0 z + g(z, z¯ ),

(3.22)

where z2 z¯ 2 z2 z¯ + g11 z¯z + g02 + g21 + ···. (3.23) 2 2 2 In the following, our motivation is to expand g in powers of z and z¯ and then obtain, from the coefﬁcients of this expansion, the values of 2 and 2 using algorithm presented by Hassard et al. [14]. Substituting (3.8) and (3.21) into g(z, z¯ ) = g20

W˙ = u˙ t − z˙ q − z˙¯ q, ¯ we have W˙ = u˙ t − z˙ q − z˙¯ q¯ = Aut + Rut − [i0 w0 z + q¯ ∗ (0) · f (z, z¯ )]q − [−i0 w0 z¯ + q ∗ (0) · f¯(z, z¯ )]q¯ = AW + 2A Re(zq) + Rut − 2 Re[q¯ ∗ (0) · f (z, z¯ )q()] − 2 Re[i0 w0 zq()] = AW − 2 Re[q¯ ∗ (0) · f (z, z¯ )q()] + Rut −1 < 0, AW − 2 Re[q¯ ∗ (0) · f (z, z¯ )q()], (3.24) = AW − 2 Re[q¯ ∗ (0) · f (z, z¯ )q()] + f, = 0. Let W˙ = AW + H (z, z¯ , ),

(3.25)

where H (z, z¯ , ) = H20 ()

z2 z¯ 2 + H11 ()z¯z + H02 () + · · · . 2 2

(3.26)

152

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

Taking the derivative of W with respect to t in (3.20), we have W˙ = Wz z˙ + Wz¯ z˙¯ .

(3.27)

Substituting (3.20) and (3.22) into (3.27), we obtain W˙ = (W20 z + W11 z¯ + · · ·)(i0 w0 z + g) + (W11 z + W02 z¯ · · ·)(−i0 w0 z¯ + g). ¯

(3.28)

Then substituting (3.20) and (3.26) into (3.25), we have the following results: z z¯ W˙ = (AW20 + H20 ) + (AW11 + H11 )z¯z + (AW02 + H02 ) + · · · . 2 2 Comparing the coefﬁcients of (3.28) with (3.29), 2

2

(3.29)

(A − 2i0 w0 )W20 () = −H20 (),

(3.30)

AW11 () = −H11 ()

(3.31)

hold. According to (3.21) and (3.22), we know 0 1 −bx 21 (t − 1) ∗ g(z, z¯ ) = q¯ (0) · f (z, z¯ ) = · , 0

¯ ¯

(3.32)

where x1 (t + ) = W (1) (t, ) + z(t)q (1) () + z¯ (t)q¯ (1) (), x2 (t + ) = W (2) (t, ) + z(t)q (2) () + z¯ (t)q¯ (2) (). From (3.32) and (3.21), we have

0 b 2 x (t − 1)

¯ 1 0 b (1) [W (t, ) + z(t)q (1) () + z¯ (t)q¯ (1) ()]2 = −

¯ 0 b (1) z2 z¯ 2 (1) (1) = − [W20 (−1) + W11 (−1)z¯z + W02 (−1)

¯ 2 2 (1) (1) 2 + z(t)q (−1) + z¯ (t)q¯ (−1)] .

g(z, z¯ ) = −

(3.33)

Comparing the coefﬁcients in (3.23) with those in (3.33), it follows that: 20 b (1) [q (−1)]2 ,

¯ 20 b (1) q (−1)q¯ (1) (−1), g11 = −

¯ 20 b (1) [q¯ (−1)]2 , g02 = −

¯ 20 b (1) W20 (−1)q¯ (1) (−1). g21 = −

¯ g20 = −

(3.34)

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

153

In the following, we focus on the computation of W20 (). Eqs. (3.24) and (3.25) imply that H (z, z¯ , ) = − 2 Re(q¯ ∗ (0) · f (z, z¯ )q()) + Rut ¯ ) + Rut = − gq() − g¯ q( 1 1 2 2 g20 z + g11 z¯z + g02 z¯ + · · · q() = − 2 2 1 1 2 2 − g¯ 20 z¯ + g¯ 11 z¯z + g¯ 02 z + · · · q( ¯ ) + Rut . (3.35) 2 2

2 From (3.33), Rut = 0 −bx 10(t−1) = ¯0g , when =0 at =0. Comparing the coefﬁcients in (3.26) with those in (3.35), we can obtain that ¯ ), −g20 q() − g¯ 02 q( −1 < 0, H20 () = (3.36)

¯ g20 −g20 q(0) − g¯ 02 q(0) ¯ + , =0 0 and H11 () =

¯ ), −g11 q() − g¯ 11 q( −1 < 0,

¯ g11 −g11 q(0) − g¯ 11 q(0) , = 0. ¯ + 0

Substituting (3.36) into (3.30) and (3.37) into (3.31) respectively, it follows that: ¯ ), W˙ 20 () = 2i0 w0 W20 () + g20 q() + g¯ 02 q( W˙ 11 () = g11 q() + g¯ 11 q( ¯ ). We can easily obtain the solutions of (3.38):  ig20 g¯ 02 −i0 w0 + E e2i0 w0 ,  q(0)ei0 w0 − q(0)e ¯  W20 () = 1 0 w 0 3i0 w0  −i0 w0 + E .  W11 () = g11 q(0)ei0 w0 − g¯ 11 q(0)e ¯ 2 i0 w0 i0 w0

(3.37)

(3.38)

(3.39)

Next we focus on the computation of E1 , from (3.30), we have AW20 (0) = 2i0 w0 W20 (0) − H20 (0), then

0

−a −1

1 W20 (−1) = 2i0 w0 W20 (0) − H20 (0). 0

Substituting (3.36) and (3.39) into (3.40), we have the following relationship: ig20 g¯ 02 −a 1 i0 w0 q(0)e−i0 w0 − q(0)e ¯ + E1 e−2i0 w0 0 −1 0 0 w0 3i0 w0 2

¯ g20 = −2g20 q(0) − g¯ 02 q(0) ¯ + 2i0 w0 E1 + g20 q(0) + g¯ 02 q(0) ¯ − , 0 3

(3.40)

154

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

2i0 w0 + 0 ae−2i0 w0 −0 e−2i0 w0 E1 2i0 w0 0 e−2i0 w0 ig20 g¯ 02 −a 1 i0 w0 + g20 q(0) q(0)e ¯ = 0 q(0)e−i0 w0 − −1 0 0 w0 3i0 w0 1

¯ g20 − g¯ 02 q(0) . ¯ + 0 3

−1 2i0 w0 + 0 ae−2i0 w0 −0 e−2i0 w0 E1 = 0 e−2i0 w0 2i0 w0 g¯ 02 ig20 −a 1 i0 w0 + g20 q(0) q(0)e−i0 w0 − q(0)e ¯ × −1 0 w0 3iw0 1

¯ g20 − g¯ 02 q(0) . ¯ + 0 3

(3.41)

Hence, we know W20 and then we can obtain g21 . The following parameters can be calculated: g21 1 i C1 (0) = g20 g11 − 2|g11 |2 − |g02 |2 + , (3.42) 2w 3 2

2 = −

Re C1 (0) , Re (0 )

2 = 2 Re C1 (0).

(3.43) (3.44)

If you want to know the detail, see appendix. As in [14], we have the following result: Theorem 3.1. Under the condition of Theorem 2.1, (I) = 0 is Hopf bifurcation value of system (3.1). (II) the direction of Hopf bifurcation is determined by the sign of 2 : if 2 > 0, the Hopf bifurcation is supercritical; if 2 < 0, the Hopf bifurcation is subcritical. (III) The stability of bifurcating periodic solutions is determined by 2 : if 2 < 0, the periodic solutions are stable; if 2 > 0, they are unstable.

4. Numerical examples In this section, some numerical results of simulating system (2.1) are presented at different data of a, b and . First, let b = 0. Then system (2.1) is linear, we investigate the Hopf bifurcation at 0 . We ﬁx a = 1.5, then 0 = arcsin(0.75). So we choose = 0.8 < 0 , 0 and 0.9 > 0 , respectively. The corresponding waveform and phase plots are shown in Figs. 1–3. By Theorem 2.1, we know in Fig. 1 its zero solution is asymptotically stable, in Fig. 2 undergoes a Hopf bifurcation at the origin, and in Fig. 3 is unstable.

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

0.6

0.6

0.4

0.4

0.2

0.2 y

0.8

x

0.8

155

0

0

-0.2

-0.2

-0.4

-0.4

-0.6

-0.6 0

50

100 t

150

200

-1

-0.5

0 x

0.5

1

-0.5

0 x

0.5

1

-100

0

100

200

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

y

x

Fig. 1. a = 1.5, b = 0, < 0 .

-0.2

-0.2

-0.4

-0.4

-0.6

-0.6

-0.8

-0.8 0

50

100 t

150

200

-1

150

150

100

100

50

50

0

0

y

x

Fig. 2. a = 1.5, b = 0, = 0 .

-50

-50

-100

-100

-150 0

50

100

150

200

-150 -200

t

x

Fig. 3. a = 1.5, b = 0, > 0 .

156

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

0.5

0.5

0.4

0.4

0.3

0.3 0.2

0.1

y

x

0.2

0.1

0 -0.1

0

-0.2

-0.1

-0.3 0

50

100

-0.2 -0.4

150

-0.2

0

t

0.2

0.4

0.6

0.2

0.4

0.6

x

0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 0

0.5 0.4 0.3 0.2 y

x

Fig. 4. a = 2.5, b = 0, < 01 .

0.1 0 -0.1 50

100 t

150

200

-0.2 -0.4

-0.2

0 x

Fig. 5. a = 2.5, b = 0, = 01 .

Next, let b = 0. Then system√(2.1) is linear, we investigate the Hopf bifurcation at 01 . We ﬁx a = 2.5, 01 = /(a + a 2 − 4). So we choose = 0.75 < 01 , 01 and 0.8 > 01 , respectively. The corresponding waveform and phase plots are shown in Figs. 4–6. By Theorem 2.1, we know in Fig. 4 its zero solution is asymptotically stable, in Fig. 5 undergoes a Hopf bifurcation at the origin, and in Fig. 6 is unstable. Finally, we ﬁx a = 1 and = 0.55 > 0 = arcsin(a/2) and we choose b = 1 and 1.5, respectively. With these parameters, 2 > 0. Hence, by Theorem 3.1, we know that the bifurcating point is supercritical. Correspondingly, 2 = −0.962 and −0.4416, and so these bifurcating periodic solutions are stable, as shown in Figs. 7 and 8.

5. Conclusions The van der Pol equation provides rich dynamical behavior. From the viewpoint of nonlinear systems, their analysis are useful in solving problems of both theoretical and practical

2.5 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5

157

1.5 1 0.5 y

x

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

0 -0.5 -1 -1.5

0

50

100

150

-4

-2

t

0

2

4

x

0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5

y

x

Fig. 6. a = 2.5, b = 0, > 01 .

0

50

100

150

0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5

200

-0.4 -0.2

t

0

0.2

0.4

0.6

x

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

y

x

Fig. 7. a = 1, b = 1, = 0.55 > 0 .

-0.1

-0.1

-0.2

-0.2

-0.3

-0.3

-0.4 0

50

100 t

150

200

-0.4 -0.4

-0.2

Fig. 8. a = 1, b = 1.5, = 0.55 > 0 .

0 x

0.2

0.4

158

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

importance. Although the systems with time delay discussed above are quite simple, they are potentially useful as the complexity found might be carried over to a general van der Pol equation with time delays. By calling the time delay as a parameter, we have shown that a Hopf bifurcation occurs when this parameter passes through a critical value. The direction of Hopf bifurcation and the stability of the bifurcating periodic orbits are also discussed.

Appendix A. In this appendix, we want to derive formula (3.42), (3.43) and (3.44). A.1. Poincaré normal form We assume we are given a 2 × 2 system in the following Poincaré norm form: X˙ = A()X +

[L/2]

Bj ()X|X|2j + o(|X||(X, )|L+1 ),

j =1

= F (X, ), where

(A.1)

() −() (() = () + i()), () () L Re cj () −Im cj () 1j , Bj () = Im cj () Re cj () 2

A() =

(A.2) (A.3)

and F (X, ) is jointly C L+2 in X and . Eq. (A.1) is equivalent to

˙ = () +

[L/2]

cj () | |2j + o(| ||( , )|L+1 ),

(A.4)

j =1

where = X1 + iX2 . We next derive the formula for the initial coefﬁcients in the MacLaurin expansions of = () and T = T (). We begin by rewrite the differential equation (A.4) as

˙ = () +

M

cj ()( ¯ )j ,

(A.5)

j =1

where is a complex variable, (0)=i0 , M 1 is arbitrary, and cj () are complex valued. This canonical is in Poincaré normal form. Observe that if is a solution, then so is ei for any real number , and trajectories of (A.5) are circles with centers at = 0. This simple geometry is reﬂected in efﬁcient computation of MacLaurin expansions of () and T ().

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

¯˙ from (A.5), we obtain Forming ¯ ˙ +   M   d ¯ Re cj ()( ¯ )j . ( ) = 2 ¯ Re () +   dt

159

(A.6)

j =1

The right-hand side of (A.6) is zero if and only if = 0 or Re () +

M

Re cj ()( ¯ )j = 0.

(A.7)

j =1

But if (A.7) holds, then (A.6) implies that ¯ = 2 0, for some 0. Setting ¯ = 2 0 and = () in (A.7) and we obtain Re (()) +

M

Re cj ()2j = 0.

(A.8)

j =1

This equation determines the coefﬁcients in the expansion

=

M j =1

j j + O(M+1 ).

In the following analysis below of the case (0) = 0, the coefﬁcients 1 , 3 , 5 , . . . are shown to vanish, which is a priori obvious from (A.8). Expanding in powers of in (A.8), we ﬁnd that 2  M M (0)  j j  + · · · + Re c1 (0)2 (0) j j + 2 j =1 j =1   M + Re c1 (0)  j j  2 + · · · + Re c2 (0)4 + · · · = 0. (A.9) j =1

At O(), (A.9) implies that (0)1 = 0. Thus,

1 = 0,

(A.10)

since (0) = 0 by hypothesis. Using this results in (A.9), we ﬁnd that at O(2 ),

(0)2 + Re c1 (0) = 0, 2 = −

Re c1 (0) . (0)

(A.11)

Given that (A.8) holds, we may rewrite (A.5) as   M ˙ = i Im () + cj ()2j  . j =1

(A.12)

160

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

Thus,

= e2it/T () , where



2 = Im () + T ( )

M

 cj ()2j  .

(A.13)

j =1

From this equation the coefﬁcients in the expansion T ( ) =

M 2 i i + O(M+1 ) 0 j =0

may be found. Explicitly at O(1), (A.13) yields

0 = 1

(A.14)

and, to high order

0

−

4

i

i

i=1

+ 0

3

2

i

i

+ · · · = (0)(2 + 4 2 )2 +

i=1

+ Im c1 (0)2 + [Im c1 (0)2 + Im c2 (0)]4 + · · · .

(0) 2 4 2 2 (A.15)

Hence, −0 1 = 0, thus

1 = 0.

(A.16)

since 0 > 0 by hypothesis. Then at O(2 ) (A.15) becomes −0 2 = (0)2 + Im c1 (0) then we obtain

2 = −

1 [Im c1 (0) + (0)2 ]. 0

(A.17)

A.2. Stability criteria We shall now apply Floquet theorem to the real 2 by 2 system (A.1). We know that () is C L+1 in ; hence P (t) = y(t, , ()) is C L+1 in t and . Since X = P (t) is a non-constant, T ()-periodic solution of (A.1), P˙ (t) is a non-trivial, T ()-periodic solution of variational system y˙ = A(t, )y, where A(t, ) = jF /jX at (P (t), ()). Following Floquet’s theorem,

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

161

one of the characteristic exponents associated with P is thus 0 mod 2i/T (). Hence, P has 0 and () as a set of characteristic exponents, where we deﬁne T () 1 tr A(s, ) ds, (A.18) () = T ( ) 0 since T () = T (, ()) is C L+1 in and A(t, ) is C L+1 jointly in t and , the function () is C L+1 in . Next we will expand (). If we write = x1 + ix2 , then x˙1 = x1 − x2 + [(Re c1 )x1 − (Im c1 )x2 ]r 2 + O(4 ), x˙2 = x1 + x2 + [(Re c1 )x2 + (Im c1 )x1 ]r 2 + O(4 ) and tr

jF (P (t, ())) = 2(()) + 4[Re c1 (())]2 + O(3 ), jX

where we have selectively used the facts that () = O(2 ) and r 2 = 2 + O(5 ), (L 2). Hence T () 1 tr A(s, ) ds = 2(()) + 4[Re c1 (())]2 + O(3 ). (A.19) T ( ) 0 But

(()) = (0)2 2 + · · · = −Re c1 (0)2 + · · · . Therefore, following Floquet’s theorem, 0 + () = 2 Re c1 (0)2 + O(3 ). Thus, () < 0 for sufﬁciently small if Re c1 (0) < 0, which is just the criterion derived earlier for asymptotic, orbital stability of P (t, ()). However, the equality () < 0 implies that the bifurcating periodic solutions of the system (A.1) are asymptotically, orbitally stable with asymptotic phase. The above computation can carried out to include terms of order 4 , provided L is at least 4. The result is

() = 2 2 + 4 4 + O(5 ), where

2 = 2 Re c1 (0),

4 = 4 Re c2 (0) + 2 Re c1 (0)2 .

(A.20)

A.3. Reduction of two-dimensional systems to Poincaré normal form Before we can apply the bifurcation formulae, derived in the preceding sections, to various model systems, we must show how general autonomous systems, satisfying the hypothesis

162

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

of Hopf Theorem, can be transformed into Poincaré normal form. We begin with the single complex equation z˙ = z + g(z, z¯ , ),

(A.21)

where

g(z, z¯ , ) =

gij ()

2 i+j L

zi z¯ j + O(|z|L+1 ) i!j !

(A.22)

and

() = () + i(). We desire to transform (A.21) by means of a transformation: j

z = + ( , ¯ , ) = +

2 i+j L

ij ()

i ¯ i!j !

(A.23)

into the Poincaré normal form

˙ = () +

[L/2]

cj () | |2j + o(| ||( , )|L+1 ) = () + ( , ¯ , ).

(A.24)

j =1

By the chain rule z˙ = ˙ + ˙ + ¯ ˙¯ or ¯ ¯ − = g( + , ¯ + ¯ ) − ( + + ¯ ¯ + ¯ ).

(A.25)

From this relation the coefﬁcients ij in (A.23) can be determined recursively. In powers of and ¯ the left-hand side of (A.25) can be written as 2 i+j L

(ij )(i + j ¯ − )

j

i ¯ . i!j !

(A.26)

Note that the expansion of the right-hand side of (A.25) to order k = 2 is independent of the ij and to order k(k = 3, . . . , L) involves exactly the coefﬁcients ij for 2 i + j k. Therefore the undetermined coefﬁcients ij with i + j = k can be found by expanding (A.25) to order k. We begin with k = 2. To order | |2 , (A.25) is

2 ¯ ¯ 2 ¯ + 11 ¯ + (2¯ − )02 = g20 + g11 ¯ + g02 . 2 2 2 2 2

20

2

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

163

Thus

20 =

g20 ,

11 =

g11 , ¯

02 =

g02 . ¯ 2 −

(A.27)

Taking this into account, we next ﬁnd that equating coefﬁcients of 2 ¯ on both sides of (A.25) implies that c1 () =

g21 g20 g11 (2 + ¯ ) |g11 |2 |g02 |2 + . + + ¯ 2||2 2 2(2 − )

Thus, for = 0, we obtain g21 i 1 2 2 g20 g11 − 2|g11 | − |g02 | + , c1 (0) = 2 0 3 2

(A.28)

(A.29)

where we have used the hypothesis that (0) = i0 . A.4. Restriction to the center manifold Writing (A.1) as X˙ = A()X + f (X, ) = F (X, ).

(A.30)

Let q() and q ∗ () be eigenvectors for A() = FX (0, ) and AT , respectively, corresponding to the simple eigenvalues

() = () + i() and ¯ () of A(); that is Aq = q,

AT q ∗ = ¯ q ∗ .

(A.31)

We normalize q ∗ relative to q so that q ∗ , q = 1, where ·, · denotes the Hermitian product u, v = For any solution x of (A.30), we deﬁne

(A.32) n

i=1 u¯i vi .

z(t) = q ∗ , x(t).

(A.33)

We shall use z and z¯ (in the directions q and q) ¯ as local coordinates. We also deﬁne ¯ ) = x(t) − 2 Re[z(t)q()]. w(t) = x(t) − z(t)q() − z¯ q(

(A.34)

Because zq and z¯ q¯ will always occur together, our choice of complex coordinates will not introduce complex-valued solutions of (A.30).

164

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

Note that q ∗ , q ¯ = 0. A slightly different way of looking at the decomposition of x(t) into z(t) and w(t) is by means of the projection matrices P = q(q¯ ∗ )T + q(q ¯ ∗ )T = 2 Re[q(q¯ ∗ )T ], and P⊥ = I − P = I − 2 Re[q(q¯ ∗ )T ]. These obey P2 = P ,

P⊥2 = P⊥ ,

P P⊥ = 0,

P⊥ P = 0.

In terms of P and P⊥ z(t)q + z¯ (t)q¯ = P x(t) and w(t) = P⊥ x(t). In the variables z and w, (A.30) becomes z˙ = ()z + G(z, z¯ , w, ), w˙ = A()w + H (z, z¯ , w, ),

(A.35)

where G(z, z¯ , w, ) = q ∗ , f (w + 2 Re[zq], ) H (z, z¯ , w, ) = f (w + 2 Re[zq], ) − 2 Re[qG].

(A.36)

Re q ∗ , w = 0

(A.37)

Since and

Im q ∗ , w = 0,

q ∗ , w = 0.

References [1] A. Buonomo, The periodic solution of van der Pol’s equation, SIAM J. Appl. Math. 59 (1998) 156–171. [2] J. Cao, T. Chen, Globally exponentially robust stability and periodicity of delayed neural networks, Chaos Solitons Fract. 22 (4) (2004) 957–963. [3] J. Cao, Daniel W.C. Ho, A general framework for global asymptotic stability analysis of delayed neural networks based on LMI approach, Chaos Solitons Fract. 24 (5) (2005) 1317–1329. [4] J. Cao, D.-S. Huang, Y. Qu, Global robust stability of delayed recurrent neural networks, Chaos Solitons Fract. 23 (1) (2005) 221–229. [5] J. Cao, J. Liang, Boundedness and stability for Cohen–Grossberg neural networks with time-varying delays, J. Math. Anal. Appl. 296 (2) (2004) 665–685. [6] J. Cao, J. Liang, J. Lam, Exponential stability of high-order bidirectional associative memory neural networks with time delays, Physica D 199 (3–4) (2004) 425–436. [7] J. Cao, J. Wang, Absolute exponential stability of recurrent neural networks with time delays and Lipschitzcontinuous activation functions, Neural Networks 17 (3) (2004) 379–390.

W. Yu, J. Cao / Nonlinear Analysis 62 (2005) 141 – 165

165

[8] J. Cao, J. Wang, Global asymptotic and robust stability of recurrent neural networks with time delays, IEEE Trans. Circuits and Systems—I 52 (2) (2005) 417–426. [9] J.E. Flaherty, F.C. Hoppensteadt, Studies in Applied Mathematics, vol. 5, Blackwell, Cambridge, MA, 1978, pp. 58–72. [10] D.R. Frey, A class of relaxation algorithms for ﬁnding the periodic steady-state solution, IEEE Trans. Circuits System—I 45 (1998) 659–663. [11] K. Gopalsamy, I. Leung, Convergence under dynamical thresholds with delays, IEEE Trans. Neural Networks 8 (1997) 341–348. [12] J. Guckenheimer, P.J. Holmes, Nonlinear Oscillations, Dynamical System and Bifurcations of Vector Fields, Springer, Berlin, 1983. [13] J. Hale, S. Verduyn Lunel, Introduction to Functional Differential Equations, Springer, New York, 1993. [14] B.D. Hassard, N.D. Kazarinoff,Y.H. Wan, Theory and Application of Hopf Bifurcation, Cambridge University Press, Cambridge, 1981. [15] G. Iooss, D.D. Joiseph, Elementary Stability and Bifurcation Theory, second ed., Springer, New York, 1989. [16] Y. Kuang, Delay Differential Equations with Applications in Population Dynamics, Academic Press, New York, 1993. [17] X. Liao, K. Wong, Z. Wu, Hopf bifurcation and stability of periodic solutions for van der Pol equation with distributed delay, Nonlinear Dyn. 26 (2001) 23–44. [18] X. Liao, Z. Wu, J. Yu, Stability switches and bifurcation analysis of a neural network with continuously delay, IEEE Trans. Systems Man Cybernet.—Part A 29 (1999) 692–696. [19] N. MacDonald, Biological Delay Systems: Linear Stability Theory, Cambridge University Press, Cambridge, 1989. [20] C.M. Marcus, R.M. Westervelt, Stability of analog neural network with delay, Phys. Rev. A 39 (1989) 347–359. [21] J.E. Marsden, M. McCracken, The Hopf Bifurcation and its Applications, Applied Mathematics Sciences, vol. 19, Springer, Heidelberg, 1976. [22] K. Murakaimi, Bifurcated periodic solutions for delayed van der Pol equation, Neural Parallel Sci. Comput. 7 (1999) 1–16. [23] H.G. Schuster, Deterministic Chaos, Physik-Verlag, Weinheim, 1984. [24] V. Venkatasubramanian, Singularity induced bifurcation and the van der Pol oscillator, IEEE Trans. Circuits System—I 41 (1994) 765–769. [25] J. Wu, Theory and Applications of Partial Functional Differential Equations, Springer, New York, 1996.

This article was originally published in a journal published by Elsevier, and the attached copy is provided by Elsevier for the author’s benefit and for the benefit of the author’s institution, for non-commercial research and educational use including without limitation use in instruction at your institution, sending it to specific colleagues that you know, and providing a copy to your institution’s administrator. All other uses, reproduction and distribution, including without limitation commercial reprints, selling or licensing copies or access, or posting on open internet sites, your personal or institution’s website or repository, are prohibited. For exceptions, permission may be sought for such use through Elsevier’s permissions site at: http://www.elsevier.com/locate/permissionusematerial

Applied Mathematical Modelling 31 (2007) 2111–2122 www.elsevier.com/locate/apm

Wangli He, Jinde Cao

*

co

py

Stability and bifurcation of a class of discrete-time neural networks Department of Mathematics, Southeast University, Nanjing 210096, China

al

Received 1 May 2001; received in revised form 1 June 2006; accepted 30 August 2006 Available online 24 October 2006

on

Abstract

pe

rs

In this paper, a class of discrete-time system modelling a network with two neurons is considered. Its linear stability is investigated and Neimark–Sacker bifurcation (also called Hopf bifurcation for map) is demonstrated by analyzing the corresponding characteristic equation. In particular, the explicit formula for determining the direction of Neimark–Sacker bifurcation and the stability of periodic solution is obtained by using the normal form method and the center manifold theory for discrete time system developed by Kuznetsov. The theoretical analysis is veriﬁed by numerical simulations. Ó 2006 Elsevier Inc. All rights reserved. Keywords: Neural network; Stability; Neimark–Sacker bifurcation; Discrete-time system

r's

1. Introduction

th o

In the past few decades, neural networks have received intensive interest due to their wide applications, such as, pattern recognition, associative memory and combinational optimization, and its dynamical behavior plays an important role. Many works [1–15] have been published to investigate the dynamics of neural networks since Hopﬁeld [16] constructed a simpliﬁed neural network model. Neural networks with one or two neurons are prototypes to understand the dynamics of large-scale networks, many results have been made for such simpliﬁed networks [17,3–7,9,10,13,14,18]. In 2003, Yuan and Huang [19] studied the asymptotical behavior of the following diﬀerence system:

Au

x1 ðn þ 1Þ ¼ bx1 ðnÞ þ a11 f ðx1 ðnÞÞ þ a12 f ðx2 ðnÞÞ; x2 ðn þ 1Þ ¼ bx2 ðnÞ þ a21 f ðx1 ðnÞÞ þ a22 f ðx2 ðnÞÞ;

n ¼ 0; 1; 2; . . . ;

where b 2 (0, 1) is a constant and f : R ! R is the activation function given by the piecewise constant McCulloch–Pitts nonlinearity *

Corresponding author. Tel.: +86 25 3792315; fax: +86 25 3792316. E-mail addresses: [email protected], [email protected] (J. Cao).

0307-904X/$ - see front matter Ó 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2006.08.006

2112

W. He, J. Cao / Applied Mathematical Modelling 31 (2007) 2111–2122

f ðuÞ ¼

1;

u > r;

þ1;

u 6 r;

where r 2 R is a constant, and acts as the threshold. In 2004, Yuan et al. [9] considered the following system: x1 ðn þ 1Þ ¼ bx1 ðnÞ þ ð1 bÞf ðax1 ðnÞÞ þ ð1 bÞf ðc1 x2 ðnÞÞ; x2 ðn þ 1Þ ¼ bx2 ðnÞ ð1 bÞf ðc2 x1 ðnÞÞ þ ð1 bÞf ðax2 ðnÞÞ;

n ¼ 0; 1; 2; . . . ;

co

py

where b 2 (0, 1) is internal decay of the neurons, the constant a > 0 and ci (i = 1, 2) denote the gain parameters, f : R ! R is a continuous transfer function and f(0) = 0. They discussed the global and local stability of the equilibrium, which gave some suﬃcient conditions to guarantee the existence of bifurcation, meanwhile got a formula to determine the direction and stability of bifurcation. In 2005, Yuan et al. [10] introduced a more general model based on [9]. They studied the following discrete-time neural network model with self-connection in the following form x1 ðn þ 1Þ ¼ bx1 ðnÞ þ a11 f1 ðx1 ðnÞÞ þ a12 f2 ðx2 ðnÞÞ; x2 ðn þ 1Þ ¼ bx2 ðnÞ þ a21 f1 ðx1 ðnÞÞ þ a22 f2 ðx2 ðnÞÞ;

n ¼ 0; 1; 2; . . . ;

ð1:0Þ

rs

on

al

where xi (i = 1, 2) denotes the state of the ith neuron, b 2 (0, 1) is internal decay of the neurons, the constants aij (i, j = 1, 2) denotes the connection weights, fi : R ! R (i = 1, 2) are continuous transfer functions and fi(0) = 0 (i = 1, 2). Some suﬃcient conditions were given to guarantee the stability of the equilibrium and the existence of bifurcation. The direction of bifurcation and the stability of bifurcating periodic solution were discussed. However, in a real neural network every neuron is usually independent and has its characteristics. In most cases, the internal decay of neurons is diﬀerent, so it is necessary to study the dynamical behavior of neural networks with diﬀerent internal decay of the neurons. In this paper, we consider the following more general model: x1 ðn þ 1Þ ¼ ax1 ðnÞ þ a11 f1 ðx1 ðnÞÞ þ a12 f2 ðx2 ðnÞÞ;

pe

x2 ðn þ 1Þ ¼ bx2 ðnÞ þ a21 f1 ðx1 ðnÞÞ þ a22 f2 ðx2 ðnÞÞ;

n ¼ 0; 1; 2; . . . ;

ð1:1Þ

where a 2 (0, 1), b 2 (0, 1) are internal decay of neurons, aij (i = 1, 2) denote the connection weights, fi : R ! R (i = 1, 2) are continuous transfer functions and fi(0) = 0 (i = 1, 2). The discrete-time system (1.1) can be regarded as a discrete form of the diﬀerential system x_ 1 ðtÞ ¼ l1 x1 ðtÞ þ c11 f1 ðx1 ðtÞÞ þ c12 f2 ðx2 ðtÞÞ;

ð1:2Þ

r's

x_ 2 ðtÞ ¼ l2 x1 ðtÞ þ c21 f1 ðx1 ðtÞÞ þ c22 f2 ðx2 ðtÞÞ; or the system with a piecewise constant arguments

th o

x_ 1 ðtÞ ¼ l1 x1 ðtÞ þ c11 f1 ðx1 ð½tÞÞ þ c12 f2 ðx2 ð½tÞÞ; x_ 2 ðtÞ ¼ l2 x1 ðtÞ þ c21 f1 ðx1 ð½tÞÞ þ c22 f2 ðx2 ð½tÞÞ;

ð1:3Þ

Au

where l1 > 0, l2 > 0 and [Æ] denotes the greatest integer function. This system has wide application in certain biomedical areas and much progress has been made in the study of such system (1.3) with the piecewise arguments [20]. For the method of discrete analogy, we refer to [7,21,19,22]. Obviously, system (1.1) includes the discrete version of systems (1.2) and (1.3), which is an improved model, compared with the model in [10]. The discussion of the stability of equilibrium is too complex and we study it by setting a proper parameter. In order to analyze bifurcation, we ﬁnd a simpler bifurcation parameter and derive more general results. Bifurcation analysis concerning the continuous dynamical systems was investigated in [13–15]. However, in this paper, based on the techniques developed by Kuznetsov [11], we not only investigate the stability of equilibrium and the existence of Neimark–Sacker bifurcation of system (1.1), but also the direction of the Neimark–Sacker bifurcation and the stability of the bifurcating periodic solution. We ﬁnd that when the bifurcation parameter exceeds a critical value, the Neimark–Sacker bifurcation will occur and its direction and stability are determined completely by an algebraic expression of a(D*).

W. He, J. Cao / Applied Mathematical Modelling 31 (2007) 2111–2122

2113

The organization of this paper is as follows: In Section 2, we will discuss the stability of the trivial solutions and the existence of Neimark–Sacker bifurcation. In Section 3, a formula for determining the direction of Neimark–Sacker bifurcation and the stability of bifurcating periodic solution will be given by using the normal form method and the center manifold theory for discrete time system developed by Kuznetsov. In Section 4, numerical simulations aimed at justifying the theoretical analysis will be reported.

py

2. Stability and existence of Neimark–Sacker bifurcation In this section, we discuss the local stability of the equilibrium (0, 0) of system (1.1). In [5,7], the transfer function f of the models which they discussed is f(u) = tanh(cu), here we only need the following assumption: and

f i ð0Þ ¼ 0;

i ¼ 1; 2:

co

ðH1 Þ f i 2 C 1 ðRÞ

For the sake of simplicity and the need of discussion, the following parameters are deﬁned: D ¼ a12 a21 f10 ð0Þf20 ð0Þ;

al

1 1 T 1 ¼ ða þ a11 f10 ð0ÞÞ; T 2 ¼ ðb þ a22 f20 ð0ÞÞ; 2 2 1 T ¼ ða11 f10 ð0Þ þ a22 f20 ð0ÞÞ: 2

rs

on

Theorem 1. The zero solution of (1.1) is asymptotically stable if (H1) is satisfied and (T1, T2, D) 2 X0, where X0 = X1 [ X2 [ X3, aþb aþb 2 3 1 þ 2ðT 1 þ T 2 Þ 4T 1 T 2 ; 2 2 aþb aþb 2 3 1 2ðT 1 þ T 2 Þ 4T 1 T 2 ; 1 2 2

pe

X 3 ¼ fðT 1 ; T 2 ; DÞ 2 R3 ; D < 1 4T 1 T 2 ; ðT 1 T 2 Þ2 < Dg:

r's

Proof. The associated characteristic matrix for the linearization of (1.1) at (0, 0) is a þ a11 f10 ð0Þ k a12 f20 ð0Þ; A¼ : a21 f10 ð0Þ b þ a22 f20 ð0Þ k Let jAj = 0, we get the associated characteristic equation of (1.1). That is k2 ða þ b þ a11 f10 ð0Þ þ a22 f20 ð0ÞÞk þ ða þ a11 f10 ð0ÞÞðb þ a22 f20 ð0ÞÞ a12 a21 f10 ð0Þf20 ð0Þ ¼ 0;

ð2:1Þ

2

th o

D ¼ ða þ b þ a11 f10 ð0Þ þ a22 f20 ð0ÞÞ 4ða þ a11 f10 ð0ÞÞðb þ a22 f20 ð0ÞÞ þ 4a12 a21 f10 ð0Þf20 ð0Þ

¼ a2 2ab þ b2 þ 2aa11 f10 ð0Þ 2aa22 f20 ð0Þ 2ba11 f10 ð0Þ þ 2ba22 f20 ð0Þ

Au

þ a211 f10 ð0Þ2 2a11 a22 f10 ð0Þf20 ð0Þ þ a222 f20 ð0Þ2 þ 4a12 a21 f10 ð0Þf20 ð0Þ 2

¼ ða b þ a11 f10 ð0Þ a22 f20 ð0ÞÞ þ 4a12 a21 f10 ð0Þf20 ð0Þ: Next we discuss it in two cases: Case 1. (T1 T2)2 > D. In this case, the roots of characteristic equation (2.1) are given by qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 k1;2 ¼ T 1 þ T 2 ðT 1 T 2 Þ D:

ð2:2Þ

Note that the eigenvalues jk1,2j < 1 if and only if ðT 1 ; T 2 ; DÞ 2 X 1 [ X 2 ;

ð2:3Þ

2114

W. He, J. Cao / Applied Mathematical Modelling 31 (2007) 2111–2122

where

aþb aþb 2 1 þ 2ðT 1 þ T 2 Þ 4T 1 T 2 ; 2 2 aþb aþb 2 1 2ðT 1 þ T 2 Þ 4T 1 T 2 ; 1 2 2 3

The modulus of the eigenvalues jk1,2j < 1 if and only if 2

co

py

Case 2. (T1 T2)2 < D. In this case the characteristic equation (2.1) has a pair of conjugate complex roots qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 k1;2 ¼ T 1 þ T 2 D ðT 1 T 2 Þ i: ð2:4Þ

ðT 1 ; T 2 ; DÞ 2 X 3 ¼ fðT 1 ; T 2 ; DÞ 2 R3 ; D < 1 4T 1 T 2 ; ðT 1 T 2 Þ < Dg:

ð2:5Þ

on

al

Combining with Cases 1 and 2, we get X0 = X1 [ X2 [ X3. Thus the eigenvalues k1,2 of characteristic equation (2.1) are inside the unit circle for (T1, T2, D) 2 X0, and thus the zero solution of (1.1) is asymptotically stable. Now, we choose D as the bifurcation parameter to study the Neimark–Sacker bifurcation of (0, 0) in this paper. In case of (T1 T2)2 < D, the characteristic equation (2.1) has a pair of conjugate complex roots, let qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 ð2:6Þ kðDÞ ¼ T 1 þ T 2 þ D ðT 1 T 2 Þ i;

rs

then the eigenvalues in (2.1) are k(D) and kðDÞ. The modulus of the eigenvalue is qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 jkj ¼ ðT 1 þ T 2 Þ þ D ðT 1 T 2 Þ ¼ D þ 4T 1 T 2 : D ¼ D ¼ 1 4T 1 T 2 : Obviously, 2

jkj < 1 for ðT 1 T 2 Þ < D < D :

pe

jkj = 1 if and only if

ð2:7Þ

ð2:8Þ

h

r's

Then we can conclude that D* is a critical value which destroy the stability of (0, 0).

th o

Lemma 1. Suppose that (H1) is satisfied and aþb < T < 1 aþb . Then 2 2 d jkðDÞj (i) > 0; dD D¼D (ii) kk(D*) 5 1, k = 1, 2, 3, 4, where k(D) and D* are given by (2.6) and (2.8), respectively.

Au

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Proof. We know jkðDÞj ¼ D þ 4T 1 T 2 , By direct calculation, we obtain that d 1 jkðDÞj ¼ > 0; dD 2 D¼D which means property (i) is true. Next we deal with the property of (ii). From aþb < T < 1 aþb , we know 2 2 k 0 < T1 + T2 < 1. Clearly, k (D*) = 1 for some k 2 {1, 2, 3, 4} if and only if the argument k(D*) 2 {0, ±p/2, ±2p/3,p}. Since qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ jkðD Þj ¼ 1; RekðD Þ ¼ T 1 þ T 2 > 0; ImkðD Þ ¼ 1 ðT 1 þ T 2 Þ2 > 0; it follows that arg kðD Þ 2 f0; p=2; 2p=3; pg. The proof is completed.

h

W. He, J. Cao / Applied Mathematical Modelling 31 (2007) 2111–2122

2115

By the Hopf bifurcation for maps [12], we have the following theorem. aþb Theorem 2. Let the assumption (H1) hold and aþb 2 < T < 1 2 . Then we have

py

(i) if (T1 T2)2 < D < D*, then the equilibrium (0, 0) of (1.1) is asymptotically stable; (ii) if (T1 T2)2 < D and D > D*, then the equilibrium (0, 0) of (1.1) is unstable; (iii) if D = D*, system (1.1) undergoes a Neimark–Sacker bifurcation, that is, system (1.1) has a periodic solution bifurcating from the equilibrium (0, 0) near D = D*, where D* is given by (2.8). 3. Direction and stability of the Neimark–Sacker bifurcation

ðH2 Þ f i 2 C 3 ðR; RÞ;

f i ð0Þ ¼ fi00 ð0Þ ¼ 0;

al

co

In the previous section, we have obtained some suﬃcient conditions to ensure that system (1.1) undergoes a Neimark–Sacker bifurcation at the equilibrium (0, 0) when D takes the critical value D*. In this section, a formula for determining the direction of Neimark–Sacker bifurcation and the stability of bifurcating periodic solutions of system (1.1) at D = D* shall be presented by employing the normal form method and the center manifold theory for discrete time system developed by Kuznetsov [11]. For most of the models in the literature [5,7,8], the transfer function f is f(u) = tanh(cu) , here we assume that the transfer functions in (1.1) satisfy fi0 ð0Þfi000 ð0Þ 6¼ 0:

ð3:1Þ

rs

on

Now system (1.1) can be rewritten as a12 f20 ð0Þ a þ a11 f10 ð0Þ x1 ðn þ 1Þ F 1 ðx; DÞ x1 ðnÞ ¼ þ ; a21 f10 ð0Þ b þ a22 f20 ð0Þ x2 ðn þ 1Þ x2 ðnÞ F 2 ðx; DÞ

pe

where x = (x1, x2)T 2 R2. From assumption (H2), we know that Fi(i = 1, 2) in (3.1) can be expanded as a11 000 a12 000 f1 ð0Þx31 þ f ð0Þx32 þ Oðkxk4 Þ; F 1 ðx; DÞ ¼ 6 6 2 a21 000 a22 000 f1 ð0Þx31 þ f ð0Þx32 þ Oðkxk4 Þ: F 2 ðx; DÞ ¼ 6 6 2 Denote B ¼ BðDÞ ¼

a þ a11 f10 ð0Þ

a12 f20 ð0Þ

a21 f10 ð0Þ

b þ a22 f20 ð0Þ

r's

:

ð3:2Þ

th o

For convenience of calculation, we introduce two notations: qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 r1 ¼ T 2 T 1 þ D ðT 1 T 2 Þ i; qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 r2 ¼ T 1 T 2 þ D ðT 1 T 2 Þ i:

ð3:3Þ ð3:4Þ

From the deﬁnition of T1, T2, we can obtain jrj j2 ¼ r1 r2 ¼ D ¼ a12 a21 f10 ð0Þf20 ð0Þ;

Au

r1 ¼ r2 ;

j ¼ 1; 2:

ð3:5Þ

Here for j = 1, 2, rj 5 0 and a12 a21 f10 ð0Þf20 ð0Þ < 0. In fact, if r1 = 0 or r2 = 0, From (3.5), we have r1r2 = D = 0, then (T1 T2)2 P D, and this is a contradiction to (T1 T2)2 < D. Hence rj 5 0(j = 1, 2) and it follows from (3.5) that a12 a21 f10 ð0Þf20 ð0Þ < 0. Let q(D) 2 C2 be the eigenvector of B(D) corresponding to k(D) given by (2.6). Then BðDÞqðDÞ ¼ kðDÞqðDÞ: Also let p(D) 2 C2 be the eigenvector of the transposed matrix BT(D) corresponding to its eigenvalue kðDÞ, BT ðDÞpðDÞ ¼ kðDÞpðDÞ:

2116

W. He, J. Cao / Applied Mathematical Modelling 31 (2007) 2111–2122

By direct computation we obtain that q

a21 f10 ð0Þ 1; r2

T

p

;

a12 f20 ð0Þ 1; r2

T ; T

where rj(j = 1, 2) are given by (3.3), (3.4). In order to normalize the q ¼ ð1; a21 f10 ð0Þ=r2 Þ and p, we obtain 1;

a12 f20 ð0Þ r2

T

py

r2 p¼ r2 r2

:

co

It is easy to see hp, qi = 1, where hÆ, Æi means the standard scalar product in C 2 : hp; qi ¼ p1 q1 þ p2 q2 . Any vector x 2 R2 can be represented for D near D* as x ¼ zqðDÞ þ zqðDÞ

al

for some complex z, obviously z ¼ hpðDÞ; xi:

on

Thus, system (3.1) can be transformed for D near D* into the following form: z_ ¼ kðDÞz þ gðz; z; DÞ;

ð3:6Þ

X kþlP2

1 g ðDÞzkzl : k!l! kl

pe

gðz; z; DÞ ¼

rs

where k(D) can be written as k(D) = (1 + u(D))eih(D)(u(D) is a smooth function with u(D*) = 0) and ð3:7Þ

We know that Fi(i = 1, 2) in (3.1) can be expanded as a11 000 a12 000 4 f1 ð0Þn31 þ f ð0Þn32 þ Oðknk Þ; 6 6 2 a21 000 a22 000 4 f1 ð0Þn31 þ f2 ð0Þn32 þ Oðknk Þ: F 2 ðn; DÞ ¼ 6 6

It follows that Bi ðx; yÞ :¼

th o

r's

F 1 ðn; DÞ ¼

2 X o2 F i ðn; D Þ on on j;k

j

k

xj xk ¼ 0;

Au

2 X o3 F i ðn; D Þ C i ðx; y; uÞ :¼ on on on j;k;l

j

k

i ¼ 1; 2;

ð3:8Þ

n¼0

l

xj xk ul n¼0

¼ ai1 f1000 ð0Þx1 y 1 u1 þ ai2 f2000 ð0Þx2 y 2 u2 ;

i ¼ 1; 2:

By (3.7)–(3.9) and the formulas g20 ðD Þ ¼ hp; Bðq; qÞi; g21 ðD Þ ¼ hp; Cðq; q; qÞi:

g11 ðD Þ ¼ hp; Bðq; qÞi;

g02 ðD Þ ¼ hp; Bðq; qÞi;

ð3:9Þ

W. He, J. Cao / Applied Mathematical Modelling 31 (2007) 2111–2122

We obtain g20 ðD Þ ¼ g11 ðD Þ ¼ g02 ðD Þ ¼ 0; qÞ þ p2 C 2 ðq; q; qÞ g21 ðD Þ ¼ p1 C 1 ðq; q; ( ) 2 2 2 0 r2 a21 f1 ð0Þ f2000 ð0Þ a12 a21 f20 ð0Þf1000 ð0Þ a221 a22 f10 ð0Þ f2000 ð0Þ 000 ¼ þ a11 f1 ð0Þ r2 f20 ð0Þ r2 r22 r2 r2 1 fa11 a12 r2 f20 ð0Þf1000 ð0Þ þ a21 a22 r2 f10 ð0Þf2000 ð0Þ a12 ðr1 þ r2 Þf20 ð0Þ 2

py

¼

2

1 ða11 a12 f20 ð0Þf1000 ð0Þ a21 a22 f10 ð0Þf2000 ð0ÞÞ 2a12 f20 ð0Þ 1 qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ þ fða11 a12 f20 ð0Þf1000 ð0Þ þ a21 a22 f10 ð0Þf2000 ð0ÞÞðT 1 T 2 Þ 2 0 2a12 D ðT 1 T 2 Þ f2 ð0Þi 2

on

¼

al

co

þ a212 a21 f20 ð0Þ f1000 ð0Þ a12 a221 f10 ð0Þ f2000 ð0Þg qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 2 0 000 qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ a11 a12 f2 ð0Þf1 ð0Þ T 1 T 2 D ðT 1 T 2 Þ i ¼ 2 0 2a12 D ðT 1 T 2 Þ f2 ð0Þi qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 þ a21 a22 f10 ð0Þf2000 ð0Þ T 1 T 2 D ðT 1 T 2 Þ i 2 000 2 000 2 0 2 0 þ a12 a21 f2 ð0Þ f1 ð0Þ a12 a21 f1 ð0Þ f2 ð0Þ

2

rs

þ a212 a21 f20 ð0Þ f1000 ð0Þ a12 a221 f10 ð0Þ f2000 ð0Þg:

2

r's

pe

Together with eihðD Þ ¼ kðD Þ and D ¼ a12 a21 f10 ð0Þf20 ð0Þ, we have ihðD Þ e g21 ð1 2eihðD Þ Þe2ihðD Þ 1 1 2 2 aðD Þ ¼ Re g20 g11 jg11 j jg02 j Re Þ ihðD 2 4 2 2ð1 e Þ ihðD Þ e g21 ¼ Re 2 1 ¼ ReðkðD Þg21 Þ 2 1 fðT 1 þ T 2 Þ½a11 a12 f20 ð0Þf1000 ð0Þ a21 a22 f10 ð0Þf2000 ð0Þ ¼ 4a12 f20 ð0Þ ða11 a12 f20 ð0Þf1000 ð0Þ þ a21 a22 f10 ð0Þf2000 ð0ÞÞðT 1 T 2 Þ 2

th o

a212 a21 f20 ð0Þ f1000 ð0Þ þ a12 a221 f10 ð0Þ f2000 ð0Þg 1 ¼ f2T 2 a11 a12 f20 ð0Þf1000 ð0Þ 2T 1 a21 a22 f10 ð0Þf2000 ð0Þ 4a12 f20 ð0Þ 2

2

Au

a212 a21 f20 ð0Þ f1000 ð0Þ þ a12 a221 f10 ð0Þ f2000 ð0Þg 1 fðb þ a22 f20 ð0ÞÞa11 a12 f20 ð0Þf1000 ð0Þ ¼ 4a12 f20 ð0Þ 2

ða þ a11 f10 ð0ÞÞa21 a22 f10 ð0Þf2000 ð0Þ a212 a21 f20 ð0Þ f1000 ð0Þ 2

þ a12 a221 f10 ð0Þ f2000 ð0Þg 1 ðb þ a22 f20 ð0ÞÞa11 a12 f20 ð0Þf1000 ð0Þ ¼ 4a12 f20 ð0Þ ða þ a11 f10 ð0ÞÞa21 a22 f10 ð0Þf2000 ð0Þ 1 4T 1 T 2 2 000 2 000 0 0 ða12 f2 ð0Þ f1 ð0Þ a21 f1 ð0Þ f2 ð0ÞÞ þ 0 f1 ð0Þf20 ð0Þ

2117

2118

W. He, J. Cao / Applied Mathematical Modelling 31 (2007) 2111–2122

¼

1 fðb þ a22 f20 ð0ÞÞa11 a12 f10 ð0Þf20 ð0Þ2 f1000 ð0Þ 0 4a12 f1 ð0Þf20 ð0Þ2 2 ða þ a11 f10 ð0ÞÞa21 a22 f10 ð0Þ f20 ð0Þf2000 ð0Þ þ ½1 ða þ a11 f10 ð0ÞÞðb þ a22 f20 ð0ÞÞða12 f20 ð0Þ2 f1000 ð0Þ

a21 f10 ð0Þ2 f2000 ð0ÞÞg

1 fba11 a12 f10 ð0Þf20 ð0Þ2 f1000 ð0Þ þ a11 a12 a22 f10 ð0Þf20 ð0Þ3 f1000 ð0Þ 0 4a12 f1 ð0Þf20 ð0Þ2 2 3 aa21 a22 f10 ð0Þ f20 ð0Þf2000 ð0Þ a11 a21 a22 f10 ð0Þ f20 ð0Þf2000 ð0Þ þ a12 f20 ð0Þ2 f1000 ð0Þð1 abÞ a21 f10 ð0Þ2 f2000 ð0Þð1 abÞ aa12 a22 f20 ð0Þ3 f1000 ð0Þ þ aa21 a22 f10 ð0Þ2 f20 ð0Þf2000 ð0Þ 2 3 ba11 a12 f10 ð0Þf20 ð0Þ f1000 ð0Þ þ ba11 a21 f10 ð0Þ f2000 ð0Þ 3 3 a11 a12 a22 f10 ð0Þf20 ð0Þ f1000 ð0Þ þ a11 a21 a22 f10 ð0Þ f20 ð0Þf2000 ð0Þg f1000 ð0Þ a21 f10 ð0Þf2000 ð0Þ 0 1 þ ab þ ba f ð0Þ þ 0 ½1 ab aa22 f20 ð0Þ: ¼ 11 1 4a12 f20 ð0Þf20 ð0Þ 4f1 ð0Þ

co

py

¼

al

ð3:10Þ

rs

on

aþb Theorem 3. Suppose that the condition (H2) hold and T 2 ð aþb 2 ; 1 2 Þ, then the direction of the Neimark– Sacker bifurcation and stability of bifurcating periodic solution can be determined by the sign of a(D*). In fact, if a(D*) < 0(>0), then the Neimark–Sacker bifurcation is supercritical (subcritical) and the bifurcating periodic solution is asymptotically stable (unstable), where D* is given by (2.8). This method is introduced by Kuznetsov in [11].

x1 ðn þ 1Þ ¼ ax1 ðnÞ þ a12 f2 ðx2 ðnÞÞ; x2 ðn þ 1Þ ¼ bx1 ðnÞ þ a21 f1 ðx2 ðnÞÞ;

pe

Remark 1. If the two-neuron network (1.1) is without self-connections, that is, a11 = a22 = 0, then system (1.1) has the following form: n ¼ 0; 1; 2; . . .

r's

From (3.10), we can obtain 1 ab f1000 ð0Þ a21 f10 ð0Þf2000 ð0Þ aðD Þ ¼ : 4 f10 ð0Þ a12 f20 ð0Þf20 ð0Þ

ð3:11Þ

ð3:12Þ

th o

From the previous analysis, we know a12 a21 f10 ð0Þf20 ð0Þ < 0. If sgnðf10 ð0Þf1000 ð0ÞÞ ¼ sgnðf20 ð0Þf2000 ð0ÞÞ, we have sgnðaðD ÞÞ ¼ sgnðf10 ð0Þf1000 ð0ÞÞ ¼ sgnðf20 ð0Þf2000 ð0ÞÞ from (3.12). Thus we obtain the following result.

Au

Corollary 1. Suppose that (H2) is satisfied and sgnðf10 ð0Þf1000 ð0ÞÞ ¼ sgnðf20 ð0Þf2000 ð0ÞÞ. Then the direction of the Neimark–Sacker bifurcation and stability of bifurcating periodic solution are determined by the sign of fk0 ð0Þfk000 ð0Þ. Indeed, if fk0 ð0Þfk000 ð0Þ < 0ð> 0Þ, the Neimark–Sacker bifurcation is supercritical (subcritical) and the bifurcating periodic solution is asymptotically stable (unstable). Remark 2. If the two-neuron network with self-connections modelled by a discrete-time system of the form (1.1) with a = b x1 ðn þ 1Þ ¼ bx1 ðnÞ þ a11 f1 ðx1 ðnÞÞ þ a12 f2 ðx2 ðnÞÞ; x2 ðn þ 1Þ ¼ bx1 ðnÞ þ a21 f1 ðx1 ðnÞÞ þ a22 f2 ðx2 ðnÞÞ;

ð3:13Þ

which is the model studied by Yuan et al. [10]. From (3.10), we obtain aðD Þ ¼

a21 f10 ð0Þf2000 ð0Þ f1000 ð0Þ 2 0 ½1 þ b ½1 b2 ba22 f20 ð0Þ: þ ba f ð0Þ þ 11 1 0 0 0 4a12 f2 ð0Þf2 ð0Þ 4f1 ð0Þ

ð3:14Þ

W. He, J. Cao / Applied Mathematical Modelling 31 (2007) 2111–2122

2119

It tallies with the result of Yuan et al. [10]. Here we ﬁnd the two bifurcation parameters are diﬀerent. In our paper, we assume D ¼ a12 a21 f10 ð0Þf20 ð0Þ as the bifurcation parameter. In [10] the bifurcation parameter is D0 ¼ ða11 a22 a12 a21 Þf10 ð0Þf20 ð0Þ. We obtain D0 ¼ D þ a11 a22 f10 ð0Þf20 ð0Þ. In fact it does not inﬂuence the sign of a(D*). The bifurcation parameter in our paper is simpler. System (1.1) includes the model of Yuan et al. [10]. Yuan et al. [10] only discussed a two-neuron model with the same internal decay of the neurons, however our model can deal with the system with different internal decay of the neurons, which is more general.

py

4. Numerical simulations

In this section, we give numerical simulations to support our theoretical analysis.

1 f20 ð0Þ ¼ ; 2

co

Example 1. Let a ¼ 14 ; b ¼ 34 ; a11 ¼ 1; a12 ¼ 1; a22 ¼ 1 and f1(u) = sin(u), f2(u) = arctan(u/2) in the system (1.1). By the simple calculation, we obtain 1 f2000 ð0Þ ¼ ; 4 1 1 3 T ¼ ða11 f10 ð0Þ þ a22 f20 ð0ÞÞ ¼ ; 1 þ 2ðT 1 þ T 2 Þ 4T 1 T 2 ¼ ; 2 4 16 45 11 1 2ðT 1 þ T 2 Þ 4T 1 T 2 ¼ ; 1 4T 1 T 2 ¼ : 16 16

11 11 It follows from (2.8) that D ¼ 16 a21 ¼ 8 , the Neimark–Sacker bifurcation occurs when D ¼ 11 . Choosing 16 a21 ¼ 54, so that D ¼ 58 < D ; ðT 1 ; T 2 ; DÞ 2 X 0 . By Theorem 1, we know the origin is asymptotically stable, the corresponding waveform and phase plots are shown in Fig. 1. Choosing a21 ¼ 32, then D ¼ 34 > D . By Theorem 2, we know that a Neimark–Sacker bifurcation occurs when D ¼ D ¼ 11 . The origin loses its stability 16 and a periodic solution bifurcates from the origin when D ¼ 34 > D . From (3.10), we have aðD Þ ¼ 131 < 0. By Theorem 3, we know the periodic solution is asymptotically stable. The corresponding 512 phase plots are shown in Figs. 2 and 3. In Fig. 2 its zero solution undergoes a Neimark–Sacker bifurcation at the origin and in Fig. 3 it is stable. f100 ð0Þ ¼ f200 ð0Þ ¼ 0;

f1000 ð0Þ ¼ 1;

pe

rs

on

al

f10 ð0Þ ¼ 1;

Example 2. Let a ¼ 13 ; b ¼ 12 ; a11 ¼ 1; a12 ¼ 1; a22 ¼ 1 and f1 ðuÞ ¼ sinðuÞ; f2 ðuÞ ¼ eu tem (1.1). By the simple calculation, we obtain

1 2

u2 þ 1 in the sys-

th o

r's

f10 ð0Þ ¼ 1; f20 ð0Þ ¼ 1; f100 ð0Þ ¼ f200 ð0Þ ¼ 0; f1000 ð0Þ ¼ 1; f2000 ð0Þ ¼ 1; 1 1 T ¼ ða11 f10 ð0Þ þ a22 f20 ð0ÞÞ ¼ 0; 1 þ 2ðT 1 þ T 2 Þ 4T 1 T 2 ¼ ; 2 2 7 5 1 2ðT 1 þ T 2 Þ 4T 1 T 2 ¼ ; 1 4T 1 T 2 ¼ : 6 3 It follows from (2.8) that D ¼ 53 ða21 ¼ 53Þ, the Neimark–Sacker bifurcation occurs when D ¼ 53. Choosing a21 = 1.6, so that D = 1.6 < D*, (T1, T2, D) 2 X0. By Theorem 1, we know the origin is asymptotically stable, 0.05

0.03

Au

0.04

0.02

0.03

0.01

0.02 0.01

0

0

–0.01

–0.01 –0.02

–0.02

–0.03 –0.03 –0.04 –0.05 –0.04

–0.03

–0.02

–0.01

0

0.01

0.02

0.03

–0.04 0

100

200

300

400

500

Fig. 1. Phase plot and waveform plot for system (1.1) with a21 ¼

5 ; 4

600

700

5 8

800

D¼ D .

0.06

0.04

0.04

0.02

0.02

0

–0.02

–0.02

–0.04

–0.02

–0.01

0

0.01

0.02

–0.04

0.03

0.04

–0.06

0

100

200

300

400

500

th o

–0.03

r's

0

–0.06 –0.04

0.6

3 4

pe

Fig. 3. a21 ¼

3 ; 2

600

700

800

Fig. 4. Phase plot and waveform plot for system (1.1) with a21 = 1.6, D = 1.6 < D*. 0.06

Au

0.04

0.02

0

–0.02

–0.04

–0.06 –0.05

–0.04

–0.03

–0.02

–0.01

Fig. 5. a21 ¼

0

5 ; 3

0.01

0.02

5 3

0.03

D¼ ¼D .

0.04

0.05

900

1000

W. He, J. Cao / Applied Mathematical Modelling 31 (2007) 2111–2122

2121

0.8 0.6 0.4 0.2 0 –0.2

–0.6 –0.8 –0.8

–0.6

–0.4

–0.2

0

0.2

0.4

0.6

0.8

co

Fig. 6. a21 = 1.7, D = 1.7 > D*.

py

–0.4

on

al

the corresponding waveform and phase plots are shown in Fig. 4. Choosing a21 = 1.7, then D = 1.7. By Theorem 2, we know that a Neimark–Sacker bifurcation occurs when D ¼ D ¼ 53. The origin loses its stability and a periodic solution bifurcates from the origin when D = 1.7 > D*. From (3.10), we have aðD Þ ¼ 18 < 0. By Theorem 3, we know the periodic solution is asymptotically stable. The corresponding phase plots are shown in Figs. 5 and 6. In Fig. 5 its zero solution undergoes a Neimark–Sacker bifurcation at the origin and in Fig. 6 it is stable. 5. Conclusions

pe

rs

The discrete-time system of neural networks provides some dynamical behaviors which enriches the theory of continuous system and has potential applications in neural networks. Although the system discussed above is quite simple, it is potentially useful as the complexity found might be carried over to the model with delay. By choosing a proper bifurcation parameter, we have shown that a Neimark–Sacker bifurcation occurs when this parameter passes through a critical value. The direction of Neimark–Sacker bifurcation and the stability of the bifurcating periodic solution are also discussed. Acknowledgement

r's

This work was jointly supported by the National Natural Science Foundation of China under Grant No. 60574043, and the Natural Science Foundation of Jiangsu Province of China under Grant No. BK2006093.

th o

References

Au

[1] M.B. D’Amico, J.L. Moiola, E.E. Paolini, Hopf bifurcation for maps: a frequency-domain approach, IEEE Trans. Circuits Syst. I 49 (2002) 281–288. [2] A.N. Michel, J.A. Farrel, W. Porod, Qualitative analysis of neural networks, IEEE Trans. Circuits Syst. 36 (1989) 229–243. [3] S. Ruan, J. Wei, Periodic solutions of planar systems with two delays, Proc. R. Soc. Edinburgh A 129 (1999) 1017–1032. [4] T. Faria, On a planar system modelling a neuron network with memory, J. Diﬀerent. Equat. 168 (2000) 129–149. [5] L. Olien, J. Belair, Bifurcations, stability, and monotonicity properties of a delayed neural network model, Physica D 102 (1997) 349– 363. [6] J. Wei, S. Ruan, Stability and bifurcation in a neural network with two delays, Physica D 130 (1999) 255–272. [7] K. Gopalsamy, I. Leung, Delay induced periodicity in a neural network of excitation and inhibition, Physica D 89 (1996) 395–426. [8] L. Sharyer, S.A. Campbell, Stability, bifurcation and multistability in a system of two coupled neurons with multiple time delays, SIAM J. Appl. Math. 61 (2000) 673–700. [9] Z. Yuan, D. Hu, L. Huang, Stability and bifurcation analysis on a discrete-time system of two neurons, Appl. Math. Letts. 17 (2004) 1239–1245. [10] Z. Yuan, D. Hu, L. Huang, Stability and bifurcation analysis on a discrete-time neural network, J. Comput. Appl. Math. 177 (2005) 89–100. [11] Y.A. Kuznetsov, Elements of Applied Bifurcation Theory, second ed., Springer-Verlag, New York, 1998. [12] J. Hale, H. Koc¸ak, Dynamics and Bifurcation, Texts in Applied Mathematics, vol. 3, Springer, New York, 1991.

2122

W. He, J. Cao / Applied Mathematical Modelling 31 (2007) 2111–2122

Au

th o

r's

pe

rs

on

al

co

py

[13] W. Yu, J. Cao, Stability and Hopf bifurcation analysis on a four-neuron BAM neural network with time delays, Phys. Lett. A 351 (2006) 64–78. [14] W. Yu, J. Cao, Stability and Hopf bifurcation on a two-neuron system with time delay in the frequency domain, Int. J. Bifurc. Chaos, in press. [15] W. Yu, J. Cao, Hopf bifurcation and stability of periodic solutions for van der Pol equation with time delay, Nonlin. Anal. 62 (2005) 141–165. [16] J. Hopﬁeld, Neurons with graded response have collective computational properties like two-state neurons, Proc. Natl. Acad. Sci. USA 81 (1984) 3088–3092. [17] Y. Chen, J. Wu, The asymptotic shapes of periodic solutions of singular delay diﬀerential system, J. Diﬀeren. Equat. 169 (2001) 614– 632. [18] J. Cao, M. Xiao, Stability and Hopf bifurcation in a simpliﬁed BAM neural network with two time delays, IEEE Trans. Neural Networks, in press. [19] Z. Yuan, L. Huang, Convergence and periodicity in a discrete-time network of two neurons with self-connections, Comput. Math. Appl. 46 (8–9) (2003) 1337–1345. [20] K.L. Cooke, J. Wiener, Retarded diﬀerential equations with piecewise constant delays, J. Math. Anal. Appl. 99 (1984) 265–297. [21] Z. Yuan, L. Huang, Y. Chen, Convergence and periodicity of solutions for a discrete-time network model of two neurons, Math. Computer Model. 35 (9–10) (2002) 941–950. [22] Z. Zhou, Periodic solutions of a class of diﬀerence systems, Fields Inst. Commun. 43 (2004) 395–402.

Recommend Documents

stability and hopf bifurcation on a two-neuron ... - Semantic Scholar

local and global hopf bifurcation in a delayed ... - Semantic Scholar

Persistence, Stability and Hopf Bifurcation in a ... - UWO Applied Math

asymptotic stability of periodic solutions for ... - Semantic Scholar

Generic Hopf bifurcation from lines of equilibria ... - Semantic Scholar

synchrony-breaking hopf bifurcation in a model of ... - Semantic Scholar