Numerical Oscillations Analysis for Nonlinear Delay Differential Equations in Physiological Control Systems

Abstract

This paper deals with the oscillations of numerical solutions for the nonlinear delay differential equations in physiological control systems. The exponential 𝜃-method is applied to 𝑝^′(𝑡)=𝛽₀𝜔^𝜇𝑝(𝑡−𝜏)/(𝜔^𝜇+𝑝^𝜇(𝑡−𝜏))−𝛾𝑝(𝑡) and it is shown that the exponential 𝜃-method has the same order of convergence as that of the classical 𝜃-method. Several conditions under which the numerical solutions oscillate are derived. Moreover, it is proven that every nonoscillatory numerical solution tends to positive equilibrium of the continuous system. Finally, the main results are illustrated with numerical examples.

1. Introduction

The nonlinear delay differential equation 𝑝^′(𝑡)=𝛽₀𝜔^𝜇𝑝(𝑡−𝜏)𝜔^𝜇+𝑝^𝜇()(𝑡),𝑡−𝜏−𝛾𝑝 (1.1) where 𝜔>0,𝛽₀>𝛾>0,𝜇,𝜏∈𝐑⁺, (1.2) has been proposed by Mackey and Glass [1] as model of hematopoiesis (blood cell production). Here, 𝑝(𝑡) denotes the density of mature cells in blood circulation, 𝜏 is the time delay between the production of immature cells in the bone marrow and their maturation for release in the circulating blood stream, and the production is a single-humped function of 𝑝(𝑡−𝜏). Equation (1.1) has been recently studied by many authors. Mackey and Heidn [2] considered the local asymptotic stability of the positive equilibrium by the well-known technique of linearization. Gopalsamy et al. [3] obtained sufficient and also necessary conditions for all positive solutions to oscillate about their positive steady states. They also obtained sufficient conditions for the positive equilibrium to be a global attractor. For more details of (1.1), we refer to Mackey [4, 5], and Su et al. [6].

Our aim in this paper is to investigate the oscillations of numerical solutions for (1.1). The oscillatory and asymptotic behavior of solutions of delay differential equations has been the subject of intensive investigations during the past decades. The strong interest in this study is motivated by the fact that it has many useful applications in some mathematical models, such as ecology, biology, and spread of some infectious diseases in humans. The general theory and basic results for this paper have been thoroughly studied in [7, 8]. In recent years, much research has been focused on the oscillations of numerical solutions for delay differential equations [9–12]. Until now, very few results dealing with the corresponding behavior for nonlinear delay differential equations have been presented in the literature except for [13]. In [13], the authors investigate the oscillations of numerical solutions for the nonlinear delay differential equation of population dynamics. Different from [13], in our paper, we will consider another nonlinear delay differential equation (1.1) in physiological control systems and obtain some new results. We not only investigate some sufficient conditions under which the numerical solutions are oscillatory but also consider the asymptotic behavior of nonoscillatory numerical solutions.

The structure of this paper is as follows. In Section 2, some necessary definitions and results for oscillations of the analytic solutions are given. In Section 3, we obtain the numerical discrete equation by applying the 𝜃-methods to the simplified form which comes from making two transformations on (1.1). Moreover, the oscillations of the numerical solutions are discussed, and conditions under which the numerical solutions oscillate are obtained. In Section 4, we investigate the asymptotic behavior of nonoscillatory solutions. In Section 5, we present numerical examples that illustrate the theoretical results for the numerical methods.

2. Preliminaries

Let us state some definitions, lemmas, and theorems that will be used throughout this paper.

Definition 2.1. A function 𝑝 of (1.1) is said to oscillate about 𝑀^∗ if 𝑝−𝑀^∗ has arbitrarily large zeros. Otherwise, 𝑝 is called non-oscillatory. When 𝑀^∗=0, we say that 𝑝 oscillates about zero or simply oscillates.

Definition 2.2. A sequence {𝑝_𝑛} is said to oscillate about {𝑧_𝑛} if {𝑝_𝑛−𝑧_𝑛} is neither eventually positive nor eventually negative. Otherwise, {𝑝_𝑛} is called non-oscillatory. If {𝑧_𝑛}={𝑧} is a constant sequence, we simply say that {𝑝_𝑛} oscillates about {𝑧}. When {𝑧_𝑛}={0}, we say that {𝑝_𝑛} oscillates about zero or simply oscillates.

Definition 2.3. We say that (1.1) oscillates if all of its solutions are oscillatory.

Theorem 2.4 (see [14]). Consider the difference equation 𝑎_𝑛+1−𝑎_𝑛+^𝑙_𝑗=−𝑘𝑞_𝑗𝑎_𝑛+𝑗=0, (2.1) and assume that 𝑘,𝑙∈𝐍 and 𝑞_𝑗∈𝐑 for 𝑗=−𝑘,…,𝑙. Then, the following statements are equivalent:(1)every solution of (2.1) oscillates;(2)the characteristic equation 𝜆−1+∑^𝑙_𝑗=−𝑘𝑞_𝑗𝜆^𝑗=0 has no positive roots.

Theorem 2.5 (see [14]). Consider the difference equation 𝑎_𝑛+1−𝑎_𝑛+𝑠𝑎_𝑛−𝑘+𝑞𝑎_𝑛=0, (2.2) where 𝑘>0,𝑝>0, and 𝑞>0. Then, the necessary and sufficient conditions for the oscillation of all solutions of (2.2) are 𝑞∈(0,1) and 𝑠(𝑘+1)^𝑘+1𝑘^𝑘>1−𝑞^𝑘+1. (2.3)

Lemma 2.6. The inequality ln(1+𝑥)>𝑥/(1+𝑥) holds for 𝑥>−1 and 𝑥≠0.

Lemma 2.7. The inequality 𝑒^𝑥<1/(1−𝑥) holds for 𝑥<−1 and 𝑥≠0.

Lemma 2.8 (see [15]). For all 𝑚≥𝑀₀, one has(1)(1+𝑎/(𝑚−𝜃𝑎))^𝑚≥𝑒^𝑎 if and only if 1/2≤𝜃≤1 for 𝑎>0, 𝜑(−1)≤𝜃≤1 for 𝑎<0;(2)(1+𝑎/(𝑚−𝜃𝑎))^𝑚<𝑒^𝑎 if and only if 0≤𝜃<1/2 for 𝑎<0, 0≤𝜃≤𝜑(1) for 𝑎>0,where 𝜑(𝑥)=1/𝑥−1/(𝑒^𝑥−1) and 𝑀₀ is a positive constant.

3. Oscillations of Numerical Solutions

3.1. Two Transformations

In order to study (1.1) conveniently, we will impose two transformations on (1.1) in this subsection.

Together with (1.1), we will consider the initial condition, 𝑝(𝑡)=𝜓(𝑡),−𝜏≤𝑡≤0, (3.1) the initial value problem (1.1), and (3.1) has a unique positive solution for all 𝑡≥0.

We introduce a similar method in [3]. The change of variables 𝑝(𝑡)=𝜔𝑥(𝑡) (3.2) turns (1.1) into the delay differential equation 𝑥^′(𝑡)=𝛽₀𝑥(𝑡−𝜏)1+𝑥^𝜇()(𝑡),𝑡−𝜏−𝛾𝑥 (3.3) with positive equilibrium 𝑀, which is denoted as 𝛽𝑀=⁰−𝛾𝛾^1/𝜇. (3.4) The following theorem gives oscillations of the analytic solution of (3.3).

Theorem 3.1 (see [3]). Assume that 𝛽𝜇>1,⁰𝛾>𝜇,𝑒𝜇−1_𝛾𝜏𝛾𝛽₀𝜇−1𝛽₀1−𝜇𝛾𝜏>𝑒, (3.5)(3.6) then every positive solution of (3.3) oscillates about its positive equilibrium 𝑀.

The following corollary is naturally obtained.

Corollary 3.2. Assume that all the conditions in Theorem 3.1 hold, then every positive solution of (1.1) oscillates about its positive equilibrium 𝑀^∗=𝜔𝑀.

Next, we introduce an invariant oscillation transformation 𝑥(𝑡)=𝑀𝑒^𝑦(𝑡), and then (3.3) can be written as 𝑦^′(𝑡)+𝛽₀1+𝑀^𝜇𝑓₁𝑦(𝑡)𝑓₂𝑦(𝑡−𝜏)−𝑓₁𝑦(𝑡)−𝑓₂𝑦(𝑡−𝜏)+2=0, (3.7) where 𝑓₁(𝑢)=1−𝑒^−𝑢,𝑓₂(𝑢)=1+1+𝑀^𝜇𝑒^𝑢1+𝑀^𝜇𝑒^𝜇𝑢. (3.8) Then, 𝑥(𝑡) oscillates about 𝑀 if and only if 𝑦(𝑡) oscillates about zero.

Moreover, since 𝛽₀1+𝑀^𝜇=𝛾, (3.9) then (3.6) and (3.7) become 𝑒^𝛾𝜏𝛾𝜇𝑀^𝜇1+𝑀^𝜇1−1𝜏>𝑒,𝑦^′(𝑡)=−𝛾𝑓₁𝑦(𝑡)𝑓₂𝑦(𝑡−𝜏)+𝛾𝑓₁𝑦(𝑡)+𝛾𝑓₂𝑦(𝑡−𝜏)−2𝛾, (3.10)(3.11) respectively.

For our convenience, denote 𝑄=𝜇𝑀^𝜇1+𝑀^𝜇−1, (3.12) then the inequality (3.10) yields 𝑒^𝛾𝜏1𝑄𝛾𝜏>𝑒. (3.13)

3.2. The Difference Scheme

Let ℎ=𝜏/𝑚 be a given stepsize with integer 𝑚>1. The adaptation of the linear 𝜃-method and the one-leg 𝜃-method to (3.11) leads to the same numerical process of the following type: 𝑦_𝑛+1=𝑦_𝑛−ℎ𝜃𝛾𝑓₁𝑦_𝑛+1𝑓₂𝑦_{𝑛+1−𝑚}−ℎ(1−𝜃)𝛾𝑓₁𝑦_𝑛𝑓₂𝑦_𝑛−𝑚+ℎ𝜃𝛾𝑓₁𝑦_𝑛+1+ℎ(1−𝜃)𝛾𝑓₁𝑦_𝑛+ℎ𝜃𝛾𝑓₂𝑦_{𝑛+1−𝑚}+ℎ(1−𝜃)𝛾𝑓₂𝑦_𝑛−𝑚−2ℎ𝛾, (3.14) where 0≤𝜃≤1,𝑦_𝑛+1 and 𝑦_{𝑛+1−𝑚} are approximations to 𝑦(𝑡) and 𝑦(𝑡−𝜏) of (3.11) at 𝑡_𝑛+1, respectively.

Letting 𝑦_𝑛=ln(𝑝_𝑛/𝑀^∗) and using the expressions of 𝑓₁ and 𝑓₂, we have 𝑝_𝑛+1=𝑝_𝑛expℎ𝛾𝜔^𝜇1+𝑀^𝜇𝜃𝑝_{𝑛+1−𝑚}𝑝_𝑛+1𝜔^𝜇+𝑝^𝜇_{𝑛+1−𝑚}+(1−𝜃)𝑝_𝑛−𝑚𝑝_𝑛𝜔^𝜇+𝑝^𝜇_𝑛−𝑚−ℎ𝛾. (3.15)

Definition 3.3. We call the iteration formula (3.15) the exponential 𝜃-method for (1.1), where 𝑝_𝑛+1 and 𝑝_{𝑛+1−𝑚} are approximations to 𝑝(𝑡) and 𝑝(𝑡−𝜏) of (1.1) at 𝑡_𝑛+1, respectively.

The following theorem, for the proof of which we refer to [16], allows us to obtain the convergence of exponential 𝜃-method.

Theorem 3.4. The exponential 𝜃-method (3.15) is convergent with order 11,𝑤ℎ𝑒𝑛𝜃≠2,12,𝑤ℎ𝑒𝑛𝜃=2. (3.16)

3.3. Oscillation Analysis

It is not difficult to know that 𝑝_𝑛 oscillates about 𝑀^∗ if and only if 𝑦_𝑛 is oscillatory. In order to study oscillations of (3.15), we only need to consider the oscillations of (3.14). The following conditions which are taken from [3] will be used in the next analysis: 𝑢𝑓₁(𝑢)>0,for𝑢≠0,lim_𝑢→0𝑓₁(𝑢)𝑢𝑓=1,₂(𝑢)>0,forevery𝑢,lim_𝑢→0𝑓₂(𝑢)𝑓=2,₁(𝑢)≤𝑢,for𝑢>0,𝑓₁(0)𝑓=0,₂(𝑢)≤2,for𝑢≥0,𝜇>2,𝑀^𝜇>1,𝑓₂(0)=2. (3.17) The linearized form of (3.14) is given by 𝑦_𝑛+1=𝑦_𝑛−ℎ𝜃𝛾𝑦_𝑛+1−ℎ(1−𝜃)𝛾𝑦_𝑛+ℎ𝜃𝛾1−𝜇𝑀^𝜇1+𝑀^𝜇𝑦_{𝑛+1−𝑚}+ℎ(1−𝜃)𝛾1−𝜇𝑀^𝜇1+𝑀^𝜇𝑦_𝑛−𝑚. (3.18) Then by (3.12), (3.18) gives 𝑦_𝑛+1=1−ℎ(1−𝜃)𝛾𝑦1+ℎ𝜃𝛾^𝑛−ℎ𝜃𝛾𝑄𝑦1+ℎ𝜃𝛾^{𝑛+1−𝑚}−ℎ(1−𝜃)𝛾𝑄𝑦1+ℎ𝜃𝛾^𝑛−𝑚. (3.19) It follows from [14] that (3.14) oscillates if (3.19) oscillates under the condition (3.17).

Definition 3.5. Equation (3.15)   is said to be oscillatory if all of its solutions are oscillatory.

Definition 3.6. We say that the exponential 𝜃-method preserves the oscillations of (1.1) if (1.1) oscillates, then there is a ℎ>0 or ℎ=∞, such that (3.15) oscillates for ℎ<ℎ. Similarly, we say that the exponential 𝜃-method preserves the nonoscillations of (1.1) if (1.1) non-oscillates, then there is a ℎ>0 or ℎ=∞, such that (3.15) nonoscillates for ℎ<ℎ.

In the following, we will study whether the exponential 𝜃-method inherits the oscillations of (1.1). Equivalently, when Corollary 3.2 holds, we will investigate the conditions under which (3.15) is oscillatory.

Lemma 3.7. The characteristic equation of (3.18) is given by 𝜉=𝑅−ℎ𝛾1+𝑄𝜉^−𝑚. (3.20)

Proof. Letting 𝑦_𝑛=𝜉^𝑛𝑦₀ in (3.18), we have 𝜉^𝑛+1𝑦₀=𝜉^𝑛𝑦₀−ℎ𝜃𝛾𝜉^𝑛+1𝑦₀−ℎ(1−𝜃)𝛾𝜉^𝑛𝑦₀+ℎ𝜃𝛾1−𝜇𝑀^𝜇1+𝑀^𝜇𝜉^{𝑛+1−𝑚}𝑦₀+ℎ(1−𝜃)𝛾1−𝜇𝑀^𝜇1+𝑀^𝜇𝜉^𝑛−𝑚𝑦₀, (3.21) that is, 𝜉=1−ℎ𝜃𝛾𝜉1−1−𝜇𝑀^𝜇1+𝑀^𝜇𝜉^−𝑚−ℎ(1−𝜃)𝛾1−1−𝜇𝑀^𝜇1+𝑀^𝜇𝜉^−𝑚, (3.22) which is equivalent to ()𝛾𝜉=1−ℎ1−𝜃1−1−𝜇𝑀^𝜇/1+𝑀^𝜇𝜉^−𝑚1+ℎ𝜃𝛾1−1−𝜇𝑀^𝜇/(1+𝑀^𝜇)𝜉^−𝑚=1−ℎ𝛾1−1−𝜇𝑀^𝜇/1+𝑀^𝜇𝜉^−𝑚1+ℎ𝜃𝛾1−1−𝜇𝑀^𝜇/(1+𝑀^𝜇)𝜉^−𝑚. (3.23) In view of [17], we know that the stability function of the 𝜃-method is 𝑅(𝑥)=1+(1−𝜃)𝑥𝑥1−𝜃𝑥=1+.1−𝜃𝑥 (3.24) By noticing (3.12), thus the characteristic equation of (3.18) is given by (3.20). The proof is completed.

Lemma 3.8. If condition (3.13) holds, then the characteristic equation (3.20) has no positive roots for 0≤𝜃≤1/2.

Proof. Let 𝑓(𝜉)=𝜉−𝑅(−ℎ𝛾(1+𝑄𝜉^−𝑚)). By Lemma 2.8, we know that 𝑅−ℎ𝛾1+𝑄𝜉^−𝑚≤exp−ℎ𝛾1+𝑄𝜉^−𝑚 (3.25) holds for 𝜉>0 and 0≤𝜃≤1/2. In the following, we will prove that 𝑔(𝜉)=𝜉−exp(−ℎ𝛾(1+𝑄𝜉^−𝑚))>0 for 𝜉>0. Suppose the opposite, that is, there exists a 𝜉₀>0 such that 𝑊(𝜉₀)≤0, then we get 𝜉₀≤exp(−ℎ𝛾(1+𝑄𝜉₀^−𝑚)), and 𝜉^𝑚₀≤exp−𝛾𝜏−𝛾𝜏𝑄𝜉₀^−𝑚. (3.26) Multiplying both sides of the inequality (3.26) by 𝑒^𝛾𝜏𝑄𝛾𝜏𝑒𝜉₀^−𝑚, we have 𝑒^𝛾𝜏𝑄𝛾𝜏𝑒≤𝑄𝛾𝜏𝜉₀^−𝑚exp1−𝑄𝛾𝜏𝜉₀^−𝑚 (3.27) Thus, we have the following two cases.
Case  1. If 1−𝑄𝛾𝜏𝜉₀^−𝑚=0, then 𝑒^𝛾𝜏𝑄𝛾𝜏𝑒≤1, which contradicts the condition (3.13).
Case  2. If 1−𝑄𝛾𝜏𝜉₀^−𝑚≠0, then in view of Lemma 2.7, we obtain exp1−𝑄𝛾𝜏𝜉₀^−𝑚<11−1−𝑄𝛾𝜏𝜉₀^−𝑚=1𝑄𝛾𝜏𝜉₀^−𝑚, (3.28) that is, 𝑄𝛾𝜏𝜉₀^−𝑚exp1−𝑄𝛾𝜏𝜉₀^−𝑚<1, (3.29) so 𝑒^𝛾𝜏𝑄𝛾𝜏𝑒<1, which is also a contradiction to (3.13).
In summary, we have, for 𝜉>0, 𝑓(𝜉)=𝜉−𝑅−ℎ𝛾1+𝑄𝜉^−𝑚≥𝜉−exp−ℎ𝛾1+𝑄𝜉^−𝑚=𝑔(𝜉)>0, (3.30) which implies that the characteristic equation (3.20) has no positive roots. The proof is complete.

Without loss of generality, in the case of 1/2<𝜃≤1, we assume that 𝑚>1.

Lemma 3.9. If condition (3.13) holds and 1/2<𝜃≤1, then the characteristic equation (3.20) has no positive roots for ℎ<ℎ, where ℎ=𝜏∞,𝑓𝑜𝑟𝑄𝛾𝜏≥1,1+𝛾𝜏+ln𝑄𝛾𝜏1+𝛾𝜏1−ln𝑄𝛾𝜏,𝑓𝑜𝑟𝑄𝛾𝜏<1. (3.31)

