Dynamics of a Discrete Internet Congestion Control Model

Abstract

We consider a discrete Internet model with a single link accessed by a single source, which responds to congestion signals from the network. Firstly, the stability of the equilibria of the system is investigated by analyzing the characteristic equation. By choosing the time delay as a bifurcation parameter, we prove that Neimark-Sacker bifurcations occur when the delay passes a sequence of critical values. Then, the explicit algorithm for determining the direction of the Neimark-Sacker bifurcations and the stability of the bifurcating periodic solutions is derived. Finally, some numerical simulations are given to verify the theoretical analysis.

1. Introduction

Congestion control algorithms and active queue management (AQM) for Internet have been the focus of intense research since the seminal work of Kelly et al. [1]. The primal algorithms they have introduced are based on explicit feedback functions, where the sum of user utilities is maximized within the capacity constraints of the links. In primal algorithms, the users adapt the source rates dynamically based on the congestion signal generated by the link, and the links select a static law to determine the link prices directly from the arrival rates at the links [2]. In [3], a stability condition was provided for a single proportionally fair congestion controller with delayed feedback. Since then, this result was extended to networks in [4, 5] and for different classes of controllers in [6].

In this paper, we consider an Internet model with a single link and single source, which can be formulated as a congestion control system with feedback delay. The model is described by [3, 7, 8] 𝑥^′(𝑡)=𝑘𝑤−𝑥(𝑡−𝜏)𝑝(𝑥(𝑡−𝜏)), (1.1) where 𝑥(𝑡) is the sending rate of the source at time 𝑡, 𝑘 is a positive gain parameter, 𝜏 is the sum of forward and return delays, 𝑤>0 is a target (set point), and the congestion indication function 𝑝(⋅) is increasing, nonnegative, and not identically zero, which can be viewed as the probability that a packet at the resource receives a feedback congestion indication signal. System (1.1) corresponds to a rate control algorithm for the user, which can be used to adjust the sending rate 𝑥(𝑡) so that the expected number of signals received by the user will tend to target 𝑤. A great deal of research has been devoted to the asymptotic behaviour, periodic solutions, global stability, and other properties of this model, for which, we refer to [3, 7–9].

But considering the need of scientific computation and real-time simulation, our interest is focused on the behaviors of discrete dynamics system corresponding to (1.1). Many authors considered the numerical approximation of a scalar delay differential equation by using different numerical methods, such as 𝜃-method, Euler method, and Trapezoidal method (see [10–15]). In this paper, we use the forward Euler scheme to make the discretization for system (1.1).

Moreover, it is also of interest to find what will happen when the congestion control system loses stability. The purpose of this paper is to discuss this version as a discrete dynamical system by using Neimark-Sacker bifurcation theory of discrete systems. We not only investigate the stability of the fixed point and the existence of the Neimark-Sacker bifurcations, but also the stability and direction of the Neimark-Sacker bifurcation of the discrete system.

The paper is organized as follows: in Section 2, we analyze the distribution of the characteristic equation associated with the discrete model and obtain the existence of the local Neimark-Sacker bifurcation. In Section 3, the direction and stability of closed invariant curve from the Neimark-Sacker bifurcation of the discrete delay model are determined by using the theories of discrete systems in [16]. In the final section, some computer simulations are performed to illustrate the analytical results found.

2. Stability Analysis

