Processing math: 63%
Research article Special Issues

On the number of the irreducible factors of xn1 over finite fields

  • Received: 19 June 2024 Revised: 25 July 2024 Accepted: 26 July 2024 Published: 05 August 2024
  • MSC : 11D04, 11T06

  • 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,,αqn1}. Some problems on normal bases can be finally reduced to the determination of the irreducible factors of the polynomial xn1 in Fq, while the latter is closely related to the cyclotomic polynomials. Denote by F(xn1) the set of all distinct monic irreducible factors of xn1 in Fq. The criteria for

    |F(xn1)|2

    have been studied in the literature. In this paper, we provide the sufficient and necessary conditions for

    |F(xn1)|=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(xn1)|=3,4,5.

    Citation: Weitao Xie, Jiayu Zhang, Wei Cao. On the number of the irreducible factors of xn1 over finite fields[J]. AIMS Mathematics, 2024, 9(9): 23468-23488. doi: 10.3934/math.20241141

    Related Papers:

    [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,,αqn1}. Some problems on normal bases can be finally reduced to the determination of the irreducible factors of the polynomial xn1 in Fq, while the latter is closely related to the cyclotomic polynomials. Denote by F(xn1) the set of all distinct monic irreducible factors of xn1 in Fq. The criteria for

    |F(xn1)|2

    have been studied in the literature. In this paper, we provide the sufficient and necessary conditions for

    |F(xn1)|=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(xn1)|=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(1PK)β1(Pm)Za1+(Pm),dZdt=β2(Pm)Za1+(Pm)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 Pm 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:

    {Pt=d1ΔP+rP(1PK)μPZα+P,Zt=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 (1m)P of the unprotected phytoplankton available for zooplankton grazing, following the Holling II functional response, i.e., f1(P)=(1m)Pa1+(1m)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(1PK)β1(1m)PZa1+(1m)P,x(0,lπ),t>0,Z(x.t)t=d2ΔZ+β2(1m)PZa1+(1m)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)
    Figure 1.  Diagram of interactions among phytoplankton, zooplankton and fish populations.

    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

    {Pt=d1ΔP+P(1P)β1(1m)PZα1+(1m)P,x(0,lπ),t>0,Zt=d2ΔZ+β2(1m)PZα1+(1m)PδP(tτ)Zα2+P(tτ)γ1ZF(a+Z)(b+F)g1Z,x(0,lπ),t>0,Ft=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(1P)β1(1m)PZα1+(1m)P,dZdt=β2(1m)PZα1+(1m)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(1m)(β2δg1)g1α1α22(1m)(β2δg1),Z2=(1P2)[α1+(1m)P2](1m)β1;

    here, B=(β2g1)(1m)α2(δ+g1)α.

    Assume that the model (2.1) has the coexistence equilibrium E(P,Z,F), where P, Z and F satisfy

    {P(1P)β1(1m)PZα1+(1m)P=0,β2(1m)PZα1+(1m)PδPZα2+Pγ1ZF(a+Z)(b+F)g1Z=0,γ2ZF(a+Z)(b+F)g2F=0. (2.2)

    By the first equation of (2.2), we can obtain

    Z=(1P)[α1+(1m)P](1m)β1. (2.3)

    Substituting Eq (2.3) into the third equation of (2.2), we have

    F=(γ2bg2)(1P)[α1+(1m)P]bg2aβ1(1m)g2[aβ1(1m)+(1P)(α1+(1m)P)]. (2.4)

    If (H2):(γ2bg2)(1P)[α1+(1m)P]bg2aβ1(1m)>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=(1m)2γ2(β2δg1),M4=γ2(1m)(1mα1)(β2δg1)+δγ2(1m)(1mα1)γ2(1m)2(1α2)(β2g1)g1γ2α1(1m),M3=[aβ1(1m)+α1]γ2(1m)(β2δg1)+(β2g1γ2(1m)(1α2)(1mα1)δγ2(1mα1)2+g1γ2α1(1mα1)+γ1(1m)2(γ2bg2)α2β2γ2(1m)2+γ2(1m)α1(δ+g1)+g1γ2α2(1m)(1mα1),M2=α1(g1γ2α2+γ1β1)(1m)+α2β2γ2(1m)(1mα2)+(1m)(1α2)α1γ2(β2g1)γ2α1(1mα1)(2δ+g1)+aβ1g1γ2α1(1m)+g1γ2α2aβ1δγ2(1m)(1mα1)+β1(1m)2(1α2)[aβ2γ2ag1γ2γ1(γ2bg2)],M1=[aβ1(1m)+α1][γ2α2(β2g1)+α1g1γ2α2γ2α1(δ+g1)](1mα1)g1γ2α1α2+γ1abg2β21(1m)2γ1β1(1m)(γ2bg2)[α2(1mα1)+α1],M0=γ1abg2β21α2(1m)2[aβ1(1m)+α1]g1γ2α1α2γ1β1(1m)(γ2bg2)α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.

    Figure 2.  Existence and uniqueness of the positive roots of f(P).

    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=12Pα1β1(1m)Z[α1+(1m)P]2,a12=β1(1m)Pα1+(1m)P,a21=α1β2(1m)Z[α1+(1m)P]2δα2Z(α2+P)2,a23=γ1bZ(a+Z)(b+F)2,a22=β2(1m)Pα1+(1m)PδPα2+Pγ1aF(a+Z)2(b+F)g1,a32=γ2aF(a+Z)2(b+F),a33=γ2bZ(a+Z)(b+F)2g2.

    The characteristic equation of the model (2.1) is

    λ3(a11+a22+a33)λ2+(a11a22+a11a33+a22a33)λ+a11a32a23a11a22a33+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(1m)α1+1mδα2+1g1<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(1m)α1+(1m)δα2+1g1)](λ+g2)=0.

    It has three roots:

    λ1=1<0,λ2=β2(1m)α1+(1m)δα2+1g1,λ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(1P2)[α1+(1m)P2]abβ1(1m)+b(1P2)[α1+(1m)P2]g2,s2=α1β1(1m)Z2+(1m)P22β2(1m)P2α1+(1m)P2,s3=1(α2+P2)[α1+(1m)P2]2[(12P2)(α1+(1m)P2)α1(1P2)][β2(1m)P2(α2+P2)(α1+(1m)P2)((δ+g1)P2+α2g1)]+α1P2(1(1m)P2)[α1+(1m)P2]2δα2P2(1P2)(α+P2)2.

    Assume that the condition (H5):s1<0, s224s3>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+A11A22A12A21A23A32,M8=A11A23A32+A12A21A33A11A22A33,A11=P+β1(1m)2Z[α1+(1m)P]2,A12=β1(1m)Pα1+(1m)P,A21=α1β2(1m)Z[α1+(1m)P]2δα2Z(α2+P)2,A22=γ1FZ(a+Z)2(b+F),A23=γ1bZ(a+Z)(b+F)2,A32=γ2aF(a+Z)2(b+F),A33=γ2FZ(a+Z)(b+F)2.

    From the Routh-Hurwitz criterion[54], if M6>0, M7>0 and M6M7M8>0, then all solutions of Eq (2.8) have negative real parts. When M6>0, M7>0 and M6M7M8<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 M6M7M8>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 M6M7M8 passes through 0; in other words, the model (2.1) undergoes a Hopf bifurcation at E when M6M7M8=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=m0,

    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(1m)2Zα1[α1+(1m)P]δ2[α2+P]2>0 and [α1+(1m)P]β1Pδβ2α1(α2+P)2β21(1m)Pγ2α1g1ab>0.

    Proof. Let (P,Z,F) be any positive solution of the model (2.1). Define a Lyapunov function

    V(t)=PPPlnPP+[α1+(1m)P]β1β2α1(ZZZlnZZ)+FFFlnFF.

    Calculating the derivative of V(t) along the solution of the model (2.1), then we have

    dV(t)dt=(PP){(PP)+β1(1m)(1m)Z(PP)[α1+(1m)P](ZZ)[α1+(1m)P][α1+(1m)P]}+[α1+(1m)P]β1β2α1(ZZ){β2(1m)α1[α1+(1m)P][α1+(1m)P](PP)+β2(1m)PZ[α1+(1m)P](ZZ)δα2(α2+P)(α2+P)(PP)Pδα2+P(ZZ)γ1b(a+Z)(b+F)(b+F)(FF)γ1aFZ(a+Z)(a+Z)(b+F)(ZZ)g1Z(ZZ)}+(FF){γ2a(a+Z)(b+F)(a+Z)(ZZ)γ2Z(b+F)(a+Z)(b+F)(FF)}{1β1(1m)2Z[α1+(1m)P][α1+(1m)P]δ2α22[α2+P]2[α2+P]2}(PP)2{2β1(1m)Pα1Z+[α1+(1m)P]β1Pδβ2α1(α2+P)+γ1aFZ(a+Z)(a+Z)(b+F)+g1Z}(ZZ)2{γ2Z(b+F)(a+Z)(b+F)+γ22a2(a+Z)2(b+F)2(a+Z)2}(FF)2{1β1(1m)2Zα1[α1+(1m)P]δ2[α2+P]2}(PP)2{2β21(1m)Pγ2α1g1ab+[α1+(1m)P]β1Pδβ2α1(α2+P)}(ZZ)2{γ2Z(b+F)(a+Z)(b+F)+γ22a2(a+Z)2(b+F)2(a+Z)2}(FF)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(1m)2Zα1[α1+(1m)P]δ2[α2+P]2>0 and 2β21(1m)Pγ2α1g1ab+[α1+(1m)P]β1Pδβ2α1(α2+P)>0 hold, then we have that the coefficients of (PP)2, (ZZ)2 and (FF)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(1P)β1(1m)PZα1+(1m)P,dZdt=β2(1m)PZα1+(1m)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(1m)[α1+P(1m)]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+b22b33b12b21b23b32,D0=b11b23b32+b12b21b33b11b22b33,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ω3D2ω3+D1iω+D0+(cosωτisinωτ)(iD3ω+D4)=0. (3.6)

    Separating the real and imaginary parts of Eq (3.6), we can obtain

    {ω3D1ω=D3ωcosωτD4sinωτ,D2ω2D0=D3ωsinωτ+D4cosωτ. (3.7)

    From Eq (3.7), we can get

    ω6+ω4(D222D1)+ω2(D212D0D2D33)+D20D24=0. (3.8)

    Let z=ω2, D5=D222D1, D6=D212D0D1D23 and D7=D20D24. 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=D253D60 and (H8):Δ1=D253D6>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 D70 and (H7) hold, Eq (3.9) has no positive root for z[0,+).

    When D70 and (H8) hold, Eq (3.10) has two roots, that is, z1 and z2, where

    z1=D5+Δ13,z2=D5Δ13.

    Furthermore, we have

    f(z1)=2D5+2Δ1+2D5=2Δ1>0,f(z2)=2D52Δ1+2D5=2Δ1<0.

    Therefore, we can obtain z1 and z2 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 D70, then Eq (3.9) has no positive root.

    (3) If the condition (H8) holds and D70, then Eq (3.9) has two positive roots when z1>0 and f(z1)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+(D2D4D1D3)ω2kD0D4D23ω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.jN0{τ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ω00.

    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ω20D12D2ω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+(2D224D1)ω40+(D212D0D2D23)ω20D23ω40+D24ω20=zkD23ω40+D24ω20[3z2k+(2D224D1)zk+(D212D0D2D23)]=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 D70 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, z1>0 and f(z1)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 ttτ; 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μ:CR3 and F:R×CR3.

    Define φ(θ)=(φ1(θ),φ2(θ),φ3(θ))TR3,θ[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+(1m)2α1β1Z[α1+P(1m)]3,C12=α1β1(1m)[α1+P(1m)]2,C21=(1m)2α1β2Z[α1+P(1m)]3,C22=γ1aF(a+Z)3(b+F),C23=γ1bZ(a+Z)(b+F)3,C24=α1β2(1m)[α1+P(1m)]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μ(φ)=01dη(θ,μ)φ(θ),φ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),01dη(μ,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),01ψ(ξ)dη(ξ,0),s=0,

    and establish the bilinear inner product

    ψ(s),φ(θ)=ˉψ(0)φ(0)01θ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)=01dη(0)q(0)=A(0)=iωτ0q(0),θ=0,

    that is,

    (τ0+μ)(b11iω0b120b21+ceiω0τ0b22iω0b230b32b33iω0)(1q1q2)=0.

    Solving it, we can get

    (1q1q2)=(1iω0b11b12b32(iω0b11)b12(iω0b33)).

    Similarly, we assume that q(s)=D(1,q1,q2)eiω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

    (1q1q2)=(1iω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)01θ0¯q(ξθ)dη(θ)q(ξ)dξ=ˉD[(1,q1,q2)(1,q1,q2)T01θ0(1,¯q1,¯q2)eiω0(ξθ)τ0dη(θ)(1,q1,q2)Teiξω0τ0dξ]=ˉD[1+q1¯q1+q2¯q2+τ0c¯q1eiω0τ0].

    Thus, we can choose

    ˉD=[1+q1¯q1+q2¯q2+τ0c¯q1eiω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 utC0 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={AW2Re{¯q(0)F0q(θ)},θ[1,0),AW2Re{¯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

    (A2iω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)Teiω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)+¯q1F2(ut)+¯q2F3(ut)]=z22[2ˉDτ0(k11+k21¯q1+k31¯q2)]+zˉz[ˉDτ0(k12+k22¯q1+k32¯q2)]+ˉz22[2ˉDτ0(k13+k23¯q1+k33¯q2)]+z2ˉz2[2ˉDτ0(k14+k24¯q1+k34¯q2)], (3.26)

    where

    k11=C11+C12q1,k12=2C11+C12(ˉq1+q1),k21=C21+C22q21+C23q22+C24q1+C25q1q2+C26e2iω0τ0+C27q1eiω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+ˉq1eiω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)+2eiω0τ0W(1)11(1))+C27(12ˉq1W(1)20(1)+q1W(1)11(1)+eiω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+ˉq1k21+ˉq2k31),g11=ˉDτ0(k12+ˉq1k22+ˉq2k32),g02=2ˉDτ0(k13+ˉq1k23+ˉq2k33),g21=2ˉDτ0(k14+ˉq1k24+ˉq2k34). (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τ0eiω0τ0θ+E1e2iω0τ0θ,W11(θ)=ig11q(0)ω0τ0eiω0τ0θ+iˉg11ˉq(0)ω0τ0eiω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

    01dη(θ)W20(θ)=2iω0τ0W20H20(0),01dη(θ)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τ0H20(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ω0b11b120ce2iω0τ0b212iω0b22b230b322iω0b33)1(k11k21k31),
    E2=(b11b120cb21b22b230b32b33)1(k11k21k31).

    Thus, all expressions of gij can be represented in full. Also, we have

    {c1(0)=i2ω0τ0(g11g202|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 nN0 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+b11b33b12b21b23b32),q0n=d1d2d3(nl)6(d1d2b33+d1d3b22+d2d3b11)(nl)4+(d1b22b33+d2b11b33+d3b11b22)(nl)2+b11b23b32+b12b21b33b11b22b33,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+b12b21b33b11b22b33,
    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+b11b33b12b21b23b32cb12)+(b11+b22+b33)(d1b22+d2b11+d2b33+d3b22+d3b11+d1b33)(d1b22b33+d2b11b33+d3b11b22)+d3cb12](nl)2+(b11+b22+b33)(cb12b11b22b22b33b11b33+b12b21+b23b32)(b11b23b32+b12b21b33b11b22b33)cb12b33.

    Assume that the following condition holds true: (H9): q2n>0, q2n(q1n+q3n)(q0n+q4n)>0, q0n+q4n>0, nN0. 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+eiω2τ(iω2q3n+q4n)=0. (4.4)

    Separating the real and imaginary parts of Eq (4.4), we can receive

    {ω32q1nω2=q3nω2cosω2τq4nsinω2τ,q2nω22q0n=q3nω3nsinω2τ+q4ncosω2τ, (4.5)

    which follows that

    ω62+(q22n2q1n)ω42+(q21n2q0nq2nq23n)ω22+(q20nq24n)=0. (4.6)

    Let p2n=q22n2q1n, p1n=q21n2q0nq2nq23n, p0n=q20nq24n 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=p22n3p1n0 and (H11):Δ2=p22n3p1n>0 are true, here y1=p1n+Δ23 and y2=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 p0n0 and the condition (H10) holds, then Eq (4.7) has no positive root.

    (3) If p0n0 and the condition (H11) holds, then Eq (4.7) has positive roots when y1>0 and f(y1)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+(q2nq4nq1nq3n)(ωk2)2q0nq4nq23n(ω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).

    Figure 3.  Positive equilibrium E_{*}(0.5043, 31.1137, 0.1191) of the model (2.1) is locally asymptotically stable when \tau = 0 . (a) P(t), (b) Z(t), (c) F(t).

    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.

    Figure 4.  When m\in(0, 1) , the dynamical behavior of the model (2.1) changes. (a) P(t), (b) Z(t), (c) F(t).
    Figure 5.  When \delta\in(0, 1) , the dynamical behavior of the model (2.1) changes. (a) P(t), (b) Z(t), (c) F(t).

    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.

    Figure 6.  Positive equilibrium E_{*}(0.2671, 35.1379, 0.1197) of the model (2.1) is locally asymptotically stable when \tau = 1 < \tau_{0} = 4.9397 . (a) P(t), (b) Z(t), (c) F(t).
    Figure 7.  Hopf bifurcation occurs at the positive equilibrium when \tau = 10 > \tau_{0} = 4.9397 . (a) P(t), (b) Z(t), (c) F(t).
    Figure 8.  Bifurcation diagrams of the model (2.1) with respect to \tau .

    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.

    Figure 9.  Positive equilibrium E_{*}(0.3794, 38.9575, 0.1202) of the model (1.6) is locally asymptotically stable when \tau = 0 . (a) P(t), (b) Z(t), (c) F(t).
    Figure 10.  Positive equilibrium E_{*}(0.3794, 38.9575, 0.1202) of the model (1.6) is locally asymptotically stable when \tau = 3.5 < \tau_{n0} = 6.2832 . (a) P(t), (b) Z(t), (c) F(t).
    Figure 11.  Positive equilibrium E_{*}(0.3794, 38.9575, 0.1202) of the model (1.6) is unstable and Hopf bifurcation occurs when \tau = 25 > \tau_{n0} = 6.2832 . (a) P(t), (b) Z(t), (c) F(t).

    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.
  • This article has been cited by:

    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
  • Reader Comments
  • © 2024 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

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

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

Metrics

Article views(977) PDF downloads(90) Cited by(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog