
Let Fq be the finite field of q elements, and Fqn its extension of degree n. A normal basis of Fqn over Fq is a basis of the form {α,αq,⋯,αqn−1}. Some problems on normal bases can be finally reduced to the determination of the irreducible factors of the polynomial xn−1 in Fq, while the latter is closely related to the cyclotomic polynomials. Denote by F(xn−1) the set of all distinct monic irreducible factors of xn−1 in Fq. The criteria for
|F(xn−1)|≤2
have been studied in the literature. In this paper, we provide the sufficient and necessary conditions for
|F(xn−1)|=s,
where s is a positive integer by using the properties of cyclotomic polynomials and results from the Diophantine equations. As an application, we obtain the sufficient and necessary conditions for
|F(xn−1)|=3,4,5.
Citation: Weitao Xie, Jiayu Zhang, Wei Cao. On the number of the irreducible factors of xn−1 over finite fields[J]. AIMS Mathematics, 2024, 9(9): 23468-23488. doi: 10.3934/math.20241141
[1] | Serap Özcan, Saad Ihsan Butt, Sanja Tipurić-Spužević, Bandar Bin Mohsin . Construction of new fractional inequalities via generalized n-fractional polynomial s-type convexity. AIMS Mathematics, 2024, 9(9): 23924-23944. doi: 10.3934/math.20241163 |
[2] | Hasan Kara, Hüseyin Budak, Mehmet Eyüp Kiriş . On Fejer type inequalities for co-ordinated hyperbolic ρ-convex functions. AIMS Mathematics, 2020, 5(5): 4681-4701. doi: 10.3934/math.2020300 |
[3] | Wei Gao, Artion Kashuri, Saad Ihsan Butt, Muhammad Tariq, Adnan Aslam, Muhammad Nadeem . New inequalities via n-polynomial harmonically exponential type convex functions. AIMS Mathematics, 2020, 5(6): 6856-6873. doi: 10.3934/math.2020440 |
[4] | Naila Mehreen, Matloob Anwar . Ostrowski type inequalities via some exponentially convex functions with applications. AIMS Mathematics, 2020, 5(2): 1476-1483. doi: 10.3934/math.2020101 |
[5] | Tekin Toplu, Mahir Kadakal, İmdat İşcan . On n-Polynomial convexity and some related inequalities. AIMS Mathematics, 2020, 5(2): 1304-1318. doi: 10.3934/math.2020089 |
[6] | Ghulam Farid, Saira Bano Akbar, Shafiq Ur Rehman, Josip Pečarić . Boundedness of fractional integral operators containing Mittag-Leffler functions via (s,m)-convexity. AIMS Mathematics, 2020, 5(2): 966-978. doi: 10.3934/math.2020067 |
[7] | M. Emin Özdemir, Saad I. Butt, Bahtiyar Bayraktar, Jamshed Nasir . Several integral inequalities for (α, s,m)-convex functions. AIMS Mathematics, 2020, 5(4): 3906-3921. doi: 10.3934/math.2020253 |
[8] | Anjum Mustafa Khan Abbasi, Matloob Anwar . Ostrowski type inequalities for exponentially s-convex functions on time scale. AIMS Mathematics, 2022, 7(3): 4700-4710. doi: 10.3934/math.2022261 |
[9] | Saad Ihsan Butt, Ahmet Ocak Akdemir, Muhammad Nadeem, Nabil Mlaiki, İşcan İmdat, Thabet Abdeljawad . (m,n)-Harmonically polynomial convex functions and some Hadamard type inequalities on the co-ordinates. AIMS Mathematics, 2021, 6(5): 4677-4690. doi: 10.3934/math.2021275 |
[10] | Andrea Aglić Aljinović, Domagoj Kovačević, Mehmet Kunt, Mate Puljiz . Correction: Quantum Montgomery identity and quantum estimates of Ostrowski type inequalities. AIMS Mathematics, 2021, 6(2): 1880-1888. doi: 10.3934/math.2021114 |
Let Fq be the finite field of q elements, and Fqn its extension of degree n. A normal basis of Fqn over Fq is a basis of the form {α,αq,⋯,αqn−1}. Some problems on normal bases can be finally reduced to the determination of the irreducible factors of the polynomial xn−1 in Fq, while the latter is closely related to the cyclotomic polynomials. Denote by F(xn−1) the set of all distinct monic irreducible factors of xn−1 in Fq. The criteria for
|F(xn−1)|≤2
have been studied in the literature. In this paper, we provide the sufficient and necessary conditions for
|F(xn−1)|=s,
where s is a positive integer by using the properties of cyclotomic polynomials and results from the Diophantine equations. As an application, we obtain the sufficient and necessary conditions for
|F(xn−1)|=3,4,5.
Plankton is divided into two groups: zooplankton and phytoplankton. Phytoplankton are the primary producers in aquatic ecological models, as well as the main supplier of dissolved oxygen to phytoplankton blooms. Phytoplankton opens up the food web of aquatic ecosystems. Zooplankton, as economic aquatic animals, constitute an important feed for fish and other economic animals in the middle and upper waters, which is of great significance to the development of fishery[1,2,3,4].
Plankton has been studied extensively in many ways. In References[1,2,3,4], the researchers found that zooplankton eat plankton and zooplankton smaller than themselves, or feed on algae, bacteria, copepods and other food scraps. Therefore, there is a predator-prey relationship between zooplankton and phytoplankton. In 1939, Fleming[5] presented the first mathematical model of plankton. Since then, researchers have done a lot of work on plankton[6,7,8,9,10], focusing on factors such as nutrients, temperature, light, viral diseases and harvest to understand the bloom and disappearance of algal blooms. However, the researchers found that toxins released by toxic-producing phytoplankton have an effect on the termination of plankton blooms, which means that toxic chemicals can act as biological control for other plankton populations[11].
When the predator captures the prey, the prey seeks shelter because of a survival instinct. Nature provides shelter to the prey, and this behavior keeps the balance of the predator-prey model[12,13,14]. For lake ecosystems, prey refuge can stabilize plankton biomass by preventing phytoplankton from being temporarily eaten by zooplankton. Scholars have found that phytoplankton shelters can be obtained through benthic sediments, which can allow phytoplankton to temporarily escape from zooplankton predation. At the same time, the water layer can also form a temporary shelter for phytoplankton, and the shelter can prevent the extinction of the prey population[15,16,17]. Li et al.[18] proposed a model with the refuge as follows:
{dPdt=rP(1−PK)−β1(P−m)Za1+(P−m),dZdt=β2(P−m)Za1+(P−m)−dZ−θPZa2+P, | (1.1) |
where P and Z represent the number of phytoplankton and zooplankton, respectively. r represents the intrinsic growth rate of phytoplankton, K represents the environmental carrying capacity of phytoplankton, β1 represents the predation rate, β2 represents the conversion rate, d represents the mortality rate of zooplankton, θ represents the toxin release rate, a1 and a2 represent half-full constants, m represents the number of protected phytoplankton when phytoplankton have the ability to shelter and P−m denotes the number of unprotected phytoplankton that can be preyed upon by zooplankton. They[18] studied the effect of refuge of phytoplankton on the phytoplankton - zooplankton model.
In nature, the population diffuse from one area to another in order to survive. The predator-prey model with diffusion can generate complex spatial patterns [19,20,21,22]. There are two types of diffusion: self-diffusion and cross - diffusion. The former refers to the diffusion of one species from the higher - density area to the lower - density area in order to survive, while the latter refers to the diffusion of one species as influenced by other species. Meanwhile, time delay in nature is an important factor affecting the predator - prey model. A dynamical model without time delay can only be an approximation [23]. Generally speaking, there are many kinds of delayed factors in the growth process, such as the digestion of food [24], the maturation of cells [25], pregnancy[26,27]. These processes are not instantaneous and take time to complete. Later, it is found that the existence of time delay would make the positive steady state of a predator - prey model lose stability, resulting in bifurcation or periodic oscillation[21,28]. Zhao et al.[19] proposed a reaction - diffusion model with mature delay:
{∂P∂t=d1ΔP+rP(1−PK)−μPZα+P,∂Z∂t=d2ΔZ+μ1PZα+P−δZ−ρP(t−τ)Zα+P(t−τ), | (1.2) |
where d1 and d2 represent the diffusion coefficients. τ is the time required for phytoplankton to form a mature cell to release toxins. They[19] analyzed the stability of the equilibrium, the existence and properties of Hopf bifurcation. They[19] also found that time delay has an effect on the model (1.2). Recently, Hopf bifurcation has also continued to be investigated in fractional-order dynamical systems [29,30,31,32] and integer-order differential systems [33,34,35,36,37].
In lake ecosystems, fish eat plankton for survival, while zooplankton eat phytoplankton. Thus, phytoplankton, zooplankton and fish form a food chain. There is a lot of work that has done by many researchers[20,38,39].
The functional response reflects the predator's predation on the prey, which is an important part of the predator-prey model. It can be divided into prey-dependent (Holling I-IV type[40], Ivlev type [41], Rosenzweig type [22]) and predator-dependent functional responses (Beddington-DeAngelis type [42,43], Crowley-Martin (C-M) type [44], Hassell-Varley type [45], ratio-dependent type[46]). A predator-prey model that takes into account interactions between predators is more realistic. In [44], Crowley and Martin first proposed the C-M functional response:
F(H,P)=aH(1+abH)(1+cP), | (1.3) |
where F(H,P) represents the predation rate per predator, H represents the density of the prey per unit of area, P represents the density of the predator per unit of area, a represents the attack coefficient, b is the handling time and c is the interference coefficient. This functional response indicates that there is interference between predators as they feed on and handle prey. In [47], the authors proposed a model with the C-M functional response and showed that the system has complex dynamical behavior. In [48], the interaction between mature prey and predator is assumed to be the C-M functional response; the authors analyzed the positivity, boundedness and existence of equilibrium points. They not only analyzed the stability behavior of the delayed and non-delayed system, but also discussed the properties of Hopf bifurcation by choosing delay as the bifurcation parameter. In recent years, the researchers have done a great deal of work on this predator-prey models with a C-M functional response[49,50,51]. Therefore, in this paper, on the basis of Eq (1.3), we will consider that fish predation on zooplankton follows the C-M functional response:
p(Z,F)=γZF(1+aZ)(1+bF), | (1.4) |
where Z and F represent the number of zooplankton and fish, respectively. a represents interference between zooplankton, b is the interference between fish and γ is the maximum rate of fish predation on zooplankton. For different values of a and b, we have that (i) Eq (1.4) has turned into a Holling II functional response if a>0 and b=0; (ii) Eq (1.4) expresses a saturation response with respect to a predator if a=0 and b>0; (iii) Eq (1.4) becomes a linear mass-action functional response if a=0 and b=0.
In this paper, we consider the self-diffusion of populations. Let P(x,t), Z(x,t) and F(x,t) be the densities of phytoplankton, zooplankton and fish at the location x and the time t, respectively. We also give the following assumptions.
(1) The phytoplankton population follows the logistic growth under the condition of no zooplankton.
(2) Zooplankton preyed upon by fish is a C-M functional response.
(3) Phytoplankton preyed upon by zooplankton is a Holling II functional response. Nature provides shelter to the prey, so a constant proportion m∈(0,1) of the phytoplankton take refuge, leaving (1−m)P of the unprotected phytoplankton available for zooplankton grazing, following the Holling II functional response, i.e., f1(P)=(1−m)Pa1+(1−m)P.
(4) The progress of release of toxins takes the Holling II functional response and considers the mature delay of toxins in the cell, i.e., f2(P)=P(t−τ)a2+P(t−τ).
Based on the above assumptions, a schematic diagram that expresses the interactions of phytoplankton, zooplankton and fish is depicted in Figure 1. The corresponding model is
{∂P(x,t)∂t=d1ΔP+r1P(1−PK)−β1(1−m)PZa1+(1−m)P,x∈(0,lπ),t>0,∂Z(x.t)∂t=d2ΔZ+β2(1−m)PZa1+(1−m)P−δP(t−τ)Za2+P(t−τ)−γ1ZF(1+aZ)(1+bF)−g1Z,x∈(0,lπ),t>0,∂F(x,t)∂t=d3ΔF+γ2ZF(1+aZ)(1+bF)−g2F,x∈(0,lπ),t>0,Px(0,t)=Zx(0,t)=Fx(0,t)=0,t>0,Px(lπ,t)=Zx(lπ,t)=Fx(lπ,t)=0,t>0,P(x,t)>0,Z(x,t)>0,F(x,t)>0,x∈[0,lπ],t∈[−τ,0], | (1.5) |
where d1, d2 and d3 represent the diffusion coefficients of each population, respectively; Δ is the Laplace operator, β1 represents the maximum predation rate of phytoplankton by zooplankton, β2 is the conversion rate, K represents the environmental carrying capacity, γ1 represents the predation rate, γ2 represents the conversion rate, m is the refuge proportion of phytoplankton, α1 and α2 are half saturation constants, δ represents the release rate of phytoplankton toxins and g1 and g2 represent the natural mortality rates of zooplankton and fish, respectively. a represents the degree of interference between zooplankton, b represents the degree of interference between fish and τ represents the mature time delay of phytoplankton toxin release. The Neumann boundary condition indicates that the area is closed and no individuals can move across this area.
We rescale the model (1.5) by
˜P=PK,˜Z=Z,˜F=F,˜t=r1t,~d1=d1r1,~β1=β1Kr1,~α1=a1K,˜a=1a,˜b=1b,~d2=d2r1,~β2=β2r1,~α2=a2K,˜δ=δr1,~γ1=γ1abr1,~d3=d3r1,~γ2=γ2abr1,~g1=g1r1,~g2=g2r1. |
For the sake of convenience, omitting the breaking line, the model (1.5) becomes
{∂P∂t=d1ΔP+P(1−P)−β1(1−m)PZα1+(1−m)P,x∈(0,lπ),t>0,∂Z∂t=d2ΔZ+β2(1−m)PZα1+(1−m)P−δP(t−τ)Zα2+P(t−τ)−γ1ZF(a+Z)(b+F)−g1Z,x∈(0,lπ),t>0,∂F∂t=d3ΔF+γ2ZF(a+Z)(b+F)−g2F,x∈(0,lπ),t>0,Px(0,t)=Zx(0,t)=Fx(0,t)=0,t>0,Px(lπ,t)=Zx(lπ,t)=Fx(lπ,t)=0,t>0,P(x,t)>0,Z(x,t)>0,F(x,t)>0,x∈[0,lπ],t∈[−τ,0], | (1.6) |
where 0<P<1, 0<m<1 and all parameters are positive.
Our paper is organized as follows. The existence and stability of equilibrium of the model (2.1) are discussed in Section 2. Meanwhile, the occurrence of Hopf bifurcation is given by choosing m as a bifurcation parameter. The existence and properties of Hopf bifurcation of the model (3.1) are discussed in Section 3. In Section 4, Hopf bifurcation of the reaction-diffusion model (1.6) at the positive equilibrium and its properties are analyzed. In Section 5, a numerical simulation is demonstrated to prove the previous theoretical results by using Matlab software. In the last section, we have a brief discussion.
In this section, we will investigate the dynamical behavior of the model (1.6) with no delay and no diffusion. That is, the model is
{dPdt=P(1−P)−β1(1−m)PZα1+(1−m)P,dZdt=β2(1−m)PZα1+(1−m)P−δPZα2+P−γ1ZF(a+Z)(b+F)−g1Z,dFdt=γ2ZF(a+Z)(b+F)−g2F. | (2.1) |
According to the existence theorem of the solution of ordinary differential equations, we can know that the solution of the model (2.1) exists. By Lemma 2.1 in Reference[52], we give the following lemma to explain the positivity of the solution of the model (2.1).
Lemma 2.1. All solutions of the model (2.1) that start positive remain positive.
Proof. Let (P(t),Z(t),F(t)) be any solution of the model (2.1). Assuming that the initial time is t0 and one solution of the model (2.1) is at least not positive, then we have the following three cases:
(i) there exists time t1 such that P(t0)>0,P(t1)=0,P′(t1)≤0,Z(t)>0,F(t)>0,t0<t<t1;
(ii) there exists time t2 such that Z(t0)>0,Z(t2)=0,Z′(t2)≤0,P(t)>0,F(t)>0,t0<t<t2;
(iii) there exists time t3 such that F(t0)>0,F(t3)=0,F′(t3)≤0,P(t)>0,Z(t)>0,t0<t<t3.
If the first case is true, then we get P′(t1)=0. This contradicts with P′(t1)≤0. Similarly, we have that Z′(t2)=0, which contradicts with Z′(t2)≤0. And, F′(t3)=0, which contradicts with F′(t3)≤0.
Because of the arbitrariness of P(t),Z(t) and F(t), all solutions of the model (2.1) remain positive for all t>t0. Thus, all solutions of the model (2.1) that start positive remain positive.
Obviously, the model (2.1) has three boundary equilibria: one trivial equilibrium E0(0,0,0) and two boundary equilibria E1(1,0,0) and E2(P2,Z2,0) if (H1):β2−δ−g1>0 holds, where
P2=−B+√B2+4(1−m)(β2−δ−g1)g1α1α22(1−m)(β2−δ−g1),Z2=(1−P2)[α1+(1−m)P2](1−m)β1; |
here, B=(β2−g1)(1−m)α2−(δ+g1)α.
Assume that the model (2.1) has the coexistence equilibrium E∗(P∗,Z∗,F∗), where P∗, Z∗ and F∗ satisfy
{P∗(1−P∗)−β1(1−m)P∗Z∗α1+(1−m)P∗=0,β2(1−m)P∗Z∗α1+(1−m)P∗−δP∗Z∗α2+P∗−γ1Z∗F∗(a+Z∗)(b+F∗)−g1Z∗=0,γ2Z∗F∗(a+Z∗)(b+F∗)−g2F∗=0. | (2.2) |
By the first equation of (2.2), we can obtain
Z∗=(1−P∗)[α1+(1−m)P∗](1−m)β1. | (2.3) |
Substituting Eq (2.3) into the third equation of (2.2), we have
F∗=(γ2−bg2)(1−P∗)[α1+(1−m)P∗]−bg2aβ1(1−m)g2[aβ1(1−m)+(1−P∗)(α1+(1−m)P∗)]. | (2.4) |
If (H2):(γ2−bg2)(1−P∗)[α1+(1−m)P∗]−bg2aβ1(1−m)>0 holds, then we have that F∗>0.
Substituting Eqs (2.3) and (2.4) into the second equation of (2.2), P∗ is the positive root of the equation
f(P)=M5P5+M4P4+M3P3+M2P2+M1P+M0=0, | (2.5) |
where
M5=(1−m)2γ2(β2−δ−g1),M4=−γ2(1−m)(1−m−α1)(β2−δ−g1)+δγ2(1−m)(1−m−α1)−γ2(1−m)2(1−α2)(β2−g1)−g1γ2α1(1−m),M3=−[aβ1(1−m)+α1]γ2(1−m)(β2−δ−g1)+(β2−g1γ2(1−m)(1−α2)(1−m−α1)−δγ2(1−m−α1)2+g1γ2α1(1−m−α1)+γ1(1−m)2(γ2−bg2)−α2β2γ2(1−m)2+γ2(1−m)α1(δ+g1)+g1γ2α2(1−m)(1−m−α1),M2=α1(g1γ2α2+γ1β1)(1−m)+α2β2γ2(1−m)(1−m−α2)+(1−m)(1−α2)α1γ2(β2−g1)−γ2α1(1−m−α1)(2δ+g1)+aβ1g1γ2α1(1−m)+g1γ2α2−aβ1δγ2(1−m)(1−m−α1)+β1(1−m)2(1−α2)[aβ2γ2−ag1γ2−γ1(γ2−bg2)],M1=[aβ1(1−m)+α1][γ2α2(β2−g1)+α1g1γ2α2−γ2α1(δ+g1)]−(1−m−α1)g1γ2α1α2+γ1abg2β21(1−m)2−γ1β1(1−m)(γ2−bg2)[α2(1−m−α1)+α1],M0=γ1abg2β21α2(1−m)2−[aβ1(1−m)+α1]g1γ2α1α2−γ1β1(1−m)(γ2−bg2)α1α2. |
Under the condition (H1), we have that M5>0. By the Descarte's rule of signs[53], Eq (2.5) has only one positive root P∗ if and only if one of the following terms is satisfied:
(1)M4>0,M3>0,M2>0,M1>0,M0<0;(2)M4>0,M3>0,M2>0,M1<0,M0<0;(3)M4>0,M3>0,M2<0,M1<0,M0<0;(4)M4>0,M3<0,M2<0,M1<0,M0<0;(5)M4<0,M3<0,M2<0,M1<0,M0<0. |
Here, we give two figures as examples of the first two cases to verify the conclusion that there is only one positive root (see Figure 2). We first give the assumption (H3): one of the conditions (1)–(5) is true. Therefore, the following conclusion can be obtained.
Theorem 2.1. Under the conditions (H1), (H2) and (H3), the model (2.1) has only one positive equilibrium E∗(P∗,Z∗,F∗), which is determined by Eqs (2.3)–(2.5).
The stability of all equilibria will be analyzed in this part.
The Jacobian matrix of the model (2.1) is
A=(a11a120a21a22a230a32a33), | (2.6) |
where
a11=1−2P−α1β1(1−m)Z[α1+(1−m)P]2,a12=−β1(1−m)Pα1+(1−m)P,a21=α1β2(1−m)Z[α1+(1−m)P]2−δα2Z(α2+P)2,a23=−γ1bZ(a+Z)(b+F)2,a22=β2(1−m)Pα1+(1−m)P−δPα2+P−γ1aF(a+Z)2(b+F)−g1,a32=γ2aF(a+Z)2(b+F),a33=γ2bZ(a+Z)(b+F)2−g2. |
The characteristic equation of the model (2.1) is
λ3−(a11+a22+a33)λ2+(a11a22+a11a33+a22a33)λ+a11a32a23−a11a22a33+a12a21a33=0. | (2.7) |
According to Eq (2.7), we have the following conclusions.
Theorem 2.2. The boundary equilibrium E0(0,0,0) of the model (2.1) is always unstable.
Proof. The characteristic equation of the model (2.1) at E0 is
(λ−1)(λ+g1)(λ+g2)=0. |
It has three roots:
λ1=1>0,λ2=−g1<0,λ3=−g2<0. |
Thus, the boundary equilibrium E0 is unstable.
Theorem 2.3. If (H4):β2(1−m)α1+1−m−δα2+1−g1<0 holds, then the boundary equilibrium E1(1,0,0) of the model (2.1) is locally asymptotically stable.
Proof. The characteristic equation of the model (2.1) at E1 is
(λ+1)[λ−(β2(1−m)α1+(1−m)−δα2+1−g1)](λ+g2)=0. |
It has three roots:
λ1=−1<0,λ2=β2(1−m)α1+(1−m)−δα2+1−g1,λ3=−g2<0. |
When (H4) is true, then λ2<0. Therefore, the boundary equilibrium E1 is locally asymptotically stable under the condition (H4).
The characteristic equation of the model (2.1) at E2(P2,Z2,0) is
(λ−s1)(λ2+s2λ+s3)=0, |
where
s1=γ2(1−P2)[α1+(1−m)P2]abβ1(1−m)+b(1−P2)[α1+(1−m)P2]−g2,s2=α1−β1(1−m)Z2+(1−m)P2−2β2(1−m)P2α1+(1−m)P2,s3=−1(α2+P2)[α1+(1−m)P2]2[(1−2P2)(α1+(1−m)P2)−α1(1−P2)][β2(1−m)P2(α2+P2)−(α1+(1−m)P2)((δ+g1)P2+α2g1)]+α1P2(1−(1−m)P2)[α1+(1−m)P2]2−δα2P2(1−P2)(α+P2)2. |
Assume that the condition (H5):s1<0, s22−4s3>0, s2<0 and s3>0 are true; we can get the following result.
Theorem 2.4. Under the assumptions (H1) and (H5), the boundary equilibrium E2(P2,Z2,0) of the model (2.1) is locally asymptotically stable.
The characteristic equation of the model (2.1) at E∗(P∗,Z∗,F∗) is
λ3+M6λ2+M7λ+M8=0, | (2.8) |
where
M6=−(A11+A22+A33),M7=A22A33+A11A33+A11A22−A12A21−A23A32,M8=A11A23A32+A12A21A33−A11A22A33,A11=−P∗+β1(1−m)2Z∗[α1+(1−m)P∗]2,A12=−β1(1−m)P∗α1+(1−m)P∗,A21=α1β2(1−m)Z∗[α1+(1−m)P∗]2−δα2Z∗(α2+P∗)2,A22=γ1F∗Z∗(a+Z∗)2(b+F∗),A23=−γ1bZ∗(a+Z∗)(b+F∗)2,A32=γ2aF∗(a+Z∗)2(b+F∗),A33=−γ2F∗Z∗(a+Z∗)(b+F∗)2. |
From the Routh-Hurwitz criterion[54], if M6>0, M7>0 and M6M7−M8>0, then all solutions of Eq (2.8) have negative real parts. When M6>0, M7>0 and M6M7−M8<0, Eq (2.8) has one negative root and a pair of complex roots with a positive real part. Assume that (H6): M6>0, M7>0 and M6M7−M8>0; then, the stability of the positive equilibrium E∗ will be obtained.
Theorem 2.5. Suppose that the conditions (H1)–(H3) are true. If the condition (H6) holds, then E∗ is locally asymptotically stable. Further, E∗ loses stability when M6M7−M8 passes through 0; in other words, the model (2.1) undergoes a Hopf bifurcation at E∗ when M6M7−M8=0.
Next, we will choose m as the bifurcation parameter to study the occurrence of Hopf bifurcation of the model (2.1) at E∗. By using the results in Reference[55], the following result can be obtained.
Theorem 2.6. If the characteristic equation of the model (2.1) at E∗(P∗,Z∗,F∗) is
λ3+M6(m)λ2+M7(m)λ+M8(m)=0, |
where M6(m), T(m)=M6(m)M7(m)−M8(m) and M8(m) are the smooth functions of m and there exists a positive number m∗ that satisfies
(1) M6(m∗)>0, T(m∗)=0 and M8(m∗)>0;
(2) dTdm|m=m∗≠0,
then Hopf bifurcation occurs at E∗(P∗,Z∗,F∗) when m=m∗.
We used Matlab software for numerical simulations to obtain this result in Section 5. Meanwhile, we can also choose δ as a bifurcation parameter to study the occurrence of Hopf bifurcation of the model (2.1) at E∗(P∗,Z∗,F∗). Since the discussion process is similar, we will only give the bifurcation diagram in Section 5.
Under some conditions, Hopf bifurcation may not take place. Thus, we will discuss the global asymptotical stability of the positive equilibrium E∗ as follows.
Theorem 2.7. Suppose that the conditions (H1)–(H3) are true. The positive equilibrium E∗ of the model (2.1) is globally asymptotically stable if 1−β1(1−m)2Z∗α1[α1+(1−m)P∗]−δ2[α2+P∗]2>0 and [α1+(1−m)P∗]β1P∗δβ2α1(α2+P∗)−2−β21(1−m)P∗γ2α1g1ab>0.
Proof. Let (P,Z,F) be any positive solution of the model (2.1). Define a Lyapunov function
V(t)=P−P∗−P∗lnPP∗+[α1+(1−m)P∗]β1β2α1(Z−Z∗−Z∗lnZZ∗)+F−F∗−F∗lnFF∗. |
Calculating the derivative of V(t) along the solution of the model (2.1), then we have
dV(t)dt=(P−P∗){−(P−P∗)+β1(1−m)(1−m)Z∗(P−P∗)−[α1+(1−m)P∗](Z−Z∗)[α1+(1−m)P∗][α1+(1−m)P]}+[α1+(1−m)P∗]β1β2α1(Z−Z∗){β2(1−m)α1[α1+(1−m)P∗][α1+(1−m)P](P−P∗)+β2(1−m)P∗Z[α1+(1−m)P∗](Z−Z∗)−δα2(α2+P)(α2+P∗)(P−P∗)−P∗δα2+P∗(Z−Z∗)−γ1b(a+Z)(b+F)(b+F∗)(F−F∗)−γ1aF∗Z(a+Z)(a+Z∗)(b+F∗)(Z−Z∗)−g1Z(Z−Z∗)}+(F−F∗){γ2a(a+Z)(b+F)(a+Z∗)(Z−Z∗)−γ2Z∗(b+F)(a+Z∗)(b+F∗)(F−F∗)}≤−{1−β1(1−m)2Z∗[α1+(1−m)P][α1+(1−m)P∗]−δ2α22[α2+P]2[α2+P∗]2}(P−P∗)2−{−2−β1(1−m)P∗α1Z+[α1+(1−m)P∗]β1P∗δβ2α1(α2+P∗)+γ1aF∗Z(a+Z)(a+Z∗)(b+F∗)+g1Z}(Z−Z∗)2−{γ2Z∗(b+F)(a+Z∗)(b+F∗)+γ22a2(a+Z)2(b+F)2(a+Z∗)2}(F−F∗)2≤−{1−β1(1−m)2Z∗α1[α1+(1−m)P∗]−δ2[α2+P∗]2}(P−P∗)2−{−2−β21(1−m)P∗γ2α1g1ab+[α1+(1−m)P∗]β1P∗δβ2α1(α2+P∗)}(Z−Z∗)2−{γ2Z∗(b+F)(a+Z∗)(b+F∗)+γ22a2(a+Z)2(b+F)2(a+Z∗)2}(F−F∗)2. |
Here, it is obvious that γ2Z∗(b+F)(a+Z∗)(b+F∗)+γ22a2(a+Z)2(b+F)2(a+Z∗)2>0. If 1−β1(1−m)2Z∗α1[α1+(1−m)P∗]−δ2[α2+P∗]2>0 and −2−β21(1−m)P∗γ2α1g1ab+[α1+(1−m)P∗]β1P∗δβ2α1(α2+P∗)>0 hold, then we have that the coefficients of (P−P∗)2, (Z−Z∗)2 and (F−F∗)2 are always negative. Thus, dV(t)dt is negative. This completes the proof.
Here, under the conditions (H1)–(H3), we choose delay τ as the bifurcation parameter and study its influence on the stability of the positive equilibrium E∗(P∗,Z∗,F∗). The model (2.1) is
{dPdt=P(1−P)−β1(1−m)PZα1+(1−m)P,dZdt=β2(1−m)PZα1+(1−m)P−δP(t−τ)Zα2+P(t−τ)−γ1ZF(a+Z)(b+F)−g1Z,dFdt=γ2ZF(a+Z)(b+F)−g2F. | (3.1) |
Let u1(t)=P(t)−P∗, u2(t)=Z(t)−Z∗, u3(t)=F(t)−F∗ and u(t)=(u1(t),u2(t),u3(t))T ∈R3. The linearized system of the model (3.1) at E∗ is
{du1(t)dt=b11u1(t)+b12u2(t),du2(t)dt=b21u1(t)+b22u2(t)+b23u3(t)+cu1(t−τ),du3(t)dt=b32u2(t)+b33u3(t), | (3.2) |
where
b11=A11,b12=A12,b21=α1β2Z∗(1−m)[α1+P∗(1−m)]2,b22=A22,b23=A23,b32=A32,b33=A33,c=−δα2Z∗(α2+P∗)2. |
Then, the model (3.2) can also be given by
du(t)dt=L1u(t)+L2u(t−τ), | (3.3) |
where
L1=(b11b120b21b22b230b32b33),L2=(000c00000). | (3.4) |
Thus, O(0, 0, 0) is the zero equilibrium of the model (3.2). The characteristic equation of the model (3.2) at O is
λ3+D2λ2+D1λ+D0+e−λτ(D3λ+D4)=0, | (3.5) |
where
D2=−(b11+b22+b33),D1=b11b22+b11b33+b22b33−b12b21−b23b32,D0=b11b23b32+b12b21b33−b11b22b33,D3=−cb12,D4=cb12b33. |
The roots of Eq (3.5) have been discussed above when τ=0. Next, we will study the effect of delay τ(τ>0) on the model (3.2).
Suppose that λ=iω(ω>0) is a pair of pure imaginary roots of Eq (3.5). Substituting it into Eq (3.5), we can obtain
−iω3−D2ω3+D1iω+D0+(cosωτ−isinωτ)(iD3ω+D4)=0. | (3.6) |
Separating the real and imaginary parts of Eq (3.6), we can obtain
{ω3−D1ω=D3ωcosωτ−D4sinωτ,D2ω2−D0=D3ωsinωτ+D4cosωτ. | (3.7) |
From Eq (3.7), we can get
ω6+ω4(D22−2D1)+ω2(D21−2D0D2−D33)+D20−D24=0. | (3.8) |
Let z=ω2, D5=D22−2D1, D6=D21−2D0D1−D23 and D7=D20−D24. Then, Eq (3.8) can be rewritten as
f(z)=z3+D5z2+D6z+D7=0, | (3.9) |
and we have
f′(z)=3z2+2D5z+D6=0. | (3.10) |
If Eq (3.9) has at least one positive root, then Hopf bifurcation takes place. Assume that (H7):Δ1=D25−3D6≤0 and (H8):Δ1=D25−3D6>0. By Lemmas 2.2 and 4.2 in Reference [56], we can get the following results.
Since limz→+∞f(z)=+∞, Eq (3.9) has at least one positive root when D7<0.
If (H7) holds, then f(z) is monotonically increasing for z∈[0,+∞); so, when D7≥0 and (H7) hold, Eq (3.9) has no positive root for z∈[0,+∞).
When D7≥0 and (H8) hold, Eq (3.10) has two roots, that is, z∗1 and z∗2, where
z∗1=−D5+√Δ13,z∗2=−D5−√Δ13. |
Furthermore, we have
f″(z∗1)=−2D5+2√Δ1+2D5=2√Δ1>0,f″(z∗2)=−2D5−2√Δ1+2D5=−2√Δ1<0. |
Therefore, we can obtain z∗1 and z∗2 as the local minimum and the local maximum of f(z), respectively.
Theorem 3.1. For Eq (3.9), the following conclusions are true.
(1) If D7<0, then Eq (3.9) has at least one positive root.
(2) If the condition (H7) holds and D7≥0, then Eq (3.9) has no positive root.
(3) If the condition (H8) holds and D7≥0, then Eq (3.9) has two positive roots when z∗1>0 and f(z∗1)≤0.
Without loss of generality, suppose that Eq (3.9) has three positive roots defined by z1, z2 and z3, respectively. Thus, Eq (3.8) has three positive roots ω1=√z1, ω2=√z2 and ω3=√z3. From Eq (3.7), we can get
τjk=1ωkarccos[D3ω4k+(D2D4−D1D3)ω2k−D0D4D23ω2k+D24]+2jπωk,k=1,2,3,j=0,1,2,⋯; | (3.11) |
thus, ±iωk is a pair of purely imaginary roots of Eq (3.5) when τ=τjk.
Define
τ0=mink=1,2,3.j∈N0{τjk},ω0=ωk|τ=τ0. |
Let λ(τ)=α(τ)+iβ(τ) be the root of Eq (3.5) near τ=τjk, and let it satisfy α(τjk)=0 and β(τjk)=ωk, j=0,1,2,⋯, k=1,2,3.
Theorem 3.2. Assume that zk=ω2k and f′(zk)≠0; then, dReλ(τ)dτ|λ=iω0≠0.
Proof. Differentiating Eq (3.5) for τ, we have
(dλ(τ)dτ)−1=(3λ2+2D2λ+D1)eλτ+D3λ(D3λ+D4)−τλ. | (3.12) |
Substituting λ(τ)=iω0 into Eq (3.12), we can get
(dλ(τ)dτ)−1|λ=iω0=(−3ω20+2D2iω0+D1)(cosω0τ+isinω0τ)+D3iω0(D3iω0+D4)−τiω0 |
=(D3ω20+iω0D4)(3ω20−D1−2D2ω0i)(cosω0τ+isinω0τ)+D3(D3ω20+iω0D4)D23ω40+D24ω20+τiω0. |
Therefore, we can get
Re(dλ(τ)dτ)−1|λ=iω0=3ω60+(2D22−4D1)ω40+(D21−2D0D2−D23)ω20D23ω40+D24ω20=zkD23ω40+D24ω20[3z2k+(2D22−4D1)zk+(D21−2D0D2−D23)]=zkf′(zk)D23ω40+D24ω20. |
Thus, we have
sign{dλ(τ)dτ|λ=iω0}=sign{Re(dλ(τ)dτ)−1|λ=iω0}=signf′(zk)≠0. |
This ends the proof of Theorem 3.2.
Theorem 3.3. Under the conditions (H1)–(H3), the dynamics of the model (3.1) at the positive equilibrium E∗ can be obtained.
(i) If D7≥0 and the condition (H7) is true, then E∗ is locally asymptotically stable for all τ>0.
(ii) If D7<0 or the condition (H8), D7>0, z∗1>0 and f(z∗1)≤0 are true, then E∗ is locally asymptotically stable when τ∈[0,τ0); but, E∗ is unstable when τ>τ0.
(iii) If the conditions in (ii) are all satisfied and f′(zk)≠0 holds, then a Hopf bifurcation occurs at E∗ when τ=τ0.
For τ=τ0, the existence of Hopf bifurcation has been discussed. The properties of the Hopf bifurcation will be discussed in consideration of the work of Hassard et al. [57].
Let τ=τ0+μ, μ∈R; we can obtain that μ=0 is a Hopf bifurcation value of the model (3.1). Let the phase space be C=C([−1,0],R3) and t→tτ; then, the model (3.1) can be expressed as a functional differential equation (FDE) in C as
˙u(t)=Lμ(ut)+F(μ,ut), | (3.13) |
where Lμ:C→R3 and F:R×C→R3.
Define φ(θ)=(φ1(θ),φ2(θ),φ3(θ))T∈R3,θ∈[−1,0] such that
Lμ(φ)=(τ0+μ)L1φ(0)+(τ0+μ)L2φ(−1) | (3.14) |
and
F(μ,φ)=(τ0+μ)(F1(φ)F2(φ)F3(φ)), | (3.15) |
where L1 and L2 are defined in Eq (3.4), and
F1(φ)=C11φ21(0)+C12φ1(0)φ2(0),F2(φ)=C21φ21(0)+C22φ22(0)+C23φ23(0)+C24φ1(0)φ2(0)+C25φ2(0)φ3(0)+C26φ21(−1)+C27φ1(−1)φ2(0),F3(φ)=C31φ22(0)+C32φ23(0)+C33φ2(0)φ3(0), |
C11=−1+(1−m)2α1β1Z∗[α1+P∗(1−m)]3,C12=−α1β1(1−m)[α1+P∗(1−m)]2,C21=(1−m)2α1β2Z∗[α1+P∗(1−m)]3,C22=γ1aF∗(a+Z∗)3(b+F∗),C23=γ1bZ∗(a+Z∗)(b+F∗)3,C24=α1β2(1−m)[α1+P∗(1−m)]2,C25=−γ1ab(a+Z∗)2(b+F∗)2,C26=δα2Z∗(α2+P∗)2,C27=−δα2(α2+P∗)2,C31=−γ2aF∗(a+Z∗)3(b+F∗),C32=−γ2bZ∗(a+Z∗)(b+F∗)3,C33=γ2ab(a+Z∗)2(b+F∗)2. |
From the Riesz representation theorem[58], there is a matrix function whose elements are functions of bounded variation η(θ,μ)∈C([−1,0],R3); also,
Lμ(φ)=∫0−1dη(θ,μ)φ(θ),φ∈C. | (3.16) |
Let
η(θ,0)=τ0L1δ(θ)−τ0L2δ(θ+1), |
where δ(θ) is a Dirac function, as follows:
δ(θ)={0,θ≠0,1,θ=1. |
For φ∈C=C([−1,0],R3), we define that
A(μ)φ={dφ(θ)dθ,θ∈[−1,0),∫0−1dη(μ,s)φ(s),θ=0,R(μ)φ={0,θ∈[−1,0),F(μ,φ),θ=0. |
Therefore, Eq (3.13) becomes
˙u(t)=Aμ(ut)+Rμut. | (3.17) |
For ψ∈C∗=C([−1,0],R3), let
A∗(μ)ψ(s)={−dψ(s)ds,s∈[0,1),∫0−1ψ(−ξ)dη(ξ,0),s=0, |
and establish the bilinear inner product
⟨ψ(s),φ(θ)⟩=ˉψ(0)φ(0)−∫0−1∫θ0ˉψ(ξ−θ)dη(θ)φ(ξ)dξ,ψ(s)∈C∗,φ(θ)∈C, |
where η(θ)=η(θ,0). Then, A(0) and A∗(0) are adjoint operators. According to Theorem 3.3, we can get that ±iω0τ0 are eigenvalues of A(0). Thus, they are also eigenvalues of A∗(0). Assume that q(θ)=(1,q1,q2)Teiω0τ0θ is the eigenvector of A(0) corresponding to iω0τ0; then, A(0)q(θ)=iω0τ0q(θ). By the definition above, we have
dq(θ)dθ=iωτ0q(θ),θ∈[−1,0), |
L0q(0)=∫0−1dη(0)q(0)=A(0)=iωτ0q(0),θ=0, |
that is,
(τ0+μ)(b11−iω0b120b21+ce−iω0τ0b22−iω0b230b32b33−iω0)(1q1q2)=0. |
Solving it, we can get
(1q1q2)=(1iω0−b11b12b32(iω0−b11)b12(iω0−b33)). |
Similarly, we assume that q∗(s)=D(1,q∗1,q∗2)e−iω0τ0s is the eigenvector of A∗(0) corresponding to −iω0τ0; then, A∗(0)q∗(s)=−iω0τ0q∗(s). By the definition above, we can get
(1q∗1q∗2)=(1−iω0+b11b12+ceiω0τ0b23(iω0+b11)b12(iω0+b33)(b21+ceiω0τ0)). |
By Eq (3.2), ⟨q∗,q⟩ can be expressed as
⟨q∗,q⟩=¯q∗(0)q(0)−∫0−1∫θ0¯q∗(ξ−θ)dη(θ)q(ξ)dξ=ˉD[(1,q∗1,q∗2)(1,q1,q2)T−∫0−1∫θ0(1,¯q∗1,¯q∗2)e−iω0(ξ−θ)τ0dη(θ)(1,q1,q2)Teiξω0τ0dξ]=ˉD[1+q1¯q∗1+q2¯q∗2+τ0c¯q∗1eiω0τ0]. |
Thus, we can choose
ˉD=[1+q1¯q∗1+q2¯q∗2+τ0c¯q∗1eiω0τ0]−1; |
then, ⟨q∗,q⟩=1.
Next, we adopt the ideas of Hassard et al.[57] to compute the coordinates describing the center manifold C0 at μ=0. We assume that ut is the solution of Eq (3.13) when μ=0. Let
z(t)=⟨q∗,ut⟩,W(t,θ)=ut(θ)−2Re{z(t)q(θ)}. | (3.18) |
On the center manifold C0, we can get
W(t,θ)=W(z(t),ˉz(t),θ), |
where
W(z(t),ˉz(t),θ)=W21(θ)z22+W11zˉz+W02(θ)ˉz22+⋯; | (3.19) |
z and ˉz express the local coordinates for the center manifold C0 in the direction of q∗ and ˉq∗. We can get that W is real when ut is real. For the real solution ut∈C0 of Eq (3.17), when μ=0, we have
˙z(t)=⟨q∗,˙u(t)⟩=⟨q∗,A(0)ut+R(0)ut⟩=⟨A∗(0)q∗,ut⟩+¯q∗(0)F(0,ut)=iω0τ0z(t)+¯q∗(0)F(0,W(t,0)+2Re{z(t)q(0)})=iω0τ0z(t)+¯q∗(0)F(0,W(z,ˉz,0)+2Re{z(t)q(0)}def=iω0τ0z(t)+¯q∗(0)F0(z,ˉz). | (3.20) |
Eq (3.20) can also be written as
˙z(t)=iω0τ0z(t)+g(z,ˉz), |
where
g(z,ˉz)=¯q∗(0)F0(z,ˉz)=g20(θ)z22+g11(θ)zˉz+g02(θ)ˉz22+g21z2ˉz2+⋯. | (3.21) |
By Eq (3.18), we can get
˙W=˙ut−˙zq−˙ˉzˉq={AW−2Re{¯q∗(0)F0q(θ)},θ∈[−1,0),AW−2Re{¯q∗(0)F0q(0)}+F0,θ=0,def=AW+H(z,ˉz,θ), | (3.22) |
where
H(z,ˉz,θ)=H20(θ)z22+H11(θ)zˉz+H02(θ)ˉz22+H21z2ˉz2+⋯. | (3.23) |
From Eqs (3.22) and (3.23), we obtain
(A−2iω0τ0)W20(θ)=−H20(θ),AW11(θ)=−H11(θ),(A+2iω0τ0)W02(θ)=−H02(θ). | (3.24) |
By Eq (3.18), we have
ut(θ)=W(t,θ)+2Re{z(t)q(θ)}=W20(θ)z22+W11(θ)zˉz+W02(θ)ˉz22+(1,q1,q2)Teiω0τ0θz+(1,¯q1,¯q2)Te−iω0τ0θˉz+⋯. | (3.25) |
Then, we have
g(z,ˉz)=¯q∗(0)F(0,ut)=ˉD(1,¯q1∗,¯q2∗)τ0(F1(ut)F2(ut)F3(ut))=ˉDτ0[F1(ut)+¯q1∗F2(ut)+¯q2∗F3(ut)]=z22[2ˉDτ0(k11+k21¯q1∗+k31¯q∗2)]+zˉz[ˉDτ0(k12+k22¯q1∗+k32¯q∗2)]+ˉz22[2ˉDτ0(k13+k23¯q1∗+k33¯q∗2)]+z2ˉz2[2ˉDτ0(k14+k24¯q1∗+k34¯q∗2)], | (3.26) |
where
k11=C11+C12q1,k12=2C11+C12(ˉq1+q1),k21=C21+C22q21+C23q22+C24q1+C25q1q2+C26e−2iω0τ0+C27q1e−iω0τ0,k31=C31q21+C32q22+C33q1q2,k33=C31ˉq21+C32ˉq22+C33ˉq1ˉq2,k22=2C21+2C22q1ˉq1+2C23q2ˉq2+C24(q1+ˉq1)+C25(q1ˉq2+ˉq1q2)+2C26+C27(q1eiω0τ0+ˉq1e−iω0τ0),k32=2C31q1ˉq1+2C32q2ˉq2+C33(q1ˉq2+ˉq1q2),k13=C11+C12ˉq1,k23=C21+C22ˉq21+C23ˉq22+C24ˉq1+C25ˉq1ˉq2+C26e2iω0τ0+C27eiω0τ0, |
k14=C11(W(1)20(0)+2W(1)11(0))+C12(12W(1)20(0)ˉq1+q1W(1)11(0)+W(2)11(0)+12W(2)20(0)),k24=C21(W(1)20(0)+2W(1)11(0))+C22(W(2)20(0)+2W(2)11(0))+C23(W(3)20(0)+2W(3)11(0))+C24(12W(1)20(0)ˉq1+q1W(1)11(0)+W(2)11(0)+12W(2)20(0))+C25(12W(2)20(0)ˉq2+W(1)11(0)q2+q1W(3)11(0)+12W(3)20(0)ˉq1)+C26(eiω0τ0W(1)20(−1)+2e−iω0τ0W(1)11(−1))+C27(12ˉq1W(1)20(−1)+q1W(1)11(−1)+e−iω0τ0W(2)11(0)+12eiω0τ0W(2)20(0)),k34=C31(W(2)20(0)+2W(1)11(0))+C33(12ˉq2W(2)11(0)+q2W(2)11(0)+q1W(3)11(0)+12W(3)20(0)ˉq1). |
Comparing the coefficients of Eqs (3.21) and (3.26), we can get
g20=2ˉDτ0(k11+ˉq∗1k21+ˉq∗2k31),g11=ˉDτ0(k12+ˉq∗1k22+ˉq∗2k32),g02=2ˉDτ0(k13+ˉq∗1k23+ˉq∗2k33),g21=2ˉDτ0(k14+ˉq∗1k24+ˉq∗2k34). | (3.27) |
Since the expression of g21 contains W20(θ) and W11(θ), we must compute W20(θ) and W11(θ). According to Eq (3.22), when θ∈[−1,0), we get
H(z,ˉz,θ)=−2Re{ˉq∗(0)F0(z,ˉz)q(θ)}=−2Re{g(z,ˉz)q(θ)}=−g(z,ˉz)q(θ)−ˉg(z,ˉz)ˉq(θ). | (3.28) |
Comparing the coefficients of Eqs (3.23) and (3.28), we can receive
H20(θ)=−g20q(θ)−ˉg02ˉq(θ),H11(θ)=−g11q(θ)−ˉg11ˉq(θ). | (3.29) |
By Eq (3.24), we have
˙W20(θ)=2iω0τ0W20(θ)+g20q(θ)+ˉg02ˉq(θ),˙W11(θ)=g11q(θ)+ˉg11ˉq(θ). | (3.30) |
Solving Eq (3.30), we can obtain
W20(θ)=ig20q(0)ω0τ0eiω0τ0θ+iˉg02ˉq(0)3ω0τ0e−iω0τ0θ+E1e2iω0τ0θ,W11(θ)=−ig11q(0)ω0τ0eiω0τ0θ+iˉg11ˉq(0)ω0τ0e−iω0τ0θ+E2, | (3.31) |
where Ei=(E(1)i,E(2)i,E(3)i)∈R3 (i=1,2) is a constant vector.
For θ=0, from Eq (3.22), we have
H(z,ˉz,0)=−2Re{ˉq∗(0)F0(z,ˉz)q(0)}+F0. |
From Eq (3.23), we can get
H20(0)=−g20q(0)−ˉg02ˉq(0)+2τ0(k11k21k31),H11(0)=−g11q(0)−ˉg11ˉq(0)+τ0(k12k22k32). | (3.32) |
According to the meaning of A(0) and Eq (3.24), we have
∫0−1dη(θ)W20(θ)=2iω0τ0W20−H20(0),∫0−1dη(θ)W11(θ)=−H11(0), | (3.33) |
where η(θ)=η(θ,0).
From Eqs (3.14), (3.16) and (3.33), we can get
τ0L1W20(0)+τ0L2W20(−1)=2iω0τ0−H20(0),τ0L1W11(0)+τ0L2W11(−1)=H11(0). | (3.34) |
Substituting Eqs (3.31) and (3.32) into Eq (3.34), we have
E1=2(2iω0−b11−b120−ce−2iω0τ0−b212iω0−b22−b230−b322iω0−b33)−1(k11k21k31), |
E2=(−b11−b120−c−b21−b22−b230−b32−b33)−1(k11k21k31). |
Thus, all expressions of gij can be represented in full. Also, we have
{c1(0)=i2ω0τ0(g11g20−2|g11|2−|g02|23)+g212,μ1=−Re{c1(0)}Re{λ′(τ0)},μ2=2Re{c1(0)},T1=−Im{c1(0)}+μ1Im{λ′(τ0)}ω0τ0, | (3.35) |
which determine the direction of Hopf bifurcation and the stability of bifurcating periodic solutions on the center manifold at τ=τ0.
Theorem 3.4. From Eq (3.35), we have the following conclusions.
(i) μ1 determines the direction of Hopf bifurcation: if μ1>0, then the Hopf bifurcation is supercritical; if μ1<0, then the Hopf bifurcation is subcritical and the bifurcating periodic solutions exist when τ>τ0;
(ii) μ2 determines the stability of the bifurcating periodic solutions: if μ2<0, then the bifurcating periodic solutions are stable; if μ2>0, then the bifurcating periodic solutions are unstable;
(iii) T1 determines the period of the bifurcating periodic solutions: if T1>0, then the period increases; if T1<0, then the period decreases.
Next, we will analyze the existence and properties of Hopf bifurcation of the model (1.6). The model (1.6) is linearized at the positive equilibrium E∗ in the phase space C=C([−τ,0],R3):
du(t)dt=DΔu(t)+L1u(t)+L2u(t−τ), | (4.1) |
where D=diag{d1,d2,d3} and L1 and L2 are defined in Eq (3.4).
We know that Δ has the eigenvalues −(nl)2 and n∈N0 under Neumann boundary conditions in [0,lπ]. Then, the characteristic equation of the model (4.1) is
λ3+q2nλ2+q1nλ+q0n+e−λτ(q3nλ+q4n)=0, | (4.2) |
where
q2n=(d1+d2+d3)(nl)2−(b11+b22+b33),q1n=(d1d2+d2d3+d1d3)(nl)4−(d1b22+d2b11+d2b33+d3b22+d3b11+d1b33)(nl)2+(b11b22+b22b33+b11b33−b12b21−b23b32),q0n=d1d2d3(nl)6−(d1d2b33+d1d3b22+d2d3b11)(nl)4+(d1b22b33+d2b11b33+d3b11b22)(nl)2+b11b23b32+b12b21b33−b11b22b33,q3n=−cb12,q4n=−d3cb12(nl)2+cb12b33. |
When τ=0, Eq (4.2) can be rewritten as
λ3+q2nλ2+(q1n+q3n)λ+(q0n+q4n)=0. | (4.3) |
Furthermore, we have
q0n+q4n=−d3cb12(nl)2+cb12b33+d1d2d3(nl)6−(d1d2b33+d1d3b22+d2d3b11)(nl)4+(d1b22b33+d2b11b33+d3b11b22)(nl)2+b11b23b32+b12b21b33−b11b22b33, |
q2n(q1n+q3n)−(q0n+q4n)=[(d1+d2+d3)(d1d2+d2d3+d1d3)−d1d2d3](nl)6−[(d1+d2+d3)(d1b22+d2b11+d2b33+d3b22+d3b11+d1b33)+(b11+b22+b33)(d1d2+d2d3+d1d3)−(d1d2b33+d1d3b22+d2d3b11)](nl)4+[(d1+d2+d3)(b11b22+b22b33+b11b33−b12b21−b23b32−cb12)+(b11+b22+b33)(d1b22+d2b11+d2b33+d3b22+d3b11+d1b33)−(d1b22b33+d2b11b33+d3b11b22)+d3cb12](nl)2+(b11+b22+b33)(cb12−b11b22−b22b33−b11b33+b12b21+b23b32)−(b11b23b32+b12b21b33−b11b22b33)−cb12b33. |
Assume that the following condition holds true: (H9): q2n>0, q2n(q1n+q3n)−(q0n+q4n)>0, q0n+q4n>0, n∈N0. Using the Routh-Hurwitz criterion[54], we can obtain the next conclusion.
Theorem 4.1. If the conditions (H1)–(H3) and (H9) hold, then all roots of Eq (4.3) have negative real parts, that is, the positive equilibrium E∗ of the model (1.6) is locally asymptotically stable when τ=0.
When τ≠0, the time delay may have some effect on the model (1.6). Therefore, we will analyze the effect of delay τ on the positive equilibrium E∗. Let λ=iω2(ω2>0) be the solution of Eq (4.2); we can get
−iω32−ω22q2n+iω2q1n+q0n+e−iω2τ(iω2q3n+q4n)=0. | (4.4) |
Separating the real and imaginary parts of Eq (4.4), we can receive
{ω32−q1nω2=q3nω2cosω2τ−q4nsinω2τ,q2nω22−q0n=q3nω3nsinω2τ+q4ncosω2τ, | (4.5) |
which follows that
ω62+(q22n−2q1n)ω42+(q21n−2q0nq2n−q23n)ω22+(q20n−q24n)=0. | (4.6) |
Let p2n=q22n−2q1n, p1n=q21n−2q0nq2n−q23n, p0n=q20n−q24n and ω22=y; we can obtain
y3+p2ny2+p1ny+p0n=0. | (4.7) |
Further, if f(y)=y3+p2ny2+p1ny+p0n, then f′(y)=3y2+2p2ny+p1n. Assume that (H10):Δ2=p22n−3p1n≤0 and (H11):Δ2=p22n−3p1n>0 are true, here y∗1=−p1n+√Δ23 and y∗2=−p1n−√Δ23 are the local minimum and the local maximum of Eq (4.7), respectively. Similarly, we have the following conclusion by Lemmas 2.2 and 4.2 in[56].
Theorem 4.2. For Eq (4.7), the following results are true.
(1) If p0n<0, then Eq (4.7) has at least one positive root.
(2) If p0n≥0 and the condition (H10) holds, then Eq (4.7) has no positive root.
(3) If p0n≥0 and the condition (H11) holds, then Eq (4.7) has positive roots when y∗1>0 and f(y∗1)≤0.
Without loss of generality, we suppose that it has three positive roots defined by y1n, y2 and y3n. So, Eq (4.6) has three positive roots:
ω12=y1n,ω22=y2n,ω32=y3n. |
Substituting ωk2(k=1,2,3) into Eq (4.5), we can get
τjkn=1ωk2arccos[q3n(ωk2)4+(q2nq4n−q1nq3n)(ωk2)2−q0nq4nq23n(ωk2)2+q24n]+2jπωk2,k=1,2,3,j=0,1,2,⋯. | (4.8) |
Thus, when τ=τjkn, we get that λ=±iωk2 is a pair of purely imaginary roots of Eq (4.2). Define
\begin{equation*} \tau_{n0} = \min\limits_{k = 1,2,3.\\ j\in N_{0}}\{\tau_{kn}^{j}\},\quad \omega_{n0} = \omega_{2}^k|_{\tau = \tau_{kn^0}}. \end{equation*} |
Assume that \lambda(\tau) = \varepsilon_{1}(\tau)+i\varepsilon_{2}(\tau) is the root of Eq (4.2) near \tau = \tau_{kn}^j , and that it {satisfies} \varepsilon_{1}(\tau_{kn}^j) = 0 and \varepsilon_{2}(\tau_{kn}^j) = \omega_{2}^k , k = 1, 2, 3, j\in N_{0}, n\in N_{0} .
Theorem 4.3. If y_{kn} = (\omega_{2}^k)^2 and f'(y_{kn})\neq0, then we have that \frac{{\rm d}Re\lambda(\tau)}{{\rm d}\tau}|_{\lambda = i\omega_{n0}}\neq0.
Proof. Differentiating Eq (4.2) for \tau , we can get
\begin{equation} (\frac{{\rm d}\lambda(\tau)}{{\rm d}\tau})^{-1} = \frac{(3\lambda^2+2q_{2n}\lambda+q_{1n})e^{\lambda\tau} +q_{3n}}{\lambda(q_{3n}\lambda+q_{4n})}-\frac{\tau}{\lambda}. \end{equation} | (4.9) |
Substituting \lambda = i\omega_{n0} into Eq (4.9), we have
(\frac{{\rm d}\lambda(\tau)}{{\rm d}\tau})^{-1}|_{\lambda = i\omega_{n0}} = \frac{(-3\omega_{n0}^2+2q_{2n}i\omega_{n0}+q_{1n})(cos\omega_{n0}\tau+isin\omega_{n0}\tau) +q_{3n}}{-q_{3n}\omega_{n0}^2+i\omega_{n0}q_{4n}}+\frac{\tau}{i\omega_{n0}}; |
then,
Re(\frac{{\rm d}\lambda(\tau)}{{\rm d}\tau})^{-1}|_{\lambda = i\omega_{n0}} = \frac{x_{kn}y'(x_{kn})}{[q_{3n}(\omega_{n)}^2]^2+(\omega_{n0}q_{4n})^2}\neq0. |
Thus, we can get
\rm{sign}\{\frac{{\rm d}Re\lambda(\tau)}{{\rm d}\tau}|_{\lambda = i\omega_{n0}}\} = \rm{sign}\{Re(\frac{{\rm d}\lambda(\tau)}{{\rm d}\tau})^{-1}|_{\lambda = i\omega_{n0}}\}\neq0. |
The proof of Theorem 4.3 is completed.
By computing as shown above, we have the following conclusion.
Theorem 4.4. Under the conditions (H_{1}) – (H_{3}) , for the the positive equilibrium E_{*} of the model (1.6), we have the following:
(i) if the condition (H_{10}) is true and p_{0n}\geq0 , then E_{*} is locally asymptotically stable for all \tau > 0;
(ii) if p_{0n} < 0 or p_{0n}\geq0 , y_{1}^{*} > 0 , f(y_{1}^*)\leq0 and the condition (H_{11}) holds, then E_{*} is locally asymptotically stable when \tau\in[0, \tau_{n0}) ; but, E_{*} is unstable when \tau > \tau_{n0};
(iii) if the conditions in (ii) are all satisfied and f'(y_{kn})\neq0 , then the spatially homogeneous Hopf bifurcation occurs at E_{*} when \tau = \tau_{0} and n = 0 ; and, the spatially inhomogeneous Hopf bifurcation occurs at E_{*} when \tau = \tau_{n0} and n > 0 .
Let \tau_{n} = \tau_{kn}^j, \omega_{n} = \omega_{2n}^k and \tau = \tau_{n}+\mu_{n}, \mu_{n}\in R . Thus, \mu_{n} = 0 is a Hopf bifurcation value of the model (1.6). Let t\rightarrow \frac{t}{\tau} ; then, the model (1.6) can be expressed as an FDE in \mathcal{C} = C([-1, 0], R^3) , as follows:
\begin{equation} \dot{u}(t) = \tau_{n}D\Delta u(t)+L(\tau_{n})u_{t}+F_{n}(\mu_{n},u_{t}), \end{equation} | (4.10) |
where L(\theta):\mathcal{C}\rightarrow X, F_{n}(\mu_{n}, u_{t}):\mathcal{C}\rightarrow X satisfies
\begin{equation*} L(\theta)(\varphi) = \theta\left({\begin{array}{ccccc} b_{11}\varphi_{1}(0)+b_{12}\varphi_{2}(0)\\ b_{21}\varphi_{1}(0)+b_{22}\varphi_{2}(0)+b_{23}\varphi_{3}(0)+c\varphi_{1}(-1)\\ b_{32}\varphi_{2}(0)+b_{33}\varphi_{3}(0)\\ \end{array}}\right) \end{equation*} |
and
\begin{equation} F_{n}(\mu_{n},\varphi) = \mu_{n}D\Delta\varphi(0)+L(\mu_{n})\varphi+F(\mu_{n},\varphi), \end{equation} | (4.11) |
where F(\mu_{n}, \varphi) is defined in Eq (3.15), L_{1} and L_{2} are defined in Eq (3.4) and \varphi = (\varphi_{1}, \varphi_{2}, \varphi_{3})^T\in \mathcal{C}.
The linear equation of Eq (4.10) at O(0, 0, 0) is
\begin{equation} \dot{u}(t) = \tau_{n}D\Delta u(t)+L(\tau_{n})u_{t}. \end{equation} | (4.12) |
Let \Lambda = \{i\omega_{n}\tau_{n}, -i\omega_{n}\tau_{n}\} and z_{t}(\theta)\in C = C([-1, 0], R^3) ; we consider the following FDE:
\begin{equation} \dot{z}(t) = L(\tau_{n})(z_{t}). \end{equation} | (4.13) |
On the basis of the Riesz representation theorem, there exists a 3\times 3 matrix function \eta_{n}(\theta, \mu) (-1\leq\theta\leq0)\in C([-1, 0], R^3) , and it satisfies
L(\tau_{n})(\varphi) = \int_{-1}^0d\eta_{n}(\theta,\mu_{n})\varphi(\theta),\quad \varphi\in C([-1,0],R^3). |
Let
\eta_{n}(\theta,\mu_{n}) = (\tau_{n}+\mu_{n})L_{1}\delta(\theta)-(\tau_{n}+\mu_{n})L_{2}\delta(\theta+1), |
where \delta(\theta) is the Dirac delta function.
We set C^* = C([0, 1], R^{3*}), and R^{3*} is the three - dimensional vector space of row vectors. The bilinear inner product is
\begin{equation} (\psi(s),\varphi(\theta)) = \psi(0)\varphi(0)-\int_{-1}^{0}\int_{0}^{\theta}\psi(\xi-\theta)d\eta_{n}(\theta)\varphi(\xi)d\xi,\quad \psi(s)\in{C^*},\quad\varphi(\theta)\in{C}. \end{equation} | (4.14) |
A_{n}(\tau_{n}) describes the infinitesimal generator of the semigroup induced by the solutions of Eq (4.13), and A_{n}^*(\tau_{n}) denotes the formal adjoint generator of A_{n}(\tau_{n}) satisfying Eq (4.14). Let V and V^* denote the center spaces of the generators A_{n}(\tau_{n}) and A_{n}^*(\tau_{n}) corresponding to \Lambda , respectively. Therefore, V^* is the adjoint space of V , {\rm dim}V = {\rm dim}V^* .
Lemma 4.1. Let
\begin{equation*} \begin{split} V_{1}& = \frac{i\omega_{n}-b_{11}}{b_{12}},\quad \quad\quad V_{2} = \frac{b_{32}(i\omega_{n}-b_{11})}{b_{12}(i\omega_{n}-b_{33})},\\ V_{1}^*& = \frac{i\omega_{n}-b_{11}}{b_{21}+ce^{-i\omega_{n}\tau_{n}}},\quad V_{2}^* = \frac{b_{23}(i\omega_{n}-b_{11})}{(i\omega_{n}-b_{33})(b_{21}+ce^{-i\omega_{n}\tau_{n}})}; \end{split} \end{equation*} |
then, p_{1}(\theta)=(1, V_{1}, V_{2})^Te^{i\omega_{n}\tau_{n}\theta} and p_{2}(\theta)=\overline{p_{1}(\theta)} , -1\leq\theta\leq0 form the basis of V associated with \Lambda ; p_{1}^*(s)=(1, V_{1}^*, V_{2}^*)e^{-i\omega_{n}\tau_{n}s} and p_{2}^*(s)=\overline{p_{1}^*(s)} , 0\leq s\leq1 form the basis of V^* associated with \Lambda .
Denote \Phi = (\Phi_{1}, \Phi_{2}) and \Psi^* = (\Psi_{1}^*, \Psi_{2}^*)^T , where
\begin{equation*} \begin{split} \Phi_{1}(\theta)& = \frac{p_{1}(\theta)+p_{2}(\theta)}{2} = (Re\{e^{i\omega_{n}\tau_{n}\theta}\},Re\{V_{1}e^{i\omega_{n}\tau_{n}\theta}\},Re\{V_{2}e^{i\omega_{n}\tau_{n}\theta}\})^T,\quad\theta\in[-1,0],\\ \Phi_{2}(\theta)& = \frac{p_{1}(\theta)-p_{2}(\theta)}{2i} = (Im\{e^{i\omega_{n}\tau_{n}\theta}\},Im\{V_{1}e^{i\omega_{n}\tau_{n}\theta}\},Im\{V_{2}e^{i\omega_{n}\tau_{n}\theta}\})^T,\quad\theta\in[-1,0],\\ \Psi_{1}(s)& = \frac{p_{1}^*(s)+p_{2}^*(s)}{2} = (Re\{e^{-i\omega_{n}\tau_{n}s}\},Re\{V_{1}^*e^{-i\omega_{n}\tau_{n}s}\},Re\{V_{2}^*e^{-i\omega_{n}\tau_{n}s}\}),\quad s\in[0,1],\\ \Psi_{2}(s)& = \frac{p_{1}^*(s)-p_{2}^*(s)}{2i} = (Im\{e^{-i\omega_{n}\tau_{n}s}\},Im\{V_{1}^*e^{-i\omega_{n}\tau_{n}s}\},Im\{V_{2}^*e^{-i\omega_{n}\tau_{n}s}\}),\quad s\in[0,1]. \end{split} \end{equation*} |
Suppose that (\Psi^*, \Phi) = (\Psi_{i}^*, \Phi_{j})(i, j = 1, 2.) is the basis \Psi of V^* , which satisfies
\Psi = (\Psi_{1},\Psi_{2})^T = (\Psi^*,\Phi)^{-1}\Psi^*. |
Thus, we have that (\Psi, \Phi) = I_{2\times2} .
Let f_{n} = (\xi_{n}^1, \xi_{n}^2, \xi_{n}^3), where \xi_{n}^1 = (cos\frac{n}{l}x, 0, 0)^T, \xi_{n}^2 = (0, cos\frac{n}{l}x, 0)^T and \xi_{n}^3 = (0, 0, cos\frac{n}{l}x)^T . \xi_{n}^{j} (j = 1, 2, 3) denotes the eigenfunctions on R^3 of the eigenvalues -(\frac{n}{l})^2 , n = 0, 1, 2\cdots . Define c_{n}\cdot f_{n} = c_{1}\xi_{n}^1+c_{2}\xi_{n}^2+c_{3}\xi_{n}^3 , c_{n} = (c_{1}, c_{2}, c_{3})^T, c_{j}\in R, j = 1, 2, 3, and the center space of Eq (4.12) is written as
P_{CN}\varphi = \Phi(\Psi,\langle\varphi,f_{n}\rangle)\cdot f_{n}, |
where \varphi\in C , C = P_{CN}C\oplus P_{s}C and P_{s}C expresses the complementary subspace of P_{CN}C .
According to [57] and [59], the center space of the linear model of (4.10) with \mu_{n} = 0 is expressed as P_{CN}C , where
\begin{equation*} P_{CN}C = \{\frac{1}{2}(p_{1}(\theta)z+p_{2}(\theta)\bar{z})\cdot f_{n},\quad z\in C\}. \label{4.16} \end{equation*} |
Thus, the solution of the model (4.10) can be written as
\begin{equation*} u_{t} = \frac{1}{2}(p_{1}(\theta)z+p_{2}(\theta)\bar{z})\cdot f_{n}+Q(z(t),\bar{z}(t))(\theta), \label{4.17} \end{equation*} |
where Q(z(t), \bar{z}(t))(\theta) = W(\frac{z+\bar{z}}{2}, i\frac{z-\bar{z}}{2}, 0) , z = x_{1}-ix_{2}.
By Wu[59], z satisfies
\begin{equation} \dot{z} = i\omega_{n}\tau_{n}z+g_{n}(z,\bar{z}), \end{equation} | (4.15) |
where
\begin{equation*} g_{n}(z,\bar{z}) = (\Psi_{1}(0)-i\Psi_{2}(0))\langle F_{n}(0,u_{t}),f{n}\rangle,\Psi(0) = (\Psi_{1}(0),\Psi_{2}(0))^T. \label{4.19} \end{equation*} |
Let
\begin{equation} Q(z,\bar{z}) = Q_{20}(\theta)\frac{z^2}{2}+Q_{11}(\theta)z\bar{z} +Q_{02}(\theta)\frac{\bar{z}^2}{2}+\cdots, \end{equation} | (4.16) |
\begin{equation*} g_{n}(z,\bar{z}) = \tilde{g}_{20}(\theta)\frac{z^2}{2}+\tilde{g}_{11}(\theta)z\bar{z} +\tilde{g}_{02}(\theta)\frac{\bar{z}^2}{2}+\tilde{g}_{21}\frac{z^2\bar{z}}{2}+\cdots \label{4.21} \end{equation*} |
and
\begin{equation*} \Psi_{1}(0)-i\Psi_{2}(0) = (\psi_{1},\psi_{2},\psi_{3}). \label{4.22} \end{equation*} |
By computing and comparing the coefficients, we have
\begin{equation*} \begin{split} \tilde{g}_{20} = &\frac{\tau_{n}}{2}\langle[(C_{11}+C_{12}V_{1})\psi_{1}+(C_{21}+C_{22}V_{1}^2+C_{23}V_{2}^2 +C_{24}V_{1}+C_{25}V_{1}V_{2}+C_{26}e^{-2i\omega_{n}\tau_{n}}\\ &+C_{27}V_{1}e^{-i\omega_{n}\tau_{n}})\psi_{2} +(C_{31}V_{1}^2+C_{32}V_{2}^2+C_{33}V_{1}V_{2})\psi_{3}]cos^2{(\frac{n}{l})x},cos(\frac{n}{l})x\rangle,\\ \end{split} \label{4.23} \end{equation*} |
\begin{equation*} \begin{split} \tilde{g}_{11} = &\frac{\tau_{n}}{4}\langle[(2C_{11}+C_{12}(V_{1}+\bar{V}_{1}))\psi_{1}+(2C_{21} +2C_{22}V_{1}\bar{V}_{1}+2C_{23}V_{2}\bar{V}_{2}+C_{24}(V_{1}+\bar{V}_{1})\\ &+C_{25}(V_{1}\bar{V}_{2}+\bar{V}_{1}V_{2})+2C_{26}+C_{27}(V_{1}e^{i\omega_{n}\tau_{n}}+\bar{V}_{1}e^{-i\omega_{n}\tau_{n}}))\psi_{2}\\ &+(2C_{31}V_{1}\bar{V}_{1}+2C_{32}V_{2}\bar{V}_{2}+C_{33}(V_{1}\bar{V}_{2}+\bar{V}_{1}V_{2}))\psi_{3}]cos^2{(\frac{n}{l})x},cos(\frac{n}{l})x\rangle,\\ \tilde{g}_{21} = &\tau_{n}\{\langle[C_{11}(Q_{20}^{(1)}(0)+2Q_{11}^{(1)}(0))+C_{12}(\frac{1}{2}Q_{20}^{(1)}(0)\bar{V}_{1} +V_{1}Q_{11}^{(1)}(0)+Q_{11}^{(2)}(0)\\ &+\frac{1}{2}Q_{20}^{(2)}(0))]cos(\frac{n}{l})x,cos(\frac{n}{l})x\rangle\psi_{1}+\langle[C_{21}(Q_{20}^{(1)}(0) +2Q_{11}^{(1)}(0))+C_{22}(Q_{20}^{(2)}(0)\\ &+2Q_{11}^{(2)}(0))+C_{23}(Q_{20}^{(3)}(0)+2Q_{11}^{(3)}(0))+C_{24}(\frac{1}{2}Q_{20}^{(1)}(0)\bar{V}_{1} +V_{1}Q_{11}^{(1)}(0)+Q_{11}^{(2)}(0)\\ &+\frac{1}{2}Q_{20}^{(2)}(0))+C_{25}(\frac{1}{2}Q_{20}^{(2)}(0)\bar{V}_{2}+Q_{11}^{(1)}(0)V_{2}+V_{1}Q_{11}^{(3)}(0)+\frac{1}{2}Q_{20}^{(3)}(0)\bar{V}_{1})\\ &+C_{26}(e^{i\omega_{n}\tau_{n}}Q_{20}^{(1)}(-1)+2e^{-i\omega_{n}\tau_{n}}Q_{11}^{(1)}(-1)) +C_{27}(\frac{1}{2}\bar{V}_{1}Q_{20}^{(1)}(-1)+V_{1}Q_{11}^{(1)}(-1)\\ &+e^{-i\omega_{n}\tau_{n}}Q_{11}^{(2)}(0)+\frac{1}{2}e^{i\omega_{n}\tau_{n}}Q_{20}^{(2)}(0))]cos(\frac{n}{l})x,cos(\frac{n}{l})x\rangle\psi_{2} +\langle[C_{31}(Q_{20}^{(2)}(0)+2Q_{11}^{(1)}(0))\\ &+C_{33}(\frac{1}{2}\bar{V}_{2}Q_{11}^{(2)}(0)+V_{2}Q_{11}^{(2)}(0)+V_{1}Q_{11}^{(3)}(0)+\frac{1}{2}Q_{20}^{(3)}(0)\bar{V}_{1})]cos(\frac{n}{l})x,cos(\frac{n}{l})x\rangle\psi_{3}\}. \end{split} \label{4.24} \end{equation*} |
We know that \int_{0}^{\pi}cos^3{(\frac{n}{l})x}dx = 0 and \tilde{g}_{02} = \overline{\tilde{g}}_{20} . Therefore, we can get that \tilde{g}_{20} = \tilde{g}_{11} = \tilde{g}_{02} = 0 when n = 1, 2, 3\cdots . When n = 0 , we have
\begin{equation} \begin{split} \tilde{g}_{20} = &\frac{\tau_{n}}{2}[(C_{11}+C_{12}V_{1})\psi_{1}+(C_{21}+C_{22}V_{1}^2+C_{23}V_{2}^2 +C_{24}V_{1}+C_{25}p_{1}V_{2}+C_{26}e^{-2i\omega_{n}\tau_{n}}\\ &+C_{27}V_{1}e^{-i\omega_{n}\tau_{n}})\psi_{2} +(C_{31}V_{1}^2+C_{32}V_{2}^2+C_{33}V_{1}V_{2})\psi_{3}],\\ \tilde{g}_{11} = &\frac{\tau_{n}}{4}[(2C_{11}+C_{12}(V_{1}+\bar{V}_{1}))\psi_{1}+(2C_{21} +2C_{22}V_{1}\bar{V}_{1}+2C_{23}V_{2}\bar{V}_{2}+C_{24}(V_{1}+\bar{V}_{1})\\ &+C_{25}(V_{1}\bar{V}_{2}+\bar{V}_{1}V_{2})+2C_{26}+C_{27}(V_{1}e^{i\omega_{n}\tau_{n}}+\bar{V}_{1}e^{-i\omega_{n}\tau_{n}}))\psi_{2}\\ &+(2C_{31}V_{1}\bar{V}_{1}+2C_{32}V_{2}\bar{V}_{2}+C_{33}(V_{1}\bar{V}_{2}+\bar{V}_{1}V_{2}))\psi_{3}]. \end{split} \end{equation} | (4.17) |
Considering the expression of \tilde{g}_{21} , it contains Q_{20}(\theta) and Q_{11}(\theta) , so we must compute Q_{20}(\theta) and Q_{11}(\theta) . Seeking the derivative on both sides of Eq (4.16), we have
\begin{equation} \dot{Q}(z,\bar{z}) = Q_{20}z+Q_{11}\dot{z}\bar{z}+Q_{11}z\dot{\bar{z}} +Q_{02}\dot{\bar{z}}\bar{z}+\cdots, \end{equation} | (4.18) |
\begin{equation} A_{\tau_{n}}Q(z,\bar{z}) = A_{\tau_{n}}Q_{20}\frac{z^2}{2}+A_{\tau_{n}}Q_{11}z\bar{z} +A_{\tau_{n}}Q_{02}\frac{\bar{z}^2}{2}+\cdots. \end{equation} | (4.19) |
From Wu [59], Q(z, \bar{z}) satisfies
\begin{equation} \dot{Q}(z,\bar{z}) = A_{\tau_{n}}Q(z,\bar{z})+S(z,\bar{z}), \end{equation} | (4.20) |
where
\begin{equation*} S(z,\bar{z}) = S_{20}\frac{z^2}{2}+S_{11}z\bar{z}+S_{02}\frac{\bar{z}^2}{2}+\cdots = X_{0}F_{n}(u_{t},0)-\Phi(\Psi,\langle X_{0}F_{n}(u_{t},0),f_{n}\rangle)\cdot f_{n}, \end{equation*} |
and
\begin{equation*} X_{0}(\theta) = \begin{cases} I,\quad &\theta = 0,\\ 0,\quad &-1\leq\theta < 0, \end{cases} \end{equation*} |
S_{ij}\in P_{S}C,\quad i+j = 2. |
From Eq (4.15) and Eqs (4.17)–(4.20), we can get
\begin{equation*} \left\{\begin{aligned} (2i\omega_{n}\tau_{n}-A_{\tau_{n}})Q_{20}& = S_{20},\\ -A_{\tau_{n}}Q_{11}& = S_{11}. \end{aligned} \right. \label{4.30} \end{equation*} |
Because A_{\tau_{n}} has only two characteristic roots with a zero real part, i.e., \pm i\omega_{n}\tau_{n} , Eq (4.20) has a unique solution Q_{ij}(i+j = 2) in P_{S}C , which satisfies
\begin{equation} \left\{ \begin{aligned} Q_{20}& = (2i\omega_{n}\tau_{n}-A_{\tau_{n}})^{-1}S_{20},\\ Q_{11}& = -A_{\tau_{n}}^{-1}S_{11}. \end{aligned} \right. \end{equation} | (4.21) |
From Eq (4.20), we can get that, for \theta\in[-1, 0) ,
\begin{equation} \begin{split} S(z,\bar{z})& = -\Phi(\theta)\Psi(0)\langle F_{n}(u_{t},0),f_{n}\rangle\cdot f_{n}\\ & = -\frac{\tau_{n}}{2}[(p_{1}(\theta)\tilde{g}_{20}+\overline{\tilde{g}_{02}}p_{2}(\theta))\cdot f_{n}\cdot\frac{z^2}{2}+(p_{1}(\theta)\tilde{g}_{11}+\overline{\tilde{g}_{11}}p_{2}(\theta))\cdot f_{n}\cdot z\bar{z}\\ &+(p_{1}(\theta)\tilde{g}_{02}+\overline{\tilde{g}_{20}}p_{2}(\theta))\cdot f_{n}\cdot\frac{\bar{z}^2}{2}]+\cdots. \end{split} \end{equation} | (4.22) |
Comparing the coefficients in Eqs (4.20) and (4.22), when \theta\in[-1, 0) , we can receive
\begin{equation} \begin{split} S_{20}(\theta)& = -\frac{\tau_{n}}{2}(p_{1}(\theta)\tilde{g}_{20}+\overline{\tilde{g}_{02}}p_{2}(\theta))cos(\frac{n}{l})x,\\ S_{11}(\theta)& = -\frac{\tau_{n}}{2}(p_{1}(\theta)\tilde{g}_{11}+\overline{\tilde{g}_{11}}p_{2}(\theta))cos(\frac{n}{l})x. \end{split} \end{equation} | (4.23) |
When \theta = 0 , we have
\begin{equation} \begin{split} S_{20}(0) = &\frac{\tau_{n}}{2} \left( \begin{array}{c} C_{11}+C_{12}p_{1}\\ C_{21}+C_{22}p_{1}^2+C_{23}p_{2}^2+C_{24}p_{1}+C_{25}p_{1}p_{2}+C_{26}e^{-2i\omega_{n}\tau_{n}}\\ +C_{27}p_{1}e^{-i\omega_{n}\tau_{n}}\\ C_{31}p_{1}^2+C_{32}p_{2}^2+C_{33}p_{1}p_{2} \end{array} \right)cos^2{(\frac{n}{l})x}\\ &-\frac{\tau_{n}}{2}(p_{1}(\theta)\tilde{g}_{20}+\overline{\tilde{g}_{02}}p_{2}(\theta))cos(\frac{n}{l})x, \end{split} \end{equation} | (4.24) |
\begin{equation} \begin{split} S_{11}(0) = &\frac{\tau_{n}}{4} \left( \begin{array}{c} 2C_{11}+C_{12}(p_{1}+\bar{p}_{1})\\ 2C_{21}+2C_{22}p_{1}\bar{p}_{1}+2C_{23}p_{2}\bar{p}_{2}+C_{24}(p_{1}+\bar{p}_{1})+C_{25}(p_{1}\bar{p}_{2}+\bar{p}_{1}p_{2})\\ +2C_{26}+C_{27}(p_{1}e^{i\omega_{n}\tau_{n}}+\bar{p}_{1}e^{-i\omega_{n}\tau_{n}})\\ 2C_{31}p_{1}\bar{p}_{1}+2C_{32}p_{2}\bar{p}_{2}+C_{33}(p_{1}\bar{p}_{2}+\bar{p}_{1}p_{2}) \end{array} \right)cos^2{(\frac{n}{l})x}\\ &-\frac{\tau_{n}}{2}(p_{1}(\theta)\tilde{g}_{11}+\overline{\tilde{g}_{11}}p_{2}(\theta))cos(\frac{n}{l})x. \end{split} \end{equation} | (4.25) |
Using Eqs (4.24) and (4.25), we can get Q_{20}(0), Q_{11}(0), Q_{20}(-1) and Q_{11}(-1) . Because p_{1}(\theta) = p_{1}(0)e^{i\omega_{n}\tau_{n}\theta}, \theta\in[-1, 0), from Eqs (4.21)–(4.25), we have
\begin{equation*} \begin{split} Q_{20}(\theta) = &\frac{i}{2}\Big[\frac{\tilde{g}_{20}}{\omega_{n}\tau_{n}}p_{1}(0)e^{i\omega_{n}\tau_{n}\theta} +\frac{\overline{\tilde{g}_{02}}}{3\omega_{n}\tau_{n}}p_{1}(0)e^{-i\omega_{n}\tau_{n}\theta}\Big]+e^{2i\omega_{n}\tau_{n}\theta}E_{3},\\ Q_{11}(\theta) = &\frac{i}{2}\Big[\frac{\tilde{g}_{11}}{\omega_{n}\tau_{n}}p_{1}(0)e^{i\omega_{n}\tau_{n}\theta} -\frac{\overline{\tilde{g}_{11}}}{\omega_{n}\tau_{n}}p_{1}(0)e^{-i\omega_{n}\tau_{n}\theta}\Big]+E_{4}, \end{split} \label{4.36} \end{equation*} |
where E_{3} = (E_{3}^{(1)}, E_{3}^{(2)}, E_{3}^{(3)})\in R^3 and E_{4} = (E_{4}^{(1)}, E_{4}^{(2)}, E_{4}^{(3)})\in R^3 satisfy
\begin{equation*} E_{3} = \frac{1}{2} \left( \begin{array}{ccc} 2i\omega_{n}-b_{11} & -b_{12}& 0\\ -ce^{-2i\omega_{n}\tau_{n}}-b_{21} & 2i\omega_{n}-b_{22} & -b_{23}\\ 0& -b_{32}& 2i\omega_{n}-b_{33} \end{array} \right)^{-1} \end{equation*} |
\begin{equation*} \times \left( \begin{array}{c} C_{11}+C_{12}V_{1}\\ C_{21}+C_{22}V_{1}^2+C_{23}V_{2}^2+C_{24}V_{1}+C_{25}V_{1}V_{2}+C_{26}e^{-2i\omega_{n}\tau_{n}}\\ +C_{27}V_{1}e^{-i\omega_{n}\tau_{n}}\\ C_{31}V_{1}^2+C_{32}V_{2}^2+C_{33}V_{1}V_{2} \end{array} \right)cos^2{(\frac{n}{l})x}, \end{equation*} |
\begin{equation*} E_{4} = \left( \begin{array}{ccc} -b_{11} & -b_{12}& 0\\ -c-b_{21} & -b_{22} & -b_{23}\\ 0& -b_{32}& -b_{33} \end{array} \right)^{-1} \end{equation*} |
\begin{equation*} \times \left( \begin{array}{c} 2C_{11}+C_{12}(V_{1}+\bar{V}_{1})\\ 2C_{21}+2C_{22}V_{1}\bar{V}_{1}+2C_{23}V_{2}\bar{V}_{2}+C_{24}(V_{1}+\bar{V}_{1})+C_{25}(V_{1}\bar{V}_{2}+\bar{V}_{1}V_{2})\\ +2C_{26}+C_{27}(V_{1}e^{i\omega_{n}\tau_{n}}+\bar{V}_{1}e^{-i\omega_{n}\tau_{n}})\\ 2C_{31}V_{1}\bar{V}_{1}+2C_{32}V_{2}\bar{V}_{2}+C_{33}(V_{1}\bar{V}_{2}+\bar{V}_{1}V_{2}) \end{array} \right)cos^2{(\frac{n}{l})x}. \end{equation*} |
Therefore, we can get the following values:
\begin{equation} \left\{\begin{array}{lr} c_{2}(0) = \frac{i}{2\omega_{n}\tau_{n}}(\tilde{g}_{11}\tilde{g}_{20}-2|\tilde{g}_{11}|^2-\frac{|\tilde{g}_{02}|^2}{3})+\frac{\tilde{g}_{21}}{2},\\ \mu_{3} = -\frac{Re\{c_{2}(0)\}}{Re\{\lambda'(\tau_{n})\}},\\ \mu_{4} = 2Re\{c_{2}(0)\},\\ T_{2} = -\frac{Im\{c_{2}(0)\}+\mu_{3}Im\{\lambda'(\tau_{n})\}}{\omega_{n}\tau_{n}}, \end{array} \right. \end{equation} | (4.26) |
which determine the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions on the center manifold at \tau = \tau_{n} .
Theorem 4.5. According to Eq (4.26), we have the following conclusions.
(i) \mu_{3} determines the direction of Hopf bifurcation: if \mu_{3} > 0 , then the Hopf bifurcation is supercritical; if \mu_{3} < 0 , then the Hopf bifurcation is subcritical and the bifurcating periodic solutions exist when \tau > \tau_{n};
(ii) \mu_{4} determines the stability of the bifurcating periodic solutions: if \mu_{4} > 0 , then the bifurcating periodic solutions are unstable; if \mu_{4} < 0 , then the bifurcating periodic solutions are stable;
(iii) T_{2} determines the period of the bifurcating periodic solutions: if T_{2} > 0 , then the bifurcating periodic solutions increase; if T_{2} < 0 , then the bifurcating periodic solutions decrease.
With the help of Matlab software, the stability of the positive equilibrium E_{*}(P_{*}, Z_{*}, F_{*}) was simulated with the given values of all parameters in order to confirm the previous theoretical results.
First, we assume that P(0) = 0.8, Z(0) = 40 and F(0) = 0.1 for the model (2.1). And, we take the values of all other parameters as follows: \beta_{1} = 0.016, \beta_{2} = 0.7, \gamma_{1} = 0.0875, \gamma_{2} = 0.075, \alpha_{1} = 0.1, \alpha_{2} = 0.2, a = 0.5, b = 0.25, g_{1} = 0.1, g_{2} = 0.2, m = 0.8 and \delta = 0.35 . According to Theorem 2.1, we can know that the model (2.1) has one unique positive equilibrium E_{*}(0.5043, 31.1137, 0.1191) . From Theorem 2.5, E_{*} is locally asymptotically stable (see Figure 3).
In the second section, we obtained Theorem 2.6 by referring to [55]. Now, we will verify some conclusions by taking m as the bifurcation parameter. When \delta = 0.2 , the critical value is m^* = 0.82 , which satisfies M_{6}(m^*) > 0 , T(m^*) = 0 , M_{8}(m^*) > 0 and \frac{{\rm d}T}{{\rm d}s}|_{m = m^*}\neq0 . That is, all conditions in Theorem 2.6 are satisfied. Therefore, Hopf bifurcation occurs at E_{*} when m = m^* . Meanwhile, we can obtain that P(t) reaches a maximum value and Z(t) and F(t) are always 0 when m\geq 0.94 . Therefore, we can get the bifurcation diagram as m changes (see Figure 4). If we choose \delta as the bifurcation parameter when m = 0.75 , we can get the critical value \delta^{*} = 0.33 . Therefore, Hopf bifurcation occurs at E_{*} when \delta = \delta^* . We can obtain that P(t) reaches a maximum value of 1 and Z(t) and F(t) are always 0 when \delta\geq 0.49 . Therefore, we can get the bifurcation diagram as \delta changes (see Figure 5). From Figure 5, when other parameter values are fixed, the density of zooplankton and fish will decrease to 0 whether the refuge capacity of phytoplankton or the probability of toxin release of phytoplankton-produced toxic substances increase to some certain value. Properly increasing the shelter capacity of phytoplankton and the rate of toxin release of by phytoplankton can stabilize the population and reach a stable state. Then, the plankton and fish populations will always exist.
For the model (3.1), we assume that P(0) = 0.5, Z(0) = 30 and F(0) = 0.115 . When m = 0.8 and \delta = 0.25 , we have that \tau_{0} = 4.9397 , and the model (2.1) has one unique positive equilibrium E_{*}(0.2671, 35.1379, 0.1197) according to Theorem 2.1. From Theorem 3.3, E_{*} is locally asymptotically stable when \tau\in[0, 4.9397] , but Hopf bifurcation occurs when \tau\in[4.9397, +\infty) . From Eq (3.35), we can know that c_{1}(0) = -512.86-540.47i < 0, \mu_{1} = 1457 > 0, \mu_{2} = -1025.7 < 0 and T_{1} = 153.5416 > 0 . Thus, the Hopf bifurcation is supercritical, the bifurcating periodic solution is stable and the period of the bifurcating periodic solutions is increasing, which can be seen in Figure 6 ( \tau = 1 ) and Figure 7 ( \tau = 10 ). Here, we give the delay bifurcation diagram (see Figure 8). This means that, if the mature delay exceeds the critical value, the model transitions to unstable from stable. At this moment, the model has a Hopf bifurcation near the equilibrium and unstable behavior occurs among populations. In other words, the presence of the mature delay can destabilize the plankton-fish population.
For the model (1.6), we choose l = 1 , that is, x\in(0, \pi) . The values of other parameters as follows: d_{1} = 0.2, d_{2} = 0.5, d_{3} = 0.1, m = 0.84 and \delta = 0.25 . The positive equilibrium is E_{*}(0.3794, 38.9575, 0.1202) . From Eq (4.8), we have that \tau_{n0} = 6.2832 . Based on Theorem 4.1, E_{*} is also locally asymptotically stable when \tau = 0 (see Figure 9), and E_{*} is locally asymptotically stable when \tau = 3.5 < \tau_{n0} = 6.2832 (see Figure 10). But, E_{*} is unstable when \tau = 25 > \tau_{n0} = 6.2832 (see Figure 11). And, we can compute that c_{2}(0) = -1141.6- 2543.9i < 0, \mu_{3} = 7257.8 > 0, \mu_{4} = -2283.3 < 0 and T_{2} = 646.8790 > 0 . From Theorem 4.5, the Hopf bifurcation is supercritical, the bifurcating periodic solution is stable and the period of the bifurcating periodic solutions is increasing. For the reaction-diffusion model, we can know that the model will transitions to unstable from stable if the mature delay exceeds the critical value. At this moment, the model has a spatially homogeneous Hopf bifurcation or spatially inhomogeneous Hopf bifurcation near the equilibrium and unstable behavior occurs between the populations. At this time, the presence of the mature delay can destabilize the plankton-fish population.
By the theoretical conclusions and numerical simulation, we not only find that the existence of delay will deteriorate the system stability under some conditions, but also that the refuge of the prey and the release of toxins will cause the stability of system be damaged in a reaction-diffusion model with delay, even causing Hopf bifurcation to occur at the positive equilibrium E_{*} .
In our paper, we establish a phytoplankton-zooplankton-fish model with mature delay and population diffusion by considering the refuge of phytoplankton, C-M functional response and Holling II functional response. In [18], the authors found that the refuge affects the stability of the positive equilibrium. However, in our paper, we not only analyzed the effect of refuge, but also studied the effects of diffusion and delay on the model; we obtained that the stability of the system may be destroyed due to the existence of delay. In [19], the authors analyzed Hopf bifurcation caused by delay. However, in our paper, we not only obtained the properties of Hopf bifurcation induced by delay, but also the influence of prey refuge on the population. We determined that the existence of prey refuge can also lead to Hopf bifurcation.
After the parameters in the model were selected, the existence and stability of the equilibrium were analyzed. First, we chose m as a bifurcation parameter to study the dynamical behavior as m changes in the model (2.1). Meanwhile, we consider the effect of the parameter \delta on the positive equilibrium in the model (2.1). Through analysis, it could be obtained that the model undergoes a Hopf bifurcation when m = m^* or \delta = \delta^* . We found that, when other parameter values are fixed, the densities of zooplankton and fish will decrease to 0 regardless whether the refuge capacity of phytoplankton or the probability of toxin release of phytoplankton-produced toxic substances increase to a certain value. And, we chose the time delay \tau as the bifurcation parameter and discussed the dynamical behavior of the model without diffusion, or with diffusion, respectively. We give the direction of the Hopf bifurcation and the stability of the bifurcating periodic solution by the center manifold theorem and normal form theory. We found that the model transitions to unstable from stable when the mature delay exceeds the critical value. At this moment, the model has a spatially homogeneous Hopf bifurcation or spatially inhomogeneous Hopf bifurcation near the positive equilibrium and unstable behavior occurs between the populations. In a word, the existence of time delay has a great influence on such a model. Meanwhile, we used Matlab software for numerical simulation to prove our theoretical results.
In this paper, we have discussed the influence of factors such as prey refuge, the disturbance between predators, time delay and diffusion on the model. However, in nature, there are external factors to influence the model, such as changing temperature, environmental pollution, human activities and noise. We did not take these influencing factors into account. Therefore, in the future work, we will introduce the influence of environmental pollution on the model and analyze the dynamical behavior of the phytoplankton-zooplankton model under the influence of environmental pollution. The model is
\begin{equation*} \begin{cases} \frac{\partial P}{\partial t} = d_{1}\Delta P+r_{1}P(1-\frac{P}{K_{1}})-\frac{\beta_{1}(1-m)PZ}{1+a_{1}(1-m)P+cZ}-m_{1}P^3,\; &x\in\Omega,t > 0,\\ \frac{\partial Z}{\partial t} = d_{2}\Delta Z+r_{2}Z(1-\frac{Z}{K_{2}})+\frac{\beta_{2}(1-m)PZ}{1+a_{1}(1-m)P+cZ}-\frac{\delta P(t-\tau)Z}{a_{2}+P(t-\tau)}-m_{2}Z^2-gZ,\; &x\in\Omega,t > 0,\\ P_{x}(x,t) = Z_{x}(x,t) = 0,\; &x\in\partial\Omega,t > 0,\\ P(x,t) > 0,Z(x,t) > 0,\; &x\in\Omega, t\in[-\tau,0], \end{cases} \end{equation*} |
where m_{1} and m_{2} are the effects of environmental pollution on phytoplankton and zooplankton, respectively. We leave this work for the future.
This work was supported by the National Natural Science Foundation of China (Grant Nos. 12161054, 11661050 and 11861044), the National Natural Science Foundation of Gansu Province (Grant No. 20JR10RA156), the Doctoral Foundation of Lanzhou University of Technology and the HongLiu First-Class Disciplines Development Program of Lanzhou University of Technology.
The authors declare that they have no conflict of interest.
[1] |
J. H. Silverman, Taxicabs and sums of two cubes, Amer. Math. Mon., 100 (1993), 331–340. http://doi.org/10.1080/00029890.1993.11990409 doi: 10.1080/00029890.1993.11990409
![]() |
[2] | Z. G. Li, P. Z. Yuan, On the number of solutions of some special simultaneous Pell equations, Acta Math. Sin. Chinese Ser., 50 (2007), 1349–1356. |
[3] |
B. Li, The solution structure of multivariate linear indeterminate equation and its application, J. Anhui Univ., 39 (2015), 6–12. http://doi.org/10.3969/j.issn.1000-2162.2015.05.002 doi: 10.3969/j.issn.1000-2162.2015.05.002
![]() |
[4] |
C. X. Zhu, Y. L. Feng, S. F. Hong, J. Y. Zhao, On the number of zeros of to the equaton f(x_{1})+\cdots+f(x_{n}) = a over finite fields, Finite Fields Their Appl., 76 (2021), 101922. https://doi.org/10.1016/j.ffa.2021.101922 doi: 10.1016/j.ffa.2021.101922
![]() |
[5] |
J. Y. Zhao, Y. Zhao, Y. J. Niu, On the number of solutions of two-variable diagonal quartic equations over finite fields, AIMS Math., 5 (2020), 2979–2991. https://doi.org/10.3934/math.2020192 doi: 10.3934/math.2020192
![]() |
[6] |
S. Perlis, Normal bases of cyclic fields of prime-power degree, Duke Math. J., 9 (1942), 507–517. http://doi.org/10.1215/S0012-7094-42-00938-4 doi: 10.1215/S0012-7094-42-00938-4
![]() |
[7] |
D. Pei, C. Wang, J. Omura, Normal basis of finite field GF(2^{m}), IEEE Trans. Inf. Theory, 32 (1986), 285–287. http://doi.org/10.1109/TIT.1986.1057152 doi: 10.1109/TIT.1986.1057152
![]() |
[8] |
Y. Chang, T. K. Truong, I. S. Reed, Normal bases over GF(q), J. Algebra, 241 (2001), 89–101. https://doi.org/10.1006/jabr.2001.8765 doi: 10.1006/jabr.2001.8765
![]() |
[9] |
H. Huang, S. M. Han, W. Cao, Normal bases and irreducible polynomials, Finite Fields Their Appl., 50 (2018), 272–278. https://doi.org/10.1016/j.ffa.2017.12.004 doi: 10.1016/j.ffa.2017.12.004
![]() |
[10] | S. H. Gao, Normal bases over finite fields, University of Waterloo, 1993. |
[11] |
W. Cao, Normal bases and factorization of x^n-1 in finite fields, Appl. Algebra Eng. Commun. Comput., 35 (2024), 167–175. http://doi.org/10.1007/s00200-022-00540-z doi: 10.1007/s00200-022-00540-z
![]() |
[12] | R. Lidl, H. Niederreiter, Finite fields, Cambridge University Press, 1997. http://doi.org/10.1017/CBO9780511525926 |
[13] | Z. Ke, Q. Sun, Lectures on number theory, Higher Education Press, 2001. |
[14] | K. Ireland, M. Rosen, A classical introduction to modern number theory, Springer-Verlag, 1990. https://doi.org/10.1007/978-1-4757-2103-4 |
[15] | A. H. Parvardi, Lifting the exponent lemma (LTE), Art of Problem Solving, 2011. |
1. | Muhammad Uzair Awan, Sadia Talib, Muhammad Aslam Noor, Yu-Ming Chu, Khalida Inayat Noor, On post quantum estimates of upper bounds involving twice (p,q)-differentiable preinvex function, 2020, 2020, 1029-242X, 10.1186/s13660-020-02496-5 | |
2. | Shu-Bo Chen, Saima Rashid, Muhammad Aslam Noor, Zakia Hammouch, Yu-Ming Chu, New fractional approaches for n-polynomial P-convexity with applications in special function theory, 2020, 2020, 1687-1847, 10.1186/s13662-020-03000-5 | |
3. | Yu-Ming Chu, Muhammad Uzair Awan, Sadia Talib, Muhammad Aslam Noor, Khalida Inayat Noor, New post quantum analogues of Ostrowski-type inequalities using new definitions of left–right (p,q)-derivatives and definite integrals, 2020, 2020, 1687-1847, 10.1186/s13662-020-03094-x | |
4. | Saad Ihsan Butt, Muhammad Umar, Saima Rashid, Ahmet Ocak Akdemir, Yu-Ming Chu, New Hermite–Jensen–Mercer-type inequalities via k-fractional integrals, 2020, 2020, 1687-1847, 10.1186/s13662-020-03093-y | |
5. | Hua Wang, Humaira Kalsoom, Hüseyin Budak, Muhammad Idrees, Ahmet Ocak Akdemir, q -Hermite–Hadamard Inequalities for Generalized Exponentially s , m ; η -Preinvex Functions, 2021, 2021, 2314-4785, 1, 10.1155/2021/5577340 | |
6. | Humaira Kalsoom, Muhammad Aamir Ali, Muhammad Idrees, Praveen Agarwal, Muhammad Arif, Adrian Neagu, New Post Quantum Analogues of Hermite–Hadamard Type Inequalities for Interval-Valued Convex Functions, 2021, 2021, 1563-5147, 1, 10.1155/2021/5529650 | |
7. | Humaira Kalsoom, Miguel Vivas-Cortez, Muhammad Idrees, Praveen Agarwal, New Parameterized Inequalities for η-Quasiconvex Functions via (p, q)-Calculus, 2021, 23, 1099-4300, 1523, 10.3390/e23111523 | |
8. | Saad Ihsan Butt, Hüseyin Budak, Muhammad Tariq, Muhammad Nadeem, Integral inequalities for n ‐polynomial s ‐type preinvex functions with applications , 2021, 44, 0170-4214, 11006, 10.1002/mma.7465 | |
9. | Humaira Kalsoom, Miguel Vivas-Cortez, Muhammad Amer Latif, Trapezoidal-Type Inequalities for Strongly Convex and Quasi-Convex Functions via Post-Quantum Calculus, 2021, 23, 1099-4300, 1238, 10.3390/e23101238 | |
10. | YONGFANG QI, GUOPING LI, FRACTIONAL OSTROWSKI TYPE INEQUALITIES FOR (s,m)-CONVEX FUNCTION WITH APPLICATIONS, 2023, 31, 0218-348X, 10.1142/S0218348X23501281 |