Let 𝑢(𝑡)=𝑥(𝑡𝜏). Then (1.1) can be rewritten as 𝑢^′(𝑡)=𝜏𝑘𝑤−𝑢(𝑡−1)𝑝(𝑢(𝑡−1)). (2.1) We consider sufficient small step size of the form ℎ=1/𝑚, where 𝑚∈ℤ₊. The Euler method applied to this equation yields the delay difference equation 𝑢_𝑛+1=𝑢_𝑛+𝜏𝑘ℎ𝑤−𝑢_𝑛−𝑚𝑝𝑢_𝑛−𝑚, (2.2) where 𝑢_𝑛 is an approximate value to 𝑢(𝑛ℎ). Let 𝑢^∗ be a positive equilibrium of (2.2), then we have𝑢^∗𝑝𝑢^∗=𝑤. (2.3) Set 𝑦_𝑛=𝑢_𝑛−𝑢^∗, then there follows: 𝑦_𝑛+1=𝑦_𝑛+𝜏𝑘ℎ𝑤−𝑦_𝑛−𝑚+𝑢^∗𝑝𝑦_𝑛−𝑚+𝑢^∗. (2.4) By introducing a new variable 𝑌_𝑛=(𝑦_𝑛,𝑦_𝑛−1,…,𝑦_𝑛−𝑚)^𝑇, we can rewrite (2.4) in the form 𝑌_𝑛+1=𝐹𝑌_𝑛,𝜏, (2.5) where 𝐹=(𝐹₀,𝐹₁,…,𝐹_𝑚)^𝑇, and 𝐹_𝑗=𝑦_𝑛+𝜏𝑘ℎ𝑤−𝑦_𝑛−𝑚+𝑢^∗𝑝𝑦_𝑛−𝑚+𝑢^∗𝑦,𝑗=0_{𝑛−𝑗+1},1≤𝑗≤𝑚. (2.6) Clearly the origin is a fixed point of (2.5), and the linear part of (2.5) is𝑌_𝑛+1=𝐴𝑌_𝑛, (2.7) where 𝑝𝑢𝐴=10⋯0−𝜏𝑘ℎ^∗+𝑢^∗𝑝^′𝑢^∗10⋯0001⋯00⋮⋮⋱⋮⋮00⋮10. (2.8) The characteristic equation of 𝐴 is given by 𝑎(𝜆)∶=𝜆^𝑚+1−𝜆^𝑚+𝜏𝑘ℎ𝑝𝑢^∗+𝑢^∗𝑝^′𝑢^∗=0. (2.9) It is well known that the stability of the zero equilibrium solution of (2.5) depends on the distribution of the zeros of the roots of (2.9). In order to proof the existence of the local Neimark-Sacker bifurcation at fixed point, we need some lemmas as follows.

Lemma 2.1. There exists a 𝜏>0 such that for 0<𝜏<𝜏 all roots of (2.9) have modulus less than one.

Proof. When 𝜏=0, (2.9) becomes 𝜆^𝑚+1−𝜆^𝑚=0. (2.10) The equation has, at 𝜏=0, an 𝑚-fold root 𝜆=0 and a simple root 𝜆=1.
Consider the root 𝜆(𝜏) such that 𝜆(0)=1. This root depends continuously on 𝜏 and is a differential function of 𝜏. From (2.9), we have 𝑑𝜆=𝑝𝑢𝑑𝜏𝑘ℎ^∗+𝑢^∗𝑝^′𝑢^∗𝑚𝜆^𝑚−1−(𝑚+1)𝜆^𝑚,𝑑𝜆=𝑝𝑢𝑑𝜏𝑘ℎ^∗+𝑢^∗𝑝^′𝑢^∗𝑚𝜆^𝑚−1−(𝑚+1)𝜆^𝑚. (2.11) Since 𝑝(⋅) is increasing and nonnegative, we have 𝑑|𝜆|²𝑑𝜏_{𝜏=0,𝜆=1}=𝜆𝑑𝜆+𝑑𝜏𝜆𝑑𝜆𝑑𝜏_{𝜏=0,𝜆=1}=−2𝑘ℎ𝑝𝑢^∗+𝑢^∗𝑝^′𝑢^∗<0. (2.12) So with the increasing of 𝜏>0, 𝜆 cannot cross 𝜆=1. Consequently, all roots of (2.9) lie in the unit circle for sufficiently small positive 𝜏>0, and the existence of the 𝜏 follows.

