The dynamics of some discrete models with delay under the effect of constant yield harvesting

Chaos, Solitons & Fractals 54 (2013) 26–38

Contents lists available at SciVerse ScienceDirect

Chaos, Solitons & Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos

The dynamics of some discrete models with delay under the effect of constant yield harvesting q Raghib Abu-Saris a,⇑, Ziyad AlSharawi b, Mohamed Ben Haj Rhouma c a

Department of Basic Sciences, College of Science and Health Professions, King Saud bin Abdulaziz University for Health Sciences, Riyadh, Saudi Arabia Department of Mathematics and Statistics, Sultan Qaboos University, P.O. Box 36, PC 123 Al-Khod, Oman c Department of Mathematics, Statistics and Physics, Qatar University, P.O. Box 2713, Doha, Qatar b

a r t i c l e

i n f o

Article history: Received 9 June 2011 Accepted 15 May 2013

a b s t r a c t In this paper, we study the dynamics of population models of the form xn+1 = xnf (xn1) under the effect of constant yield harvesting. Results concerning stability, boundedness, persistence and oscillations of solutions are given. Also, some regions of persistence and extinction are characterized. Pielous equation was considered as an example on these models, and a connection with a Lyness type equation has been established at certain harvesting level, which is used to give an explicit description of a persistent set. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Difference equations of the form xn+1 = xnf (xn), n = 0, 1,2,. . . are used to model single species with non-overlapping generations [19], where xn denotes the population density at discrete time n and f (xn) represents the per-capita growth rate of the population. The appropriate form of f (x) is usually chosen to reflect known observation about the studied species. For instance, in fish populations with high fertility rates and low survivorship to adulthood, Beverton and Holt [5] considered f ðxÞ ¼

lK , Kþðl1Þx

where

l > 1 is interpreted as the growth rate per generation and K is the carrying capacity of the environment. Straightforward re-scaling gives the simpler one parameter equation bxn xnþ1 ¼ 1þx . Interestingly, the classical Beverton–Holt (BH) n model can be derived using several arguments and approaches [39,8,13,10]. The BH model has been generalized to several more sophisticated models such as the Hassell model and the Maynard Smith–Slatkin model in order to make the model better accommodate certain observations

q This work is partially supported by SQU/ Internal Grants IG/SCI/ DOMS/09/05 and IG/SCI/DOMS/11/14 ⇑ Corresponding author. Tel.: +17 807572044. E-mail address: [email protected] (R. Abu-Saris).

0960-0779/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.chaos.2013.05.008

on real populations. Another prototype model of single species with non-overlapping generations is the Ricker model [29,27], which is obtained when f ðxÞ ¼ exp   rð1  kxÞ ; r; k > 0. In populations with substantial maturation time to sexual maturity, certain delay effect must be included in the function f (x), which motivates considering difference equations with delay [22] of the form

xnþ1 ¼ xn f ðxnm Þ;

n ¼ 0; 1; 2; . . . ; m P 1:

ð1:1Þ

For instance, Pielou [28, (page 80)] suggested taking lK f ðxnm Þ ¼ Kþðl1Þx in order to account for oscillations in nm certain populations. When the delay is one and after rescaling, Pielou’s equation takes the form

ynþ1 ¼

byn : 1 þ yn1

ð1:2Þ

Eq. (1.2) was studied in [17,20] where it was shown that when b > 1, every positive solution converges to the posi ¼ b  1. The main objective of this paper tive equilibrium y is to study the effect of harvesting on the dynamics of Eq. (1.1) with delay one, i.e., we consider the difference equation

xnþ1 ¼ xn f ðxn1 Þ  h;

ð1:3Þ

R. Abu-Saris et al. / Chaos, Solitons & Fractals 54 (2013) 26–38

where h is a fixed harvesting quota and f (x) obeys certain assumptions. In particular, f (x) is C2([0,1)) and satisfies the following conditions: (C1) f(0) = b > 1. (C2) f (x) is strictly decreasing on [0,1) and  1x 6 f 0 ðxÞ 6 0 for all x > 0. (C3) xf(x) is an increasing concave down function and xf(x) 6 M for some positive constant M. Observe that without condition (C1), there is no long term survival for any population regardless of the h value. Also, condition (C3) implies that limx?1 f (x) = 0. The effect of harvesting on the dynamics of populations governed by differential equations has been studied extensively by several authors [30–37]. In contrast, the effect of harvesting on populations governed by discrete models has received less attention. Much of the literature regarding discrete model involves, game-theoretic models in fishery, see [15] for a nice brief literature survey starting with the work of Levhari and Mirman [23] in the early 80’s. For example, in [6,7], it was shown that harvesting in the presence of a reserve area leads to higher levels of sustainability in exploiting fish stocks. In [14,15] dynamical models of commercial exploitations of renewable resources were studied assuming that two agents behave adaptively. In particular, they obtained results about the existence and stability of the positive equilibrium characterizing the sustainable use of the renewable resource. Recently, Al-Sharawi and Rhouma have studied the effect of both constant and periodic harvesting for the BH equation in constant and fluctuating capacity environment and have shown that constant rate harvesting in constant capacity environment can be the optimal strategy [1,2]. Under certain conditions, harvesting was also shown to lead to coexistence in an exclusive competitive environment with species governed by a discrete model [3]. This paper is structured as follows: We investigate local stability, boundedness and persistence of solutions in sections two and three. In section four, we investigate the global behavior of solutions. As a particular case of Eq. (1.3), we study the effect of constant harvesting on Pielou’s equation bxn with delay one in section five, i.e., xnþ1 ¼ 1þx  h. We charn1 acterize the persistent set, and use a 8-periodic solution of the Lyness type equation to construct an invariant region for Pielou’s equation when h < 1. 2. Equilibrium solutions and local stability Without harvesting (h = 0), Eq. (1.3) has two equilibrium solutions, namely, the origin  x1;0 ¼ 0 and another positive one, say  x2;0 ¼ c ¼ f 1 ð1Þ. A positive h shifts  x1;0 upward and  x2;0 downward till they collide and disappear when h > hmax, where

hmax :¼ maxfxðf ðxÞ  1Þ : x 2 ½0; cg:

ð2:1Þ

The equilibrium points are determined by the intersections between the curve y = x(f (x)  1) and the line y = h. Let 0 6 h 6 hmax. If we assume that we have more than two equilibrium points, then we can use the Mean Value

27

Theorem to obtain a contradiction with concavity assumption in (C3). To stress the role of h in the existence of the equilibrium points, we denote them by  x1;h and  x2;h , where  x1;h 6  x2;h . Thus, we have the following: Proposition 2.1. Consider hmax as defined in Eq. (2.1), and let 0 6 h < hmax. Eq. (1.3) has two equilibrium solutions. The two equilibrium solutions are equal when h = hmax. The linearized equation associated with Eq. (1.3) at a fixed point  x is given by

ynþ1  f ðxÞyn  xf 0 ðxÞyn1 ¼ 0: Þ, where  Let p :¼ f ð xÞ and q :¼  xf 0 ðx x 2 f x1;h ;  x2;h g. The linearized stability theorem gives the stability of each fixed point depending on the values of p and q. Notice that because  xf ð xÞ ¼  x þ h, then p ¼ 1 þ h= x and thus

1 6 p 6 b: Also, the fact that f is decreasing makes q negative which combined with assumption (C2) leads to the inequality

1 6 q 6 0: Condition (C3) implies that p + q > 0 for both fixed points  x1;h and  x2;h . However, since the graph of xf(x) crosses the line y = x with a slope greater than 1 at  x1;h , then p + q > 1 and thus  x1;h is a saddle. On the other hand the graph of xf(x) crosses the line y = x with a slope less than 1 at  x2;h , giving p + q < 1 and thus  x2;h is locally asymptotically stable when 0 < h < hmax. We summarize these facts in the following proposition. Proposition 2.2. Consider hmax as defined in Eq. (2.1), and let p ¼ f ð xÞ; q ¼  xf 0 ð xÞ. Each of the following holds true for Eq. (1.3). (i)  x1;h is a saddle and  x2;h is locally asymptotically stable for all 0 6 h < hmax. (ii) If h = hmax, then  x1;hmax ¼  x2;hmax is nonhyperbolic (one of the eigenvalues is equal to 1.) (iii)  x2;h starts as stable at h = 0, then as h varies, it can go to a repeller or nonhyperbolic before its disappearance. 3. Persistence and boundedness Persistence in population models is of paramount importance in general, and harvesting adds another factor of significance to its characterization in our equation. One can find several definitions pertaining to persistence (persistence, permanence, uniform persistence, strong persistence,. . . etc.) in the literature (see for instance [16,38]); however, we give ourselves the liberty to use the following definition throughout this paper. Definition 3.1. A solution of Eq. (1.3) is called persistent if the corresponding population survives indefinitely. We call a persistent solution strongly persistent if lim inf xn > 0. A 2 set D :¼ fðx; yÞ : ðx; yÞ 2 Rþ g is persistent if each solution of Eq. (1.3) with ðx1 ; x0 Þ 2 D is persistent. It is called the persistence set if it is the largest persistent set. We use Dh for the persistence set at a harvesting level h.

