Next Article in Journal
Efficient Approaches for Solving Systems of Nonlinear Time-Fractional Partial Differential Equations
Next Article in Special Issue
Comparative Numerical Study of Spline-Based Numerical Techniques for Time Fractional Cattaneo Equation in the Sense of Caputo–Fabrizio
Previous Article in Journal
On Starlike Functions of Negative Order Defined by q-Fractional Derivative
Previous Article in Special Issue
A New Three-Step Root-Finding Numerical Method and Its Fractal Global Behavior
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Discretization, Bifurcation, and Control for a Class of Predator-Prey Interactions

1
Department of Basic Sciences and Humanities, College of Computer and Information Sciences, Majmaah University, Al-Majmaah 11952, Saudi Arabia
2
Department of Mathematics, University of the Poonch Rawalakot, Azad Kashmir 12350, Pakistan
*
Author to whom correspondence should be addressed.
Fractal Fract. 2022, 6(1), 31; https://doi.org/10.3390/fractalfract6010031
Submission received: 30 November 2021 / Revised: 22 December 2021 / Accepted: 27 December 2021 / Published: 6 January 2022

Abstract

:
The present study focuses on the dynamical aspects of a discrete-time Leslie-Gower predator-prey model accompanied by a Holling type III functional response. Discretization is conducted by applying a piecewise constant argument method of differential equations. Moreover, boundedness, existence, uniqueness, and a local stability analysis of biologically feasible equilibria were investigated. By implementing the center manifold theorem and bifurcation theory, our study reveals that the given system undergoes period-doubling and Neimark-Sacker bifurcation around the interior equilibrium point. By contrast, chaotic attractors ensure chaos. To avoid these unpredictable situations, we establish a feedback-control strategy to control the chaos created under the influence of bifurcation. The fractal dimensions of the proposed model are calculated. The maximum Lyapunov exponents and phase portraits are depicted to further confirm the complexity and chaotic behavior. Finally, numerical simulations are presented to confirm the theoretical and analytical findings.

1. Introduction

predator-prey models have a wide range of applications in ecological and biological fields. Although various fundamental aspects of the nonlinear dynamics of predator-prey population models related to continuous dynamical systems have been studied, the characteristics of discrete dynamical systems remain comparatively unknown. A discrete dynamical structure possesses a solitary dynamical nature as compared to a continuous system. There are several critical and practical problems in daily life that can be characterized with the help of a discrete dynamical system. To consider the analytical aspects of a solution that is difficult to calculate, various schemes can be implemented to discretize a continuous system and discuss the numerical solution. Therefore, detailed critical inspections of discrete-time dynamical systems have contributed immensely to various fields such as engineering, physics, chemistry, and mathematics. There have been numerous studies conducted that are related to the dynamics of predator-prey models.
Chen et al. [1] applied the Euler scheme and center manifold theorem to a ratio-dependent prey-predator model and scrutinized the dynamic characteristics of the model. Ghaziani et al. [2] studied a prey-predator system with a Holling functional response and discussed the resonance and bifurcation analyses. Jana [3] found extremely powerful dynamical conditions through numerical and theoretical investigations of discrete-time prey-predator models, such as stability conditions, flip, and hopf-bifurcation. Misra et al. [4] studied a predator-prey model based on age predation and discussed the dynamic behavior of the models. Zhang et al. [5] presented a biological economic system related to the predator-prey model of a differential algebraic system by applying a new normal form. Hu and Cao [6] investigated the Holling and Leslie type predator-prey model and discussed a chaos and bifurcation analysis. Wang and Li [7] proposed a lemma that is extremely meaningful for discussing the stability and bifurcation of the systems. The fundamental finding in the dynamics of prey-predator species is the classical Lotka-Volterra prey-predator model, which exhibits unrealistic behavior (see, Murdoch et al. [8]). To remove such unrealistic behavior, Holling introduced three different types of functional responses (see, Holling [9]). Rosenzweig and MacArthur [10] implemented a functional response to modify the predator-prey model. An investigation into population interaction focused on the continuous dynamical system of two species [11,12,13]. By contrast, a recent study led to the discrete dynamical system becoming more suitable than a continuous version when the population is non-overlapping (e.g., see, Jing et al. [14], Liu et al. [15], Lopez-Ruiz and Fournier-Prunaret [16], Neubert and Kot [17]). Furthermore, multiple existing studies related to the dynamics of predator-prey models are described in [18,19,20,21,22,23,24,25,26]. In [27], the Holling type-III functional response was introduced in both populations (prey and predator). The stability conditions around biologically suitable equilibria were further discussed. Diagrams of the phase portraits, bifurcation, and time series were plotted. It was shown that the system is sensitive to the initial conditions, which means that the system is chaotic. A two-dimensional continuous model with a Holling-III functional response in both prey and predator was presented [28]. Furthermore, Euler’s scheme was used to discretize the model and study the complex behavior of the system. Elettreby et al. [29] discussed a discrete-time prey-predator model with predator and prey populations having Holling type I and III functional responses, respectively. Moreover, they described a fascinating dynamical nature of the model, including stability, bifurcation, and chaos, which ensure the rich dynamics of discrete-time models.
In this study, we evaluate the specific prey –predator model discussed by Murray [30]:
d x d t = x r 1 x k a x y b 2 + x 2 , d y d t = y s 1 h y x ,
where x t and y t denote the densities of prey and predator species at any time t , respectively; the carrying capacity of prey in the absence of predator is k , and r ,     b ,     a ,     s ,   and   h are positive constants. Moreover, the carrying capacity is proportional to the prey population size and population of prey attacked by predators, as specified by the Holling type III functional response. He and Lai [31] examined the bifurcation and chaos control of the discrete-time version of model (1) by implementing Euler’s forward scheme with step size h as the bifurcation parameter. The numerical results in [31] show that period-doubling bifurcation occurs when a large step size is considered in Euler’s method; this fact contravenes the precision of the numerical method for discretization. To overcome this deficiency, the following discretization method was implemented. Considering the regular time interval for the average growth rate in both populations, by resorting to piecewise constant arguments for solving nonlinear differential equations, system (1) can then be rewritten as follows:
1 x t d x t d t = r 1 x t k a x t y t b 2 + x [ t ] 2 , 1 y t d y t d t = s 1 h y t x t ,
where the integer part of t is given by t within the interval 0 < t < 1 . In addition, by integrating system (1) for t n : n + 1 , ( n = 0 , 1 , 2 , ), we have the following system:
x t = x n exp r 1 x n k a x n y n b 2 + x n 2 t n , y t = y n exp s 1 h y n x n t n .
Applying t n + 1 , we obtain the following prey-predator system:
x n + 1 = x n exp r 1 x n k a x n y n b 2 + x n 2 , y n + 1 = y n exp s 1 h y n x n .
The key contributions and findings of the current study are as follows for model (4):
  • The existence and uniqueness of biologically feasible equilibria and their stability analysis are discussed.
  • Our findings indicate that model (4) undergoes periodic doubling as well as a Neimark-Sacker bifurcation at its unique positive equilibrium.
  • The direction and existance criteria for both types of bifurcation are examined under interior equilibrium.
  • A hybrid control strategy is applied to control the chaos in model (4).
