Dynamics of a Single Species in a Fluctuating Environment under Periodic Yield Harvesting

Abstract

We discuss the effect of a periodic yield harvesting on a single species population whose dynamics in a fluctuating environment is described by the logistic differential equation with periodic coefficients. This problem was studied by Brauer and Sánchez (2003) who attempted the proof of the existence of two positive periodic solutions; the flaw in their argument is corrected. We obtain estimates for positive attracting and repelling periodic solutions and describe behavior of other solutions. Extinction and blow-up times are evaluated for solutions with small and large initial data; dependence of the number of periodic solutions on the parameter 𝜎 associated with the intensity of harvesting is explored. As 𝜎 grows, the number of periodic solutions drops from two to zero. We provide bounds for the bifurcation parameter whose value in practice can be efficiently approximated numerically.

1. Introduction

Environmental conditions like weather or food availability change significantly throughout the year and influence directly the growth of populations. Responding to seasonal environmental fluctuations, population density can alter quite fast during relatively brief periods, reflecting changes in the living conditions that become less favorable or converse. Since in many cases environmental fluctuations have a clearly pronounced seasonal character, they can be efficiently modeled with the help of nonautonomous differential equations with periodic coefficients. A striking example of a positive effect of a periodically fluctuating environment on the dynamics of a species has been reported by Jillson [1] who observed that total population numbers in the flour beetle population in the periodically fluctuating environment were more than twice those in the constant environment. On the other hand, Walters and Bandy [2] demonstrated positive effect of periodic harvesting concluding that periodic harvest of some big game populations may increase the total yield by 10 to 20 percent with the best interval between harvests in 2 to 4 years.

Although importance of the systematic study of the effect of environmental changes on the dynamics of populations has been emphasized in the monographs of MacArtur and Wilson [3], Nisbet and Gurney [4], Renshaw [5], Thieme [6], and other authors, many important problems remain open even for simple cases. As mentioned by Rosenblat [7, page 23], “seasonal and circadian changes in the surrounding conditions... can have a significant effect on birth and death rates, availability of resources, and so on. In spite of this, the question of the influence of these variations has received surprisingly little attention, certainly by comparison with the massive literature devoted to the analysis of systems in constant environments, and even by comparison with the studies of ecosystems in randomly fluctuating environments.” In the introduction to the special issue of the journal Theoretical Population Biology “Understanding the role of environmental variation in population and community dynamics” (volume 64 (2003), issue 3), its editor Chesson [8] stressed that “the dominant focus in theoretical models of population and community dynamics has not been on how populations change in response to the physical environment, but on how populations depend on their own population densities or the population densities of other organisms.”

In this paper, we investigate the effect of a periodic yield harvesting on the dynamics of a population in a fluctuating environment described by 𝑑𝑥(𝑡)(𝑡)𝑥(𝑡)𝑥(𝑡)𝑑𝑡=𝑟1−𝐾(𝑡), (1) where the intrinsic growth rate 𝑟 and the carrying capacity of the environment 𝐾 are positive, continuous functions that vary periodically with time, 𝑟(𝑡+𝑇)=𝑟(𝑡) and 𝐾(𝑡+𝑇)=𝐾(𝑡), for all 𝑡∈ℝ. Logistic equation (1) is widely used by ecologists, although its appropriateness as a model has been questioned, see Gabriel et al. [9]. Equation (1) does not describe correctly behavior of solutions if the condition 𝑟=^def(1/𝑇)∫^𝑇₀𝑟(𝑡)𝑑𝑡>0 fails to hold, see Rogovchenko and Rogovchenko [10]. Nevertheless, as Gabriel et al. [9, page 147] fairly noticed, “independently of the status that one gives to this model, it has been and remains a corner-stone of empirical and theoretical ecology.”

The dynamics of harvested populations in a fluctuating environment has been addressed by several authors. We mention papers by Benardete et al. [11], Brauer and Sánchez [12], and Campbell and Kaplan [13] that stimulated the interest of the authors to the topic, as well as contributions by Lazer [14], Lazer and Sánchez [15, 16], and Liu et al. [17]. Contrary to proportional harvesting, the case where both 𝑟 and 𝐾 in (1) are periodic along with the harvesting term 𝐻 has been studied only by Brauer and Sánchez [12].

A bifurcation problem for a differential equation 𝑑𝑥(𝑡)(𝑡(𝑡((𝑑𝑡=𝑘𝑥)(1−𝑥))−𝐻1+sin2𝜋𝑡)) (2) has been discussed by Campbell and Kaplan [13] and, in more detail, by Benardete et al. [11], cf. [16, Example 1, page 157]. Unfortunately, argument developed by Benardete et al. [11] uses symmetry of the differential equation (2). Additional difficulties arise in case of variable coefficients because, as mentioned by Nkashama [18, page 2], “unlike the constant-coefficient case, the nonlinearity might have a string of non-zero 𝑥-intercepts in time.” Recently, using Crandall-Rabinowitz saddle-node bifurcation theorem, Liu et al. [17] established existence of periodic solutions for 𝑑𝑥(𝑡)/𝑑𝑡=𝑓(𝑥(𝑡))−𝜎ℎ(𝑡). However, neither they provide estimates for periodic solutions nor describe behavior of other solutions. As stressed by Padhi et al. [19, page 2617], “it would be interesting to develop results that identify the exact number of positive periodic solutions admitted by the considered model and study their stability nature. Such study becomes imperative from resource management perspective.”

As Brauer and Sánchez [12, page 243] pointed out, “a general theory of the qualitative behavior of periodic population models, both single species and interacting species, would have many applications.” In this paper, we obtain estimates for positive attracting and repelling periodic solutions to (1) in case of periodic yield harvesting, describe behavior of other solutions, and derive estimates for extinction and blow-up times. This information is important for ecologists who can predict asymptotic behavior of solutions and evaluate their “life span.” We also perform a detailed bifurcation analysis providing bounds for the bifurcation parameter 𝜎_bif; these bounds can be tightened numerically. We are not, however, concerned with optimal harvesting policies; the reader is referred, for example, to the papers by Braverman and Mamdani [20], Castilho and Srinivasu [21], Fan and Wang [22], or Xu et al. [23]. Finally, we note that environmental fluctuations may be also modeled by including deviated arguments in logistic differential equations in a variety of ways, see, for instance, Gopalsamy [24], Zhang and Gopalsamy [25], Gopalsamy et al. [26, 27], and the references cited therein.

