Processing math: 100%

On Lennard-Jones systems with finite range interactions and their asymptotic analysis

  • The aim of this work is to provide further insight into the qualitative behavior of mechanical systems that are well described by Lennard-Jones type interactions on an atomistic scale. By means of $Γ$-convergence techniques, we study the continuum limit of one-dimensional chains of atoms with finite range interactions of Lennard-Jones type, including the classical Lennard-Jones potentials. So far, explicit formula for the continuum limit were only available for the case of nearest and next-to-nearest neighbour interactions. In this work, we provide an explicit expression for the continuum limit in the case of finite range interactions. The obtained homogenization formula is given by the convexification of a Cauchy-Born energy density.

    Furthermore, we study rescaled energies in which bulk and surface contributions scale in the same way. The related discrete-to-continuum limit yields a rigorous derivation of a one-dimensional version of Griffith' fracture energy and thus generalizes earlier derivations for nearest and next-to-nearest neighbors to the case of finite range interactions.

    A crucial ingredient to our proofs is a novel decomposition of the energy that allows for refined estimates.

    Citation: Mathias Schäffner, Anja Schlömerkemper. 2018: On Lennard-Jones systems with finite range interactions and their asymptotic analysis, Networks and Heterogeneous Media, 13(1): 95-118. doi: 10.3934/nhm.2018005

    Related Papers:

    [1] San-Xing Wu, Xin-You Meng . Dynamics of a delayed predator-prey system with fear effect, herd behavior and disease in the susceptible prey. AIMS Mathematics, 2021, 6(4): 3654-3685. doi: 10.3934/math.2021218
    [2] Binfeng Xie, Na Zhang . Influence of fear effect on a Holling type III prey-predator system with the prey refuge. AIMS Mathematics, 2022, 7(2): 1811-1830. doi: 10.3934/math.2022104
    [3] Fatao Wang, Ruizhi Yang, Yining Xie, Jing Zhao . Hopf bifurcation in a delayed reaction diffusion predator-prey model with weak Allee effect on prey and fear effect on predator. AIMS Mathematics, 2023, 8(8): 17719-17743. doi: 10.3934/math.2023905
    [4] Jie Liu, Qinglong Wang, Xuyang Cao, Ting Yu . Bifurcation and optimal harvesting analysis of a discrete-time predator–prey model with fear and prey refuge effects. AIMS Mathematics, 2024, 9(10): 26283-26306. doi: 10.3934/math.20241281
    [5] A. Q. Khan, Ibraheem M. Alsulami . Complicate dynamical analysis of a discrete predator-prey model with a prey refuge. AIMS Mathematics, 2023, 8(7): 15035-15057. doi: 10.3934/math.2023768
    [6] Xiaoyan Zhao, Liangru Yu, Xue-Zhi Li . Dynamics analysis of a predator-prey model incorporating fear effect in prey species. AIMS Mathematics, 2025, 10(5): 12464-12492. doi: 10.3934/math.2025563
    [7] Xiaoming Su, Jiahui Wang, Adiya Bao . Stability analysis and chaos control in a discrete predator-prey system with Allee effect, fear effect, and refuge. AIMS Mathematics, 2024, 9(5): 13462-13491. doi: 10.3934/math.2024656
    [8] Kottakkaran Sooppy Nisar, G Ranjith Kumar, K Ramesh . The study on the complex nature of a predator-prey model with fractional-order derivatives incorporating refuge and nonlinear prey harvesting. AIMS Mathematics, 2024, 9(5): 13492-13507. doi: 10.3934/math.2024657
    [9] Weili Kong, Yuanfu Shao . Bifurcations of a Leslie-Gower predator-prey model with fear, strong Allee effect and hunting cooperation. AIMS Mathematics, 2024, 9(11): 31607-31635. doi: 10.3934/math.20241520
    [10] Yaping Wang, Yuanfu Shao, Chuanfu Chai . Dynamics of a predator-prey model with fear effects and gestation delays. AIMS Mathematics, 2023, 8(3): 7535-7559. doi: 10.3934/math.2023378
  • The aim of this work is to provide further insight into the qualitative behavior of mechanical systems that are well described by Lennard-Jones type interactions on an atomistic scale. By means of $Γ$-convergence techniques, we study the continuum limit of one-dimensional chains of atoms with finite range interactions of Lennard-Jones type, including the classical Lennard-Jones potentials. So far, explicit formula for the continuum limit were only available for the case of nearest and next-to-nearest neighbour interactions. In this work, we provide an explicit expression for the continuum limit in the case of finite range interactions. The obtained homogenization formula is given by the convexification of a Cauchy-Born energy density.

    Furthermore, we study rescaled energies in which bulk and surface contributions scale in the same way. The related discrete-to-continuum limit yields a rigorous derivation of a one-dimensional version of Griffith' fracture energy and thus generalizes earlier derivations for nearest and next-to-nearest neighbors to the case of finite range interactions.

    A crucial ingredient to our proofs is a novel decomposition of the energy that allows for refined estimates.



    In population ecology, understanding how predators and primary producers influence nutrient flow relative to each other is important. Ecosystem interactions and predator-prey relationships are governed by predation and the delivery of resource processes. The identification of ecological factors that can alter or control dynamic behavior requires theoretical and experimental research. One way to study these questions is by means of experimental control, and another useful way is via mathematical modeling as well as computer simulations. Over decades of theoretical ecology and biomathematics development, mathematical modeling has become an indispensable tool for scientists in related fields to study ecosystems. Since Lotka [1] and Volterra [2], as cornerstones of theoretical ecology, published the first study of predator-prey dynamics, any species in nature can be a predator or prey, and due to its prevalence, it has become one of the most popular topics for researchers to study [3,4,5]. Besides, because biological resources are renewable and have the most unique development mechanisms, the over-utilization of biological resources and the destruction of the environment by humans will directly affect the balance of the ecosystem. Maintaining ecological balance and meeting humans material needs have attracted the most attention from researchers focused on the scientific management of renewable resource development [6,7,8].

    Shelter serves as a defense strategy. It refers broadly to a series of behaviors by prey to avoid predators in order to increase their survival rate. The concept of sanctuary was first developed by Maynard-Smith [9] and Gause et al. [10], and its popularity has been very high, garnering widespread attention from many scholars [11,12,13,14,15]. Sih et al.[16] investigated the effects of prey refuge in a three-species model and concluded that the system's stability is strongly related to the refuge. Also, similar findings can be displayed in [17,18,19,20,21,22]. The two modes of refuge analyzed by Gonzalez-Olivares et al. [17] have diverse stability domains in terms of the parameter space. Qi et al.[21] ensure the stability of the system by varying the strength of the refuge.

    Through reviewing a large amount of literature, we begin to consider [23,24] as a basis for the two prey and one predator species that will be modeled in this article. We assume that at a certain time $ t $, the populations of the two prey and one predator are $ x_1(t) $, $ x_2(t) $, and $ y(t) $, respectively. Based on the above, we construct the following model:

    $ {dx1dt=r1x1(1x1K1)a1x1x2c1(1m1)x1yq1E1x1,dx2dt=r2x2(1x2K2)a2x1x2c2(1m2)x2yq2E2x2,dydt=e1(1m1)x1y+e2(1m2)x2ydyq3E3y.
    $
    (1.1)

    The significance of the full parameters is annotated in Table 1.

    Table 1.  Biological meaning of parameters.
    Parameters Biological meaning
    $ r_1, r_2 $ Growth rates of prey $ x_1 $ and prey $ x_2 $
    $ K_1, K_2 $ Carrying capacity of prey $ x_1 $ and prey $ x_2 $
    $ a_1, a_2 $ Interspecific competition between prey $ x_1 $ and prey $ x_2 $
    $ c_1, c_2 $ Predation coefficients for prey $ x_1 $ and prey $ x_2 $
    $ m_1, m_2 $ Refuge rates of prey $ x_1 $ and prey $ x_2 $
    $ e_1, e_2 $ Conversion factors for prey $ x_1 $ and prey $ x_2 $
    $ q_1, q_2, q_3 $ Captureability factors for prey $ x_1 $, prey $ x_2 $ and predator $ y $
    $ E_1, E_2, E_3 $ Harvesting efforts for prey $ x_1 $, prey $ x_2 $ and predator $ y $
    $ d $ Predator $ y $ mortality rate

     | Show Table
    DownLoad: CSV

    Most species in nature, including humans, are influenced by fear. Fear may cause an abnormal state and behavior to arise. As usual, prey have an innate fear of predators. The ecology of fear is related to combining the optimal behavior of prey and predators with their population densities [25,26]. In view of reality, it is a fact that prey fear predators, which is seen as a psychological effect that can have a lasting impact on prey populations. This psychological influence is often easy to overlook, but it is necessary to consider it in the context of practical ecology [27]. Wang et al. [28] first considered the effect of the fear factor on the model and first proposed the fear of prey $ F(k, y) $. Afterwards, some researchers have investigated the effects of the fear effect and predator interferences in some three-dimensional systems as well as explored the generation of Hopf bifurcation conditions in the presence of a fear parameter as a bifurcation parameter [29,30,31,32]. Zanette et al. [33] observed that prey will reduce reproducing because of fear of being killed by predators, thus decreasing the risk of being killed after giving birth, which also leads directly to a decline in prey birth rates. According to the above discussion, our paper considers the different fears $ k_i $ caused by predators for the two prey species.

    In reality, when prey feel the crisis of being hunted, they will reproduce less and increase their survival rate. These conditions about the fear factor $ F(k_i, y) $ $ (i = 1, 2) $ are listed as follows:

    $ \mathbf{1}) $ $ F(0, y) = 1 $: prey production does not decrease when the prey does not fear the predator;

    $ \mathbf{2}) $ $ F(k_i, 0) = 1 $: even though the prey will develop a fear of predators and there will be no predators, prey production will still not decline;

    $ \mathbf{3}) $ $ \lim\limits_{k_i\rightarrow \infty} F(k_i, y) = 0 $: when the prey's fear of the predator is very high, this will result in the prey production tending to zero;

    $ \mathbf{4}) $ $ \lim\limits_{y\rightarrow \infty} F(k_i, y) = 0 $: prey have a fear of predators, and when predator numbers are too large, this can also lead to prey production tending to zero;

    $ \mathbf{5}) $ $ \frac{\partial F(k_i, y)}{\partial k_i} < 0 $: the greater the prey's fear of predators, the less productive it will be;

    $ \mathbf{6}) $ $ \frac{\partial F(k_i, y)}{\partial y} < 0 $: predators are inversely proportional to their prey.

    For ease of analysis, we draw on Wang et al. [28] to consider the fear effect:

    $ F(ki,y)=11+kiy(i=1,2),
    $
    (1.2)

    obviously, $ F(k_i, y) $ $ (i = 1, 2) $ in (1.2) satisfies conditions $ \mathbf{1}) $–$ \mathbf{6}) $. Based on the above conditions, this study will consider the effect of fear on system (1.1) to obtain system (1.3).

    $ {dx1dt=r1x11+k1y(1x1K1)a1x1x2c1(1m1)x1yq1E1x1,dx2dt=r2x21+k2y(1x2K2)a2x1x2c2(1m2)x2yq2E2x2,dydt=e1(1m1)x1y+e2(1m2)x2ydyq3E3y.
    $
    (1.3)

    Notably, most biological parameters in much of the literature are fixed constants. However, in reality, the survival of species is full of unknowns, and all data are not always constant, which can lead to deviations from the ideal model with fixed parameters. In order to make the model more relevant and the results more accurate, we cannot just consider fixed parameters. Therefore, to make the study more convincing, it is necessary to target imprecise parameters. Professor Zadeh [34], who first proposed the fuzzy set theory, also argued that the application of fuzzy differential equations is a more accurate method for modeling biological dynamics in the absence of accurate data conditions [35]. Moreover, the first introduction of the idea of fuzzy derivatives came from Chang and Zadeh [36]. Further, Kaleva [37] studied the generalized fuzzy derivatives based on Hukuhara differentiability, the Zadeh extension principle, and the strong generalized differentiability concept. Bede et al. [38] employed the notion of strongly generalized differentiability to investigate fuzzy differential equations. Khastan and Nieto [39] solved the margin problem for fuzzy differential equations in their article. Motivated by the method of Pal [13] and Wang [23], we assume that the imprecise parameters $ \widetilde{r_1} $, $ \widetilde{r_2} $, $ \widetilde{a_1} $, $ \widetilde{a_2} $, $ \widetilde{c_1} $, $ \widetilde{c_2} $, $ \widetilde{e_1} $, $ \widetilde{e_1} $ and $ \widetilde{d} $ represent all triangular fuzzy numbers (the relevant theories of fuzzy sets are detailed in Appendix A), then the system (1.3) can be written as

    $ {~dx1dt=~r1x11+k1y(1x1K1)~a1x1x2~c1(1m1)x1yq1E1x1,~dx2dt=~r2x21+k2y(1x2K2)~a2x1x2~c2(1m2)x2yq2E2x2,~dydt=~e1(1m1)x1y+~e2(1m2)x2y˜dyq3E3y,
    $
    (1.4)

    we cut these imprecise parameters $ \widetilde{r_1} $, $ \widetilde{r_2} $, $ \widetilde{a_1} $, $ \widetilde{a_2} $, $ \widetilde{c_1} $, $ \widetilde{c_2} $, $ \widetilde{e_1} $, $ \widetilde{e_1} $, and $ \widetilde{d} $ by using $ \alpha $-level. System (1.4) can be expressed as follows:

    $ {(dx1dt)αL=rα1Lx11+k1yrα1R1+k1yx21K1aα1Rx1x2cα1R(1m1)x1yq1E1x1,(dx1dt)αR=rα1Rx11+k1yrα1L1+k1yx21K1aα1Lx1x2cα1L(1m1)x1yq1E1x1,(dx2dt)αL=rα2Lx21+k2yrα2R1+k2yx22K2aα2Rx1x2cα2R(1m2)x2yq2E2x2,(dx2dt)αR=rα2Rx21+k2yrα2L1+k2yx22K2aα2Lx1x2cα2L(1m2)x2yq2E2x2,(dydt)αL=eα1L(1m1)x1y+eα2L(1m2)x2ydαRyq3E3y,(dydt)αR=eα1R(1m1)x1y+eα2R(1m2)x2ydαLyq3E3y.
    $
    (1.5)

    Introducing weighted sum, we change (1.5) to (1.6)

    $ {dx1dt=w1(dx1dt)αL+w2(dx1dt)αR,dx2dt=w1(dx2dt)αL+w2(dx2dt)αR,dydt=w1(dydt)αL+w2(dydt)αR,
    $
    (1.6)

    where $ w_1 $ and $ w_2 $ are satisfied with $ w_1+w_2 = 1 $, and $ w_1, w_2\geq0 $. Simplifying the system (1.6), we obtain

    $ {dx1dt=A11+k1yx1A21+k1yx21K1A3x1x2A4(1m1)x1yq1E1x1,dx2dt=B11+k2yx2B21+k2yx22K2B3x1x2B4(1m2)x2yq2E2x2,dydt=C1(1m1)x1y+C2(1m2)x2yC3yq3E3y,
    $
    (1.7)

    where

    $ A1=w1rα1L+w2rα1R,  A2=w1rα1R+w2rα1L,  A3=w1aα1R+w2aα1L,A4=w1cα1R+w2cα1L,  B1=w1rα2L+w2rα2R,  B2=w1rα2R+w2rα2L,B3=w1aα2R+w2aα2L,  B4=w1cα2R+w2cα2L,  C1=w1eα1L+w2eα1R,C2=w1eα2L+w2eα2R,  C3=w1dαR+w2dαL.
    $

    The rest of the paper is shown below: In Section 2, we first prove the nonnegativity and boundedness of the system (1.7). Sections 3 and 4 discuss all possible equilibria and give conditions for the local asymptotic stability and global asymptotic stability of the equilibria. Immediately after that, in Section 5, we analyze the Hopf bifurcation by using the normal form theory. In Section 6, we numerically simulate the theoretical results of Sections 4 and 5. Finally, the article ends with detailed conclusions.

    In this section, we give the following theorem to ensure the boundedness and nonnegativity of the solutions of the system (1.7).

    Theorem 2.1. Provided that the initial values $ x_1(0) > 0 $, $ x_2(0) > 0 $, and $ y(0) > 0 $, all solutions of system (1.7) are nonnegative.

    Proof. It is not difficult to find that the right half of the system (1.7) fulfills the local Lipschitzian condition. Integrating both sides of the system (1.7) at the same time yields

    $ x1(t)=x1(0)[expt0(A11+k1yA21+k1yx1K1A3x2A4(1m1)yq1E1)ds]>0,x2(t)=x2(0)[expt0(B11+k2yB21+k2yx2K2B3x1B4(1m2)yq2E2)ds]>0,y(t)=y(0)[expt0(C1(1m1)x1C2(1m2)x2C3q3E3)ds]>0.
    $
    (2.1)

    If the solution curve starts at any internal point of $ \mathbb{R}^{3}_{+} = \{(x_1(t), x_2(t), y(t))\in\mathbb{R}^{3}:x_1(t)\geq0 $, $ x_2(t)\geq0 $, $ y(t)\geq0\} $, then $ x_1(t) $, $ x_2(t) $, and $ y(t) $ will always be nonnegative.

    Theorem 2.2. Assume that the initial values $ x_1(0) $, $ x_2(0) $, and $ y(0) $ are all greater than zero. The feasible region $ \Omega $ is a positive invariant set of the system (1.7) defined by

    $ Ω={(x1(t),x2(t),y(t))R3+:C1A4x1(t)+C2B4x2(t)+y(t)ϕμ},
    $
    $ (i.e  Ω={(x1(t),x2(t),y(t))R3+:w1eα1L+w2eα1Rw1cα1R+w2cα1Lx1(t)+w1eα2L+w2eα2Rw1cα2R+w2cα2Lx2(t)+y(t)ϕμ}),
    $

    where $ \mu = \mathit{\mbox{min}}\{q_1E_1, q_2E_2, C_3+q_3E_3\} $.

    Proof. Define a function

    $ W(t)=C1A4x1(t)+C2B4x2(t)+y(t).
    $
    (2.2)

    After taking the derivative on both sides of (2.2), we obtain

    $ dWdt=C1A4dx1dt+C2B4dx2dt+dydt.
    $
    (2.3)

    Furthermore, we can obtain

    $ dWdt+μW=C1A4(1+k1y)(A1x1A2x21K1)C1A3A4x1x2C1(1m1)x1ydWdt+μW=+C2B4(1+k2y)(B1x2B2x22K2)C2B3B4x1x2C2(1m2)x2ydWdt+μW=+C1(1m1)x1y+C2(1m2)x2yC3yq3E3y+C1μA4x1dWdt+μW=C1A4q1E1x1C2B4q2E2x2+C2μB4x2+μy,dWdt+μW=C1A4(1+k1y)(A1x1A2x21K1)+C2B4(1+k2y)(B1x2B2x22K2)dWdt+μW=+C1A4x1(μq1E1)+C2B4x2(μq2E2)+y(μC3q3E3)dWdt+μW=(C1A3A4+C2B3B4)x1x2,
    $
    (2.4)

    where $ \mu = \mbox{min}\{q_1E_1, q_2E_2, C_3+q_3E_3\} $. Let $ \phi_1 = \frac{A_1^{2}K_1}{4A_2}, \phi_2 = \frac{B_1^{2}K_2}{4B_2} $, $ \phi\equiv\frac{C_1}{A_4}\phi_1+\frac{C_2}{B_4}\phi_2 $, we have

    $ dWdt+μWC1A4ϕ1+C2B4ϕ2=ϕ.
    $
    (2.5)

    Therefore, it can be deduced that

    $ Wϕμ+Neμt,
    $
    (2.6)

    where $ N $ is a positive constant. Then we can further obtain

    $ lim suptWϕμ,
    $
    (2.7)

    which indicates that the feasible domain $ \Omega $ is a positive invariant set.

    In this section, we discuss the existence of all equilibria in the system (1.7). All equilibria for system (1.7) are provided by

    (1) Trivial equilibrium $ P_1 = (0, 0, 0) $.

    (2) Axial equilibrium $ P_2 = (x_1^{\varsigma}, 0, 0) $ exists if $ A_1 > q_1E_1 $, where $ x_1^{\varsigma} = \frac{K_1(A_1-q_1E_1)} {A_2} $.

    (3) Axial equilibrium $ P_3 = (0, x_2^{\Upsilon}, 0) $ exists if $ B_1 > q_2E_2 $, where $ x_2^{\Upsilon} = \frac{K_2(B_1-q_2E_2)}{B_2} $.

    (4) Axial equilibrium $ P_4 = (x_1^{\Psi}, 0, y^{\Psi}) $ exists if $ A_1\geq\frac{A_2x_1^{\Psi}}{K_1}+q_1E_1 $ and $ \sqrt{\Delta_1} > A_4(1-m_1)+k_1q_1E_1 $, where

    $ xΨ1=C3+q3E3C1(1m1),yΨ=Δ1(A4(1m1)+k1q1E1)2k1A4(1m1),Δ1=4(k1A4(1m1))(A1A2xΨ1K1q1E1)+(A4m1A4k1q1E1)2.
    $
    (3.1)

    (5) Axial equilibrium $ P_5 = (0, x_2^{\nu}, y^{\nu}) $ exists if $ B_1\geq\frac{B_2x_2^{\nu}}{K_2}+q_2E_2 $ and $ \sqrt{\Delta_2} > B_4(1-m_2)+k_2q_2E_2 $, where

    $ xν2=C3+q3E3C2(1m2),yν=Δ2(B4(1m2)+k2q2E2)2k2B4(1m2),Δ2=4(k2B4(1m2))(B1B2xν2K2q2E2)+(B4m2B4k2q2E2)2.
    $
    (3.2)

    (6) Axial equilibrium $ P_6 = (x_1^{\iota}, x_2^{\iota}, 0) $ exists if $ B_1 > B_3x_1^{\iota}+q_2E_2 $, $ K_1K_2A_3B_3 > A_2B_2 $ and $ A_3B_1K_2+B_2q_1E_1 > A_3K_2q_2E_2+A_1B_2 $, where

    $ xι1=A3K1K2(B1q2E2)+K1B2(q1E1A1)K1K2A3B3A2B2,xι2=K2(B1B3xι1q2E2)B2.
    $
    (3.3)

    (7) Internal equilibrium $ P_7 = (x_1^{*}, x_2^{*}, y^{*}) $ exists, and its value will be given in the proof of Theorem 3.1.

    Theorem 3.1. When $ g_3 > 0 $ and $ g_4g_5 < 0 $ are met, there is an internal equilibrium $ P_7 $.

    Proof. We derive that from the second equation of the system (1.7)

    $ g1y2+g2y+g3=0,
    $
    (3.4)

    where

    $ g_1 = -k_2B_4(1-m_2),\ \ g_2 = -k_2B_3x_1-B_4(1-m_2)-q_2E_2k_2,\ \ g_3 = \left(B_1-\frac{x_2}{K_2}B_2-B_3x_1-q_2E_2\right). $

    It follows from the Descartes law of signs that Eq (3.4) has one and only one solution $ y^{*} $ greater than zero if and only if $ g_3 > 0 $, i.e., $ B_1 > \frac{x_2}{K_2}B_2+B_3x_1+q_2E_2 $. Substituting $ y^{*} $ into the algebra expression on the right side of the first equation of the system (1.7) equals zero; furthermore, we obtain

    $ x2=A1A3+A3k1yA2x1K1A3+K1A3k1yA4(1m1)yA3q1E1A3.
    $
    (3.5)

    Introduce (3.5) into the right side of the third equation of the system (1.7), which meets zero, it simplifies to obtain

    $ g4x1g5=0,
    $
    (3.6)

    where

    $ g4=[C1(1m1)A2C2(1m2)K1A3+K1A3k1y],g5=[A1C2(1m2)A3+A3k1yA4C2(1m1)(1m2)yA3C2q3E3(1m2)A3C3q3E3].
    $

    Reusing the Descartes law of signs, we can assert that there exists at least one positive solution $ x_1^{*} $ of Eq (3.6) if and only if $ g_4g_5 < 0 $. And then we can deduce that

    $ x_2^{*} = \frac{A_1}{A_3+A_3k_1y^{*}}-\frac{A_2x_1^{*}}{K_1A_3+K_1A_3k_1y^{*}}-\frac{A_4(1-m_1)y^{*}}{A_3}-\frac{q_1E_1}{A_3}. $

    then the interior equilibrium $ P_7(x_1^{*}, x_2^{*}, y^{*}) $ exists.

    In this section, the Jocabian matrix will be used to prove the local stability of all equilibria. Moreover, we prove the global stability of the internal equilibrium $ P_7 $ by constructing a Lyapunov function.

    The Jocabian matrix for system (1.7) is given below:

    $ M=(M11M12M13M21M22M23M31M32M33),
    $
    (4.1)

    where

    $ M11=A11+k1y2A21+k1yx1K1A3x2A4(1m1)yq1E1,M12=A3x1,M13=k1A1x1(1+k1y)2+k1A2(1+k1y)2x21K1A4(1m1)x1,M21=B3x2,M22=B11+k2y2B21+k2yx2K2B3x1B4(1m2)yq2E2,M23=k2B1x2(1+k2y)2+k2B2(1+k2y)2x22K2B4(1m2)x2,M31=C1(1m1)y,M32=C2(1m2)y,  M33=C1(1m1)x1+C2(1m2)x2C3q3E3.
    $
    (4.2)

    Through simple calculation, we directly draw the conclusion that trivial and axial equilibria are locally asymptotically stable:

    $ \mathbf{(1)} $ $ P_1(0, 0, 0) $ is locally asymptotically stable if

    $ A1q1E1<0andB1q2E2<0.
    $
    (4.3)

    $ \mathbf{(2)} $ $ P_2(x_1^{\varsigma}, 0, 0) $ is locally asymptotically stable if

    $ B1q2E2<B3xς1q2<B3(C3+q3E3)C1(1m1)q2.
    $
    (4.4)

    $ \mathbf{(3)} $ $ P_3(0, x_2^{\Upsilon}, 0) $ is locally asymptotically stable if

    $ A1q1E1<A3xΥ2q1<A3(C3+q3E3)C2(1m2)q1.
    $
    (4.5)

    $ \mathbf{(4)} $ $ P_4(x_1^{\Psi}, 0, y^{\Psi}) $ is locally asymptotically stable if

    $ B1q2(1+k2yΨ)E2<B3xΨ1+B4(1m2)yΨq2.
    $
    (4.6)

    $ \mathbf{(5)} $ $ P_5(0, x_2^{\nu}, y^{\nu}) $ is locally asymptotically stable if

    $ A1q1(1+k1yν)E1<A3xν2+A4(1m1)yνq1.
    $
    (4.7)

    $ \mathbf{(6)} $ $ P_6(x_1^{\iota}, x_2^{\iota}, 0) $ is locally asymptotically stable if

    $ C1(1m1)xι1+C2(1m2)xι2<C3+q3E3.
    $
    (4.8)

    We draw the conclusion that the internal equilibrium $ P_7(x_1^{*}, x_2^{*}, y^{*}) $ is locally asymptotically stable from the proof of Theorem 4.1.

    Theorem 4.1. The internal equilibrium $ P_7 $ is locally asymptotically stable if it exists and the following conditions are fulfilled:

    $ ψ1>0,ψ1ψ2>0,ψ3>0.
    $
    (4.9)

    Proof. The Jocabian matrix of system (1.7) at $ (x_1^{*}, x_2^{*}, y^{*}) $ is

    $ (L1L2L3L4L5L6L7L8L9),
    $
    (4.10)

    where

    $ L1=A2(1+k1y)x1K1<0,L2=A3x1<0,L3=k1A1x1(1+k1y)2A4(1m1)x1+k1A2(1+k1y)2(x1)2K1,L4=B3x2<0,L5=B2(1+k2y)x2K2<0,L6=k2B1x2(1+k2y)2B4(1m2)x2+k2B2(1+k2y)2(x2)2K2,L7=C1(1m1)y>0,L8=C2(1m2)y>0,L9=0.
    $
    (4.11)

    Therefore, the characteristic equation at $ P_7 $ can be expressed as

    $ η3+ψ1η2+ψ2η+ψ3=0,
    $
    (4.12)

    where

    $ ψ1=L1L5,ψ2=L1L5L6L8L3L7L2L4,ψ3=L8(L1L6L3L4)+L7(L3L5+L2L6).
    $
    (4.13)

    The Routh-Hurwitz criterion shows that the internal equilibrium $ P_7 $ is locally asymptotically stable; the following conditions need to be met: $ \psi_1 > 0 $, $ \psi_1\psi_2 > 0 $, and $ \psi_3 > 0 $.

    This subsection studies the global asymptotic stability of interior equilibrium $ P_7 $.

    Theorem 4.2. If condition $ 4\Gamma_1\Gamma_2l_1l_2A_2B_2(1+k_1y)(1+k_2y)$ > $(l_1A_3+l_2B_3)^{2} $ (i.e. $ 4\Gamma_1\Gamma_2l_1l_2(w1rα1R+w2rα1L)

    (w_1r_{2R}^{\alpha}+w_2r_{2L}^{\alpha})(1+k1y)(1+k2y)$>$(l1(w1aα1R+w2aα1L)
    +l_2(w_1a_{2R}^{\alpha}$$+w_2a_{2L}^{\alpha}))^{2} $) holds, then $ P_7 $ is globally asymptotically stable.

    Proof. We construct a Lyapunov function:

    $ V(x1,x2,y)=l1[x1x1x1ln(x1x1)]+l2[x2x2x2ln(x2x2)]  +yyyln(yy).
    $
    (4.14)

    Obviously, $ x_i-x_i^{*}-x_i^{*}\ln(\frac{x_i}{x_i^{*}})\geq0\ (i = 1, 2) $ and $ y-y^{*}-y^{*}\ln(\frac{y}{y^{*}})\geq0 $, thus $ V\geq0 $. Taking the derivative of $ V(x_1, x_2, y) $ over $ t $, one has

    $ dVdt=l1(x1x1x1)dx1dt+l2(x2x2x2)dx2dt+yyydydt,
    $
    (4.15)

    where

    $ x1x1x1dx1dt=k1A1(1+k1y)(1+k1y)(x1x1)(yy)A2(x1x1)2K1(1+k1y)(1+k1y)x1x1x1dx1dt=A2k1(x1yx1y)K1(1+k1y)(1+k1y)(x1x1)A3(x1x1)(x2x2)x1x1x1dx1dt=A4(1m1)(x1x1)(yy),x2x2x2dx2dt=k2B1(1+k2y)(1+k2y)(x2x2)(yy)B2(x2x2)2K2(1+k2y)(1+k2y)x2x2x2dx2dt=B2k2(x2yx2y)K2(1+k2y)(1+k2y)(x2x2)B3(x1x1)(x2x2)x2x2x2dx2dt=B4(1m2)(x2x2)(yy),yyydydt=C1(1m1)(x1x1)(yy)+C2(1m2)(x2x2)(yy).
    $
    (4.16)

    To simplify the calculation, let

    $ x1yx1y=y(x1x1)x1(yy),  x2yx2y=y(x2x2)x2(yy),Γ1=1K1(1+k1y)(1+k1y),  Γ2=1K2(1+k2y)(1+k2y),l1=C1(1m1)Γ1k1(K1A1+A2x1)+A4(1m1),  l2=C2(1m2)Γ2k2(K2B1+B2x2)+A4(1m1).
    $
    (4.17)

    We obtain

    $ dVdt={Γ1l1A2(1+k1y)(x1x1)2+(l1A3+l2B3)(x1x1)(x2x2)+Γ2l2B2(1+k2y)(x2x2)2}dVdt  =YGY,
    $
    (4.18)

    where

    $ Y^{'} = \left[(x_1-x_1^{*}),(x_2-x_2^{*})\right],\ \ G = (Γ1l1A2(1+k1y)l1A3+l2B32l1A3+l2B32Γ2l2B2(1+k2y))
    . $

    Therefore, $ \frac{dV}{dt} < 0 $ if and only if $ 4\Gamma_1\Gamma_2l_1l_2A_2B_2(1+k_1y)(1+k_2y) > (l_1A_3+l_2B_3)^{2} $.

    In this section, we will use the normal form theory introduced by Hassard et al.[40] and the central manifold theory [41] to study the Hopf bifurcation of the system (1.7). When the system (1.7) undergoes Hopf bifurcation, the corresponding characteristic equation must have a pair of conjugate pure imaginary roots, that is,

    $ η1,2=±iω,i=1.
    $
    (5.1)

    Consider the parameter $ k_1 $ as a bifurcation parameter. When the value of parameter $ k_{1} $ changes near the critical point $ k_{1}^{\Xi} $ of Hopf bifurcation, the pure imaginary roots $ \pm i\omega $ will become a complex eigenvalue $ \eta = \rho+i\widetilde{\omega} $. Substituting $ \eta = \rho+i\widetilde{\omega} $ into Eq (4.12), we need to separate the imaginary and real parts to get

    $ ρ3+ψ3+ρψ2+ρ2ψ13ρ˜ω2ψ1˜ω2=0,
    $
    (5.2)
    $ 3ρ2˜ω+ψ2˜ω+2ρψ1˜ω˜ω3=0.
    $
    (5.3)

    By simplifying Eqs (5.2) and (5.3), we obtain

    $ ψ38ρ32ρψ28ρ2ψ1ψ1ψ22ρψ21=0,
    $
    (5.4)

    at $ k_1 = k_1^{\Xi} $, taking the derivative of Eq (5.4) over $ k_1 $ yields

    $ dρdk1|k1=kΞ1=12(dψ3dk1ψ1dψ2dk1ψ2dψ1dk1)/(ψ2+ψ21).
    $
    (5.5)

    If it satisfies $ \frac{d\rho}{dk_1}\big|_{k_1 = k_1^{\Xi}}\neq0 $, the system (1.7) will generate Hopf bifurcation, which indicates that when parameter $ k_1 $ crosses the bifurcation critical point $ k_1^{\Xi} $, the population state evolves from stable equilibrium to periodic oscillation over time.

    When the system (1.7) undergoes Hopf bifurcation at $ k_1 = k_1^{\Xi} $, the final decision condition is also met. Considering that the characteristic roots of Eq (4.12) are $ \eta_{1, 2} = \pm i\omega $ and $ \eta_3 = -\psi_1 $, in order to obtain this condition, we introduce

    $ z1=x1x1,z2=x2x2,z3=yy.
    $
    (5.6)

    Substituting (5.6) into the system (1.7) and separating the linear and nonlinear parts, it can be obtained that

    $ (˙z1˙z2˙z3)=J(P7)(z1z2z3)+(F1(z1,z2,z3)F2(z1,z2,z3)F3(z1,z2,z3)),
    $
    (5.7)

    where

    $ F1(z1,z2,z3)=2j1+j2+j33tj1j2j3zj11zj22zj33+O((|z1|+|z2|+|z3|)4),F2(z1,z2,z3)=2j1+j2+j33nj1j2j3zj11zj22zj33+O((|z1|+|z2|+|z3|)4),F3(z1,z2,z3)=2j1+j2+j33lj1j2j3zj11zj22zj33+O((|z1|+|z2|+|z3|)4),
    $
    (5.8)

    where $ O((|z_1|+|z_2|+|z_3|)^{4}) $ is a fourth-order polynomial function about variables $ (|z_1|, |z_2|, |z_3|) $, while $ t_{j_1j_2j_3} $, $ n_{j_1j_2j_3} $, and $ l_{j_1j_2j_3} $ can be obtained through calculation:

    $ t101=A1k12(1+k1y)2+A2k1(1+k1y)2x1K1A4(1m1)2,  t002=A1x1k21(1+k1y)3A2k21(1+k1y)3(x1)2K1,t102=A1k213(1+k1y)3A2k213(1+k1y)32x1K1,  t003=A1x1k31(1+k1y)4+A2k31(1+k1y)4(x1)2K1,t110=A32,  t200=A21+k1y1K1,  t201=A2k13(1+k1y)21K1,t011=t020=t030=t012=t021=t111=t120=t210=t300=0,n011=B1k22(1+k2y)2+B2k2(1+k2y)2x2K2B4(1m2)2,  n002=B1x2k22(1+k2y)3B2k22(1+k2y)3(x2)2K2,n012=B1k223(1+k2y)3B2k223(1+k2y)32x2K2,  n003=B1x2k32(1+k2y)4+B2k32(1+k2y)4(x2)2K2,n110=B32,  n020=B21+k2y1K2,  n021=B2k23(1+k2y)21K2,n101=n200=n030=n111=n120=n102=n210=n201=n300=0,l101=C1(1m1)2,  l011=C2(1m2)2,l110=l200=l020=l002=l030=l003=l012=l021=l111=l120=l102=l210=l201=l300=0.
    $
    (5.9)

    By introducing a reversible transformation

    $ (z1z2z3)=(101q21q22q23q31q32q33)(y1y2y3),
    $
    (5.10)

    which $ q_{21}q_{32} - q_{22}q_{31} + q_{22}q_{33} - q_{23}q_{32}\neq0 $, the expression for the coefficient $ q_{ij}(i = 2, 3;j = 1, 2, 3) $ is

    $ q21=L2L3L4L6+L1L3L5L6L1L2L26+L3L6ω2L23L4L5L23L252L2L3L5L6+L22L26+L23ω2,q22=(L23L4L1L3L6+L3L5L6L2L6)ωL23L252L2L3L5L6+L22L26+L23ω2,q23=L1L6L6η3L3L4L3L5L2L6L3η3,q31=L2L3L4L5L1L3L25L22L4L6+L1L2L5L6L1L3ω2L2L6ω2L23L252L2L3L5L6+L22L26+L23ω2,q32=ω(L1L2L6+L2L5L6L2L3L4L3L25L3ω2)L23L252L2L3L5L6+L22L26+L23ω2,q33=L2L4L1L5+L1η3+L5η3η23L3L5L2L6L3η3,
    $
    (5.11)

    $ L_{i} (i = 1, 2, \cdots, 9) $ is defined in (4.11). So the standard type of system (5.7) can be written as

    $ (˙y1˙y2˙y3)=(0ω0ω0000η3)(y1y2y3)+(~F1(y1,y2,y3)~F2(y1,y2,y3)~F3(y1,y2,y3)),
    $
    (5.12)

    where

    $ (~F1(y1,y2,y3)~F2(y1,y2,y3)~F3(y1,y2,y3))=(101q21q22q23q31q32q33)1(F1(z1,z2,z3)F2(z1,z2,z3)F3(z1,z2,z3)).
    $
    (5.13)

    In Eq (5.13), the coefficients of polynomial $ \widetilde{F}(y_1, y_2, y_3) $ are $ \widetilde{t}_{j_1j_2j_3} $, $ \widetilde{n}_{j_1j_2j_3} $ and $ \widetilde{l}_{j_1j_2j_3} $.

    Based on the central manifold theory, use the center manifold $ W^{c}(0, 0, 0) $ existing at the origin to reduce the dimension of the system (5.12), that is

    $ Wc(0,0,0)={(y1,y2,y3)R3|y3=h(y1,y2),h(0,0)=0,Dh(0,0)=0},
    $
    (5.14)

    where

    $ h(y1,y2)=1v1+v23hv1v2y1y2+O((|y1|+|y2|)4).
    $
    (5.15)

    Through calculation, we can obtain

    $ h10=h01=h11=0,h20=˜l200+˜l020ω2η3(1+ω4),h02=˜l020˜l200ω21+ω4,h30=1(1+ω4)(η23+ω6)(˜l200˜n101˜l300η32˜l020˜n200η3ω+˜l020˜n101ω2+˜l011˜l020ω3h30=˜l030ω3+2˜l200˜n200η3ω3˜l300η3ω4˜l011˜l200ω5˜l030ω7),h21=(˜l200+˜l020ω2)(˜l011η3+2˜t200η3ω+˜n110η23ω32˜t110ω4+˜n101η3ω52˜n020η3ω6)η3(1+ω4)(η23+ω6)h21=˜l020(2˜n110η3ω+˜n101ω32˜n020ω4)η23+ω6,h12=(˜l200+˜l020ω2)(2˜t110η3ω˜n101η23ω2+˜l011ω3+2˜n020η23ω3+2˜t200ω4+2˜n110η3ω6)η3(1+ω4)(η23+ω6)h12=˜l020(2˜n020η3ω+2˜n110ω4˜n101η3)η23+ω6,h03=(˜l200+˜l020ω2)(˜l011η23+ω2+˜n101ω3+2˜n200η3ω6)η3(1+ω4)(η23+ω6)h03=+˜l011˜l020η3˜l030η3+˜l300ω3+2˜l020˜n200ω4η23+ω6.
    $
    (5.16)

    Correspondingly, the dynamic properties of the system are limited to the central flow $ W^{c}(0, 0, 0) $, and in conjunction with Eq (5.14), system (5.12) can be simplified as

    $ (˙y1˙y2)=(U(y1,y2)N(y1,y2)),
    $
    (5.17)

    where

    $ U(y1,y2)=ωy2+1y21+2y1y2+3y22+4y31+5y21y2+6y1y22+7y32,N(y1,y2)=ωy1+ȷ1y21+ȷ2y1y2+ȷ3y22+ȷ4y31+ȷ5y21y2+ȷ6y1y22+ȷ7y32,
    $
    (5.18)

    in the formula, we have

    $ 1=˜t200,2=˜t110,3=˜t020,4=h20˜t101+˜t300,5=˜t210,6=h02˜t101,7=0,ȷ1=˜n200,ȷ2=˜n110,ȷ3=˜n020,ȷ4=˜n300,ȷ5=h20˜n011,ȷ6=0,ȷ7=h02˜n011+˜n030.
    $
    (5.19)

    We introduce the partial derivative sign

    $ Uy1(ykΞ1)=Uy1,  3Uy21y2(ykΞ1)=Uy1y1y2,  2Ny22(ykΞ1)=Ny2y2,  ,
    $
    (5.20)

    where subscripts $ y_1 $ and $ y_2 $ indicate partial derivatives for the first and second variable, respectively. Based on Eq (5.18), it can be obtained that $ U_{y_1} = 0 $, $ U_{y_2}\neq0 $, $ N_{y_1}\neq0 $, $ N_{y_2} = 0 $, and $ U_{y_2}N_{y_1}\neq0 $. In addition, it ensures that the system (5.18) has pure virtual feature roots $ \pm i \sqrt{|U_{y_2}N_{y_1}|} $. Thus, it can be determined that system (1.7) produces Hopf bifurcation; the direction of the bifurcation is determined by the following equation:

    $ QkΞ1=116ω(3+5+ȷ5+ȷ7)+116ω(13ȷ2ȷ3ȷ1ȷ2ȷ11).
    $
    (5.21)

    Theorem 5.1. If $ \frac{d\rho}{dk_1}\big|_{k_1 = k_1^{\Xi}}\neq0 $, then system (1.7) will generate Hopf bifurcation at interior equilibrium $ P_7 $. In addition, when $ \frac{d\rho}{dk_1}\big|_{k_1 = k_1^{\Xi}} < 0 $, if $ Q_{k_1^{\Xi}} < 0 $ and $ 0 < k_1-k_1^{\Xi}\ll1 $, then system (1.7) will generate supercritical Hopf bifurcation and form a stable periodic orbit, or if $ Q_{k_1^{\Xi}} > 0 $ and $ 0 < k_1-k_1^{\Xi}\ll1 $, then system (1.7) will generate subcritical Hopf bifurcation and form a stable periodic orbit.

    In this section, we first discussed equilibria $ P_1 $ to $ P_7 $ of system (1.7) with distinct values of $ \alpha $, $ w_1 $, and $ w_2 $. Consider the parameter values as follows: $ \widetilde{r_1} = (2.8, 3, 3.2) $, $ \widetilde{r_2} = (2.8, 3, 3.2) $, $ \widetilde{c_1} = (0.1, 0.2, 0.3) $, $ \widetilde{c_2} = (0.5, 0.6, 0.7) $, $ \widetilde{a_1} = (0.1, 0.2, 0.3) $, $ \widetilde{a_2} = (0.2, 0.3, 0.4) $, $ \widetilde{e_1} = (0.2, 0.3, 0.4) $, $ \widetilde{e_2} = (0.3, 0.4, 0.5) $, and $ \widetilde{d} = (0.1, 0.2, 0.3) $. Tables 28 showed that the trivial equilibrium $ P_1 $ retained constant at (0, 0, 0), the values of prey $ x_1 $, prey $ x_2 $, and predator $ y $ always maintained at $ 0 $; the values of prey $ x_1 $ in $ P_2 $ and prey $ x_2 $ in $ P_3 $ severally decreased with increasing $ w_1 $ under the same $ \alpha $; the values of prey $ x_1 $ and predator $ y $ in $ P_4 $ increased with increasing $ w_1 $, and for $ P_5 $ the value of prey $ x_2 $ and predator $ y $ rose with growing $ w_1 $; the values of prey $ x_1 $ and $ x_2 $ in $ P_6 $ decreased with growing $ w_1 $; and for the same $ \alpha $, considering interior equilibrium $ P_7 $, the values of prey $ x_1 $, prey $ x_2 $, and predator $ y $ decreased with growing $ w_1 $.

    Table 2.  The trivial equilibrium $ P_1 $ for $ k_1 = 0.1$, $k_2 = 0.7$, $q_1 = 0.7$, $q_2 = 0.5$, $q_3 = 0.7 $, $ E_1 = 0.3$, $E_2 = 0.2$, $E_3 = 0.2$, $K_1 = 5$, $K_2 = 5$, $m_1 = 0.9$, $m_2 = 0.3 $.
    $ w_1 $ $ w_2 $ $ P_1 $ at $ \alpha=0 $ $ P_1 $ at $ \alpha=0.3 $ $ P_1 $ at $ \alpha=0.6 $ $ P_1 $ at $ \alpha=0.9 $
    0 1 $ (0, 0, 0) $ $ (0, 0, 0) $ $ (0, 0, 0) $ $ (0, 0, 0) $
    0.2 0.8 $ (0, 0, 0) $ $ (0, 0, 0) $ $ (0, 0, 0) $ $ (0, 0, 0) $
    0.4 0.6 $ (0, 0, 0) $ $ (0, 0, 0) $ $ (0, 0, 0) $ $ (0, 0, 0) $
    0.6 0.4 $ (0, 0, 0) $ $ (0, 0, 0) $ $ (0, 0, 0) $ $ (0, 0, 0) $
    0.8 0.2 $ (0, 0, 0) $ $ (0, 0, 0) $ $ (0, 0, 0) $ $ (0, 0, 0) $
    1 0 $ (0, 0, 0) $ $ (0, 0, 0) $ $ (0, 0, 0) $ $ (0, 0, 0) $

     | Show Table
    DownLoad: CSV
    Table 3.  The axial equilibrium $ P_2 $ for $ k_1 = 0.1$, $k_2 = 0.7$, $q_1 = 0.7$, $q_2 = 0.5$, $q_3 = 0.7 $, $ E_1 = 0.3$, $E_2 = 0.2$, $E_3 = 0.2$, $K_1 = 5$, $K_2 = 5$, $m_1 = 0.9$, $m_2 = 0.3 $.
    $ w_1 $ $ w_2 $ $ P_2 $ at $ \alpha=0 $ $ P_2 $ at $ \alpha=0.3 $ $ P_2 $ at $ \alpha=0.6 $ $ P_2 $ at $ \alpha=0.9 $
    0 1 $ (5.3393, 0, 0) $ $ (5.1224, 0, 0) $ $ (4.9144, 0, 0) $ $ (4.7148, 0, 0) $
    0.2 0.8 $ (5.0521, 0, 0) $ $ (4.9280, 0, 0) $ $ (4.8069, 0, 0) $ $ (4.6888, 0, 0) $
    0.4 0.6 $ (4.7804, 0, 0) $ $ (4.7409, 0, 0) $ $ (4.7017, 0, 0) $ $ (4.6629, 0, 0) $
    0.6 0.4 $ (4.5230, 0, 0) $ $ (4.5608, 0, 0) $ $ (4.5988, 0, 0) $ $ (4.6372, 0, 0) $
    0.8 0.2 $ (4.2788, 0, 0) $ $ (4.3872, 0, 0) $ $ (4.4980, 0, 0) $ $ (4.6116, 0, 0) $
    1 0 $ (4.0469, 0, 0) $ $ (4.2197, 0, 0) $ $ (4.3994, 0, 0) $ $ (4.5861, 0, 0) $

     | Show Table
    DownLoad: CSV
    Table 4.  The axial equilibrium $ P_3 $ for $ k_1 = 0.7$, $k_2 = 0.1$, $q_1 = 0.7$, $q_2 = 0.5$, $q_3 = 0.7 $, $ E_1 = 0.3$, $E_2 = 0.2$, $E_3 = 0.2$, $K_1 = 5$, $K_2 = 5$, $m_1 = 0.3$, $m_2 = 0.9 $.
    $ w_1 $ $ w_2 $ $ P_3 $ at $ \alpha=0 $ $ P_3 $ at $ \alpha=0.3 $ $ P_3 $ at $ \alpha=0.6 $ $ P_3 $ at $ \alpha=0.9 $
    0 1 $ (0, 5.5357, 0) $ $ (0, 5.3147, 0) $ $ (0, 5.1027, 0) $ $ (0, 4.8993, 0) $
    0.2 0.8 $ (0, 5.2431, 0) $ $ (0, 5.1166, 0) $ $ (0, 4.9932, 0) $ $ (0, 4.8728, 0) $
    0.4 0.6 $ (0, 4.9662, 0) $ $ (0, 4.9260, 0) $ $ (0, 4.8861, 0) $ $ (0, 4.8465, 0) $
    0.6 0.4 $ (0, 4.7039, 0) $ $ (0, 4.7424, 0) $ $ (0, 4.7812, 0) $ $ (0, 4.8202, 0) $
    0.8 0.2 $ (0, 4.4551, 0) $ $ (0, 4.5655, 0) $ $ (0, 4.6785, 0) $ $ (0, 4.7942, 0) $
    1 0 $ (0, 4.2187, 0) $ $ (0, 4.3949, 0) $ $ (0, 4.5779, 0) $ $ (0, 4.7682, 0) $

     | Show Table
    DownLoad: CSV
    Table 5.  The axial equilibrium $ P_4 $ for $ k_1 = 0.3$, $k_2 = 0.7$, $q_1 = 0.7$, $q_2 = 0.5$, $q_3 = 0.7 $, $ E_1 = 0.3$, $E_2 = 0.2$, $E_3 = 0.2$, $K_1 = 20$, $K_2 = 20$, $m_1 = 0.2$, $m_2 = 0.6 $.
    $ w_1 $ $ w_2 $ $ P_4 $ at $ \alpha=0 $ $ P_4 $ at $ \alpha=0.3 $ $ P_4 $ at $ \alpha=0.6 $ $ P_4 $ at $ \alpha=0.9 $
    0 1 $ (0.4800, 0, 0.1496) $ $ (0.5838, 0, 0.1767) $ $ (0.7059, 0, 0.2040) $ $ (0.8516, 0, 0.2317) $
    0.2 0.8 $ (0.6222, 0, 0.1858) $ $ (0.6971, 0, 0.2022) $ $ (0.7802, 0, 0.2188) $ $ (0.8732, 0, 0.2354) $
    0.4 0.6 $ (0.8000, 0, 0.2224) $ $ (0.8306, 0, 0.2280) $ $ (0.8623, 0, 0.2335) $ $ (0.8954, 0, 0.2391) $
    0.6 0.4 $ (1.0286, 0, 0.2595) $ $ (0.9902, 0, 0.2539) $ $ (0.9534, 0, 0.2484) $ $ (0.9181, 0, 0.2428) $
    0.8 0.2 $ (1.3333, 0, 0.2969) $ $ (1.1845, 0, 0.2801) $ $ (1.0551, 0, 0.2633) $ $ (0.9415, 0, 0.2465) $
    1 0 $ (1.7600, 0, 0.3343) $ $ (1.4261, 0, 0.3062) $ $ (1.1692, 0, 0.2782) $ $ (0.9655, 0, 0.2502) $

     | Show Table
    DownLoad: CSV
    Table 6.  The axial equilibrium $ P_5 $ for $ k_1 = 0.3$, $k_2 = 0.7$, $q_1 = 0.7$, $q_2 = 0.5$, $q_3 = 0.7 $, $ E_1 = 0.3$, $E_2 = 0.2$, $E_3 = 0.2$, $K_1 = 20$, $K_2 = 20$, $m_1 = 0.8$, $m_2 = 0.6 $.
    $ w_1 $ $ w_2 $ $ P_5 $ at $ \alpha=0 $ $ P_5 $ at $ \alpha=0.3 $ $ P_5 $ at $ \alpha=0.6 $ $ P_5 $ at $ \alpha=0.9 $
    0 1 $ (0, 0.1920, 0.3638) $ $ (0, 0.2298, 0.3833) $ $ (0, 0.2727, 0.4029) $ $ (0, 0.3220, 0.4226) $
    0.2 0.8 $ (0, 0.2435, 0.3899) $ $ (0, 0.2697, 0.4016) $ $ (0, 0.2981, 0.4134) $ $ (0, 0.3291, 0.4252) $
    0.4 0.6 $ (0, 0.3048, 0.4160) $ $ (0, 0.3150, 0.4199) $ $ (0, 0.3255, 0.4239) $ $ (0, 0.3363, 0.4278) $
    0.6 0.4 $ (0, 0.3789, 0.4422) $ $ (0, 0.3668, 0.4383) $ $ (0, 0.3551, 0.4343) $ $ (0, 0.3437, 0.4304) $
    0.8 0.2 $ (0, 0.4706, 0.4683) $ $ (0, 0.4268, 0.4566) $ $ (0, 0.3872, 0.4448) $ $ (0, 0.3513, 0.4330) $
    1 0 $ (0, 0.5867, 0.4942) $ $ (0, 0.4970, 0.4748) $ $ (0, 0.4222, 0.4553) $ $ (0, 0.3590, 0.4357) $

     | Show Table
    DownLoad: CSV
    Table 7.  The axial equilibrium $ P_6 $ for $ k_1 = 0.3$, $k_2 = 0.4$, $q_1 = 0.7$, $q_2 = 0.5$, $q_3 = 0.7 $, $ E_1 = 0.6$, $E_2 = 0.1$, $E_3 = 0.2$, $K_1 = 5$, $K_2 = 5$, $m_1 = 0.9$, $m_2 = 0.4 $.
    $ w_1 $ $ w_2 $ $ P_6 $ at $ \alpha=0 $ $ P_6 $ at $ \alpha=0.3 $ $ P_6 $ at $ \alpha=0.6 $ $ P_6 $ at $ \alpha=0.9 $
    0 1 $ (4.9826, 3.8455, 0) $ $ (4.6419, 3.5356, 0) $ $ (4.3385, 3.2568, 0) $ $ (4.0671, 3.0043, 0) $
    0.2 0.8 $ (4.5369, 3.4395, 0) $ $ (4.3577, 3.2745, 0) $ $ (4.1901, 3.1191, 0) $ $ (4.0331, 2.9724, 0) $
    0.4 0.6 $ (4.1543, 3.0858, 0) $ $ (4.1016, 3.0366, 0) $ $ (4.0500, 2.9883, 0) $ $ (3.9995, 2.9408, 0) $
    0.6 0.4 $ (3.8232, 2.7740, 0) $ $ (3.8700, 2.8184, 0) $ $ (3.9177, 2.8636, 0) $ $ (3.9665, 2.9097, 0) $
    0.8 0.2 $ (3.5351, 2.4958, 0) $ $ (3.6599, 2.6172, 0) $ $ (3.7926, 2.7447, 0) $ $ (3.9338, 2.8789, 0) $
    1 0 $ (3.2834, 2.2448, 0) $ $ (3.4690, 2.4307, 0) $ $ (3.6742, 2.6311, 0) $ $ (3.9017, 2.8485, 0) $

     | Show Table
    DownLoad: CSV
    Table 8.  The interior equilibrium $ P_7 $ for $ k_1 = 0.4$, $k_2 = 0.5$, $q_1 = 0.6$, $q_2 = 0.4$, $q_3 = 0.2 $, $ E_1 = 0.2$, $E_2 = 0.3$, $E_3 = 0.2$, $K_1 = 100$, $K_2 = 100$, $m_1 = 0.9$, $m_2 = 0.3 $.
    $ w_1 $ $ w_2 $ $ P_7 $ at $ \alpha=0 $ $ P_7 $ at $ \alpha=0.3 $
    0 1 $ (1.7240, 3.6276, 1.1986) $ $ (1.5910, 3.4597, 1.0271) $
    0.2 0.8 $ (1.6168, 3.5123, 1.0583) $ $ (1.5042, 3.3223, 0.9139) $
    0.4 0.6 $ (1.5417, 3.4296, 0.9337) $ $ (1.4211, 3.1793, 0.8023) $
    0.6 0.4 $ (1.4309, 3.2886, 0.7930) $ $ (1.3913, 3.1413, 0.7269) $
    0.8 0.2 $ (1.4291, 3.2451, 0.6884) $ $ (1.3195, 3.0467, 0.6350) $
    1 0 $ (1.3764, 2.9890, 0.5245) $ $ (1.2797, 2.9766, 0.5529) $
    $ w_1 $ $ w_2 $ $ P_7 $ at $ \alpha=0.6 $ $ P_7 $ at $ \alpha=0.9 $
    0 1 $ (1.7066, 3.7444, 1.0113) $ $ (1.5325, 3.4231, 0.7999) $
    0.2 0.8 $ (1.5291, 3.4203, 0.8678) $ $ (1.4899, 3.3773, 0.7762) $
    0.4 0.6 $ (1.5273, 3.4215, 0.8306) $ $ (1.4789, 3.3677, 0.7638) $
    0.6 0.4 $ (1.4894, 3.3737, 0.7779) $ $ (1.4685, 3.3587, 0.7516) $
    0.8 0.2 $ (1.4594, 3.3320, 0.7279) $ $ (1.4590, 3.3504, 0.7397) $
    1 0 $ (1.4328, 3.2552, 0.6687) $ $ (1.4124, 3.2841, 0.7105) $

     | Show Table
    DownLoad: CSV

    Considering four sets of different initial values, it could be seen from Figure 1 that different orbits eventually converged to the same value, which concluded that the interior equilibrium of the system (1.7) fulfills the character of globally asymptotical stability. Figure 2 plotted the bifurcation graph of system (1.7) with the horizontal coordinates $ k_1 $, and the Hopf bifurcation of the system occurred with $ k_1 $ taking values in the range of $ 0.01\leq k_1\leq0.7 $. When $ 0.01\leq k_1 < 0.384 $, the system oscillates periodically, while it maintains a stable steady-state when $ 0.384 < k_1\leq0.7 $. Therefore, based on Figure 2, it could be concluded that the fear of prey $ x_1 $ for predator y affected the stability of the system. We further observed that as $ k_1 $ increased, the prey $ x_1 $ density continued to decrease while the predator $ y $ density kept increasing. Thus, the result also suggested that greater fear of predators had a negative impact on prey populations while having a positive impact on predator populations. Correspondingly, Figures 3 and 4 showed the waveform plots and phase diagram at $ k_1 = 0.1 $ and $ k_1 = 0.7 $, respectively.

    Figure 1.  Global stability of the internal equilibrium $ P_7 $ = (5.665, 1.668, 2.047) of system (1.7) is given by the following parameter values: $ \alpha = 1$, $w_{1}+w_{2} = 1$, $A_1 = 2.0$, $A_2 = 2.0 $, $ B_1 = 2.0$, $B_2 = 2.0$, $k_1 = 0.2$, $k_2 = 0.1$, $q_1 = 0.4$, $q_2 = 0.4$, $q_3 = 0.2$, $E_1 = 0.2$, $E_2 = 0.2$, $E_3 = 0.2$, $A_3 = 0.1$, $B_3 = 0.1$, $A_4 = 0.3$, $B_4 = 0.6$, $K_1 = 10$, $K_2 = 10$, $m_1 = 0.4$, $m_2 = 0.4$, $C_1 = 0.1 $, $ C_2 = 0.2$, $C_3 = 0.5 $.
    Figure 2.  Hopf bifurcation occurs as a bifurcation parameter $ k_1 $, and the remaining parameters take the following values: $ \alpha = 1$, $w_{1}+w_{2} = 1 $, $ A_1 = 3.0$, $A_2 = 3.0$, $B_1 = 3.0$, $B_2 = 3.0$, $k_2 = 0.4$, $q_1 = 0.6$, $q_2 = 0.4$, $q_3 = 0.2$, $E_1 = 0.2$, $E_2 = 0.2$, $E_3 = 0.2$, $A_3 = 0.2$, $B_3 = 0.1$, $A_4 = 0.3$, $B_4 = 0.6$, $K_1 = 10$, $K_2 = 70$, $m_1 = 0.9$, $m_2 = 0.3$, $C_1 = 0.1$, $C_2 = 0.2$, $C_3 = 0.5 $.
    Figure 3.  Waveform plots and phase diagram of system (1.7) with $ k_1 = 0.1 $, and $ \alpha = 1$, $w_{1}+w_{2} = 1 $.
    Figure 4.  Waveform plots and phase diagram of system (1.7) with $ k_1 = 0.7 $, and $ \alpha = 1$, $w_{1}+w_{2} = 1 $.

    In addition, Figure 5 also plots the bifurcation graph with changing $ m_1 $. As can be seen in Figure 5, $ m_1 $ took values from 0.3 to 1, in which the system also underwent a Hopf bifurcation. When the value $ m_1 $ ranged from 0.3 to 0.657, the system (1.7) was stable; nevertheless, it would become unstable at $ 0.657 < m_1\leq 1 $. Correspondingly, Figures 6 and 7 showed the waveform plots and phase diagram at $ m_1 = 0.6 $ and $ m_1 = 0.9 $, respectively.

    Figure 5.  Hopf bifurcation occurs as a bifurcation parameter of system (1.7) parameter $ m_1 $, and the remaining parameters take the following values: $ \alpha = 1$, $w_{1}+w_{2} = 1 $, $ A_1 = 2.0$, $A_2 = 2.0$, $B_1 = 2.0$, $B_2 = 2.0$, $k_1 = 0.1$, $k_2 = 0.4$, $q_1 = 0.7$, $q_2 = 0.4$, $q_3 = 0.2$, $E_1 = 0.2$, $E_2 = 0.2$, $E_3 = 0.3$, $A_3 = 0.1$, $B_3 = 0.1$, $A_4 = 0.3$, $B_4 = 0.6$, $K_1 = 10$, $K_2 = 70$, $m_2 = 0.4$, $C_1 = 0.1$, $C_2 = 0.2$, $C_3 = 0.5 $.
    Figure 6.  Waveform plots and phase diagram of system (1.7) with $ m_1 = 0.6 $ and $ \alpha = 1$, $w_{1}+w_{2} = 1 $.
    Figure 7.  Waveform plots and phase diagram of system (1.7) with $ m_1 = 0.9 $ and $ \alpha = 1$, $w_{1}+w_{2} = 1 $.

    Further, we find an interesting dynamic phenomenon through some numerical simulations. System (1.7) appears as a chaotic phenomenon, as shown in Figure 8.

    Figure 8.  Waveform plots and phase diagram of chaotic phenomena with the following parameter values: $ \alpha = 1$, $w_{1}+w_{2} = 1 $, $ A_1 = 2.0$, $A_2 = 2.0$, $B_1 = 3.0$, $B_2 = 3.0$, $k_1 = 0.2$, $k_2 = 0.5$, $q_1 = 0.6$, $q_2 = 0.4$, $q_3 = 0.2$, $E_1 = 0.2$, $E_2 = 0.3$, $E_3 = 0.2$, $A_3 = 0.2$, $B_3 = 0.3$, $A_4 = 0.3$, $B_4 = 0.6$, $K_1 = 10$, $K_2 = 70$, $m_1 = 0.9$, $m_2 = 0.3$, $C_1 = 0.1$, $C_2 = 0.2$, $C_3 = 0.5 $.

    In this work, we develop a model of one-predator and two-prey interactions in a fuzzy environment, examine the effects of fear and prey refuge on the system, and provide insight into the dynamic complexity. The proofs of the theoretical parts of this paper are based on system (1.7). It has been proven that all equilibria in system (1.7) are locally asymptotically stable, and interior equilibrium $ P_7 $ is also globally asymptotically stable. We have been further concerned about the appearance and direction of Hopf bifurcation. With the support of theoretical research, our numerical simulations have been able to display a wealth of charts and graphs.

    First of all, different equilibria are displayed from Tables 28 with different $ \alpha, w_1, w_2 $, respectively. Throughout Figure 1, we have verified the global asymptotical stability of interior equilibrium $ P_7 $, and find that the system is from unstable to stable with the increase of fear $ k_1 $, which demonstrates that the fear effect may be an important factor influencing the stability of the system (see Figures 24). Furthermore, it has also been observed that an increase in prey refuge $ m_1 $ leads to oscillatory phenomena (see Figures 57). Finally, through studying the Hopf bifurcation, we have discovered some interesting biological phenomena, namely that system (1.7) appears to be in a chaotic state (see Figure 8).

    Xuyang Cao: Conceptualization, Investigation, Methodology, Validation, Writing-original draft, Formal analysis, Software; Qinglong Wang: Conceptualization, Methodology, Formal analysis, Writing-review and editing, Supervision; Jie Liu: Validation, Visualization, Data curation.

    The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

    The authors thank the editor and referees for their careful reading and valuable comments.

    The work is supported by the Natural Science Foundation of Hubei Province (No. 2023AFB1095) and the National Natural Science Foundation of China (No. 12101211) and the Program for Innovative Research Team of the Higher Education Institution of Hubei Province (No. T201812) and the Teaching Research Project of Education Department of Hubei Province (No. 2022367) and the Graduate Education Innovation Project of Hubei Minzu University (Nos. MYK2024071, MYK2023042).

    The authors declare that there are no conflicts of interest regarding the publication of this paper.

    Definition 1. [34] Fuzzy set: A fuzzy set $ \widetilde{\hbar} $ in a universe of discourse $ S $ is denoted by the set of pairs

    $ \widetilde{\hbar} = \{(s,\mu_{\widetilde{\hbar}}(s)):s\in S\}, $

    where the mapping $ \mu_{\widetilde{\hbar}}: S\rightarrow [0, 1] $ is the membership function of the fuzzy set $ \widetilde{\hbar} $ and $ \mu_{\widetilde{\hbar}} $ is the membership value or degree of membership of $ s\in S $ in the fuzzy set $ \widetilde{\hbar} $.

    Definition 2. [42] $ \alpha $-cut of fuzzy set: For any $ \alpha\in(0, 1] $, the $ \alpha $-cut of fuzzy set $ \widetilde{\hbar} $ defined by $ \hbar_{\alpha} = \{s:\mu_{\widetilde{\hbar}}(s))\geq \alpha\} $ is a crisp set. For $ \alpha = 0 $ the support of $ \widetilde{\hbar} $ is defined as $ \hbar_{0} = {\rm Supp}(\widetilde{\hbar}) = \overline{\{s\in \mathbb{R}, \mu_{\widetilde{\hbar}}(s) > 0\}} $.

    Definition 3. [43] Fuzzy number: A fuzzy number satisfying the property $ S = \mathbb{R} $ is called a convex fuzzy set.

    Definition 4. [44] Triangular fuzzy number: A triangular fuzzy number (TFN) $ \widetilde{\hbar}\equiv (b_{1}, b_{2}, b_{3}) $ represent fuzzy set of the real line $ \mathbb{R} $ satisfying the property that the membership function $ \mu_{\widetilde{\hbar}}:\mathbb{R} $$ \rightarrow [0, 1] $ can be espressed by

    $ \mu_{\widetilde{\hbar}} = \left\{sb1b2b1  if  b1sb2,b3sb3b2  if  b2sb3,0otherwise.
    \right. $

    Hence, the $ \alpha $-cut of triangular fuzzy number meets boundedness and encapsulation on $ [\hbar_{L}(\alpha), \hbar_{R}(\alpha)] $, in which $ \hbar_{L}(\alpha) = {\rm inf} {s:\mu_{\widetilde{\hbar}}(s)\geq \alpha} = b_{1}+\alpha(b_{2}-b_{1}) $ and $ \hbar_{R}(\alpha) = {\rm sup} \{s:\mu_{\widetilde{\hbar}}(s)\geq \alpha\} = b_{3}+\alpha(b_{3}-b_{2}) $.

    Lemma 1. [45] In weighted sum method, $ w_{j} $ stands for the weight of $ j $th objective. $ w_{j}g_{j} $ represent a utility function for $ j $th objective, and the total utility function $ \pi $ is represented by

    $ \pi = \sum\limits_{j}^l w_{j}g_{j},\; \; \; j = 1,2,\cdots,l, $

    where $ w_{j} > 0 $ and $ \sum^{l}_{j} w_{j} = 1 $ are satisfied.

    [1] (2000) Functions of bounded variation and free discontinuity problems. Oxford Univ. Press.
    [2]

    J. M. Ball, Some open problems in elasticity, In Geometry, Mechanics and Dynamics, Springer, New York, 2002, 3-59

    [3] From molecular models to continuum mechanics. Arch. Ration. Mech. Anal. (2002) 164: 341-381.
    [4] (2002) $Γ$-Convergence for Beginners. Oxford Univ. Press.
    [5] Surface energies in nonconvex discrete systems. Math. Models Methods Appl. Sci. (2007) 17: 985-1037.
    [6] Variational formulation of softening phenomena in fracture mechanics: the one-dimensional case. Arch. Rational Mech. Anal. (1999) 146: 23-58.
    [7] Continuum limits of discrete systems without convexity hypotheses. Math. Mech. Solids (2002) 7: 41-66.
    [8]

    A. Braides and M. S. Gelli, The passage from discrete to continuous variational problems: A nonlinear homogenization process, in Nonlinear Homogenization and its Applications to Composites, Polycrystals and Smart Materials, NATO Sci. Ser. II Math. Phys. Chem., Kluwer Acad. Publ., Dordrecht, 170 (2004), 45-63.

    [9] From discrete systems to continuous variational problems: An introduction. Lect. Notes Unione Mat. Ital. (2006) 2: 3-77.
    [10] Asymptotic analysis of microscopic impenetrability constraints for atomistic systems. J. Mech. Phys. Solids (2016) 96: 235-251.
    [11] The passage from nonconvex discrete systems to variational problems in Sobolev spaces: The one-dimensional case. Proc. Steklov Inst. Math. (2002) 236: 395-414.
    [12] Effective cohesive behavior of layers of interatomic planes. Arch. Rational Mech. Anal. (2006) 180: 151-182.
    [13] Asymptotic analysis of Lennard-Jones systems beyond the nearest-neighbour setting: A one-dimensional prototypical case. Math. Mech. Solids (2016) 21: 915-930.
    [14] Asymptotic expansions by $Γ$-convergence. Cont. Mech. Thermodyn. (2008) 20: 21-62.
    [15] An atomistic-to-continuum analysis of crystal cleavage in a two-dimensional model problem. J. Nonlinear Sci. (2014) 24: 145-183.
    [16] An analysis of crystal cleavage in the passage from atomistic models to continuum theory. Arch. Ration. Mech. Anal. (2015) 217: 263-308.
    [17] Asymptotic behaviour of a pile-up of infinite walls of edge dislocations. Arch. Ration. Mech. Anal. (2013) 209: 495-539.
    [18] Atomistic-to-continuum coupling. Acta Numer. (2013) 22: 397-508.
    [19] Boundary layer energies for nonconvex discrete systems. Math. Models Methods Appl. Sci. (2011) 21: 777-817.
    [20] Towards uniformly $Γ$-equivalent theories for nonconvex discrete systems. Discrete Contin. Dyn. Syst. Ser. B (2012) 17: 661-686.
    [21] On a $Γ$-convergence analysis of a quasicontinuum method. Multiscale Model. Simul. (2015) 13: 132-172.
    [22]

    M. Schäffner, Multiscale analysis of non-convex discrete systems via $Γ$-convergence, Ph. D thesis, University of Würzburg, 2015.

    [23] Quasicontinuum analysis of defects in solids. Phil. Mag. A (2006) 73: 1529-1563.
    [24] Fracture as a phase transition. Contemporary Research in the Mechanics and Mathematics of Marterials (1996) 322-332.
  • This article has been cited by:

    1. Yuan Tian, Hua Guo, Wenyu Shen, Xinrui Yan, Jie Zheng, Kaibiao Sun, Dynamic analysis and validation of a prey-predator model based on fish harvesting and discontinuous prey refuge effect in uncertain environments, 2025, 33, 2688-1594, 973, 10.3934/era.2025044
  • Reader Comments
  • © 2018 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(5270) PDF downloads(309) Cited by(6)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog