Oscillations of Numerical Solutions for Nonlinear Delay Differential Equations in the Control of Erythropoiesis

Abstract

We consider the oscillations of numerical solutions for the nonlinear delay differential equations in the control of erythropoiesis. The exponential 𝜃-method is constructed and some conditions under which the numerical solutions oscillate are presented. Moreover, it is proven that every nonoscillatory numerical solution tends to the equilibrium point of the continuous system. Numerical examples are given to illustrate the main results.

1. Introduction

The oscillatory and asymptotic behavior of solutions of delay differential equations has been the subject of intensive investigations during the past decades. A large number of articles has appeared in the literature, we refer to [1–4] and the references therein. 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, spread of some infectious diseases in humans, and so on. For more information on this investigation, the reader can see [5, 6] and the references therein.

By contrast with the research on the oscillations of the analytic solutions, much studies have been focused on the oscillations of the numerical solutions for delay differential equations. In [7, 8], oscillations of numerical solutions in 𝜃-methods and Runge-Kutta methods for a linear differential equation with piecewise constant arguments (EPCA, a special type of Delay Differential Equations) 𝑥^′(𝑡)+𝑎𝑥(𝑡)+𝑎₁𝑥([𝑡−1])=0 were considered, respectively. More recently, Wang et al. [9] studied numerical oscillations of alternately advanced and retarded linear EPCA, the conditions of oscillations for the 𝜃-methods are obtained. To the best of our knowledge, until now less attention had been paid for the oscillations of the numerical solutions for nonlinear delay differential equations except for [10]. Differently from [10], in our paper, we will investigate another nonlinear delay differential equation in the control of erythropoiesis and obtain some new results.

Consider the following nonlinear delay differential equation: 𝑥^′(𝑡)=𝛽₀𝜔^𝜇𝜔^𝜇+𝑥^𝜇()(𝑡),𝑡−𝜏−𝛾𝑥 (1) with conditions 𝜇,𝜔,𝛽₀,𝛾>0,𝜏≥0. (2) Mackey and Glass [11] have proposed (1) as model of hematopoiesis (blood cell production). In (1), 𝑥(𝑡) 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 monotonic decreasing function of 𝑥(𝑡−𝜏). Equation (1) has been recently studied by many authors. Gopalsamy et al. [12] obtained sufficient and also necessary and sufficient 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. Using the linearization method, Zaghrout et al. [13] considered (1) and gave a sufficient condition for oscillations of all solutions about the positive steady state 𝑥^∗ and proved that every nonoscillatory positive solution of (1) tends to 𝑥^∗ as 𝑡→∞. For more details of (1), we refer to Mackey and Milton [14] and Mackey [15]. Up to now, few results on the properties of numerical solutions for (1) were obtained. In the present paper, our main goal is to investigate some sufficient conditions under which the numerical solutions are oscillatory. We also consider the asymptotic behavior of nonoscillatory numerical solutions.

The contents of this paper are as follows. In Section 2, some necessary definitions and results for oscillations of the analytic solutions are given. In Section 3, we obtain a recurrence relation by applying the 𝜃-methods to the simplified form which comes from making two transformations on (1). Moreover, the oscillations of the numerical solutions are discussed and conditions under which the numerical solutions oscillate are obtained. In Section 4, we study the asymptotic behavior of nonoscillatory solutions. In Section 5, we present numerical examples that illustrate the theoretical results for the numerical methods. Finally, Section 6 gives conclusions and issues for future research.

2. Preliminaries

In this section, we start by introducing some definitions, lemmas and theorems that will be employed throughout the work.

Definition 1. A function 𝑥(𝑡) of (1) is said to oscillate about 𝐾^∗ if 𝑥(𝑡)−𝐾^∗ has arbitrarily large zeros. Otherwise, 𝑥(𝑡) is called nonoscillatory. When 𝐾^∗=0, we say that 𝑥(𝑡) oscillates about zero or simply oscillates.