Remark 1. For obvious reasons, in population biology, only solutions that take on positive values should be taken into consideration. However, for completeness of mathematical analysis of the problem, we also investigate behavior of solutions that satisfy negative initial conditions or become negative at some instant 𝑡_∗. In the former case, such analysis is completely irrelevant for applications, whereas in the latter case the phrase “solutions decay to −∞” should be interpreted in biological terms as “the population goes extinct”; we provide useful estimates for extinction times.

2. Periodic Solutions and Harvesting

2.1. Constant Yield Harvesting versus Proportional

We start by providing an introductory information regarding harvesting of a single species. In general case, harvesting of a population can be modeled by a differential equation 𝑑𝑥(𝑡)((𝑡((𝑡,𝑑𝑡=𝑓𝑡,𝑥))−ℎ𝑡,𝑥)) (3) where the function 𝑓 describes the growth of unharvested population and the function ℎ provides a law according to which members of the population are removed. Two main harvesting options are described in the literature. A commonly used and widely studied type of harvesting where ℎ is a linear function of population size, ℎ(𝑡,𝑥(𝑡))=𝐻₀𝑥(𝑡), is known as proportional or constant effort harvesting. It often arises in mathematical models of fisheries under the assumption that the catch is proportional to the fishing effort 𝐸, see, for instance, the fundamental monograph of Clark [28]. The principal assumption that the catch is proportional to effort appears to be reasonable in many practical situations yet, it may be questionable for small or exhausted fisheries where much higher fishing effort should be required.

The type of harvesting where members of the population are removed at the constant rate per unit time, that is, ℎ(𝑡,𝑥(𝑡))=𝐻₀, is called constant rate or constant yield harvesting. It arises in situations when a certain quota is specified (fishing or hunting licenses, etc.) and can be also described as regular harvesting at a stock-independent rate. It is reasonable to consider the case where 𝐻₀ is a function of time, 𝐻₀=𝐻₀(𝑡); it can also be a periodic function. In a pioneering paper by Brauer and Sánchez [29], the case of logistic growth with a constant harvesting rate was considered. Since then, quite a few papers dealing with the harvesting of single and competing species have been published. However, as recently mentioned by the same authors [12, pages 233-234], “a plausible situation which has received little attention is when 𝐻(𝑡) is periodic, corresponding to seasonal harvesting such as seasonal open hunting or fishing seasons or crop spraying for parasites.”

In addition to theoretical importance of the study undertaken in this paper, we also stress its practical importance. In fact, a logistic growth model with periodic harvesting 𝑑𝑥(𝑡)(𝑡)𝑥(𝑡)𝑥(𝑡)𝑑𝑡=𝑟1−𝐾(𝑡)−𝐻(𝑡), (4) where 𝐻(𝑡) is a certain piecewise constant function with the period 12 has been used by Laham et al. [30] for the mathematical analysis of the best harvesting strategy for tilapia fish farming at selected fish farms in Malaysia. Furthermore, in the recent report by Keesom et al. [31], a differential equation similar to (2), 𝑑𝑥(𝑡)(𝑡)𝑥(𝑡)𝑑𝑡=𝑘𝑥1−𝐶−𝑎(1+sin(𝑏𝑡)), (5) has been used for determination of an optimal harvesting frequency of fishing cycles 𝑏 required for maintaining a steady population of Alaskan salmon. Since both Keesom et al. [31] and Laham et al. [30] provide only numerical analysis for the models mentioned above, importance of a comprehensive theoretical analysis for this class of equations becomes obvious.

In what follows, we employ concepts of lower and upper fences, also termed lower and upper solutions (subsolutions and supersolutions). Basic facts regarding fences and funnels can be found in Hubbard and West [32, Chapters 1 and 4]. We interpret attractors and repellers using forward and pullback convergence, see Wiggins [33, page 112] and use the definition in Berger and Siegmund [34, Definition 3, page 3792], described by the authors as a “tailor-made specialization of more general concepts.”

2.2. Proportional Harvesting

The case of periodic proportional harvesting, 𝑑𝑥(𝑡)(𝑡)𝑥(𝑡)𝑥(𝑡)𝑑𝑡=𝑟1−𝐾(𝑡)−𝐻(𝑡)𝑥(𝑡), (6) where 𝐻 is a continuous positive periodic function has been studied by Castilho and Srinivasu [21, Proposition 3.2, page 6] and Fan and Wang [22, Theorem 4.1, page 172], cf. also Hale and Koçak [35, Exercise 4.23 on pages 128-129] and Sánchez [36, discussion on pages 884-885]. Clearly, (6) can be recast in the form (1) with a modified intrinsic growth rate 𝑟₁(𝑡)=𝑟(𝑡)−𝐻(𝑡) and carrying capacity 𝐾₁(𝑡)=(𝐾(𝑡)(𝑟(𝑡)−𝐻(𝑡)))/𝑟(𝑡). Consequently, an application of results reported by Coleman et al. [37], Fan and Wang [22, Theorem 2.1, pages 167-168], or S. P. Rogovchenko and Y. V. Rogovchenko [10, Theorem 13, page 1176] yields existence of a unique positive 𝑇-periodic solution 𝑥_𝑝(𝑡) to (6) which is asymptotically stable with the domain of attraction containing positive initial data provided that the time average 𝑟₁ of the function 𝑟₁(𝑡) is positive, 𝑟₁>0.

2.3. Periodic Yield Harvesting

Compared to proportional harvesting, the case of periodic yield harvesting, (4) received much less attention. Results where at least one of the three functions 𝑟,  𝐾, or 𝐻 is constant are known, see, Braverman and Mamdani [20], Fan and Wang [22], Xu et al. [23]. To the best of our knowledge, (4) with all three periodic coefficients has been studied only by Brauer and Sánchez [12].

Contrary to (6), Riccati differential equation (4) does not have the trivial solution 𝑥_triv(𝑡)≡0, unless 𝐻(𝑡)≡0. This distinguishes dynamics of solutions with small positive initial data. Although any two solutions of (4) with 0<𝑥₁(𝑡₀)<𝑥₂(𝑡₀) satisfy 𝑥₁(𝑡)<𝑥₂(𝑡), for all 𝑡≥𝑡₀, for any solution 𝑥_small(𝑡) of this equation with 0<𝑥_small(𝑡₀)<𝜀, for a small 𝜀>0, there exists a 𝑡_ext, termed the “extinction time,” such that 𝑥_small(𝑡_ext)=0 and 𝑥_small(𝑡)<0, for all 𝑡>𝑡_ext, see Section 2.5. This cannot happen for (6), for which the unstable solution 𝑥_triv(𝑡) acts as a nonporous lower fence.

Sánchez [38, page 959] mentioned that a minor modification of the result due to Pliss [39, Theorem 9.6, pages 102-103] yields existence of at most two periodic solutions to equations with the right-hand side quadratic with respect to 𝑥, cf. [40, Theorem on page 30]. Consequence of the celebrated Massera Theorem [41], see Brauer and Sánchez [12, Theorem 1, page 234] or Sánchez [40, Theorem on page 34], is often used for establishing existence of periodic solutions to equations with a continuous, 𝑇-periodic right-hand side. As Brauer and Sánchez [12, page 234] pointed out, this 𝑇-periodic solution is asymptotically stable. Thus, we concentrate our efforts on establishing bounds for positive periodic solutions and analyze behavior of other solutions, cf. Rizaner and Rogovchenko [42].

In the sequel, we use notation 𝑔_min=min_{0≤𝑡≤𝑇}𝑔(𝑡) and 𝑔_max=max_{0≤𝑡≤𝑇}𝑔(𝑡) and assume that at least one of the inequalities 𝑟_min≤𝑟_max,𝐾_min≤𝐾_max,𝐻_min≤𝐻_max is strict. Completing the square on the right-hand side of (4) as in Brauer and Sánchez [12, Section 4, pages 241-242], one concludes that the slope 𝑑𝑥(𝑡)/𝑑𝑡 is negative for all 𝑡∈ℝ provided that 𝛾(𝑡)^def(𝑡)−𝑟(𝑡)𝐾(𝑡)=𝐻4>0, (7) in which case 𝑥(𝑡)→−∞ as 𝑡→+∞. Thus, (4) has no periodic solutions and the population goes extinct.

Passing to the case when (7) fails to hold, Brauer and Sánchez [12, Section 4, pages 241-242] reasoned as follows. Assuming that the inequality 𝛾(𝑡)<0 (8) holds for all 𝑡∈ℝ, they argued that, for all 𝑡,𝑥^′|_𝑥=0<0,𝑥^′|_𝑥=𝐾max<0 and 𝑥^′|_{𝑥=𝐾max/2}>0. However, one can easily construct counter examples where the latter inequality does not hold, although (8) is satisfied. In fact, consider (4) with 𝑟(𝑡)=4+sin(2𝜋𝑡),𝐾(𝑡)=5+2sin(2𝜋𝑡) and 𝐻(𝑡)=3+sin(2𝜋𝑡). Then, for all 𝑡∈ℝ,𝛾(𝑡)≤−1/4, and (8) is satisfied, but the derivative of solution of the given differential equation, 𝑥^′_{𝑥=𝐾max/2}=7(4+sin(2𝜋𝑡))2×71−2((−((,5+2sin2𝜋𝑡))3+sin2𝜋𝑡)) (9) changes sign infinitely many times on ℝ.

Theorem 2. Assume that condition 𝐻_max<𝑟_min𝐾_min4 (10) is satisfied. Then (4) has two positive periodic solutions, 𝑥₊(𝑡) and 𝑥₋(𝑡), which are the forward and pullback attractors, respectively, and, for all 𝑡∈ℝ, 0<𝑥₋(𝑡)<𝐾_min2<𝑥₊(𝑡)<𝐾_max. (11)

Proof. Observe first that [40, Theorem on page 30] yields that (4) cannot have more than two periodic solutions. Furthermore, for all 𝑡∈ℝ, 𝑑𝑥(𝑡)𝑑𝑡_{𝑥(𝑡)=𝐾max}=𝑟(𝑡)𝐾_max𝐾1−^max𝐾(𝑡)(𝑡)(𝑡)𝐾−𝐻≤𝑟_max𝐾1−^max𝐾_max(𝑡)(𝑡)−𝐻=−𝐻<0. (12) On the other hand, by virtue of (10), 𝑑𝑥(𝑡)𝑑𝑡_{𝑥(𝑡)=𝐾min/2}(𝑡)𝐾=𝑟^min2𝐾1−^min(𝑡)(𝑡)2𝐾−𝐻≥𝑟_min𝐾_min2𝐾1−^min2𝐾_min−𝐻_max=𝑟_min𝐾_min4−𝐻_max>0. (13) By [40, Theorem on page 34], there exists a periodic solution 𝑥₊(𝑡) of (4) satisfying 𝐾_min/2<𝑥₊(𝑡)<𝐾_max. Direction field and uniqueness arguments imply that this solution is a forward attractor for all solutions of (4) with initial data satisfying 𝐾_min/2<𝑥(𝑡₀)<𝐾_max.
Keeping in mind that a pullback attractor is a forward repeller, see Rasmussen [43, Remark 3.4, page 272], introduce a new variable 𝜏=−𝑡. If 𝑥(𝑡) is a solution to (4), then 𝑥^∼(𝜏)=𝑥(−𝑡) satisfies 𝑑𝑥^∼(𝜏)(𝜏)𝑑𝜏=−𝑟𝑥^∼(𝜏)1−^∼𝑥(𝜏)𝐾(𝜏)+𝐻(𝜏). (14) Note first that, for 𝑥^∼(𝜏)=0, the slope 𝑑𝑥^∼(𝜏)/𝑑𝜏 is positive for all 𝜏∈ℝ, 𝑑𝑥^∼(𝜏)/𝑑𝜏|^∼_𝑥(𝜏)=0=𝐻(𝜏)>0. Furthermore, by virtue of (10), 𝑑𝑥^∼(𝜏)𝑑𝜏^∼_{𝑥(𝜏)=𝐾min/2}(𝜏)𝐾=−𝑟^min2𝐾1−^min(𝜏)(𝜏)2𝐾+𝐻≤−𝑟_min𝐾_min2𝐾1−^min2𝐾_min+𝐻_max=𝐻_max−𝑟_min𝐾_min4<0. (15) Therefore, another application of [40, Theorem on page 34] yields the existence of a unique forward attractor of (14) as 𝜏→+∞. Correspondingly, there exists a unique positive pullback attractor 𝑥₋(𝑡) satisfying (11).

2.4. Sharper Estimates for Periodic Solutions

Rough preliminary estimates (11) for the attractor-repeller pair are further improved in this section. For 𝑖=1,2, let 𝜑_𝑖 and 𝜓_𝑖 be equilibrium solutions to autonomous differential equations 𝑑𝜑(𝑡)𝑑𝑡=𝑟_min𝜑(𝑡)𝜑(𝑡)1−𝐾_min−𝐻_max,𝑑𝜓(𝑡)𝑑𝑡=𝑟_max𝜓(𝑡)𝜓(𝑡)1−𝐾_max−𝐻_min. (16)(17)

Theorem 3. Assume that (10) holds for all 𝑡∈ℝ. Suppose also that 2𝐻_min𝐻_max≤𝑟_max𝑟_min. (18) Then, for the attractor-repeller pair 𝑥₊(𝑡),𝑥₋(𝑡) in Theorem 2, one has 0<𝜓₁<𝑥₋(𝑡)<𝜑₁,𝜑₂<𝑥₊(𝑡)<𝐾_max. (19) If, in addition, 𝐾_max<2𝐾_min, (20) then 𝜑₂<𝑥₊(𝑡)<𝜓₂<𝐾_max. (21)

Proof. By (11), it suffices to consider only the values of 𝑥 between 0 and 𝐾_max. However, one can completely control the behavior of the right-hand side of (4) only for 0≤𝑥≤𝐾_min, in which case the following estimates hold: 𝑟_min𝑥(𝑡)𝑥(𝑡)1−𝐾_min−𝐻_max(𝑡)𝑥(𝑡)𝑥(𝑡)≤𝑟1−𝐾(𝑡)−𝐻(𝑡)≤𝑟_max𝑥(𝑡)𝑥(𝑡)1−𝐾_max−𝐻_min. (22) Equilibrium solutions to differential equations (16) and (17) are given by 𝜑_1,2=𝐾_min21∓1−4𝐻_max𝑟_min𝐾_min,𝜓_1,2=𝐾_max21∓1−4𝐻_min𝑟_max𝐾_max. (23) Condition (18) ensures that these four equilibria are ordered as follows: 𝜓₁<𝜑₁<𝜑₂<𝜓₂. (24) Indeed, note first that 𝜑₁<𝜑₂<𝐾_min. Then, by virtue of 𝐻_min/(𝑟_max𝐾_max)≤𝐻_min/(𝑟_min𝐾_min)≤𝐻_max/(𝑟_min𝐾_min)<1/4, we conclude that 𝜑₂<𝜓₂. To demonstrate that 𝜓₁<𝜑₁, we keep 𝐻_max,𝐻_min,𝑟_max, and 𝑟_min fixed and analyze the behavior of 𝜑₁(𝐾_min) and 𝜓₁(𝐾_max) described by the function 𝜒(𝑥)=𝑥(1−1−𝜈/𝑥)/2,𝑥≥𝜈. One can verify that 𝜒 is strictly decreasing on [𝜈,+∞), attains its maximal value 𝜈/2 at 𝑥=𝜈, and decays to 𝜈/4 as 𝑥→+∞. Correspondingly, 𝜑₁(𝐾_min)∈(𝐻_max/𝑟_min,2𝐻_max/𝑟_min] and 𝜓₁(𝐾_max)∈(𝐻_min/𝑟_max,2𝐻_min/𝑟_max]. Finally, condition (18) guarantees that max_𝐾max𝜓₁(𝐾_max)=2𝐻_min/𝑟_max<𝐻_max/𝑟_min=inf_𝐾min𝜑₁(𝐾_min), which completes the proof of (24).
Consider now the functions 𝛼₁(𝑡)=𝜑₁ and 𝛽₁(𝑡)=𝜓₁. One has 𝑓𝑡,𝛼₁(𝑡)=𝑟(𝑡)𝜑₁1−𝜑₁𝐾(𝑡)−𝐻(𝑡)≥𝑟_min𝜑₁1−𝜑₁𝐾_min−𝐻_max=𝛼^′₁(𝑡),𝑓𝑡,𝛽₁(𝑡)=𝑟(𝑡)𝜓₁1−𝜓₁𝐾(𝑡)−𝐻(𝑡)≤𝑟_max𝜓₁1−𝜓₁𝐾_max−𝐻_min=𝛽^′₁(𝑡). (25) Hence, 𝛼₁(𝑡)=𝜑₁ and 𝛽₁(𝑡)=𝜓₁ are lower and upper fences, whereas the horizontal strip bounded by 𝜑₁ and 𝜓₁ is an antifunnel, and there exists a periodic solution to (4) located between 𝛼₁(𝑡) and 𝛽₁(𝑡). For the repeller 𝑥₋(𝑡), one has, for all 𝑡∈ℝ,𝜓₁<𝑥₋(𝑡)<𝜑₁. If (20) is satisfied, one can show that both estimates for the attractor 𝑥₊(𝑡) are also “tightened” to (21). Otherwise, only a lower bound can be improved, which leads to (19).

Numerical values in the following example, as well as in the rest of the paper are truncated to four decimal places.

Example 4. Consider a 1-periodic differential equation 𝑑𝑥(𝑡)=((𝑥(𝑡)(𝑡)𝑑𝑡5+sin2𝜋𝑡))1−2𝑥()−123+5cos2𝜋𝑡2(11−3cos(2𝜋𝑡)). (26) In this case, 𝜑_1,2=323∓√2,𝜓_1,2√=71∓357.21 (27) Condition (18) is satisfied because 2𝐻_min/𝐻_max=2⋅4/7≤6/4=𝑟_max/𝑟_min. By (19), 7(1−√357/21)<𝑥₋(𝑡)<3(3−√2)/2<3(3+√2)/2<𝑥₊(𝑡)<7(1+√357/21), see Figure 1. The improvement achieved in comparison with rough estimates provided by Theorem 2 is seen as a light-colored “corridor” separating the funnel and antifunnel containing two periodic solutions. Corridor's width is given by 𝜑₂−𝜑₁=𝐾_min1−4𝐻_max/(𝑟_min𝐾_min); it equals 3√2 for (26). In addition, one can see two light-colored “corridors” of the same width 14−7(1+√357/21)=7(1−√357/21)=0.7019 corresponding to improved upper and lower estimates for 𝑥₊(𝑡) and 𝑥₋(𝑡).

Figure 1 

Tighter estimates for both periodic solutions of (26).

Remark 5. Exact solutions 𝜑(𝑡) and 𝜓(𝑡) to (16) and (17) give tighter estimates for the solution 𝑥(𝑡) of (4) satisfying the same initial condition, provided that 0≤𝑥(𝑡₀)=𝜑(𝑡₀)=𝜓(𝑡₀)≤𝐾_min or (20) holds. With the growth of 𝑡, both 𝜑(𝑡) and 𝜓(𝑡) approach equilibria quite fast; even for relatively small values of 𝑡, the difference becomes hardly visible, see Figure 2. Thus, from the practical point of view, one can use estimates (19) and (21) excluding a small interval where solutions of (4) approach periodic ones as 𝑡→∓∞.

Figure 2 

Unstable and stable periodic solutions embraced by upper and lower solutions.

Remark 6. To describe behavior of solutions 𝑥(𝑡) with initial data satisfying 𝑥(𝑡₀)<0 or 𝑥(𝑡₀)>𝐾_max, one needs estimates 𝑟_max𝑥(𝑡)𝑥(𝑡)1−𝐾_min−𝐻_max(𝑡)𝑥(𝑡)𝑥(𝑡)≤𝑟1−𝐾(𝑡)−𝐻(𝑡)≤𝑟_min𝑥(𝑡)𝑥(𝑡)1−𝐾_max−𝐻_min (28) rather than (22), because the first term on the right-hand side of (4) takes on negative values. Consequently, a pair of autonomous differential equations 𝑑𝜂(𝑡)𝑑𝑡=𝑟_max𝜂(𝑡)𝜂(𝑡)1−𝐾_min−𝐻_max,𝑑𝜉(𝑡)𝑑𝑡=𝑟_min𝜉(𝑡)𝜉(𝑡)1−𝐾_max−𝐻_min, (29)(30) is used instead of (16) and (17).

2.5. Extinction Times

We start with an upper bound for the extinction time 𝑡_ext for solutions with initial data 𝑥(𝑡₀)∈(0,𝐾_min]. By (22), one has to study the behavior of solutions to differential equation (17). Suppose that 𝐻_min>𝐻_crit=𝑟_max𝐾_max4. (31) Note that (31) yields (7). Therefore, the slope defined by (17) is negative and all solutions of (4) decay to −∞. Integrating (4), one derives the formula for the extinction time, 𝑡^𝑠_ext=𝑡₀+2𝑟⁻¹_max𝐻_min𝐻_crit−1^−1/2×tan⁻¹𝐻_min𝐻_crit−1^−1/2+tan⁻¹2𝐾_max𝜓𝑡₀−𝐾_max2×𝐻_min𝐻_crit−1^−1/2. (32) In particular, if (20) holds, the extinction time for the solution 𝑥_∗(𝑡) of (4) satisfying 𝑥_∗(𝑡₀)=𝐾_max/2 is determined by a simpler expression, 𝑡^𝑠_ext=𝑡₀+2𝑟⁻¹_max𝐻_min𝐻_crit−1^−1/2tan⁻¹𝐻_min𝐻_crit−1^−1/2. (33)

For solutions with initial data 𝑥(𝑡₀)>𝐾_max, the situation is similar. Assume that (31) holds. By virtue of (28), differential equation (30) is used. Then, 𝐻_max≥𝐻_min>𝐻_crit=𝑟_max𝐾_max/4≥𝑟_min𝐾_max/4=𝐻_∗. Integrating (30), letting 𝜉(𝑡)=0 and solving the resulting equation for 𝑡^𝑙_ext, one has 𝑡^𝑙_ext=𝑡₀+2𝑟⁻¹_min𝐻_min𝐻_∗−1^−1/2×tan⁻¹𝐻_min𝐻_∗−1^−1/2+tan⁻¹2𝐾_max𝜉𝑡₀−𝐾_max2×𝐻_min𝐻_∗−1^−1/2. (34)

Example 7. Consider a 1-periodic differential equation 𝑑𝑥(𝑡)()𝑑𝑡=exp1+cos2𝜋𝑡2(𝑡)𝑥(𝑡)×𝑥1−((−(((.6+coshcos2𝜋𝑡))7+sinhsin2𝜋𝑡))) (35) Condition (31) holds since 5.8248=𝐻_min>𝑟_max𝐾_max/4=5.1260. Consequently, all solutions of (35) decay to −∞. By (32), the extinction time for the solution satisfying 𝑥(0)=1<𝐾_min is 𝑡^𝑠_ext=0.2229, see Figure 3. For a solution of (35) satisfying 𝑥(0)=8>𝐾_max=6+cosh1=7.5431, one has 8.1752=𝐻_max>𝑟_min𝐾_max/4=1.8857. By (34), an upper bound for the extinction time is 𝑡^𝑙_ext=1.7506, see Figure 4.

