
Citation: Enrico Troiani. Analytical evaluation of the Stress Intensity Factor in stiffened sheets with multiple side damage[J]. AIMS Materials Science, 2016, 3(4): 1615-1622. doi: 10.3934/matersci.2016.4.1615
[1] | Sze-Bi Hsu, Feng-Bin Wang, Xiao-Qiang Zhao . Mathematical modeling and analysis of harmful algal blooms in flowing habitats. Mathematical Biosciences and Engineering, 2019, 16(6): 6728-6752. doi: 10.3934/mbe.2019336 |
[2] | Konstantin E. Starkov, Giovana Andres Garfias . Dynamics of the tumor-immune-virus interactions: Convergence conditions to tumor-only or tumor-free equilibrium points. Mathematical Biosciences and Engineering, 2019, 16(1): 421-437. doi: 10.3934/mbe.2019020 |
[3] | Lin Zhang, Yongbin Ge, Xiaojia Yang . High-accuracy positivity-preserving numerical method for Keller-Segel model. Mathematical Biosciences and Engineering, 2023, 20(5): 8601-8631. doi: 10.3934/mbe.2023378 |
[4] | Guo Lin, Shuxia Pan, Xiang-Ping Yan . Spreading speeds of epidemic models with nonlocal delays. Mathematical Biosciences and Engineering, 2019, 16(6): 7562-7588. doi: 10.3934/mbe.2019380 |
[5] | Jizhi Huang, Guoyuan Xu, Yu Wang, Xiaowei Ouyang . Equivalent deformation modulus of sandy pebble soil—Mathematical derivation and numerical simulation. Mathematical Biosciences and Engineering, 2019, 16(4): 2756-2774. doi: 10.3934/mbe.2019137 |
[6] | Huanan Yu, Yutang Gao, Guoping Qian, Chao Zhang, Changyun Shi, Jinguo Ge, Wan Dai . Study on synergistic effect of multiple physical fields on hot mix asphalt during compaction process. Mathematical Biosciences and Engineering, 2024, 21(4): 5181-5206. doi: 10.3934/mbe.2024229 |
[7] | Rong Liu, Fengqin Zhang . A hierarchical age-structured model of optimal vermin management by contraception. Mathematical Biosciences and Engineering, 2023, 20(4): 6691-6720. doi: 10.3934/mbe.2023288 |
[8] | Xiawei Yang, Wenya Li, Yaxin Xu, Xiurong Dong, Kaiwei Hu, Yangfan Zou . Performance of two different constitutive models and microstructural evolution of GH4169 superalloy. Mathematical Biosciences and Engineering, 2019, 16(2): 1034-1055. doi: 10.3934/mbe.2019049 |
[9] | K. E. Starkov, Svetlana Bunimovich-Mendrazitsky . Dynamical properties and tumor clearance conditions for a nine-dimensional model of bladder cancer immunotherapy. Mathematical Biosciences and Engineering, 2016, 13(5): 1059-1075. doi: 10.3934/mbe.2016030 |
[10] | Edwin Duque-Marín, Alejandro Rojas-Palma, Marcos Carrasco-Benavides . A soil water indicator for a dynamic model of crop and soil water interaction. Mathematical Biosciences and Engineering, 2023, 20(8): 13881-13899. doi: 10.3934/mbe.2023618 |
The dynamical analysis of complex ecosystem, such as food chain, is based on the interaction among species or between two species, especially the dynamical relationship between predator and prey [1,2,3,4]. The current theory of predator-prey dynamics must depend on the study of nonlinear mathematical model [5]. With the continuous improvement of ecological knowledge like theoretical research, empirical research and observational research, etc., there are more and more basic elements of predation to be considered. Therefore, modelers add some complexities to their abstraction in order to obtain authenticity from the emergence of the far-reaching Lotka-Volterra model[6] and the modifications introduced by Volterra[5], taking into account the self-interference of prey populations.
Allee effect, which affects the number of prey, is one of these factors. It changes the qualitative stability and quantitative aspects of predator-prey model dynamics. Because the interaction between predator and prey is naturally prone to vibration, it is obvious to study this phenomenon as a potential mechanism for the generation of population cycles. Lots of researches about predator-prey models are done with the Allee effect[7,8,9].
The most popular framework for modeling expert predator-prey interaction has the following structure:
{dxdt=xg(x,k)−yh(x),dydt=(ph(x)−c)y, | (1.1) |
where x(t) and y(t) are the prey and predator population sizes in the time t, respectively, p,c>0 indicate the birth rate and background mortality rate, respectively, g(x,k) describes the specific growth rate of the prey in the absence of predator, and h(x) describes the predator functional response.
Any mechanism leading to a positive correlation between the components of individual fitness and the number or density of similar individuals can be regarded as Allee effect; it describes a scenario in which low population size is affected by the positive correlation between population growth rate and density, increasing the possibility of their extinction.
Recent ecological studies have shown that two or more Allee effects can lead to mechanisms acting on a population at the same time. The combined effects of some of these phenomena are called multiple (double) Allee effects. The author's analysis in[10] showed that the results of strong and weak Allee effects on the dynamics of Volterra predator-prey model are similar, which originate from the limit cycle of the model.
In this paper, we continue to consider the following predator-prey model with double Allee effect functional response raised by[10].
{dXdT=rX(1−XK)X−MX+N−qXY,dYdT=(pX−C)Y, | (1.2) |
where r scales the prey growth rate, K is the intrinsic carrying capacity for the environment to the prey in the absence of predator, M is the Allee threshold, and the auxiliary parameter N satisfies N>0, q is the prey captured rate by the predator, p,C>0 indicate the birth rate and background mortality rate, respectively.
One can see that the first equation in the model (1.2) includes double Allee effects, expressed by the first factor m(X)=X−M, and a second term r(X)=rXX+N, This can be interpreted as an approximation of population dynamics, in which the difference between fertile and non-fertile are not clearly modeled.
In order to simplify the analysis of system (1.2), we make a topologically equivalent change of variables and time rescaling as in[11,12,13,14], defining the function ϕ, such that ϕ(x,y)=(Kx,rqy)=(X,Y), rx+ndT=dt. Then, system (1.2) is transformed into
{dxdt=((1−x)(x−m)−(x+n)y)x,dydt=b(x−c)(x+n)y, | (1.3) |
where b=pKr, CpK, m=MK, n=NK, and K>X>M is obtained from equation (1.2), so 1>m>0.
We now use the semidiscretization method, which has been applied in many studies[15,16,17,18], to study the discrete model of system (1.3). For this, suppose that [t] denotes the greatest integer not exceeding t. Consider the following semidiscretization version of (1.3).
{1x(t)dx(t)dt=(1−x([t]))(x([t])−m)−(x([t])+n)y([t]),1y(t)dy(t)dt=b(x([t])−c)(x([t])+n). | (1.4) |
It is easy to see that the system (1.4) has piecewise constant arguments, and that a solution (x(t),y(t)) of the system (1.4) for t∈[0,+∞) possesses the following natures:
1. on the interval [0,+∞), x(t) and y(t) are continuous;
2. when t∈[0,+∞) except for the points t∈{0,1,2,3,⋯}, dx(t)dt and dy(t)dt exist everywhere.
The following system can be obtained by integrating (1.4) over the interval [n,t] for any t∈[n,n+1) and n=0,1,2,⋯
{x(t)=xne(1−xn)(xn−m)−(xn+n)yn(t−n),y(t)=yneb(xn−c)(xn+n)(t−n), | (1.5) |
where xn=x(n) and yn=y(n). Letting t→(n+1)− in (1.5) produces
{xn+1=xne(1−xn)(xn−m)−(xn+n)yn,yn+1=yneb(xn−c)(xn+n), | (1.6) |
where b,c,n>0, 1>m>0 are the same as in (1.3). We mainly study the properties of system (1.6) in the sequel.
The rest of the paper is organized as follows: In Section 2, we investigate the existence and stability of fixed points of the system (1.6). In Section 3, we derive the sufficient conditions for the transcritical bifurcation and the Neimark-Sacker bifurcation of the system (1.6) to occur. In Section 4, numerical simulations are performed to verify the above obtained theoretical results and reveal some new dynamical properties.
Before we analyze the fixed points of the system (1.6), we recall the following lemma see [16,19].
Lemma 1.1. Let F(λ)=λ2+Pλ+Q, where P and Q are two real constants. Suppose λ1 and λ2 are two roots of F(λ)=0. Then the following statements hold.
(ⅰ) If F(1)>0, then
(ⅰ.1) |λ1|<1 and |λ2|<1 if and only if F(−1)>0 and Q<1;
(ⅰ.2) λ1=−1 and λ2≠−1 if and only if F(−1)=0 and P≠2;
(ⅰ.3) |λ1|<1 and |λ2|>1 if and only if F(−1)<0;
(ⅰ.4) |λ1|>1 and |λ2|>1 if and only if F(−1)>0 and Q>1;
(ⅰ.5) λ1 and λ2 are a pair of conjugate complex roots with |λ1|=|λ2|=1 if and only if −2<P<2 and Q=1;
(ⅰ.6) λ1=λ2=−1 if and only if F(−1)=0 and P=2.
(ii) If F(1)=0, namely, 1 is one root of F(λ)=0, then the another root λ satisfies |λ|=(<,>)1 if and only if |Q|=(<,>)1.
(iii) If F(1)<0, then F(λ)=0 has one root lying in (1,∞). Moreover,
(iii.1) the other root λ satisfies λ<(=)−1 if and only if F(−1)<(=)0;
(iii.2) the other root −1<λ<1 if and only if F(−1)>0.
In this section, we first consider the existence of fixed points and then analyze the local stability of each fixed point.
The fixed points of the system (1.6) satisfy
x=xe(1−x)(x−m)−(x+n)y, |
y=yeb(x−c)(x+n). |
Considering the biological meanings of the system (1.6), one only takes into account nonnegative fixed points. Thereout, one finds that the system (1.6) has and only has four nonnegative fixed points E0=(0,0), E1=(1,0), E2=(m,0) and E3=(x0,y0) for m<c<1, where
x0=c,y0=(1−c)(c−m)c+n. |
The Jacobian matrix of the system (1.6) at any fixed point E(x,y) takes the following form
J(E)=([1+x(−2x−y+m+1)]eA−x(x+n)eAb(2x+n−c)yeb(x−c)(x+n)eb(x−c)(x+n)). |
where A=e(1−x)(x−m)−(x+n)y.
The characteristic polynomial of Jacobian matrix J(E) reads
F(λ)=λ2−pλ+q, |
where
p=Tr(J(E))=[1+x(−2x−y+m+1)]eA+eb(x−c)(x+n), q=Det(J(E))=[1+x(−2x−y+m+1)][bx(x+n)(2x+n−c)]eA+b(x−c)(x+n).
For the stability of fixed points E0, E1, E2 and E3, we can easily get the following Theorems 2.1-2.4 respectively.
Theorem 2.1. The fixed point E0=(0,0) of the system (1.6) is a sink.
Proof. The Jacobian matrix J(E0) of the system (1.6) at the fixed point E0=(0,0) is given by
J(E0)=(e−m00e−bcn). |
Obviously, |λ1|=e−m<1 and |λ2|=e−bcn<1, so E0=(0,0) is a sink.
Theorem 2.2. The following statements about the fixed point E1=(1,0) of the system (1.6) are true.
1. If c<1, then E1 is a saddle.
2. If c=1, then E1 is non-hyperbolic.
3. If c>1, then E1 is a stable node.
Proof. The Jacobian matrix of the system (1.6) at E1=(1,0) is
J(E1)=(m−(1+n)0eb(1−c)(1+n)). |
Obviously, λ1=m and λ2=eb(1−c)(1+n).
Note |λ1|<1 is always true. If c<1, then |λ2|>1, so E1 is a saddle; if c=1, then |λ2|=1, therefore E1 is non-hyperbolic; if c>1, implying |λ2|<1, then E1 is a stable node, namely, a sink. The proof is complete.
Theorem 2.3. The following statements about the fixed point E2=(m,0) of the system (1.6) are true.
1. If c<m, then E2 is a source.
2. If c=m, then E2 is non-hyperbolic.
3. If c>m, then E2 is a saddle.
Proof. The Jacobian matrix of the system (1.6) at E2=(m,0) is
J(E2)=(−m2+m+1−m(m+n)0eb(m−c)(m+n)). |
Obviously, λ1=−m2+m+1 and λ2=eb(m−c)(m+n).
Note 0<m<1, so |λ1|>1 is always true. If c<m, then |λ2|>1, so E2 is a source; if c=m, then |λ2|=1, therefore E2 is non-hyperbolic; if c>m, implying |λ2|<1, then E2 is a saddle. The proof is finished.
Theorem 2.4. When (1−c)(c−m)>0, namely, 0<m<c<1, the fixed point E3=(c,(1−c)(c−m)c+n) is a positive fixed point of the system (1.6). Let b0=c2+2cn−mn−m−n(c+n)2(1−c)(c−m). Then the following statements are true about the positive fixed point E3.
Ⅰ. If m<c2+2cn−nn+1, then,
1. for 0<b<b0, E3 is a stable node;
2. for b=b0, E3 is non-hyperbolic;
3. for b>b0, E3 is an unstable node.
Ⅱ. If m≥c2+2cn−nn+1, then, E3 is an unstable node.
Proof. The Jacobian matrix of the system (1.6) at E3 can be simplified as follows
J(E3)=(1+−c(c2+2cn−mn−m−n)c+n−c(c+n)b(1−c)(c−m)1). |
The characteristic polynomial of Jacobian matrix J(E3) reads as
F(λ)=λ2−pλ+q, | (2.1) |
where
p=Tr(J(E3))=2−b0c(1−c)(c−m)(c+n),
q=Det(J(E3))=1+(b−b0)c(1−c)(c−m)(c+n). By calculating we get
F(1)=bc(1−c)(c−m)(c+n)>0,
and
F(−1)=2(2−c2)+2c(n(1−c)+mn+n)c+n+bc(1−c)(c−m)(c+n)>0. |
Ⅰ. If m<c2+2cn−nn+1, then b0>0. So, when 0<b<b0, q<1. By Lemma 1.1 (i.1), |λ1|<1 and |λ2|<1, therefore E3 is a stable node, i.e., a sink. When b=b0, q=1, −2<p<2. By Lemma 1.1 (i.5), Eq. (2.1) has a pair of conjugate complex roots λ1 and λ2 with |λ1|=|λ2|=1, implying E3 is non-hyperbolic. When b>b0<, q>1. Lemma 1.1 (i.4) tells us that |λ1|>1 and |λ2|>1, so E3 is an unstable node, i.e., a source.
Ⅱ. If m≥c2+2cn−nn+1, then b0≤0. So, b>0≥b0. Hence q>1. Lemma 1.1 (i.4) reads that E3 is an unstable node. The proof is complete.
In this section, we are in a position to use the Center Manifold Theorem and bifurcation theorem to analyze the local bifurcation problems of the fixed points E1, E2 and E3. For related work, refer to[20,21,22,23,24,25].
Theorem 2.2 shows that a bifurcation of E1 may occur in the space of parameters (b,c,m,n)∈SE+={(b,c,m,n)∈R4+|b>0,c>0,1>m>0,n>0}.
Theorem 3.1. Set the parameters (b,c,m,n)∈SE+={(b,c,m,n)∈R4+|b>0,c>0,1>m>0,n>0}. Let c0=1, then the system (1.6) undergoes a transcritical bifurcation at E1 when the parameter c varies in a small neighborhood of c0.
Proof. In order to show the detailed process, we proceed according to the following steps.
The first step. Let un=xn−1,vn=yn−0, which transforms the fixed point E1=(1,0) to the origin O(0,0), and the system (1.6) to
{un+1=(un+1)e−un(un−m+1)−(un+n+1)vn−1,vn+1=vneb(un−c+1)(un+n+1). | (3.1) |
The second step. Giving a small perturbation c∗ of the parameter c, i.e., c∗=c−c0, with 0<|c∗|≪1, the system (3.1) is perturbed into
{un+1=(un+1)e−un(un−m+1)−(un+n+1)vn−1,vn+1=vneb(un−c∗)(un+n+1). | (3.2) |
Letting c∗n+1=c∗n=c∗, the system (3.2) can be written as
{un+1=(un+1)e−un(un−m+1)−(un+n+1)vn−1,vn+1=vneb(un−c∗n)(un+n+1),c∗n+1=c∗n. | (3.3) |
The third step. Taylor expanding of the system (3.3) at (un,vn,c∗n)=(0,0,0) takes the form
{un+1=a100un+a010vn+a200u2n+a020v2n+a110unvn+a300u3n+a030v3n+a210u2nvn+a120unv2n+o(ρ31),vn+1=b100un+b010vn+b001c∗n+b200u2n+b020v2n+b002c∗n2+b110unvn+b101unc∗n+b011vnc∗n+b300u3n+b030v3n+b003c∗n3+b210u2nvn+b120unv2n+b201u2nc∗n+b102unc∗n2+b021v2nc∗n+b012vnc∗n2+b111unvnc∗n+o(ρ31),c∗n+1=c∗n, | (3.4) |
where ρ1=√u2n+v2n+(c∗n)2,
a100=m,a010=−(n+1),a200=12(m−1)2+m−2,a020=12(n+1)2,a110=−mn−m−1,a300=16(m−1)3+12(m−1)2−m,a030=−16(n+1)3,a120=12m(n+1)2+n+1,a210=−12(m−1)2(n+1)−mn−2m−2n+3, |
b100=b001=b200=b020=b002=b101=b300=b030=b003=b120=b201=b102=b021=0,b010=1,b110=b(n+1),b011=−b(n+1),b210=b+12b2(n+1)2,b012=12b2(n+1)2,b111=−b−b2(n+2)2. |
Let
J(E1)=(a100a0100b100b0100001) |
i.e.,J(E1)=(m−(n+1)0010001). Then, we derive the three eigenvalues of J(E1) to be
λ1=m,λ2,3=1, |
and the corresponding eigenvectors
(ξ1,η1,φ1)T=(1,0,0)T,(ξ2,η2,φ2)T=(n+1,m−1,0)T,(ξ3,η3,φ3)T=(0,0,1)T. |
The fourth step. Let T=(ξ1ξ2ξ3η1η2η3φ1φ2φ3), namely,
T=(1n+100m−10001), |
then T−1=(11+n1−m001m−10001). |
Taking the following transformation
(un,vn,c∗n)T=T(Xn,Yn,δn)T, |
the system (3.4) is changed into the following form
{Xn+1=mXn+F(Xn,Yn,δn)+o(ρ32),Yn+1=Yn+G(Xn,Yn,δn)+o(ρ32),δn+1=δn, | (3.5) |
where ρ2=√X2n+Y2n+δ2n,
F(Xn,Yn,δn)=m200X2n+m020Y2n+m002δn2+m110XnYn+m101Xnδn+m011Ynδn+m300X3n+m030Y3n+m003δn3+m210X2nYn+m120XnY2n+m201X2nδn+m102Xnδn2+m021Y2nδn+m012Ynδn2+m111XnYnδn, |
G(Xn,Yn,δn)=l200X2n+l020Y2n+l002δn2+l110XnYn+l101Xnδn+l011Ynδn+l300X3n+l030Y3n+l003δn3+l210X2nYn+l120XnY2n+l201X2nδn+l102Xnδn2+l021Y2nδn+l012Ynδn2+l111XnYnδn, |
m200=a200,m300=a300,m002=m101=m003=m201=m102=0,m020=(a200−b110)(1+n)2+a020(m−1)2+a110(m−1)(n+1),m110=(2a200−b110)(1+n)+a110(m−1),m011=−b011(1+n),m030=(a300−b210)(1+n)3+a030(m−1)3+a210(1+n)2(m−1)+a120(1+n)(m−1)2,m210=(3a300−b210)(1+n)+a210(m−1),m120=(3a300−2b210)(1+n)2+2a210(1+n)(m−1)+a120(m−1)2,m021=−b111(1+n)2,m012=−b012(1+n),m111=−b111(1+n),l200=l002=0,l020=b110(1+n),l110=b110,l101=0,l011=b011,l300=0,l030=b210(1+n)2,l003=0,l210=b210,l120=2b210(1+n),l201=l102=0,l021=b111(1+n),l012=b012,l111=b111. |
The fifth step. Suppose on the center manifold
Xn=h(Yn,δn)=h20Y2n+h11Ynδn+h02δ2n+o(ρ23), |
where ρ3=√Y2n+δ2n, then, according to
Xn+1=mh(Yn,δn)+F(h(Yn,δn),Yn,δn)+o(ρ23), |
h(Yn+1,δn+1)=h20Y2n+1+h11Yn+1δn+1+h02δ2n+1+o(ρ23)=h20(Yn+G(h(Yn,δn),Yn,δn))2+h11(Yn+G(h(Yn,δn),Yn,δn))δn+h02δ2n+o(ρ23). |
and Xn+1=h(Yn+1,δn+1), we obtain the center manifold equation to satisfy the following relation
mh(Yn,δn)+F(h(Yn,δn),Yn,δn)=h20(Yn+G(h(Yn,δn),Yn,δn))2+h11(Yn+G(h(Yn,δn),Yn,δn))δn+h02δ2n+o(ρ23).
Comparing the corresponding coefficients of terms with the same orders in the above center manifold equation, we get
h20=(m−2)(1+n)2−(1+n)(2+n)(m−1)−b(1+n)31−m,
h11=b(1+n)21−m,h02=0.
So, the system (3.5) restricted to the center manifold takes as
Yn+1=f(Yn,δn):=Yn+G(h(Yn,δn),Yn,δn)+o(ρ23)=Yn+b(1+n)2Y2n−b(1+n)Ynδn+b(1+n)21−m(1−2m−n−12b(1+n)2(1+m))Y3n+(mb2(1+n)31−m−b(1+n))Yn2δn+12b2(1+n)2Ynδ2n+o(ρ33).
Therefore one has
f(Yn,δn)|(0,0)=0,∂f∂Yn|(0,0)=1,∂f∂δn|(0,0)=0,∂2f∂Yn∂δn|(0,0)=−b(1+n)≠0,∂2f∂Y2n|(0,0)=2b(1+n)2≠0. |
According to (21.1.42)-(21.1.46) in the literature ([26], pp. 507), all the conditions for the occurrence of the transcritical bifurcation are established, hence, it is valid for the {occurrence} of transcritical bifurcation in the fixed point E1. The proof is over.
According to Theorem 2.3, the fixed point E2(m,0) is non-hyperbolic, the system (1.6) may undergo a bifurcation (the correspond eigenvalue are λ1=−m2+m+1, λ2=1). By using the same method as that in Section 3.2, we get the following result.
Theorem 3.2. Set the parameters (b,c,m,n)∈SE+={(b,c,m,n)∈R4+|b>0,c>0,1>m>0,n>0}. Let c1=m, then the system (1.6) undergoes a transcritical bifurcation at E2 when the parameter c varies in a small neighborhood of c1.
When m<c2+2cn−nn+1, b=b0=c2+2cn−mn−n−m(c+n)2(1−c)(c−m), Theorem 2.4 with Lemma 1.2 (i.5) shows that F(1)>0, F(−1)>0, −2<p<2 and q=1, so λ1 and λ2 are a pair of conjugate complex roots with |λ1|=|λ2|=1. At this time we derive that the system (1.6) at the fixed point E3 can undergo a Neimark-Sacker bifurcation in the space of parameters (b,c,m,n)∈SE+={(b,c,m,n)∈R4+|b>0,1>c>m,0<m<c2+2cn−nn+1}.
In order to show the process clearly, we carry out the following steps.
The first step. Take the changes of variables un=xn−x0,vn=yn−y0, which transform fixed point E3=(x0,y0) to the origin O(0,0), and the system (1.6) into
{un+1=(un+x0)e(−un−x0+1)(un+x0−m)−(un+x0+n)vnyn−x0,vn+1=(vn+y0)eb(un+x0−c)(un+x0+n)−y0. | (3.6) |
The second step. Give a small perturbation b∗ of the parameter b, i.e., b∗=b−b0, then the perturbation of the system (3.6) can be regarded as follows
{un+1=(un+x0)e(−un−x0+1)(un+x0−m)−(un+x0+n)vnyn−x0,vn+1=(vn+y0)e(b∗+b0)(un+x0−c)(un+x0+n)−y0. | (3.7) |
The corresponding characteristic equation of the linearized equation of the system (3.7) at the equilibrium point (0, 0) can be expressed as
F(λ)=λ2−p(b∗)λ+q(b∗)=0, |
where
p(b∗)=2+c(−c2−2cn+mn+m+n)c+n, |
and
q(b∗)=c[(c2+2cn−mn−m−n)−(b∗+b0)(c+n)2(1−c)(c−m)]−(c+n)+1.
It is easy to derive p2(b∗)−4q(b∗)<0 when b∗=0, and 0<p(b∗)<2, then the two roots of F(λ)=0 are
λ1,2(b∗)=p(b∗)±√p2(b∗)−4q(b∗)2=p(b∗)±i√4q(b∗)−p2(b∗)2, |
which implies
(|λ1,2(b∗)|)|b∗=0=√q(b∗)|b∗=0=1, |
and
(d|λ1,2(b∗)|db∗)|b∗=0=12c(c+n)(1−c)(c−m)>0. |
The occurrence of Neimark-Sacker bifurcation requires the following conditions to be satisfied
(H.1)(d|λ1,2(b∗)|db∗)|b∗=0≠0; |
(H.2)λi1,2(0)≠1,i=1,2,3,4. |
Since p(b∗)|b∗=0=2+c(−c2−2cn+mn+m+n)c+n and q(b∗)|b∗=0=1, we have λ1,2(0)=2(c+n)+c(−c2−2cn+mn+m+n)±i√4(c+n)2−[2(c+n)+c(−c2−2cn+mn+m+n)]22(c+n), then it is easy to derive λi1,2(0)≠1 for all i=1,2,3,4. According to ([2], pp517-522), they satisfy all of the conditions for Neimark-Sacker bifurcation to occur.
The third step. In order to derive the normal form of the system (3.7), we expand the system (3.7) into power series up to the following third-order form around the origin
{un+1=c10un+c01vn+c20u2n+c11unvn+c02v2n+c30u3n+c21u2nvn+c12unv2n+c03v3n+o(ρ34),vn+1=d10un+d01vn+d20u2n+d11unvn+d02v2n+d30u3n+d21u2nvn+d12unv2n+d03v3n+o(ρ34), | (3.8) |
where ρ4=√u2n+v2n,
c10=1+c(−c2−2cn+mn+m+n)c+n,c01=−c(c+n),c20=−c+−c2−2cn+mn+m+nc+n+c(−c2−2cn+mn+m+n)22(c+n)2,c02=c(c+n)22,c11=−[2c+n+c(−c2−2cn+mn+m+n)],c30=−1−c(−c2−2cn+mn+m+n)c+n+(−c2−2cn+mn+m+n)2c(c+n)2+c(−c2−2cn+mn+m+n)36(c+n)3,c03=−c(c+n)36,c21=−1−(−c2−2cn+mn+m+n)+c[c+n−−c2−2cn+mn+m+nc+n−(−c2−2cn+mn+m+n)22(c+n)],c12=(c+n)22+c(c+n)[1+−c2−2cn+mn+m+n2],d10=−−c2−2cn+mn+m+n(c+n)2,d01=1,d20=(−c2−2cn+mn+m+n)22(c+n)3(1−c)(c−m)−2(−c2−2cn+mn+m+n)(1−c)(c−m)2(c+n)3(1−c)(c−m),d11=−−c2−2cn+mn+m+n(c+n)(1−c)(c−m),d02=d03=d12=0,d30=(−c2−2cn+mn+m+n)2(c+n)4(1−c)(c−m)[1−−c2−2cn+mn+m+n6(1−c)(c−m)],d21=(−c2−2cn+mn+m+n)(c+n)2(1−c)(c−m)[1−−c2−2cn+mn+m+n2(1−c)(c−m)]. |
Let
J(E3)=(c10c01d10d01),namely, |
J(E3)=(1+K−c(c+n)−Kc(c+n)1). |
It is easy to derive the two eigenvalues of the matrix J(E3) are
λ1,2=(1+12K)±βi, |
where K=c(−c2−2cn+mn+m+n)c+n,
β=√−c(−c2−2cn+mn+m+n)[4(c+n)+c(−c2−2cn+mn+m+n)]2(c+n),
with the corresponding eigenvectors
v1,2=(−c(c+n)−12K)±i(0β).
Let
T=(0−c(c+n)β−12K),then, |
T−1=(−K2c(c+n)β1β−1c(c+n)0). |
Make a change of variables
(u,v)T=T(X,Y)T, |
then, the system (3.8) is transformed into the following form
{X→(1+12K)X−βY+¯F(X,Y)+o(ρ35),Y→βX+(1+12K)Y+¯G(X,Y)+o(ρ35), | (3.9) |
where ρ5=√X2+Y2,
¯F(X,Y)=e20X2+e11XY+e02Y2+e30X3+e21X2Y+e12XY2+e03Y3,¯G(X,Y)=f20X2+f11XY+f02Y2+f30X3+f21X2Y+f12XY2+f03Y3, |
e20=c02βK2c01,e11=c01c11K+2c201d11−c02K22c01,e02=4c201(c20K+2c01d20−d11K)8c01β+K2(c02K−2c01c11)8c01β,e30=c03β2K2c01,e21=(2c01c12−3c03K)βK4c01,e12=c01c21K−c12K2+2d21c2012+3c03K38c01,e03=8c301(c30K+2c01d30−d21K)16c01β−K2(4c201c21+c03K2−2c01c12K)16c01β,f20=c02c01β2,f11=c11β−c02c01βK,f02=c01c20−12c11K+c024c01K2,f30=c03c01β3,f21=c12β2−3c032c01β2K,f12=c01c21β−c12βK+3c034c01βK2,f03=c30c201−12c01c21K+14c12K2−c038c01K3. |
Furthermore,
¯FXX=c02βKc01,¯FXY=c01c11K+2c201d11−c02K22c01,¯FXXX=3c03β2Kc01,¯FYY=4c201(c20K+2c01d20−d11K)4c01β+K2(c02K−2c01c11)4c01β,¯FXXY=c12βK−3c03βK22c01,¯FXYY=c01c21K−c12K2+2d21c201+3c03K34c01,¯FYYY=3c301(c30K+2c01d30)−d21Kc01β−3K2(4c201c21+c03K2−2c01c12K)8c01β,¯GXX=2c02β2c01,¯GXY=c11β−c02βKc01,¯GYY=2c01c20−c11K+c02K22c01,¯GXXX=6c03β3c01,¯GXXY=2c12β2−3c03β2Kc01,¯GXYY=2c01c21β−2c12βK+3c03βK22c01,¯GYYY=6c30c201−3c01c21K+32c12K2−3c03K34c01. |
The {fourth} step. In order to ensure that the system (3.9) has a Neimark-Sacker bifurcation occurring, we need to calculate the discriminating quantity
L=−Re((1−2λ1)λ221−λ1ζ20ζ11)−12|ζ11|2−|ζ02|2+Re(λ2ζ21), | (3.10) |
and L is required not to be zero, where
ζ20=18[¯FXX−¯FYY+2¯GXY+i(¯GXX−¯GYY−2¯FXY)],ζ11=14[¯FXX+¯FYY+i(¯GXX+¯GYY)],ζ02=18[¯FXX−¯FYY−2¯GXY+i(¯GXX−¯GYY+2¯FXY)],ζ21=116[¯FXXX+¯FXYY+¯GXXY+¯GYYY+i(¯GXXX+¯GXYY−¯FXXY−¯FYYY)]. |
By calculation we get
ζ20=18(−4c201(c20K+2c01d20−d11K)4c01β−K2(c02K−2c01c11)4c01β+(2c01c11−c02K)βc01)+18(c02(K2+4β2)2c01−2c01(c20+d11))i,ζ11=14(c02βKc01+4c201(c20K+2c01d20−d11K)4c01β+K2(c02K−2c01c11)4c01β)+14(c02(4β2+K2)2c01+2c01c20−c11K)i,ζ02=18(−4c201(c20K+2c01d20−d11K)4c01β−K2(c02K−2c01c11)4c01β+(3c02K−2c01c11)βc01)+14(c02(4β2−3K2)4c01+c11K+c01(d11−c20))i,ζ21=116(2c01(3c30c01−c21K+c01d21)+c12(12K2+2β2))+116(3c03β(K2+2β2)c01+3K2(4c201c21+c03K2−2c01c12K)8c01β+β(2c01c21−3c12K)−3c301(c30K+2c01d30)−d21Kc01β)i. |
Theorem 3.3. Assume the parameters b, c, m, n in the space SE+={(b,c,m,n)∈R4+|b>0,1>c>m,0<m<c2+2cn−nn+1}. Let b0=c2+2cn−mn−m−n(c+n)2(1−c)(c−m) and L be defined as above (3.10). If L≠0 holds and the parameter a varies in the small neighborhood of b0, then the system (1.6) at the fixed point E3 undergoes a Neimark-Sacker bifurcation. In addition, if L<(or>)0, then an attracting (or repelling) invariant closed curve bifurcates from the fixed point E3 for b<(or>)b0.
In this section, we use the bifurcation diagrams, phase portraits and Lyapunov exponents of the system (1.6) to verify our theoretical results and further reveal some new dynamical behaviors to occur as the parameters vary by Matlab software.
Fix the parameter values c=0.8,m=0.3,n=0.8, let b∈(1.5,3.0) and take the initial values (x0,y0)=(0.55,0.25),(0.80,0.05) in Fig. 2 and Fig. 3 respectively. Figure 1(a) shows the {bifurcation} diagram of (b,x)-plane, from which the fixed point E3 is stable when b<b0=2.266. Moreover, the fixed point E3 is unstable when b>b0. Hence, the Neimark-Sacker bifurcation occurs at the fixed point E3=(0.800,0.625) when b=b0, whose multipliers are λ1,2=0.855±0.519i with |λ1,2|=1.
The corresponding maximum Lyapunov exponent diagram of the system (1.6) is plotted in Figure 1(b). Figures 2(a)-2(f) and Figures 3(a)-3(d) show that the dynamical properties of the fixed point E3 change from stable to unstable as the value of the parameter b decreases and there is an occurrence of invariant closed curve around E3 when b=b0, which agrees to the result of Theorem 3.3.
From the phase portraits in Figs 2 and 3, we infer the stability of E3. Figures 2(d)-2(f) show that the closed curve is stable outside, while Figures 3(a)-3(d) indicate that the closed curve is stable inside for the fixed point E3 as long as the assumptions of Theorem 3.2 hold.
In this paper, we discuss the dynamical behaviors of a predator-prey model (1.6) of Gause-type with double Allee effect affecting the prey population. Under the given parametric conditions, we completely show the existence and stability of four nonnegative equilibria E0=(0,0), E1=(1,0), E2=(m,0) and E3=(c,(1−c)(c−m)c+n). Then we derive the sufficient conditions for its transcritical bifurcation and Neimark-Sacker bifurcation to occur. Meanwhile, it is clear that the positive equilibrium E3=(x0,y0) is asymptotically stable when b<b0=c2+2cn−mn−m−n(c+n)2(1−c)(c−m) and unstable when b>b0 under the condition m<c2+2cn−nn+1. Hence, the system (1.6) undergoes a bifurcation which has been shown to be a Neimark-Sacker bifurcation when the parameter b goes through the critical value b0. Finally, numerical simulations illustrate the theoretical analysis results of the system (1.6).
The perturbations of different parameters in this system may lead to different bifurcations. This demonstrates that this system is sensitive to its parameters. Especially, the occurrence of Neimark-Sacker bifurcation implies that the predator and the prey can coexist under such parametric conditions.
All authors contributed equally and significantly in writing this paper. All authors read and approved the final manuscript.
This work is partly supported by the National Natural Science Foundation of China (61473340), the Distinguished Professor Foundation of Qianjiang Scholar in Zhejiang Province, and the National Natural Science Foundation of Zhejiang University of Science and Technology (F701108G14).
The authors declare that they have no competing interests.
[1] | Goranson U (2007) Damage tolerance facts and fiction. International conference on damage tolerance of aircraft structures, Delft University of Technology, Delft, the Netherlands. |
[2] | Nesterenko B, Nesterenko G (2011) Analysis of requirements on fatigue and damage tolerance for civil transport airplanes. In: Komorowski J editor, Proceedings of the 26th Symposium of the international committee on aeronautical fatigue 39–59. |
[3] | Bellinger NC, Komorowski JP (1997) Corrosion pillowing stresses in fuselage lap joints. AIAA J 35: 317–320. |
[4] | Federal Aviation Administration (2011) Damage Tolerance and Fatigue Evaluation of Structure. Advisory Circular AC 25.571-1D. |
[5] |
Wang W, Rans C, Alderliesten RC, et al. (2015) Predicting the influence of discretely notched layers on fatigue crack growth in fiber metal laminates. Eng Fract Mech 145: 1–14. doi: 10.1016/j.engfracmech.2015.06.062
![]() |
[6] | Wang L, Chow WT, Kawai H, et al. (1996) Prediction of widespread fatigue damage threshold. Proceeding of the FAA/NASA Symposium on Continued Airworthiness of Aircraft Structures, Atlanta (GA), USA. |
[7] | Gruber ML, Wilkins KE, Worden RE (1996) Investigations of fuselage structure subject to widespread fatigue damage. In: Bigelow C editor. Proceedings from the FAA/NASA Symposium on the Continued Airworthiness of Aircraft Structures 439–460. |
[8] |
Galatolo R, Nilsson KF (2001) An experimental and numerical analysis of residual strength of butt-joints panels with multiple site damage. Eng Fract Mech 68: 1437–1461. doi: 10.1016/S0013-7944(01)00027-3
![]() |
[9] | Schijve J (1994) Some calculation on the stress distribution in an infinite sheet with a single crack and with periodic collinear cracks. Report LR-770, TUDelft, the Netherlands. |
[10] | Molinari G, Meneghin I, Melega M, et al. (2012) Parametric damage tolerance design of metallic aeronautical stiffened panels. Aeronaut J 1182: 815–831. |
[11] |
Bloom JM, Sanders JL (1966) The effect of a riveted stringer on the stress of a cracked sheet. ASME J Appl Mech 33: 561–570. doi: 10.1115/1.3625122
![]() |
[12] |
Tada H (1970) Westergaard stress functions for several periodic crack problems. Eng Fract Mech 2: 177–180. doi: 10.1016/0013-7944(70)90023-8
![]() |
[13] | Poe CC (1971) Stress intensity factor for a cracked sheet with riveted and uniformly spaced stringers. NASA TR R-358, Langley Research Center, Hampton (VA), USA. |
[14] | Westergaard HM (1939) Bearing pressures and cracks. J Appl Mech 6: 49–53. |
[15] | Eftis J, Leibowitz H (1972) On the modified Westergaard equation for certain plane crack problems. Int J Fracture 8: 383–392. |
[16] | Carthwrigth DJ, Rooke DP (1976) Compendium of Stress Intensity Factors, HMSO, London, UK. |
1. | Chen Zhang, Xianyi Li, Dynamics of a Discrete Leslie–Gower Model with Harvesting and Holling-II Functional Response, 2023, 11, 2227-7390, 3303, 10.3390/math11153303 | |
2. | Chirodeep Mondal, Ritwika Mondal, Dipak Kesh, Debasis Mukherjee, Dynamics of predator-prey system with the consequences of double Allee effect in prey population, 2025, 51, 0092-0606, 10.1007/s10867-025-09670-0 | |
3. | Debjit Pal, Ritwika Mondal, Dipak Kesh, Debasis Mukherjee, Non-spatial Dynamics and Spatiotemporal Patterns Formation in a Predator–Prey Model with Double Allee and Dome-shaped Response Function, 2025, 87, 0092-8240, 10.1007/s11538-025-01411-7 |