A Neimark-Sacker bifurcation occurs when a complex conjugate pair of eigenvalues of 𝐴 cross the unit circle as 𝜏 varies. We have to find values of 𝜏 such that there are roots on the unit circle. Denote the roots on the unit circle by 𝑒^𝑖𝜔, 𝜔∈(−𝜋,𝜋]. Since we are dealing with complex roots of a real polynomial, we only need to look for 𝜔∈(0,𝜋]. For 𝜔∈(0,𝜋], 𝑒^𝑖𝜔 is a root of (2.9) if and only if 𝑒^𝑖𝜔−1+𝜏𝑘ℎ𝑝𝑢^∗+𝑢^∗𝑝^′𝑢^∗𝑒^{−𝑖𝑚𝜔}=0. (2.13) Separating the real part and imaginary part from (2.13), there are cos𝜔+𝜏𝑘ℎ𝑝𝑢^∗+𝑢^∗𝑝^′𝑢^∗𝑝𝑢cos𝑚𝜔=1,sin𝜔−𝜏𝑘ℎ^∗+𝑢^∗𝑝^′𝑢^∗sin𝑚𝜔=0. (2.14)(2.15) So 1cos𝜔=1−2(𝜏𝑘ℎ)²𝑝𝑢^∗+𝑢^∗𝑝^′𝑢^∗2. (2.16) From (2.14) and (2.15), we obtain that sin𝑚𝜔cos𝜔+cos𝑚𝜔sin𝜔=sin𝑚𝜔, (2.17) that is, sin(𝑚+1)𝜔=sin𝑚𝜔, (2.18) then there exists a sequence of 𝜔_𝑗∈2𝑗𝜋𝑚,2𝑗+1𝜋𝑚,𝑗=0,1,2,…,𝑚−12, (2.19) where [⋅] denotes the greatest integer function. It is clear that there exists a sequence of the time delay parameters 𝜏_𝑗 satisfying (2.14) and (2.15), according to 𝜔=𝜔_𝑗.

From (2.15), we have 𝜏_𝑗=𝜏𝜔_𝑗=sin𝜔_𝑗𝑘ℎ𝑝(𝑢^∗)+𝑢^∗𝑝^′(𝑢^∗)sin𝑚𝜔_𝑗. (2.20) Differentiating (2.16) with respect to 𝜔, we obtain thatsin𝜔=𝜏𝑘ℎ𝑝𝑢^∗+𝑢^∗𝑝^′𝑢^∗2𝑑𝜏.𝑑𝜔 (2.21) This gives𝑑𝜏𝑑𝜔_{𝜔=𝜔𝑗,𝜏=𝜏𝑗}=sin𝜔_𝑗𝜏_𝑗𝑘ℎ𝑝(𝑢^∗)+𝑢^∗𝑝^′(𝑢^∗)²>0. (2.22) Therefore, 0<𝜏₀<𝜏₁<⋯<𝜏_{[(𝑚−1)/2]}. (2.23)

Lemma 2.2. Let 𝜆_𝑗(𝜏)=𝑟_𝑗(𝜏)𝑒^{𝑖𝜔𝑗(𝜏)} be a root of (2.9) near 𝜏=𝜏_𝑗 satisfying 𝑟_𝑗(𝜏_𝑗)=1 and 𝜔_𝑗(𝜏_𝑗)=𝜔_𝑗. Then 𝑑𝑟²_𝑗(𝜏)𝑑𝜏_{𝜏=𝜏𝑗,𝜔=𝜔𝑗}>0. (2.24)

Proof. From (2.11), we have 𝑑𝑟²_𝑗𝑑𝜏_{𝜏=𝜏𝑗,𝜔=𝜔𝑗}=𝑑𝜆_𝑗²𝑑𝜏_{𝜏=𝜏𝑗,𝜔=𝜔𝑗}=𝜆𝑑𝜆+𝑑𝜏𝜆𝑑𝜆𝑑𝜏_{𝜏=𝜏𝑗,𝜔=𝜔𝑗}=2𝑘ℎ𝑝𝑢^∗+𝑢^∗𝑝^′𝑢^∗𝑚cos𝑚𝜔_𝑗−(𝑚+1)cos(𝑚+1)𝜔_𝑗𝑚𝑒^{𝑖(𝑚−1)𝜔𝑗}−(𝑚+1)𝑒^{𝑖𝑚𝜔𝑗2}. (2.25) From (2.14) and (2.15), we obtain that cos𝑚𝜔_𝑗=1−cos𝜔_𝑗𝜏_𝑗𝑘ℎ𝑝(𝑢^∗)+𝑢^∗𝑝^′(𝑢^∗,)sin𝑚𝜔_𝑗=sin𝜔_𝑗𝜏_𝑗𝑘ℎ𝑝(𝑢^∗)+𝑢^∗𝑝^′(𝑢^∗.) (2.26) It is easy to see that cos(𝑚+1)𝜔_𝑗=cos𝑚𝜔_𝑗cos𝜔_𝑗−sin𝑚𝜔_𝑗sin𝜔_𝑗=cos𝜔_j−1𝜏𝑘ℎ𝑝(𝑢^∗)+𝑢^∗𝑝^′(𝑢^∗.) (2.27)
Using (2.25)–(2.27), we have 𝑑𝜆_𝑗²𝑑𝜏_{𝜏=𝜏𝑗,𝜔=𝜔𝑗}=2(2𝑚+1)1−cos𝜔_𝑗𝜏_𝑗𝑚𝑒^{𝑖(𝑚−1)𝜔𝑗}−(𝑚+1)𝑒^{𝑖𝑚𝜔𝑗2}>0. (2.28) This completes the proof.

Lemma 2.3. Equation (2.9) has a pair of simple roots 𝑒^{±𝑖𝜔𝑗} on the unit circle when 𝜏=𝜏_𝑗,𝑗=0,1,2,…,[(𝑚−1)/2]. Furthermore, if 𝜏∈[0,𝜏₀), then all the roots of (2.9) have modulus less than one; If 𝜏=𝜏₀, then all roots of (2.9) except 𝑒^±𝑖𝜔0 have modulus less than one. But if 𝜏∈(𝜏_𝑗,𝜏_𝑗+1], for 𝑗=0,1,2,…,[(𝑚−1)/2], (2.9) has 2(𝑗+1), roots have modulus more than one.

Proof. From (2.18)–(2.20), we have that (2.9) has toots 𝑒^{±𝑖𝜔𝑗} if and only if 𝜏=𝜏_𝑗 and 𝜔=𝜔_𝑗 given in (2.19) and (2.20). By Lemmas 2.1 and 2.2, we know that if 𝜏∈[0,𝜏₀), then all roots of (2.9) have modulus less than one; if 𝜏=𝜏₀, then all roots of (2.9) except 𝑒^{±𝑖𝜔𝑗} have modulus less than one; furthermore, by the Rouchs theorem (Dieudonné [17, Theorem 9.17.4]), the statement on the number of eigenvalues with modulus more than one is as follows.

Lemma 2.3 immediately leads to stability of the zero solution of (2.4), and equivalently, of the positive fixed point 𝑢=𝑢^∗ of (2.2). So we have the following results on stability and bifurcation in system (2.2).

Theorem 2.4. There exists a sequence of values of the time delay parameter 0<𝜏₀<𝜏₁<⋯<𝜏_{[(𝑚−1)/2]} such that the positive fixed point 𝑢=𝑢^∗ of (2.2) is asymptotically stable for 𝜏∈[0,𝜏₀) and unstable for 𝜏>𝜏₀. Equation (2.2) undergoes a Neimark-Sacker bifurcation at the 𝑢^∗ when 𝜏=𝜏_𝑗,𝑗=0,1,2,…,[(𝑚−1)/2], where 𝜏_𝑗 satisfies (2.20).

3. Direction and Stability of the Neimark-Sacker Bifurcation in Discrete Model

In the previous section, we obtain the conditions under which a family of periodic solutions bifurcate from the steady state at the critical value 𝜏=𝜏_𝑗,𝑗=0,1,2,…,[(𝑚−1)/2]. Without loss of generality, denote the critical value 𝜏=𝜏_𝑗 by 𝜏^∗. For continuous systems, as pointed out in Hassard et al. [18], it is interesting and important to determine the direction, stability, and the period of these periodic solutions bifurcating from the steady state. In this section, following the idea of Hassard et al., we shall study the direction, stability, and the period of the bifurcating periodic solution when 𝜏=𝜏₀ in the discrete Internet congestion control model. The method we used is based on the theories of discrete system by Kuznetsov [16].