Figure 3 

The upper solution (blue) and solution (orange) to (35) satisfying the initial condition 𝑥(0)=1. An upper bound for the extinction time for this solution is 𝑡𝑠ext=0.2229.

Figure 4 

The upper solution (blue) and solution (orange) to (35) satisfying the initial condition 𝑥(0)=8. An upper bound for the extinction time for this solution is 𝑡𝑙ext=1.7506.

2.6. Forward and Backward Blow-Up Time

In this section, we assume that (10) holds for all 𝑡∈ℝ. First, we estimate forward blow-up time for “small” solutions with initial data 𝑥(𝑡₀)∈(0,𝜓₁), that is, for solutions located below the repeller 𝑥₋(𝑡). These solutions decay rapidly to −∞ and have vertical asymptotes to the right of 𝑡=𝑡₀. Taking into account estimates (22) that hold for small values of 𝑥(𝑡₀), we integrate differential equation (17) from 𝑡₀ to 𝑡 obtaining 𝜓(𝑡)=𝜓₁𝑐₀𝑒^{𝜘(𝑡−𝑡0)}−𝜓₂𝑐₀𝑒^{𝜘(𝑡−𝑡0)},−1 (36) where 𝑐₀=(𝜓(𝑡₀)−𝜓₂)/(𝜓(𝑡₀)−𝜓₁),  𝜘=−𝑟_max(𝜓₂−𝜓₁)/𝐾_max<0, and equilibrium solutions 𝜓₁ and 𝜓₂ are defined above. By (36), solutions with small initial data 0<𝑥(𝑡₀)<𝜓₁ blow up in the future, 𝑥(𝑡)→−∞ as 𝑡→𝑡⁻_forw; an estimated value for a forward blow-up time 𝑡_forw is given by 𝑡_forw=𝑡₀+𝐾_maxln𝑐₀⁻¹𝑟_max𝜓₁−𝜓₂. (37) Equation 𝑡=𝑡_forw defines a vertical asymptote for solutions to (17); 𝑡_forw>𝑡₀ because, for 𝑥(𝑡₀)<𝜓(𝑡₀)<𝜓₁, one always has 𝑐₀⁻¹<1.

Example 8. To estimate a forward blow-up time for the solution to (26) satisfying the initial condition 𝑥(0)=0.5, note that equilibria of differential equation (17) associated with (26) are given by (27), 𝑐₀=63.4036. By virtue of (37), an upper bound for a forward blow-up time for the solution of (26) passing through the point (0,0.5) is 𝑡_forw=0.7686, see Figure 5.

Figure 5 

The upper solution (blue), lower (purple), and solution (orange) to (26), all satisfying the initial condition, 𝑥(0)=0.5, and the vertical asymptote for the upper solution.

To estimate backward blow-up time for solutions to (4) with “large” initial data 𝑥(𝑡₀)>𝐾_max located above the attractor 𝑥₊(𝑡), we analyze asymptotic behavior of solutions using inequalities (28) rather than (22). Integrating (30) from 𝑡₀ to 𝑡, one arrives at the formula 𝑡_back=𝑡₀+𝐾_maxln𝑐₀⁻¹𝑟_min𝜉₁−𝜉₂, (38) where 𝑐₀=(𝜉(𝑡₀)−𝜉₂)/(𝜉(𝑡₀)−𝜉₁),  𝜘=−𝑟_min(𝜉₂−𝜉₁)/𝐾_max<0, 𝜉_1,2=𝐾_max(1∓1−4𝐻_min/(𝑟_min𝐾_max))/2, and 𝑡_back<𝑡₀ since 𝑐₀⁻¹>1 for “large” initial data. Correspondingly, 𝑥(𝑡)→+∞ as 𝑡→𝑡⁺_back.

Example 9. To find a backward blow-up time for a solution to (26) satisfying 𝑥(0)=16>𝐾_max=14, note that the equilibria of (30) are 𝜉₁=7(1−√35/7)=1.0839 and 𝜉₂=7(1+√35/7)=12.9161. Then, 𝑐₀=0.2067, and a lower bound for a backward blow-up time is estimated as 𝑡_back=−0.4662, see Figure 6. For an upper estimate for a backward blow-up time, note that equilibria for (29) are 𝜂₁=9(1−√39/9)/2=1.3775 and 𝜂₂=9(1+√39/9)/2=7.6225. Then, 𝑐^∗₀=(𝜂(𝑡₀)−𝜂₂)/(𝜂(𝑡₀)−𝜂₁)=0.5729, and an upper bound for a backward blow-up time for a solution to (26) with 𝑥(0)=16 is 𝑡^∗_back=−0.1337, see Figure 7.

Figure 6 

The upper (blue), lower (purple), and solution to (26) (orange), all satisfying the initial condition 𝑥(0)=16; a red dashed vertical asymptote for the upper solution is used for estimating a backward blow-up time.

Figure 7 

Red dashed lines represent two-sided estimates for backward blow-up time for a solution 𝑥(𝑡) (orange) of (26) with 𝑥(0)=16 squeezed between the upper (blue) and lower (purple) solutions.

Remark 10. In a similar manner, one can obtain two-sided estimates for extinction times derived in Section 2.5. Such estimates are useful for the evaluation of a “life span” of a given solution.

Remark 11. Nkashama [18, Theorem 2.1] established that solutions to 𝑑𝑥(𝑡)(𝑡𝑎(𝑡)(𝑡)𝑥(𝑡𝑑𝑡=𝑥)[−𝑏)] (39) located above the attractor 𝑥₊(𝑡) and below the repeller 𝑥₋(𝑡) blow up to +∞ and −∞, respectively, in finite time backward and forward. In our case, there exists a positive repeller for (4) instead of the trivial solution for (39). In addition, there exist positive solutions with small initial data that decay rapidly to −∞, but do not have a vertical asymptote.

