- Email: [email protected]
Controlling predator–prey discrete dynamics utilizing a threshold policy with hysteresis
Applied Mathematics and Computation 217 (2011) 7874–7886
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Controlling predator–prey discrete dynamics utilizing a threshold policy with hysteresis Magno Enrique Mendoza Meza a,⇑, Amit Bhaya b a
Center of Engineering, Modeling and Applied Social Sciences (CECS), Federal University of ABC (UFABC), Rua Catequese 242, 5 Andar, SP 09090-400, Jardim, Santo André, São Paulo, Brazil b Dept. of Electrical Engineering, COPPE/NACAD, Federal University of Rio de Janeiro, P.O. Box 68516, RJ 21941-972, Brazil
a r t i c l e
i n f o
Keywords: Discrete models Nonstandard scheme discretization Threshold policy with hysteresis
a b s t r a c t This paper introduces a threshold policy with hysteresis (TPH) for the control of the logistic one-species model, the Lotka–Volterra and Rosenzweig–MacArthur two species densitydependent predator–prey models. A nonstandard scheme is used for the discretization of the models since it results in preservation of the qualitative characteristics of the continuous-time models. Two theorems that establish the global stability of the discrete logistic model subject to the threshold policy (TP) and the TPH are proved. The proposed policy (TPH) is more realistic than a pure threshold policy (TP) proposed earlier in the literature and changes the dynamics of the system in such a way that a low amplitude bounded oscillation, far from the extinction region, is achieved. Furthermore, it can be designed by a suitable choice of so called virtual equilibrium points in a simple and intuitive manner. 2011 Elsevier Inc. All rights reserved.
1. Introduction In population dynamics, there are two kinds of mathematical models: continuous-time models, described by differential equations or dynamical systems and the discrete-time models, described by difference equations, discrete dynamical systems or iterative maps [1]. Investigations have shown that many economic, physical and biological phenomena are best represented via difference equations. However, there are many situations for which continuous models, i.e., differential equations, are the best fit. One may then invoke some numerical scheme to transform the given differential equation into a difference equation. The resulting difference equation is said to be dynamically consistent with its continuous counterpart if both exhibit the same qualitative behavior, such as stability, bifurcation and chaos [2]. A major difficulty in the numerical integration of ordinary differential equations (ODE) is the existence of numerical instabilities [3]. These are solutions to the discrete equations that do not correspond to any of the solutions to the original ODEs. Conventional numerical techniques, such as the forward Euler and higher-order Runge–Kutta schemes, are extensively used in solving systems of nonlinear differential equations. However, these numerical procedures may lead to chaotic or spurious solutions which are strongly scheme dependent. A nonstandard scheme was introduced by Mickens [3] to alleviate these problems and lead to difference equations that preserve the qualitative behavior of their continuous counterparts. The effects of harvesting on dynamic behavior in discrete-time population models have received increased attention [4–6, and references therein]. This paper studies the response of three discrete-time predator–prey models to a sustained perturbation (removal of species through harvesting). The contribution of this paper is to show that a continuous-time harvesting ⇑ Corresponding author. E-mail addresses: [email protected] (M.E. Mendoza Meza), [email protected] (A. Bhaya). 0096-3003/$ - see front matter 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2011.02.082
7875
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
or control policy [7] called a threshold policy with hysteresis (TPH) can be used successfully in the discrete-time case, provided that nonstandard discretization referred to above is used, instead of the standard discretization schemes. In this paper, we use a nonstandard scheme [3] to discretize the logistic model, the Lotka–Volterra model as well as the Rosenzweig–MacArthur model, which are subjected to TPH. 1.1. Preliminary definitions The definition of the real and virtual equilibria concepts, as well as the definition of the threshold policy (TP) and the TPH, are briefly summarized in this section for the following discrete time system
zkþ1 ¼ f ðzk ; tÞ u;
ð1Þ
T
where zk = [xk yk] is the state vector and the control policy u can be defined as the TP uTP = ezkw(s) or as the TPH uTPH = ezkwhys(s), where e is the effort (harvesting effort), w(s) is defined as follows:
wðsÞ ¼
1; if 0;
if
s > 0; s < 0;
ð2Þ
where s is the threshold that should be chosen according to the problem to be solved (Fig. 1(a)). A graphical representation of the uTP is shown in Fig. 1(b). The threshold policy with hysteresis (TPH) is defined as follows:
whys ðsÞ ¼
0;
if
1; if
s < r; or s 6 r and sðk þ 1Þ > sðkÞ; s > r; or s P r and sðk þ 1Þ < sðkÞ;
ð3Þ
where s is the threshold that should be chosen adequately, depending on the problem to be solved, r > 0 is the hysteresis parameter and s(k + 1), s(k) are the values of s at instants k + 1 and k (Fig. 2(a)). A graphical representation of the TPH can be seen in Fig. 2(b). i
i
Definition 1. Let zG be the stable equilibrium point of the dynamics of region Gi for i = {1, 2}. Then zG is called a real equilibrium if it belongs to Gi and a virtual equilibrium if it belongs to Gj, j – i, for i = {1, 2}. Remark 1. In Fig. 2(b), let the switching curves be denoted as s(xk) = r and s(xk) = r and virtual equilibrium points denoted 1 2 1 as zG and zG . Considering, for instance, a trajectory starting in G1 and seeking a stable equilibrium zG located in G2. This G1 trajectory will never attain z , because when it crosses the line s(xk) = r, the dynamics change and the trajectory now 2 2 seeks the stable equilibrium zG located in G1. Once again, it will never attain zG , because when the trajectory crosses the 1 switching line s(xk) = r the dynamics change, and the trajectory once again seeks the stable equilibrium zG and so on. 2. Logistic discrete model subject to TP and TPH The one-species discrete models are obtained according to nonstandard discretization method [3] by making the followx x ing replacements: (i) dx ! kþ1/ k where / is such that /(h) = h + O(h2) and h = Dt, e.g., /(h) = eh 1, (ii) x ? xk, (iii) x2 ? xk+1xk. dt The resulting one-species discrete models are shown in Table 1 where r is the intrinsic growth rate, K is the carrying capacity,
(a)
(b) 1
2
Fig. 1. (a) Graphical representation of w(s). (b) Graphical representation of the TP in a phase plane. zG and zG are the stable equilibrium points of the i i dynamics in region G1 and G2, respectively, and both are virtual. f G is the vector field in region Gi for i = {1, 2}. uG is the control applied to the system in i region G for i = {1, 2}.
7876
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
(a)
(b) 1
Fig. 2. (a) Graphical representation of whys(s). (b) Graphical representation of the TPH in a phase plane. The gray region is the hysteresis region G3. zG and 2 i i zG are the stable equilibrium points of the dynamics in region G1 and G2, respectively, and both are virtual. f G is the vector field in region Gi for i = {1, 2}. uG is the control applied to the system in region Gi for i = {1, 2}.
Table 1 Discrete logistic models with different consumption terms and their nonstandard discretizations. Continuous model x_ ¼ rx 1 Kx cx x_ ¼ rx 1 Kx xþd x_ ¼ rx 1 Kx
c 2 x2 2 x2 þd
Discrete model k xkþ1 ¼ ð1þr/Þx 1þ r /x K k
xkþ1 ¼ xkþ1 ¼
1þr/x c/þd xk k 1þKr /xk
ð1þr/Þxk 1þ
r c2 K þx2 þd2 k
/xk
c is the endogenous consumption rate and d relates to the prey (x) density at which predator satiation occurs. The discrete logistic model has been studied in [3], whereas the discrete logistic model with endogenous consumption Holling Type II and Holling Type III and nonstandard discretization are introduced in this paper. Another nonstandard discretization scheme for the models above can be seen in Appendix A. As an application example we take the continuous logistic model, which although it assumes a drastic simplification of the real-world complexities inherent in predator–prey systems, can still provide useful insight into the underlying behavior of the system. This type of model has been studied by [8–11] in the continuous-time case with an emphasis on the productive aspects of renewable-resource management. In the simplest form of the model, changes in the species abundance are described by
dx x ¼ rx 1 ; dt K
ð4Þ
where x denotes the density population, r is the intrinsic growth rate and K is the carrying capacity. Under nonstandard discretization [3] model (4) becomes
xkþ1 ¼
ð1 þ r/Þxk : 1 þ Kr /xk
ð5Þ
2.1. Logistic discrete model subject to TP The model (5) under TP is as follows:
xkþ1 ¼
ð1 þ r/Þxk uTP ; 1 þ Kr /xk
where the control is given by
uTP ¼ exk wðsÞ and w(s) is defined as in (2), where s is the threshold and e is an effort parameter. The simplest choice is
s ¼ xk xth ; where xth is the threshold level of population density.
ð6Þ
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
7877
2.1.1. Equilibrium points under TP The equilibrium point of the system without harvesting (without control) is xVNC = K and the equilibrium point of the sysr/e tem with harvesting (proportional control1) is xVC ¼ K r/ð1þ eÞ. If the effort parameter satisfies the restriction e < r/, then the equilibrium point xVC is positive. If e does not satisfy this restriction, then the system with control goes to extinction. In this paper, we consider values of e slightly greater than r/, i.e., e P r/. We choose the threshold level such that the equilibria xVNC and xVC are virtual, i.e.,
xVC < xth < xVNC : Note that the subscript VNC refers to virtual equilibrium of the no control region, and VC refers to virtual equilibrium of the control region. 2.2. Logistic discrete model subject to TPH Using the nonstandard discretization method [3] the logistic model (4) under TPH is as follows:
xkþ1 ¼
ð1 þ r/Þxk uTPH 1 þ Kr /xk
ð7Þ
where the control is given by
uTPH ¼ exk whys ðsÞ and whys(s) is defined as in (3), where s is the threshold that should be chosen according to the problem to be solved, and e is an effort parameter. The simplest choice of s is
s ¼ xk xth ; where xth is the threshold level of population density. We chose the threshold level such that
xVC < xth r < xth < xth þ r < xVNC : Simulation of model (6) under TP is shown in Fig. 3 and simulation of model (7) under TPH is shown in Fig. 4. Theorem 1. Consider a system of the following type
xkþ1 ¼ f ðxk ; /Þ uTP ; uTP ¼ exk wðsÞ;
ð8Þ
where w(s) is defined as in Eq. (2). Assume that: 1. the function f(xk, /) xk is nonnegative on the interval [0, xVNC), in which xVNC is the equilibrium point of the system without control, and 2. the effort parameter e satisfies
e<
1 : 1 þ r/
ð9Þ
Under these assumptions, the system stabilizes in a neighborhood of the threshold level xth, i.e., in an invariant interval, if the threshold level is chosen such that the equilibrium points of the system with control and without control are virtual. Proof 1. Details of the proof are given in Appendix B. The proof of Theorem 1 is similar to the proof of Theorem 2 and is referenced within the text of the proof of Theorem 2. h Theorem 2. Consider a system of the following type
(
xkþ1 ¼ f ðxk ; /Þ uTPH ; uTPH ¼ exk whys ðsÞ;
ð10Þ
where whys(s) is defined as in Eq. (3). Assume that: 1. the function f(xk, /) xk is nonnegative on the interval [0, xVNC), in which xVNC is the equilibrium point of the system without control, and 2. the effort parameter e satisfies
1
The proportional control is defined as u = exk.
7878
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
e<
1 : 1 þ r/
ð11Þ
Under these assumptions, the system stabilizes in the invariant interval (xL, xR), if the threshold levels xLL = xth r and xRR = xth + r are chosen such that the equilibrium points of the system with control and without control are virtual. Proof 2. Details of the proof are given in Appendix B. h
3. Lotka–Volterra model subject to TPH The Lotka–Volterra model is a simple model of a predator–prey interaction:
x_ ¼ r 1 x axy; y_ ¼ r 2 y þ bxy;
ð12Þ
where r1 is the specific growth rate for prey, r2 is the specific death rate for predator; a is a coefficient denoting the intensity of predation and b is a coefficient denoting the contribution to the growth of predators as a result of the predation. The state variable x denotes the prey population density and the state variable y denotes the predator population density. The Lotka– Volterra model has been of considerable value in posing highly relevant questions and is a jumping-off place for more realistic models; this is the main motivation for studying it here. The difficulty with obtaining the correct behavior of numerical solutions to discretizations of the Lotka–Volterra system is that all the trajectories in the relevant part of the phase-plane (x > 0, y > 0) are closed. Consequently, perturbations of the Lotka–Volterra system can change the formerly closed curves into another closed curves. In other words, this system is not structurally stable. A nonstandard finite-difference scheme, as generated from the rules given by Mickens [3,12], can be applied to a structurally unstable planar dynamical system, the Lotka–Volterra equations, to obtain numerical solutions having the correct behavior. The discretization of the Lotka–Volterra model (12) according to [2] makes the following replacements: in the first equax x y y tion of (12), (i) dx ! kþ1/ k , (ii) x ? xk, (iii) xy ? xk+1yk; and in the second equation of (12): (i) dy ! kþ1/ k , (ii) y ? yk, (iii) dt dt xy ? xkyk+1. The Lotka–Volterra discrete model is then as follows:
8 < xkþ1 ¼ ð1þr1 /1 ðhÞÞxk ; 1þa/1 ðhÞy k
ð13Þ
: ykþ1 ¼ ð1r2 /2 ðhÞÞyk ; 1b/2 ðhÞx k
where /1(h) = eh 1 and /2(h) = eh 1. The discrete time Lotka–Volterra model subject to a TPH is as follows:
8 < xkþ1 ¼ ð1þr1 /1 ðhÞÞxk e1 xk whys ðsÞ; 1þa/1 ðhÞy k
: ykþ1 ¼ ð1r2 /2 ðhÞÞyk e2 yk whys ðsÞ: 1b/2 ðhÞx
ð14Þ
k
The equilibrium point of the system without harvesting (without control) is
r
r1 b a
ðxVCN ; yVNC Þ ¼
2
;
and the equilibrium point of the system with harvesting (proportional control) is
ðxVC ; yVC Þ ¼
r 2 / þ e2 r 1 / e1 : ; b/ðe2 þ 1Þ a/ðe1 þ 1Þ
The simulation of the system (14) under TPH is shown in Fig. 5. Another way to discretize the Lotka–Volterra model (12), according to [12] is to make the following replacements in the x x y y first equation of (12): (i) dx ! kþ1/ k , (ii) x ? 2xk xk+1, (iii) xy ? xk+1yk; and in the second equation of (12): (i) dy ! kþ1/ k , (ii) dt dt y ? yk+1, (iii) xy ? 2xk+1yk xk+1yk+1. The resulting Lotka–Volterra discrete model subject to TPH is as follows:
(
1 /1 ðhÞÞxðkÞ xðk þ 1Þ ¼ 1þrð1þ2r e1 xðkÞwhys ðsÞ; 1 /1 ðhÞþa/1 ðhÞyðkÞ
ð1þ2b/1 ðhÞxðkþ1ÞÞyðkÞ e1 yðkÞwhys ðsÞ: xðk þ 1Þ ¼ 1þr 2 /1 ðhÞþb/1 ðhÞxðkþ1Þ
The equilibrium point of the system without harvesting (without control) is
zVNC ¼ ðxVNC ; yVNC Þ ¼
r r 2 1 ; b a
and the equilibrium point of the system with harvesting (proportional control) is
zVC ¼ ðxVC ; yVC Þ ¼
r 2 /ð1 þ e2 Þ þ e2 r 1 /ð1 e1 Þ e1 ; : b/ð1 e2 Þ a/ðe1 þ 1Þ
ð15Þ
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
7879
For both discretization schemes, the threshold is chosen from a weighted (or linear) combination of prey and predator densities, first proposed in [13],
s ¼ aa xk þ a2 yk Sd ; where a1 and a2 are attributed population weights, and Sd is a constant. The value of a1, a2, Sd, e1 and e2 are chosen such that zVC and zVNC are virtual and satisfied a slope condition defined in [13]. The simulation of the system (15) under TPH is shown in Fig. 6. Dynamics of the Lotka–Volterra discrete models (14) and (15) have a similar behavior under the TPH, as shown in Figs. 5 and 6. 4. Rosenzweig–MacArthur model The Rosenzweig–MacArthur model can be regarded as one of the simplest nontrivial paradigms that can be proposed by refining the more classical but biologically unrealistic Lotka–Volterra model. The Rosenzweig–MacArthur model is as follows:
(
xy _ xðtÞ ¼ rx 1 Kx xþA ; sAðxJÞ _ y; yðtÞ ¼ ðJþAÞðxþAÞ
ð16Þ
where r is the intrinsic growth rate of the prey, K is the carrying capacity of the environment, A is the half saturation constant, s conversion efficiency of predator, and J is the minimum prey population for which the predator can survive below the carrying capacity K; J < K. The state variable x denotes the prey population density and the state variable y denotes the predator population density. The discretization of the Rosenzweig–MacArthur model (16) is carried out according to [2]. The folx x lowing replacements are made in the first equation of (16): (i) dx ! kþ1/ k , (ii) x ? xk, (iii) x2 ? xk+1xk, (iii) xy ? xk+1yk; and in dt ykþ1 yk dy the second equation of (16): (i) dt ! / , (ii) y ? yk, (iii) xy ? xkyk+1. The resulting Rosenzweig–MacArthur discrete model is as follows:
8 k xkþ1 ¼ rð1þr/1 ðhÞÞx > / ðhÞy ; > 1þK /1 ðhÞxk 1x þA k > < k sAJ/2 ðhÞ 1 ðJþAÞx yk > > k > ; : ykþ1 ¼ sA/2 ðhÞxk 1ðJþAÞðx
ð17Þ
k þAÞ
where /1(h) = exp(h) 1 and /2(h) = exp(h) 1. The Rosenzweig–MacArthur model subject to a TPH is as follows:
8 k xkþ1 ¼ rð1þr/1 ðhÞÞx > / ðhÞy ; > 1þK /1 ðhÞxk 1x þA k > < k sAJ/2 ðhÞ 1 ðJþAÞx yk > > k > e2 yk whys ðsÞ: : ykþ1 ¼ sA/2 ðhÞxk 1ðJþAÞðx
ð18Þ
k þAÞ
The equilibrium point of the system without harvesting (without control) is
r zVNC ¼ ðxVNC ; yVNC Þ ¼ J; ðxVNC þ AÞðxVNC KÞ : K The equilibrium point of the system with harvesting (proportional control) is
zVC ¼ ðxVC ; yVC Þ ¼
AðsJ/ þ e2 ðJ þ AÞÞ r ; ðxVC þ AÞðxVC KÞ : e2 ðJ þ A sA/Þ þ sA/ K
sA/ sA/ If e2 < JþAsA/ the equilibrium point of the system is positive. We consider in this paper values of e2 > JþAsA/ . We choose the threshold level as follows:
s ¼ yk yth such that yVC < yth < yVNC, i.e., zVNC and zVC are virtual. Simulation of the model (18) under TPH is shown in Fig. 7. 5. Conclusions The policies TP and TPH have been shown to be appropriate for the control of the discrete logistic model, stabilizing the population density in a bounded region with positive values, i.e., avoiding extinction of species. The comparison between a nonstandard scheme and a standard scheme discretization for the logistic and Lotka–Volterra discrete models are shown in Fig. 8. A discrete logistic model with a standard scheme is
r xkþ1 ¼ xk 1 þ ðK xk Þ K
7880
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
which shows a behavior quite different from its continuous counterpart [3,14]. A Lotka–Volterra discrete model with a Euler’s scheme is as follows:
xkþ1 ¼ xk ð1 þ r 1 ayk Þ; ykþ1 ¼ yk ð1 r2 þ bxk Þ;
Fig. 3. Simulation of the discrete logistic model subject to a TP. Parameter values are: r = 1, K = 90, h = 0.01, / = 0.0101, xth = 45, e = 0.0201, initial conditions 15 and 75. The system stabilizes in a small neighborhood of xth.
Fig. 4. Simulation of the discrete logistic model subject to a TPH. Parameter values are: r = 1, K = 90, h = 0.01, / = 0.0101, xth = 45, e = 0.0201, r = 5, initial conditions 15 and 75. The system stabilizes in a small neighborhood around the hysteresis region.
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
7881
Fig. 5. Simulation of the discrete Lotka–Volterra model (14) subject to a TPH, see [2]. Parameter values a = 1, b = 1, r1 = 1, r2 = 1, e1 = 0.05, e2 = 0.05, h = 0.01, / = 0.01005, Sd = 1, r = 0.1, a1 = 0.2, a2 = 1. The system stabilizes in a low amplitude bounded oscillation.
Fig. 6. Simulation of the Lotka–Volterra model (15) subject to a TPH. Parameter values a = 1, b = 1, r1 = 1, r2 = 1, e1 = 0.5, e2 = 0.5, h = 0.01, / = 0.01005, Sd = 1, r = 0.1, a1 = 0.2, a2 = 1. The system stabilizes in a low amplitude bounded oscillation.
7882
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
Fig. 7. Simulation of the Rosenzweig–MacArthur model (18) subject to a TPH. Parameter values: r = 2, K = 70, s = 1, A = 10, J = 20, yth = 30, r = 5, / (h) = 0.1052, 2 = 0.05. Initial conditions: (10, 10), (50, 60), (70, 20), (30, 60), (70, 70). The system stabilizes in a low amplitude bounded oscillation.
which is plagued with deficiencies and lacks dynamical consistency. Another discrete version of the Lotka–Volterra model is
xkþ1 ¼ xk expðr 1 ayk Þ; ykþ1 ¼ yk expðr 2 þ bxk Þ
this system can show behavior quite different from that of the ordinary differential equation as can be seen in Fig. 8(b). The TPH allows a period of overexploitation, but avoids extinction of the species by an appropriate switching policy between periods of overexploitation and no exploitation. The important novel characteristic of a TPH is that it ensures that, even though the system is subjected to a period of overexploitation, it eventually stabilizes in low amplitude bounded oscillations in a desired safe region of the state space. In contrast, the commonly used constant harvesting policy (proportional control) cannot ensure this and will often lead to extinction while the system is being subjected to overexploitation. Under TP and TPH, the system is stabilized in a neighborhood of a point or in a hysteresis region, respectively, i.e., a bounded oscillation is achieved (Figs. 3–7). It is well known that the discrete logistic model with a standard discretization, e.g., forward Euler scheme, has solutions with various periods, as well as chaotic solutions [3,14]. The nonstandard scheme was introduced by [3] to alleviate these problems and lead to difference equations that preserve the qualitative behavior of their continuous counterparts. It must be stressed that when a threshold policy with hysteresis effect is considered, we are implicitly considering errors and delays in the implementation of the policy, i.e., errors and delays in the measurement of the species density. The hysteresis loop around the threshold level takes such occurrences into account explicitly and stabilizes the discrete logistic model in low amplitude bounded oscillations within the region delimited by the hysteresis itself. Finally, from the management standpoint, it is important to stress that the proposed strategy does not interfere directly in the harvesting intensities. Acknowledgement This research was partially financed by Project E-26/100.455/2007 (CNE) of FAPERJ, and Project No. 154631/2006-0 of CNPq, and also by the agency CAPES and FAPESP. We thank the reviewer for his constructive and detailed criticism and helpful suggestions that, in our opinion, have improved the manuscript.
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
7883
Appendix A. Another possibility for nonstandard discretization The one-species discrete models are obtained according to nonstandard discretization method [12], by making the folx x lowing replacements: (i) dx ! kþ1/ k where / is such that /(h) = h + O(h2) and h = Dt, e.g., /(h) = eh 1, (ii) dt x = 2x x ? 2xk xk+1, (iii) x2 ? xk+1xk. The resulting one-species discrete models are as follows: Continuous model x_ ¼ rx 1 Kx cx x_ ¼ rx 1 Kx xþd x_ ¼ rx 1 Kx
Discrete model ð1þ2r/Þxk xkþ1 ¼ 1þr/þ r K /xk
xkþ1 ¼
1þ2r/x2c/ xk þd k
1þr/þKr /xk x c/þd k
c 2 x2 x2 þd2
xkþ1 ¼
ð1þ2r/Þxk c2 /x 1þr/þKr /xk þ 2 k2 x þd k
Appendix B. Proof of the Theorems 2 and 1 Trajectories initiating in the intervals I1 = [0, xL], I3 = [xR, 1) (bold lines, Fig. B.1) will be shown to converge to the interval I2 = (xL, xR) which contains the points xLL = xth r and xRR = xth + r which determine the thresholds, in the sense that convergence occurs to a neighborhood of the interval (xLL, xRR). I2 is shown to be an invariant interval. The proof of Theorem 1 is similar to the proof of Theorem 2. For the proof of Theorem 1 is necessary to consider r = 0. The difference is that there are two different switching points (thresholds), one switching point for trajectories, which initiate in the interval [0, xLL] and switch when crossing the point xRR = xth + r and another switching point for trajectories, which initiate in the interval [xRR, 1) and switch when crossing the point xLL = xth r. All trajectories remain inside the interval I2. B.1. Equilibrium points under TPH and TP The system considered to prove the theorem is a logistic discrete model, which is subjected to the TPH,
xkþ1 ¼
ð1 þ r/Þxk uTPH ; 1 þ Kr /xk
where uTPH = exkwhys(s) and whys(s) is defined as in Eq. (3). The function
f ðxk ; /Þ xk ¼
ð1 þ r/Þxk xk 1 þ Kr /xk
is positive for values of xk < xVNC. The equilibrium point of the system without harvesting (without control) is
xVNC ¼ K and the equilibrium point of the system with harvesting (proportional control) is
xVC ¼ K
r/ e : r/ð1 þ eÞ
If the effort parameter satisfies the restriction e < r/ then the equilibrium point xVC is positive. If e does not satisfy this restriction, then the system without control goes to extinction. In this paper, we consider values of e slightly greater than r/, i.e.,
e P r/:
ðB:1Þ
We choose the threshold level such that the equilibria xVNC and xVC are virtual, see Fig. B.1, i.e.,
xVC < xLL < xth < xRR < xVNC : B.2. Determining the bounds of the globally attractive invariant interval We first calculate an interval that contains xLL = xth r and xRR = xth + r (the threshold values) by reasoning as follows. Suppose that the state xk ¼ xþ RR (i.e., infinitesimally to the left of xRR), then it is in the region without harvesting and the next state, which will denote xR, can be written as
xR ¼
ð1 þ r/ÞxRR : 1 þ Kr /xRR
For the proof of Theorem 1 is necessary to consider r = 0, then xRR = xth.
ðB:2Þ
7884
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
(a)
(b) Fig. 8. The comparison between a nonstandard scheme and the standard scheme discretization for the logistic and Lotka–Volterra discrete models. (a) Simulation of the discrete logistic model for a nonstandard and standard scheme discretization. Parameter values for both simulations are the same. r = 3, K = 1, initial condition 0.1. (b) Simulation of the Lotka–Volterra discrete model for a nonstandard and standard scheme discretization. Parameter values for both simulations are the same. a = 1, b = 1, r1 = 1, r2 = 1, initial condition (2.1, 2.1).
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
7885
Similarly, if xþ LL (i.e., infinitesimally to the right of xLL) and therefore in the region with harvesting, the next state, denoted as xL, can be calculated as
xL ¼
ð1 þ r/ÞxLL exLL : 1 þ Kr /xLL
ðB:3Þ
For the proof of Theorem 1 is necessary to consider r = 0, then xLL = xth. We will show that xL and xR are the left and right endpoints of an interval, denoted I2, that will be proved to be globally attractive and invariant under the TPH and TP. Eqs. (B.2) and (B.3) define the endpoints of the interval I2. B.3. Proof of the global attractivity and invariance of the interval I2 A suitable reaching phase Lyapunov function, adapted from [15], is
Vðxk Þ ¼ xk 2xth þ
x2th ðxk xth Þ2 ¼ xk xk
ðB:4Þ
for which
x2 DVðxk Þ ¼ Vðxkþ1 Þ Vðxk Þ ¼ ðxkþ1 xk Þ 1 th : xkþ1 xk
ðB:5Þ
This function is used to prove: B.3.1. Global attractivity of interval I2 We now use V, DV to carry out an interval analysis of the system (10). Interval I1: Assume xk 2 I1. If xk+1 2 I1, then xk+1 6 xth, and from the dynamics in this interval, we have
xkþ1 ¼
ð1 þ r/Þxk : 1 þ Kr /xk
ðB:6Þ
Since xk < xVNC = K, "k, the factor ð1 þ r/Þ= 1 þ Kr /xk is greater than 1, thus from (B.6), xk+1 > xk and xk+1 6 xth. Therefore, from (B.5), D V < 0 and all trajectories with initial condition in the interval I1 converges to the interval I2. Interval I3: Assume xk 2 I3. If xk+1 2 I3, then xk+1 P xth, and from the dynamics in this interval, we have
xkþ1
1 þ r/ e 1 þ Kr /xk xk ð1 þ r/Þxk ¼ exk ¼ : 1 þ Kr /xk 1 þ Kr /xk
K r/e . The latter holds by (11) since xk is a nonnegative The factor ð1 þ r/Þ= 1 þ Kr /xk e is smaller than 1 for xk > r/ 1þe variable. Thus xk+1 < xk and xk+1 P xth, "xk 2 I3. Therefore, DV < 0 and all trajectories with initial condition in the interval I3 converge to the interval I2. Therefore, trajectories initiating in the intervals I1 and I3 converge to the interval I2, proving that it is globally attractive. B.3.2. Invariance of interval I2 Inside the interval I2, the plan of proof is to show that points in I2 to the left of xRR cannot jump, under the dynamics (10), to the right of xR and, similarly, points to the right of xLL cannot jump to the left of xL. This completes the proof of invariance of interval I2.
Fig. B.1. Graphical representation of the phase interval. xVC and xVNC are the equilibrium point of the system with harvesting and without harvesting, respectively. xL and xR are the left bounded and the right bounded of the interval I2, respectively. xLL = xth r and xRR = xth + r are the threshold levels.
7886
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
Fig. B.2. Graphical representation inside the interval I2. xl is the initial condition in the region without harvesting. xr is the initial condition in the region with harvesting.
1. Consider an initial condition xk in the region without harvesting (Fig. B.2) that belongs to the interval (xL, xRR) and suppose that the next state is on the right of xR, i.e., xk+1 2 I3, and xk+1 > xR,
ð1 þ r/Þxk ð1 þ r/ÞxRR > : 1 þ Kr /xk 1 þ Kr /xRR
ðB:7Þ
Inequality (B.7) is satisfied if xk > xRR, which contradicts xk 2 I2. In other words, if xk 2 I2, then xk+1 is smaller than xR, i.e., the trajectory remains inside the interval I2. For the proof of Theorem 1 is necessary to consider r = 0, then xRR = xth. 2. Consider an initial condition xk in the region with harvesting (Fig. B.2) that belongs to the interval (xLL, xR) and suppose that next state is on the left of xL, i.e., xk+1 2 I1. Let a = 1 + r/, b = r//K for notational simplicity. Then xk+1 < xL can be written as
axk axLL exk < exLL ; 1 þ bxk 1 þ bxLL which implies
ðxr xLL Þða eð1 þ bxr Þð1 þ bxLL ÞÞ < 0:
ðB:8Þ
We will prove that the second term of (B.8) is positive. In fact,
a eð1 þ bxr Þð1 þ bxLL Þ > 0 ()
e<
1 þ r/ : x 1 þ r/ x 1 þ r/ K r K LL
By assumption 2 and for values of xk and xLL in the interval [0, xVNC], we have
e<
1 1 þ r/ < 1 þ r/: < 1 þ r/ x 1 þ r/ x 1 þ r/ K r K LL
Thus, Eq. (B.8) satisfied for xr < xLL, which contradicts xk 2 (xLL, xR). In other words, if xk 2 (xLL, xR), then xk+1 is greater than xL, i.e., the trajectory remains inside the interval I2. For the proof of Theorem 1 is necessary to consider r = 0, then xLL = xth. The proof of the analogs of Theorems 1 and 2 for the logistic discrete model with endogenous consumption of Holling Type II and Holling Type III, models in Table 1, can be carried out similarly. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]
X. Liu, D. Xiao, Complex dynamics behaviors of discrete-time predator–prey system, Chaos Solutions Fractals 32 (2007) 80–94. P. Liu, S.N. Elaydi, Discrete competitive and cooperative models of Lotka–Volterra type, J. Comput. Anal. Appl. 3 (1) (2001) 53–73. R.E. Mickens, Nonstandard Finite Difference Models of Differential Equations, World Scientific, Singapore, 1994. M. Basson, M.J. Fogarty, Harvesting in discrete-time predator–prey systems, Math. Biosci. 141 (1) (1997) 41–74. H. Kokko, J. Lindström, Seasonal density dependence, timing of mortality, and sustainable harvesting, Ecol. Model. 110 (3) (1998) 293–304. S.Y. Tang, L.S. Chen, The effect of seasonal harvesting on stage-structured population models, J. Math. Biol. 48 (4) (2004) 357–374. Magno Enrique Mendoza Meza, Amit Bhaya, Eugenius Kaszkurewicz, Michel Iskin Silveira Costa, On-off policy and hysteresis on-off policy control of the herbivore-vegetation dynamics in a semi-arid grazing system, Ecol. Eng. 28 (2) (2006) 114–123. I. Noy-Meir, Stability of grazing systems: an application of predator–prey graphs, J. Ecol. 63 (2) (1975) 459–481. C.W. Clark, Mathematical Bioeconomics: The Optimal Management of Renewable Resources, John Wiley & Sons, USA, 1976. C.W. Clark, Bioeconomics Modelling and Fisheries Management, John Wiley & Sons, USA, 1985. M. Kot, Elements of Mathematical Ecology, Cambridge University Press, Cambridge, 2001. R.E. Mickens, A nonstandard finite-difference scheme for the Lotka–Volterra system, Appl. Numer. Math. 45 (2003) 309–314. M.I.S. Costa, E. Kaszkurewicz, A. Bhaya, L. Hsu, Achieving global convergence to an equilibrium population in predator–prey systems by the use of discontinuous harvesting policy, Ecol. Model. 128 (2000) 89–99. J.D. Murray, Mathematical Biology I. An Introduction, Springer-Verlag, New York, 2002. B.-S. Goh, Management and analysis of biological populations, Developments in Agricultural and Managed-Forest Ecology, vol. 8, Elsevier, 1980.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Controlling predator–prey discrete dynamics utilizing a threshold policy with hysteresis Magno Enrique Mendoza Meza a,⇑, Amit Bhaya b a
Center of Engineering, Modeling and Applied Social Sciences (CECS), Federal University of ABC (UFABC), Rua Catequese 242, 5 Andar, SP 09090-400, Jardim, Santo André, São Paulo, Brazil b Dept. of Electrical Engineering, COPPE/NACAD, Federal University of Rio de Janeiro, P.O. Box 68516, RJ 21941-972, Brazil
a r t i c l e
i n f o
Keywords: Discrete models Nonstandard scheme discretization Threshold policy with hysteresis
a b s t r a c t This paper introduces a threshold policy with hysteresis (TPH) for the control of the logistic one-species model, the Lotka–Volterra and Rosenzweig–MacArthur two species densitydependent predator–prey models. A nonstandard scheme is used for the discretization of the models since it results in preservation of the qualitative characteristics of the continuous-time models. Two theorems that establish the global stability of the discrete logistic model subject to the threshold policy (TP) and the TPH are proved. The proposed policy (TPH) is more realistic than a pure threshold policy (TP) proposed earlier in the literature and changes the dynamics of the system in such a way that a low amplitude bounded oscillation, far from the extinction region, is achieved. Furthermore, it can be designed by a suitable choice of so called virtual equilibrium points in a simple and intuitive manner. 2011 Elsevier Inc. All rights reserved.
1. Introduction In population dynamics, there are two kinds of mathematical models: continuous-time models, described by differential equations or dynamical systems and the discrete-time models, described by difference equations, discrete dynamical systems or iterative maps [1]. Investigations have shown that many economic, physical and biological phenomena are best represented via difference equations. However, there are many situations for which continuous models, i.e., differential equations, are the best fit. One may then invoke some numerical scheme to transform the given differential equation into a difference equation. The resulting difference equation is said to be dynamically consistent with its continuous counterpart if both exhibit the same qualitative behavior, such as stability, bifurcation and chaos [2]. A major difficulty in the numerical integration of ordinary differential equations (ODE) is the existence of numerical instabilities [3]. These are solutions to the discrete equations that do not correspond to any of the solutions to the original ODEs. Conventional numerical techniques, such as the forward Euler and higher-order Runge–Kutta schemes, are extensively used in solving systems of nonlinear differential equations. However, these numerical procedures may lead to chaotic or spurious solutions which are strongly scheme dependent. A nonstandard scheme was introduced by Mickens [3] to alleviate these problems and lead to difference equations that preserve the qualitative behavior of their continuous counterparts. The effects of harvesting on dynamic behavior in discrete-time population models have received increased attention [4–6, and references therein]. This paper studies the response of three discrete-time predator–prey models to a sustained perturbation (removal of species through harvesting). The contribution of this paper is to show that a continuous-time harvesting ⇑ Corresponding author. E-mail addresses: [email protected] (M.E. Mendoza Meza), [email protected] (A. Bhaya). 0096-3003/$ - see front matter 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2011.02.082
7875
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
or control policy [7] called a threshold policy with hysteresis (TPH) can be used successfully in the discrete-time case, provided that nonstandard discretization referred to above is used, instead of the standard discretization schemes. In this paper, we use a nonstandard scheme [3] to discretize the logistic model, the Lotka–Volterra model as well as the Rosenzweig–MacArthur model, which are subjected to TPH. 1.1. Preliminary definitions The definition of the real and virtual equilibria concepts, as well as the definition of the threshold policy (TP) and the TPH, are briefly summarized in this section for the following discrete time system
zkþ1 ¼ f ðzk ; tÞ u;
ð1Þ
T
where zk = [xk yk] is the state vector and the control policy u can be defined as the TP uTP = ezkw(s) or as the TPH uTPH = ezkwhys(s), where e is the effort (harvesting effort), w(s) is defined as follows:
wðsÞ ¼
1; if 0;
if
s > 0; s < 0;
ð2Þ
where s is the threshold that should be chosen according to the problem to be solved (Fig. 1(a)). A graphical representation of the uTP is shown in Fig. 1(b). The threshold policy with hysteresis (TPH) is defined as follows:
whys ðsÞ ¼
0;
if
1; if
s < r; or s 6 r and sðk þ 1Þ > sðkÞ; s > r; or s P r and sðk þ 1Þ < sðkÞ;
ð3Þ
where s is the threshold that should be chosen adequately, depending on the problem to be solved, r > 0 is the hysteresis parameter and s(k + 1), s(k) are the values of s at instants k + 1 and k (Fig. 2(a)). A graphical representation of the TPH can be seen in Fig. 2(b). i
i
Definition 1. Let zG be the stable equilibrium point of the dynamics of region Gi for i = {1, 2}. Then zG is called a real equilibrium if it belongs to Gi and a virtual equilibrium if it belongs to Gj, j – i, for i = {1, 2}. Remark 1. In Fig. 2(b), let the switching curves be denoted as s(xk) = r and s(xk) = r and virtual equilibrium points denoted 1 2 1 as zG and zG . Considering, for instance, a trajectory starting in G1 and seeking a stable equilibrium zG located in G2. This G1 trajectory will never attain z , because when it crosses the line s(xk) = r, the dynamics change and the trajectory now 2 2 seeks the stable equilibrium zG located in G1. Once again, it will never attain zG , because when the trajectory crosses the 1 switching line s(xk) = r the dynamics change, and the trajectory once again seeks the stable equilibrium zG and so on. 2. Logistic discrete model subject to TP and TPH The one-species discrete models are obtained according to nonstandard discretization method [3] by making the followx x ing replacements: (i) dx ! kþ1/ k where / is such that /(h) = h + O(h2) and h = Dt, e.g., /(h) = eh 1, (ii) x ? xk, (iii) x2 ? xk+1xk. dt The resulting one-species discrete models are shown in Table 1 where r is the intrinsic growth rate, K is the carrying capacity,
(a)
(b) 1
2
Fig. 1. (a) Graphical representation of w(s). (b) Graphical representation of the TP in a phase plane. zG and zG are the stable equilibrium points of the i i dynamics in region G1 and G2, respectively, and both are virtual. f G is the vector field in region Gi for i = {1, 2}. uG is the control applied to the system in i region G for i = {1, 2}.
7876
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
(a)
(b) 1
Fig. 2. (a) Graphical representation of whys(s). (b) Graphical representation of the TPH in a phase plane. The gray region is the hysteresis region G3. zG and 2 i i zG are the stable equilibrium points of the dynamics in region G1 and G2, respectively, and both are virtual. f G is the vector field in region Gi for i = {1, 2}. uG is the control applied to the system in region Gi for i = {1, 2}.
Table 1 Discrete logistic models with different consumption terms and their nonstandard discretizations. Continuous model x_ ¼ rx 1 Kx cx x_ ¼ rx 1 Kx xþd x_ ¼ rx 1 Kx
c 2 x2 2 x2 þd
Discrete model k xkþ1 ¼ ð1þr/Þx 1þ r /x K k
xkþ1 ¼ xkþ1 ¼
1þr/x c/þd xk k 1þKr /xk
ð1þr/Þxk 1þ
r c2 K þx2 þd2 k
/xk
c is the endogenous consumption rate and d relates to the prey (x) density at which predator satiation occurs. The discrete logistic model has been studied in [3], whereas the discrete logistic model with endogenous consumption Holling Type II and Holling Type III and nonstandard discretization are introduced in this paper. Another nonstandard discretization scheme for the models above can be seen in Appendix A. As an application example we take the continuous logistic model, which although it assumes a drastic simplification of the real-world complexities inherent in predator–prey systems, can still provide useful insight into the underlying behavior of the system. This type of model has been studied by [8–11] in the continuous-time case with an emphasis on the productive aspects of renewable-resource management. In the simplest form of the model, changes in the species abundance are described by
dx x ¼ rx 1 ; dt K
ð4Þ
where x denotes the density population, r is the intrinsic growth rate and K is the carrying capacity. Under nonstandard discretization [3] model (4) becomes
xkþ1 ¼
ð1 þ r/Þxk : 1 þ Kr /xk
ð5Þ
2.1. Logistic discrete model subject to TP The model (5) under TP is as follows:
xkþ1 ¼
ð1 þ r/Þxk uTP ; 1 þ Kr /xk
where the control is given by
uTP ¼ exk wðsÞ and w(s) is defined as in (2), where s is the threshold and e is an effort parameter. The simplest choice is
s ¼ xk xth ; where xth is the threshold level of population density.
ð6Þ
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
7877
2.1.1. Equilibrium points under TP The equilibrium point of the system without harvesting (without control) is xVNC = K and the equilibrium point of the sysr/e tem with harvesting (proportional control1) is xVC ¼ K r/ð1þ eÞ. If the effort parameter satisfies the restriction e < r/, then the equilibrium point xVC is positive. If e does not satisfy this restriction, then the system with control goes to extinction. In this paper, we consider values of e slightly greater than r/, i.e., e P r/. We choose the threshold level such that the equilibria xVNC and xVC are virtual, i.e.,
xVC < xth < xVNC : Note that the subscript VNC refers to virtual equilibrium of the no control region, and VC refers to virtual equilibrium of the control region. 2.2. Logistic discrete model subject to TPH Using the nonstandard discretization method [3] the logistic model (4) under TPH is as follows:
xkþ1 ¼
ð1 þ r/Þxk uTPH 1 þ Kr /xk
ð7Þ
where the control is given by
uTPH ¼ exk whys ðsÞ and whys(s) is defined as in (3), where s is the threshold that should be chosen according to the problem to be solved, and e is an effort parameter. The simplest choice of s is
s ¼ xk xth ; where xth is the threshold level of population density. We chose the threshold level such that
xVC < xth r < xth < xth þ r < xVNC : Simulation of model (6) under TP is shown in Fig. 3 and simulation of model (7) under TPH is shown in Fig. 4. Theorem 1. Consider a system of the following type
xkþ1 ¼ f ðxk ; /Þ uTP ; uTP ¼ exk wðsÞ;
ð8Þ
where w(s) is defined as in Eq. (2). Assume that: 1. the function f(xk, /) xk is nonnegative on the interval [0, xVNC), in which xVNC is the equilibrium point of the system without control, and 2. the effort parameter e satisfies
e<
1 : 1 þ r/
ð9Þ
Under these assumptions, the system stabilizes in a neighborhood of the threshold level xth, i.e., in an invariant interval, if the threshold level is chosen such that the equilibrium points of the system with control and without control are virtual. Proof 1. Details of the proof are given in Appendix B. The proof of Theorem 1 is similar to the proof of Theorem 2 and is referenced within the text of the proof of Theorem 2. h Theorem 2. Consider a system of the following type
(
xkþ1 ¼ f ðxk ; /Þ uTPH ; uTPH ¼ exk whys ðsÞ;
ð10Þ
where whys(s) is defined as in Eq. (3). Assume that: 1. the function f(xk, /) xk is nonnegative on the interval [0, xVNC), in which xVNC is the equilibrium point of the system without control, and 2. the effort parameter e satisfies
1
The proportional control is defined as u = exk.
7878
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
e<
1 : 1 þ r/
ð11Þ
Under these assumptions, the system stabilizes in the invariant interval (xL, xR), if the threshold levels xLL = xth r and xRR = xth + r are chosen such that the equilibrium points of the system with control and without control are virtual. Proof 2. Details of the proof are given in Appendix B. h
3. Lotka–Volterra model subject to TPH The Lotka–Volterra model is a simple model of a predator–prey interaction:
x_ ¼ r 1 x axy; y_ ¼ r 2 y þ bxy;
ð12Þ
where r1 is the specific growth rate for prey, r2 is the specific death rate for predator; a is a coefficient denoting the intensity of predation and b is a coefficient denoting the contribution to the growth of predators as a result of the predation. The state variable x denotes the prey population density and the state variable y denotes the predator population density. The Lotka– Volterra model has been of considerable value in posing highly relevant questions and is a jumping-off place for more realistic models; this is the main motivation for studying it here. The difficulty with obtaining the correct behavior of numerical solutions to discretizations of the Lotka–Volterra system is that all the trajectories in the relevant part of the phase-plane (x > 0, y > 0) are closed. Consequently, perturbations of the Lotka–Volterra system can change the formerly closed curves into another closed curves. In other words, this system is not structurally stable. A nonstandard finite-difference scheme, as generated from the rules given by Mickens [3,12], can be applied to a structurally unstable planar dynamical system, the Lotka–Volterra equations, to obtain numerical solutions having the correct behavior. The discretization of the Lotka–Volterra model (12) according to [2] makes the following replacements: in the first equax x y y tion of (12), (i) dx ! kþ1/ k , (ii) x ? xk, (iii) xy ? xk+1yk; and in the second equation of (12): (i) dy ! kþ1/ k , (ii) y ? yk, (iii) dt dt xy ? xkyk+1. The Lotka–Volterra discrete model is then as follows:
8 < xkþ1 ¼ ð1þr1 /1 ðhÞÞxk ; 1þa/1 ðhÞy k
ð13Þ
: ykþ1 ¼ ð1r2 /2 ðhÞÞyk ; 1b/2 ðhÞx k
where /1(h) = eh 1 and /2(h) = eh 1. The discrete time Lotka–Volterra model subject to a TPH is as follows:
8 < xkþ1 ¼ ð1þr1 /1 ðhÞÞxk e1 xk whys ðsÞ; 1þa/1 ðhÞy k
: ykþ1 ¼ ð1r2 /2 ðhÞÞyk e2 yk whys ðsÞ: 1b/2 ðhÞx
ð14Þ
k
The equilibrium point of the system without harvesting (without control) is
r
r1 b a
ðxVCN ; yVNC Þ ¼
2
;
and the equilibrium point of the system with harvesting (proportional control) is
ðxVC ; yVC Þ ¼
r 2 / þ e2 r 1 / e1 : ; b/ðe2 þ 1Þ a/ðe1 þ 1Þ
The simulation of the system (14) under TPH is shown in Fig. 5. Another way to discretize the Lotka–Volterra model (12), according to [12] is to make the following replacements in the x x y y first equation of (12): (i) dx ! kþ1/ k , (ii) x ? 2xk xk+1, (iii) xy ? xk+1yk; and in the second equation of (12): (i) dy ! kþ1/ k , (ii) dt dt y ? yk+1, (iii) xy ? 2xk+1yk xk+1yk+1. The resulting Lotka–Volterra discrete model subject to TPH is as follows:
(
1 /1 ðhÞÞxðkÞ xðk þ 1Þ ¼ 1þrð1þ2r e1 xðkÞwhys ðsÞ; 1 /1 ðhÞþa/1 ðhÞyðkÞ
ð1þ2b/1 ðhÞxðkþ1ÞÞyðkÞ e1 yðkÞwhys ðsÞ: xðk þ 1Þ ¼ 1þr 2 /1 ðhÞþb/1 ðhÞxðkþ1Þ
The equilibrium point of the system without harvesting (without control) is
zVNC ¼ ðxVNC ; yVNC Þ ¼
r r 2 1 ; b a
and the equilibrium point of the system with harvesting (proportional control) is
zVC ¼ ðxVC ; yVC Þ ¼
r 2 /ð1 þ e2 Þ þ e2 r 1 /ð1 e1 Þ e1 ; : b/ð1 e2 Þ a/ðe1 þ 1Þ
ð15Þ
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
7879
For both discretization schemes, the threshold is chosen from a weighted (or linear) combination of prey and predator densities, first proposed in [13],
s ¼ aa xk þ a2 yk Sd ; where a1 and a2 are attributed population weights, and Sd is a constant. The value of a1, a2, Sd, e1 and e2 are chosen such that zVC and zVNC are virtual and satisfied a slope condition defined in [13]. The simulation of the system (15) under TPH is shown in Fig. 6. Dynamics of the Lotka–Volterra discrete models (14) and (15) have a similar behavior under the TPH, as shown in Figs. 5 and 6. 4. Rosenzweig–MacArthur model The Rosenzweig–MacArthur model can be regarded as one of the simplest nontrivial paradigms that can be proposed by refining the more classical but biologically unrealistic Lotka–Volterra model. The Rosenzweig–MacArthur model is as follows:
(
xy _ xðtÞ ¼ rx 1 Kx xþA ; sAðxJÞ _ y; yðtÞ ¼ ðJþAÞðxþAÞ
ð16Þ
where r is the intrinsic growth rate of the prey, K is the carrying capacity of the environment, A is the half saturation constant, s conversion efficiency of predator, and J is the minimum prey population for which the predator can survive below the carrying capacity K; J < K. The state variable x denotes the prey population density and the state variable y denotes the predator population density. The discretization of the Rosenzweig–MacArthur model (16) is carried out according to [2]. The folx x lowing replacements are made in the first equation of (16): (i) dx ! kþ1/ k , (ii) x ? xk, (iii) x2 ? xk+1xk, (iii) xy ? xk+1yk; and in dt ykþ1 yk dy the second equation of (16): (i) dt ! / , (ii) y ? yk, (iii) xy ? xkyk+1. The resulting Rosenzweig–MacArthur discrete model is as follows:
8 k xkþ1 ¼ rð1þr/1 ðhÞÞx > / ðhÞy ; > 1þK /1 ðhÞxk 1x þA k > < k sAJ/2 ðhÞ 1 ðJþAÞx yk > > k > ; : ykþ1 ¼ sA/2 ðhÞxk 1ðJþAÞðx
ð17Þ
k þAÞ
where /1(h) = exp(h) 1 and /2(h) = exp(h) 1. The Rosenzweig–MacArthur model subject to a TPH is as follows:
8 k xkþ1 ¼ rð1þr/1 ðhÞÞx > / ðhÞy ; > 1þK /1 ðhÞxk 1x þA k > < k sAJ/2 ðhÞ 1 ðJþAÞx yk > > k > e2 yk whys ðsÞ: : ykþ1 ¼ sA/2 ðhÞxk 1ðJþAÞðx
ð18Þ
k þAÞ
The equilibrium point of the system without harvesting (without control) is
r zVNC ¼ ðxVNC ; yVNC Þ ¼ J; ðxVNC þ AÞðxVNC KÞ : K The equilibrium point of the system with harvesting (proportional control) is
zVC ¼ ðxVC ; yVC Þ ¼
AðsJ/ þ e2 ðJ þ AÞÞ r ; ðxVC þ AÞðxVC KÞ : e2 ðJ þ A sA/Þ þ sA/ K
sA/ sA/ If e2 < JþAsA/ the equilibrium point of the system is positive. We consider in this paper values of e2 > JþAsA/ . We choose the threshold level as follows:
s ¼ yk yth such that yVC < yth < yVNC, i.e., zVNC and zVC are virtual. Simulation of the model (18) under TPH is shown in Fig. 7. 5. Conclusions The policies TP and TPH have been shown to be appropriate for the control of the discrete logistic model, stabilizing the population density in a bounded region with positive values, i.e., avoiding extinction of species. The comparison between a nonstandard scheme and a standard scheme discretization for the logistic and Lotka–Volterra discrete models are shown in Fig. 8. A discrete logistic model with a standard scheme is
r xkþ1 ¼ xk 1 þ ðK xk Þ K
7880
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
which shows a behavior quite different from its continuous counterpart [3,14]. A Lotka–Volterra discrete model with a Euler’s scheme is as follows:
xkþ1 ¼ xk ð1 þ r 1 ayk Þ; ykþ1 ¼ yk ð1 r2 þ bxk Þ;
Fig. 3. Simulation of the discrete logistic model subject to a TP. Parameter values are: r = 1, K = 90, h = 0.01, / = 0.0101, xth = 45, e = 0.0201, initial conditions 15 and 75. The system stabilizes in a small neighborhood of xth.
Fig. 4. Simulation of the discrete logistic model subject to a TPH. Parameter values are: r = 1, K = 90, h = 0.01, / = 0.0101, xth = 45, e = 0.0201, r = 5, initial conditions 15 and 75. The system stabilizes in a small neighborhood around the hysteresis region.
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
7881
Fig. 5. Simulation of the discrete Lotka–Volterra model (14) subject to a TPH, see [2]. Parameter values a = 1, b = 1, r1 = 1, r2 = 1, e1 = 0.05, e2 = 0.05, h = 0.01, / = 0.01005, Sd = 1, r = 0.1, a1 = 0.2, a2 = 1. The system stabilizes in a low amplitude bounded oscillation.
Fig. 6. Simulation of the Lotka–Volterra model (15) subject to a TPH. Parameter values a = 1, b = 1, r1 = 1, r2 = 1, e1 = 0.5, e2 = 0.5, h = 0.01, / = 0.01005, Sd = 1, r = 0.1, a1 = 0.2, a2 = 1. The system stabilizes in a low amplitude bounded oscillation.
7882
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
Fig. 7. Simulation of the Rosenzweig–MacArthur model (18) subject to a TPH. Parameter values: r = 2, K = 70, s = 1, A = 10, J = 20, yth = 30, r = 5, / (h) = 0.1052, 2 = 0.05. Initial conditions: (10, 10), (50, 60), (70, 20), (30, 60), (70, 70). The system stabilizes in a low amplitude bounded oscillation.
which is plagued with deficiencies and lacks dynamical consistency. Another discrete version of the Lotka–Volterra model is
xkþ1 ¼ xk expðr 1 ayk Þ; ykþ1 ¼ yk expðr 2 þ bxk Þ
this system can show behavior quite different from that of the ordinary differential equation as can be seen in Fig. 8(b). The TPH allows a period of overexploitation, but avoids extinction of the species by an appropriate switching policy between periods of overexploitation and no exploitation. The important novel characteristic of a TPH is that it ensures that, even though the system is subjected to a period of overexploitation, it eventually stabilizes in low amplitude bounded oscillations in a desired safe region of the state space. In contrast, the commonly used constant harvesting policy (proportional control) cannot ensure this and will often lead to extinction while the system is being subjected to overexploitation. Under TP and TPH, the system is stabilized in a neighborhood of a point or in a hysteresis region, respectively, i.e., a bounded oscillation is achieved (Figs. 3–7). It is well known that the discrete logistic model with a standard discretization, e.g., forward Euler scheme, has solutions with various periods, as well as chaotic solutions [3,14]. The nonstandard scheme was introduced by [3] to alleviate these problems and lead to difference equations that preserve the qualitative behavior of their continuous counterparts. It must be stressed that when a threshold policy with hysteresis effect is considered, we are implicitly considering errors and delays in the implementation of the policy, i.e., errors and delays in the measurement of the species density. The hysteresis loop around the threshold level takes such occurrences into account explicitly and stabilizes the discrete logistic model in low amplitude bounded oscillations within the region delimited by the hysteresis itself. Finally, from the management standpoint, it is important to stress that the proposed strategy does not interfere directly in the harvesting intensities. Acknowledgement This research was partially financed by Project E-26/100.455/2007 (CNE) of FAPERJ, and Project No. 154631/2006-0 of CNPq, and also by the agency CAPES and FAPESP. We thank the reviewer for his constructive and detailed criticism and helpful suggestions that, in our opinion, have improved the manuscript.
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
7883
Appendix A. Another possibility for nonstandard discretization The one-species discrete models are obtained according to nonstandard discretization method [12], by making the folx x lowing replacements: (i) dx ! kþ1/ k where / is such that /(h) = h + O(h2) and h = Dt, e.g., /(h) = eh 1, (ii) dt x = 2x x ? 2xk xk+1, (iii) x2 ? xk+1xk. The resulting one-species discrete models are as follows: Continuous model x_ ¼ rx 1 Kx cx x_ ¼ rx 1 Kx xþd x_ ¼ rx 1 Kx
Discrete model ð1þ2r/Þxk xkþ1 ¼ 1þr/þ r K /xk
xkþ1 ¼
1þ2r/x2c/ xk þd k
1þr/þKr /xk x c/þd k
c 2 x2 x2 þd2
xkþ1 ¼
ð1þ2r/Þxk c2 /x 1þr/þKr /xk þ 2 k2 x þd k
Appendix B. Proof of the Theorems 2 and 1 Trajectories initiating in the intervals I1 = [0, xL], I3 = [xR, 1) (bold lines, Fig. B.1) will be shown to converge to the interval I2 = (xL, xR) which contains the points xLL = xth r and xRR = xth + r which determine the thresholds, in the sense that convergence occurs to a neighborhood of the interval (xLL, xRR). I2 is shown to be an invariant interval. The proof of Theorem 1 is similar to the proof of Theorem 2. For the proof of Theorem 1 is necessary to consider r = 0. The difference is that there are two different switching points (thresholds), one switching point for trajectories, which initiate in the interval [0, xLL] and switch when crossing the point xRR = xth + r and another switching point for trajectories, which initiate in the interval [xRR, 1) and switch when crossing the point xLL = xth r. All trajectories remain inside the interval I2. B.1. Equilibrium points under TPH and TP The system considered to prove the theorem is a logistic discrete model, which is subjected to the TPH,
xkþ1 ¼
ð1 þ r/Þxk uTPH ; 1 þ Kr /xk
where uTPH = exkwhys(s) and whys(s) is defined as in Eq. (3). The function
f ðxk ; /Þ xk ¼
ð1 þ r/Þxk xk 1 þ Kr /xk
is positive for values of xk < xVNC. The equilibrium point of the system without harvesting (without control) is
xVNC ¼ K and the equilibrium point of the system with harvesting (proportional control) is
xVC ¼ K
r/ e : r/ð1 þ eÞ
If the effort parameter satisfies the restriction e < r/ then the equilibrium point xVC is positive. If e does not satisfy this restriction, then the system without control goes to extinction. In this paper, we consider values of e slightly greater than r/, i.e.,
e P r/:
ðB:1Þ
We choose the threshold level such that the equilibria xVNC and xVC are virtual, see Fig. B.1, i.e.,
xVC < xLL < xth < xRR < xVNC : B.2. Determining the bounds of the globally attractive invariant interval We first calculate an interval that contains xLL = xth r and xRR = xth + r (the threshold values) by reasoning as follows. Suppose that the state xk ¼ xþ RR (i.e., infinitesimally to the left of xRR), then it is in the region without harvesting and the next state, which will denote xR, can be written as
xR ¼
ð1 þ r/ÞxRR : 1 þ Kr /xRR
For the proof of Theorem 1 is necessary to consider r = 0, then xRR = xth.
ðB:2Þ
7884
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
(a)
(b) Fig. 8. The comparison between a nonstandard scheme and the standard scheme discretization for the logistic and Lotka–Volterra discrete models. (a) Simulation of the discrete logistic model for a nonstandard and standard scheme discretization. Parameter values for both simulations are the same. r = 3, K = 1, initial condition 0.1. (b) Simulation of the Lotka–Volterra discrete model for a nonstandard and standard scheme discretization. Parameter values for both simulations are the same. a = 1, b = 1, r1 = 1, r2 = 1, initial condition (2.1, 2.1).
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
7885
Similarly, if xþ LL (i.e., infinitesimally to the right of xLL) and therefore in the region with harvesting, the next state, denoted as xL, can be calculated as
xL ¼
ð1 þ r/ÞxLL exLL : 1 þ Kr /xLL
ðB:3Þ
For the proof of Theorem 1 is necessary to consider r = 0, then xLL = xth. We will show that xL and xR are the left and right endpoints of an interval, denoted I2, that will be proved to be globally attractive and invariant under the TPH and TP. Eqs. (B.2) and (B.3) define the endpoints of the interval I2. B.3. Proof of the global attractivity and invariance of the interval I2 A suitable reaching phase Lyapunov function, adapted from [15], is
Vðxk Þ ¼ xk 2xth þ
x2th ðxk xth Þ2 ¼ xk xk
ðB:4Þ
for which
x2 DVðxk Þ ¼ Vðxkþ1 Þ Vðxk Þ ¼ ðxkþ1 xk Þ 1 th : xkþ1 xk
ðB:5Þ
This function is used to prove: B.3.1. Global attractivity of interval I2 We now use V, DV to carry out an interval analysis of the system (10). Interval I1: Assume xk 2 I1. If xk+1 2 I1, then xk+1 6 xth, and from the dynamics in this interval, we have
xkþ1 ¼
ð1 þ r/Þxk : 1 þ Kr /xk
ðB:6Þ
Since xk < xVNC = K, "k, the factor ð1 þ r/Þ= 1 þ Kr /xk is greater than 1, thus from (B.6), xk+1 > xk and xk+1 6 xth. Therefore, from (B.5), D V < 0 and all trajectories with initial condition in the interval I1 converges to the interval I2. Interval I3: Assume xk 2 I3. If xk+1 2 I3, then xk+1 P xth, and from the dynamics in this interval, we have
xkþ1
1 þ r/ e 1 þ Kr /xk xk ð1 þ r/Þxk ¼ exk ¼ : 1 þ Kr /xk 1 þ Kr /xk
K r/e . The latter holds by (11) since xk is a nonnegative The factor ð1 þ r/Þ= 1 þ Kr /xk e is smaller than 1 for xk > r/ 1þe variable. Thus xk+1 < xk and xk+1 P xth, "xk 2 I3. Therefore, DV < 0 and all trajectories with initial condition in the interval I3 converge to the interval I2. Therefore, trajectories initiating in the intervals I1 and I3 converge to the interval I2, proving that it is globally attractive. B.3.2. Invariance of interval I2 Inside the interval I2, the plan of proof is to show that points in I2 to the left of xRR cannot jump, under the dynamics (10), to the right of xR and, similarly, points to the right of xLL cannot jump to the left of xL. This completes the proof of invariance of interval I2.
Fig. B.1. Graphical representation of the phase interval. xVC and xVNC are the equilibrium point of the system with harvesting and without harvesting, respectively. xL and xR are the left bounded and the right bounded of the interval I2, respectively. xLL = xth r and xRR = xth + r are the threshold levels.
7886
M.E. Mendoza Meza, A. Bhaya / Applied Mathematics and Computation 217 (2011) 7874–7886
Fig. B.2. Graphical representation inside the interval I2. xl is the initial condition in the region without harvesting. xr is the initial condition in the region with harvesting.
1. Consider an initial condition xk in the region without harvesting (Fig. B.2) that belongs to the interval (xL, xRR) and suppose that the next state is on the right of xR, i.e., xk+1 2 I3, and xk+1 > xR,
ð1 þ r/Þxk ð1 þ r/ÞxRR > : 1 þ Kr /xk 1 þ Kr /xRR
ðB:7Þ
Inequality (B.7) is satisfied if xk > xRR, which contradicts xk 2 I2. In other words, if xk 2 I2, then xk+1 is smaller than xR, i.e., the trajectory remains inside the interval I2. For the proof of Theorem 1 is necessary to consider r = 0, then xRR = xth. 2. Consider an initial condition xk in the region with harvesting (Fig. B.2) that belongs to the interval (xLL, xR) and suppose that next state is on the left of xL, i.e., xk+1 2 I1. Let a = 1 + r/, b = r//K for notational simplicity. Then xk+1 < xL can be written as
axk axLL exk < exLL ; 1 þ bxk 1 þ bxLL which implies
ðxr xLL Þða eð1 þ bxr Þð1 þ bxLL ÞÞ < 0:
ðB:8Þ
We will prove that the second term of (B.8) is positive. In fact,
a eð1 þ bxr Þð1 þ bxLL Þ > 0 ()
e<
1 þ r/ : x 1 þ r/ x 1 þ r/ K r K LL
By assumption 2 and for values of xk and xLL in the interval [0, xVNC], we have
e<
1 1 þ r/ < 1 þ r/: < 1 þ r/ x 1 þ r/ x 1 þ r/ K r K LL
Thus, Eq. (B.8) satisfied for xr < xLL, which contradicts xk 2 (xLL, xR). In other words, if xk 2 (xLL, xR), then xk+1 is greater than xL, i.e., the trajectory remains inside the interval I2. For the proof of Theorem 1 is necessary to consider r = 0, then xLL = xth. The proof of the analogs of Theorems 1 and 2 for the logistic discrete model with endogenous consumption of Holling Type II and Holling Type III, models in Table 1, can be carried out similarly. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]
X. Liu, D. Xiao, Complex dynamics behaviors of discrete-time predator–prey system, Chaos Solutions Fractals 32 (2007) 80–94. P. Liu, S.N. Elaydi, Discrete competitive and cooperative models of Lotka–Volterra type, J. Comput. Anal. Appl. 3 (1) (2001) 53–73. R.E. Mickens, Nonstandard Finite Difference Models of Differential Equations, World Scientific, Singapore, 1994. M. Basson, M.J. Fogarty, Harvesting in discrete-time predator–prey systems, Math. Biosci. 141 (1) (1997) 41–74. H. Kokko, J. Lindström, Seasonal density dependence, timing of mortality, and sustainable harvesting, Ecol. Model. 110 (3) (1998) 293–304. S.Y. Tang, L.S. Chen, The effect of seasonal harvesting on stage-structured population models, J. Math. Biol. 48 (4) (2004) 357–374. Magno Enrique Mendoza Meza, Amit Bhaya, Eugenius Kaszkurewicz, Michel Iskin Silveira Costa, On-off policy and hysteresis on-off policy control of the herbivore-vegetation dynamics in a semi-arid grazing system, Ecol. Eng. 28 (2) (2006) 114–123. I. Noy-Meir, Stability of grazing systems: an application of predator–prey graphs, J. Ecol. 63 (2) (1975) 459–481. C.W. Clark, Mathematical Bioeconomics: The Optimal Management of Renewable Resources, John Wiley & Sons, USA, 1976. C.W. Clark, Bioeconomics Modelling and Fisheries Management, John Wiley & Sons, USA, 1985. M. Kot, Elements of Mathematical Ecology, Cambridge University Press, Cambridge, 2001. R.E. Mickens, A nonstandard finite-difference scheme for the Lotka–Volterra system, Appl. Numer. Math. 45 (2003) 309–314. M.I.S. Costa, E. Kaszkurewicz, A. Bhaya, L. Hsu, Achieving global convergence to an equilibrium population in predator–prey systems by the use of discontinuous harvesting policy, Ecol. Model. 128 (2000) 89–99. J.D. Murray, Mathematical Biology I. An Introduction, Springer-Verlag, New York, 2002. B.-S. Goh, Management and analysis of biological populations, Developments in Agricultural and Managed-Forest Ecology, vol. 8, Elsevier, 1980.