
In this paper, a generalized full orthogonalization method (GFOM) based on weighted inner products is discussed for computing PageRank. In order to improve convergence performance, the GFOM algorithm is accelerated by two cheap methods respectively, one is the power method and the other is the extrapolation method based on Ritz values. Such that two new algorithms called GFOM-Power and GFOM-Extrapolation are proposed for computing PageRank. Their implementations and convergence analyses are studied in detail. Numerical experiments are used to show the efficiency of our proposed algorithms.
Citation: Yu Jin, Chun Wen, Zhao-Li Shen. Acceleration of the generalized FOM algorithm for computing PageRank[J]. Electronic Research Archive, 2022, 30(2): 732-754. doi: 10.3934/era.2022039
[1] | Yichao Shao, Hengguo Yu, Chenglei Jin, Jingzhe Fang, Min Zhao . Dynamics analysis of a predator-prey model with Allee effect and harvesting effort. Electronic Research Archive, 2024, 32(10): 5682-5716. doi: 10.3934/era.2024263 |
[2] | Yujia Xiang, Yuqi Jiao, Xin Wang, Ruizhi Yang . Dynamics of a delayed diffusive predator-prey model with Allee effect and nonlocal competition in prey and hunting cooperation in predator. Electronic Research Archive, 2023, 31(4): 2120-2138. doi: 10.3934/era.2023109 |
[3] | Jianzhi Cao, Xuhong Ji, Pengmiao Hao, Peiguang Wang . Bifurcation analysis of an SIS model with a modified nonlinear incidence rate. Electronic Research Archive, 2025, 33(6): 3901-3930. doi: 10.3934/era.2025173 |
[4] | Zhuo Ba, Xianyi Li . Period-doubling bifurcation and Neimark-Sacker bifurcation of a discrete predator-prey model with Allee effect and cannibalism. Electronic Research Archive, 2023, 31(3): 1405-1438. doi: 10.3934/era.2023072 |
[5] | Xin Du, Quansheng Liu, Yuanhong Bi . Bifurcation analysis of a two–dimensional p53 gene regulatory network without and with time delay. Electronic Research Archive, 2024, 32(1): 293-316. doi: 10.3934/era.2024014 |
[6] | Fengrong Zhang, Ruining Chen . Spatiotemporal patterns of a delayed diffusive prey-predator model with prey-taxis. Electronic Research Archive, 2024, 32(7): 4723-4740. doi: 10.3934/era.2024215 |
[7] | Jiange Dong, Xianyi Li . Bifurcation of a discrete predator-prey model with increasing functional response and constant-yield prey harvesting. Electronic Research Archive, 2022, 30(10): 3930-3948. doi: 10.3934/era.2022200 |
[8] | Mengting Sui, Yanfei Du . Bifurcations, stability switches and chaos in a diffusive predator-prey model with fear response delay. Electronic Research Archive, 2023, 31(9): 5124-5150. doi: 10.3934/era.2023262 |
[9] | Mengxin He, Zhong Li . Dynamic behaviors of a Leslie-Gower predator-prey model with Smith growth and constant-yield harvesting. Electronic Research Archive, 2024, 32(11): 6424-6442. doi: 10.3934/era.2024299 |
[10] | Yanhe Qiao, Hui Cao, Guoming Xu . A double time-delay Holling Ⅱ predation model with weak Allee effect and age-structure. Electronic Research Archive, 2024, 32(3): 1749-1769. doi: 10.3934/era.2024080 |
In this paper, a generalized full orthogonalization method (GFOM) based on weighted inner products is discussed for computing PageRank. In order to improve convergence performance, the GFOM algorithm is accelerated by two cheap methods respectively, one is the power method and the other is the extrapolation method based on Ritz values. Such that two new algorithms called GFOM-Power and GFOM-Extrapolation are proposed for computing PageRank. Their implementations and convergence analyses are studied in detail. Numerical experiments are used to show the efficiency of our proposed algorithms.
Predator-prey interactions have been an interesting and challenging issue that is frequently discussed in marine ecosystems, especially in fish populations. Predator-prey interactions are the most important component of ecology, determining various factors such as community composition, species behavior and dynamics. Mathematical modeling helps to provide insights into the dynamics of the system, which was investigated in the early pioneering work of Lotka[1] and Volterra[2]. In dynamical systems, continuity, equilibrium stability, bifurcation, and control problems are also often studied[3,4,5,6,7]. Traditionally, predators could be distinguished as specialists or generalists based on whether they ate only one or more types of prey. In general, a predator-prey model can be described by the following ordinary differential equation for both specialists and generalist predators[8]:
˙P=F(P)P−G(P,N)N,˙N=γυ(G(P,N))N+H(N)N, |
where P and N denote the prey and predator's population densities at moment t, respectively. F(⋅) and H(⋅) represent the growth of the species in the absence of the other one. G(⋅) is known as the functional response, which characterizes the average individual of prey consumed by a predator; γ represents the conversion rate from prey consumption to predator, υ(⋅) is a monotonically increasing function. For specialist predators, there is H(N)=−d<0; and for generalist predators, it is required that γυ(G(0,N))+H(N)>0. In literature, different type of functional responses were adopted to model different species, which can be prey dependent[9,10] or prey and predator dependent[11,12,13,14,15].
The classical view of population dynamics claims that the higher the population density, the lower the overall growth rate due to the competition for resources. The lower the density, the higher the overall growth rate. However, when the population density is low, Allee[16] introduced the opposite view that the lower the population density, the lower the overall growth rate, that is, the Allee effect. The Allee effect is a common phenomenon in marine populations[17,18,19]. When the population density is low, it may affect population development due to pairing restrictions, dispersal, habitat changes, cooperative foraging, cooperative defense, and predator saturation. Therefore, the study of the Allee effect on ecosystems has attracted the attention of many scholars. In general, the Allee effect can be represented by a multiplier of the form P-L[20,21,22], where L is the threshold for the Allee effect. When L<0, it is a weak Allee effect, and the Allee effect is always positive no matter how much the prey growth rate decreases. When L>0, it is a strong Allee effect. The strong Allee effect indicates that in order for the population to grow, the population size or density must be higher than L; otherwise, the population will die out. Scholars have analyzed the dynamics of the system with the Allee effect of the form P-L, and discussed the existence of the equilibrium of the system and various bifurcation phenomena such as saddle-node bifurcation and B-T bifurcation[23,24,25].
Fish is a kind of important ecological resources. In view of the fish resources development problems, scholars studied population behavior by adding harvest items on the basis of continuous systems. However, fishing activities are not continuous, so continuous dynamic systems cannot accurately describe the actual fishing process. In the process of fish harvesting, the periodic harvesting of fish is a kind of human-controlled behavior that can be described by an impulse differential equation, and it has been found that the impulse differential equation is more accurate in describing and portraying the dynamical behavior of the population[26]. The theory of semi-continuous dynamic systems has been widely used in modeling research on pest management modeling[27,28,29,30,31]. Based on the analysis of the literature, due to the fact that impulsive differential systems (semi-continuous dynamical systems) have the characteristics of both continuous and discrete dynamical systems, there are some studies applying the theory to the development of fish populations in deterministic environments[32,33,34,35,36] and uncertain environment[37,38,39]. In addition, most of the studies considered fish harvesting or fish protection only (unilateral control); in this study, a bilateral state control[40,41,42] is considered, and a predator-prey model for conservation and harvesting of two fish species is developed by constructing a semi-continuous dynamical system. When designing the state feedback control strategy, the number of objective fish was used as the state variable for feedback control. On the one hand, when the number of prey fish is small, the Allee effect will lead to their extinction, which will lead to the lack of enough food for prey fish and destroy the balance of the ecosystem. Therefore, it is necessary to release a certain amount of prey fish when the number of prey fish decreases to a certain threshold. When the number of prey fish is high, it is necessary to catch prey fish from the economic point of view. Since the fishing behavior will also lead to the harvesting of some predator fish, in order to maintain the ecological balance and avoid the extinction of predator fish caused by the fishing behavior, it is necessary to release a certain amount of predator fish larvae at the same time. Based on the above two aspects, we propose a bilateral control strategy to maintain the population size of the two species in a suitable range.
This paper considers a predator-prey model in which the predator is a generalist and the growth of the prey is affected by the Allee effect. The structure is as follows: In Section Ⅱ, we describe the fishery model, non-negativity, persistent survivability and discuss the existence and stability of equilibria of the model. In Section Ⅲ, the bifurcation dynamics of the model are discussed using bifurcation theory. In Section Ⅳ, based on this model, we analyze the model of marine fish harvesting and conservation with bilateral controls. The existence and stability of the order-1 and order-2 periodic solutions of the system are analyzed by using the geometry theory of differential equations. In the fifth section, we performed numerical simulations using MATLAB to verify the correctness of the results.
In this paper, we present a predator-prey model in which predators are generalists, prey growth rates are logical and subject to strong Allee effects, and the functional response is a Holling-Ⅰ type, while the conversion from prey consumption to predator species is saturated,
{dPdT=rP(1−PK)(P−L)−APN,dNdT=e(AP1+BP)N+(s1+fN−d)N, | (2.1) |
where p(T), N(T) denote the prey and predator's densities at the moment of T; K denotes the prey's environmental holding capacity; L is the threshold of the prey's Allee effect; r and s represent the intrinsic growth rates of prey and predator, respectively; A is the capture rate of prey by the predator; e is the efficiency with which prey is converted to predator; B is the half-saturation constant; f is the intensity of predator density dependence, and d is the predator mortality rate. Since the predators are generalist, it requires that eA≤Bd, s>d, and all parameters of model (2.1) are positive.
To facilitate the analysis, let x=PK, y=N, t=rKT, α=ArK, β=LK, γ=eAr, δ=BK, s1=srK, d1=drK. Then system (2.1) is simplified to system (2.2):
{dxdt=x(1−x)(x−β)−αxy,dydt=γxy1+δx+(s11+fy−d1)y, | (2.2) |
and from the biological point of consideration, the model (2.2) is limited in the region
Ω={(x,y)∈R2+|0≤x≤1,y≥0}. |
On the other hand, as a renewable resource, fish species are closely related to human life. In order to maintain the balance of the ecosystem during the fishing process, we consider a fish stock control strategy with a combination of fishing and investment. First, to avoid the distinction of prey fish caused by the Allee effect, a quantity (η) of juvenile prey fish is released when the prey density declines to the level x=h1. But at a higher level x=h2, a proportion a of prey fish together with a proportion b of predator fish will also be caught for economic purposes, and simultaneously, a quantity τ of juvenile predator fish is released into the system to maintain the level of predator fish. Based on the control measures, the model can be described as follows:
{dxdt=x(1−x)(x−β)−αxydydt=γxy1+δx+(s11+fy−d1)y}h1<x<h2,Δx=ηΔy=0}x=h1,Δx=−axΔy=−by+τ}x=h2, | (2.3) |
where η, a, b, τ are all positive, and a,b∈(0,1).
For a given planar model
{dxdt=f1(x,y),dydt=f2(x,y)ω(x,y)≠0,Δx=I1(x,y),Δy=I2(x,y)ω(x,y)=0, | (2.4) |
Definition 1 (Order-k periodic solution[32,33,36]) The solution ˜z(t)=(˜x(t),˜y(t)) is called periodic if there exists m(⩾1) satisfying ˜zm=˜z0. Furthermore, ˜z is an order-k T-periodic solution with k≜min{j|1≤j≤m,˜zj=˜z0}.
Lemma 1 (Analogue of Poincaré Criterion[32,33,36]). The order-k T-periodic solution z(t)=(ξ(t),η(t))T is orbitally asymptotically stable if |μq|<1, where
μk=k∏j=1Δjexp(∫T0[∂f1∂x+∂f2∂y](ξ(t),η(t))dt), |
with
Δj=f+1(∂I2∂y∂ω∂x−∂I2∂x∂ω∂y+∂ω∂x)+f+2(∂I1∂x∂ω∂y−∂I1∂y∂ω∂x+∂ω∂y)f1∂ω∂x+f2∂ω∂y, |
f+1=f1(ξ(θ+j),η(θ+j)), f+2=f2(ξ(θ+j),η(θ+j)) and f1, f2, ∂I1∂x, ∂I1∂y, ∂I2∂x, ∂I2∂y, ∂ω∂x, ∂ω∂y are calculated at (ξ(θj),η(θj)).
In this section, the bounded-ness of the solution for Model (2.2) is discussed. Moreover, the existence, type, and local stability of the equilibrium as well as the bifurcation properties are verified.
Define
g1(x,y)=(1−x)(x−β)−αy,f1(x,y)=xg1(x,y); |
g2(x,y)=γx1+δx+(s11+fy−d1),f2(x,y)=yg2(x,y). |
Theorem 1. The solution of Model (2.2) with non-negative initial values will remain non-negative for all time and is bounded on R2+.
Proof. By Eq (2.2), it can be obtained that
x(t)=x(0)exp[∫t0g1(x(s),y(s))ds],y(t)=y(0)exp[∫t0g2(x(s),y(s))ds], |
for all t≥0, as long as x(0) and y(0) are non-negative, then x(t) and y(t) are also non-negative.
Next, we define a function u(t)=γαx(t)+y(t). Then
˙u=γα˙x+˙y=γα[x(1−x)(x−β)−αxy]+γxy1+δx+(s11+fy−d1)y≤γα[x(1−x)(x−β)−αxy]+γxy+(s11+fy−d1)y≤γ(1−x)(x−β)αx+s1f−d1y≤[γα((1−β)24+d1)+s1f]−d1u, |
which implies that
u(t)≤1d1[γα((1−β)24+d1)+s1f]+(u0−1d1[γα((1−β)24+d1)+s1f])e−d1t, |
so long as u0=γαx0+y0 is bounded, u(t) is bounded in Ω. To sum up, any solution of Model (2.2) starting with a non-negative bounded initial condition is non-negative and bounded in Ω.
Obviously, four equilibria always exist
O(0,0),E1(0,1f(s1−d1d1)),E2(β,0),E3(1,0). |
Define
y1(x)=1α(1−x)(x−β), | (3.1) |
y2(x)=1f[s1δ(1δ+x)d1−(γ−d1δ)x−1] | (3.2) |
and denote ¯γ1≜d1δ. Due to the assumptions eA≤Bd and s>d, then γ<¯γ1, i.e., y2(x)<s1δ¯γ1−γ. The positional relationship between y1(x) and y2(x) for different cases is shown in Figure 1.
Let y1(x)=y2(x). Then it has
a1x3+a2x2+a3x+a4=:g(x)=0, | (3.3) |
where
a1=f(δd1−γ)>0,a2=f[(γ−δd1)(β+1)+d1],a3=α[δ(s1−d1)+γ]+f[β(δd1−γ)−d1(β+1)],a4=α(s1−d1)+βfd1>0. |
Define
¯α≜fβ,¯δ1≜d1(1+β)s1β,¯δ2≜fd1(1+β)fβd1+α(s1−d1),¯δ3≜(1+β)fd1αs1, |
¯γ2≜αδ(s1−d1)+fd1(βδ−(1+β))fβ−α,xd≜√a22−3a1a3−a23a1,ρ≜g(xd). |
Theorem 2. For any of the following cases: C1) α<¯α, δ<¯δ2, γ<¯γ1; C2) α<¯α, ¯δ2<δ<¯δ3, ¯γ2<γ<¯γ1; C3) α>¯α, δ<¯δ2, 0<γ<−¯γ2, if ρ<0 holds, there exists two interior equilibria in Model (2.2); if ρ=0, there exists a unique interior equilibrium in Model (2.2); if ρ>0, there doesn't exists interior equilibrium in Model (2.2).
Proof. Clearly, the existence of an interior equilibrium is equivalent to that of a positive root of Eq (3.3) in the interval (0, 1). Since
g(0)=a4>0,g(1)=a1+a2+a3+a4=α[(δ+1)(s1−d1)+γ]>0, |
and
g′(x)=3a1x2+2a2x+a3, | (3.4) |
then for any of the cases C1)–C3), there is a3<0. Moreover, g′(0)<0, g′(xd)=0, g″(xd)>0 and
g′(1)=(1−β)f(δd1−γ)+(1−β)d1f+α[δ(s1−d1)+γ]>0, |
so that for x∈(0,xd), g′(x)<0, for x∈(xd,0), g′(x)>0. If ρ<0, Eq (3.3) has two distinct positive roots x∗i∈(0,1), i=1,2. Denote y∗i=y1(x∗i), i=1,2. Then two interior equilibria exist in Model (2.2), denoted by E∗1(x∗1,y∗1) and E∗2(x∗2,y∗2). If ρ=0, Eq (3.3) has a unique positive root x∗=xd∈(0,1), then a unique positive equilibrium exists in Model (2.2), denoted by E∗(x∗,y1(x∗); when ρ>0, Eq (3.3) does not have positive root, and thus there does not exist interior equilibrium in Model (2.2).
For any equilibrium ¯E(¯x,¯y), there is
J(¯E)=[−3¯x2+2(β+1)¯x−β−α¯y−α¯xγ¯y(1+δ¯x)2γ¯x1+δ¯x+s1(1+f¯y)2−d1], |
its characteristic equation is
λ2−TR(J(¯E))λ+DET(J(¯E))=0, |
where
TR(J(¯E))=−3¯x2+2(β+1)¯x−β−α¯y+γ¯x1+δ¯x+s1(1+f¯y)2−d1, |
DET(J(¯E))=[−3¯x+2(β+1)¯x−β−α¯y][γ¯x1+δ¯x+s1(1+f¯y)2−d1]+α¯x(γ¯y(1+δ¯x)2). |
1) Boundary equilibria
At O(0,0), there is
J(O)=[−β00s1−d1], |
since λ1=−β<0, λ2=s1−d1>0, then O is an unstable higher-order singularity.
At E1(0,1f(s1−d1d1)), there is
J(E1)=[−β−αf(s1−d1d1)0γf(s1−d1d1)0], |
since λ1=0, λ2=−αf(s1−d1d1)<0, so that E1 is locally stable.
At E2(β,0), there is
J(E2)=[β(1−β)−αβ0γβ1+δβ+s1−d1], |
since λ1=β(1−β)>0, λ2=γβ1+δβ+s1−d1>0, then E2 is unstable.
At E3(1,0), there is
J(E3)=[β−1−α0γ1+δ+s1−d1], |
since λ1=β−1<0, λ2=γ1+δ+s1−d1>0, then E3 is unstable.
2) Interior equilibrium
At E∗(x∗,y∗), there are g1(x∗,y∗)=0 and g2(x∗,y∗)=0, then
J(E∗)=[(1+β)x∗−2(x∗)2−αx∗γy∗(1+δx∗)2−s1fy∗(1+fy∗)2], |
thus,
DET(J(E))=[x∗y∗∂g1∂y∂g2∂y(dy1dx−dy2dx)](x∗,y∗), |
where x∗y∗∂g1∂y∂g2∂y|x∗,y∗>0, then the sign of DET(J(E∗)) is identical to that of dy1dx|x=x∗−dy2dx|x=x∗. Next, it will discuss the sign of dy1dx|x=x∗−dy2dx|x=x∗ for different cases in Theorem 2.
ⅰ) When ρ<0 holds, two interior equilibria E∗1(x∗1,y∗1) and E∗2(x∗2,y∗2) with 0<x∗1<x∗2<1 exists in Model (2.2), as illustrated in Figure 2(a). It can be easily checked that
SIGN(J(E∗1))=[+−+−],SIGN(J(E∗2))=[−−+−]. |
Besides, at E∗1, there is dy1dx|x=x∗1>dy2dx|x=x∗1, thus
DET(J(E∗1))=[xy∂f1∂y∂f2∂y(dy1dx−dy2dx)](x∗1,y∗1)<0, |
i.e., E∗1 is unstable. Similarly, at E∗2, there is dy1dx|x=x∗2<dy2dx|x=x∗2, thus
(λ1+λ2)|(x∗2,y∗2)=TR(J(E∗2))<0, |
(λ1λ2)|(x∗2,y∗2)=DET(J(E∗2))=[xy∂f1∂y∂f2∂y(dy1dx−dy2dx)](x∗2,y∗2)>0, |
i.e., E∗2 is a locally asymptotically stable node.
ⅱ) When ρ=0 holds, system (2.2) has a unique interior equilibrium E∗(x∗,y∗), as shown in Figure 2(b). Let X=x−x∗,Y=y−y∗, then E∗ is converted to the origin O(0,0), and the model is written as
{dXdt=a11X+a12Y+A1X2+A2XY,dYdt=a21X+a22Y+B1X2+B2XY+B3Y2+P3(X,Y), | (3.5) |
where,
a11=(1+β)x∗−2(x∗)2,a12=−αx∗,a21=γy∗(1+δx∗)2,a22=−s1fy∗(1+fy∗)2,A1=β+1−3x∗,A2=−α2,B1=−γδy∗(1+δx∗)3,B2=γ(1+δx∗)2,B3=−s1f(1+fy∗)3 |
and P3(X,Y) is a function of (X,Y) with degree of three or higher. The Jacobian matrix at E∗ is
J(E∗)=[a11a12a21a22], |
thus
DET(J(E∗))=a11a22−a12a21=0,TR(J(E∗))=a11+a22. |
a) If TR(J(E∗))=a11+a22=0, then λ1=λ2=0. Make the transformation
(XY)=(a110a211)(x1y1), |
then Model (3.5) is converted into the following standard form:
{dx1dt=ˉa12y1+ˉA1x21+ˉA2x1y1,dy1dt=ˉB1x21+ˉB2x1y1+ˉB3y21+ˉP3(x1,y1), |
where
ˉa12=a12a11,ˉA1=a11A1+a21A2,ˉA2=A2, |
ˉB1=a211B1+a11a21(B2−A1)+a221(B3−A2),ˉB2=a11B2+a21(2B3−A2),ˉB3=B3, |
and ˉP3(X,Y) is a function of three or more degrees about (X,Y).
Let τ=ˉa12t. For convenience, t is still used to represent τ, then it has
{dx1dt=y1+˜A1x21+˜A2x1y1,dy1dt=˜B1x21+˜B2x1y1+˜B3y21+˜P3(x1,y1), | (3.6) |
where
˜Ai=ˉAiˉa12,i=1,2,˜Bj=ˉBjˉa12,j=1,2,3. |
and ˜P3(X,Y) is a function of (X,Y) with three or higher degree. Model (3.6) can be transformed into the following form[43]:
{dx1dt=y1,dy1dt=˜B1x21+(˜B2+2˜A1)x1y1+˜P′3(x1,y1), |
and if ˜B1≠0, then E∗ is a cusp point.
Meanwhile, if ˜B2+2˜A1≠0, E∗ is a cusp of codimension two. If ˜B2+2˜A1=0, E∗ is a cusp with at least codimension three.
b) If TR(J(E∗))=a11+a22≠0, then λ1=0, λ2≠0. Make the transformation
(XY)=(a22a11−a21a21)(x1y1), |
system (3.5) is converted into the following standard form
{dx1dt=A′1x21+A′2x1y1+A′y21,dy1dt=a′22y1+B′1x21+B′2x1y1+B′3y21+P′3(x1,y1), |
where
A′1=a21a222A1−a221a22A2−a11a222B1+a11a21a22B2−a11a221B3a21(a11+a22),A′2=2a211a22+a221a22A2−a11a221A2−2a211a22B1−a11a21a22B2+a11a221B2+2a11a221B3a21(a11+a22),A′3=a211a21A1+a11a221A2−a311B1−a211a21B2−a11a221B3a21(a11+a22),B′1=a21a222A1−a221a22A2+a322B1−a21a222B2+a221a22B3a21(a11+a22),B′2=2a11a21a22A2+a221a22A2−a11a221A2+2a11a222B1+a21a222B2−a11a21a22B2−2a221a22B3a21(a11+a22),B′3=a211a21A1+a11a221A2+a11a222B1+a11a21a22B2+a221a22B3a21(a11+a22). |
Next, introduce a new variable τ=a′22t (for convenience, is represented by t), then
{dx1dt=ˆA1x21+ˆA2x1y1+ˆA3y21,dy1dt=y1+ˆB1x21+ˆB2x1y1+ˆB3y21+ˆP3(x1,y1), |
where, ˆAi=A′ia′22,ˆBi=B′ia′22,i=1,2,3.
If ˆA1≠0 then E∗ is unstable. Meanwhile, E∗ is a saddle node of attraction.
To sum up, the following result hold.
Theorem 3. For Model (2.2), 1) O(0,0) is unstable, E1(0,1f(s1−d1d1)) is locally stable, E2(β,0), E3(1,0) is unstable; 2) For the interior equilibrium, when ρ>0, E∗1(x∗1,y∗1) is a saddle (unstable), and E∗2(x∗2,y∗2) is a stable node; when ρ=0, if TR(J(E∗))=0 and ˜B1≠0, E∗(x∗,y∗) is a cusp of codimension two in case of ˜B2+2˜A1≠0, and a cusp with at least codimension three in case of ˜B2+2˜A1=0; If TR(J(E∗))≠0 and ˆA1≠0, E∗(x∗,y∗) is an attractive saddle node.
Let α=α0 satisfy
DET(J(E∗))=(λ1λ2)|(x∗,y∗)=0,TR(J(E∗))=−(λ1+λ2)|(x∗,y∗)≠0. |
Denote
ξ1=2(x∗)2−(1+β)x∗,ξ2=−γy∗(1+δx∗)2,ω1=αx∗,ω2=s1fy∗(1+fy∗)2 |
and define
Φ2Δ=2ω21ξ2x∗−2γδy∗(1+δx∗)3ω21ξ1−2s1f2(1+fy∗)3ξ31. |
Theorem 4. Let the parameters of Model (2.2) satisfy TR(J(E∗))≠0 and ˆA1≠0. If Φ2≠0, system (2.2) undergoes a saddle node bifurcation near E∗(x∗,y∗) when α=α0.
Proof. At E∗(x∗,y∗), there is
J(E∗)=(−ξ1−ω1−ξ2−ω2). |
Let V=(V1,V2)T (W=(W1,W2)T) be the eigenvector corresponding to the zero eigenvalue of J(E∗) ((J(E∗))T). Then V=(ω1,−ξ1)T, W=(−ξ2,ξ1)T. Let g=(g1,g2)T. Then
gα(E∗,α0)=∂g∂α(E∗,α0)=(−y∗0). |
and
D2g(E∗,α0)(V,V)=(∂2g1∂x2V21+2∂2g1∂x∂yV1V2+∂2g1∂y2V22∂2g2∂x2V21+2∂2g2∂x∂yV1V2+∂2g2∂y2V22)=(−2ω21−2γδ(1+δx∗)3ω21−2s1f2(1+fy∗)3ξ21). |
Since
Φ1=WTgα(E∗,α0)=−x∗ξ1≠0 |
and
Φ2=WT[D2g(E∗,α0)(V,V)]=2ω21ξ2−2γδ(1+δx∗)3ω21ξ1−2s1f2(1+fy∗)3ξ31≠0, |
then according to Sotomayor's theorem[44], Model (2.2) undergoes a saddle-node bifurcation near E∗(x∗,y∗) when α=α0.
Remark 1. For Model (2.2), it can be concluded that for a lager α, the interior equilibrium does not exist. When α decreases to α=α0, a unique equilibrium exists in the system. As α decreases below α0, the system undergoes a saddle-node bifurcation at the interior equilibrium E∗, giving rise to two interior equilibria E∗1 and E∗2.
According to Theorem 3, when ρ=0, TR(J(E∗))=0, ˜B1≠0 and ˜B2+2˜A1≠0, Model (2.2) to undergo a Bogdanov-Takens bifurcation of codimension two near E∗ when (α,β)=(α0,β0). Next, it will show the universal unfolding of the Bodmanov-Takens bifurcation of codimension two under parameter perturbation when α and β are taken as bifurcation parameters.
Define
b01=x∗(1−x∗)(x∗−ϵ2)−ϵ1x∗y∗,b11=[1+(β+ϵ2)]x∗−2(x∗)2,b12=−(α+ϵ1)x∗,b10=0,b21=γy∗(1+δx∗)2,b21=−s1fy∗(1+fy∗)2,E1=(β+ϵ2)+1−3x∗,E2=−(α+ϵ1)2,F1=−γδy∗(1+δx∗)3,F2=γ(1+δx∗)2,F3=−s1f(1+fy∗)3,H1=−b212b10+b01b12b22−b201F3b12,H2=b201E22+b212b10E22+b312b21−b11b212b22+b01b212F2b212,H3=b11b12−b01E2+b12b22−2b10F3b12,H4=−b01b12E1E2−b01b11E22−b212b22E1−b212F1+b11b212F2−b211b12F3b212,H5=2b212E21−b11b12E2−2b01E22+b212F2−2b11b12F3b212,H6=E2+F3b12,J1=H1,J2=H2−2H1H6,J3=H3J4=H1H26−2H2H6+H4,J5=H5−H3H6, |
K1=−J1J4,K2=−J2J4,K3=−J3√−J4,K4=−J5√−J4,L1=K1+K224,L2=K3+K2K42,L3=K4,O1=−L1L43,O2=−L2L3, |
Let (ϵ1,ϵ2) be a parameter vector in a small neighborhood of (0,0). Then
Theorem 5. Let the parameters of Model (2.2) satisfy ρ=0, TR(J(E∗))=0, ˜B1≠0, ˜B2+2˜A1≠0 and |∂(O1,O2)∂(ϵ1,ϵ2)|≠0. When (α,β) varies in the neighborhood of (α0,β0), Model (2.2) changes in the small neighborhood of E∗(x∗,y∗), and a codimensional 2 Bogdanov-Takens branching occurs.
Proof. Consider the perturbation system
{dxdt=x(1−x)[x−(γ+ϵ2)]−(α+ϵ1)xy≜F(x,y),dydt=γxy1+δx+(s11+fy−d1)y≜G(x,y), |
For (α,β)=(α0,β0), there is DET(J(E∗))=0,TR(J(E∗))=0. With the transformation x1=x−x∗, y1=y−y∗, we can obtain
{dx1dt=b01+b11x1+b12y1+E1x21+E2x1y1,dy1dt=b10+b21x1+b22y1+F1x21+F2x1y1+F3y21+N3(x1,y1), | (4.1) |
where N3(x1,y1,ϵ1,ϵ2)∈C∞.
Make the following transformation:
{x2=x1,y2=b01+b11x1+b12y1+E1x21+E2x1y1, |
then Model (4.1) is converted to
{dx2dt=y2,dy2dt=H1+H2x2+H3y2+H5x22+H5x2y2+H6y22+N′3(x2,y2,ϵ1,ϵ2), |
where N′3(x2,y2,ϵ1,ϵ2)∈C∞ with coefficients depending smoothly on ϵ1,ϵ2.
Next, introduce the variable τ, denoted as dt=(1−H6x2)dτ (still denote τ as t), then
{dx2dt=(1−H6x2)y2,dy2dt=(1−H6x2)(H1+H2x2+H3y2+H4x22+H5x2y2+H6y22+N′3(x2,y2,λ1,λ2)), |
Let x3=x2,y3=(1−H6x2)y2, then the above system of equations is transformed into
{dx3dt=y3,dy3dt=J1+J2x3+J3y3+J4x23+J5x3y3+˜N3(x3,y3,ϵ1,ϵ2), |
where ˜N3(x3,y3,ϵ1,ϵ2)∈C∞ with coefficients depending smoothly on ϵ1,ϵ2.
(ⅰ) When J4<0, the following transformations are applied to the variables:
x4=x3,y4=y3√−J4,τ=√−J4t, |
still denote τ as t, there is
{dx4dt=y4,dy4dt=K1+K2x4+K3y4−x24+K4x4y4+M3(x4,y4,ϵ1,ϵ2), |
where M3(x4,y4,ϵ1,ϵ2)∈C∞ with at least third order.
Let x5=x4−K22, y5=y4, and obtain
{dx5dt=y5,dy5dt=L1+L2y5−x25+L3x5y5+M′3(x5,y5,ϵ1,ϵ2), |
where M′3(x5,y5,ϵ1,ϵ2)∈C∞ with at least third order.
If set J5≠0, then L3≠0. Define new variables: x6=−L23x5, y6=L33y5, τ=−tL3, and denote x6 by x, y6 by y, and τ by t, which yields that
{dxdt=y,dydt=O1+O2y+x2+xy+˜M3(x,y,λ1,λ2), | (4.2) |
where ˜N3(x2,y2,ϵ1,ϵ2)∈C∞ with at least third order, and O1, O2 can be represented by ϵ1 and ϵ2.
(ⅱ) When J4>0, the following transformations are applied
x4=x3,y4=y3√J4,τ=√J4t, |
still denote τ by t, it has
{dx4dt=y4,dy4dt=K′1+K′2x4+K′3y4+x24+K′4x4y4+N3(x4,y4,ϵ1,ϵ2), |
where N3(x4,y4,ϵ1,ϵ2)∈C∞ with at least third order, and
K′1=J1J4,K′2=J2J4,K′3=J3√J4,K′4=J5√J4. |
Let x5=x4+K′2/2,y5=y4. Then it has
{dx5dt=y5,dy5dt=L′1+L′2y5+x25+L′3x5y5+N′3(x5,y5,ϵ1,ϵ2), |
where N′3(x5,y5,ϵ1,ϵ2)∈C∞ with at least third order, and
L′1=K′1−K′224,L′2=K′3−K′2K′42,L′3=K′4. |
If set J5≠0, then L′3≠0. Define new variables: x6=L′32x5, y6=L′33y5, τ=t/L′3, and still denote x6 by x, y6 by y, τ by t, which yields that
{dxdt=y,dydt=O′1+O′2y+x2+xy+˜N3(x2,y2,ϵ1,ϵ2), | (4.3) |
where O′1=L′1L′34,O′2=L′2L′3, and ϵ1, ϵ2 can be represented by O′1, O′2.
For convenience of discussion, O′1, O′2 is still denoted by O1,O2. When |∂(O1,O2)∂(ϵ1,ϵ2)|≠0, Models (4.2) and (4.3) are the cardinal folds of the Bogdanov-Takens bifurcation[44], when (α,β) varies in the vicinity of (α0,β0), Model (2.2) undergoes a codimension 2 Bogdanov-Takens bifurcation in a small neighborhood of E∗(x∗,y∗).
It only focuses on the case of ρ<0 in Theorem 3. For Model (2.3), there are
M1={(x,y)∈R2+:x=h1,y≥0},M2={(x,y)∈R2+:x=h2,y≥0}, |
N1={(x,y)∈R2+:x=h1+η,y≥0},N2={(x,y)∈R2+:x=(1−a)h2,y≥τ}. |
Denote l1=1α(1−x)(x−β) as the prey isocline, l2=(s1−d1)(1+δx)−γxf[d1+(d1δ−γ)x] as the predator's isocline, l3, l4, l5, l6 as the saddle point separatrix of E∗1 in different directions. The intersection point between l1 and N2 is denoted by A3. The intersection point between l4 and N2 is denoted by M. The trajectory from A3 intersects M2 at the point B3. Define ¯τ2≜yM−(1−b)yB3. Denote
Ω1={(x,y)|0≤x≤x∗1,y≥0},Ω2={(x,y)|x∗1≤x≤1,y≥0}. |
Definition 2 (Successor function). For a point S∈N1, if the trajectory from S intersects M1, then denote the intersection point by S−∈M1. Under the impulse effect, the point S− is mapped to S+∈N1. In such a case, we can define fISOR1: N1→R, S→fISOR1(S)≜yS+−yS. If the trajectory from S intersects M2, then denote the intersection point by S−. Under the impulse effect, the point S− is mapped to S+∈N2. If the trajectory from S+ intersects M1, denote the intersection point by S+−∈M1. Under the impulse effect, the point S+− is mapped to S++∈N1. In such cases, we can define fIISOR1: N1→R, S→fIISOR1(S)≜yS++−yS. Similarly, For a point S′∈N2, define fISOR2: N2→R, S′→fISOR2(S′)≜yS′+−yS′.
Theorem 6. For system (2.3) with model parameters satisfying any one of C1)–C3) and ρ>0, if D-1) h1+η<x∗1<(1−a)h2<h2<x∗2, y2(h1)≥y1(h1+η) and τ<¯τ2, an order-1 periodic solution exists in each Ωi, i=1,2; if D-2) h1<(1−a)h2<h1+η<x∗1<h2<x∗2 and y2(h1)≥y1(h1+η), an order-1 periodic solution exists in Ω1; if D-3) h1<x∗1<h1+η<(1−a)h2<h2<x∗2 and τ<¯τ2, an order-1 periodic solution exists in Ω2.
Proof. For case D-1) h1+η<x∗1<(1−a)h2<h2<x∗2 and y2(h1)≥y1(h1+η), denote the intersection point between l2 and N1 by A0. Select a point A1∈N1 above A0; the trajectory from A1 intersects M1 at B1. Under the impulse effect, the point B1 is mapped to A+1∈N1. Since g1(xA1,yA1)<0, g2(xA1,yA1)<0, then we have fISOR1(A1)=yA+1−yA1=yB1−yA1<0. Besides, denote the intersection between l1 and N1 by A2. Since g1(xA2,yA2)=0, g2(xA2,yA2)>0 and g1(xS,yS)<0 for S∈U(A1,ϵ) with xS<h1+η and yS<yA1, then we have fISOR1(A2)=yA+2−yA2=yB2−yA2>0 due to y2(h1)≥y1(h1+η). The continuity of fISOR1 implies that a point S∈¯A1A2 exists so that fISOR1(S)=0, i.e., an order-1 periodic solution exists in Ω1 (Figure 3(a)). Similarly, it can be proved that an order-1 periodic solution exists in Ω1 for case D-2).
Under the impulse effect, B3 is mapped to A3+∈N2, where yA3+=(1−b)yB3+τ. Define ¯τ1≜yA3−(1−b)yB3.
ⅰ) If τ=¯τ1, then fISOR2(A3)=yA3+−yA3=0, i.e., an order-1 periodic solution exists in Ω1 (Figure 3(a)).
ⅱ) If τ<¯τ1, then fISOR2(A3)=yA3+−yA3<0. On the other hand, for A4((1−a)h2,τ), there is fISOR2(A4)=yA4+−yA4>0. The continuity of fISOR2 implies that a point S′∈¯A3A4 exists so that fISOR2(S′)=0, i.e., an order-1 periodic solution exists in Ω2 (Figure 3(b)).
ⅲ) If ¯τ1<τ<¯τ2, then fISOR2(A3)=yA3+−yA3>0. According to the trend of the trajectory and the fact that any two trajectories cannot be intersected, select D∈N2 above and sufficiently close to the A3, then fISOR2(D)=yD+−yD>0. On the other hand, for C=A+3, there is fISOR2(C)=yC+−yC<yA+3−yC=0. The continuity of fISOR2 implies that a point S′∈¯CD exists so that fISOR2(S′)=0, i.e., an order-1 periodic solution exists in Ω2 (Figure 3(c)).
Similarly, system (2.3) has an order-1 periodic solution in Ω2 for case D-3).
To sum up, an order-1 periodic solution exists in Model (2.3) for case D-1).
Let zi(t)=(ϕi(t),φi(t)) (k−1)Ti≤t≤kTi, k∈N be the order-1 periodic solution in Ωi, i=1,2. For z1(t)=(ϕ1(t),φ1(t)) (k−1)T1≤t≤kT1, denote
ϕ1(T+1)=ϕ1(T1)+η=h1+η,φ1(T+1)=φ1(T1)=δ1. |
For z2(t)=(ϕ2(t),φ2(t)) (k−1)T2≤t≤kT2, denote
ϕ2(T+2)=(1−a)ϕ2(T2)=(1−a)h2,φ2(T+2)=(1−b)φ2(T2)+τ=δ2. |
Theorem 7. For the model parameters with any one of C1)–C3), ρ>0 and D-1), zi(t)=(ϕi(t),φi(t)) (k−1)Ti≤t≤kTi is orbitally asymptotically stable if μi<1, where
μ1=|[1−(h1+η)][(h1+η)−β]−αδ1(1−h1)(h1−β)−αδ1|exp{∫T10+[ϕ1(t)(β+1−2ϕ1(t))−s1fφ1(t)(1+fφ1(t))2]dt},μ2=|{[1−(1−a)h2][(1−a)h2−β]−αδ2}(δ2−τ)δ2[(1−b)(1−h2)(h2−β)−α(δ2−τ)]|Λ(t) |
with Λ(t)=exp{∫T20+[ϕ2(t)(β+1−2ϕ2(t))−s1fφ2(t)(1+fφ2(t))2]dt}.
Proof. For system (2.3), there are
f1(x,y)=x(1−x)(x−β)−αxy,f2(x,y)=γxy1+δx+(s11+fy−d1)y, |
and ω1(x,y)=x−h1,I11(x,y)=η,I21(x,y)=0, So it can be concluded that
∂f1(x,y)∂x=˙xx+x(β+1−2x),∂f2(x,y)∂x=˙yy−s1fy(1+fy)2, |
∂ω1(x,y)∂x=1,∂ω1(x,y)∂y=∂I11(x,y)∂x=∂I11(x,y)∂y=∂I21(x,y)∂x=∂I21(x,y)∂y=0. |
Denote f1+≜f1(ϕ1(T+1),φ1(T+1)),f2+≜f2(ϕ1(T+1),φ1(T+1)). Then by Lemma 1, there are
κ1=(∂I21∂x⋅∂ω1∂x−∂β1∂x⋅∂ω1∂y+∂ω1∂x)f1++(∂I11∂x⋅∂ω1∂y−∂I11∂y⋅∂ω1∂x+∂ω1∂y)f2+∂ω1∂xf1+∂ω1∂yf2=(h1+η)[1−(h1+η)][(h1+η)−β]−α(h1+η)δ1h1(1−h1)(h1−β)−αh1δ1 |
and
μ1=|κ1|exp[∫T10+(∂f1∂x+∂f2∂y)dt]=|[1−(h1+η)][(h1+η)−β]−αδ1(1−h1)(h1−β)−αδ1|exp{∫T10+[ϕ1(t)(β+1−2ϕ1(t))−s1fφ1(t)(1+fφ1(t))2]dt}. |
If μ1<1, the order-1 periodic solution z1(t)=(ϕ1(t),φ1(t)) (k−1)T1≤t≤kT1 is orbitally asymptotically stable.
Similarly, for the order-1 periodic solution z2(t)=(ϕ2(t),φ2(t)) (k−1)T2≤t≤kT2, there is
μ2=|{[1−(1−a)h2][(1−a)h2−β]−αδ2}(δ2−τ)δ2[(1−b)(1−h2)(h2−β)−α(δ2−τ)]|Λ(t), |
where Λ(t)≜exp{∫T20+[ϕ2(t)(β+1−2ϕ2(t))−s1fφ2(t)(1+fφ2(t))2]dt}. By Lemma 1, if μ2<1, the order-1 periodic solution z2(t)=(ϕ2(t),φ2(t)) (k−1)T2≤t≤kT2 is orbitally asymptotically stable. The proof is completed.
Here only the case for h1<(1−a)h2<x∗1<h1+η<h2<x∗2 is considered. Let G1(h1+η,yG1) be the intersection point between l1 and N1, G2(h1,yG2) be the intersection point between l2 and M1. The trajectory that starts with G1 intersects M2 at B1(h2,yB1). Define ¯τ3≜yG1−(1−b)yB1.
Theorem 8. For the model parameters with any one of C1)-C3), ρ>0 and D-4) h1<(1−a)h2<x∗1<h1+η<h2<x∗2, if byG2≤τ≤min{yG2,¯τ3} holds, an order-2 periodic solution exists in system (2.3).
Proof. For G1∈N1, the trajectory starting from G1 intersects with M2 at B1. Then B1 is mapped to B″1∈N2, and next intersects with M1 at ˆB1, and then ˆB1 is mapped to G+1 under impulse effect. Since τ<¯τ3, then yB″1=(1−b)yB1+τ<yG1. Since g1(B″1)<0 and g2(B″1)<0, it has yˆB1<yB″1, then yG+1=yˆB1<yB″1=(1−b)yB1+τ<yG1, i.e., fIISOR1(G1)=yG+1−yG1<0.
On the other hand, since τ≤yG2, then G2∈N1. The trajectory starting from A1∈N1 with yA1=yG2 intersects with M2 at A−1, and then A−1 is mapped to A+1. Then it intersects with M1 at A+−1, and next A+−1 is mapped to A++1. Since yA+1=(1−b)yA−1+τ>(1−b)yG2+τ>yG2, so fIISOR1(A1)=yA++1−yA1=yA+−1−yG2>0.
The continuity of fIISOR1 implies that a point S∈¯A1G1⊂N1 exists so that fIISOR1(S)=0, i.e., an order-2 periodic solution exists (Figure 4(a)).
Similarly, a point S′∈¯A2B2⊂N2 exists so that fIISOR2(S′)=0, i.e., an order-2 periodic solution exists (Figure 4(b)).
Let z3(t)=(ϕ3(t),φ3(t)) (k−1)(T1+T2)≤t≤k(T1+T2), k∈N be the order-2 periodic solution. Denote
ϕ3(0)=h1+η,ϕ3(T1)=h2,ϕ3(T+1)=(1−a)h2,ϕ3(T1+T2)=h1,ϕ3((T1+T2)+)=h1+η, |
φ3(0)=δ3,φ3(T1)=δ4,φ3(T+1)=(1−b)δ4+τ,φ3(T1+T2)=δ3,φ3((T1+T2)+)=δ3. |
Define
Θ0≜δ4(1−b)δ4+τγ1γ3γ2γ4, |
where
γ1≜(1−h1−η)(h1+η−β)−αδ3,γ2≜(1−h2)(h2−β)−αδ4,γ3≜(1−(1−a)h2)((1−a)h2−β)−α((1−b)δ4+τ),γ4≜(1−h1)(h1−β)−αδ3. |
Theorem 9. For the model parameters with any one of C1)–C3), ρ>0 and D-4) h1<(1−a)h2<x∗1<h1+η<h2<x∗2, and byG2≤τ≤min{yG2,¯τ3}, if
μ3=Θ0exp{∫T1+T20+[ϕ3(t)(β+1−2ϕ3(t))−s1fφ3(t)(1+fφ3(t))2]dt}<1, |
then z3(t)=(ϕ3(t),φ3(t)) (k−1)(T1+T2)≤t≤k(T1+T2) is orbitally asymptotically stable.
Proof. For convenience, denote the intersection point between z3(t)=(ϕ3(t),φ3(t)) and N1 (M2, N2, ) by (, , ). Then, according to analogue of Poincaré Criterion, there are
and
Then
If , then is orbitally asymptotically stable.
For system (2.2) with model parameters , , , , , , , there is , then two interior equilibria exist in the system, and the phase diagram is presented in Figure 5(a); For the model parameters , , , , , , , there is , then a unique equilibrium exists in system (2.2), which is a sharp point. The phase diagram of the system (2.2) for such a case is presented in Figure 5(b). Moreover, system (2.2) undergoes a Bogdanov-Takens bifurcation of codimension two in a very small neighborhood of the unique interior equilibrium.
While for system (2.2) with model parameters , , , , , , , there is also , in such case a unique equilibrium exists in system (2.2), which is a saddle node, and the phase diagram is presented in Figure 5(c). System (2.2) undergoes a saddle-node bifurcation of codimension one in a very small neighborhood of the interior equilibrium.
For system (2.2) with model parameters , , , , , , there is . In such a case, system (2.2) does not have interior equilibrium, and the phase diagram of the system (2.2) is presented in Figure 5(d).
For system (2.2) with parameters , , , , , , the bifurcation diagrams of the residual dimension 1 for the system (2.2) are shown in Figure 6(a), (b) when is selected as the bifurcation parameter. The result shows that for larger , the interior equilibrium does not exist. When decreases to , a unique equilibrium exists in the system. As decreases below , the system undergoes a saddle-node bifurcation at the interior equilibrium , giving rise to two interior equilibria and .
For system (2.3) with model parameters , , , , , , , the control parameters: , , , , , , an order-1 periodic solution can be formed in both and (Theorem 6), as presented in Figure 7(a); For model parameters , , , , , , and control parameters , , , , , , an order-1 periodic solution can be formed in (Theorem 6), as presented in Figure 7(b); while for control parameters , , , , , , an order-1 periodic solution can be formed (Theorem 6), as presented in Figure 7(c).
For system (2.3) with given model parameters , , , , , , and control parameters , , , , , , an order-2 periodic solution can be formed in system (2.3) (Theorem 3.2), as presented in Figure 8(a); while for the control parameters , , , , , , a different order-2 periodic solution can be formed in system (2.3) (Theorem 8), as presented in Figure 8(b).
Considering that the Allee effect is an important mechanism in ecosystems and a realistic description of the interaction between species, we presented a model of prey-predator system with prey's Allee effect and generalist predator in the context of fish resources (Models (2.1) or (2.2)). We investigated the dynamic properties of Model (2.2) such as the type and stability of the boundary equilibria as well as the existence and stability of the interior equilibrium in detail (Theorems 1–3, Figure 5).
To show the influence of the parameters on the dynamics of Model (2.2), we analyzed the bifurcations in the predation system by selecting the capture rate of prey by predator and Allee threshold as key parameters. We showed that Model (2.2) will undergo a saddle-node bifurcation as changing of the capture rate (Figure 6), and undergo a Bogdanov-Takens bifurcation of codimension at least 2 and 3 as changing of .
To achieve sustainable and efficient exploitation of fish stocks, we adopted a bilateral intervention strategy, i.e., to avoid the distinction of prey fish caused by the Allee effect, releasing juvenile prey fish is adopted at a lower-level ; while for economic purposes, capturing both prey and predator fish is adopted at a higher-level . We obtained the conditions for the existence and stability of the order-1 periodic solution (Theorems 6, 7, Figure 7) and order-2 periodic solution (Theorems 8, 9, Figure 8) of the control system (2.3). The results showed that the extinction can be prevented by control even when the prey density is low, while in the case of the prey density increasing to a certain extent, fishing activities can be taken in a periodic way (periodic solution) to obtain the fish resources. Therefore, as long as the fish stocks are properly managed, the number of fish stocks can be controlled within an appropriate range, and the sustainable development and utilization of biological resources can be realized.
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
The work was supported by the National Natural Science Foundation of China (No. 11401068).
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
[1] | L. Page, S. Brin, R. Motwami, T. Winograd, The PageRank citation ranking: Bringing order to the web, Stanford Digital Library Technol. Proj., 1998. https://doi.org/10.1007/978-3-319-08789-4-10 |
[2] |
I. C. Ipsen, T. M. Selee, PageRank computation, with special attention to dangling nodes, SIAM J. Matrix Anal. Appl., 29 (2008), 1281–1296. https://doi.org/10.1137/060664331 doi: 10.1137/060664331
![]() |
[3] |
A. Langville, C. Meyer, A survey of eigenvector methods for web information retrieval, SIAM Rev., 47 (2005), 135–161. https://doi.org/10.1137/S0036144503424786 doi: 10.1137/S0036144503424786
![]() |
[4] | A. Langville, C. Meyer, Deeper inside PageRank, Internet Math., 1 (2004), 335–380. https://doi.org/10.1080/15427951.2004.10129091 |
[5] | G. H. Golub, C. F. Van Loan, Matrix Computations, 3 edition, The Johns Hopkins University Press, Baltimore, London, 1996. https://doi.org/10.1007/978-1-4612-5118-7-5 |
[6] | S. Kamvar, T. Haveliwala, C. Manning, G. Golub, Extrapolation methods for accelerating PageRank computations, in Proceedings of the Twelfth Internatinal World Wide Web Conference, ACM Press, New York, (2003), 261–270. https://doi.org/10.1145/775152.775190 |
[7] |
A. Sidi, Vector extrapolation methods with applications to solution of large systems of equations and to PageRank computations, Comput. Appl. Math., 56 (2008), 1–24. https://doi.org/10.1016/j.camwa.2007.11.027 doi: 10.1016/j.camwa.2007.11.027
![]() |
[8] |
X. Y. Tan, A new extrapolation method for PageRank computations, J. Comput. Appl. Math., 313 (2017), 383–392. https://doi.org/10.1016/j.cam.2016.08.034 doi: 10.1016/j.cam.2016.08.034
![]() |
[9] |
S. Cipolla, M. Redivo-Zaglia, F. Tudisco, Extrapolation methods for fixed-point multilinear PageRank computations, Numer. Linear Algebra Appl., 27 (2020), e2280. https://doi.org/10.1002/nla.2280 doi: 10.1002/nla.2280
![]() |
[10] |
D. Gleich, A. Gray, C. Greif, T. Lau, An inner-outer iteration for computing PageRank, SIAM J. Sci. Comput., 32 (2010), 349–371. https://doi.org/10.1137/080727397 doi: 10.1137/080727397
![]() |
[11] |
Z. Z. Bai, On convergence of the inner-outer iteration method for computing PageRank, Numer. Algebra Control Optim., 2 (2012), 855–862. https://doi.org/10.3934/naco.2012.2.855 doi: 10.3934/naco.2012.2.855
![]() |
[12] |
C. Q. Gu, F. Xie, K. Zhang, A two-step matrix splitting iteration for computing PageRank, J. Comput. Appl. Math., 278 (2015), 19–28. https://doi.org/10.1016/j.cam.2014.09.022 doi: 10.1016/j.cam.2014.09.022
![]() |
[13] |
C. Wen, T. Z. Huang, Z. L. Shen, A note on the two-step matrix splitting iteration for computing PageRank, J. Comput. Appl. Math., 315 (2017), 87–97. https://doi.org/10.1016/j.cam.2016.10.020 doi: 10.1016/j.cam.2016.10.020
![]() |
[14] |
Z. L. Tian, Y. Liu, Y. Zhang, Z. Y. Liu, M. Y. Tian, The general inner-outer iteration method based on regular splittings for the PageRank problem, Appl. Math. Comput., 356 (2019), 479–501. https://doi.org/10.1016/j.amc.2019.02.066 doi: 10.1016/j.amc.2019.02.066
![]() |
[15] |
J. F. Yin, G. J. Yin, M. Ng, On adaptively accelerated Arnoldi method for computing PageRank, Numer. Linear Algebra Appl., 19 (2012), 73–85. https://doi.org/10.1002/nla.789 doi: 10.1002/nla.789
![]() |
[16] |
C. Wen, Q. Y. Hu, G. J. Yin, X. M. Gu, Z. L. Shen, An adaptive Power-GArnoldi algorithm for computing PageRank, J. Comput. Appl. Math., 386 (2021), 113209. https://doi.org/10.1016/j.cam.2020.113209 doi: 10.1016/j.cam.2020.113209
![]() |
[17] |
C. Wen, Q. Y. Hu, B. Y. Pu, Y. Y. Huang, Acceleration of an adaptive generalized Arnoldi method for computing PageRank, AIMS Math., 6 (2021), 893–907. https://doi.org/10.3934/math.2021053 doi: 10.3934/math.2021053
![]() |
[18] |
H. D. Sterck, T. A. Manteuffel, S. F. McCormick, Q. Nguyen, J. Ruge, Multilevel adaptive aggregation for Markov chains, with application to web ranking, SIAM J. Sci. Comput., 30 (2008), 2235–2262. https://doi.org/10.1137/070685142 doi: 10.1137/070685142
![]() |
[19] |
Z. L. Shen, T. Z. Huang, B. Carpentieri, C. Wen, X. M. Gu, Block-accelerated aggregation multigrid for Markov chains with application to PageRank problems, Commun. Nonlinear. Sci. Numer. Simulat., 59 (2018), 472–487. https://doi.org/10.1016/j.cnsns.2017.11.031 doi: 10.1016/j.cnsns.2017.11.031
![]() |
[20] |
G. H. Golub, C. Greif, An Arnoldi-type algorithm for computing PageRank, BIT Numer. Math., 46 (2006), 759–771. https://doi.org/10.1007/s10543-006-0091-y doi: 10.1007/s10543-006-0091-y
![]() |
[21] |
Z. X. Jia, Refined iterative algorithms based on Arnoldi's process for large unsymmetric eigenproblems, Linear Algebra Appl., 259 (1997), 1–23. https://doi.org/10.1016/S0024-3795(96)00238-8 doi: 10.1016/S0024-3795(96)00238-8
![]() |
[22] |
G. Wu, Y. M. Wei, A Power-Arnoldi algorithm for computing PageRank, Numer. Linear Algebra Appl., 14 (2007), 521–546. https://doi.org/10.1002/nla.531 doi: 10.1002/nla.531
![]() |
[23] |
R. B. Morgan, M. Zeng, A harmonic restarted Arnoldi algorithm for calculating eigenvalues and determining multiplicity, Linear Algebra Appl., 415 (2006), 96–113. https://doi.org/10.1016/j.laa.2005.07.024 doi: 10.1016/j.laa.2005.07.024
![]() |
[24] |
Q. Y. Hu, C. Wen, T. Z. Huang, Z. L. Shen, X. M. Gu, A variant of the Power-Arnoldi algorithm for computing PageRank, J. Comput. Appl. Math., 381 (2021), 113034. https://doi.org/10.1016/j.cam.2020.113034 doi: 10.1016/j.cam.2020.113034
![]() |
[25] |
G. Wu, Y. M. Wei, An Arnoldi-Extrapolation algorithm for computing PageRank, J. Comput. Appl. Math., 234 (2010), 3196–3212. https://doi.org/10.1016/j.cam.2010.02.009 doi: 10.1016/j.cam.2010.02.009
![]() |
[26] |
H. F. Zhang, T. Z. Huang, C. Wen, Z. L. Shen, FOM accelerated by an extrapolation method for solving PageRank problems, J. Comput. Appl. Math., 296 (2016), 397–409. https://doi.org/10.1016/j.cam.2015.09.027 doi: 10.1016/j.cam.2015.09.027
![]() |
[27] |
C. Q. Gu, X. L. Jiang, C. C. Shao, Z. B. Chen, A GMRES-Power algorithm for computing PageRank problems, J. Comput. Appl. Math., 343 (2018), 113–123. https://doi.org/10.1016/j.cam.2018.03.017 doi: 10.1016/j.cam.2018.03.017
![]() |
[28] |
Y. Saad, M. H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput., 7 (1986), 857–869. https://doi.org/10.1137/0907058 doi: 10.1137/0907058
![]() |
[29] |
Z. L. Shen, T. Z. Huang, B. Carpentieri, X. M. Gu, C. Wen, An efficient elimination strategy for solving PageRank problems, Appl. Math. Comput., 298 (2017), 111–122. https://doi.org/10.1016/j.amc.2016.10.031 doi: 10.1016/j.amc.2016.10.031
![]() |
[30] |
Z. L. Shen, T. Z. Huang, B. Carpentieri, X. M. Gu, C. Wen, X. Y. Tan, Off-diagonal low-rank preconditioner for difficult PageRank problems, J. Comput. Appl. Math., 346 (2019), 456–470. https://doi.org/10.1016/j.cam.2018.07.015 doi: 10.1016/j.cam.2018.07.015
![]() |
[31] |
B. Y. Pu, T. Z. Huang, C. Wen, A preconditioned and extrapolation-accelerated GMRES method for PageRank, Appl. Math. Lett., 37 (2014), 95–100. https://doi.org/10.1016/j.aml.2014.05.017 doi: 10.1016/j.aml.2014.05.017
![]() |
[32] |
C. Q. Miao, X. Y. Tan, Accelerating the Arnoldi method via Chebyshev polynomials for computing PageRank, J. Comput. Appl. Math., 377 (2020), 112891. https://doi.org/10.1016/j.cam.2020.112891 doi: 10.1016/j.cam.2020.112891
![]() |
[33] | X. M. Gu, S. L. Lei, K. Zhang, Z. L. Shen, C. Wen, B. Carpentieri, A Hessenberg-type algorithm for computing PageRank problems, Numer. Algorithms, 2021. https://doi.org/10.1007/s11075-021-01175-w |
[34] |
Z. L. Shen, H. Yang, B. Carpentieri, X. M. Gu, C. Wen, A preconditioned variant of the refined Arnoldi method for computing PageRank eigenvectors, Symmetry, 13 (2021), 1327. https://doi.org/10.3390/sym13081327 doi: 10.3390/sym13081327
![]() |
[35] |
Z. L. Tian, Y. Zhang, J. X. Wang, C. Q. Gu, Several relaxed iteration methods for computing PageRank, J. Comput. Appl. Math., 388 (2021), 113295. https://doi.org/10.1016/j.cam.2020.113295 doi: 10.1016/j.cam.2020.113295
![]() |
[36] | Z. L. Tian, Z. Y. Liu, Y. H. Dong, The coupled iteration algorithms for computing PageRank, Numer. Algor., (2021), 1–15. https://doi.org/10.1007/s11075-021-01166-x |
[37] | Y. H. Feng, J. X. You, Y. X. Dong, An extrapolation iteration and its lumped type iteration for computing PageRank, Bull. Iran. Math. Soc., (2021), 1–8. https://doi.org/10.1007/s41980-021-00656-x |
[38] | Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, 2003. |
[39] | SuiteSparse Matrix Collection, Available from: https://sparse.tamu.edu/. |
[40] | T. Haveliwala, S. Kamvar, The second eigenvalue of the Google matrix, in Proceedings of the Twelfth International World Wide Web of Conference, 2003. |
[41] |
R. Horn, S. Serra-Capizzano, A general setting for the parametric Google matrix, Internet Math., 3 (2008), 385–411. https://doi.org/10.1080/15427951.2006.10129131 doi: 10.1080/15427951.2006.10129131
![]() |
1. | Xinrui Yan, Yuan Tian, Kaibiao Sun, Dynamic analysis of additional food provided non-smooth pest-natural enemy models based on nonlinear threshold control, 2025, 1598-5865, 10.1007/s12190-024-02318-7 | |
2. | Yuan Tian, Xinlu Tian, Xinrui Yan, Jie Zheng, Kaibiao Sun, Complex dynamics of non-smooth pest-natural enemy Gomportz models with a variable searching rate based on threshold control, 2025, 33, 2688-1594, 26, 10.3934/era.2025002 | |
3. | Xinrui Yan, Yuan Tian, Kaibiao Sun, Hybrid Effects of Cooperative Hunting and Inner Fear on the Dynamics of a Fishery Model With Additional Food Supplement, 2025, 0170-4214, 10.1002/mma.10805 |