3. Saddle-Node Bifurcation

Consider now a 𝑇-periodic differential equation 𝑑𝑦(𝑡)(𝑡)𝑦(𝑡)𝑦(𝑡)𝑑𝑡=𝑟1−𝐾(𝑡)−𝜎𝐻(𝑡), (40) where 𝜎>0 is a parameter that characterizes the intensity of harvesting. An averaged system associated with (40) is 𝑑𝑦(𝑡)=𝑑𝑡(𝑡)𝑦(𝑡)𝑟𝑦1−𝐾−𝜎𝐻. (41) The change of variable 𝑥=𝑦−𝐾/2 reduces (41) to the form 𝑑𝑥(𝑡)𝑑𝑡=𝑎−𝑏𝑥²(𝑡), (42) where 𝑎=(𝑟𝐾)/4−𝜎𝐻 and 𝑏=𝑟/𝐾. Differential equation (42) with real parameters 𝑎,  𝑏 is the simplest canonical example of a saddle-node bifurcation with a nonhyperbolic equilibrium point at the origin. Consequently, one expects that (40) undergoes a so-called nonautonomous saddle-node bifurcation.

We know from Section 2.4 that, for 𝑦∈[0,𝐾_min], (22) holds; for 𝑦<0 or 𝑦>𝐾_max, (28) is satisfied, whereas for 𝑦∈(𝐾_min,𝐾_max), the sign of 𝑟(𝑡)𝑦(1−𝑦/𝐾(𝑡)) cannot be controlled. Let 𝜆(𝑡)=𝑟(𝑡)𝑦(𝑡)(1−𝑦(𝑡)/𝐾(𝑡)). Denote by 𝐼_neg the subset of points 𝑡∈𝐼_𝑇=[𝑡₀,𝑡₀+𝑇] such that 𝜆(𝑡)<0 for 𝑡∈𝐼_neg. Then, for 𝑦∈(0,𝐾_max), one has 𝜆(𝑡)(𝑡)≤−𝜎𝐻0,if𝑡∈𝐼_neg,𝑟_max𝑦(𝑡)𝑦(𝑡)1−𝐾_max−𝜎𝐻_min,if𝑡∈𝐼_𝑇⧵𝐼_neg. (43)

Lemma 12. For 𝜎>𝜎^∗=𝑟_max𝐾_max4𝐻_min, (44) (40) has no periodic solutions; all solutions diverge to −∞.

Proof. Let 𝑦(𝑡) be an arbitrary solution to (40). Taking into account that the function 𝑔(𝑦)=𝑟_max𝑦(1−𝑦/𝐾_max) attains its maximum value 𝑟_max𝐾_max/4 at 𝑦=𝐾_max/2, we have, for all 𝑡∈[𝑡₀,𝑡₀+𝑇], 𝑦𝑡₀+𝑇−𝑦𝑡₀=^𝑡0_𝑡^+𝑇₀𝑟(𝑡)𝑦(𝑡)𝑦(𝑡)1−𝐾(𝑡)(𝑡≤−𝜎𝐻)𝑑𝑡_{𝐼𝑇⧵𝐼neg}𝑟_max𝑦(𝑡𝑦(𝑡))1−𝐾_max−𝜎𝐻_min𝑑𝑡≤^𝑡0_𝑡^+𝑇₀𝑟_max𝐾_max4−𝜎𝐻_min=𝑟𝑑𝑡_max𝐾_max4−𝜎𝐻_min𝑇. (45) If (44) is satisfied, one has 𝑦(𝑡₀+𝑇)<𝑦(𝑡₀), and solution 𝑦(𝑡) cannot be periodic. Similarly, for any 𝑛∈ℕ, (44) yields 𝑦𝑡₀+(𝑛+1)𝑇−𝑦𝑡₀≤𝑟+𝑛𝑇_max𝐾_max4−𝜎𝐻_min𝑇<0. (46) Therefore, for any solution 𝑦(𝑡) to (40), one has 𝑦(𝑡)→−∞ as 𝑡→+∞.

Remark 13. Condition (44) perfectly agrees with the assumption (7) that forces all solutions to decay to −∞.

Theorem 2 assures the existence of the attractor-repeller pair for (40), for all 𝜎<𝜎_∗=𝑟_min𝐾_min4𝐻_max. (47) To explore transition from the attractor-repeller pair to the case with no periodic solutions, we use nullclines and generalized nullclines.

Example 14. Consider a 1-periodic differential equation 𝑑𝑦(𝑡)=((𝑦(𝑡)𝑦(𝑡)𝑑𝑡4+cos2𝜋𝑡))1−()((.8+cos2𝜋𝑡−𝜎2+sin2𝜋𝑡)) (48) Condition (47) holds for all 𝜎<7/4. Equilibria for (16) and (17) are defined, respectively, by 𝜑_1,2=7(1∓7(7−4𝜎)/7)/2 and 𝜓_1,2=9(1∓5(45−4𝜎)/15)/2. Taking in (48) 𝜎=1, we obtain bounds for generalized nullclines 𝜁₁(𝑡) and 𝜁₂(𝑡),  9(1−√205/15)/2<𝜁₁(𝑡)<7(1−√21/7)/2<7(1+√21/7)/2<𝜁₂(𝑡)<9(1+√205/15)/2, see Figure 8.

Figure 8 

Nullclines and constant fences for (48) for 𝜎=1.

Increasing the value of parameter 𝜎 beyond 𝜎_∗, one observes that, at certain point, nullclines for (40) become pinched together, cf. Benardete et al. [11] or Campbell and Kaplan [13]. As the slope takes on negative values in regions between “trapping regions,” some solutions can escape to −∞.