Proof. Since 𝑅(−ℎ𝛾(1+𝑄𝜉^−𝑚)) is an increasing function of 𝜃 when 𝜉>0, then, for 𝜉>0 and 1/2<𝜃≤1, 𝑅−ℎ𝛾1+𝑄𝜉^−𝑚=1−ℎ(1−𝜃)𝛾1+𝑄𝜉^−𝑚1+ℎ𝜃𝛾1+𝑄𝜉^−𝑚≤11+ℎ𝛾1+𝑄𝜉^−𝑚. (3.32) In the following, we will prove that the inequality 1𝜉−1+ℎ𝛾1+𝑄𝜉^−𝑚>0 (3.33) holds under certain conditions.
The left side of inequality (3.33) can be rewritten as 1𝜉−1+ℎ𝛾1+𝑄𝜉^−𝑚=1+ℎ𝛾𝜉^1−𝑚1+ℎ𝛾1+𝑄𝜉^−𝑚Γ(𝜉), (3.34) where Γ(𝜉)=𝜉^𝑚−1𝜉1+ℎ𝛾^𝑚−1+ℎ𝛾𝑄,1+ℎ𝛾 (3.35) so we only need to prove that Γ(𝜉)>0 for 𝜉>0. It is easy to know that Γ(𝜉)=0 is the characteristic equation of the following difference equation𝑦_𝑛+1−𝑦_𝑛+ℎ𝛾𝑄𝑦1+ℎ𝛾^{𝑛+1−𝑚}+ℎ𝛾𝑦1+ℎ𝛾^𝑛=0. (3.36) By Theorems 2.4 and 2.5, we know that Γ(𝜉) has no positive roots if and only if ℎ𝛾𝑄𝑚1+ℎ𝛾(𝑚−1^𝑚)^𝑚−1>1−ℎ𝛾1+ℎ𝛾^𝑚, (3.37) which is equivalent to ln𝑄𝛾𝜏+(𝑚−1)ln1+1+𝛾𝜏𝑚−1>0. (3.38) We examine two cases depending on the position of 𝑄𝛾𝜏: either 𝑄𝛾𝜏≥1 or 𝑄𝛾𝜏<1.
Case  1. If 𝑄𝛾𝜏≥1, by 𝑚>1, (3.38) holds true.
Case  2. If 𝑄𝛾𝜏<1 and 𝜏ℎ<1+𝛾𝜏+ln𝑄𝛾𝜏1+𝛾𝜏1−ln𝑄𝛾𝜏, (3.39) then by Lemma 2.6, we have ln𝑄𝛾𝜏+(𝑚−1)ln1+1+𝛾𝜏()/()𝑚−1>ln𝑄𝛾𝜏+𝑚−11+𝛾𝜏𝑚−1/()()1+1+𝛾𝜏𝑚−1=ln𝑄𝛾𝜏+𝑚−11+𝛾𝜏𝑚+𝛾𝜏>0. (3.40) Therefore, the inequality (3.33) holds for ℎ<ℎ, where ℎ is defined in (3.31). So, we get that the following inequality 𝑓(𝜉)=𝜉−𝑅−ℎ𝛾1+𝑄𝜉^−𝑚1≥𝜉−1+ℎ𝛾1+𝑄𝜉^−𝑚>0 (3.41) holds for ℎ<ℎ and 𝜉>0, which implies that the characteristic equation (3.20) has no positive roots. This completes the proof.

Remark 3.10. By inequality (3.38) and condition 𝑄𝛾𝜏<1, we have 𝜏1+𝛾𝜏+ln𝑄𝛾𝜏1+𝛾𝜏1−ln𝑄𝛾𝜏>0, (3.42) thus ℎ is meaningful.

In view of (3.17), Lemmas 3.8 and 3.9, and Theorem 2.4, we present the first main theorem of this paper.

Theorem 3.11. If condition (3.13) holds, then (3.15) is oscillatory for 1ℎ<ℎ=∞,𝑤ℎ𝑒𝑛0≤𝜃≤2,1ℎ,𝑤ℎ𝑒𝑛2<𝜃≤1, (3.43) where ℎ is defined in Lemma 3.9.

4. Asymptotic Behavior of Nonoscillatory Solutions

In this section, we will study the asymptotic behavior of non-oscillatory solutions of (3.15). We first recall the following result about asymptotic behavior of (3.3).

Lemma 4.1 (see [3]). Assume that 𝜏>0,𝜇>1,𝛽₀>𝛾, (4.1) then every solution 𝑥(𝑡) of the initial value problem is 𝑥^′(𝑡)=𝛽₀𝑥(𝑡−𝜏)1+𝑥^𝜇()(𝑡),𝑥(𝑡)𝑡−𝜏−𝛾𝑥=𝜑,−𝜏≤𝑡≤0, (4.2) which is non-oscillatory about 𝑀 satisfies lim_𝑡→∞𝑥(𝑡)=𝑀. (4.3)

From (3.3) and (3.7), we know that the non-oscillatory solution of (3.7) satisfies lim_𝑡→∞𝑦(𝑡)=0 if Lemma 4.1 is satisfied. Furthermore, lim_𝑡→∞𝑝(𝑡)=𝑀^∗ is also obtained. In the following, we will prove that the numerical solution of (1.1) can inherit this property.

Lemma 4.2. Let 𝑦_𝑛 be a non-oscillatory solution of (3.14), then lim_𝑛→∞𝑦_𝑛=0.

Proof. Without loss of generality, we assume that 𝑦_𝑛>0 for sufficiently large 𝑛. Then by condition (3.17), we have 𝑓₁𝑦_𝑖>0,𝑓₂𝑦_𝑖−1>0,𝑓₂𝑦_𝑖−2<0, (4.4) for sufficiently large 𝑖. Moreover, it is can be seen from (3.14) that 𝑦_𝑛+1−𝑦_𝑛=−ℎ𝜃𝛾𝑓₁𝑦_𝑛+1𝑓₂𝑦_{𝑛+1−𝑚}−1−ℎ(1−𝜃)𝛾𝑓₁𝑦_𝑛𝑓₂𝑦_𝑛−𝑚𝑓−1+ℎ𝜃𝛾₂𝑦_{𝑛+1−𝑚}−2+ℎ(1−𝜃)𝛾𝑓₂𝑦_𝑛−𝑚−2, (4.5) which gives 𝑦_𝑛+1−𝑦_𝑛−ℎ𝜃𝛾𝑓₂𝑦_{𝑛+1−𝑚}−2−ℎ(1−𝜃)𝛾𝑓₂𝑦_𝑛−𝑚−2<0, (4.6) thus we have 𝑦_𝑛+1−𝑦_𝑛<ℎ𝜃𝛾𝑓₂𝑦_{𝑛+1−𝑚}−2+ℎ(1−𝜃)𝛾𝑓₂𝑦_𝑛−𝑚−2<0, (4.7) then the sequence {𝑦_𝑛} is decreasing, and, therefore, lim_𝑛→∞𝑦_𝑛=𝜂∈[0,∞). (4.8) Next, we prove that 𝜂=0. If 𝜂>0, then there exist𝑁∈𝐍 and 𝜀>0 such that 0<𝜂−𝜀<𝑦_𝑛<𝜂+𝜀 for 𝑛−𝑚>𝑁. Thus, 𝑦_𝑛−𝑚>𝜂−𝜀 and 𝑦_{𝑛−𝑚+1}>𝜂−𝜀. So, inequality (4.6) yields 𝑦_𝑛+1−𝑦_𝑛−ℎ𝜃𝛾𝑓₂𝜂−𝜀−2−ℎ(1−𝜃)𝛾𝑓₂𝜂−𝜀−2<0, (4.9) which implies that 𝑦_𝑛+1−𝑦_𝑛<𝐶<0, where 𝐶=ℎ𝛾𝑓₂𝑒𝜂−𝜀−2=ℎ𝛾^𝜂−𝜀−1+𝑀^𝜇𝑒^𝜂−𝜀1−𝑒^{(𝜇−1)(𝜂−𝜀)}1+𝑀^𝜇𝑒^{𝜇(𝜂−𝜀)}. (4.10) Thus, 𝑦_𝑛→−∞ as 𝑛→∞, which is a contradiction to (4.8). This completes the proof.

As a consequence, the second main theorem of this paper is as follows.

Theorem 4.3. Letting 𝑝_𝑛 be a positive solution of (3.15), which does not oscillate about 𝑀^∗, then lim_𝑛→∞𝑝_𝑛=𝑀^∗, where 𝑀^∗ is the positive equilibrium of the continuous system (1.1).

5. Numerical Examples

In order to verify our results, three numerical examples are examined in this section.

Example 5.1. Consider the following equation: 𝑝^′(𝑡)=3×2.6⁵𝑝(𝑡−2.5)2.6⁵+𝑝⁵()(𝑡),𝑡−2.5−1.5𝑝 (5.1) subject to the initial condition 𝑝(𝑡)=1.7,for𝑡≤0. (5.2) In (5.1), it can be seen that condition (3.6) holds true and 𝑄𝛾𝜏≈5.6250>1. That is, the analytic solutions of (5.1) are oscillatory. In Figures 1–3, we draw the figures of the analytic solutions and the numerical solutions of (5.1), respectively. Set 𝑚=25, 𝜃=0.25 in Figure 2 and 𝑚=10, 𝜃=0.6 in Figure 3, respectively. From the two figures, we can see that the numerical solutions of (5.1) oscillate about 𝑀^∗≈2.6, which are in agreement with Theorem 3.11.

Figure 1 

The analytic solution of (5.1).

Figure 2 

The numerical solution of (5.1) with 𝑚=25 and 𝜃=0.25.

Figure 3 

The numerical solution of (5.1) with 𝑚=10 and 𝜃=0.6.

Example 5.2. Let us consider the equation 𝑝^′(𝑡)=1.77×1.98³𝑝(𝑡−2)1.98³+𝑝³()(𝑡),𝑡−2−𝑝 (5.3) with the initial value 𝑝(𝑡)=1.9 for 𝑡≤0. In (5.3), it is easy to see that condition (3.6) is fulfilled and 𝑄𝛾𝜏≈0.6102<1. That is, the analytic solutions of (5.3) are oscillatory. In Figures 4–6, we draw the figures of the analytic solutions and the numerical solutions of (5.3), respectively. Set 𝑚=20, 𝜃=0.4 in Figure 5 and 𝑚=40, 𝜃=0.7 in Figure 6, respectively. We can see from the three figures that the numerical solutions of (5.3) oscillate about 𝑀^∗≈1.8149, which are consistent with Theorem 3.11. On the other hand, by direct calculation, we get ℎ≈1.2568, so the stepsize ℎ is not optimal.

Figure 4 

The analytic solution of (5.3).

Figure 5 

The numerical solution of (5.3) with 𝑚=20 and 𝜃=0.4.

Figure 6 

The numerical solution of (5.3) with 𝑚=40 and 𝜃=0.7.

Example 5.3. Consider the following equation: 𝑝^′(𝑡)=1.5×3.7²𝑝(𝑡−1)3.7²+𝑝²()(𝑡),𝑡−1−0.5𝑝 (5.4) subject to the initial condition 𝑝(𝑡)=9,for𝑡≤0. (5.5) For (5.4), it is easy to see that 𝑄𝛾𝜏𝑒^𝛾𝜏+1≈0.7469<1, so the condition (3.6) is not satisfied. That is, the analytic solutions of (5.4) are non-oscillatory. In Figures 7–9, we draw the figures of the analytic solutions and the numerical solutions of (5.4), respectively. In Figure 7, we can see that 𝑝(𝑡)→𝑀^∗≈5.2326 as 𝑡→∞. From Figures 8 and 9, we can also see that the numerical solutions of (5.4) satisfy 𝑝_𝑛→𝑀^∗≈5.2326 as 𝑛→∞. That is, the numerical method inherits the asymptotic behavior of non-oscillatory solutions of (5.4), which coincides with Theorem 4.3.

Figure 7 

The analytic solution of (5.4).

Figure 8 

The numerical solution of (5.4) with 𝑚=15 and 𝜃=0.3.

Figure 9 

The numerical solution of (5.4) with 𝑚=20 and 𝜃=0.8.

Furthermore, according to Definition 3.6, we can see from these figures that the exponential 𝜃-method preserves the oscillations of (5.1) and (5.3) and the non-oscillations of (5.4), respectively.

Acknowledgments

Q. Wang’s work was supported by the National Natural Science Foundation of China (no. 11201084). The authors would like to thank Professor Mingzhu Liu, Minghui Song, and Dr. Zhanwen Yang for their useful suggestions.