Set 𝜏=𝜏₀+𝜈,𝜈∈ℝ. Then 𝜈=0 is a Neimark-Sacker bifurcation value for (2.4). Rewrite (2.4) as 𝑦_𝑛+1=𝑦_𝑛−𝜏𝑘ℎ𝑝𝑢^∗+𝑢^∗𝑝^′𝑢^∗𝑦_𝑛−𝑚+122𝑝^′𝑢^∗+𝑢^∗𝑝^′′𝑢^∗𝑦²_𝑛−𝑚+163𝑝^′′𝑢^∗+𝑢^∗𝑝^′′′𝑢^∗𝑦³_𝑛−𝑚+𝑂𝑦_𝑛−𝑚⁴. (3.1) So system (2.5) is turned into 𝑌_𝑛+1=𝐴𝑌_𝑛+12𝐵𝑌_𝑛,𝑌_𝑛+16𝐶𝑌_𝑛,𝑌_𝑛,𝑌_𝑛+𝑂𝑌_𝑛⁴, (3.2) where𝐵𝑌_𝑛,𝑌_𝑛=𝑏₀𝑌_𝑛,𝑌_𝑛,0,…,0,𝐶𝑌_𝑛,𝑌_𝑛,𝑌_𝑛=𝑐₀𝑌_𝑛,𝑌_𝑛,𝑌_𝑛,𝑏,0,…,0₀𝜙,𝜓=−𝜏𝑘ℎ2𝑝^′𝑢^∗+𝑢^∗𝑝^′′𝑢^∗𝜙_𝑚𝜓_𝑚,𝑐₀𝜙,,𝜓,𝜂=−𝜏𝑘ℎ3𝑝^′′𝑢^∗+𝑢^∗𝑝^′′′𝑢^∗𝜙_𝑚𝜓_𝑚𝜂_𝑚. (3.3)(3.4) Let 𝑞=𝑞(𝜏₀)∈ℂ^𝑚+1 be an eigenvector of 𝐴 corresponding to 𝑒^𝑖𝜔0, then 𝐴𝑞=𝑒^𝑖𝜔0𝑞,𝐴𝑞=𝑒^−𝑖𝜔0𝑞. (3.5) We also introduce an adjoint eigenvector 𝑞^∗=𝑞^∗(𝜏)∈ℂ^𝑚+1 having the properties 𝐴^𝑇𝑞^∗=𝑒^−𝑖𝜔0𝑞^∗,𝐴^𝑇𝑞^∗=𝑒^𝑖𝜔0𝑞^∗, (3.6) and satisfying the normalization ⟨𝑞^∗,𝑞⟩=1, where ⟨𝑞^∗,𝑞⟩=∑^𝑚_𝑖=0𝑞^∗_𝑖𝑞_𝑖.

Lemma 3.1 (see [19]). Define a vector-valued function 𝑞∶ℂ→ℂ^𝑚+1 by 𝑝(𝜉)=𝜉^𝑚,𝜉^𝑚−1,…,1^𝑇. (3.7) If 𝜉 is an eigenvalue of 𝐴, then 𝐴𝑝(𝜉)=𝜉𝑝(𝜉).

In view of Lemma 3.1, we have 𝑞=𝑝𝑒^𝑖𝜔0=𝑒^{𝑖𝑚𝜔0},𝑒^{𝑖(𝑚−1)𝜔0},…,𝑒^𝑖𝜔0,1^𝑇. (3.8)

Lemma 3.2. Suppose that 𝑞^∗=(𝑞^∗₀,𝑞^∗₁,…,𝑞^∗_𝑚)^𝑇 is the eigenvector of 𝐴^𝑇 corresponding to eigenvalue 𝑒^−𝑖𝜔0, and ⟨𝑞^∗,𝑞⟩=1. Then 𝑞^∗=𝐾1,𝑎₀𝑒^{𝑖𝑚𝜔0},𝑎₀𝑒^{𝑖(𝑚−1)𝜔0},…,𝑎₀𝑒^𝑖2𝜔0,𝑎₀𝑒^{𝑖𝜔0𝑇}, (3.9) where 𝑎₀=−𝜏𝑘ℎ(2𝑝^′(𝑢^∗)+𝑢^∗𝑝^′′(𝑢^∗)) and 𝐾=𝑒^{𝑖𝑚𝜔0}+𝑚𝑎₀𝑒^{−𝑖𝜔0−1}. (3.10)

Proof. Assign 𝑞^∗ satisfies 𝐴^𝑇𝑞^∗=𝑧𝑞^∗ with 𝑧=𝑒^−𝑖𝜔0, then the following identities hold 𝑞^∗₀+𝑞^∗₁=𝑒^−𝑖𝜔0𝑞^∗₀,𝑞^∗_𝑘=𝑒^−𝑖𝜔0𝑞^∗_𝑘−1𝑎,𝑘=2,…,𝑚,₀𝑞^∗₀=𝑒^−𝑖𝜔0𝑞^∗_𝑚. (3.11) Let 𝑞^∗_𝑚=𝑎₀𝑒^𝑖𝜔0𝐾, then 𝑞^∗=𝐾1,𝑎₀𝑒^{𝑖𝑚𝜔0},𝑎₀𝑒^{𝑖(𝑚−1)𝜔0},…,𝑎₀𝑒^𝑖2𝜔0,𝑎₀𝑒^{𝑖𝜔0𝑇}. (3.12) From normalization ⟨𝑞^∗,𝑞⟩=1 and computation, (3.10) is held.

Let 𝑎(𝜆) be characteristic polynomial of 𝐴 and 𝜆=𝑒^𝑖𝜔0. Following the algorithms in [16] and using a computation process similar to that in [19], we can compute an expression for the critical coefficient 𝑐₁(𝜏₀), 𝑐₁𝜏₀=𝑔₂₀𝑔₁₁1−2𝜆₀2𝜆²₀−𝜆₀+𝑔₁₁²1−𝜆₀+𝑔₀₂²2𝜆²₀−𝜆₀+𝑔₂₁2, (3.13) where 𝑔₂₀=𝑞^∗,𝐵𝑞,𝑞,𝑔₁₁=𝑞^∗,𝐵𝑞,𝑞,𝑔₀₂=𝑞^∗,𝐵𝑞,𝑞,𝑔₂₁=𝑞^∗,𝐵𝑞,𝜔₂₀+2𝑞^∗,𝐵𝑞,𝜔₁₁+𝑞^∗,𝐶𝑞,𝑞,𝑞,𝜔₂₀=𝑏₀𝑞,𝑞𝑎𝜆²₀𝑝𝜆²₀−𝑞^∗,𝐵𝑞,𝑞𝜆²₀−𝜆₀𝑞−𝑞^∗,𝐵𝑞,𝑞𝜆²₀−𝜆₀𝜔𝑞,₁₁=𝑏₀𝑞,𝑞𝑎(1)𝑝(1)−𝑞^∗,𝐵𝑞,𝑞1−𝜆₀𝑞−𝑞^∗,𝐵𝑞,𝑞1−𝜆₀𝑞. (3.14) By (3.4), (3.8), and Lemma 3.2, we get 𝑏₀𝑞,𝑝𝑒^𝑖2𝜔0=−𝜏𝑘ℎ2𝑝^′𝑢^∗+𝑢^∗𝑝^′′𝑢^∗,𝑏₀𝑞,𝑞=−𝜏𝑘ℎ2𝑝^′𝑢^∗+𝑢^∗𝑝^′′𝑢^∗,𝑏₀𝑞,𝑞=−𝜏𝑘ℎ2𝑝^′𝑢^∗+𝑢^∗𝑝^′′𝑢^∗,𝑐₀𝑞,𝑞,𝑞=−𝜏𝑘ℎ3𝑝^′′𝑢^∗+𝑢^∗𝑝^′′′𝑢^∗,𝑎𝑒^𝑖2𝜔0=𝑒^{𝑖2(𝑚+1)𝜔0}−𝑒^{𝑖2𝑚𝜔0}+𝜏𝑘ℎ2𝑝^′𝑢^∗+𝑢^∗𝑝^′′𝑢^∗,𝑎(1)=𝜏𝑘ℎ2𝑝^′𝑢^∗+𝑢^∗𝑝^′′𝑢^∗,𝑏₀𝑞,𝑝(1)=−𝜏𝑘ℎ2𝑝^′𝑢^∗+𝑢^∗𝑝^′′𝑢^∗. (3.15) Substituting these into (3.13), we have 𝑐₁𝜏₀=𝐾2𝜏𝑘ℎ2𝑝^′𝑢^∗+𝑢^∗𝑝^′′𝑢^∗2𝑎𝑒^𝑖2𝜔0+2𝜏𝑘ℎ2𝑝^′𝑢^∗+𝑢^∗𝑝^′′𝑢^∗2𝑎(1)−𝜏𝑘ℎ3𝑝^′′𝑢^∗+𝑢^∗𝑝^′′′𝑢^∗. (3.16)