The remainder of this paper is organized as follows. After presenting some related preliminaries in Section 2, the boundedness of the steady state is analyzed in Section 3. In Section 4, the dynamics of system (4), including the existence of equilibria and local stability, are presented. Section 5 describes an investigation of the birfurcation analysis at the interior fixed point of system (4). In Section 6, we study a hybrid control method to control the chaos. The fractal dimensions are calculated in Section 7. Finally, numerical simulations are provided in Section 8 to verify our analytical approach. Conclusions related to these results are presented in Section 9 and the future directions are providing in Section 10.
Furthermore, a detailed investigation of some charismatic population models and their qualitative behavior are provided (see, Din and Din et al. [18,19,20,21,22,23,24,25,26] and the references therein).

2. Preliminaries

Definition 1.
([32]) A point x * is said to be a fixed point of the map for an equilibrium point if f x * = x * .
Theorem 1.
([32]) Let f   :   I     I be a continuous map, where I = a ,   b is a closed interval in R . Then,   f has a fixed point.
Theorem 2.
([32]) Let f   :   I = a ,   b     R be a continuous map such that f I     I . Then, f has a fixed point in I .
Definition 2.
([32]) Let f   :   I     I be a map and x * be a fixed point of f , where I is an interval in the set of real numbers R . Then, the following conditions hold true:
  • x * is said to be stable if for any ε > 0 , there exists δ > 0 such that for all x 0 I with x 0 x * < δ we have f n x 0 x * < ε for all n     Z + . Otherwise, the fixed point x * is unstable.
  • x * is said to be attractive if η > 0 exists, such that x 0 x * < η implies lim n f n x 0 = x * .
  • x * is asymptotically stable if it is both stable and attractive. If in (2), η = , then x * is said to be globally asymptotically stable.
Definition 3.
([32]) A fixed point x * of a map f is said to be hyperbolic if f x * = 1 . Otherwise, it is non-hyperbolic.
Theorem 3.
([32]) Let x * be a hyperbolic fixed point of a map f , where f is continuously differentiable at x * . The following statements hold true:
  • If f x * < 1 , then x * is asymptotically stable.
  • If f x * > 1 , then x * is unstable.

3. Boundedness

The boundedness of system (4) is based on the following Remark.
Remark 1.
([25]) Assuming that x 0 > 0 for every x t and x t + 1 x t e x p A 1 B x t for every t t 1 , , where B > 0 is constant. Then,
lim n Sup x t 1 A B exp A 1
Using Remark 1, we state the following theorem for the uniform boundedness of system (4).
Theorem 4.
Any positive solution x n , y n of model (4) is uniformly bounded.
Proof. 
Assuming that x n , y n is any positive solution of system (4), we then have
x n + 1 x n exp r 1 x n k ,   for   all   n = 0 , 1 , 2 , .
Let x 0 > 0 . Using Remark 1, we obtain the following result.
lim n Sup x n k r exp r 1 = l 1 .
Furthermore, from the second part of system (4), we obtain the following:
y n + 1 y n exp s 1 h y n l 1 .
Let y 0 > 0 . Applying Remark 1, we obtain the following result:
lim n Sup y n l 1 s h s 1 = l 2
Thus, it follows that lim n Sup x n , y n l , where l = max l 1 , l 2 . The proof is completed. □

4. Existence of a Positive Fixed Point and Local Stability

To explore the existence of a fixed point of model (4), suppose that x , y is any arbitrary fixed point of (4). Then, x , y must satisfy the following algebraic system of equations:
x = x   exp r 1 x k a x y b 2 + x 2 , y = y   exp s 1 h y x .
Then, (7) has a boundary equilibrium point k , 0 . In addition, we also explore the existence and uniqueness of the solution of system (4) because the positive fixed points are not in a closed form. For this purpose, the following computation, using Theorem 4, exhibits the existence and uniqueness of the solution to model (4).
Theorem 5.
There exists a unique positive steady-state x * ,     y * 0 ,     l 1 × 0 ,     l 2 of system (4).
Proof. 
To attain the fixed point by solving system (7), we have
r 1 x k = a x y x 2 + b 2 , x = y h .
Suppose that
F x = r 1 x k a x 2 h x 2 + b 2
for all x 0 ,     l 1 . Then, we can see that F 0 = r > 0 and
F l 1 = a exp 2 r 1 k 2 h b 2 + exp 2 r 1 k 2 r 2 r 2 + r 1 exp r 1 r < 0
for all a ,     b ,     s ,   r ,     k ,     and   h > 0 . Hence, there exists at least one root of F x = 0 , for x 0 ,     l 1 . In addition,
F ' x = r k 2 a b 2 x h b 2 + x 2 2 < 0
for all x 0 ,     l 1 . Therefore, the system (4) has a unique positive fixed point x * ,     y * 0 ,     l 1 × 0 ,     l 2 .
Initially, we explored the stability analysis of the boundary equilibrium k , 0 . The Jacobian matrix F J evaluated at k , 0 , is expressed as
F J k , 0 = 1 r k 2 a b 2 + k 2 0 exp s ,  
and the characteristic equation computed at k , 0 is given by
F η = η 2 1 r + exp s η + 1 r exp s  
Hence, F η = 0 has two roots, namely, η 1 = exp s and η 2 = 1 r . In addition, k , 0 is the source if r > 2 , and is the saddle point if 0 < r < 2 . Next, we explored the stability analysis of the fixed points. To investigate the stability of the equilibria, we calculated the Jacobian F J of system (4) at any point x , y as follows:
F J x , y = b 11 b 12 b 21 b 22
The characteristic polynomial of F J at x , y is given by
η = η 2 T 1 η + D 1 ,
where
T 1 = b 11 + b 22 ,
and
D 1 = b 11 b 22 b 12 b 21
The following Lemma is extremely useful to examine the stability of the equilibria. □
Lemma 1.
Let η = η 2 T 1 η + D 1 and 1 > 0 . Moreover, η 1 , η 2 are the roots of equation η = 0 , and thus
 (i)
η 1 < 1 & η 2 < 1     1 > 0 and D 1 < 1 ;
 (ii)
η 1 < 1 & η 2 > 1 or ( η 1 > 1 and η 2 < 1 )   1 < 0 ;
 (iii)
η 1 > 1 & η 2 > 1     1 > 0 and D 1 > 1 ;
 (iv)
η 1 = 1 & η 2 1   1 = 0 and T 1 0 , 2 ;
 (v)
η 1 and η 2 are complex and η 1 = 1 and η 2 = 1     T 1 2 4 D 1 < 0 and D 1 = 1 .
Because η 1 and η 2 are the eigenvalues of (9), the following topological results are obtained.
The equilibrium x , y is known as a sink if η 1 < 1   and   η 2 < 1 , which is locally asymptotically stable, and as a source if η 1 > 1   and   η 2 > 1 ; thus, the nature of the source is always unstable. Moreover, the equilibrium point x , y is always known as the saddle point if η 1 < 1   and   η 2 > 1 or ( η 1 > 1   and   η 2 < 1 ) . In the case of a non-hyperbolic equilibrium x , y , either η 1 = 1   or   η 2 = 1 .
Our next aim is to discuss the local stability of the unique positive equilibrium x * , y * of system (4). Let (9) be the characteristic polynomial of the variational matrix evaluated at x * , y * , such that
T 1 = 2 r x * k Ω s )     and     D 1 = 1 r x * k Ω 1 s + s Φ h
where Ω = a x * y * b 2 x * 2 b 2 + x * 2 2 and Φ = a x * 2 b 2 + x * 2 . Thus, by applying Lemma 1, we discuss the local stability of system (4) around x * , y * by stating the following proposition.
Proposition 1.
The interior equilibrium point x * , y * of system (4) satisfies the following results:
 (i)
