In this paper, some integral inequalities of Hermite-Hadamard type for multiplicatively preinvex functions are established. Also, new inequalities involving multiplicative integrals are obtained for product and quotient of preinvex and multiplicatively preinvex functions.
1.
Introduction
The well-known Leslie-Gower model with a Holling Ⅱ functional response is expressed as follows:
where the biological interpretation of each parameter is provided in Table 1. It is noteworthy that the Leslie-Gower model incorporates a response function indicating the predator's carrying capacity, which is proportional to the prey population, a feature absent in the Lotka-Volterra model. This distinction has been extensively investigated in many references, such as [1,2,3,4,5].
Previous studies predominantly focused on direct predation effects. In 2016, Wang et al. [6] demonstrated through experiments that predator-induced fear (indirect effects) in prey leads to a reduction in prey birth rates. They introduced the fear factor F(k,y)=11+ky, where k represents the intensity of fear driving antipredator behavior. To further investigate the influence of fear on population dynamics, we enhance the system (1.1) by incorporating the fear factor F(k,y). The modified model is formulated as follows:
For additional significant findings on the impact of fear in predator-prey models, refer to [7,8,9].
There has been increasing recognition that traditional integer-order differential equations may not sufficiently capture the complexity of biological systems [10,11,12,13,14]. In reality, the behaviors of most organisms in nature are influenced by their historical context. Fractional-order derivatives, which extend the concept of integer-order differentiation, offer a more flexible and accurate framework for modeling memory and hereditary properties in population dynamics. Therefore, the authors of this paper aim to apply the Caputo fractional derivative to (1.2), thereby extending it into a fractional model. The Caputo fractional derivative of a function u(t) of order α∈(0,1] is given in [15] as follows
where Γ is the Gamma function. By replacing the integer-order derivative with the Caputo fractional derivative, the following model is obtained:
Specifically, when α=1, the fractional-order model (1.3) reduces to the integer-order model (1.2), demonstrating that the fractional-order model serves as a generalization of the integer-order model. Moreover, in fractional calculus, the rate of change at any given moment, expressed by the fractional-order derivative, depends on the population density over a specified time interval. This feature gives the fractional-order model (1.3) a distinct advantage in capturing memory effects within populations.
When studying populations with nonoverlapping generations or small sizes, mathematical models are often expressed in discrete terms. Although some variables may change continuously in real-world scenarios, the recording of these changes typically happens at specific time intervals during data collection. Thus, employing discrete systems to examine the dynamic behavior of biological populations is highly practical and significant. In references [16,17], the authors employed the piecewise constant approximation method to discretize continuous fractional predator-prey models and explored their dynamical properties. Singh and Sharma [18] examine a discrete prey-predator model with Holling type Ⅱ functional response and prey refuge, identifying bifurcations and controlling chaos through state feedback, pole placement, and hybrid techniques. Berkal and Almatrafi [19] used the exponential piecewise constant argument to discretize continuous fractional activator-inhibitor system. Their analysis includes stability assessments, investigations of Neimark-Sacker and period-doubling bifurcations, and numerical simulations that validate the theoretical findings on the system's dynamics. For a more comprehensive exploration of discrete model studies, readers are advised to consult references [20,21,22,23,24,25] and the related literature cited therein. However, current literature lacks studies on fractional discrete-time predator-prey Leslie-Gower systems that incorporate the fear effect in prey population. Using the same method as in references [16,17] to discretize model (1.3), we can obtain the following discrete model:
Here h>0 represents the time interval of production.
Furthermore, the key contributions and findings of this study are summarized as follows:
● The Caputo fractional derivative of order (0,1] is utilized to incorporate the memory effect into the dynamical behavior of the proposed model.
● The existence and stability of fixed points are investigated.
● Conditions for the occurrence and direction of period-doubling bifurcation and Neimark–Sacker bifurcation at the positive fixed point are established.
● State feedback and hybrid control strategies are employed to manage bifurcations and chaotic behavior in the model.
● To validate the accuracy of our theoretical findings, numerical examples for the fractional-order discrete-time Leslie-Gower model with a fear factor are provided.
The remainder of this paper is structured as follows: In Section 2, we investigate the existence and stability of the equilibrium points of model (1.4). Section 3 analytically demonstrates that model (1.4), under specific parametric conditions, undergoes period-doubling or Neimark-Sacker bifurcation. Section 4 explores the control of chaos toward an unstable equilibrium point using feedback control or hybrid control approaches. Section 5 presents a quantitative analysis of the dynamics of model (1.4) to validate our analytical findings. Finally, Section 6 offers brief conclusions.
2.
Existence and stability of the equilibrium points
This section explores the existence of equilibrium points in model (1.4) and assesses their stability by evaluating the eigenvalues of the Jacobian matrix at these points. The following definition and lemma are introduced to assist in this stability analysis of equilibrium points.
Definition 1. [26] Let λ1 and λ2 denote the two roots of the characteristic equation F(λ)=λ2+pλ+q=0 associated with the Jacobian matrix J(x,y). The equilibrium point (x,y) is termed
(1) sink if |λ1|<1 and |λ2|<1, and the sink is locally asymptotically stable;
(2) source if |λ1|>1 and |λ2|>1, and the source is locally unstable;
(3) saddle if |λ1|>1 and |λ2|<1 (or |λ1|<1 and |λ2|>1);
(4) non-hyperbolic if either |λ1|=1 or |λ2|=1.
Lemma 1. [26] Let F(1)>0 in F(λ)=λ2−Mλ+N, where λ1,λ2 are the two roots of F(λ)=0. Then, the following results hold true:
(1) |λ1|<1 and |λ2|<1 if, and only if, F(−1)>0 and N<1;
(2) |λ1|<1 and |λ2|>1 (or |λ1|>1 and |λ2|<1) if, and only if, F(−1)<0;
(3) |λ1|>1 and |λ2|>1 if, and only if, F(−1)>0 and N>1;
(4) λ1=−1 and |λ2|≠1 if, and only if, F(−1)=0 and M≠0,−2;
(5) λ1 and λ2 are conjugate complex and |λ1,2|=1 if, and only if, M2<4N and N=1.
2.1. Existence of equilibrium points
It is evident that the equilibrium points of model (1.4) satisfy the following equations:
This algebraic system is satisfied if x=K and y=0, indicating that the model (1.4) has a boundary equilibrium point E0(K,0) for all model parameters. To find the interior equilibrium point E1, we will solve the following system simultaneously:
From the second equation, we find y=x. Substituting y=x into the first equation yields:
Given that all model parameters are positive, we obtain the following results through direct calculations.
Theorem 1. The model (1.4) always has a boundary equilibrium point E0(K,0) and a positive equilibrium point E1=(x∗,y∗), where
2.2. Stability of equilibrium points
The Jacobian matrix of model (1.4) evaluated at any point (x,y) is as follows:
Theorem 2. The boundary equilibrium point E0(K,0) exhibits the following behaviors:
(1) It is a source point if r>2Γ(α+1)hα;
(2) It is a saddle point if 0<r<2Γ(α+1)hα;
(3) It is a non-hyperbolic point if r=2Γ(α+1)hα.
Proof. The Jacobian matrix (2.1) evaluated at the boundary equilibrium point E0(K,0) is given by:
The eigenvalues of J(E0) are λ1=1−hαrΓ(α+1) and λ2=1+hαsΓ(α+1). Clearly, |λ2|>1 and:
By applying Definition 1, the proof is complete. □
Next, we analyze the local dynamics of model (1.4) at the positive equilibrium point E1(x∗,y∗). The Jacobian matrix (2.1) evaluated at E1(x∗,y∗) is expressed as:
where
The characteristic equation of J(E1) is given by:
where
Let F(λ)=λ2−Mλ+N, then:
As a result, we establish the following theorem.
Theorem 3. Let E1 be the unique positive equilibrium point of model (1.4). Then:
(1) E1 is a sink point if |M|<1+N<2;
(2) E1 is a source point if |M|<|1+N| and |N|>1;
(3) E1 is a saddle point if M2>4N and |M|>|1+N|;
(4) E1 is non-hyperbolic if |M|=|1+N| or N=1 and |M|<2.
Proof. (1) According to Lemma 1, E1 is a sink point if and only if F(1)>0, F(−1)>0, and N<1. This condition is satisfied when |M|<1+N<2. Similarly, the proofs for Theorem 3 (2)–(4) follow straightforwardly from the definitions and properties of M and N. □
3.
Bifurcation analysis
In this section, we will analyze the period-doubling and Neimark-Sacker bifurcation behaviors of model (1.4) at the positive equilibrium point E1(x∗,y∗). The reason for not analyzing the boundary equilibrium point is that, at E0(K,0), the predators have become extinct, leaving only the prey.
3.1. Period-doubling bifurcation
Assume M2>4N and
and we can ascertain that the characteristic equation (2.3) satisfies
Therefore, the eigenvalues of the characteristic equation (2.3) are
Let
The condition |λ2|≠1 implies that s1≠s2,s3.
Based on the above analysis, we can conclude that when the parameters vary within the set
model (1.4) will undergo period-doubling bifurcation at E1(x∗,y∗).
Select parameters (h,α,r,k,K,b,a,s)∈P.D and consider s∗ as a small perturbation of s, i.e., s∗=s−s1, where |s∗|≪1. After this perturbation, model (1.4) can be represented as follows:
Let un=xn−x∗ and vn=yn−y∗, transforming the equilibrium point E1(x∗,y∗) into the origin O(0,0). The model (3.2) can thus be expressed as:
The Taylor expansion of model (3.3) around (un,vn)=(0,0) yields the following form:
where
and
Define
then
Using the transformation:
model (3.4) is transformed into the following form:
where
and
Next, we apply the center manifold theorem [26] to analyze the dynamics around the equilibrium point (˜un,˜vn)=(0,0) at s∗=0. According to the theorem, the model (3.5) has a center manifold, which can be represented as:
Assume that
Then, the center manifold must satisfy
By comparing the coefficients, it can be obtained that
Thus, the model (3.5), when restricted to the center manifold Wc(0,0,0), is given by:
where
In order for Eq (3.6) to undergo period-doubling bifurcation, it is necessary that the following two quantities possess nonzero values.
Summarize the above analysis into the following theorem.
Theorem 4. If β1β2≠0, then model (1.4) undergoes a period-doubling bifurcation at the positive equilibrium point E1(x∗,y∗) when parameters vary in a small neighborhood of P.D. Additionally, when β2>0 (respectively, β2<0), model (1.4) bifurcates from the equilibrium point E1(x∗,y∗) to a stable (respectively, unstable) 2-periodic orbit.
3.2. Nerimark-Sacker bifurcation
Assume M2<4N and
and we can determine that the eigenvalues of the characteristic equation (2.3) satisfy
Namely, Eq (2.3) has two complex conjugate roots with unit modulus. It is clear from the above discussion that
In addition, it is crucial that when s=s4, λθ1,2(s4)≠1(θ=1,2,3,4), which is equivalent to M(s4)≠−2,−1,0,2. Since M2<4N, we deduce M(s4)≠−2,2. Additionally, we necessitate that M(s4)≠0,−1, which leads to
Based on the preceding analysis, we conclude that when the parameters vary within a small neighborhood of the set
model (1.4) undergoes a Neimark-Sacker bifurcation at E1(x∗,y∗).
Next, assuming (h,α,r,k,K,b,a,s)∈N.S. and |s∗|≪1 represents a small perturbation of s4, model (1.4) can be described as follows:
Let un=xn−x∗ and vn=yn−y∗, then we have
where
and c13, c14, c15, c16, c17, c18, c19, c23, c24, c25, c26, c27, c28, c29 are given in (2.3) by substituting s1 for s4+s∗.
In order to obtain the normal form of model (3.9) at s∗=0, we use the following transformation:
with
Using this transformation, model (3.9) will transform as follows:
where
and
Next, a nonzero real number is defined as follows:
where
Through some complicated calculations, we get
and
Based on the aforementioned calculations, we establish the theorem regarding the existence and direction of Neimark-Sacker bifurcation.
Theorem 5. If (h,α,r,k,K,b,a,s)∈N.S. and Ω≠0, then model (1.4) undergoes a Neimark-Sacker bifurcation at the equilibrium point E1(x∗,y∗) when the parameter s varies in the vicinity of s4. Moreover, if Ω<0 (respectively, Ω>0), then an attracting (respectively, repelling) closed invariant curve bifurcates from the equilibrium point for s∗>0 (respectively, s∗<0).
4.
Chaos control
Chaos often has detrimental effects on biological systems, disrupting the ecological balance of populations and directly influencing long-term population growth projections. Implementing effective control policies not only safeguards the size of ecological populations but also establishes a strong foundation for the sustainable exploitation of ecological resources [27,28]. This section explores two control methods aimed at effectively managing the chaos generated by model (1.4).
4.1. State feedback control
In this subsection, the state feedback control method [22] will be employed to regulate the chaos exhibited by model (1.4). To achieve this, we introduce the following controlled model.
which corresponds to model (1.4). The feedback controlling force is defined as
where (x∗,y∗) represents the positive equilibrium point of the model (1.4), and p1, p2 stand for the feedback gains. The Jacobian matrix of the controlled model (4.1) evaluated at the positive equilibrium point E1(x∗,y∗) is given by
where the variables A, a11, and a12 are defined in Eq (2.2). The corresponding characteristic equation of the Jacobian matrix J1(x∗,y∗) is
Let λ1 and λ2 be the roots of the Eq (4.3), then
The lines of marginal stability l1, l2, and l3 are derived by solving λ1λ2=1, λ1=1 and λ2=±1, respectively. These conditions ensure that |λ1,2|=1. Then, we derive the marginal stability lines as follows:
Therefore, l1, l2, and l3 in the (p1,p2)-plane form a triangular region which leads to |λ1,2|<1.
Theorem 6. If p1 and p2 lie within a triangular region bounded by the lines l1, l2, and l3, it can be concluded that the model (4.1) is stable.
4.2. Hybrid control
Next, we apply the hybrid control approach proposed by [23] to control chaos. The controlled model of (1.4) with the hybrid control approach is depicted below:
where 0<ρ<1, and the controlled strategy in (4.9) combines feedback control and parameter perturbation. By appropriately choosing the controlled parameter ρ, the chaotic behaviors of the equilibrium point (x∗,y∗) of the controlled model (4.9) can be accelerated (delayed) or even entirely eliminated. The Jacobian matrix of the controlled model (4.9), evaluated at the positive equilibrium point (x∗,y∗), is given by
where the variables A, a11, and a12 are defined in Eq (2.2). Then, the positive equilibrium point (x∗,y∗) of the controlled model (4.9) is locally asymptotically stable if the roots of the characteristic polynomial of (4.10) lie within the open unit disk. According to the Jury condition, the equilibrium point of the model remains stable if, and only if, the following conditions are met:
Theorem 7. If |2+Aρa11−Aρs|<1+(1+Aρa11)(1−Aρs)−(Aρa12)(Aρs)<2 can hold, it can be concluded that the model (4.9) is stable.
5.
Numerical experiments
In this section, by exploring specific cases of model (1.4), we confirm the theoretical analysis above and discover new intriguing complex dynamic behaviors. Additionally, we validate the effectiveness of linear feedback techniques and hybrid control strategies for chaos control through numerical simulations.
5.1. Period-doubling bifurcation at positive equilibrium point
Take
After straightforward calculations, we determine the positive equilibrium point E1=(1.7498,1.7498) and the critical bifurcation value s1=2.8472. The Jacobian matrix of model (1.4) evaluated at E1 with s=s1 is given by
whose characteristic equation is
The roots of (5.2) are λ1=−1 and λ2=0.4768. Moreover, we have β1=−0.7959≠0 and β2=1.0080>0. According to Theorem 4, model (1.4) undergoes period-doubling bifurcation at E1 as s passes through s1. This behavior, verified through corresponding bifurcation diagrams shown in Figure 1, utilizes initial conditions (x0,y0)=(1.74,1.74) and varies s in the range [2.7,3.7].
5.2. Neimark-Sacker bifurcation at positive equilibrium point
Take
After straightforward calculations, we determine the positive equilibrium point E1=(2.1746,2.1746) and the critical bifurcation value s4=1.1872. The Jacobian matrix of model (1.4) evaluated at E1 with s=s4 is given by
whose characteristic equation is
The roots of (5.4) are λ1,2=0.02930±0.9996i. Moreover, we have d=0.3912≠0 and Ω=−1392.0464<0. According to Theorem 5, model (1.4) undergoes Neimark-Sacker bifurcation at E1 as s passes through s4. This behavior, verified through corresponding bifurcation diagrams shown in Figure 2, utilizes initial conditions (x0,y0)=(2.2,2.2) and varies s in the range [1.1,1.6].
5.3. Chaos control
We chose the parameter values as follows:
By Figure 3, we can get that the variables xn and yn in the model (1.4) are in a chaotic state.
5.3.1. State feedback control
We use the same parameter values as in (5.5). In Figure 4, the triangular region defined by Theorem 6 bounds the parameters p1 and p2. Inside this region, the chaotic behavior produced by model (1.4) is effectively managed, resulting in asymptotic convergence toward the equilibrium point E1=(2.1746,2.1746).
In this case, with feedback gains set to p1=0.57 and p2=−0.38 and the controller activated at the 3000th iteration of the model, Figure 5 demonstrates the control effect. A chaotic trajectory is successfully stabilized at the equilibrium point E1=(2.1746,2.1746). This indicates that the feedback control approach is effective in mitigating bifurcation and chaos.
5.3.2. Hybrid control
Finally, using the same parameter values as in (5.5), the Jacobian matrix of the controlled model (4.9), evaluated at E1, is given by:
The characteristic polynomial of (5.6) is given by
The roots of (5.7) lie within the open unit disk if, and only if, 0<ρ<0.881058. Additionally, the plots for xn and yn of the controlled model (4.9) are shown in Figure 6 with ρ=0.87. From Figure 6, it is clear that the positive equilibrium point E1 is stable. Therefore, it can be concluded that employing the hybrid control approach is effective in mitigating bifurcation and chaos.
6.
Conclusions
In this study, we proposed a fractional Leslie-Gower model with a Holling type Ⅱ functional response and antipredator behavior. By employing the piecewise constant approximation method, we derived the discrete model (1.4) and analyzed its dynamical behavior, including the existence and stability of equilibrium points and the possibility of local bifurcations. The following conclusions can be drawn from our research:
(1) The discrete model (1.4) has two equilibrium points: E0(K,0) and E1(x∗,y∗). The only positive coexistence equilibrium point, E1(x∗,y∗), reflects the coexistence of predators and prey.
(2) Our theoretical analysis and numerical simulations of the positive equilibrium point E1(x∗,y∗) indicate that the model undergoes period-doubling bifurcation and Neimark-Sacker bifurcation under specific parameter conditions. Figures 1 and 2 illustrate how these bifurcations can lead to chaotic behavior at E1.
(3) By applying state feedback control and hybrid control methods, we effectively managed the chaotic behavior generated by the model (1.4), as shown in Figures 3–6. These interventions mitigated the adverse effects of chaos and bifurcations, consequently enhancing ecosystem resilience.
This study advances our understanding of the complex dynamics in ecological models influenced by fear effects and provides practical techniques for controlling chaotic behavior in such models. Future work could explore different discretization methods for the model, and new parameters could be chosen to study the influence of various ecological effects on population dynamics.
Author contributions
Yao Shi: Formal analysis, validation, writing-original draft & editing, visualization, methodology; Zhenyu Wang: Conceptualization, investigation, writing-review & editing, software, supervision. All authors have read and approved the final version of the manuscript for publication.
Acknowledgments
This work was supported by the Innovation Foundation of Hebei University of Engineering (Grant Nos. SJ2401002097), the National Natural Science Foundation of China (Grant Nos. 12401519), and the Central Guidance on Local Science and Technology Development Fund of Hebei Province (Grant Nos. 246Z1825G).
Conflict of interest
The authors declare no conflict of interest.