Lemma 3.3 (see [20]). Given the map (2.5) and assume that(1)𝜆(𝜏)=𝑟(𝜏)𝑒^{𝑖𝜔(𝜏)}, where 𝑟(𝜏^∗)=1,𝑟^′(𝜏^∗)≠0 and 𝜔(𝜏^∗)=𝜔^∗;(2)𝑒^{𝑖𝑘𝜔∗}≠1 for 𝑘=1,2,3,4;(3)Re[𝑒^{−𝑖𝜔∗}𝑐₁(𝜏^∗)]≠0.Then an invariant closed curve, topologically equivalent to a circle, for map (2.5) exists for 𝜏 in a one-side neighborhood of 𝜏^∗. The radius of the invariant curve grows like 𝑂(|𝜏−𝜏^∗|). One of the four cases below applies:(1)𝑟^′(𝜏^∗)>0,Re[𝑒^{−𝑖𝜔∗}𝑐₁(𝜏^∗)]<0. The origin is asymptotically stable for 𝜏<𝜏^∗ and unstable for 𝜏>𝜏^∗. An attracting invariant closed curve exists for 𝜏>𝜏^∗.(2)𝑟^′(𝜏^∗)>0,Re[𝑒^{−𝑖𝜔∗}𝑐₁(𝜏^∗)]>0. The origin is asymptotically stable for 𝜏<𝜏^∗ and unstable for 𝜏>𝜏^∗. An repelling invariant closed curve exists for 𝜏<𝜏^∗.(3)𝑟^′(𝜏^∗)<0,Re[𝑒^{−𝑖𝜔∗}𝑐₁(𝜏^∗)]<0. The origin is asymptotically stable for 𝜏>𝜏^∗ and unstable for 𝜏<𝜏^∗. An attracting invariant closed curve exists for 𝜏<𝜏^∗.(4)𝑟^′(𝜏^∗)<0,Re[𝑒^{−𝑖𝜔∗}𝑐₁(𝜏^∗)]>0. The origin is asymptotically stable for 𝜏>𝜏^∗ and unstable for 𝜏<𝜏^∗. An repelling invariant closed curve exists for 𝜏>𝜏^∗.

From the discussion in Section 2, we know that 𝑟^′(𝜏^∗)>0; therefore, by Lemma 3.3, we have the following result.

Theorem 3.4. For (2.2), 𝑢=𝑢^∗ is asymptotically stable for 𝜏∈[0,𝜏₀) and unstable for 𝜏>𝜏₀. An attracting (repelling) invariant closed curve exists for 𝜏>𝜏₀ if Re[𝑒^−𝑖𝜔0𝑐₁(𝜏₀)]<0(>0).

4. Computer Simulation

In this section, we will confirm our theoretical analysis by numerical simulation. We choose 𝑝(𝑥)=𝑥/(20−3𝑥) as in [8] and set 𝑘=1.5,𝑚=4,𝑤=1. Then 𝜏₀=1.0159 is the Neimark-Sacker bifurcation value.

Figures 1 and 2 are about delay difference equation (2.2) when step size ℎ=0.25.

(a)

(b)

Figure 1 

(a) The equilibrium 𝑢∗ of (2.2) is asymptotically stable for 𝜏=0.9<𝜏0=1.0159, (b) Phase plot of (2.2) on the plane (𝑢(𝑛),𝑢(𝑛−4)) for 𝜏=0.9<𝜏0=1.0159.

(a)

(b)

Figure 2 

(a) Waveform plot of (2.2) for 𝜏=1.1>𝜏0=1.0159, (b) A bifurcating periodic solution appears on the plane (𝑢(𝑛),𝑢(𝑛−4)) for 𝜏=1.1>𝜏0=1.0159.

In Figure 1, we show the waveform plot and phase plot for (2.2) with initial values 𝑢_𝑗=3(𝑗=0,1,…,4) for 𝜏=0.9<𝜏₀=1.0159. The equilibrium 𝑢^∗=3.2170 of (2.2) is asymptotically stable. In Figure 2, we show the waveform plot and phase plot for (2.2) with initial values 𝑢_𝑗=3(𝑗=0,1,…,4). The equilibrium 𝑢^∗=3.2170 of (2.2) is unstable for 𝜏=1.1>𝜏₀=1.0159. When 𝜏 varies and passes through 𝜏₀=1.0159, the equilibrium loses its stability and a periodic solution bifurcates from the equilibrium for 𝜏=1.1>𝜏₀=1.0159. That is the delay difference equation (2.2) which has a Neimark-Sacker bifurcation at 𝜏₀.

Acknowledgment

The authors are very grateful to the referees for their valuable suggestions, which help to improve the paper significantly.