The interior equilibrium x * , y * is stable iff:
2 r x * k Ω s < 1 + 1 r x * k Ω 1 s + s Φ h ,
and
s Φ h + 1 s 1 r x * k Ω > 1
 (ii)
The positive fixed point x * , y * is a saddle point if and only if
2 r x * k Ω s 2 > 4 1 s 1 r x * k Ω + s Φ h ,
and
 (iii)
The interior fixed point x * , y * is non-hyperbolic if and only if
2 r x * k Ω s = 1 + 1 r x * k Ω 1 s + s Φ h
or
1 r x * k Ω 1 s + s Φ h = 1     a n d     2 r x * k Ω s < 2 .
To explore the local stability criteria for x * , y * of model (4), we have the following theorem:
Theorem 6.
If neither (10) nor (11) is satisfied, then the positive steady-state x * , y * of system (4) is locally asymptotically stable if and only if the following condition is satisfied.
2 r x * k Ω s < 1 + 1 r x * k Ω 1 s + s Φ h < 2

5. Bifurcation Analysis

In this section, we discuss the period-doubling and Neimark-Sacker bifurcations of system (4) around the interior equilibrium. Initially, we explored the period-doubling bifurcation at a positive fixed point x * , y * of system (4). To study the period-doubling bifurcation, assume that T 1 2 > 4 D 1 , that is,
2 r x * k Ω s > 4 1 r x * k Ω 1 s + s Φ h
and T 1 + D 1 + 1 = 0 . It then follows that
s : = 2 h r x * 2 Ω k h Ω   k + r x * + k Φ 2   h
Then, η 1 = 1 and η 2 1 if
1 r x * k Ω 1 s + s Φ h ± 1 .
Consider the map T P B = { a , b , k , r , s R + 5 for which (12)–(14) are thus satisfied. Then, the equilibrium x * , y * of system (4) sustains period-doubling bifurcation whenever the parameters deviate within the small neighborhood of T P B . Thus, system (4) along with parameters a , b , k , r , s 1 T P B , can be written as follows:
x y x e r 1 x k a x y b 2 + x 2 y e s 1 1 h y x
The following perturbation of system (15) can be obtained by taking s ¯ as a bifurcation parameter:
x y x e r 1 x k a x y b 2 + x 2 y e s 1 + s ¯ 1 h y x
where s ¯ < < 1 denotes the least perturbation parameter. Assuming that N = x x * ,     P = y y * , system (16) is reduced to the following form:
N P b 11 b 12 b 21 b 22 N P + f 1 N , P , s ¯ f 2 N , P , s ¯
Here,
f 1 N , P , s ¯ = b 13 N 2 + b 14 N P + b 15 P 2 + a 1 N 3 + a 2 N 2 P + a 3 N P 2 + a 4 P 3 + O ( N + P + s ¯ ) 4 , f 2 N , P , s ¯ = b 23 N 2 + b 24 N P + b 25 P 2 + d 1 N 3 + d 2 N 2 P + d 3 N P 2 + d 4 P 3 + c 1 s ¯ N + c 2 r ¯ P + c 3 r ¯ 2 + c 4 r ¯ N P + c 5 r ¯ N 2 + c 6 r ¯ P 2 + c 7 r ¯ 2 N + c 8 s ¯ 2 P + c 9 r ¯ 3 + O ( N + P + s ¯ ) 4
Whereas the descriptions and computations of the involved coefficients are given in Appendix A.
The canonical form of (17) at s 1 = 0 can be obtained by assuming the following map:
N P = b 12 b 12 1 b 11 η 2 b 11 u v .
The normal form of system (17) under translation (18) can be expressed as
u v 1 0 0 η 2 u v + f ˜ u , v , s ¯ g ˜ u , v , s ¯ ,
the descriptions and computations of the involved functions and parameters leading to the following expression are provided in Appendix B.
F 1 : = 2 f ˜ u s ¯ + 1 2 F r ¯ 2 F u 2 0 , 0 = c 2 1 + b 11 η 2 + 1 c 1 b 12 η 2 + 1 , F 2 : = 1 6 3 F u 3 + 1 2 2 F u 2 2 0 , 0 = t 1 2 + t 5
Hence, we arrive at the following conclusions based on the aforementioned calculations.
Theorem 7.
There exists a period-doubling bifurcation at x * , y * of system (4), whenever F 2 0 and r deviates within a small neighboring point of s 1 . In addition, if F 2 > 0 ,   F 2 < 0 , the orbit is period-2 stable (unstable).
Next, we investigated the Neimark-Sacker bifurcation around x * , y * of system (4). For identical results, we referred to the studies by Din [24,25], Shabbir et al. [20], and Jing et al. [14]. Furthermore, the equilibrium point moves around the close invariant curve, owing to the Neimark-Sacker bifurcation. To explore the Neimark-Sacker bifurcation, we find the conditions for which x * , y * is a non-hyperbolic point with a complex conjugate root of the characteristic equation of the unit modulus. Thus, if the following results hold true, then η = 0 has two complex conjugate roots with a unit modulus.
s : = Ω   k + r x * h h Ω   k + r x * + k Φ h ,
and
2 r x * k Ω s < 2 .
Consider
T N S = { a , b , k , r , s R + 5 : s = Ω   k + r x * h h Ω   k + r x * + k Φ h     and     2 r x * k Ω s < 2 } .
Assuming that s 2 = Ω   k + r x * h h Ω   k + r x * + k Φ h , the fixed point x * , y * ensures the Neimark-Sacker bifurcation when the parameters fluctuate in the least neighborhood of T N S . Thus, system (4) along with parameters a , b , k , r 2 , s can be expressed as follows:
x y x e r 1 x k a x y b 2 + x 2 y e s 2 1 h y x .
The following perturbation of system (20) can be obtained by taking s ˜ as the bifurcation parameter, i.e.,
x y x e r 1 x k a x y b 2 + x 2 y e s 2 + s ˜ 1 h y x
where s ˜ < < 1 denotes the least perturbation. Assuming that N = x x * ,     P = y y * , then system (21) takes the following modified form:
N P a 11 a 12 a 21 a 22 N P + g 1 N , P g 2 N , P ,
where
g 1 N , P = b 13 N 2 + b 14 N P + b 15 P 2 + a 1 N 3 + a 2 N 2 P + a 3 N P 2 + a 4 P 3 + O ( N + P + r ¯ ) 4 , g 2 N , P = b 23 N 2 + b 24 N P + b 25 P 2 + d 1 N 3 + d 2 N 2 P + d 3 N P 2 + d 4 P 3 + O ( N + P + r ¯ ) 4 .
Here, b 11 ,   b 12 ,   b 21 ,   b 22 , b 13 ,   b 14 ,   b 15 ,   a 1 ,   a 2 ,   a 3 ,   a 4 ,   b 23 ,   b 24 ,   b 25 ,   d 1 ,   d 2 ,   d 3 ,   and   d 4 , are defined in (17) by replacing s 1 by s 2 + s ˜ . Let
η 2 T 1 s ˜ η + D 1 s ˜ = 0 ,
be the characteristic equation of the variational matrix of (22) evaluated at 0 , 0 , where
T 1 s ˜ = 2 r x * k Ω s 2 + s ˜     and     D 1 s ˜ = 1 r x * k Ω 1 s 2 + s ˜ + s 2 + s ˜ Φ h
where Ω = a x * y * b 2 x * 2 b 2 + x * 2 2 and Φ = a x * 2 b 2 + x * 2 . Because a , b , k , r , s 2 T N S , η 1 = η 2 such that η 1 and η 2 are the complex conjugate roots of (23), it follows that
η 1 , η 2 = T 1 s ˜ 2 ± i 2 4 D 1 s ˜ T 1 2 s ˜
We then obtain
η 1 = η 2 = D 1 s ˜ , d D 1 s ˜ d s ˜ s ˜ = 0 = Ω   h + Φ h k + h r x 2 Ω 1 k + r x s 1 h + Φ   k s 0
Moreover, T 1 0 = 2 r x * k Ω s 2 0 ,   1 . Because a , b , k , r , s 2 T N S , it follows that 2 < T 1 0 = 2 r x * k Ω s 2 < 2 . Thus, we have η 1 m , η 2 m 1 for all m = 1 ,   2 ,   3 ,   4 at s ˜ = 0 , for T 1 0 0 , 1 , ± 2 . Hence, for s ˜ = 0 , zeros of (23) do not belong to the intersection of the unit circle with coordinate axes if the following condition is satisfied:
2 Ω s 2 r x * k ,     3 Ω s 2 r x * k
The canonical form of (22) at s ˜ = 0 can be obtained by taking γ = T 1 0 2 , δ = 1 2 4 D 1 0 T 1 2 0 and assuming
N P = b 12 0 γ b 11 δ u v
Using transformation (25), we obtain the following canonical form of system (22):
u v γ δ δ γ u v + f ˜ u , v g ˜ u , v
where
f ˜ u , v = a 1 N 3 b 12 + a 2 N 2 P b 12 + b 13 N 2 b 12 + a 3 N P 2 b 12 + b 14 N P b 12 + a 4 P 3 b 12 + b 15 P 2 b 12 + O ( u + v ) 4
g ˜ u , v = γ b 11 a 1 b 12 δ d 1 δ N 3 + γ b 11 a 2 b 12 δ d 2 δ N 2 P
+ γ b 11 b 13 b 12 δ b 23 δ N 2 + γ b 11 a 3 b 12 δ d 3 δ N P 2
+ γ b 11 b 14 b 12 δ b 24 δ N P + γ b 11 a 4 b 12 δ d 4 δ P 3
+ γ b 11 b 15 b 12 δ b 25 δ P 2 + O ( u + v ) 4
N = a 12 u   and   < P = γ b 11 u δ v
Owing to the aforementioned computation, we state a nonzero real number
L : = R e 1 2 η 1 η 2 2 1 η 1 ξ 20 ξ 11 1 2 ξ 11 | 2 ξ 02 | 2 + R e η 2 ξ 21 c ˜ = 0
where,
ξ 11 = 1 4 f ˜ u u + f ˜ v v + i g ˜ u u + g ˜ v v ,  
ξ 02 = 1 8 f ˜ u u f ˜ v v 2 g ˜ u v + i g ˜ u u g ˜ v v + 2 f ˜ u v
ξ 20 = 1 8 f ˜ u u f ˜ v v + 2 g ˜ u v + i g ˜ u u g ˜ v v 2 f ˜ u v ,  
ξ 21 = 1 16 f ˜ u u u + f ˜ u v v + g ˜ u u v + g ˜ v v v + i g ˜ u u u + g ˜ u v v f ˜ u u v f ˜ v v v .
Ultimately, we deduced the following conclusions for the direction and existence of the Neimark-Sacker bifurcation, based on the aforementioned calculation:
Theorem 8.
There exists a Neimark-Sacker bifurcation around x * , y * whenever s deviates wtihin the neighborhood of s 2 = Ω   k + r x * h h Ω   k + r x * + k Φ h . In addition, if L < 0   ( L > 0 ) , then an attracting (or repelling) invariant closed curve fluctuates in the range x * , y * for s > s 2   ( o r   s < s 2 ) .

6. Chaos Control

In this section, we implement the hybrid control method for controlling the chaos caused by the period-doubling bifurcation and for controlling the Neimark-Sacker bifurcation in (4). Such strategies have been discussed elsewhere in [21,33,34,35,36,37]. We assume the following controlled system corresponding to model (4):
x n + 1 = ϵ x n exp r 1 x n k a x n y n b 2 + x n 2 + 1 ϵ x n , y n + 1 = ϵ y n exp s 1 h y n x n + 1 ϵ y n ,
where 0 < ϵ < 1 . Furthermore, both types of bifurcations can be controlled by choosing an appropriate value of parameter ϵ . The controlled system (27) and the original system (4) have the same equilibrium point; the Jacobian matrix of the controlled system (27) at x * , y * is expressed by
1 ϵ   x r k ϵ   Ω ϵ   Φ ϵ   s h 1 ϵ   s
Consequently, the necessary and sufficient condition for local stability around x * , y * of the controlled system (27) yields the following result.
Theorem 9.
The interior fixed point x * , y * of (27) is locally asymptotically stable if
2 ε r x * k ε Ω ε s < 1 + 1 ε r x * k ε Ω 1 ε s + s Φ h < 2 .

7. Fractal Dimension

The fractal dimension that describes the strange attractors of discrete-time models is defined as follows [38,39]:
D 𝓁 = 𝓂 + j = 1 𝓂 𝓀 𝓂 𝓀 𝓂 ,
where 𝓀 1 ,   𝓀 2 , , 𝓀 𝓂 are Lyapunov exponents, and 𝓂 is the largest integer such that j = 1 𝓂 𝓀 𝓂 0 and j = 1 𝓂 + 1 𝓀 𝓂 < 0 . For the discrete-time model (4), the fractal dimension takes the following form:
D 𝓁 = 1 + 𝓀 1 𝓀 2 ,         𝓀 1 > 0 > 𝓀 2
Furthermore, for the values of parameters a ,   b ,   r ,   k ,   h , and s , the two Lyapunov exponents F 1 and F 2 are computed numerically. If b = 3.3 ,     a = 0.8 ,     r = 1.3 ,     h = 2.7 ,   and     k = 1.8 , then F 1 and F 2 corresponding to the values of the bifurcation (period-doubling) parameter s from the chaotic region, with the help of Mathematica software, are shown in Table 1.
The strange attractors for fixed parametric values illustrate that the discrete model (4) has a complex dynamical behavior as parameter s increases by s > 2.1894756175566834 . Similarly, for the Neimark-Sacker bifurcation, the Lyapunov exponents and fractal dimension can be calculated for the values of parameter s from the chaotic region, that is, s = 1.66 ,   1.68 ,   1.89 , and so on. The strange attractors corresponding to these values are also shown in Figure 1. In particular, Figure 1g,h,k below, demonstrate that the discrete time model (4) has a complex dynamical nature when parameter s > 1.3874082082631611 .

8. Numerical Simulation

This section verifies the aforementioned theoretical discussion. The first example is related to the existence and direction of the Neimark-Sacker bifurcation. The second example shows that for a suitable choice of parameters, system (4) undergoes period-doubling bifurcation. Moreover, to confirm the control of flip and Neimark-Sacker bifurcation, we provide two examples for different choices of parameters defined in T P B and T N S .
Example 1.
Let b = 1.7 ,   a = 2.26 ,   r = 2.4 ,   h = 1.8 ,   k = 1.4 , s 1 , 1.8 , and initial conditions x 0 , y 0 = 1.0219 , 0.567721 . Then, both species undergo a Neimark-Sacker bifurcation, as shown in Figure 2. To confirm the chaotic behavior of model (4) MLE is shown in Figure 2c.
Furthermore, Figure 1a–l shows the interesting behavior of system (4). Figure 1f–l shows the chaotic behavior of system (4). Assuming that b = 1.7 , we have a positive fixed point x * , y * = 1.165750001 , 0.6476388892 , which loses stability and undergoes a Neimark-Sacker bifurcation. Thus, for the aforementioned parameters, we have the following control system:
x n + 1 = ϵ   x n e 2.4 1.714285714   x n 2.26   x n y n x n 2 + 2.89 + 1 ϵ x n , y n + 1 = ϵ   y n e 1.5 2.70   y n x n + 1 ϵ y n .
It can be clearly observed that the controlled system (28) has a unique positive equilibrium point x * , y * = 1.165750001 ,   0.6476388892 , which is similar to the original system (4). In addition, the Jacobian at equilibrium x * , y * = 1.165750001 , 0.6476388892 has the following form:
2.143126283   ϵ + 1 0.7228285706   ϵ 0.8333333325   ϵ 1.500000000   ϵ + 1 .
Example 2.
Assuming the parameters b = 3.3 , a = 0.8 , r = 1.3 , h = 2.7 , k = 1.8 , and s 2 ,   3.5 , and the initial conditions x 0 , y 0 = 1.74693 , 0.647013 , both species then undergo period-doubling bifurcation when the bifurcation parameter passes through s = 2.1894756175566834 , as shown in Figure 3. In particular this fact is obvious in Figure 3a,b. Moreover, to confirm the chaotic behavior of model (4) MLE is shown in Figure 3c.
Furthermore, if s = 2.4 , then the equilibrium point x * , y * = 1.712924751 , 0.6344165743 loses its stability and undergoes periodic doubling (see Figure 4).
Thus, for the aforementioned parameters, we present the following control system:
x n + 1 = ϵ   x n e 1.3 0.7222222223   x n 0.8   x n y n x n 2 + 10.89 + 1 ϵ x n , y n + 1 = ϵ   y n e 2.4 6.48   y n x n + 1 ϵ y n ,
The fixed point x * , y * = 1.712924751 ,   0.6344165743 was preserved in the case of a controlled system (29). Furthermore, the variational matrix of the aforementioned controlled system computed at a fixed point x * , y * = 1.712924751 ,   0.6344165743 is given by
1.273304693   ϵ + 1 0.1697967362   ϵ 0.8888888882   ϵ 2.399999999   ϵ + 1
The characteristic polynomial of the aforementioned Jacobian matrix is given by
η 2 + 3.673304692   ϵ 2 η + 3.206861694   ϵ 2 3.673304692   ϵ + 1 = 0 .
According to Lemma 1, the control system is locally asymptotically stable, if 0 < ϵ < 0.8910230450195268 and bifurcation is controlled for 0 < ϵ < 0.8910230450195268 (see Figure 4c,d).
Finally, some local implications of the MLE diagrams, shown in Figure 1c and Figure 2c for the Neimark-Sacker bifurcation and period-doubling bifurcation, respectively, are plotted in Figure 5a,b, respectively. It has also been verified that the system undergoes Neimark-Sacker bifurcation at s = 1.3874082082631611 , where the phase portrait at this point shows a closed invariant curve, as already shown in Figure 4c.

9. Concluding Remarks

In this study, we examined the qualitative and dynamical analyses of a discrete-time predator-prey model. Piecewise constant arguments have been applied to achieve the discrete-time counterpart of a continuous model. Thus, a comprehensive analysis of model (4) was presented. In particular, we investigated the boundedness, local stability of the boundary, and positive equilibrium points, which seem to present more challenging cases of Euler’s discretization scheme in [31]. Moreover, it was proved that the population sustains both period-doubling bifurcation and Neimark-Sacker bifurcation near the interior equilibrium. The parametric conditions were obtained for the direction and existence of both types of bifurcations using the theory of bifurcation and center manifold theorem. Moreover, the chaotic attractors shown in Figure 5 ensures chaos in the system. To control the chaotic behavior of system (4), a hybrid control method was implemented. Hence, by applying a control strategy, both types of bifurcations can be controlled for a maximum range of control parameters. We also presented the fractal dimension of model (4), which characterizes the strange attractors provided in Figure 5 thereby illustrating the complexity and rich dynamics of discrete model (4). Finally, numerical simulations were conducted to verify the analytical and theoretical approaches.

10. Future Direction

Our future research will include the Leslie–Gower predator-prey model with the functional response of Holling type-II. In this case, we aim to conduct stability, bifurcation, and chaos-control analyses of the model. A comparison of both functional responses will be conducted in a future study.

Author Contributions

Conceptualization, A.T., M.S.S. and Q.D.; methodology, A.T., M.S.S. and Q.D.; validation, A.T., M.S.S., Q.D. and H.N.; formal analysis, A.T., M.S.S., Q.D. and H.N.; investigation, A.T., M.S.S., Q.D. and H.N.; data curation, A.T., M.S.S., Q.D. and H.N.; writing—original draft preparation, A.T., M.S.S., Q.D. and H.N.; writing—review and editing, A.T., M.S.S., Q.D. and H.N. All authors have read and agreed to the published version of the manuscript.

Funding

This research did not receive any specific external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

Asifa Tassaddiq would like to thank the Deanship of Scientific Research at Majmaah University for supporting this work under Project No. R-2021-274. The authors are also thankful to the worthy reviewers and editors for their useful and valuable suggestions for the improvement of this paper which led to a better presentation.

Conflicts of Interest

None of the authors have any competing interests in this manuscript.

Appendix A

b 11 = 1 + x * B 1 ,     b 12 = a x * 2 b 2 + x * 2 ,     b 21 = s 1 h , b 22 = 1 s 1 ,     b 15 = a 2 x * 3 2 b 2 + x * 2 2 b 13 = B 1 + 1 2 x * y * B 3 + 1 2 x * B 1 2 , b 14 = x * a x * B 1 b 2 B 2 x * 2 B 2 + a b 2 + x * 2 , a 1 = y * B 3 + 1 2 B 1 2 + x * a y * b 4 6 b 2 x * 2 + x * 4 b 2 + x * 2 4 + 1 2 x * y * B 3 B 1 + 1 6 x * B 1 3 , a 2 = B 2 B 1 a x * b 2 + x * 2 + 1 2 x * B 3 1 2 x * 2 y * B 3 a b 2 + x * 2 + x * B 1 B 2 1 2 x * 2 B 1 2 a b 2 + x * 2 , a 3 = 1 2 a 2 x * 2 b 2 + x * 2 2 x * 2 B 2 a b 2 + x * 2 + 1 2 x * 3 B 1 a 2 b 2 + x * 2 2 , a 4 = 1 6 x * 4 a 3 b 2 + x * 2 3 , b 23 = h s 1 y * 2 h s 1 y * 2 x * 2 x * 4 ,         b 24 = h s 1 y * 2 x * h s 1 y * x * 3 , b 25 = h s 1 h s 1 y * 2 x * 2 x * 2 ,         d 1 = h s 1 y * 2 6 x * 2 6 h s 1 x * y * + h 2 s 1 2 y * 2 6 x * 6 , d 2 = h s 1 y * h s 1 y * 4 x * h s 1 y * x * 2 x * 5 , d 3 = h s 1 2 x * 2 4 h s 1 x * y * + h 2 s 1 2 y * 2 2 x * 4 , d 4 = h 2 s 1 2 h s 1 y * 3 x * 6 x * 3 ,         c 1 = h y * 2 1 + s 1 x * h s 1 y * x * 3 , c 2 = x * 2 h 2 + s 1 x * y * + h 2 s 1 y * 2 x * 2 , c 3 = h y * 2 1 + s 1 x * 2 h s 1 4 + s 1 x * y * + h 2 s 1 2 y * 2 x * 4 , c 4 = h y * 2 2 1 + s 1 x * 2 h s 1 4 + s 1 x * y * + h 2 s 1 2 y * 2 2 x * 5 , c 5 = h 2 1 + s 1 x * 2 h s 1 4 + s 1 x * y * + h 2 s 1 2 y * 2 2 x * 3 .
where B 1 = r 1 k + a y * b 2 x * 2 b 2 + x * 2 2 , B 2 = a b 2 x * 2 b 2 + x * 2 2 , B 3 = 2 a x * 3 b 2 x * 2 b 2 + x * 2 3 , B 4 = y * B 3 .

Appendix B

f ˜ u , v , s ¯ = b 13 η 2 b 11 b 12 1 + η 2 b 23 1 + η 2 N 2 b 24 1 + η 2 + b 14 b 11 η 2 b 12 1 + η 2 N P b 25 1 + λ 2 + b 15 b 11 η 2 b 12 1 + η 2 P 2 d 1 1 + λ 2 + a 1 b 11 η 2 b 12 1 + η 2 N 3 d 2 1 + λ 2 + a 2 b 11 η 2 b 12 1 + η 2 N 2 P d 3 1 + η 2 + a 3 b 11 η 2 b 12 1 + η 2 N P 2 d 4 1 + η 2 + a 4 b 11 η 2 b 12 1 + η 2 P 3 c 1 1 + η 2 s ¯ N c 2 1 + η 2 r ¯ P c 3 1 + η 2 N P c 4 1 + η 2 r ¯ N 2 c 5 1 + η 2 r ¯ P 2 + O ( u + v + s ¯ ) 4 , g ˜ u , v , s ¯ = b 23 1 + η 2 + 1 + b 11 b 13 b 12 1 + η 2 N 2 + b 24 1 + η 2 + 1 + b 11 b 14 b 12 1 + η 2 N P + b 25 1 + η 2 + 1 + b 11 b 15 b 12 1 + η 2 P 2 + a 1 1 + b 11 b 12 1 + η 2 + d 1 1 + η 2 N 3 + a 2 1 + b 11 b 12 1 + η 2 + d 2 1 + η 2 N 2 P + a 3 1 + b 11 b 12 1 + η 2 + d 3 1 + η 2 N P 2 + a 4 1 + b 11 b 12 1 + η 2 + d 4 1 + η 2 P 3 + c 1 1 + η 2 s ¯ N + c 2 1 + η 2 s ¯ P + c 3 1 + η 2 s ¯ N P + c 4 1 + η 2 s ¯ N 2 + c 5 1 + η 2 r ¯ P 2 + O ( u + v + s ¯ ) 4 ,
where N = b 12 u + v and P = 1 b 11 u + η 2 b 11 v .
Thus, the approximation of the center manifold W c 0 , 0 , 0 of (19) within the neighborhood of s ¯ = 0 evaluated at the origin can be expressed as
W c 0 , 0 , 0 = { u , v , s 1 R 3 = M 3 s 2 + M 2 s u + m 1 u 2 + ( O u , s 1 ) 4 } ,
where
M 1 = b 12 2 b 23 1 + η 2 + 1 + b 11 b 13 b 12 1 + η 2 1 η 2 1 + b 11 b 12 b 24 1 + η 2 + 1 + b 11 b 14 b 12 1 + η 2 1 η 2 + 1 + b 11 2 b 25 1 + η 2 + 1 + b 11 b 15 b 12 1 + η 2 1 η 2 , M 2 = c 2 1 + b 11 c 1 b 12 η 2 2 1 ,   M 3 = 0 .
Consequently, the restricted map to center manifold W c 0 , 0 , 0 is expressed as follows:
F : u u + t 1 u 2 + t 2 u s + t 3 u 2 s + t 4 u s 2 + t 5 u 3 + ( O u , s 1 ) 4
where
t 1 = b 12 2 b 23 1 + η 2 + b 13 b 11 η 2 b 12 1 + η 2 1 + b 11 2 b 25 1 + η 2 + b 15 b 11 η 2 b 12 1 + η 2 + 1 + b 11 b 12 b 24 1 + η 2 + b 14 b 11 η 2 b 12 1 + η 2 , t 2 = c 2 1 + b 11 η 2 + 1 c 1 b 12 η 2 + 1 , t 3 = 2 η 2 b 11 b 13 b 12 η 2 + 1 b 23 η 2 + 1 b 12 2 M 2 + c 3 b 12 1 + b 11 η 2 + 1 c 4 b 12 2 η 2 + 1 + η 2 b 11 b 14 b 12 η 2 + 1 b 24 η 2 + 1 b 12 η 2 b 11 M 2 c 2 η 2 b 11 M 1 η 2 + 1 + b 11 η 2 b 14 b 12 η 2 + 1 + b 24 η 2 + 1 b 12 M 2 1 + b 11 c 5 1 + b 11 2 η 2 + 1 + 2   b 11 η 2 b 15 b 12 η 2 + 1 + b 25 η 2 + 1 1 + b 11 η 2 b 11 M 2 c 1 b 12 M 1 η 2 + 1 , t 4 = 2 η 2 b 11 b 13 b 12 η 2 + 1 b 23 η 2 + 1 b 12 2 M 3 c 2 η 2 b 11 M 2 η 2 + 1 + η 2 b 11 b 14 b 12 η 2 + 1 b 24 η 2 + 1 b 12 η 2 b 11 M 3 c 1 b 12 M 2 η 2 + 1 + 2   b 11 η 2 b 15 b 12 η 2 + 1 + b 25 η 2 + 1 1 + b 11 η 2 b 11 M 3 + b 11 η 2 b 14 b 12 η 2 + 1 + b 24 η 2 + 1 b 12 M 3 1 + b 11 t 5 = η 2 b 11 b 14 b 12 η 2 + 1 b 24 η 2 + 1 b 12 M 1 η 2 b 11 + b 11 b 11 η 2 b 14 b 12 η 2 + 1 + b 24 η 2 + 1 b 12 M 1 1 + b 11 + b 11 η 2 a 4 b 12 η 2 + 1 + d 4 η 2 + 1 1 + b 11 3 + 2   b 11 η 2 b 15 b 12 η 2 + 1 + b 25 η 2 + 1 1 + b 11 η 2 b 11 M 1 + 2   η 2 b 11 b 13 b 12 η 2 + 1 b 23 η 2 + 1 b 12 2 M 1 + η 2 b 11 a 3 b 12 η 2 + 1 d 3 η 2 + 1 b 12 1 + b 11 2 + η 2 b 11 a 1 b 12 η 2 + 1 d 1 η 2 + 1 b 12 3 + b 11 η 2 a 2 b 12 η 2 + 1 + d 2 η 2 + 1 b 12 2 1 + b 11 .

References

  1. Chen, B.S.; Chen, J.J. Bifurcation and chaotic behavior of a discrete singular biological economic system. Appl. Math. Comput. 2012, 219, 2371–2386. [Google Scholar] [CrossRef]
  2. Ghaziani, R.K.; Govaerts, W.; Sonck, C. Resonance and bifurcation in a discrete-time predator-prey system with Holling functional response. Nonlinear Anal. RWA 2012, 13, 1451–1465. [Google Scholar] [CrossRef] [Green Version]
  3. Jana, D. Chaotic dynamics of a discrete predator-prey system with prey refuge. Appl. Math. Comput. 2013, 224, 848–865. [Google Scholar] [CrossRef]
  4. Misra, O.P.; Sinha, P.; Singh, C. Stability and bifurcation analysis of a prey-predator model with age based predation. Appl. Math. Model. 2013, 37, 6519–6529. [Google Scholar] [CrossRef]
  5. Zhang, G.D.; Shen, Y.; Chen, B.S. Bifurcation analysis in a discrete differential-algebraic predator-prey system. Appl. Math. Model. 2014, 38, 4835–4848. [Google Scholar] [CrossRef]
  6. Hu, D.P.; Cao, H.J. Bifurcation and chaos in a discrete-time predator-prey system of Holling and Leslie type. Commun. Nonlinear Sci. Numer. Simulat. 2015, 22, 702–715. [Google Scholar] [CrossRef]
  7. Wang, C.; Li, X.Y. Further investigations into the stability and bifurcation of a discrete predator-prey model. J. Math. Anal. Appl. 2015, 422, 920–939. [Google Scholar] [CrossRef]
  8. Murdoch, W.; Briggs, C.; Nisbet, R. Consumer-Resource Dynamics; Princeton University Press: New York, NY, USA, 2003. [Google Scholar]
  9. Holling, C.S. The functional response of predator to prey density and its role in mimicry and population regulation. Mem. Entomol. Soc. Can. 1965, 97, 5–60. [Google Scholar] [CrossRef] [Green Version]
  10. Rosenzweig, M.L.; MacArthur, R.H. Graphical representation and stability conditions of predator-prey interactions. Am. Nat. 1963, 97, 209–223. [Google Scholar] [CrossRef]
  11. Freedman, H.I.; Mathsen, R.M. Persistence in predator-prey systems with ratio-dependent predator influence. Bull. Math. Biol. 1993, 55, 817–827. [Google Scholar] [CrossRef]
  12. Hastings, A. Multiple limit cycles in predator-prey models. J. Math. Biol. 1981, 11, 51–63. [Google Scholar] [CrossRef]
  13. Lindstrom, T. Qualitative analysis of a predator-prey systems with limit cycles. J. Math. Biol. 1993, 31, 541–561. [Google Scholar] [CrossRef]
  14. Jing, Z.J.; Yang, J. Bifurcation and chaos in discrete-time predator-prey system. Chaos Solitons Fractals 2006, 27, 259–277. [Google Scholar] [CrossRef]
  15. Liu, X.; Xiao, D. Complex dynamic behaviors of a discrete-time predator-prey system. Chaos Solitons Fractals 2007, 32, 80–94. [Google Scholar] [CrossRef]
  16. Lopez-Ruiz, R.; Fournier-Prunaret, D. Indirect Allee effect, bistability and chaotic oscillations in a predator-prey discrete model of logistic type. Chaos Solitons Fractals 2005, 24, 85–101. [Google Scholar] [CrossRef]
  17. Neubert, M.G.; Kot, M. The subcritical collapse of predator populations in discrete-time predator-prey models. Math. Biosci. 1992, 110, 45–66. [Google Scholar] [CrossRef]
  18. Shabbir, M.S.; Din, Q.; Ahmad, K.; Tassaddiq, A.; Soori, A.H.; Khan, M.A. Stability, bifurcation and chaos control of a novel discrete-time model involving Allee effect and cannibalism. Adv. Differ. Equ. 2020, 2020, 379. [Google Scholar] [CrossRef]
  19. Din, Q.; Shabbir, M.S.; Khan, M.A.; Ahmad, K. Bifurcation analysis and chaos control for a plant-herbivore model with weak predator functional response. J. Biol. Dyn. 2019, 13, 481–501. [Google Scholar] [CrossRef] [Green Version]
  20. Shabbir, M.S.; Din, Q.; Safeer, M.; Khan, M.A.; Ahmad, K.A. A dynamically consistent nonstandard finite difference scheme for a predator-prey model. Adv. Differ. Equ. 2019, 2019, 381. [Google Scholar] [CrossRef] [Green Version]
  21. Samaddar, S.; Dhar, M.; Bhattacharya, P. Effect of fear on prey-predator dynamics: Exploring the role of prey refuge and additional food. Chaos Interdiscip. J. Nonlinear Sci. 2020, 30, 63129. [Google Scholar] [CrossRef]
  22. Anacleto, M.; Vidal, C. Dynamics of a delayed predator-prey model with Allee effect and Holling type II functional response. Math. Methods Appl. Sci. 2020, 43, 5708–5728. [Google Scholar] [CrossRef]
  23. Tang, B. Dynamics for a fractional-order predator-prey model with group defense. Sci. Rep. 2020, 10, 4906. [Google Scholar] [CrossRef]
  24. Sarwardi, S.; Haque, M.M.; Hossain, S. Analysis of Bogdanov-Takens bifurcations in a spatiotemporal harvested-predator and prey system with Beddington–DeAngelis-type response function. Nonlinear Dyn. 2020, 100, 1755–1778. [Google Scholar] [CrossRef]
  25. Din, Q. Complexity and chaos control in a discrete-time prey-predator model. Commun. Nonlinear. Sci. Numer. Simulat. 2017, 49, 113–134. [Google Scholar] [CrossRef]
  26. Shabbir, M.S.; Din, Q.; Alabdan, R.; Tassaddiq, A.; Ahmad, K. Dynamical complexity in a class of novel discrete-time predator-prey interaction with cannibalism. IEEE Access 2020, 8, 100226–100240. [Google Scholar] [CrossRef]
  27. Selvam, A.G.; Vianny, D.A.; Jacob, S.B. Dynamical behaviour of discrete time prey-predator model with Holling type III functional response. Cikitusi J. Multidiscip. Res. 2019, 6, 75–81. [Google Scholar]
  28. Jiangang, Z.; Tian, D.; Yandong, C.; Shuang, Q.; Wenju, D.; Hongwei, L. Stability and bifurcation analysis of a discrete predator-prey model with Holling type III functional response. J. Nonlinear Sci. Appl. 2016, 2016, 6228–6243. [Google Scholar]
  29. Elettreby, M.F.; Nabil, T.; Khawagi, A. Stability and bifurcation analysis of a discrete predator-prey model with mixed Holling interaction. Comput. Modeling Eng. Sci. 2020, 122, 907–921. [Google Scholar] [CrossRef]
  30. Murray, J.D. Mathematical Biology, 2nd ed.; Springer: Berlin, Germany, 1993. [Google Scholar]
  31. He, Z.; Lai, X. Bifurcation and chaotic behavior of a discrete-time predator-prey system. Nonlinear Anal. Real World Appl. 2011, 12, 403–417. [Google Scholar] [CrossRef]
  32. Saber, N. Elaydi, Discrete Chaos; Chapman & Hall/CRC: Boca Raton, FL, USA, 2007. [Google Scholar]
  33. Chen, Y.U. Controlling and anti-controlling Hopf bifurcations in discrete maps using polynomial functions. Chaos Solitons Fractals 2005, 26, 1231–1248. [Google Scholar] [CrossRef]
  34. ELabbasy, E.M.; Agiza, H.N.; Metwally, H.E.L.; Elsadany, A.A. Bifurcation analysis, chaos and control in the Burgers mapping. Int. J. Nonlinear Sci. 2007, 4, 171–185. [Google Scholar]
  35. Tassaddiq, A.; Shabbir, M.S.; Din, Q.; Ahmad, K. A ratio-dependent nonlinear predator-prey model with certain dynamical results. IEEE Access 2020, 8, 195074–195088. [Google Scholar] [CrossRef]
  36. Din, Q.; Saleem, N.; Shabbir, M.S. A class of discrete predator-prey interaction with bifurcation analysis and chaos control. Math. Model. Nat. Phenom. 2020, 15, 1–27. [Google Scholar] [CrossRef]
  37. Din, Q.; Shabbir, M.S. A Cubic autocatalator chemical reaction model with limit cycle analysis and consistency preserving discretization. MATCH Commun. Math. Comput. Chem. 2022, 87, 441–462. [Google Scholar] [CrossRef]
  38. Cartwright, J.H.E. Nonlinear stiffness Lyapunov exponents and attractor dimension. Phys. Lett. A 1999, 264, 298–304. [Google Scholar] [CrossRef] [Green Version]
  39. Kaplan, J.L.; Yorke, J.A. Preturbulence: A regime observed in a fluid flow model of Lorenz. Commun. Math. Phys. 1979, 67, 93–108. [Google Scholar] [CrossRef]
Figure 1. a   s = 1.34 , b   s = 1.3811 ,   c   s = 1.3874082082631611 , d   s = 1.43 , e   s = 1.607 , f   s = 1.625 , g   s = 1.66 ,   h   s = 1.68 , i   s = 1.735 , j   s = 1.83 , k   s = 1.89 , l   s = 1.99 . (a)–(l) Phase portraits of system (4) for different values of s 1 , 2 with b = 1.7 , a = 2.26 , r = 2.4 , h = 1.8 , and   k = 1.4 under the initial conditions x 0 = 1.0219 , y 0 = 0.567721 .
Figure 1. a   s = 1.34 , b   s = 1.3811 ,   c   s = 1.3874082082631611 , d   s = 1.43 , e   s = 1.607 , f   s = 1.625 , g   s = 1.66 ,   h   s = 1.68 , i   s = 1.735 , j   s = 1.83 , k   s = 1.89 , l   s = 1.99 . (a)–(l) Phase portraits of system (4) for different values of s 1 , 2 with b = 1.7 , a = 2.26 , r = 2.4 , h = 1.8 , and   k = 1.4 under the initial conditions x 0 = 1.0219 , y 0 = 0.567721 .
Fractalfract 06 00031 g001aFractalfract 06 00031 g001b
Figure 2. Bifurcation diagrams and maximum Lypunov exponents (MLE) for system (4) with parametric values b = 1.7 ,   a = 2.26 ,   r = 2.4 ,   h = 1.8 ,   k = 1.4 , and s 1 , 1.8 and initial conditions x 0 , y 0 = 1.0219 , 0.567721 : (a) bifurcation for x n , (b) bifurcation for y n , and (c) MLE.
Figure 2. Bifurcation diagrams and maximum Lypunov exponents (MLE) for system (4) with parametric values b = 1.7 ,   a = 2.26 ,   r = 2.4 ,   h = 1.8 ,   k = 1.4 , and s 1 , 1.8 and initial conditions x 0 , y 0 = 1.0219 , 0.567721 : (a) bifurcation for x n , (b) bifurcation for y n , and (c) MLE.
Fractalfract 06 00031 g002
Figure 3. Bifurcation diagrams and MLE for system (4) with parameters b = 3.3 ,   a = 0.8 ,   r = 1.3 ,   h = 2.7 ,   k = 1.8 , s 2 ,   3 , and x 0 , y 0 = 1.74693 , 0.647013 : (a) period-doubling bifurcation for x n , (b) period-doubling bifurcation for y n , and (c) MLE.
Figure 3. Bifurcation diagrams and MLE for system (4) with parameters b = 3.3 ,   a = 0.8 ,   r = 1.3 ,   h = 2.7 ,   k = 1.8 , s 2 ,   3 , and x 0 , y 0 = 1.74693 , 0.647013 : (a) period-doubling bifurcation for x n , (b) period-doubling bifurcation for y n , and (c) MLE.
Fractalfract 06 00031 g003
Figure 4. (a) Bifurcation diagrams of x n for controlled system (28), (b) bifurcation diagrams of y n for controlled system (28), (c) bifurcation diagrams of x n for controlled system (30), and (d) bifurcation diagrams of y n for controlled system (30).
Figure 4. (a) Bifurcation diagrams of x n for controlled system (28), (b) bifurcation diagrams of y n for controlled system (28), (c) bifurcation diagrams of x n for controlled system (30), and (d) bifurcation diagrams of y n for controlled system (30).
Fractalfract 06 00031 g004
Figure 5. (a) Local implication for the Neimark-Sacker bifurcation and (b) local implication for period-doubling bifurcation.
Figure 5. (a) Local implication for the Neimark-Sacker bifurcation and (b) local implication for period-doubling bifurcation.
Fractalfract 06 00031 g005
Table 1. Fractal dimension of model (4).
Table 1. Fractal dimension of model (4).
Values   of   s 1 st   Lyapunov   Exponents   F 1 2 nd   Lyapunov   Exponents   F 2 Fractal   Dimension   D 𝓁
2.85 0.08462596943938297 1.1368798813345231 1.0744370366903186
2.90 0.22741225178613755 1.2160641186798002 1.187006793714976
3.0 0.31895100399320747 1.0790255038850196 1.2955917194216706
3.1 0.22493177216760443 1.244489487648406 1.1807422034497348
3.2 0.4025124673527987 1.1944261491614452 1.3369923436751492
3.3 0.3894200849244259 1.2328810276916689 1.3158618521801246
3.4 0.47745428811163265 1.2528556268034083 1.3810928233844715
3.5 0.47582043180971084 1.2925566246939908 1.3681234715131803
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Tassaddiq, A.; Shabbir, M.S.; Din, Q.; Naaz, H. Discretization, Bifurcation, and Control for a Class of Predator-Prey Interactions. Fractal Fract. 2022, 6, 31. https://doi.org/10.3390/fractalfract6010031

AMA Style

Tassaddiq A, Shabbir MS, Din Q, Naaz H. Discretization, Bifurcation, and Control for a Class of Predator-Prey Interactions. Fractal and Fractional. 2022; 6(1):31. https://doi.org/10.3390/fractalfract6010031

Chicago/Turabian Style

Tassaddiq, Asifa, Muhammad Sajjad Shabbir, Qamar Din, and Humera Naaz. 2022. "Discretization, Bifurcation, and Control for a Class of Predator-Prey Interactions" Fractal and Fractional 6, no. 1: 31. https://doi.org/10.3390/fractalfract6010031

APA Style

Tassaddiq, A., Shabbir, M. S., Din, Q., & Naaz, H. (2022). Discretization, Bifurcation, and Control for a Class of Predator-Prey Interactions. Fractal and Fractional, 6(1), 31. https://doi.org/10.3390/fractalfract6010031

Article Metrics

Back to TopTop