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

Authors: Mustafa Hasanbulli; Svitlana P. Rogovchenko; Yuriy V. Rogovchenko
Journal: Journal of Applied Mathematics
 (2013)

---

## 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.


---

## Body

## 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)dx(t)dt=r(t)x(t)(1-x(t)K(t)),
where the intrinsic growth rate r and the carrying capacity of the environment K are positive, continuous functions that vary periodically with time, r(t+T)=r(t) and K(t+T)=K(t), for all t∈ℝ. 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 r¯=def(1/T)∫0T‍r(t)dt>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 r and K in (1) are periodic along with the harvesting term H has been studied only by Brauer and Sánchez [12].A bifurcation problem for a differential equation(2)dx(t)dt=kx(t)(1-x(t))-H(1+sin(2πt))
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 x-intercepts in time.” Recently, using Crandall-Rabinowitz saddle-node bifurcation theorem, Liu et al. [17] established existence of periodic solutions for dx(t)/dt=f(x(t))-σh(t). 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 instantt*. 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)dx(t)dt=f(t,x(t))-h(t,x(t)),
where the function f describes the growth of unharvested population and the function h 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 h is a linear function of population size, h(t,x(t))=H0x(t), 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 E, 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,h(t,x(t))=H0, is called constant rate or constant yieldharvesting. 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 H0 is a function of time, H0=H0(t); 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 H(t) 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(4)dx(t)dt=r(t)x(t)(1-x(t)K(t))-H(t),
where H(t) 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),
(5)dx(t)dt=kx(t)(1-x(t)C)-a(1+sin(bt)),
has been used for determination of an optimal harvesting frequency of fishing cycles b 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,(6)dx(t)dt=r(t)x(t)(1-x(t)K(t))-H(t)x(t),
where H 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 r1(t)=r(t)-H(t) and carrying capacity K1(t)=(K(t)(r(t)-H(t)))/r(t). 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 T-periodic solution xp(t) to (6) which is asymptotically stable with the domain of attraction containing positive initial data provided that the time average r¯1 of the function r1(t) is positive, r¯1>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 r,  K, or H 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 xtriv(t)≡0, unless H(t)≡0. This distinguishes dynamics of solutions with small positive initial data. Although any two solutions of (4) with 0<x1(t0)<x2(t0) satisfy x1(t)<x2(t), for all t≥t0, for any solution xsmall(t) of this equation with 0<xsmall(t0)<ε, for a small ε>0, there exists a text, termed the “extinction time,” such that xsmall(text)=0 and xsmall(t)<0, for all t>text, see Section 2.5. This cannot happen for (6), for which the unstable solution xtriv(t) 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 x, 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, T-periodic right-hand side. As Brauer and Sánchez [12, page 234] pointed out, this T-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 notationgmin=min0≤t≤Tg(t) and gmax=max0≤t≤Tg(t) and assume that at least one of the inequalities rmin≤rmax,Kmin≤Kmax,Hmin≤Hmax 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 dx(t)/dt is negative for all t∈ℝ provided that
(7)γ(t)=defH(t)-r(t)K(t)4>0,
in which case x(t)→-∞ as t→+∞. 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
(8)γ(t)<0
holds for all t∈ℝ, they argued that, for all t,x′|x=0<0,x′|x=Kmax<0 and x′|x=Kmax/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 r(t)=4+sin(2πt),K(t)=5+2sin(2πt) and H(t)=3+sin(2πt). Then, for all t∈ℝ,γ(t)≤-1/4, and (8) is satisfied, but the derivative of solution of the given differential equation,
(9)x′|x=Kmax/2=7(4+sin(2πt))2×(1-72(5+2sin(2πt)))-(3+sin(2πt)),
changes sign infinitely many times on ℝ.Theorem 2.
Assume that condition(10)Hmax<rminKmin4
is satisfied. Then (4) has two positive periodic solutions, x+(t) and x-(t), which are the forward and pullback attractors, respectively, and, for all t∈ℝ,
(11)0<x-(t)<Kmin2<x+(t)<Kmax.Proof.
Observe first that [40, Theorem on page 30] yields that (4) cannot have more than two periodic solutions. Furthermore, for all t∈ℝ,
(12)dx(t)dt|x(t)=Kmax=r(t)Kmax(1-KmaxK(t))-H(t)≤r(t)Kmax(1-KmaxKmax)-H(t)=-H(t)<0.
On the other hand, by virtue of (10),
(13)dx(t)dt|x(t)=Kmin/2=r(t)Kmin2(1-Kmin2K(t))-H(t)≥rminKmin2(1-Kmin2Kmin)-Hmax=rminKmin4-Hmax>0.
By [40, Theorem on page 34], there exists a periodic solution x+(t) of (4) satisfying Kmin/2<x+(t)<Kmax. Direction field and uniqueness arguments imply that this solution is a forward attractor for all solutions of (4) with initial data satisfying Kmin/2<x(t0)<Kmax.
Keeping in mind that a pullback attractor is a forward repeller, see Rasmussen [43, Remark 3.4, page 272], introduce a new variable τ=-t. If x(t) is a solution to (4), then x~(τ)=x(-t) satisfies
(14)dx~(τ)dτ=-r(τ)x~(τ)(1-x~(τ)K(τ))+H(τ).
Note first that, for x~(τ)=0, the slope dx~(τ)/dτ is positive for all τ∈ℝ, dx~(τ)/dτ|x~(τ)=0=H(τ)>0. Furthermore, by virtue of (10),
(15)dx~(τ)dτ|x~(τ)=Kmin/2=-r(τ)Kmin2(1-Kmin2K(τ))+H(τ)≤-rminKmin2(1-Kmin2Kmin)+Hmax=Hmax-rminKmin4<0.
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 x-(t) 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 i=1,2, let φ¯i and ψ¯i be equilibrium solutions to autonomous differential equations
(16)dφ(t)dt=rminφ(t)(1-φ(t)Kmin)-Hmax,(17)dψ(t)dt=rmaxψ(t)(1-ψ(t)Kmax)-Hmin.Theorem 3.
Assume that (10) holds for all t∈ℝ. Suppose also that
(18)2HminHmax≤rmaxrmin.
Then, for the attractor-repeller pair x+(t),x-(t) in Theorem 2, one has
(19)0<ψ¯1<x-(t)<φ¯1,φ¯2<x+(t)<Kmax.
If, in addition,
(20)Kmax<2Kmin,
then
(21)φ¯2<x+(t)<ψ¯2<Kmax.Proof.
By (11), it suffices to consider only the values of x between 0 and Kmax. However, one can completely control the behavior of the right-hand side of (4) only for 0≤x≤Kmin, in which case the following estimates hold:
(22)rminx(t)(1-x(t)Kmin)-Hmax≤r(t)x(t)(1-x(t)K(t))-H(t)≤rmaxx(t)(1-x(t)Kmax)-Hmin.
Equilibrium solutions to differential equations (16) and (17) are given by
(23)φ¯1,2=Kmin2(1∓1-4HmaxrminKmin),ψ¯1,2=Kmax2(1∓1-4HminrmaxKmax).
Condition (18) ensures that these four equilibria are ordered as follows:
(24)ψ¯1<φ¯1<φ¯2<ψ¯2.
Indeed, note first that φ¯1<φ¯2<Kmin. Then, by virtue of Hmin/(rmaxKmax)≤Hmin/(rminKmin)≤Hmax/(rminKmin)<1/4, we conclude that φ¯2<ψ¯2. To demonstrate that ψ¯1<φ¯1, we keep Hmax,Hmin,rmax, and rmin fixed and analyze the behavior of φ¯1(Kmin) and ψ¯1(Kmax) described by the function χ(x)=x(1-1-ν/x)/2,x≥ν. One can verify that χ is strictly decreasing on [ν,+∞), attains its maximal value ν/2 at x=ν, and decays to ν/4 as x→+∞. Correspondingly, φ¯1(Kmin)∈(Hmax/rmin,2Hmax/rmin] and ψ¯1(Kmax)∈(Hmin/rmax,2Hmin/rmax]. Finally, condition (18) guarantees that maxKmaxψ¯1(Kmax)=2Hmin/rmax<Hmax/rmin=infKminφ¯1(Kmin), which completes the proof of (24).
Consider now the functionsα1(t)=φ¯1 and β1(t)=ψ¯1. One has
(25)f(t,α1(t))=r(t)φ¯1(1-φ¯1K(t))-H(t)≥rminφ¯1(1-φ¯1Kmin)-Hmax=α1′(t),f(t,β1(t))=r(t)ψ¯1(1-ψ¯1K(t))-H(t)≤rmaxψ¯1(1-ψ¯1Kmax)-Hmin=β1′(t).
Hence, α1(t)=φ¯1 and β1(t)=ψ¯1 are lower and upper fences, whereas the horizontal strip bounded by φ¯1 and ψ¯1 is an antifunnel, and there exists a periodic solution to (4) located between α1(t) and β1(t). For the repeller x-(t), one has, for all t∈ℝ,ψ¯1<x-(t)<φ¯1. If (20) is satisfied, one can show that both estimates for the attractor x+(t) 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(26)dx(t)dt=(5+sin(2πt))x(t)(1-2x(t)23+5cos(2πt))-12(11-3cos(2πt)).
In this case,
(27)φ¯1,2=32(3∓2),ψ¯1,2=7(1∓35721).
Condition (18) is satisfied because 2Hmin/Hmax=2·4/7≤6/4=rmax/rmin. By (19), 7(1-357/21)<x-(t)<3(3-2)/2<3(3+2)/2<x+(t)<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 φ¯2-φ¯1=Kmin1-4Hmax/(rminKmin); it equals 32 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 x+(t) and x-(t).Figure 1
Tighter estimates for both periodic solutions of (26).Remark 5.
Exact solutionsφ(t) and ψ(t) to (16) and (17) give tighter estimates for the solution x(t) of (4) satisfying the same initial condition, provided that 0≤x(t0)=φ(t0)=ψ(t0)≤Kmin or (20) holds. With the growth of t, both φ(t) and ψ(t) approach equilibria quite fast; even for relatively small values of t, 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 t→∓∞.Figure 2
Unstable and stable periodic solutions embraced by upper and lower solutions.Remark 6.
To describe behavior of solutionsx(t) with initial data satisfying x(t0)<0 or x(t0)>Kmax, one needs estimates
(28)rmaxx(t)(1-x(t)Kmin)-Hmax≤r(t)x(t)(1-x(t)K(t))-H(t)≤rminx(t)(1-x(t)Kmax)-Hmin
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
(29)dη(t)dt=rmaxη(t)(1-η(t)Kmin)-Hmax,(30)dξ(t)dt=rminξ(t)(1-ξ(t)Kmax)-Hmin,
is used instead of (16) and (17).