28

R. Abu-Saris et al. / Chaos, Solitons & Fractals 54 (2013) 26–38

If h < Mb and {xn} is a persistent solution of Eq. (1.3), then we have

xnþ2 ¼ xnþ1 f ðxn Þ  h ¼ xn f ðxn Þf ðxn1 Þ  hf ðxn Þ  h 6 Mb  hf ðxn Þ  h < Mb  h:

With the convention that empty sum is set to 0 and empty product is set to 1, we use induction on Eq. (1.3) to obtain

xn ¼

n2 Y j¼1

Also,

Thus, the following result becomes clear.

f ðxj Þ þ

j¼0

n2 Y

þ

ðxnþ1 þ hÞ 1 h xn ¼ > ðxnþ1 þ hÞ > : f ðxn1 Þ b b

f ðxj Þx0  h

" n2 Y

n2 Y

f ðxj Þ þ   

j¼1

#

f ðxj Þ þ f ðxn2 Þ þ 1 :

ð3:1Þ

j¼n3

Now, factor the coefficient of x0 and use the summation notation to obtain n2 Y

"

n1 Y k1 X

# 1

Lemma 3.1. Persistent solutions of Eq. (1.3) are bounded and strongly persistent.

xn ¼

The next results shows that the harvesting quota cannot exceed hmax.

for all n 2 N. Consider {xn} to be a persistent solution, then it is bounded, and the partial sums of the series in the numerator form an increasing bounded sequence, then

Theorem 3.1. If h > hmax, then Dh is empty.

Proof. Consider h > hmax and start with an initial condition (x1,x0) such that x0 6 x1, then

x1 ¼ x0 f ðx1 Þ  h 6 x0 f ðx0 Þ  h 6 M  h: If h > M then x1 is negative and Dh is empty. If M P h > hmax, then the function F (t) = tf (t)  h has no fixed points. In particular, x1 = x0f (x0)  h < x0 and x1 < x0 6 x1. and thus by induction, we obtain a decreasing sequence. This decreasing sequence either becomes negative and thus is non persistent or it converges to a positive fixed point, which is a contradiction. On the other hand, if x0 > x1, then either we obtain an increasing sequence or xm 6 xm1 for some positive integer m. If the sequence {xn} is increasing, then it is bounded because xf(x) 6 M by assumption and thus must converge to a positive fixed point which is a contradiction. If xm 6 xm1, we obtain that xm+1 6 xm 6 xm1 and again induction leads to a scenario similar to the first case, and hence, no population persists. h The next result shows that if a population starting above  x1;h goes below  x1;h at a certain time, then it is doomed to extinction. Proposition 3.1. Let 0 < h < hmax, and define the set 2 2 S 0 :¼ fðx; yÞ 2 Rþ : y < x1;h < xg [ fðx; yÞ 2 Rþ

:y<

h ; x 6 x1;h g: f ðxÞ  1

Then S 0 \ Dh ¼ /. Proof. Let ðx1 ; x0 Þ 2 S 0 , we show that ðx1 ; x0 Þ R Dh . If x0 < x1, then

f ðxj Þ x0  h

j¼1

1 Y k X

ð3:2Þ

ðf ðxj ÞÞ

k¼0 j¼1

1

ðf ðxj ÞÞ

k¼1j¼1

must converge. Thus, the limit of the nth term is zero, and consequently 1 Y k X x0 1 ¼ ðf ðxj ÞÞ : h k¼1j¼1

ð3:3Þ

We close this section by giving a result concerning the case when the two fixed points collide to form a nonhyperbolic fixed point. Theorem 3.2. Consider Eq. (1.3) with h = hmax. The nonhyperbolic fixed point  x1;h ¼  x2;h attracts all persistent solutions. Proof. Consider h = hmax and suppose that fxn g1 n¼1 is a persistent solution. Define An = f (xn1)(xn + h), then

Anþ1 ¼ gðxn Þ; An

where gðtÞ ¼

t f ðtÞ: tþh

Observe that g(t) 6 1 for all t 2 [0,1]. Indeed, if t > c, then

gðtÞ ¼

t f ðtÞ < f ðtÞ < f ðcÞ ¼ 1: tþh

However, if t 2 [0,c], then

gðtÞ ¼

1 1 1 ½tðf ðtÞ  1Þ þ t 6 ½hmax þ t ¼ 1: tf ðtÞ ¼ tþh tþh tþh

That said, we obtain

Anþ1 ¼ gðxn Þ < 1 whenever xn An

is not an equilibrium

x1 ¼ x0 f ðx1 Þ  h < x0 f ðx0 Þ  h < x0 :

Thus, An is decreasing. Since any persistent solution {xn} is bounded away from zero, the sequence An is also bounded and must converge to a positive value, say c0 > 0. Now,

By induction, we obtain a strictly decreasing sequence less than  x1 , which implies that ðx1 ; x0 Þ R Dh . If x1 < x0, then

Anþ1 ¼ An gðxn Þ

x0 <

h f ðx1 Þ  1

implies

implies x1 ¼ x0 f ðx1 Þ  h < x0 :

Thus, x1 < x0 and we apply the first scenario on (x0,x1). Hence, the proof is complete. h

lim gðxn Þ ¼ 1 and consequently xn is attracted to the nonhyperbolic equilibrium. h

29

R. Abu-Saris et al. / Chaos, Solitons & Fractals 54 (2013) 26–38

Theorem 4.2. If t 0 ¼ hb in Lemma 4.1, then any persistent solution fxn g1 1 of Eq. (1.3) satisfies xn P tn for all n 2 N.

4. A comparison principle In this section, we develop a useful comparison principle. Consider the difference equations

xnþ1 ¼ Fðxn ; xn1 Þ

ð4:1Þ

and

xnþ2 ¼ Fðxn ; xn1 Þ;

ð4:2Þ

where F (x,y) is increasing in x and decreasing in y. The next result gives a comparison principle for equations of the form (4.1) and (4.2) in a non-autonomous settings. Theorem 4.1. Suppose there exist two sequences an and bn such that an 6 xn 6 bn for all n P 1, and define fn+1 (x) = F (x,an) and gn+1(x) = F (x,bn), then each of the following holds true:

Proof. By definition, ‘‘if fxn g1 n¼1 is a persistent solution, then (xn1,xn) 2 Dh for all n P 0.’’ Indeed, if fxn g1 n¼1 is a persistent solution, then by Proposition 3.1 (xn1,xn) R S0 for all n P 0. So, for all n P 0, either   h   ; ðx xn P  P and x > Þ, x1;h and xn P f ðxn1 x x n n1 1;h 1;h Þ1 h Þ. As such, for all n P 0, x1;h and xn P f ðxn1 or ðxn1 6  Þ1

 xn P min x1;h ;

 h : f ðxn1 Þ  1

x1;h Þ ¼  x1;h . But h/(f (x)  1) is increasing with f(0) = b and f ð Therefore, xn P h/(b  1) for all n P 0. In other words, we established the fact that xn P h/(b  1) > h/b = t0 for all n P 0. Now, we assume xn P tk for all n P k, and use induction on k to show that xn P tk+1 for all n P k + 1. Indeed,

(i) The solution xn of Eq. (4.1) satisfies

ðg n  g n1     g 0 Þðx0 Þ 6 xnþ1 6 ðfn  fn1     f0 Þðx0 Þ (ii) The solution xn of Eq. (4.2) satisfies

ðg 2k g 2k2 g 0 Þðx0 Þ 6 x2k 6 ðf2k2 f2k4 f0 Þðx0 Þ and ðg 2kþ1  g 2k1     g 1 Þðx1 Þ 6 x2kþ1 6 ðf2k1  f2k3    f1 Þðx1 Þ

for all k P 1.

Proof. We use mathematical induction on n to prove part (i). For n = 0,

g 0 ðx0 Þ ¼ Fðx0 ; b1 Þ 6 x1 ¼ Fðx0 ; x1 Þ 6 Fðx0 ; a1 Þ ¼ f0 ðx0 Þ:

1 ðt k þ hÞ ¼ t kþ1 ; f ðt k Þ ðt k þ hÞ ðt k þ hÞ t k 6 xkþ3 ¼ xkþ2 f ðxkþ1 Þ  h 6 xkþ2 f ðt kþ1 Þ  h ) xkþ2 P ¼ t kþ1 P f ðt kþ1 Þ f ðt k Þ

