Table Of Content

Nonstandard finite difference schemes for a general predator-prey system Quang A Danga, Manh Tuan Hoangb a Center for Informatics and Computing, Vietnam Academy of Science and Technology (VAST), 18 Hoang Quoc Viet, Cau Giay, Hanoi, Vietnam b Institute of Information Technology, Vietnam Academy of Science and Technology (VAST), 7 1 18 Hoang Quoc Viet, Cau Giay, Hanoi, Vietnam 0 2 n a J Abstract 0 2 In this paper we transform a continuous-time predator-prey system with general functional response and recruitment for both species into a discrete-time model by nonstandard finite difference scheme ] A (NSFD). The NSFD model shows complete dynamic consistency with its continuous counterpart for N any step size. Especially, the global stability of a non-hyperbolic equilibrium point in a particular . case of parameters is proved by the Lyapunov stability theorem. The performed numerical simulations h confirmed the validity of the obtained theoretical results. t a m Keywords: Predator-Prey system; Nonstandard finite-difference scheme; Dynamically consistent; [ Lyapunov stability theorem; Global stability. 1 v 3 1. Introduction 6 6 Predator-prey systems are among the most interested topics in mathematical biology and 5 0 ecology. Their dynamics continues to attract attention from both applied mathematicians and . 1 ecologists because of its universal existence and importance [1, 5, 6, 18]. The majority of the 0 results in this topic is mainly concentrated on the study of the qualitative aspects of the contin- 7 1 uous systems described by systems of differential equations. It is worth to mention some recent : typical works such as [16, 17, 20, 29, 33, 34]...In that time the conversion of the continuous sys- v i tems to discrete systems preserving the properties of the original continuous systems is of great X importance. This problem has attracted attention from many researchers and a lot of results r a are reached for some dynamical systems. Nevertheless, some continuous models are studied fully from theoretical point of view but their discrete counterparts still are not investigated. Biological systems including predator-prey systems often are described by ordinary or partial differential equations. There are many ways for converting continuous models to discrete counterparts. The most popular way for this purpose is to use standard difference methods such as Euler, Runge-Kutta schemes. However, in many nonlinear problems the standard difference schemes reveal a serious drawback which is called ”numerical instabilities” [23, 24, 26]. Under this concept we have in mind the phenomena when the discrete models, for example, the difference schemes, do not preserve properties of the corresponding differential equations. In Email addresses: [email protected](Quang A Dang), [email protected](Manh Tuan Hoang) Preprint submitted to Elsevier January 23, 2017 [23, 24, 25, 26] Mickens showed many examples and analysed the numerical instabilities when using standard difference schemes. In general, standard difference schemes preserve the properties of the differential equations only in the case if the discretization parameter h is sufficiently small. Therefore, when studying dynamical models in large time intervals the selection of small time steps will requires very large computational effort, so these discrete models are inefficient. Besides, for some special dynamical problems standard difference schemes cannot preserve the properties of the problems for any step sizes. In order to overcome the numerical instabilities phenomena in 1989 Mickens [22] introduced theconcept Nonstandard Finite Difference (NSDF)schemes andafterthat hasdeveloped NSDF methodsin many works, such as[23, 24, 25, 26]. According to Mickens, NSDFschemes arethose constructed following a set of five basic rules. The NSDF schemes preserve main properties of the differential counterparts, such as positivity, monotonicity, periodicity, stability and some other invariants including energy and geometrical shapes. It should be emphasized that NSFD schemes can preserve all properties of the continuous models for any discretization parameters. The discrete models with these properties are called dynamically consistent [2, 3, 8, 10, 11, 12, 13, 14, 26, 27, 30]. Up to date NSFD schemes became a power and efficient tool for simulating dynamical systems, especially inconverting continuous modelstodynamically consistent discrete counterparts [2, 3, 7, 9, 10, 11, 12, 13, 25, 28, 31, 32, 35]. The majority of these models are met in physics, mechanics, chemistry, biology, epidemiology, finance, ...with complicated dynamical behaviour. For predator-prey systems, some NSFD schemes, which are dynamically consistent with them, are constructed, such as [4, 9, 12, 14]. Below we mention some models. In 2006 Dimitrov and Kojouharov [12] used NSFD schemes for constructing discrete model for the general Rosenzweig-MacArthur predator-prey model with a logistic intrinsic growth of the prey population. The model has the form dx dy = bx(x−1)−ag(x)xy, = g(x)xy −dy, dt dt where x and y represent the prey and predator population sizes, respectively, b > 0 represents theintrinsicgrowthrateoftheprey, a > 0standsforthecapturingrateandd > 0isthepredator death rate. Under some assumptions for the function g(x) the obtained discrete model preserves thepositivity ofsolutions andlocallyasymptotical stability of theset of equilibrium points. The finite difference method is called positive and elementary stable nonstandard (PESN) method. Next, in 2008 Dimitrov and Kojouharov [14] constructed a PESN scheme for predator-prey models with general functional response of the form dx dy = p(x)−af(x,y)y, = f(x,y)y−µ(y), dt dt where functions p(x) and µ(y) describe the intrinsic growth rate of the prey and the mortality rate of the predator, respectively, the function f(x,y) is called ”functional response” and represents theper capita predator”feeding rate”per unit time. Recently, in2016, BairagiandBiswas [4] used PESN scheme for converting a predator-prey model with Beddington-DeAngelis Func- tional Response to a dynamically consistent discrete model. The model under consideration is of the form dx αxy dy Exy = x(x−1)− , = −Dy, dt 1+βx+µy dt 1+βx+µy 2 where α,β,µ,E,D are positive constants. This model is a particular case of the model considered by Dimitrov and Kojouharov in [14]. An another interesting predator-prey model, which should be mentioned is the harvesting Leslie-Gower predator-prey model. In the recent paper [9] the authors also constructed NSFD scheme preserving the positivity of solutions and asymptotical stability for this model. Itshouldbeemphasized thatalltheequilibriumpointsinthementionedabovepredator-prey models are hyperbolic, therefore, for establishing the stability, it suffices to consider eigenvalues of the Jacobian matrix of the linearized system around equilibrium points. Aiming at the conversion of continuous predator-prey models studied fully from the theoretical point of view to dynamically consistent discrete models, in the present paper we consider a mathematical model for a predator-prey system with general functional response and recruitment, which also includes capture on both species [21]. The results of qualitative aspects of the model are given in [21]. One important difference of this model in comparison with the above models is that the model [21] has one non-hyperbolic equilibrium point. The paper is organized as follows. In Section 2 we recall from [21] the predator-prey system under consideration with the theoretical results of the existence of equilibrium points and their stability properties. In Section 3 we propose NSFD schemes and study their positivity and existence of equilibrium points. Next, Section 4 is devoted to the stability analysis of the equilibrium points. Some numerical simulations for demonstrating the validity of the theoretical results obtained in the previous section are given in Section 5. Finally, Section 6 is Conclusion. 2. Mathematical model We consider a mathematical model for a predator-prey system with general functional response and recruitment for both species [21]. The model is described by the system of nonlinear differential equations x˙(t) = x(t)f(x(t),y(t)) = x(t) r(x(t))−y(t)φ(x(t))−m , 1 (cid:2) (cid:3) (1) y˙(t) = y(t)g(x(t),y(t)) = y(t) s(y(t))+cx(t)φ(x(t))−m , 2 (cid:2) (cid:3) where: • x(t) and y(t) are prey population and predator population, respectively; • r(x) and s(y) are the per capita recruitment rates of prey and predators, respectively; • xyφ(x) is the predator response, and xφ(x) is the number of prey consumed per predator in a unit of time; • c is a constant named conversion efficiency of prey into predators, generally 0 < c < 1 , and cxyφ(x) is the predator numerical response; • m > 0 and m > are the total mortality rates of prey and predators, respectively. 1 2 From the biological significance we have ′ ′ ∀x ≥ 0, r(x) > 0, r (x) < 0, [xr(x)] ≥ 0, and lim r(x) = 0, x→∞ ′ ′ ∀y ≥ 0, s(y) > 0, s(y) < 0, [ys(y)] ≥ 0, and lim s(y) = 0, (2) y→∞ ′ ′ ∀x ≥ 0, φ(x) > 0, φ(x) ≤ 0, and [xφ(x)] ≥ 0. 3 It is easy to check that the region Ω = R2 is a positive invariant set for the system (1). + Other qualitative properties of the model (1) including the existence of equilibrium points, their stability, are studied in detail in [21]. In order to easily track the paper, below we mention these properties. Theorem 1. (Existence of equilibrium points) [21, Proposition 1] System (1) has four distinct kinds of possible equilibrium points in Ω: (i) A trivial equilibrium point P∗ = (x∗,y∗) = (0,0), for all the values of the parameter. 0 0 0 (ii) An equilibrium point of the form P∗ = (x∗,y∗) = (K,0), with r(K) = m , if and only if 1 1 1 1 m < r(0). 1 (iii) An equilibrium point of the form P∗ = (x∗,y∗) = (0,M), with s(M) = m , if and only if 2 2 2 2 m < s(0). 2 (iv) An equilibrium point of the form P∗ = (x∗,y∗) = (x∗,y∗), where x∗ satisfies the equation 3 3 3 r(x∗)−m ∗ ∗ 1 cx φ(x )+s −m = 0, (cid:16) φ(x∗) (cid:17) 2 and y∗ is given, as a function of x∗, by r(x∗)−m ∗ 1 y = , φ(x∗) if and only if (m ,m ) verifies m < r(0) − Mφ(0) and m < s(0) or m < r(0) and 1 2 1 2 1 s(0) < m < s(0)+cKφ(K). 2 Theorem 2. (Stability analysis) (i) If m > r(0) and m > s(0), then the extinction equilibrium point P∗ = (0,0) is locally 1 2 0 asymptotically stable, and unstable otherwise. (ii) If m ≥ r(0) and m ≥ s(0) then the extinction equilibrium point P∗ is globally asymptot- 1 2 0 ically stable. (iii) If m < r(0) and m > s(0)+cKφ(K), then the equilibrium point of the form P∗ = (K,0) 1 2 1 is locally asymptotically stable, and unstable otherwise.The equilibrium point (K,0) shall be called the equilibrium point of extinction of the predator species. (iv) If m > r(0)−Mφ(0) and m < s(0), then the equilibrium point of the form P∗ = (0,M) 1 2 2 is locally asymptotically stable, and unstable otherwise. The equilibrium point (0,M) shall be called the equilibrium point of extinction of the prey species. (v) If an equilibrium point of the form P∗ = (x∗,y∗) belongs to Ω, then it is locally asymptot- 3 ically stable. This equilibrium point (x∗,y∗) shall be named the ecological stability equilibrium. 4 In general, the property of stability of the set of equilibria of differential equations plays an essential role in the study of asymptotical behavior of the solutions of differential equations. The construction of difference schemes, which preserve the stability of the equilibrium points, is important in numerical simulation of differential equations. The difference schemes with this stability property is called elementary stable schemes [2, 10, 11, 35]. There are many works concerning the elementary stable schemes. The typical results are for general dynamical systems [10, 11, 19] and for other specific systems [14, 30, 31, 32, 35] .... One popular approach to the elementary stability is the investigation of Jacobian matrices of the discrete models at the equilibria, namely, determination of conditions ensuring that all eigenvalues of Jacobian matriceshavemodulilessorequal 1. Thisisthenecessary andsufficient conditionforhyperbolic equilibrium points to be locally stable [1, 18]. The mentioned above approach has the following weaknesses and limitations: 1. It is applicable when all the equilibrium points are hyperbolic. To our best knowledge, at present no results on NSFDschemes preserving the stability of non-hyperbolic equilibrium points are available. 2. The consideration of Jacobian only guarantees the local stability meanwhile many models have the global stability. Return to the model (1). A difficulty is that when m = r(0) or m = s(0) the equilibrium point 1 2 P∗ = (0,0) becomes non-hyperbolic. Therefore, it is impossible to study its stability via the 0 eigenvalues of Jacobian J(P∗). Consequently, we cannot directly apply the results concerning 0 the elementary stability [10, 11, 19] of the system. For the continuous model, the Lyapunov stability theorem should be used for investigating the global stability [21, Theorem 1]. The models considered in [4, 9, 12, 14] have equilibrium points, which are all hyperbolic. Meanwhile the system (1) has a non-hyperbolic equilibrium point in a particular case of the parameters. Moreover, it is globally asymptotical stable. Therefore, for the corresponding discrete system we shall use the Lyapunov stability theorem for proving the stability of this equilibrium point. This is the important contribution of our paper. Besides, we propose a more general model with many iterative parameters in the discretization of the right-hand sides. The combination of appropriate selection of denominator with these parameters will give sufficient conditions for dynamical consistency of the discrete model with the continuous system. 3. Construction of NSFD scheme In this section we construct NSFD scheme for the system (1) preserving all dynamic properties of the original continuous model for any discretization parameter or step size h > 0. Recall that, according to Mickens, a finite difference scheme is called nonstandard if at least one of the following conditions is satisfied [23, 24, 25, 26]: • A nonlocal approximation is used. • The discretization of the derivative is not traditional and uses a function 0 < ϕ(h) = h+O(h2). 5 We propose a general NSFD scheme for the model (1) in the form x −x k+1 k = α x r(x )+α x r(x )−α x y φ(x )−α x y φ(x )−α m x −α m x , 1 k k 2 k+1 k 3 k k k 4 k+1 k k 5 1 k 6 1 k+1 ϕ(h) y −y k+1 k = β y s(y )+β y s(y )+cβ x y φ(x )+cβ x y φ(x )−β m y −β m y , 1 k k 2 k+1 k 3 k k k 4 k k+1 k 5 2 k 6 2 k+1 ϕ(h) α +α = β +β = 1, j = 1,3,5; ϕ = h+O(h2), h → 0. j j+1 j j+1 (3) Concerning theset ofequilibrium pointsandpositivityofthescheme (3)thereholdthefollowing assertions. Proposition 1. The region Ω = (x,y) ∈ R2 x ≥ 0, y ≥ 0 is a positive invariant set for (cid:8) (cid:12) (cid:9) the scheme (3) if (cid:12) α ≥ 0, α ≤ 0, α ≤ 0, α ≥ 0, α ≤ 0, α ≥ 0, 1 2 3 4 5 6 (4) β ≥ 0, β ≤ 0, β ≥ 0, β ≤ 0, β ≤ 0, β ≥ 0. 1 2 3 4 5 6 Proof It is easy to reduce the scheme (3) to the explicit form x +ϕα x r(x )−ϕα x y φ(x )−ϕα m x k 1 k k 3 k k k 5 1 k x = F(x ,y ) := , k+1 k k 1−ϕα r(x )+ϕα y φ(x )+ϕα m 2 k 4 k k 6 1 (5) y +ϕβ y s(y )+ϕβ cx y φ(x )−ϕβ m y k 1 k k 3 k k k 5 2 k y = G(x ,y ) := . k+1 k k 1−ϕβ s(y )−ϕβ cx φ(x )+ϕβ m 2 k 4 k k 6 2 Since the parameters α and β (j = 1,6) satisfy (4) from the above formulas we obtain the fact j j to be proved. Proposition 2. The scheme (3) preserves the set of equilibrium points of system (3). Proof The equilibrium (x∗,y∗) of the scheme (3) is the solution of the equation ∗ ∗ ∗ ∗ ∗ ∗ x = F(x ,y ), y = G(x ,y ), where F(x,y) and G(x,y) are defined by (5). Since the parameters α ,β satisfy (3), if the i i scheme (3) is defined then it is equivalent to ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ x r(x )−y φ(x )−m = 0, y s(y )+cx φ(x )−m = 0. 1 2 (cid:2) (cid:3) (cid:2) (cid:3) It is the system for determining equilibrium points of the model (1). Thus, the proposition is proved. 4. Stability analysis In this section we give sufficient conditions for the scheme (3) to preserve the stability properties of equilibrium points of the model (1). For the purpose of easy tracking we recall the following results of the necessary and sufficient conditions for an equilibrium point to be locally asymptotically stable [1, Theorem 2.10], [18]. 6 Lemma 1. Assume the functions f(x,y) and g(x,y) have continuous first-order partial deriva- tives in x and y on some open set in R2 that contains the point (x∗,y∗). Then the equilibrium point (x∗,y∗) of the nonlinear system x = f(x ,y ), y = g(x ,y ), k+1 k k k+1 k k is locally asymptotically stable if the eigenvalues of the Jacobian matrix J evaluated at the equilibrium satisfy |λ | < 1 iff i (i)det(J) < 1, (ii)1−Tr(J)+det(J) > 0, (iii)1+Tr(J)+det(J) > 0. The equilibrium is unstable if some |λ | > 1, that is, if any one of three inequalities is satisfied, i (i)det(J) > 1, (ii)1−Tr(J)+det(J) < 0, (iii)1+Tr(J)+det(J) < 0. 4.1. The extinction equilibrium point P∗ = (0,0) 0 Proposition 3. Forthe casem > r(0) and m > s(0) considerthe differencescheme (3) under 1 2 the assumptions of Proposition 1. Then the extinction equilibrium point P∗ = (x∗,y∗) = (0,0) 0 0 0 is locally asymptotically stable if m > r(0) and m > s(0), and unstable otherwise. 1 2 Proof Computing the Jacobian matrix of system (3), evaluated at the extinction equilibrium point P∗ = (x∗,y∗) = (0,0) , one obtains 0 0 0 1+ϕα r(0)−ϕα m 1 5 1  0  J(P∗) = 1−ϕα2r(0)+ϕm1α6 . 0  1+ϕβ1s(0)−ϕβ5m2  0   1−ϕβ2s(0)+ϕβ6m2 In this case, the eigenvalues are 1+ϕα r(0)−ϕα m 1+ϕβ s(0)−ϕβ m 1 5 1 1 5 2 λ = , λ = . 1 2 1−ϕα r(0)+ϕm α 1−ϕβ s(0)+ϕβ m 2 1 6 2 6 2 Since α and β satisfy Proposition 1 we have λ > 0, λ > 0. On the other hand j j 1 2 r(0)−m s(0)−m 1 2 |λ |−1 = ϕ , |λ |−1 = ϕ , 1 2 1−ϕα r(0)+ϕm α 1−ϕβ s(0)+ϕβ m 2 1 6 2 6 2 Therefore, |λ |,|λ | < 1 iff m > r(0) and m > s(0). By Lemma 1 Proposition is proved. 1 2 1 2 Next, consider scheme (3) under the assumptions of Proposition 1 in the case m ≥ r(0) and 1 m ≥ s(0). Notice that when m = r(0) or m = s(0) the equilibrium point P∗ = (0,0) becomes 2 1 2 0 non-hyperbolic one. Therefore, it is impossible to use Lemma 1 for proving the locally asymptotical stability. Moreover, the equilibrium point P∗, as shown in [21, Theorem 1], is globally 0 asymptotically stable. As for the continuous system we shall also use Lyapunov’s stability theorem to prove that the point P∗ = (0,0) is a globally asymptotically stable equilibrium point 0 of the scheme (3). For this purpose we consider a family of functions V(x ,y ) := αx y +βx2 +γx +δy , (x ,y ) ∈ R2, (6) k k k k k k k k k + 7 where α,β,γ,δ > 0 are parameters, which are selected so that the function V(x ,y ) satisfies k k the conditions of the Lyapunov’s stability theorem [15, Theorem 4.20]. Obviously, the function V(x ,y ) defined by (6) is continuous on R2 and V(x ,y ) → ∞ as k k + k k ||(x ,y )|| → ∞. Moreover, V(P∗) = 0andV(x ,y ) > 0forany(x ,y ) ∈ R2,(x ,y ) 6= (0,0). y k 0 k k k k + k k Therefore, in order to show that V(x ,y ) satisfies [15, Theorem 4.20] it suffices to determine k k conditions for ∆V(x ,y ) = V(x ,y )−V(x ,y ) < 0 ∀(x ,y ) ∈ R2\{(0,0)}. k k k+1 k+1 k k k k + We have ∆V(x ,y ) = α(x y −x y )+β(x2 −x2)+γ(x −x )+δ(y −y ). (7) k k k+1 k+1 k k k+1 k k+1 k k+1 k Now we rewrite (5) in the form x [r(x )−m ]−x y φ(x ) k k 1 k k k x = x +ϕ , k+1 k 1−ϕα r(x )+ϕα y φ(x )+ϕα m 2 k 4 k k 6 1 (8) y [s(y )−m ]+cx y φ(x ) k k 2 k k k y = y +ϕ . k+1 k 1−ϕβ s(y )−ϕβ cx φ(x )+ϕβ m 2 k 4 k k 6 2 According to the estimates in Theorem 1 in [21] we have: if m ≥ r(0) and m ≥ s(0), then 1 2 r(x)−m ≤ 0 for all x ≥ 0 and s(y)−m ≤ 0 for all y ≥ 0. Therefore, from (8) we see that if 1 2 m ≥ r(0), then x ≤ x for all k ≥ 0. Consequently, from (7) it implies that if y −y < 0 1 k+1 k k+1 k then ∆V(x ,y ) < 0. Hence, it is sufficient to consider only the case y ≥ y . k k k+1 k Using the mean value theorem for two variables function u(x ,y ) = x y we have k k k k x y −x y = ξ (x −x )+ξ (y −y ) ≤ y (x −x )+x (y −y ), (9) k+1 k+1 k k yk k+1 k xk k+1 k k k+1 k k k+1 k where ξ ∈ (x ,x ) and ξ ∈ (y ,y ). On the other hand, since x ≤ x we have xk k+1 k yk k k+1 k+1 k x2 −x2 ≤ x x −x2 = x (x −x ). (10) k+1 k k k+1 k k k+1 k From (9), (10), (7) and (8) we obtain the estimate ∆V(x ,y ) ≤ (αy +βx +γ)(x −x )+(αx +δ)(y −y ) k k k k k+1 k k k+1 k x [r(x )−m ]−x y φ(x ) k k 1 k k k = (αy +βx +γ)ϕ k k 1−ϕα r(x )+ϕα y φ(x )+ϕα m 2 k 4 k k 6 1 (11) y [s(y )−m ]+cx y φ(x ) k k 2 k k k +(αx +δ)ϕ k 1−ϕβ s(y )−ϕβ cx φ(x )+ϕβ m 2 k 4 k k 6 2 ≤ ϕx y φ(x )Q(x ,y ), k k k k k where Q(x ,y ) is defined by k k c(αx +δ) αy +βx +γ k k k Q(x ,y ) := − . k k 1−ϕβ s(y )−ϕβ cx φ(x )+ϕβ m 1−ϕα r(x )+ϕα y φ(x )+ϕα m 2 k 4 k k 6 2 2 k 4 k k 6 1 From here it follows that if Q(x ,y ) < 0 for all (x ,y ) ∈ R2\{(0,0)} then ∆V(x ,y ) < 0 for k k k k + k k all (x ,y ) ∈ R2\{(0,0)}. On the other hand we have k k + Q(x ,y )[1−ϕα r(x )+ϕα y φ(x )+ϕα m ][1−ϕβ s(y )−ϕβ cx φ(x )+ϕβ m ] k k 2 k 4 k k 6 1 2 k 4 k k 6 2 = τ x +τ ϕx +τ ϕcαx y φ(x )+τ +ϕτ +τ ϕy −Q (x ,y ), 1 k 2 k 3 k k k 4 5 6 k 1 k k 8 where τ := cα−β, τ := −α cαr(x )+α cαm −ββ m , τ := α +β , 1 2 2 k 6 1 6 2 3 4 4 τ := cδ −γ, τ := −α cδr(x )+cδα m −β γm , τ := α cδφ(x )−αβ m , 4 5 2 k 6 1 6 2 6 4 k 6 2 Q (x ,y ) := βx [−ϕβ s(y )−ϕβ cx φ(x )]+αy [1−ϕβ s(y )]+γ[−ϕβ s(y )−ϕβ cx φ(x )]. 1 k k k 2 k 4 k k k 2 k 2 k 4 k k (12) Notice that Q (x ,y ) > 0 for all (x ,y ) ∈ R2\{(0,0)}. Therefore, Q(x ,y ) < 0 for all 1 k k k k + k k (x ,y ) ∈ R2\{(0,0)} if τ < 0 (i = 1,6). Taking into account (2) we have 0 < r(x ) ≤ r(0) k k + i k and 0 < φ(x ) ≤ φ(0) for all x ≥ 0. Therefore from (12) it follows k k β 1. τ < 0 if c < , 1 α 2. τ = −α cαr(x )+α cαm −ββ m ≤ −α cαr(0)+α cαm −ββ m < 0 2 2 k 6 1 6 2 2 6 1 6 2 −α cr(0)+α cm β 2 6 1 if < , β m α 6 2 3. τ < 0 if α +β < 0, 3 4 4 γ 4. τ < 0 if c < , 4 δ 5. τ = −α cδr(x )+cδα m −β γm ≤ −α cδr(0)+cδα m −β γm < 0 5 2 k 6 1 6 2 2 6 1 6 2 −α cr(0)+cα m γ 2 6 1 if < , β m δ 6 2 6. τ = α cδφ(x )−αβ m ≤ α cδφ(0)−αβ m < 0 6 4 k 6 2 4 6 2 cα φ(0) α 4 if < . β m δ 6 2 In summary the function V(x ,y ) defined by (6) satisfies ∆V(x ,y ) < 0 for all (x ,y ) ∈ k k k k k k R2\{(0,0)} if + −α cr(0)+α cm β −α cr(0)+cα m γ 2 6 1 2 6 1 max c, < , max c, < , n β m o α n β m o δ 6 2 6 2 (13) cα φ(0) α 4 < , α +β < 0. 4 4 β m δ 6 2 Once, the scheme (3) is fixed the selection of the parameters α,β,γ,δ > 0 satisfying the above relations is completely possible. Thus, we obtain the following theorem of the global stability of the equilibrium point P∗ = (0,0). 0 Theorem 3. For the case m ≥ r(0) and m ≥ s(0) consider the difference scheme (3) under 1 2 the assumptions of Proposition 1. If additionally assume that α +β < 0, (14) 4 4 then the extinction equilibrium point P∗ = (0,0) is globally asymptotically stable. 0 9 4.2. The equilibrium point P∗ = (K,0) (the equilibrium point of extinction of the predator 1 species) Proposition 4. For the case m < r(0) and m > s(0) + cKφ(K) consider the difference 1 2 scheme (3) under the assumptions of Proposition 1. If additionally assume that ′ T := 2α m −2α r(K)+Kr (K) > 0, 1 6 1 2 (15) T := s(0)−m +cKφ(K)−2β s(0)−2β cKφ(K)+2β m > 0, 2 2 2 4 6 2 then the equilibrium point of the form P∗ = (K,0) is locally asymptotically stable if m < r(0) 1 1 and m > s(0)+cKφ(K), and unstable otherwise. 2 Proof Recall that m < r(0) is necessary and sufficient condition for the existence of the 1 equilibrium point P∗ (see Theorem 1). Computing the Jacobian matrix of system (3), evaluated 1 at the extinction equilibrium point P∗ = (K,0), one obtains 1 ϕKr′(K) −ϕKφ(K) 1+  J(P∗) = 1−ϕα2r(K)+ϕα6m1 1−ϕα2r(K)+ϕα6m1 . 1  s(0)−m2 +cKφ(K)   0 1+ϕ   1−ϕβ2s(0)−ϕβ4cKφ(K)+ϕβ6m2 In this case, the eigenvalues are ϕKr′(K) s(0)−m +cKφ(K) 2 λ = 1+ , λ = 1+ϕ . 1 2 1−ϕα r(K)+ϕα m 1−ϕβ s(0)−ϕβ cKφ(K)+ϕβ m 2 6 1 2 4 6 2 Therefore, λ < 1 iff m > s(0)+cKφ(K). Besides, due to r′(K) < 0 there holds λ < 1. On 2 2 1 the other hand we have 2+ϕT 2+ϕT 1 2 λ +1 = , λ +1 = , 1 2 1−ϕα r(K)+ϕα m 1−ϕβ s(0)−ϕβ cKφ(K)+ϕβ m 2 6 1 2 4 6 2 Therefore, if (15) is fulfilled then λ > −1, λ > −1. In result we have |λ | < 1 and |λ | < 1. 1 2 1 2 Hence, by Lemma 1 the proposition is proved. 4.3. The equilibrium point P∗ = (0,M) (the equilibrium point of extinction of the prey species) 2 Proposition 5. For the case m > r(0)−Mφ(0) and m < s(0) consider the difference scheme 1 2 (3) under the assumptions of Proposition 1. If additionally assume that T := r(0)−Mφ(0)−m −2α r(0)+2α Mφ(0)+2α m > 0, 3 1 2 4 6 1 (16) ′ T := Ms(M)−2β s(M)+2β m > 0, 4 2 6 2 then the equilibrium point of the form P∗ = (0,M) is locally asymptotically stable if m > 2 1 r(0)−Mφ(0) and m < s(0), and unstable otherwise. 2 10