### 2.5. Extinction Times

We start with an upper bound for the extinction timetext for solutions with initial data x(t0)∈(0,Kmin]. By (22), one has to study the behavior of solutions to differential equation (17). Suppose that
(31)Hmin>Hcrit=rmaxKmax4.
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,
(32)texts=t0+2rmax-1(HminHcrit-1)-1/2×{(HminHcrit-1)-12tan-1(HminHcrit-1)-1/2+tan-1((HminHcrit-1)-122Kmax(ψ(t0)-Kmax2)×(HminHcrit-1)-1/2(HminHcrit-1)-12)}.
In particular, if (20) holds, the extinction time for the solution x*(t) of (4) satisfying x*(t0)=Kmax/2 is determined by a simpler expression,
(33)texts=t0+2rmax-1(HminHcrit-1)-1/2tan-1(HminHcrit-1)-1/2.For solutions with initial datax(t0)>Kmax, the situation is similar. Assume that (31) holds. By virtue of (28), differential equation (30) is used. Then, Hmax≥Hmin>Hcrit=rmaxKmax/4≥rminKmax/4=H*. Integrating (30), letting ξ(t)=0 and solving the resulting equation for textl, one has
(34)textl=t0+2rmin-1(HminH*-1)-1/2×{(HminHcrit-1)-12tan-1(HminH*-1)-1/2+tan-1((HminHcrit-1)-122Kmax(ξ(t0)-Kmax2)×(HminH*-1)-1/2(HminHcrit-1)-12)}.Example 7.
Consider a 1-periodic differential equation(35)dx(t)dt=exp(1+cos(2πt)2)×x(t)(1-x(t)6+cosh(cos(2πt)))-(7+sinh(sin(2πt))).
Condition (31) holds since 5.8248=Hmin>rmaxKmax/4=5.1260. Consequently, all solutions of (35) decay to -∞. By (32), the extinction time for the solution satisfying x(0)=1<Kmin is texts=0.2229, see Figure 3. For a solution of (35) satisfying x(0)=8>Kmax=6+cosh1=7.5431, one has 8.1752=Hmax>rminKmax/4=1.8857. By (34), an upper bound for the extinction time is textl=1.7506, see Figure 4.Figure 3
The upper solution (blue) and solution (orange) to (35) satisfying the initial condition x(0)=1. An upper bound for the extinction time for this solution is texts=0.2229.Figure 4
The upper solution (blue) and solution (orange) to (35) satisfying the initial condition x(0)=8. An upper bound for the extinction time for this solution is textl=1.7506.