t k 6 xkþ2 ¼ xkþ1 f ðxk Þ  h 6 xkþ1 f ðt k Þ  h ) xkþ1 P

and similarly xk+j P tk+1 for all j = 3,4, . . . . Thus, xn P tk+1 for all n P k + 1 and the proof is complete. h Proposition 3.1 shows that it is necessary for a persish tent solution {xn} to have x0 > b1 ; however, it is not necessarily sufficient. Consider t0 ¼ hb, and use Lemma 4.1 to obtain a sequence {tk}. Now, in conjunction with Eq. (3.3), define

sn :¼ h

n Y i1 X

1

ðf ðtj ÞÞ :

i¼1 j¼0

Assume the statement is true for n = k  1, i.e.,

The sequence {sn} is a strictly increasing sequence bounded above by f ðx h Þ1 ¼  x1;h . Thus, we define

ðg k1  g k2      g 0 Þðx0 Þ 6 xk 6 ðfk1  fk2      f0 Þðx0 Þ;

c :¼ n!1 lim sn

then for n = k, we obtain

xkþ1 ¼ Fðxk ; xk1 Þ 6 Fðxk ; ak1 Þ ¼ fk ðxk Þ 6 ðfk  fk1      f0 Þðx0 Þ

1;h

and give the following result. 1

Proposition 4.1. Define S 1 :¼ fðx; yÞ : y < ðc þ hÞðf ðxÞÞ S 1 \ Dh ¼ /.

g.

and

xkþ1 ¼ Fðxk ; xk1 Þ P Fðxk ; bk1 Þ ¼ g k ðxk Þ P ðg k  g k1      g 0 Þðx0 Þ: The proof of (ii) and is therefore omitted. h Lemma 4.1. Consider tnþ1 ¼ ðtfnðtþhÞ such that t0 > 0. The nÞ x1;h in an increasing fashion if sequence {tn} converges to  x1;h and on a decreasing fashion if  x1;h < t 0 <  x2;h . 0 < t0 <  ðtþhÞ , f ðtÞ

Proof. Define GðtÞ :¼ then G(t) is increasing with the same fixed points as Eq. (1.3). Since Gð0Þ ¼ hb, then G(t) > t x1;h ; GðtÞ < t for  x1;h < t <  x2;h and G(t) > t for for 0 < t <  t >  x2;h . Now, the fact that tn = Gn(t0) and a simple stair step diagram completes the proof.

Proof. Suppose ðx1 ; x0 Þ 2 S 1 \ Dh , and let fxn g1 n¼1 be the associated persistent solution. We use the comparison principle of Theorem 4.1 to obtain a contradiction. Consider the sequence an to be the sequence tn obtained in Lemma 4.1 with t0 ¼ hb, and define fn+1(x) = f(an)x  h. By part (i) of Theorem 4.1, xn+1 < fnfn1    f0(x0). However, bx0 when ðx1 ; x0 Þ 2 S 1 ; x1 ¼ 1þx  h < c and fnfn1    f0(x0) 1

goes negative for sufficiently large n. This contradicts our assumption that {xn} is persistent. h Theorem 4.2 and Proposition 3.1 show that persistent solutions with initial conditions below  x1;h (if any) either converge to  x1;h or go up and stay above  x1;h . Now, we proceed to establish certain bounds on the asymptotic behavior of persistent solutions.

30

R. Abu-Saris et al. / Chaos, Solitons & Fractals 54 (2013) 26–38

Lemma 4.2. Let 0 6 h < hmax. The function gðtÞ ¼ h f ðtÞff ðtÞþ1 x1;h Þ1 ð is increasing and has exactly two fixed points in the interval ½0; f 1 ð1=f ð x1;h ÞÞ. One of the fixed points is  x1;h and the other one (say x2;h ) is larger than  x2;h .

Proof. It is clear by observing that g 0 ðtÞ > 0; gð x1;h Þ ¼  x1;h ; gð x2;h Þ <  x2;h and g(t) ? 1 as t approaches f 1 ð1=f ð x1;h ÞÞ from the left. Theorem 4.3. Let 0 6 h < hmax. and x2;h as in Lemma 4.2. Any persistent solution fxn g1 n¼1 of Eq. (1.3) satisfies

x1;h 6 lim inf xn 6 lim sup xn 6 x2;h < f 1





1 : f ðx1;h Þ

to observe that persistent solutions that are non-monotonic must oscillate about the curves of y  x = 0 and y  xf(x) + h = 0. Otherwise, if xn+1 < xnf (xn)  h for all n, then xnf (xn1)  h < xnf (xn)  h for all n and, consequently, xn1 > xn for all n, i.e., {xn} is monotonic which is a contradiction. Other cases can be handled similarly. Proposition 5.1. A persistent solution of Eq. (1.3) either converges to an equilibrium point or oscillates about the curve y ¼ ðxþhÞ . f ðxÞ Proof. Let {xn} be a persistent solution. If xn eventually stays on one side of the curve, say xnþ1 < ðxf nðxþhÞ for all nÞ n P n0, then Eq. (1.3) gives

xnþ2 ¼ xnþ1 f ðxn Þ  h < xn

lim inf xn P x1;h :

for all n P n0. This leads to convergence to a periodic solution of period 2. If this is the case then there exists positive values a – b such that

Next, write Eq. (1.3) in the form of Eq. (4.2) by taking

b ¼ af ðbÞ  h and a ¼ bf ðaÞ  h:

Proof. From Theorem 4.2 and Lemma 4.1, we obtain

Fðxn ; xn1 Þ ¼ xn f ðxn Þf ðxn1 Þ  hf ðxn Þ  h; h b

then again, consider an to be the sequence tn with t 0 ¼ as given in Lemma 4.1. Now, we take fn(x) = f(an)f (x)x  h f(x)  h and use part (ii) of Theorem 4.1 to obtain

x2k 6 f2k2  f2k4      f0 ðx0 Þ and

x2kþ1 6 f2k1  f2k3      f1 ðx1 Þ for all k P 1. In both cases, the right hand side is eventually smaller than the larger fixed point of f1(x), which is given by x2;h as clarified in Lemma 4.2. Thus, we obtain

lim sup xn 6 x2;h < f 1





1 f ðx1;h Þ

and the proof is complete. h Theorem 4.3 shows that a persistent solution is either converging to  x1;h monotonically or eventually larger than  x1;h . Since  x1;h is a saddle all the time, then orbits that converge to  x1;h are the ones assured by the stable manifold theorem. Those orbits show that harvesting can alter the oscillatory behavior of populations; however, those populations make no significant contribution to the persistent set, and therefore, we are more concerned about populations that are larger than  x1;h . Thus, without further mention, we limit our attention to solutions that satisfy xn P  x1;h for all n. 5. Oscillation of solutions In this section, we discuss the oscillatory character of solutions of Eq. (1.3). If a persistent solution {xn} of Eq. (1.3) is neither eventually less than nor larger than  x2;h , then we call it oscillatory about  x2;h , or oscillatory for short. Notice that because of Theorem 4.2, we cannot have solutions oscillating about  x1;h . We call a solution {xn} of Eq. (1.3) oscillatory about a curve H(x,y) = 0 if Xn = (xn,xn+1) does not eventually stay on one side of the curve. It is easy

Without loss of generality, if we assume that a < b, then we have

f ðaÞ ¼

aþh b

<

bþh

a

¼ f ðbÞ

which contradicts the assumption that f is decreasing. Since Eq. (1.3) has no periodic solutions of minimal period 2, then the convergence must be to an equilibrium solution. The other side is similar and omitted. h We associate orbits of Eq. (1.3) in the positive quadrant by the map

Tðx; yÞ ¼ ðy; yf ðxÞ  hÞ:

ð5:1Þ