Example 15. Consider a 1-periodic differential equation 𝑑𝑦(𝑡)=((𝑦(𝑡)𝑦(𝑡)𝑑𝑡7−3sin2𝜋𝑡))1−()((.10+2sin2𝜋𝑡−3.8⋅5+2cos2𝜋𝑡)) (49) Observe that 𝑟_min𝐾_min/(4𝐻_max)=8/7<3.8=𝜎. Generalized nullclines for (49) are shown in Figure 9.

Figure 9 

Pinched together pieces of upper (blue) and lower (purple) nullclines and one of solutions to (49) escaping the trapping region for 𝜎=3.8.

Theorem 16. There exists a bifurcation value 𝜎_𝑏𝑖𝑓 satisfying the inequality 𝑟_min𝐾_min/(4𝐻_max)≤𝜎_𝑏𝑖𝑓≤𝑟_max𝐾_max/(4𝐻_min) such that (40) has (i) no periodic solutions if 𝜎>𝜎_𝑏𝑖𝑓;  (ii) exactly one periodic solution if 𝜎=𝜎_𝑏𝑖𝑓;  (iii) two periodic solutions if 𝜎<𝜎_𝑏𝑖𝑓 .

Proof. Define the Poincaré map ℎ(𝑦₀,𝜎) for (40) on (0,𝑇). Equation ℎ(𝑦₀,𝜎)=𝑦₀ is quadratic, see, for instance, [11, page 211]. Since the second derivative ℎ_𝑦𝑦(𝑦₀,𝜎) is negative, the graph of ℎ(𝑦₀,𝜎) is concave down. For 𝜎<𝑟_min𝐾_min/4𝐻_max, it crosses the line 𝑔(𝑦₀)=𝑦₀ at two points corresponding to two periodic solutions of (40). On the other hand, the graph has no intersections with this line for 𝜎>𝑟_max𝐾_max/4𝐻_min, and thus, in this case there are no periodic solutions to (40). The derivative ℎ_𝜎(𝑦₀,𝜎) is negative. Hence, for each fixed 𝑦₀, the function 𝛿(𝜎)=ℎ(𝑦₀,𝜎) decreases as 𝜎 increases from 𝑟_min𝐾_min/4𝐻_max to 𝑟_max𝐾_max/4𝐻_min, and there exists a unique value 𝜎=𝜎_bif  such that the graph of ℎ(𝑦₀,𝜎) is tangent to the line 𝑔(𝑦₀)=𝑦₀.

Remark 17. We stress that Benardete et al. [11, Section 5, pages 212–215] established bounds for the bifurcation parameter 𝜎_bif for a much simpler differential equation (2) using extensively its symmetry about the point (1/4,1/2), whereas Theorem 16 requires no symmetry properties of (40) at all.

Example 18. Consider a 1-periodic differential equation 𝑑𝑦(𝑡)=((𝑦(𝑡)(𝑡)𝑑𝑡5+sin2𝜋𝑡))1−2𝑦()−𝜎23+5cos2𝜋𝑡2(11−3cos(2𝜋𝑡)). (50) For 𝜎=1, (50) turns into (26) for which the attractor-repeller pair exists, see Example 4. One has 𝑟_min𝐾_min/(4𝐻_max)=9/7 and 𝑟_max𝐾_max/(4𝐻_min)=21/4. By Theorem 16, there exists a bifurcation value 𝜎_bif∈[9/7,21/4]. Numerical experiments can be used to approximate the value of 𝜎_bif. In fact, for 𝜎=2.395, (50) has two periodic solutions plotted in Figure 10. On the other hand, for 𝜎=2.396, (50) does not have periodic solutions anymore, see Figure 11. Consequently, 𝜎_bif∈(2.395,2.396).

Figure 10 

Several solutions to (50) for 𝜎=2.395.

Figure 11 

Several solutions to (50) for 𝜎=2.396.

4. Conclusions

Differential equation (4) exhibits quite interesting dynamics. Existence of the attractor-repeller pair is assured by condition (10); efficient estimates for periodic solutions are derived in Section 2.4. In the presence of the positive attractor-repeller pair, all other solutions to (4) fall into one of the three groups. Namely, as time 𝑡 advances, (i) solutions with small initial data 𝑥(𝑡₀)∈(0,𝑥₋(𝑡₀)) move away from 𝑥₋(𝑡), decay rapidly to −∞, and blow up in the future as 𝑡→𝑡⁻_forw (approach the repeller as 𝑡→−∞); (ii) solutions with initial data 𝑥(𝑡₀)∈(𝑥₋(𝑡₀),𝑥₊(𝑡₀)) leave the vicinity of the repeller 𝑥₋(𝑡) and approach the attractor 𝑥₊(𝑡) (correspondingly, leave the vicinity of the attractor 𝑥₊(𝑡) and approach the repeller 𝑥₋(𝑡) as 𝑡→−∞); these are so-called heteroclinic orbits; (iii) solutions with initial data 𝑥(𝑡₀)>𝑥₊(𝑡₀) approach the attractor 𝑥₊(𝑡), and blow up back in time as 𝑡→𝑡⁺_back. If condition (7) holds, (4) has no periodic solutions. All solutions decay to −∞; extinction time is estimated for solutions with small and large initial data.

Qualitative properties of (40) vary with the parameter 𝜎. As 𝜎 increases, the number of periodic solutions changes from two to zero; differential equation (40) undergoes a nonautonomous saddle-node bifurcation. Estimates for the bifurcation parameter 𝜎_bif can be significantly tightened by using numerical methods. Contrary to [11], our bifurcation analysis does not require symmetry of (40). We conclude by noting that although all numerical experiments in this paper were performed using Wolfram Mathematica 7.0, any computer algebra system can be used instead.

Acknowledgments

This research started when S. P. Rogovchenko and Yu. V. Rogovchenko visited the Abdus Salam International Centre for Theoretical Physics, Trieste, Italy, whose warm hospitality, excellent research facilities, and financial support are gratefully acknowledged. Yu. V. Rogovchenko also acknowledges the research grant from the Faculty of Science and Technology of Umeå University. The authors thank an anonymous referee for useful remarks that helped us to accentuate the importance of the study undertaken in this paper.