### 2.6. Forward and Backward Blow-Up Time

In this section, we assume that (10) holds for all t∈ℝ. First, we estimate forward blow-up time for “small” solutions with initial data x(t0)∈(0,ψ¯1), that is, for solutions located below the repeller x-(t). These solutions decay rapidly to -∞ and have vertical asymptotes to the right of t=t0. Taking into account estimates (22) that hold for small values of x(t0), we integrate differential equation (17) from t0 to t obtaining
(36)ψ(t)=ψ¯1c0eϰ(t-t0)-ψ¯2c0eϰ(t-t0)-1,
where c0=(ψ(t0)-ψ¯2)/(ψ(t0)-ψ¯1),  ϰ=-rmax(ψ¯2-ψ¯1)/Kmax<0, and equilibrium solutions ψ¯1 and ψ¯2 are defined above. By (36), solutions with small initial data 0<x(t0)<ψ¯1 blow up in the future, x(t)→-∞ as t→tforw-; an estimated value for a forward blow-up time tforw is given by
(37)tforw=t0+Kmaxln(c0-1)rmax(ψ¯1-ψ¯2).
Equation t=tforw defines a vertical asymptote for solutions to (17); tforw>t0 because, for x(t0)<ψ(t0)<ψ¯1, one always has c0-1<1.Example 8.
To estimate a forward blow-up time for the solution to (26) satisfying the initial condition x(0)=0.5, note that equilibria of differential equation (17) associated with (26) are given by (27), c0=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 tforw=0.7686, see Figure 5.Figure 5
The upper solution (blue), lower (purple), and solution (orange) to (26), all satisfying the initial condition, x(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 x(t0)>Kmax located above the attractor x+(t), we analyze asymptotic behavior of solutions using inequalities (28) rather than (22). Integrating (30) from t0 to t, one arrives at the formula
(38)tback=t0+Kmaxln(c0-1)rmin(ξ¯1-ξ¯2),
where c0=(ξ(t0)-ξ¯2)/(ξ(t0)-ξ¯1),  ϰ=-rmin(ξ¯2-ξ¯1)/Kmax<0, ξ¯1,2=Kmax(1∓1-4Hmin/(rminKmax))/2, and tback<t0 since c0-1>1 for “large” initial data. Correspondingly, x(t)→+∞ as t→tback+.Example 9.
To find a backward blow-up time for a solution to (26) satisfying x(0)=16>Kmax=14, note that the equilibria of (30) are ξ¯1=7(1-35/7)=1.0839 and ξ¯2=7(1+35/7)=12.9161. Then, c0=0.2067, and a lower bound for a backward blow-up time is estimated as tback=-0.4662, see Figure 6. For an upper estimate for a backward blow-up time, note that equilibria for (29) are η¯1=9(1-39/9)/2=1.3775 and η¯2=9(1+39/9)/2=7.6225. Then, c0*=(η(t0)-η¯2)/(η(t0)-η¯1)=0.5729, and an upper bound for a backward blow-up time for a solution to (26) with x(0)=16 is tback*=-0.1337, see Figure 7.Figure 6
The upper (blue), lower (purple), and solution to (26) (orange), all satisfying the initial condition x(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 solutionx(t) (orange) of (26) with x(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 Section2.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)dx(t)dt=x(t)[a(t)-b(t)x(t)]
located above the attractor x+(t) and below the repeller x-(t) 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.


## 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)dx(t)dt=f(t,x(t))-h(t,x(t)),
where the function f describes the growth of unharvested population and the function h 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 h is a linear function of population size, h(t,x(t))=H0x(t), 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 E, 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,h(t,x(t))=H0, is called constant rate or constant yieldharvesting. 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 H0 is a function of time, H0=H0(t); 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 H(t) 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(4)dx(t)dt=r(t)x(t)(1-x(t)K(t))-H(t),
where H(t) 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),
(5)dx(t)dt=kx(t)(1-x(t)C)-a(1+sin(bt)),
has been used for determination of an optimal harvesting frequency of fishing cycles b 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,(6)dx(t)dt=r(t)x(t)(1-x(t)K(t))-H(t)x(t),
where H 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 r1(t)=r(t)-H(t) and carrying capacity K1(t)=(K(t)(r(t)-H(t)))/r(t). 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 T-periodic solution xp(t) to (6) which is asymptotically stable with the domain of attraction containing positive initial data provided that the time average r¯1 of the function r1(t) is positive, r¯1>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 r,  K, or H 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 xtriv(t)≡0, unless H(t)≡0. This distinguishes dynamics of solutions with small positive initial data. Although any two solutions of (4) with 0<x1(t0)<x2(t0) satisfy x1(t)<x2(t), for all t≥t0, for any solution xsmall(t) of this equation with 0<xsmall(t0)<ε, for a small ε>0, there exists a text, termed the “extinction time,” such that xsmall(text)=0 and xsmall(t)<0, for all t>text, see Section 2.5. This cannot happen for (6), for which the unstable solution xtriv(t) 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 x, 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, T-periodic right-hand side. As Brauer and Sánchez [12, page 234] pointed out, this T-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 notationgmin=min0≤t≤Tg(t) and gmax=max0≤t≤Tg(t) and assume that at least one of the inequalities rmin≤rmax,Kmin≤Kmax,Hmin≤Hmax 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 dx(t)/dt is negative for all t∈ℝ provided that
(7)γ(t)=defH(t)-r(t)K(t)4>0,
in which case x(t)→-∞ as t→+∞. 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
(8)γ(t)<0
holds for all t∈ℝ, they argued that, for all t,x′|x=0<0,x′|x=Kmax<0 and x′|x=Kmax/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 r(t)=4+sin(2πt),K(t)=5+2sin(2πt) and H(t)=3+sin(2πt). Then, for all t∈ℝ,γ(t)≤-1/4, and (8) is satisfied, but the derivative of solution of the given differential equation,
(9)x′|x=Kmax/2=7(4+sin(2πt))2×(1-72(5+2sin(2πt)))-(3+sin(2πt)),
changes sign infinitely many times on ℝ.Theorem 2.
Assume that condition(10)Hmax<rminKmin4
is satisfied. Then (4) has two positive periodic solutions, x+(t) and x-(t), which are the forward and pullback attractors, respectively, and, for all t∈ℝ,
(11)0<x-(t)<Kmin2<x+(t)<Kmax.Proof.
Observe first that [40, Theorem on page 30] yields that (4) cannot have more than two periodic solutions. Furthermore, for all t∈ℝ,
(12)dx(t)dt|x(t)=Kmax=r(t)Kmax(1-KmaxK(t))-H(t)≤r(t)Kmax(1-KmaxKmax)-H(t)=-H(t)<0.
On the other hand, by virtue of (10),
(13)dx(t)dt|x(t)=Kmin/2=r(t)Kmin2(1-Kmin2K(t))-H(t)≥rminKmin2(1-Kmin2Kmin)-Hmax=rminKmin4-Hmax>0.
By [40, Theorem on page 34], there exists a periodic solution x+(t) of (4) satisfying Kmin/2<x+(t)<Kmax. Direction field and uniqueness arguments imply that this solution is a forward attractor for all solutions of (4) with initial data satisfying Kmin/2<x(t0)<Kmax.
Keeping in mind that a pullback attractor is a forward repeller, see Rasmussen [43, Remark 3.4, page 272], introduce a new variable τ=-t. If x(t) is a solution to (4), then x~(τ)=x(-t) satisfies
(14)dx~(τ)dτ=-r(τ)x~(τ)(1-x~(τ)K(τ))+H(τ).
Note first that, for x~(τ)=0, the slope dx~(τ)/dτ is positive for all τ∈ℝ, dx~(τ)/dτ|x~(τ)=0=H(τ)>0. Furthermore, by virtue of (10),
(15)dx~(τ)dτ|x~(τ)=Kmin/2=-r(τ)Kmin2(1-Kmin2K(τ))+H(τ)≤-rminKmin2(1-Kmin2Kmin)+Hmax=Hmax-rminKmin4<0.
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 x-(t) 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 i=1,2, let φ¯i and ψ¯i be equilibrium solutions to autonomous differential equations
(16)dφ(t)dt=rminφ(t)(1-φ(t)Kmin)-Hmax,(17)dψ(t)dt=rmaxψ(t)(1-ψ(t)Kmax)-Hmin.Theorem 3.
Assume that (10) holds for all t∈ℝ. Suppose also that
(18)2HminHmax≤rmaxrmin.
Then, for the attractor-repeller pair x+(t),x-(t) in Theorem 2, one has
(19)0<ψ¯1<x-(t)<φ¯1,φ¯2<x+(t)<Kmax.
If, in addition,
(20)Kmax<2Kmin,
then
(21)φ¯2<x+(t)<ψ¯2<Kmax.Proof.
By (11), it suffices to consider only the values of x between 0 and Kmax. However, one can completely control the behavior of the right-hand side of (4) only for 0≤x≤Kmin, in which case the following estimates hold:
(22)rminx(t)(1-x(t)Kmin)-Hmax≤r(t)x(t)(1-x(t)K(t))-H(t)≤rmaxx(t)(1-x(t)Kmax)-Hmin.
Equilibrium solutions to differential equations (16) and (17) are given by
(23)φ¯1,2=Kmin2(1∓1-4HmaxrminKmin),ψ¯1,2=Kmax2(1∓1-4HminrmaxKmax).
Condition (18) ensures that these four equilibria are ordered as follows:
(24)ψ¯1<φ¯1<φ¯2<ψ¯2.
Indeed, note first that φ¯1<φ¯2<Kmin. Then, by virtue of Hmin/(rmaxKmax)≤Hmin/(rminKmin)≤Hmax/(rminKmin)<1/4, we conclude that φ¯2<ψ¯2. To demonstrate that ψ¯1<φ¯1, we keep Hmax,Hmin,rmax, and rmin fixed and analyze the behavior of φ¯1(Kmin) and ψ¯1(Kmax) described by the function χ(x)=x(1-1-ν/x)/2,x≥ν. One can verify that χ is strictly decreasing on [ν,+∞), attains its maximal value ν/2 at x=ν, and decays to ν/4 as x→+∞. Correspondingly, φ¯1(Kmin)∈(Hmax/rmin,2Hmax/rmin] and ψ¯1(Kmax)∈(Hmin/rmax,2Hmin/rmax]. Finally, condition (18) guarantees that maxKmaxψ¯1(Kmax)=2Hmin/rmax<Hmax/rmin=infKminφ¯1(Kmin), which completes the proof of (24).
Consider now the functionsα1(t)=φ¯1 and β1(t)=ψ¯1. One has
(25)f(t,α1(t))=r(t)φ¯1(1-φ¯1K(t))-H(t)≥rminφ¯1(1-φ¯1Kmin)-Hmax=α1′(t),f(t,β1(t))=r(t)ψ¯1(1-ψ¯1K(t))-H(t)≤rmaxψ¯1(1-ψ¯1Kmax)-Hmin=β1′(t).
Hence, α1(t)=φ¯1 and β1(t)=ψ¯1 are lower and upper fences, whereas the horizontal strip bounded by φ¯1 and ψ¯1 is an antifunnel, and there exists a periodic solution to (4) located between α1(t) and β1(t). For the repeller x-(t), one has, for all t∈ℝ,ψ¯1<x-(t)<φ¯1. If (20) is satisfied, one can show that both estimates for the attractor x+(t) 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(26)dx(t)dt=(5+sin(2πt))x(t)(1-2x(t)23+5cos(2πt))-12(11-3cos(2πt)).
In this case,
(27)φ¯1,2=32(3∓2),ψ¯1,2=7(1∓35721).
Condition (18) is satisfied because 2Hmin/Hmax=2·4/7≤6/4=rmax/rmin. By (19), 7(1-357/21)<x-(t)<3(3-2)/2<3(3+2)/2<x+(t)<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 φ¯2-φ¯1=Kmin1-4Hmax/(rminKmin); it equals 32 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 x+(t) and x-(t).Figure 1
Tighter estimates for both periodic solutions of (26).Remark 5.
Exact solutionsφ(t) and ψ(t) to (16) and (17) give tighter estimates for the solution x(t) of (4) satisfying the same initial condition, provided that 0≤x(t0)=φ(t0)=ψ(t0)≤Kmin or (20) holds. With the growth of t, both φ(t) and ψ(t) approach equilibria quite fast; even for relatively small values of t, 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 t→∓∞.Figure 2
Unstable and stable periodic solutions embraced by upper and lower solutions.Remark 6.
To describe behavior of solutionsx(t) with initial data satisfying x(t0)<0 or x(t0)>Kmax, one needs estimates
(28)rmaxx(t)(1-x(t)Kmin)-Hmax≤r(t)x(t)(1-x(t)K(t))-H(t)≤rminx(t)(1-x(t)Kmax)-Hmin
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
(29)dη(t)dt=rmaxη(t)(1-η(t)Kmin)-Hmax,(30)dξ(t)dt=rminξ(t)(1-ξ(t)Kmax)-Hmin,
is used instead of (16) and (17).


## 2.5. Extinction Times

We start with an upper bound for the extinction timetext for solutions with initial data x(t0)∈(0,Kmin]. By (22), one has to study the behavior of solutions to differential equation (17). Suppose that
(31)Hmin>Hcrit=rmaxKmax4.
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,
(32)texts=t0+2rmax-1(HminHcrit-1)-1/2×{(HminHcrit-1)-12tan-1(HminHcrit-1)-1/2+tan-1((HminHcrit-1)-122Kmax(ψ(t0)-Kmax2)×(HminHcrit-1)-1/2(HminHcrit-1)-12)}.
In particular, if (20) holds, the extinction time for the solution x*(t) of (4) satisfying x*(t0)=Kmax/2 is determined by a simpler expression,
(33)texts=t0+2rmax-1(HminHcrit-1)-1/2tan-1(HminHcrit-1)-1/2.For solutions with initial datax(t0)>Kmax, the situation is similar. Assume that (31) holds. By virtue of (28), differential equation (30) is used. Then, Hmax≥Hmin>Hcrit=rmaxKmax/4≥rminKmax/4=H*. Integrating (30), letting ξ(t)=0 and solving the resulting equation for textl, one has
(34)textl=t0+2rmin-1(HminH*-1)-1/2×{(HminHcrit-1)-12tan-1(HminH*-1)-1/2+tan-1((HminHcrit-1)-122Kmax(ξ(t0)-Kmax2)×(HminH*-1)-1/2(HminHcrit-1)-12)}.Example 7.
Consider a 1-periodic differential equation(35)dx(t)dt=exp(1+cos(2πt)2)×x(t)(1-x(t)6+cosh(cos(2πt)))-(7+sinh(sin(2πt))).
Condition (31) holds since 5.8248=Hmin>rmaxKmax/4=5.1260. Consequently, all solutions of (35) decay to -∞. By (32), the extinction time for the solution satisfying x(0)=1<Kmin is texts=0.2229, see Figure 3. For a solution of (35) satisfying x(0)=8>Kmax=6+cosh1=7.5431, one has 8.1752=Hmax>rminKmax/4=1.8857. By (34), an upper bound for the extinction time is textl=1.7506, see Figure 4.Figure 3
The upper solution (blue) and solution (orange) to (35) satisfying the initial condition x(0)=1. An upper bound for the extinction time for this solution is texts=0.2229.Figure 4
The upper solution (blue) and solution (orange) to (35) satisfying the initial condition x(0)=8. An upper bound for the extinction time for this solution is textl=1.7506.


## 2.6. Forward and Backward Blow-Up Time

In this section, we assume that (10) holds for all t∈ℝ. First, we estimate forward blow-up time for “small” solutions with initial data x(t0)∈(0,ψ¯1), that is, for solutions located below the repeller x-(t). These solutions decay rapidly to -∞ and have vertical asymptotes to the right of t=t0. Taking into account estimates (22) that hold for small values of x(t0), we integrate differential equation (17) from t0 to t obtaining
(36)ψ(t)=ψ¯1c0eϰ(t-t0)-ψ¯2c0eϰ(t-t0)-1,
where c0=(ψ(t0)-ψ¯2)/(ψ(t0)-ψ¯1),  ϰ=-rmax(ψ¯2-ψ¯1)/Kmax<0, and equilibrium solutions ψ¯1 and ψ¯2 are defined above. By (36), solutions with small initial data 0<x(t0)<ψ¯1 blow up in the future, x(t)→-∞ as t→tforw-; an estimated value for a forward blow-up time tforw is given by
(37)tforw=t0+Kmaxln(c0-1)rmax(ψ¯1-ψ¯2).
Equation t=tforw defines a vertical asymptote for solutions to (17); tforw>t0 because, for x(t0)<ψ(t0)<ψ¯1, one always has c0-1<1.Example 8.
To estimate a forward blow-up time for the solution to (26) satisfying the initial condition x(0)=0.5, note that equilibria of differential equation (17) associated with (26) are given by (27), c0=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 tforw=0.7686, see Figure 5.Figure 5
The upper solution (blue), lower (purple), and solution (orange) to (26), all satisfying the initial condition, x(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 x(t0)>Kmax located above the attractor x+(t), we analyze asymptotic behavior of solutions using inequalities (28) rather than (22). Integrating (30) from t0 to t, one arrives at the formula
(38)tback=t0+Kmaxln(c0-1)rmin(ξ¯1-ξ¯2),
where c0=(ξ(t0)-ξ¯2)/(ξ(t0)-ξ¯1),  ϰ=-rmin(ξ¯2-ξ¯1)/Kmax<0, ξ¯1,2=Kmax(1∓1-4Hmin/(rminKmax))/2, and tback<t0 since c0-1>1 for “large” initial data. Correspondingly, x(t)→+∞ as t→tback+.Example 9.
To find a backward blow-up time for a solution to (26) satisfying x(0)=16>Kmax=14, note that the equilibria of (30) are ξ¯1=7(1-35/7)=1.0839 and ξ¯2=7(1+35/7)=12.9161. Then, c0=0.2067, and a lower bound for a backward blow-up time is estimated as tback=-0.4662, see Figure 6. For an upper estimate for a backward blow-up time, note that equilibria for (29) are η¯1=9(1-39/9)/2=1.3775 and η¯2=9(1+39/9)/2=7.6225. Then, c0*=(η(t0)-η¯2)/(η(t0)-η¯1)=0.5729, and an upper bound for a backward blow-up time for a solution to (26) with x(0)=16 is tback*=-0.1337, see Figure 7.Figure 6
The upper (blue), lower (purple), and solution to (26) (orange), all satisfying the initial condition x(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 solutionx(t) (orange) of (26) with x(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 Section2.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)dx(t)dt=x(t)[a(t)-b(t)x(t)]
located above the attractor x+(t) and below the repeller x-(t) 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 aT-periodic differential equation
(40)dy(t)dt=r(t)y(t)(1-y(t)K(t))-σH(t),
where σ>0 is a parameter that characterizes the intensity of harvesting. An averaged system associated with (40) is
(41)dy(t)dt=r¯y(t)(1-y(t)K¯)-σH¯.
The change of variable x=y-K¯/2 reduces (41) to the form
(42)dx(t)dt=a-bx2(t),
where a=(r¯K¯)/4-σH¯ and b=r¯/K¯. Differential equation (42) with real parameters a,  b 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 Section2.4 that, for y∈[0,Kmin], (22) holds; for y<0 or y>Kmax, (28) is satisfied, whereas for y∈(Kmin,Kmax), the sign of r(t)y(1-y/K(t)) cannot be controlled. Let λ(t)=r(t)y(t)(1-y(t)/K(t)). Denote by Ineg the subset of points t∈IT=[t0,t0+T] such that λ(t)<0 for t∈Ineg. Then, for y∈(0,Kmax), one has
(43)λ(t)-σH(t)≤{0,ift∈Ineg,rmaxy(t)(1-y(t)Kmax)-σHmin,ift∈IT∖Ineg.Lemma 12.
For(44)σ>σ*=rmaxKmax4Hmin,
(40) has no periodic solutions; all solutions diverge to -∞.Proof.
Lety(t) be an arbitrary solution to (40). Taking into account that the function g(y)=rmaxy(1-y/Kmax) attains its maximum value rmaxKmax/4 at y=Kmax/2, we have, for all t∈[t0,t0+T],
(45)y(t0+T)-y(t0)=∫t0t0+T‍[r(t)y(t)(1-y(t)K(t))-σH(t)]dt≤∫IT∖Ineg‍[rmaxy(t)(1-y(t)Kmax)-σHmin]dt≤∫t0t0+T‍(rmaxKmax4-σHmin)dt=(rmaxKmax4-σHmin)T.
If (44) is satisfied, one has y(t0+T)<y(t0), and solution y(t) cannot be periodic. Similarly, for any n∈ℕ, (44) yields
(46)y(t0+(n+1)T)-y(t0+nT)≤(rmaxKmax4-σHmin)T<0.
Therefore, for any solution y(t) to (40), one has y(t)→-∞ as t→+∞.Remark 13.
Condition (44) perfectly agrees with the assumption (7) that forces all solutions to decay to -∞.Theorem2 assures the existence of the attractor-repeller pair for (40), for all
(47)σ<σ*=rminKmin4Hmax.
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(48)dy(t)dt=(4+cos(2πt))y(t)(1-y(t)8+cos(2πt))-σ(2+sin(2πt)).
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 ζ1(t) and ζ2(t),  9(1-205/15)/2<ζ1(t)<7(1-21/7)/2<7(1+21/7)/2<ζ2(t)<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(49)dy(t)dt=(7-3sin(2πt))y(t)(1-y(t)10+2sin(2πt))-3.8·(5+2cos(2πt)).
Observe that rminKmin/(4Hmax)=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σbif satisfying the inequality rminKmin/(4Hmax)≤σbif≤rmaxKmax/(4Hmin) such that (40) has (i) no periodic solutions if σ>σbif;  (ii) exactly one periodic solution if σ=σbif;  (iii) two periodic solutions if σ<σbif .Proof.
Define the Poincaré maph(y0,σ) for (40) on (0,T). Equation h(y0,σ)=y0 is quadratic, see, for instance, [11, page 211]. Since the second derivative hyy(y0,σ) is negative, the graph of h(y0,σ) is concave down. For σ<rminKmin/4Hmax, it crosses the line g(y0)=y0 at two points corresponding to two periodic solutions of (40). On the other hand, the graph has no intersections with this line for σ>rmaxKmax/4Hmin, and thus, in this case there are no periodic solutions to (40). The derivative hσ(y0,σ) is negative. Hence, for each fixed y0, the function δ(σ)=h(y0,σ) decreases as σ increases from rminKmin/4Hmax to rmaxKmax/4Hmin, and there exists a unique value σ=σbif  such that the graph of h(y0,σ) is tangent to the line g(y0)=y0.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(50)dy(t)dt=(5+sin(2πt))y(t)(1-2y(t)23+5cos(2πt))-σ2(11-3cos(2πt)).
For σ=1, (50) turns into (26) for which the attractor-repeller pair exists, see Example 4. One has rminKmin/(4Hmax)=9/7 and rmaxKmax/(4Hmin)=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 t advances, (i) solutions with small initial data x(t0)∈(0,x-(t0)) move away from x-(t), decay rapidly to -∞, and blow up in the future as t→tforw- (approach the repeller as t→-∞); (ii) solutions with initial data x(t0)∈(x-(t0),x+(t0)) leave the vicinity of the repeller x-(t) and approach the attractor x+(t) (correspondingly, leave the vicinity of the attractor x+(t) and approach the repeller x-(t) as t→-∞); these are so-called heteroclinic orbits; (iii) solutions with initial data x(t0)>x+(t0) approach the attractor x+(t), and blow up back in time as t→tback+. 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.



---

*Source: 167671-2013-03-27.xml*