It is straightforward to verify that T is one to one on the positive quadrant (the x axis is not included), and therefore, the action of Eq. (1.3) on a region can be visualized by testing the action of T on the boundary of that region. For a fixed constant c > 0, points above y = c are mapped to points on the right of x = c. Also, since T(c,t) = (t,tf(c)  h) and x1 = x0f (x1)  h < 0 if and only if x0 < h/f (x1), then one can square the persistent set inside the box Sc ¼ fðx; yÞ : 0 6 x < c; 0 < y < cg provided that f (x) < h/ (xf(c)  h). It is worth stressing again that Theorem 4.2 and Lemma 3.1 limit our interest to the region fðx; yÞ : x; y P  x1;h g as persistent solutions with initial conditions below  x1;h (if any) either converge to  x1;h or go up and stay above  x1;h . We define the following subregions:

R1 :¼ fðx; yÞ 2 R2 : x1;h < x < x2;h ; y P x2;h g; R2 :¼ fðx; yÞ 2 R2 : x P x2;h ; y > x2;h g; R2i :¼ fðx; yÞ 2 R2 : y > xg; R2ii :¼ fðx; yÞ 2 R2 : y 6 xg; R3 :¼ fðx; yÞ 2 R2 : x1;h < y 6 x2;h ; x > x2;h g; R4 :¼ fðx; yÞ 2 R2 : x1;h < x 6 x2;h ; x1;h < y < x2;h g; R4i :¼ fðx; yÞ 2 R4 : y < xg; R4ii :¼ fðx; yÞ 2 R4 : y P xg:

31

R. Abu-Saris et al. / Chaos, Solitons & Fractals 54 (2013) 26–38

For our convenience, if (x,y) is in one of the above regions but yf ðxÞ  h 6  x1;h , then the population with initial conditions (x1,x0) = (x,y) goes extinct (see Proposition 3.1), and therefore, we say T(x,y) 2 {0}. Now, by tracing the action of T on the defined boundaries of these regions, we conclude the following result: Proposition 5.2. Let 0 < h < hmax. Each of the following holdstrue: (i) (ii) (iii) (iv) (v) (vi)

TðR1 Þ  R2i . TðR2i Þ  R2ii if h = 0 and TðR2i Þ  R2i [ R2ii . TðR3 Þ  R4i [ f0g. T 2 ðR2i Þ  R2ii [ R3 [ f0g. TðR4ii Þ  R4ii [ R1 . TðR4i Þ  R4i [ R4ii [ f0g.

We recall from [20] that given a fixed point  x for a sequence {xn}, a positive semi-cycle is subsequence of consecutive terms that are all greater or equal to  x. A negative semi-cycle is subsequence of consecutive terms that are all smaller or equal than  x. The next result becomes straightforward. Proposition 5.3. Consider Eq. (1.3) with 0 6 h < hmax. Each of the following holds true: (i) Non-equilibrium persistent solutions that do not converge to  x1;h monotonically must oscillate about  x2;h . (ii) When h is zero, a positive semi-cycle has more than two terms (except possibly the first one) with a maximum in the second term, and decreasing afterwards. (iii) When h is zero, a negative semi-cycle has at least two terms (except possibly the first one) with a minimum in the second term, and increasing afterwards. For h > 0 and by checking the effect of T, T2 on the regions R1 and R3 , one can observe that an increase in the harvesting quota increases the length of semi-cycles. Therefore, solutions overshoot or undershoot the stable equilibrium when they oscillate, and consequently, the risk of extinction for populations governed by Eq. (1.3) becomes high. However, obtaining exact values of h that do not lead to population extinction remains an extremely challenging task. Nevertheless, one can use the extreme values of semi-cycles and the approach used in [17,20,9] to prove that certain persistent solutions are attracted to  x2;h under some mild conditions on h. Define g1(t) = t(f(t)  1) and g 2 ðtÞ ¼ ðh=f ð x2;h ÞÞf ðtÞ, then g1 = g2 at t¼ x2;h . Furthermore, for small values of h, they intersect at a second positive point smaller than  x2;h . Thus, g1(t) 6 g2(t) for all t P  x2;h and for sufficiently small values of h. As we increase h, we continue to obtain g1(t) 6 g2(t) for all t P  x2;h as long as  x2;h is the largest solution of g1 = g2. This fact gives us a constraint on h, which we use to obtain the following result: Lemma 5.1. Assume that tf(t) is increasing. Let 0 6 h < hmax x2;h is the largest solution of g1 = g2. A positive semisuch that  cycle of a persistent and oscillatory solution of Eq. (1.3) has at least three terms (except possibly the first one) with a maximum in the second or third term, and decreasing afterwards.

Proof. We trace the boundary of R1 . Since

T 2 ðx2;h ; tÞ ¼ Tðt; tf ðx2;h Þ  hÞ ¼ ðtf ðx2;h Þ  h; ðtf ðx2;h Þ  hÞf ðtÞ  hÞ: x2;h ; tf ð x2;h Þ  h > t and T(t,t) = (t,tf(t)  h) implies For t >  tf ðtÞ  h >  x2;h . Thus, a positive semi-cycle has at least three terms. On the other hand,

tf ðx2;h Þ  h > ðtf ðx2;h Þ  hÞf ðtÞ  h () g 1 ðtÞ < g 2 ðtÞ and therefore, the condition on h makes the maximum takes place at the second or third term. Next, define the functions

f0 ðtÞ ¼ x2;h f ðtÞ  h; G0 ðtÞ ¼ f1 ðf0 ðtÞÞ;

f 1 ðtÞ ¼ f ðx2;h Þt  h; and G1 ðtÞ ¼ f1 ðG0 ðtÞÞ:

The motivation behind taking these functions will be clear in the proof of Lemma 5.2 and Theorem 5.1; however, it is easy to observe that G0 and G1 are both decreasing, and Gj ð x2;h Þ ¼  x2;h . For j = 0, . . . ,3, let C rþj :¼ fxkj þ1 ; xkj þ2 ; . . . ; xkjþ1 g be four consecutive semi-cycles starting with the negative semi-cycle Cr. Also, let cr+j, j = 0, . . . ,3 be the extreme values of the semi-cycles respectively. Next, assume that xf(x) is increasing. If cr+3 and cr+2 lie in the second term, then

crþ3 ¼ xk3 þ2 ¼ xk3 f ðxk3 Þf ðxk3 1 Þ  hf ðxk3 Þ  h 6 x2;h f ðx2;h Þf ðxk3 1 Þ  hf ðx2;h Þ  h 6 G0 ðcrþ2 Þ and

crþ2 ¼ xk2 þ2 ¼ xk2 f ðxk2 Þf ðxk2 1 Þ  hf ðxk2 Þ  h P x2;h f ðx2;h Þf ðxk2 1 Þ  hf ðx2;h Þ  h P G0 ðcrþ1 Þ: Which implies

crþ3 6 G0 ðcrþ2 Þ 6 G0 ðG0 ðcrþ1 ÞÞ ¼ G20 ðcrþ1 Þ: If cr+3 and cr+2 lie in the third term, then

crþ3 ¼ xk3 þ3 ¼ xk3 þ2 f ðxk3 þ1 Þ  h 6 f1 ðxk3 þ2 Þ 6 f1 ðG0 ðcrþ2 ÞÞ and

crþ2 ¼ xk2 þ3 ¼ xk2 þ2 f ðxk2 þ1 Þ  h P f1 ðxk2 þ2 Þ P f1 ðG0 ðcrþ1 ÞÞ: Which implies

crþ3 6 G1 ðcrþ2 Þ 6 G1 ðG1 ðcrþ1 ÞÞ ¼ G21 ðcrþ1 Þ: Similarly, if cr+2 and cr+3 lie in the third and second terms respectively, then we obtain

crþ3 6 G0 ðcrþ2 Þ 6 G0 ðG1 ðcrþ1 ÞÞ and if cr+2 and cr+3 lie in the second and third terms respectively, then we obtain

crþ3 6 G1 ðcrþ2 Þ 6 G1 ðG0 ðcrþ1 ÞÞ: Now, we can summarize in the following result: Lemma 5.2. Suppose that x f (x) is increasing. Let cr+j, j = 1,2,3 be the extreme values of three consecutive semi-cycles starting with a positive one. Each of the following holds true: (i) If cr+2 and cr+3 lie in the second term, then crþ3 6 G20 ðcrþ1 Þ.

32

R. Abu-Saris et al. / Chaos, Solitons & Fractals 54 (2013) 26–38

(ii) If cr+2 and cr+3 lie in the third term, then crþ3 6 G21 ðcrþ1 Þ. (iii) If cr+2 and cr+3 lie in the second and third terms respectively, then cr+3 6 G1(G0((cr+1))). (iv) If cr+2 and cr+3 lie in the third and second terms respectively, then cr+3 6 G0(G1((cr+1))).

pffiffiffi 2 satisfied if b 6 4. Now, we have hmax :¼ ð b  1Þ , and when 0 6 h < hmax, the two positive equilibria are given by

xi;h ¼

b  1  h þ ð1Þi

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ðb  1  hÞ  4h 2

;

i ¼ 1; 2: The next result is a classical one, and the proof is straightforward. Lemma 5.3. Let g 2 Cð½0; cÞ be decreasing such that g(c) = 0 and 0 < g(0) < c. The map g has a unique fixed point in (0,c). Moreover, the fixed point is a global attractor for the interval [0,c] if and only if g has no periodic points of period two. Now, we give the main result of this section. Theorem 5.1. Suppose that tf(t) is increasing, and let 0 6 h < hmax such that  x2;h is the largest solution of g1 = g2. If G21 ð0Þ > 0 and G1 has no periodic points of period two, then for any persistent solution {xn} of Eq. (1.3), we have lim sup xn 6  x2;h .

Proof. From Lemma 5.1, the extreme values of positive semi-cycles take place at the second or third term. When t> x2;h , we have

> G1 G0 ðtÞ and G0 G1 ðtÞ > G20 ðtÞ: x2;h , then G21 ðtÞ ¼ f1 G0 G1 ðtÞ > G0 Also, since f1(t) > t for t >  G1 ðtÞ. Thus, Lemma 5.2 narrows down to crþ3 < G21 ðcrþ1 Þ. Now, for a persistent and oscillatory solution fxn g1 n¼1 take fcm g1 m¼0 to be the sequence of extreme values of the positive semi-cycles (ordered as they appear in {xn}). Then we have cm < G2m 1 ðc 0 Þ and the conditions on G1 together with Lemma 5.3 give us the conclusion. h Corollary 5.1. Suppose the hypotheses of Theorem 5.1 are satisfied, then persistent and oscillatory solutions are attracted to  x2;h . We close this section by the following remark: Remark 5.1. In Theorem 5.1, the conditions on G1 could be stronger than the condition on  x2;h to be the largest solution of g1 = g2, and therefore, the latter one could be redundant. However, keeping the two conditions saved us some deep technicalities in the proof.

6. Pielou’s equation with harvesting In this section, we focus on a particular case of Eq. (1.3) and consider Pielou’s Equation with harvesting

xnþ1

It is worth mentioning that Eq. (6.1) was studied when h = 0 by Kuruklis and Ladas in [21]. In Section 2, we discussed the local stability of these equilibria in general; however, we can be more specific here. If h < 1, then  x2;h is stable. At h ¼ 1;  x2;h is a nonhyperbolic fixed point and Eq. (6.1) reduces to a Lyness type equation. If h > 1, then   x2;h is a repeller. pffiffiffi x1;h is a saddle for all h < hmax. At h = hmax,  x1;h ¼  x2;h ¼ b  1 is a nonhyperbolic point. It is stable when b < 4 and unstable when b > 4. In the remaining of this section, we first study the special case h = 1 which, as mentioned earlier, reduces Eq. (6.1) to Lyness equation allowing us to use previous results to establish both the persistence set and the general behavior of the solutions. Then, in subSection (6.2) we find a range of parameters (b,h) where solutions converge to the larger equilibrium  x2;h . The last subsection is dedicated to narrowing down the persistence set Dh. While this is a challenging task requiring tedious calculations, we feel that it is necessary to attempt to describe the set of initial populations that will survive a certain level of harvesting h. 6.1. A connection with a Lyness type equation

G1 ðtÞ ¼ f1 ðG0 ðtÞÞ < G0 ðtÞ implies G21 ðtÞ

bxn ¼  h; 1 þ xn1

ð6:2Þ

When h = 1, Eq. (6.1) reduces to the well-known Lyness equation (or its conjugate form) [24–26]. Indeed, we can use the change of variables b yn = 1 + xn, to obtain

ynþ1 ¼

yn  1b ; yn1

n 2 N:

ð6:3Þ

Eq. (6.3) has been extensively studied [4,12,11,18]. Notice that when h = 1 and 1 < b < 4, then Theorem 3.1 gives an empty persistent set. The case h = 1 and b = 4 has been discussed in [11] where certain solutions can be written in explicit form. The case h = 1 is of particular interest to us, since we use one of the 8-periodic solutions of Eq. (6.3) to introduce a trapping region for Eq. (6.1). Theorem 6.1. Consider Eq. (6.1) with h = 1 and b P 4. Each of the following holds true. (i) The curves given by

 I b ðx; yÞ :¼ 1 þ

b 1þx

 1þ

b ð1 þ x þ yÞ ¼ C; 1þy

define the invariants of Eq. (6.1), where C is a constant given by I b ðx1 ; x0 Þ. (ii) A population persists if and only if the constant C in (i) is bounded by 2

C 1 :¼ 2 þ ðb þ 1Þ  x2 ðb  4Þ 6 C 6 b > 1:

ð6:1Þ

b In this equation, we have f ðxÞ ¼ 1þx . Furthermore, it is straightforward to check that conditions (C1)-(C3) are

2

C 2 :¼ 2 þ ðb þ 1Þ  x1 ðb  4Þ: (iii) The invariant I b ðx; yÞ ¼ C 2 defines the boundary of the persistence set D1 .

33

R. Abu-Saris et al. / Chaos, Solitons & Fractals 54 (2013) 26–38

(iv)  x2 is a center and the positive solutions inside D1 oscillate about  x2 such that every semi-cycle (except possibly the first one) has at least three elements.

f1 ðbx2;h  hx1;h  2h þ hbÞ < f01

hð2 þ x1;h Þ ð1 þ x1;h Þ2

!

and consequently

ð1 þ x1;h Þðbx2;h  hx1;h  2h þ hbÞ  h Proof. (i) and (ii) can be extracted from [18,12]. Since b x ¼ ð1 þ  xÞ2 , by direct computations, which can be checked easily using a symbolic computational language such as MAPLE,

<

bð1 þ x1;h Þ2  1: x1;h ð2 þ x1;h þ ð1 þ x1;h Þ2 Þ

ð6:5Þ

Again, we rewrite this inequality as

2

C 1 ¼ Ib ðx2 ; x2 Þ ¼ 2 þ ðb þ 1Þ  x2 ðb  4Þ

2 2 ðb  h þ 2hÞx2;h þ hðh þ 2h þ bh  1Þ þ 1

and

<

2 C 2 ¼ Ib ðx1 ; x1 Þ ¼ 2 þ ðb þ 1Þ  x1 ðb  4Þ:

To prove (iii), notice that when b = 4 we have  x1 ¼  x2 and thus C1 = C2. Accordingly, using (ii) we obtain the persistent set

b 2

ððb þ 2Þh þ ðb þ b þ 1Þx1;h þ hx2;h Þ

or 2 4

x2;h >

I 2ðb1 2Þ ðx; yÞ ¼ 18b1  24;

I b1 ðx; yÞ ¼ 9ðb1  1Þ;

I 5 ðx; yÞ ¼ 36; pffiffiffi where b1 :¼ 4 þ 7 are 7-periodic, 7-periodic, 8-periodic respectively. It is worth mentioning that it is possible for a fixed b > 4 to have two invariants with different periods. For instance, if we fix b = 5, then I 5 ðx; yÞ ¼ 36 and I 5 ðx; yÞ ¼ 75 are 8-periodic and 10-periodic respectively. 2

2

3

3

3

2

2 3

2

4

2

2

2

b h þ ð3b þ b þ 1Þh  ð2 þ 2b þ 3b Þh þ b þ b þ 1

:

ð6:6Þ

Graph Fig. 1 shows the solution of this inequality in the (b,h)-plane. 6.3. Characterizing the persistent set Finding the exact persistent set when h – 1 is a noneasy task; however, in this section, we characterize the persistent set through bounds, trapping regions and numerical psimulations. We start by the case ffiffiffi 2 h ¼ hmax ¼ ð b  1Þ . (See Figs. 2–4). pffiffiffi 2 Theorem 6.2. Consider h ¼ ð b  1Þ . Each of the following holds true:

(i) If b > 4, then

Dh  fðx; yÞ : 0 6 x 6 pffiffiffi 6 y 6 b  1g:

6.2. Convergence to the large equilibrium In this subsection, we apply the results of Section 5 on Pielou’s equation. For Eq. (6.1), the condition g1(t) 6 g2(t) in Lemma 5.1 simplifies to

2

b h  b ð3 þ 2bÞh þ ð2b þ 3b  1 þ b Þh þ ðb þ 2Þh  1

D1 ¼ fðx; yÞ : C 1 6 Ib ðx; yÞ 6 C 2 g and the boundary of this region is given by the invariant I b ðx; yÞ ¼ C 2 . Finally, (iv) follows from (i) and Proposition 5.2. h The map in (5.1) rotates the invariants I b ðx; yÞ ¼ C; C 1 6 C 6 C 2 clockwise such that each orbit along the invariant is periodic of the same period or dense in the invariant. For instance, the orbits along the invariants

3

pffiffiffi pffiffiffi pffiffiffi b  1; ð b  1Þðx þ 2  bÞ

(ii) If b = 4, then Dh  fðx; yÞ : 3 þ 16x  x2 þ

t 2 þ ð1  bÞt þ hð1 þ x2;h Þ P 0;

y¼

and that is satisfied as long as

0 6 x 6 1g:

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 49  148x þ 150x2  52x3 þ x4 ; 2ðx þ 5Þ

2

h6

ðb  1Þ : 2ðb þ 1Þ

ð6:4Þ

This constraint on the harvesting level guarantees that the extreme values of positive semi-cycles in an oscillatory and persistent solution take place at the second or third term. Next, we test the conditions of Theorem 5.1. The condition G21 ð0Þ > 0 implies G1 ð0Þ < f01 f11 f11 ð0Þ, or equivalently

  h : f1 ðf1 ðbx2;h  hÞÞ < f01 f11 ð1 þ x1;h Þ

b Keep in mind that f ðxÞ ¼ 1þx , which implies f ð x2;h Þ ¼ 1 þ  x1;h and  x1;h  x2;h ¼ h. Further computations give us,

hþxn Proof. From the proof of Theorem 3.2, An ¼ 1þx conn1 pffiffiffi  verges to A ¼ b  1 on a decreasing fashion. Thus, pffiffiffi which implies h þ xnþ1 > ð b  1Þð1 þ xn Þ, pffiffiffi pffiffiffi xnþ1 P ð b  1Þðxn þ 2  bÞ, and we obtain (i). (ii) can be obtained from the previous section. In particular, when h = 1 and b = 4, we obtain one invariant curve given by I 4 ðx; yÞ ¼ 27. If ðxn1 ; xn Þ 2 I 4 where xn < xn1, then

xnþ1 ¼

4xn 4xn 1<  1 < xn : 1 þ xn1 1 þ xn

Thus the population decreases to extinction along the invariant curve. So, the lower branch of the invariant curve

34

R. Abu-Saris et al. / Chaos, Solitons & Fractals 54 (2013) 26–38

2

1.5

1

0.5

0 0

1

2

3

4

5

6

pffiffiffi 2 2 Fig. 1. This figure shows the curves of hmax ¼ ð b  1Þ ; h ¼ ðb1Þ ; h ¼ 1. Also, the red dashed curve gives the upper boundary of the region given by 2ðbþ1Þ Inequality (6.6). (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)

1.4

1

1.2 0.8 1

y

0.8

y

0.6

0.6

0.4

0.4 0.2 0.2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

x

0

0.2

0.4

x

0.6

0.8

1

Fig. 2. These figures show numerical simulations of the persistent set at the nonhyperbolic fixed point. The fixed point is indicated by a bold black dot.

does not belong to the persistent set. If ðxn1 ; xn Þ 2 I 4 where xn > xn1, then

xnþ1 ¼

4xn 4xn 1>  1 > xn : 1 þ xn1 1 þ xn

This inequality is obvious due to the fact that xn > 1þx4n1 . Hence, the persistent set is the upper branch of I 4 , which is what we have in (ii). h

Theorem 6.3. Consider S 0 as defined in Proposition 3.1. If pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 h > 13 ðb  2 þ 2 b þ b þ 1Þ, then

Dh  fðx; yÞ : 0 < x < b  1; 0 < y < b  1g n S 0 : Proof. Since few iterates of the map T in (5.1) are sufficient to take initial conditions (x1,x0) with x0 P b  1 to ðxn0 1 ; xn0 Þ with xn0 < b  1 and xn0 1 P b  1, then we start with ðxn0 1 ; xn0 Þ and show that we have no persistence. We have

35

R. Abu-Saris et al. / Chaos, Solitons & Fractals 54 (2013) 26–38

Now, the given condition on h implies

4

bðxn0  hÞ h  h < ð1 þ xn0 Þ; b ð1 þ xn0 Þ and Proposition 4.1 shows that there is no persistence. h

3

y

In the rest of this section, we use an 8-periodic solution of the case h = 1, to define a trapping region when h < 1. To achieve this task and avoid complicated computations, we transform Eq. (6.1) to a new form where the origin becomes the small equilibrium. Take zn ¼ xn   x1;h to obtain

2

znþ1 ¼

1

b 1þx1;h

hþx

zn  1þx1;h zn1 1;h

zn1 1 þ 1þ x

ð6:7Þ

:

1;h

Now, define 0

yn :¼ 1

2

3

4

x Fig. 3. The shaded region in this figure shows the set of initial conditions that survived 500 iterations at h = 0.90 and b = 4.25, where x = xn1 and y = xn. The blue circles show an orbit oscillating and converging to  x2 . The fading loop shows the largest invariant curve at h = 1,b = 4.25. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)

bxn0  h 6 xn 0  h 1 þ xn0 1  bxn0 þ1 bxn0 bðxn0  hÞ b ¼ h¼ h h6  h: 1 þ xn0 1 þ xn0 1 ð1 þ xn0 Þ 1 þ xn 0

xn0 þ1 ¼ xn0 þ2

zn ; 1 þ x1;h

a :¼

b ; 1 þ x1;h

and b :¼

h þ x1;h ; 1 þ x1;h

to obtain

ynþ1 ¼

ayn  byn1 1 þ yn1

;

n 2 N:

ð6:8Þ

The positive equilibrium of Eq. (6.8) is given by  ¼ a  b  1 and that should not be confused with  y x2;h . Recall the map T introduced in (5.1), which takes the ybx form Tðx; yÞ ¼ ðy; a1þx Þ. Next, define qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A :¼ a ða4  2a3 þ 4ab  4b2 þ a2 þ 4b2 a  2ba2  2a3 b þ b2 a2 Þ

s :¼

1 ða3  a2 b  2ab  a2 þ 4b þ AÞ 4bða  1Þ

Fig. 4. This figure shows the trapping region inside the red bold curve when h < 1. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)

36

R. Abu-Saris et al. / Chaos, Solitons & Fractals 54 (2013) 26–38

and

s⁄ > a  b  1 if and only if

ða  bÞt bt ; f 2 ðtÞ :¼ 1þt a1t tða2  ab  b  atÞ tða2  a  b  atÞ g 1 ðtÞ :¼ ; g 2 ðtÞ :¼ : bða  1Þ ab

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð3h  x Þx2 2;h 2;h 2 ðx42;h  6hx22;h þ h Þ >  h; ðh þ x2;h Þ

f1 ðtÞ :¼

To assure that A is real, we need to impose the following constraint on a and b:

1 1 a P a0 :¼ ðb þ 1Þ þ 2 2

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi b2 þ ð10 þ 8 2Þb þ 1:

ð6:9Þ

To keep things in perspective of the previous section, we  x have a ¼ 1 þ  and Inequality (6.9) is equivalent x2;h ; ab ¼ 2;h h to

x2;h

pffiffiffi 2ð2 þ 2Þh P ðb  h  1Þ

()

pffiffiffiffiffiffi b P h þ 1 þ 2 2h:

pffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h < b þ 3  2 2ðb þ 1Þ; 1 < b 6 2ð1 þ 2Þ pffiffiffi bðb2Þ ; b > 2ð1 þ 2Þ: h < 2ðbþ2Þ

ð6:14Þ

The solution of Inequality (6.14) is given by Condition (6.11). (iii) We show that T4(s⁄,s⁄) lies along the diagonal y = x. Computations show that T

T

fðx; yÞ : x ¼ g 2 ðyÞg ! fðx; yÞ : y ¼ f2 ðxÞg !fðx; yÞ : y ¼ xg and T

ð6:10Þ

Here, computations are becoming tedious. So, we keep our writing concise and use the MAPLE computer algebra system to perform computations. Condition (6.9) assures that A is real, and to assure that s⁄ is larger than the positive  ¼ a  b  1, we impose the following equilibrium y constraint:

(

which can be written as

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4h 2 ð1 þ x2;h Þ  ðb  h  1Þ: ðb  h  1Þ  8h > b

T

fðx; yÞ : y ¼ xg ! fðx; yÞ : y ¼ f1 ðxÞg ! fðx; yÞ : y ¼ g 1 ðxÞg: In particular,



 T   T T g 2 f21 ðs Þ; f21 ðs Þ ! f21 ðs Þ; s ! ðs ; s Þ ! ðs ; f1 ðs ÞÞ

and T

T

ðs ; f1 ðs ÞÞ ! ðf1 ðs Þ; g 1 f1 ðs ÞÞ ! ðg 1 f1 ðs Þ; f2 g 1 f1 ðs ÞÞ T

 ! ðf2 g 1 f1 ðs Þ; f2 g 1 f1 ðs ÞÞ: ð6:11Þ

Proposition 6.1. Each of the following holds true: (i) If Condition (6.10) is satisfied, then A is real.  ¼ a  b  1. (ii) Condition (6.11) assures that s > y (iii) The choice of s⁄ assures that T4(s⁄,s⁄) lies along the diagonal y = x. (iv) If h < 1 and Condition (6.11) is satisfied, then Tj(s⁄,s⁄),j = 1,2,3,4 stay in the positive quadrant.

We alert the reader that computations here were tedious; however, if you take y1 = y0 = s⁄ as given in Eq. (6.12) and take x2 for  x2;h , then use the MAPLE commands > y[-1]:¼s⁄; y[0]:¼s⁄; alpha:¼x[2]+1; beta:¼h⁄(x[2]+1)/(h + x[2]); > For j from 0 to 3 do > y[j + 1]:¼factor (simplify (expand (rationalize ((alpha⁄y[j]-beta⁄y[j-1])/ (1 + y[j-1]))))): > end do;

to obtain x2;h and Proof. (i) Write a ¼ 1 þ 

b¼

y1 ¼

ðh þ x1;h Þx2;h hðx2;h þ 1Þ ¼ ; ðx2;h þ hÞ ð1 þ x1;h Þx2;h

y2 ¼

then we obtain

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1 þ x2;h Þ2 2 ðx42;h  6hx22;h þ h Þ A¼ ðh þ x2;h Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2;h ð1 þ x2;h Þ2 2 ðb  h  1Þ  8h ¼ ðh þ x2;h Þ

y3 ¼

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 h  3hx2;h þ x22;h þ x32;h þ ð1 þ x2;h Þ ðx42;h  6hx22;h þ h Þ 4hx2;h ð6:12Þ

¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 h þ x42  6hx2

; 2ðh þ x2 Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 x22 þ 2x2 h  h  h þ x42  6hx2 2ðh þ x2 Þ

;

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 2 ð1 þ x2 Þðx32 þ hx2  3x2 h þ h  ðx2 þ hÞ h þ x42  6hx2 Þ 4ðh þ x2 Þx2

;

y4 ¼ y3 :

and Condition (6.10) makes the radicand positive. (ii) Write s⁄ in terms of  x1;h ;  x2;h and h as s ¼

x22 þ 2x2  h þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  x1;h  3h þ x2;h þ x22;h þ ð1 þ x2;h Þ ðx22;h þ x21;h  6hÞ : 4h ð6:13Þ

Thus T4(s⁄,s⁄) lies along the diagonal y = x. Finally, to prove (iv), assume conditions (6.10) and (6.11) are satisfied. Because h ¼  x1;h  x2;h , then it is obvious that y1 > 0. Also, since

2x1;h < b () 2h < bx2;h () ðx2;h þ 1Þðh þ x2;h Þ () ðx22;h þ ð1 þ hÞx2;h  hÞ > 0 x22;h þ ð1 þ hÞ x2;h  hÞ > 0 is sufficient condition for and ð y2,y3 > 0, then the proof is complete. The next result shows that our choice of s⁄ is based on an 8-periodic solution of the case h = 1.

37

R. Abu-Saris et al. / Chaos, Solitons & Fractals 54 (2013) 26–38

pffiffiffi Theorem 6.4. Let h = 1 and assume that b > 2ð1 þ 2Þ. The ⁄ ⁄ initial conditions (y1,y0) = (s ,s ) define an 8-periodic solution of Eq. (6.8).

the positive factors, then substitute  x1;h  x2;h in place of h and factor again. Neglect the positive factors and collect the terms with a radical to obtain 3

Proof. The proof is computational. When h = 1, we obtain b = 1,

q1 ðtÞ ¼ q2 ðtÞ ¼

bð1  bÞtðt  s2 Þ ðb þ tÞ2

y3 ¼ y4 ¼ yconj 0 ;

ð6:15Þ ð6:16Þ

ðb þ tÞ2

ð6:17Þ

The next lemma is used in the proof of Theorem 6.5. Lemma 6.1. Assume that Condition (6.11) is satisfied, h 6 1 and s⁄ < s1. Then q2(s⁄) P q1(s⁄). Proof. First observe that s1  s2 ¼ 1b ð1  bÞðb  a2 Þ < 0 for  < s < s1 < s2 . Also, all h < 1. Thus, y

bða  bÞ3

; ða  1Þ2  x2;h x2;h þ h 2 > 0: q2 ð0Þ  q1 ð0Þ ¼ 2 x2;h þ 1 h ÞðAt 2 þ Bt þ CÞ ðb  1Þðt  y 2

ðt þ bÞ ðs1  tÞ

2 3

ð6:18Þ

Proof. Since the map T is one-to-one, then it is sufficient to show that T does not map the boundary of the closed curve outside. Since T

T

T

T

T

C0 ! C1 ! C2 ! C3 and C4 ! C5 ! C6 ! C7 ; 

C3 :¼ ðx; yÞ : y ¼

  a bs x ; f 1 ðs Þ 6 x 6 s ;   1þs 1þs

then ( TðC3 Þ ¼ ðx;yÞ : y ¼

) ða2  bðs þ 1ÞÞx  b2 s af1 ðs Þ  bs  : ; 6 x 6 f ðs Þ 1 a þ bs þ ðs þ 1Þx 1 þ s

Since T(C3) defines a concave function, then it is above C4. Also, since T(x,y) is below the diagonal for all points (x,y) such that y < x and x P  x2 , then T(C3) is within the region bounded by the curves Cj, j = 0, . . . ,8. Next, we check T(C8). From the one-to-one property of T and Proposition 5.2, it is sufficient to show that x1 6 x7 for all h 6 1. Observe that x1 6 x7 () g 2 ðf21 ðs ÞÞ 6

ða  bÞf2 ðg 1 ðf1 ðs ÞÞÞ 1 þ f2 ðg 1 ðf1 ðs ÞÞÞ

() ½ða  1Þ þ ðb  1Þg 1 ðf1 ðs ÞÞg 2 ðf21 ðs ÞÞ 6 ða  bÞbg 1 ðf1 ðs ÞÞ

Next, we write

q1 ðtÞ  q2 ðtÞ ¼

2

then we need to check T(C3) and T(C8). Since

ð1 þ tÞ2

ða  bÞb2 ðt  s2 Þ : ðt  s1 Þ

Þ ¼ q2 ðy Þ ¼ q1 ðy

3

2

Theorem 6.5. Let h 6 1 and assume that Condition (6.11) is satisfied. The region bounded by the closed curve defined by the line segments Cj, j = 0,1 . . . ,8 (C8 is a point when h = 1) defines an invariant region for Eq. (6.8). Moreover, when b < 1, removing the segment C8 makes the region a trapping region where orbits can enter through the removed segment and do not exit.

T

ða2  bÞ a b

þ ða  1Þ

3

which must be positive. A graphing utility shows that it is indeed positive when h < 1 and Condition (6.11) is satisfied. h

where ‘‘conj’’ is used to mean that we p replace A with  A. ffiffiffi Finally, the strict inequality b > 2ð1 þ 2Þ makes the radicand in A positive and the period a minimal period. h Consider ðx1 ; x0 Þ ¼ ðg 2 f21 ðs Þ; f21 ðs ÞÞ and define Cj to be the line segment connecting the point (xj1,xj) with (xj,xj+1) for all j = 0,1, . . . ,7. For j = 8 and b < 1, we define C8 to be the line segment connecting (x7,x8) with (x1,x0). Observe that when b = 1, (x7,x8) = (x1,x0), and consequently C8 shrinks to a point. Define

s1 :¼ a2  a  b and s2 ¼

2

 ð1  hÞð7h þ 11h  11h þ 1Þb þ 2hð3h  1Þð1  hÞ

and consequently

y2 ¼ yconj 1 ;

2

3

1 s ¼ ða3 þ 4  2a2  2a þ AÞ 4ða  1Þ

1 ða3  2a þ AÞ; 2a2 y5 ¼ y2 ; y6 ¼ y1 ;

2

þ ð3h þ 4h  3Þb þ ð3h þ 5h  13h þ 3Þb

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A ¼ a ða2  2Þða2  4a þ 2Þ and

y1 ¼

2

½ð1  hÞb þ 2ðh þ 2h  3Þb þ ð1  hÞðh  6h þ 1Þb qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 2 4 þ 2hð1  hÞ  ðb  h  1Þ  8h þ ð1  hÞb

() ½ða  1Þ þ ðb  1Þg 1 ðf1 ðs ÞÞ

;

s ðs1  s Þ ðb þ s Þ2

6 ða  bÞbg 1 ðf1 ðs ÞÞ;

where s1 as defined in (6.17). If s⁄ P s1, then it is an obvious case since the intersection of the region bounded by the nine-gone defined by Cj with the positive quadrant  < s < s1 and forms an invariance. So, we proceed with y obtain

where

 þ 2b; A :¼ ðb þ 1Þy ½bða þ 1Þðb þ 1Þ þ a þ b2 ða þ 2Þ þ ab; C :¼ y B :¼ a½2ða þ 1Þðb þ 1Þ  a2  3  b2 ða þ b þ 1Þ  b  1: 2

The task will be achieved if we show that p2(t):¼At + Bt + C is positive at t = s⁄. Write A,B and C in terms of  x2;h and h, then substitute s⁄ from Eq. (6.13) and factor p(s⁄). Neglect

½ð1 þ s Þ2 ða  1Þ þ ðb  1Þbs ðs2  s Þ

ðs1  s Þ ðb þ s Þ2

 2

ðb  1Þbs ðs2  sÞ þ ða  1Þð1 þ s Þ ðb þ s Þ2

6 ða  bÞb2 ðs2  s Þ

ð6:19Þ

2

6

ða  bÞb ðs2  s Þ ; ðs1  s Þ

ð6:20Þ

38

R. Abu-Saris et al. / Chaos, Solitons & Fractals 54 (2013) 26–38

where s2 as defined in (6.17). Now, Lemma 6.1 completes the proof. h Finally, we remark that when Condition (6.6) is satisfied, then all solutions inside the trapping region as defined above are attracted to the positive equilibrium a  b  1. However, when Condition (6.6) is not satisfied but h < 1, then we conjecture that the same scenario happens. Acknowledgment The authors thank an anonymous referee for his thorough reading and detailed suggestions. The second author would like to thank Eduardo Liz for the kind hospitality and stimulating discussion we had during a research visit to the University of Vigo. References [1] AlSharawi Z, Rhouma M. The Beverton–Holt model with periodic and conditional harvesting. J Biol Dyn 2009;3:463–78. [2] AlSharawi Z, Rhouma M. The discrete Beverton–Holt model with periodic harvesting in a periodically fluctuating environment. Adv Difference Equ 2010. http://dx.doi.org/10.1155/2010/215875. Article ID 215875. [3] AlSharawi Z, Rhouma M. Coexistence and extinction in a competitive exclusion Leslie/Gower model with harvesting and stocking. J Difference Equ Appl 2009;15:1031–53. [4] Barbeau E, Gelbord B, Tanny S. Periodicities of solutions of the generalized Lyness recursion. J Difference Equ Appl 1995;1:291–306. [5] Beverton R, Holt SJ. On the dynamics of exploited fish populations. New Jersey: The Blackburn Press; 2004. [6] Bischi GI, Lamantia F. Harvesting dynamics in protected and unprotected areas. J Econ Behaviour Org 2007;62:348–70. [7] Bischi GI, Lamantia F, Sbragia L. Strategic interaction and imitation dynamics in patch differentiated exploitation of fisheries. Ecol Complexity 2009;6:353–62. [8] Brannstrom A, Sumpter DJT. The role of competition and clustering in population dynamics. Proc R Soc B 2005;272:2065–72. [9] Camouzis E, Ladas G. Dynamics of third-order rational difference equations with open problems and conjectures. New York: Chapman & Hall; 2008. [10] Eskolaa H, Geritz S. On the mechanistic derivation of various discrete-time population models. Bull Math Biol 2007;69:329–46. [11] Feuer J, Janowski EJ, Ladas G. Lyness-type equations in the third quadrant. Nonlinear Anal 1997;30:1183–9. [12] Gardini L, Bischi GI, Mira C. Invariant curves and focal points in a Lyness iterative process. Int J Bifurcation Chaos 2003;13:1841–52. [13] Geritz S, Kisdi E. On the mechanistic underpinning of discrete-time population models with complex dynamics. J Theor Biol 2004;228:261–9.

[14] Getz WM, Haight RG. Population harvesting. New Jersey: Princeton University Press; 1989. [15] Gu En-Guo, Tian Fan. Complex dynamics analysis for a duopoly model of common fishery resource. Nonlinear Dyn 2010;61:579–90. [16] Hutson V, Moran W. Persistence of species obeying difference equations. J Math Biol 1982;15:203–13. [17] Kocic VL, Ladas G. Global asymptotic behavior of nonlinear difference equations of higher order with applications. Dordrecht: Kluwer Academic Publishers; 1993. [18] Kocic VL, Ladas G, Rodrigues IW. On rational recursive sequences. J Math Anal Appl 1993;173:127–57. [19] Kot M. Elements of mathematical ecology. UK: Cambridge University Press; 2003. [20] Kulenovic MRS, Ladas G. Dynamics of second order rational difference equations; with open problems and applications. Chapman & Hall/CRC Press; 2001. [21] Kuruklis SA, Ladas G. Oscillations and global attractivity in a discrete delay logistic model. Quart Appl Math 1992;1:227–33. [22] Leven SA, May RM. A note on difference delay equations. Theor Popul Biol 1976;9:178–87. [23] Levhari D, Mirman LJ. The great fish war: an example using a dynamics Gournot–Nash solution. Bell J Econ 1980;11:322–34. [24] Lyness RC. Note 1581. Math Gazette 1942;26:62. [25] Lyness RC. Note 1847. Math Gazette 1945;29:231–3. [26] Lyness RC. Note 2952. Math Gazette 1961;45:207–9. [27] May RM. Simple mathematical models with very complicated dynamics. Nature 1976;261:459–67. [28] Pielou EC. Population and community ecology. New York: Gordon and Breach; 1974. [29] Ricker WE. Stock and recruitment. J Fisheries Res Board Can 1954;11:559–623. [30] Beddington JR, May RM. harvesting natural population in a randomly fluctuating environment. Science 1977;197:463–5. [31] Brauer F, Sanchez DA. Constant rate population harvesting: equilibrium and stability. Theor Popul Biol 1975;8:12–30. [32] Brauer F, Sanchez DA. Periodic environments and periodic harvesting. Nat Res Model 2003;16:233–44. [33] Clark CW. Mathematical bioeconomics. New York, NY, USA: John Wiley and Sons Inc; 1990. [34] Fan M, Wang K. Optimal harvesting policy for single population with periodic coefficients. Math Biosci 1998;152:165–78. [35] Tang S, Chen L. The effect of seasonal harvesting on stage-structured population models. J Math Biol 2004;48:357–74. [36] Xu C, Boyce MS, Daley DJ. Harvesting in seasonal environments. J Math Biol 2005;50:663–82. [37] Zhang X, Shuai Z, Wang K. Optimal impulsive harvesting policy for single population. Nonlinear Anal: Real World Appl 2003;40:639–51. [38] Schriber S. Chaos and population disappearances in simple ecological models. J Math Biol 2001;42:239–60. [39] Walters CJ, Korman J. Linking recruitment to trophic factors: revisiting the Beverton–Holt recruitment model from a life history and multispecies perspective. Rev Fish Biol Fisheries 1999;9:187–202.