Definition 2. A sequence {𝑥_𝑛} is said to oscillate about {𝑦_𝑛} if {𝑥_𝑛−𝑦_𝑛} is neither eventually positive nor eventually negative. Otherwise, {𝑥_𝑛} is called nonoscillatory. If {𝑦_𝑛}={𝑦} is a constant sequence, we simply say that {𝑥_𝑛} oscillates about {𝑦}. When {𝑦_𝑛}={0}, we say that {𝑥_𝑛} oscillates about zero or simply oscillates.

Definition 3. We say (1) oscillates if all of its solutions are oscillatory.

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

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

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

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

Lemma 8 (see [17]). For all 𝑚≥𝑀,(i)(1+𝑎/(𝑚−𝜃𝑎))^𝑚≥𝑒^𝑎 if and only if 1/2≤𝜃≤1 for 𝑎>0,  𝜑(−1)≤𝜃≤1 for 𝑎<0;(ii)(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. Transformation

For (1), we take an initial condition of the form 𝑥(𝑡)=𝜓(𝑡),−𝜏≤𝑡≤0, (6) where 𝜓∈𝐶([−𝜏,0],(0,∞)),𝜓(0)>0.

In order to simplify (1), we introduce a similar method in [13]. The change of variables 𝑥(𝑡)=𝜔𝑦(𝑡) (7) transforms (1) to the delay differential equation 𝑦^′(𝑡)=𝑎1+𝑦^𝜇()(𝑡),𝑡−𝜏−𝛾𝑦 (8) where 𝑎=𝛽₀/𝜔. One can see that (8) has a unique equilibrium 𝐾 and that 𝑎1+𝐾^𝜇=𝛾𝐾. (9) The following result concerning oscillations of the analytic solution of (8) is given in [13].

Theorem 9. Assume that 𝑎𝜇𝐾^𝜇−1𝜏(1+𝐾^𝜇)²𝑒^{𝑎𝜏/𝐾(1+𝐾𝜇)}>1𝑒, (10) then every positive solution of (8) oscillates about its positive equilibrium 𝐾.

Therefore, we obtain the following corollary naturally.

Corollary 10. Assume that the condition 𝑎𝜇𝐾^𝜇−1𝜏(1+𝐾^𝜇)²𝑒^{𝑎𝜏/𝐾(1+𝐾𝜇)}>1𝑒 (11) holds, then every positive solution of (1) oscillates about its positive equilibrium 𝐾^∗=𝜔𝐾.

Next, we also introduce an invariant oscillation transformation 𝑦(𝑡)=𝐾𝑒^𝑧(𝑡), then (8) can be written as 𝑧^′(𝑡)+𝑎𝐾(1+𝐾^𝜇)𝑓₁(𝑧(𝑡))𝑓₂(𝑧(𝑡−𝜏))−𝑓₁(𝑧(𝑡))−𝑓₂(𝑧(𝑡−𝜏))+2=0, (12) where 𝑓₁(𝑢)=1−𝑒^−𝑢,𝑓₂(𝑢)=1+1+𝐾^𝜇1+𝐾^𝜇𝑒^𝜇𝑢. (13) Then 𝑦(𝑡) oscillates about 𝐾 if and only if 𝑧(𝑡) oscillates about zero. Moreover, for the sake of brevity, let 𝑎𝑇=𝐾(1+𝐾^𝜇), (14) then (12) becomes 𝑧^′(𝑡)=−𝑇𝑓₁(𝑧(𝑡))𝑓₂(𝑧(𝑡−𝜏))+𝑇𝑓₁(𝑧(𝑡))+𝑇𝑓₂(𝑧(𝑡−𝜏))−2𝑇. (15) For convenience, denote 𝐿=𝑎𝜇𝐾^𝜇−1(1+𝐾^𝜇)², (16) then the inequality (11) yields 𝐿𝜏𝑒^𝑇𝜏>1𝑒. (17)

3.2. The Difference Scheme

In this subsection we consider the adaptation of the 𝜃-methods. Let ℎ=𝜏/𝑚 be a given step size with integer 𝑚>1. The adaptation of the linear 𝜃-method and the one-leg 𝜃-method to (15) leads to the same numerical process of the following type: 𝑧_𝑛+1=𝑧_𝑛−𝑇ℎ𝜃𝑓₁𝑧_𝑛+1𝑓₂𝑧_{𝑛+1−𝑚}+(1−𝜃)𝑓₁𝑧_𝑛𝑓₂𝑧_𝑛−𝑚+𝑇ℎ𝜃𝑓₁𝑧_𝑛+1+(1−𝜃)𝑓₁𝑧_𝑛+𝑇ℎ𝜃𝑓₂𝑧_{𝑛+1−𝑚}+(1−𝜃)𝑓₂𝑧_𝑛−𝑚−2𝑇ℎ, (18) where 0≤𝜃≤1,  𝑧_𝑛+1 and 𝑧_{𝑛+1−𝑚} are approximations to 𝑧(𝑡) and 𝑧(𝑡−𝜏) of (15) at 𝑡_𝑛+1, respectively.

Let 𝑧_𝑛=ln(𝑥_𝑛/𝐾𝜔), and take into account of the expressions of 𝑓₁ and 𝑓₂ we have 𝑥_𝑛+1=𝑥_𝑛expℎ𝑇𝐾𝜔^𝜇+11+𝐾^𝜇×𝜃𝑥_𝑛+1𝜔^𝜇+𝑥^𝜇_{𝑛+1−𝑚}+1−𝜃𝑥_𝑛𝜔^𝜇+𝑥^𝜇_𝑛−𝑚−ℎ𝑇. (19)

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

The following theorem gives the convergence of exponential 𝜃-method. We can easily prove it by the method of steps which is used in [18].

Theorem 12. The exponential 𝜃-method (19) is convergent with order 11,𝑖𝑓𝜃≠2,12,𝑖𝑓𝜃=2. (20)

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 (19), we only need to consider the oscillations of (18). The following conditions which are taken from [13] will be used in the next analysis: 𝑢𝑓₁(𝑢)>0,for𝑢≠0,lim_𝑢→0𝑓₁(𝑢)𝑢𝑓=1,₂(𝑢)>0,any𝑢,lim_𝑢→0𝑓₂(𝑢)𝑓=2,₁(𝑢)≤𝑢for𝑢>0,𝑓₁(0)𝑓=0,₂(𝑢)≤2𝑢≥0,𝑓₂(0)=2. (21) For (18), its linearized form is given by 𝑧_𝑛+1=𝑧_𝑛−ℎ𝜃𝑇𝑧_𝑛+1−ℎ(1−𝜃)𝑇𝑧_𝑛−ℎ𝜃𝐿𝑧_{𝑛+1−𝑚}−ℎ(1−𝜃)𝑧_𝑛−𝑚, (22) which is equivalent to 𝑧_𝑛+1=1−ℎ(1−𝜃)𝑇𝑧1+ℎ𝜃𝑇^𝑛−ℎ𝜃𝐿𝑧1+ℎ𝜃𝑇^{𝑛+1−𝑚}−ℎ(1−𝜃)𝐿𝑧1+ℎ𝜃𝑇^𝑛−𝑚. (23) It follows from [16] that (18) oscillates if (23) oscillates under the condition (21).

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

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

In the following, we will study whether the exponential 𝜃-method preserves the oscillations of (1). That is, when Corollary 10 holds, we will investigate the conditions under which (19) is oscillatory.

Lemma 15. The characteristic equation of (22) is given by 𝜉=𝑅−ℎ𝑇+𝐿𝜉^−𝑚. (24)

Proof. Let 𝑧_𝑛=𝜉^𝑛𝑧₀ in (22), we have 𝜉^𝑛+1𝑧₀=𝜉^𝑛𝑧₀−ℎ𝜃𝑇𝜉^𝑛+1𝑧₀−ℎ(1−𝜃)𝑇𝜉^𝑛𝑧₀−ℎ𝜃𝐿𝜉^{𝑛+1−𝑚}𝑧₀−ℎ(1−𝜃)𝐿𝜉^𝑛−𝑚𝑧₀, (25) that is 𝜉=1−ℎ𝜃𝑇+𝐿𝜉^−𝑚𝜉−ℎ(1−𝜃)𝑇+𝐿𝜉^−𝑚, (26) which is equivalent to ()𝜉=1−ℎ1−𝜃𝑇+𝐿𝜉^−𝑚1+ℎ𝜃𝑇+𝐿𝜉^−𝑚ℎ=1−𝑇+𝐿𝜉^−𝑚1+ℎ𝜃𝑇+𝐿𝜉^−𝑚. (27) In view of [19], we know that the stability function of the 𝜃-method is 𝑅(𝑥)𝑥=1+,1−𝜃𝑥 (28) thus the characteristic equation of (22) is given by (24). This completes the proof of the lemma.

Lemma 16. If 𝐿𝜏𝑒^𝑇𝜏>1/𝑒, then the characteristic equation (24) has no positive roots for 0≤𝜃≤1/2.

Proof. Let 𝑉(𝜉)=𝜉−𝑅(−ℎ(𝑇+𝐿𝜉^−𝑚)). By Lemma 8, we have 𝑅−ℎ𝑇+𝐿𝜉^−𝑚≤exp−ℎ𝑇+𝐿𝜉^−𝑚,𝜉>0,0≤𝜃≤1/2. (29) Now we will prove that 𝑊(𝜉)=𝜉−exp(−ℎ(𝑇+𝐿𝜉^−𝑚))>0 for 𝜉>0. Suppose the opposite, that is, there exists a 𝜉₀>0 such that 𝑊(𝜉₀)≤0, then we have 𝜉₀≤exp(−ℎ𝑇−ℎ𝐿𝜉₀^−𝑚), and 𝜉^𝑚₀≤exp−𝑇𝜏−𝐿𝜏𝜉₀^−𝑚. (30) Multiplying both sides of the inequality (30) by 𝐿𝜏𝑒^𝑇𝜏𝑒𝜉₀^−𝑚, we obtain 𝐿𝜏𝑒^𝑇𝜏𝑒𝜉₀^−𝑚𝜉^𝑚₀≤𝐿𝜏𝑒^𝑇𝜏𝑒𝜉₀^−𝑚exp−𝑇𝜏−𝐿𝜏𝜉₀^−𝑚, (31) which gives 𝐿𝜏𝑒^𝑇𝜏𝑒≤𝐿𝜏𝜉₀^−𝑚exp1−𝐿𝜏𝜉₀^−𝑚, (32) therefore we have the following two cases.
Case 1. If 1−𝐿𝜏𝜉₀^−𝑚=0, then 𝐿𝜏𝑒^𝑇𝜏𝑒≤1, which contradicts the condition 𝐿𝜏𝑒^𝑇𝜏>1/𝑒.
Case 2. If 1−𝐿𝜏𝜉₀^−𝑚≠0, then in view of Lemma 7, we get exp1−𝐿𝜏𝜉₀^−𝑚<11−1−𝐿𝜏𝜉₀^−𝑚=1𝐿𝜏𝜉₀^−𝑚, (33) that is, 𝐿𝜏𝜉₀^−𝑚exp1−𝐿𝜏𝜉₀^−𝑚<1, (34) so 𝐿𝜏𝑒^𝑇𝜏𝑒<1, which is also a contradiction to 𝐿𝜏𝑒^𝑇𝜏>1/𝑒.
Combining both the cases, by (29) we obtain that for 𝜉>0𝑉(𝜉)=𝜉−𝑅−ℎ𝑇+𝐿𝜉^−𝑚≥𝜉−exp−ℎ𝑇+𝐿𝜉^−𝑚=𝑊(𝜉)>0, (35) which implies that the characteristic equation (24) has no positive roots. The proof of the lemma is complete.

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

Lemma 17. If 𝐿𝜏𝑒^𝑇𝜏>1/𝑒 and 1/2<𝜃≤1, then the characteristic equation (24) has no positive roots for ℎ<ℎ₀, where ℎ₀=𝜏()∞,𝑓𝑜𝑟𝐿𝜏≥1,1+𝑇𝜏+ln𝐿𝜏1+𝑇𝜏(1−ln𝐿𝜏),𝑓𝑜𝑟𝐿𝜏<1. (36)

Proof. Since 𝑅(−ℎ(𝑇+𝐿𝜉^−𝑚)) is an increasing function of 𝜃 when 𝜉>0, then for 𝜉>0 and 1/2<𝜃≤1𝑅−ℎ𝑇+𝐿𝜉^−𝑚=1−ℎ(1−𝜃)𝑇+𝐿𝜉^−𝑚1+ℎ𝜃𝑇+𝐿𝜉^−𝑚≤11+ℎ𝑇+𝐿𝜉^−𝑚. (37) Next, we will prove that the inequality 1𝜉−1+ℎ𝑇+𝐿𝜉^−𝑚>0,𝜉>0 (38) holds under certain conditions.
From (38), it follows that 1𝜉−1+ℎ𝑇+𝐿𝜉^−𝑚=(1+ℎ𝑇)𝜉^1−𝑚1+ℎ𝑇+𝐿𝜉^−𝑚𝜆(𝜉), (39) where 𝜆(𝜉)=𝜉^𝑚−1𝜉1+ℎ𝑇^𝑚−1+ℎ𝐿,1+ℎ𝑇 (40) so we only need to prove 𝜆(𝜉)>0 for 𝜉>0. It is not difficult to know that 𝜆(𝜉) is the characteristic polynomial of the following difference scheme 𝑧_𝑛+1−𝑧_𝑛+ℎ𝐿𝑧1+ℎ𝑇^{𝑛+1−𝑚}+ℎ𝑇𝑧1+ℎ𝑇^𝑛=0. (41) In view of Theorems 4 and 5, we have that 𝜆(𝜉) has no positive roots if and only if ℎ𝐿𝑚1+ℎ𝑇(𝑚−1^𝑚)^𝑚−1>1−ℎ𝑇1+ℎ𝑇^𝑚, (42) which is equivalent to ln𝐿𝜏+(𝑚−1)ln1+1+𝑇𝜏𝑚−1>0. (43) We examine two cases depending on the position of 𝐿𝜏: Either 𝐿𝜏≥1 or 𝐿𝜏<1.
Case 1. If 𝐿𝜏≥1, by 𝑚>1, then (43) holds true.
Case 2. If 𝐿𝜏<1 and 𝜏()ℎ<1+𝑇𝜏+ln𝐿𝜏(),1+𝑇𝜏1−ln𝐿𝜏 (44) then according to Lemma 6 we have ln𝐿𝜏+(𝑚−1)ln1+1+𝑇𝜏()()/()𝑚−1>ln𝐿𝜏+𝑚−11+𝑇𝜏𝑚−1()/()()1+1+𝑇𝜏𝑚−1=ln𝐿𝜏+𝑚−1)(1+𝑇𝜏𝑚+𝑇𝜏>0. (45) Therefore the inequality (38) holds for ℎ<ℎ₀, where ℎ₀=𝜏()∞,for𝐿𝜏≥1,1+𝑇𝜏+ln𝐿𝜏1+𝑇𝜏(1−ln𝐿𝜏),for𝐿𝜏<1. (46) So we get that the following inequality: 𝑉(𝜉)=𝜉−𝑅−ℎ𝑇+𝐿𝜉^−𝑚1≥𝜉−1+ℎ𝑇+𝐿𝜉^−𝑚>0, (47) holds for ℎ<ℎ₀ and 𝜉>0, which implies that the characteristic equation (24) has no positive roots. This completes the proof.

Remark 18. By inequality (43) and condition 𝐿𝜏<1, we have that 𝜏(1+𝑇𝜏+ln𝐿𝜏)1+𝑇𝜏(1−ln𝐿𝜏)>0, (48) thus ℎ₀ is meaningful.

In view of (21), Lemmas 16 and 17, and Theorem 4, we have the first main theorem of this paper.

Theorem 19. If 𝐿𝜏𝑒^𝑇𝜏>1/𝑒, then (19) is oscillatory for ℎℎ<∞,𝑤ℎ𝑒𝑛0≤𝜃≤1/2,₀,𝑤ℎ𝑒𝑛1/2<𝜃≤1, (49) where ℎ₀ is defined in Lemma 17.

4. Asymptotic Behavior of Nonoscillatory Solutions

In this section, we will investigate the asymptotic behavior of nonoscillatory solutions of (19). The following lemma is an important result about asymptotic behavior of (8).

Lemma 20 (see [13]). Let 𝑦(𝑡) be a positive solution of (8), which does not oscillate about 𝐾. Then lim_𝑡→∞𝑦(𝑡)=𝐾.

From the relationship between (8) and (12), we know that the nonoscillatory solution of (12) satisfies lim_𝑡→∞𝑧(𝑡)=0 if Lemma 20 holds. Furthermore, lim_𝑡→∞𝑥(𝑡)=𝐾^∗ is also obtained. Next, we will prove that the numerical solution of (1) can inherit this property.

Lemma 21. Let 𝑧_𝑛 be a nonoscillatory solution of (18), then lim_𝑛→∞𝑧_𝑛=0.

Proof. Without loss of generality, we may assume that 𝑧_𝑛>0 for sufficiently large 𝑛. Then by condition (21) we know that the following inequalities holds true: 𝑓₁𝑧_𝑖>0,𝑓₂𝑧_𝑖−1>0,𝑓₂𝑧_𝑖−2<0 (50) for sufficiently large 𝑖. Moreover, it is can be seen from (18) that 𝑧_𝑛+1−𝑧_𝑛=−ℎ𝜃𝑇𝑓₁𝑧_𝑛+1𝑓₂𝑧_{𝑛+1−𝑚}()−1−ℎ1−𝜃𝑇𝑓₁𝑧_𝑛𝑓₂𝑧_𝑛−𝑚𝑓−1+ℎ𝜃𝑇₂𝑧_{𝑛+1−𝑚}()𝑇𝑓−2+ℎ1−𝜃₂𝑧_𝑛−𝑚−2, (51) which gives 𝑧_𝑛+1−𝑧_𝑛−ℎ𝜃𝑇𝑓₂𝑧_{𝑛+1−𝑚}()𝑇𝑓−2−ℎ1−𝜃₂𝑧_𝑛−𝑚−2<0. (52) Thus 𝑧_𝑛+1−𝑧_𝑛<ℎ𝜃𝑇𝑓₂𝑧_{𝑛+1−𝑚}()𝑇𝑓−2+ℎ1−𝜃₂𝑧_𝑛−𝑚−2<0, (53) then the sequence {𝑧_𝑛} is decreasing, and therefore lim_𝑛→∞𝑧_𝑛=𝜂∈[0,∞). (54) Now we will prove that 𝜂=0. If 𝜂>0, then there exists 𝑁∈𝐍 and 𝜀>0 such that for 𝑛−𝑚>𝑁, 0<𝜂−𝜀<𝑧_𝑛<𝜂+𝜀. Hence 𝜂−𝜀<𝑧_𝑛−𝑚 and 𝜂−𝜀<𝑧_{𝑛−1+𝑚}. So inequality (52) yields 𝑧_𝑛+1−𝑧_𝑛−ℎ𝜃𝑇1+𝐾^𝜇1+𝐾^𝜇𝑒^{𝜇(𝜂−𝜀)}()𝑇−1−ℎ1−𝜃1+𝐾^𝜇1+𝐾^𝜇𝑒^{𝜇(𝜂−𝜀)}−1<0, (55) which implies that 𝑧_𝑛+1−𝑧_𝑛<𝐵<0, where 𝐵=ℎ𝑇𝐾^𝜇1−𝑒^{𝜇(𝜂−𝜀)}1+𝐾^𝜇𝑒^{𝜇(𝜂−𝜀)}. (56) Thus 𝑧_𝑛→−∞ as 𝑛→∞, which is a contradiction to (54). This completes the proof.

Therefore, the second main theorem of this paper is as follows.

Theorem 22. Let 𝑥_𝑛 be a positive solution of (19), which does not oscillate about 𝐾^∗, then lim_𝑛→∞𝑥_𝑛=𝐾^∗.

5. Numerical Experiments

In this section, we will give some numerical examples to illustrate our results.

Firstly, we consider the equation 𝑥^′(𝑡)=21+𝑥⁷()(𝑡),𝑡−2−𝑥 (57) with initial value 𝑥(𝑡)=0.5 for 𝑡≤0. In (57), it is easy to see that condition (11) holds true and 𝐿𝜏≈7>1. That is, the analytic solutions of (57) are oscillatory. In Figures 1–3, we draw the figures of the analytic solutions and the numerical solutions of (57), respectively. The parameters 𝑚=40,  𝜃=0.2 in Figure 2 and 𝑚=20, 𝜃=0.8 in Figure 3. From the two figures, we can see that the numerical solutions of (57) oscillate about 𝐾^∗=1, which are in agreement with Theorem 19.

Figure 1 

The analytic solution of (57).

Figure 2 

The numerical solution of (57) with 𝑚=40 and 𝜃=0.2.

Figure 3 

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

Secondly, we consider 𝑥^′(𝑡)=11+𝑥⁶()(𝑡),𝑡−1−2𝑥 (58) with initial value 𝑥(𝑡)=0.6 for 𝑡≤0. In (58), it is not difficult to see that condition (11) is fulfilled. That is, the analytic solutions of (58) are oscillatory. In Figures 4–7, we draw the figures of the analytic solutions and the numerical solutions of (58), respectively. The parameters 𝑚=25, 𝜃=0.3 in Figure 5, 𝑚=50,  𝜃=0.6 in Figure 6 and 𝑚=5,  𝜃=0.6 in Figure 7. We can see from the three figures that the numerical solutions of (58) oscillate about 𝐾^∗≈0.4929, which are consistent with Theorem 19. On the other hand, by direct calculation, we get ℎ₀≈0.1882. We notice that ℎ=0.02<ℎ₀ and ℎ=0.2>ℎ₀ in Figures 6 and 7, respectively, so the stepsize ℎ₀ is not optimal.

Figure 4 

The analytic solution of (58).

Figure 5 

The numerical solution of (58) with 𝑚=25 and 𝜃=0.3.

Figure 6 

The numerical solution of (58) with 𝑚=50 and 𝜃=0.6.

Figure 7 

The numerical solution of (58) with 𝑚=5 and 𝜃=0.6.

Thirdly, we consider another equation, 𝑥^′(𝑡)=0.2×2.1⁸2.1⁸+𝑥⁸()(𝑡),𝑡−0.15−0.1𝑥 (59) with initial value 𝑥(𝑡)=39 for 𝑡≤0. For (59), it is easy to see that 𝐿𝜏e^𝑇𝜏+1≈0.0505<1, so the condition (11) is not satisfied. That is, the analytic solutions of (59) are nonoscillatory. In Figures 8–10, we draw the figures of the analytic solutions and the numerical solutions of (59), respectively. In Figure 8, we can see that 𝑥(𝑡)→𝐾^∗≈1.6949 as 𝑡→∞. From Figures 9 and 10, we can also see that the numerical solutions of (59) satisfy 𝑥_𝑛→𝐾^∗≈1.6949 as 𝑛→∞. That is, the numerical method preserves the asymptotic behavior of nonoscillatory solutions of (59), which coincides with Theorem 22.

Figure 8 

The analytic solution of (59).

Figure 9 

The numerical solution of (59) with 𝑚=10 and 𝜃=0.4.

Figure 10 

The numerical solution of (59) with 𝑚=15 and 𝜃=0.75.

Finally, according to Definition 14, we can see from these figures that the exponential 𝜃-method preserves the oscillations of (57) and (58) and the nonoscillations of (59), respectively.

All the above numerical examples confirm our theoretical findings.

6. Conclusions

In this paper, we discuss the oscillations of the numerical solutions of a nonlinear delay differential equation in the control of erythropoiesis. The convergent exponential 𝜃-method, namely the linear 𝜃-method and the one-leg 𝜃-method in exponential form, is constructed. We establish some conditions under which the numerical solutions oscillate in the case of oscillations of the analytic solutions. We also prove that nonoscillatory numerical solutions can inherit the corresponding properties of analytic solutions. It is pointed out that the stepsize ℎ₀ in Lemma 17 is not optimal. Therefore, our future work will be devoted to investigating this problem.

Acknowledgments

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