Dynamical Analysis in a Delayed Predator-Prey Model with Two Delays

Abstract

A class of Beddington-DeAngelis functional response predator-prey model is considered. The conditions for the local stability and the existence of Hopf bifurcation at the positive equilibrium of the system are derived. Some explicit formulae for determining the stability and the direction of the Hopf bifurcation periodic solutions bifurcating from Hopf bifurcations are obtained by using the normal form theory and center manifold theory. Some numerical simulations for justifying the theoretical analysis are also provided. Finally, main conclusions are given.

1. Introduction

In recent years, population dynamics (including stable, unstable, persistent, and oscillatory behavior) has become very popular since Vito Volterra and James Lotka proposed the seminal models of predator-prey models in the mid-1920s. Great attention has been paid to the dynamics properties of the predator-prey models which have significant biological background. Many excellent and interesting results have been obtained [1–21]. In 2009, Gakkhar et al. [4] investigated the local stability and Hopf bifurcation of the autonomous delayed predator-prey system with Beddington-DeAngelis functional response:𝑢̇₁(𝑡)=𝑢₁(𝑡)1−𝑢₁𝑡−𝜏₁−𝑎₁𝑢₁(𝑡)𝑢₂(𝑡)𝑎+𝑢₁(𝑡)+𝑏𝑢₂(𝑡),𝑢̇₂(𝑡)=𝑑₁𝑢₂(𝑡𝑎)−𝑑+¹𝑢₁𝑡−𝜏₂𝑎+𝑢₁𝑡−𝜏₂+𝑏𝑢₂𝑡−𝜏₂, (1.1) where 𝑢₁(𝑡),𝑢₂(𝑡) represent the prey density and the predator density, respectively. The delay terms occur in growth as well as interaction terms. For this, it means that the prey takes time 𝜏₁ to convert the food into its growth [11], whereas the predator takes time 𝜏₂ for the same [22]. All the parameters in the model take positive values, that is, 𝑎₁>0,𝑑>0,𝑑₁>0,𝑎>0,𝑏>0. The more detail biological meaning of the coefficients of system (1.1), one can see [11] or [22].

We would like to point out that Gakkhar et al. [4] studied the local stability and Hopf bifurcation of system (1.1) under the assumption: 𝜏₁=𝜏₂=𝜏 and obtained some excellent results. While in most cases, 𝜏₁≠𝜏₂. Considering the factor, we further investigate the model (1.1) with 𝜏₁≠𝜏₂ as a complementarity.

In this paper, we go on to study the stability, the local Hopf bifurcation for system (1.1). To the best of our knowledge, it is the first time to deal with the research of Hopf bifurcation for model (1.1) under the assumption 𝜏₁≠𝜏₂.

The remainder of the paper is organized as follows. In Section 2, we investigate the stability of the positive equilibrium and the occurrence of local Hopf bifurcations. In Section 3, the direction and stability of the local Hopf bifurcation are established. In Section 4, numerical simulations are carried out to illustrate the validity of the main results. Some main conclusions are drawn in Section 5.

2. Stability of the Positive Equilibrium and Local Hopf Bifurcations

In this section, we shall study the stability of the positive equilibrium and the existence of local Hopf bifurcations.

Since time delay does not change the equilibrium of system and according to [4], we know that the delayed prey predator model (1.1) has three equilibrium points: two boundary equilibrium 𝐸₁(0,0) and 𝐸₂(1,0), and a nontrivial equilibrium point 𝐸₀(𝑢^∗₁,𝑢^∗₂), where 𝑢^∗₁,𝑢^∗₂ are the positive solutions of the following quadratic equations:𝑢²₁+𝛼₁𝑢₁+𝛽₁𝑢=0,₂₂+𝛼₂𝑢₂+𝛽₂=0, (2.1) where𝛼₁=𝑎₁−𝑑𝑎₁−𝑏𝑏,𝛽₁=−𝑎𝑑𝑎₁𝑑,𝛼₂=𝑎₁(𝑑−1)²+𝑏(2𝑎𝑑+𝑑−1)𝑏²𝑑,𝛽₂𝑎()=−𝑎𝑑+𝑑−1𝑏²𝑑. (2.2) Since 𝛽₁<0, (2.1) in 𝑢₁ admits a unique positive solution. If one of the following conditions:(a)𝛼₂>0,𝛽₂<0,(b)𝛼₂<0,𝛽₂=𝛼²₂4,(c)𝛼₂<0,0<𝛽₂<𝛼²₂4, (2.3) holds, then system (2.1) has at least one positive equilibrium point 𝐸₀(𝑢^∗₁,𝑢^∗₂).

Let 𝑢₁(𝑡)=𝑢₁(𝑡)−𝑢^∗₁,𝑢₂(𝑡)=𝑢₂(𝑡)−𝑢^∗₂ and still denote 𝑢_𝑖(𝑡)(𝑖=1,2) by 𝑢_𝑖(𝑡)(𝑖=1,2), respectively, then (1.1) becomes𝑢̇₁(𝑡)=𝑚₁𝑢₁(𝑡)+𝑚₂𝑢₂(𝑡)+𝑚₃𝑢₁𝑡−𝜏₁+𝐹₁,𝑢̇₂(𝑡)=𝑛₁𝑢₂(𝑡)+𝑛₂𝑢₁𝑡−𝜏₂+𝑛₃𝑢₂𝑡−𝜏₂+𝐹₂, (2.4) where 𝑚_𝑖,𝑛_𝑖(𝑖=1,2,3) and 𝐹_𝑗(𝑗=1,2) are defined by Appendix A.

The linearization of (2.4) at (0,0) is𝑢̇₁(𝑡)=𝑚₁𝑢₁(𝑡)+𝑚₂𝑢₂(𝑡)+𝑚₃u₁𝑡−𝜏₁,𝑢̇₂(𝑡)=𝑛₁𝑢₂(𝑡)+𝑛₂𝑢₁𝑡−𝜏₂+𝑛₃𝑢₂𝑡−𝜏₂, (2.5) whose characteristic equation is𝜆²−𝑚₁+𝑛₁𝜆+𝑚₁𝑛₁−𝑛₃𝜆−𝑚₁𝑛₃+𝑚₂𝑛₂𝑒^−𝜆𝜏2−𝑚₃𝜆+𝑚₃𝑛₁𝑒^−𝜆𝜏1+𝑚₃𝑛₃𝑒^{−𝜆(𝜏1+𝜏2)}=0. (2.6) In order to investigate the distribution of roots of the transcendental equation (2.6), the following Lemma is useful.

Lemma 2.1 (see [23]). For the transcendental equation 𝑃𝜆,𝑒^−𝜆𝜏1,…,𝑒^{−𝜆𝜏𝑚}=𝜆^𝑛+𝑝₁⁽⁰⁾𝜆^𝑛−1+⋯+𝑝⁽⁰⁾_𝑛−1𝜆+𝑝_𝑛⁽⁰⁾+𝑝₁⁽¹⁾𝜆^𝑛−1+⋯+𝑝⁽¹⁾_𝑛−1𝜆+𝑝_𝑛⁽¹⁾𝑒^−𝜆𝜏1+𝑝+⋯_1(𝑚)𝜆^𝑛−1+⋯+𝑝^(𝑚)_𝑛−1𝜆+𝑝_𝑛^(𝑚)𝑒^{−𝜆𝜏𝑚}=0, (2.7) as (𝜏₁,𝜏₂,𝜏₃,…,𝜏_𝑚) vary, the sum of orders of the zeros of 𝑃(𝜆,𝑒^−𝜆𝜏1,…,𝑒^{−𝜆𝜏𝑚}) in the open right half plane can change, and only a zero appears on or crosses the imaginary axis.

In the sequel, we consider three cases.

Case a. 𝜏₁=𝜏₂=0, (2.6) becomes 𝜆²−𝑚₁+𝑚₃+𝑛₁+𝑛₃𝜆+𝑚₁𝑛₁+𝑚₁𝑛₃+𝑚₃𝑛₃−𝑚₂𝑛₂−𝑚₃𝑛₁=0. (2.8) A set of necessary and sufficient conditions for all roots of (2.8) to have a negative real part are given in the following form: (H1)𝑚₁+𝑚₃+𝑛₁+𝑛₃<0,𝑚₁𝑛₁+𝑚₁𝑛₃+𝑚₃𝑛₃−𝑚₂𝑛₂−𝑚₃𝑛₁>0. (2.9) Then, the equilibrium point 𝐸₀(𝑢^∗₁,𝑢^∗₂) is locally asymptotically stable when the condition (H1) holds.

Case b. 𝜏₁=0,𝜏₂>0, (2.6) becomes 𝜆²+𝑝𝜆+𝑟+𝑠𝜆+𝑞𝑒^−𝜆𝜏2=0, (2.10) where 𝑝=−𝑚₁+𝑚₃+𝑛₁,𝑟=𝑚₁𝑛₁−𝑚₃𝑛₁,𝑠=−𝑛₃,𝑞=𝑚₁𝑛₃+𝑚₃𝑛₃−𝑚₂𝑛₂. (2.11) For 𝜔>0,𝑖𝜔 be a root of (2.10), then it follows that 𝑞cos𝜔𝜏₂+𝑠𝜔sin𝜔𝜏₂=𝜔²−𝑟,𝑠𝜔cos𝜔𝜏₂−𝑞sin𝜔𝜏₂=−𝑝𝜔, (2.12) which leads to 𝜔⁴+𝑝²−𝑠²−2𝑟𝜔²+𝑟²−𝑞²=0. (2.13) It is easy to see that if the condition (H2)𝑝²−𝑠²−2𝑟>0,𝑟²−𝑞²>0 (2.14) holds, then (2.13) has no positive roots. Hence, all roots of (2.10) have negative real parts when 𝜏₂∈[0,+∞) under the conditions (H1) and (H2).
If (H1) and (H3)𝑟²−𝑞²<0 (2.15) hold, then (2.13) has a unique positive root 𝜔²₀. Substituting 𝜔²₀ into (2.12), we obtain 𝜏_2𝑛=1𝜔₀𝑞𝜔arccos²⁰−𝑟−𝑝𝑠𝜔²₀𝑠²𝜔²₀+𝑞²+2𝑛𝜋,𝑛=0,1,2,…. (2.16) If (H1) and (H4)𝑝²−𝑠²−2𝑟<0,𝑟²−𝑞²>0,𝑝²−𝑠²−2𝑟²>4𝑟²−𝑞² (2.17) hold, then (2.10) has two positive roots 𝜔²₊ and 𝜔²₋. Substituting 𝜔²_± into (2.12), we obtain 𝜏^±_2𝑘=1𝜔_±𝑞𝜔arccos^2±−𝑟−𝑝𝑠𝜔²_±𝑠²𝜔²_±+𝑞²+2𝑘𝜋,𝑘=0,1,2,…. (2.18) Let 𝜆(𝜏₂)=𝛼(𝜏₂)+𝑖𝜔(𝜏₂) be a root of (2.10) near 𝜏₂=𝜏_2𝑛 and 𝛼(𝜏_2𝑛)=0,𝜔(𝜏_2𝑛)=𝜔₀. Due to functional differential equation theory, for every 𝜏_2𝑛,𝑛=0,1,2,…, there exists 𝜀>0 such that 𝜆(𝜏₂) is continuously differentiable in 𝜏₂ for |𝜏₂−𝜏_2𝑛|<𝜀. Substituting 𝜆(𝜏₂) into the left-hand side of (2.10) and taking derivative with respect to 𝜏₂, we have 𝑑𝜆𝑑𝜏₂⁻¹=2𝜆+𝑝𝑒^𝜆𝜏2𝜆+𝑠𝑠𝜆+𝑞𝜆−𝜏𝑠𝜆+𝑞𝜆², (2.19) which leads to 𝑑(Re𝜆(𝜏))𝑑𝜏_2𝜏⁻¹_2=𝜏2𝑛𝑒=Re2𝜆+𝑝^𝜆𝜏2𝜆𝑠𝜆+𝑞_{𝜏2=𝜏2𝑛}𝑠+Re𝜆𝑠𝜆+𝑞_{𝜏2=𝜏2𝑛}=Re𝑝cos𝜔₀𝜏_2𝑛−2𝜔₀sin𝜔₀𝜏_2𝑛+𝑖2𝜔₀cos𝜔₀𝜏_2𝑛+𝑝sin𝜔₀𝜏_2𝑛−𝑠𝜔²₀+𝑖𝑞𝜔₀𝑠+Re−𝑠𝜔²₀+𝑖𝑞𝜔₀=1Λ−𝑠𝜔²₀𝑝cos𝜔₀𝜏_2𝑛−2𝜔₀sin𝜔₀𝜏_2𝑛+𝑞𝜔₀2𝜔₀cos𝜔₀𝜏_2𝑛+𝑝sin𝜔₀𝜏_2𝑛−𝑠²𝜔²₀=1Λ𝑝𝜔₀𝑞sin𝜔₀𝜏_2𝑛−𝑠𝜔₀cos𝜔₀𝜏_2𝑛+2𝜔²₀𝑞cos𝜔₀𝜏_2𝑛+𝑠𝜔₀sin𝜔₀𝜏_2𝑛−𝑠²𝜔²₀=1Λ𝑝²𝜔²₀+2𝜔⁴₀−2𝑟𝜔²₀−𝑠²𝜔²₀=𝜔²₀Λ2𝜔²₀+𝑝²−𝑠²=1−2𝑟Λ−𝑝²+𝑠²+2𝑟+Δ_∗+𝑝²−𝑠²=𝜔−2𝑟²⁰ΛΔ_∗>0, (2.20) where Λ=𝑠²𝜔⁴₀+𝑞²𝜔²₀>0,Δ_∗=𝑠²−𝑝²+2𝑟²−4𝑟²−𝑞². (2.21) Noting that 𝑑()signRe𝜆𝑑𝜏_{2𝜏2=𝜏2𝑛}=signRe𝑑𝜆𝑑𝜏_{2𝜏2=𝜏2𝑛}=1, (2.22) we have 𝑑(Re𝜆)𝑑𝜏_{2𝜏2=𝜏2𝑛}>0. (2.23) Similarly, we can obtain 𝑑(Re𝜆)𝑑𝜏_{2𝜏2=𝜏+2𝑘}>0,𝑑(Re𝜆)𝑑𝜏_{2𝜏2=𝜏−2𝑘}<0. (2.24) According to above analysis the Corollary  2.4 in Ruan and Wei [23], we have the following results.

Lemma 2.2. For𝜏₁=0, assume that one of the conditions (a), (b), (c), and (d) holds and (H1) is satisfied. Then, the following conclusions hold.(i)If (H2) holds, then the positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂) of system (1.1) is asymptotically stable for all 𝜏₂≥0.(ii) If (H3) holds, then the positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂) of system (1.1) is asymptotically stable for 𝜏₂<𝜏₂₀, and unstable for 𝜏₂<𝜏₂₀. Furthermore, system (1.1) undergoes a Hopf bifurcation at the positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂) when 𝜏₂=𝜏₂₀.(iii) If (H4) holds, then there is a positive integer 𝑚 such that the positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂) is stable when 𝜏₂∈[0,𝜏⁺₂₀)∪(𝜏⁻₂₀,𝜏⁺₂₁)∪⋯∪(𝜏⁻_2𝑚−1,𝜏⁺_2𝑚), and unstable when 𝜏₂∈[𝜏⁺₂₀,𝜏⁻₂₀)∪(𝜏⁺₂₁,𝜏⁻₂₁)∪⋯∪(𝜏⁺_2𝑚,𝜏⁻_2𝑚)∪(𝜏⁺_2𝑚,∞). Furthermore, system (1.1) undergoes a Hopf bifurcation at the positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂) when 𝜏₂=𝜏^±_2𝑘,𝑘=0,1,2,….

Case c (𝜏₁>0,𝜏₂>0)
We consider (2.6) with 𝜏₂ in its stable interval. Regarding 𝜏₁ as a parameter, without loss of generality, we consider system (1.1) under the assumptions (H1) and (H3). Let 𝑖𝜔(𝜔>0) be a root of (2.6), then we can obtain 𝜔⁴+𝑘₁𝜔³+𝑘₂𝜔²+𝑘₃𝜔+𝑘₄=0, (2.25) where 𝑘₁=2𝑛₃sin𝜔𝜏₂,𝑘₂=2𝑚₂𝑛₂−𝑚₁𝑛₃cos𝜔𝜏₂−𝑚₁𝑛₁+𝑚₁+𝑛₁+𝑛₃cos𝜔𝜏₂²−𝑚²₃,𝑘₃=2𝑚₂𝑛₂−𝑚₁𝑛₃cos𝜔𝜏₂𝑛₃sin𝜔𝜏₂−2𝑚²₃𝑛₃sin𝜔𝜏₂+2𝑚₁+𝑛₁+𝑛₃cos𝜔𝜏₂𝑚₂𝑛₂−𝑚₁𝑛₃sin𝜔𝜏₂,𝑘₄=𝑚₂𝑛₂−𝑚₁𝑛₃cos𝜔𝜏₂−𝑚₁𝑛₁²+𝑚₂𝑛₂−𝑚₁𝑛₃sin𝜔𝜏₂²−𝑚₃𝑛₁−𝑚₃𝑛₃cos𝜔𝜏₂²−𝑚₃𝑛₃sin𝜔𝜏₂². (2.26) Denote 𝐻(𝜔)=𝜔⁴+𝑘₁𝜔³+𝑘₂𝜔²+𝑘₃𝜔+𝑘₄. (2.27) Assume that (H5)𝑚₂𝑛₂−𝑚₁𝑛₃−𝑚₁𝑛₁²<𝑚₃𝑛₁−𝑚₃𝑛₃². (2.28) It is easy to check that 𝐻(0)<0 if (H5) holds and lim_𝜔→+∞𝐻(𝜔)=+∞. We can obtain that (2.25) has finite positive roots 𝜔₁,𝜔₂,…,𝜔_𝑛. For every fixed 𝜔_𝑖,𝑖=1,2,3,…,𝑘, there exists a sequence {𝜏^𝑗_1𝑖∣𝑗=1,2,3,…}, such that (2.25) holds. Let 𝜏₁₀=min𝜏^𝑗_1𝑖∣𝑖=1,2,…,𝑘;𝑗=1,2,…. (2.29) When 𝜏₁=𝜏₁₀, (2.6) has a pair of purely imaginary roots ±𝑖𝜔^∗ for 𝜏₂∈[0,𝜏₂₀).
In the following, we assume that (H6)𝑑(Re𝜆)𝑑𝜏_{1𝜆=𝑖𝜔∗}≠0. (2.30) Thus, by the general Hopf bifurcation theorem for FDEs in Hale [24], we have the following result on the stability and Hopf bifurcation in system (1.1).

Theorem 2.3. For system (1.1), assume that one of the conditions (a), (b), (c), and (d) holds and suppose (H1),(H3), and (H5) are satisfied, and 𝜏₂∈[0,𝜏₂₀), then the positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂) is asymptotically stable when 𝜏₁∈(0,𝜏₁₀), and system (1.1) undergoes a Hopf bifurcation at the positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂) when 𝜏₁=𝜏₁₀.

3. Direction and Stability of the Hopf Bifurcation

In the previous section, we obtained conditions for Hopf bifurcation to occur when 𝜏₁=𝜏₁₀. In this section, we shall derive the explicit formulae determining the direction, stability, and period of these periodic solutions bifurcating from the positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂) at this critical value of 𝜏₁, by using techniques from normal form and center manifold theory [7]. Throughout this section, we always assume that system (1.1) undergoes Hopf bifurcation at the positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂) for 𝜏₁=𝜏₁₀, and then ±𝑖𝜔^∗ is corresponding purely imaginary roots of the characteristic equation at the positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂).

Without loss of generality, we assume that 𝜏^∗₂<𝜏₁₀, where 𝜏^∗₂∈(0,𝜏₂₀). For convenience, let 𝑢_𝑖(𝑡)=𝑢_𝑖(𝜏𝑡)(𝑖=1,2) and 𝜏₁=𝜏₁₀+𝜇, where 𝜏₁₀ is defined by (2.28) and 𝜇∈𝑅, drop the bar for the simplification of notations, then system (1.1) can be written as an FDE in 𝐶=𝐶([−1,0],𝑅²) as𝑢̇(𝑡)=𝐿_𝜇𝑢_𝑡+𝐹𝜇,𝑢_𝑡, (3.1) where 𝑢(𝑡)=(𝑢₁(𝑡),𝑢₂(𝑡))^𝑇∈𝐶 and 𝑢_𝑡(𝜃)=𝑢(𝑡+𝜃)=(𝑢₁(𝑡+𝜃),𝑢₂(t+𝜃))^𝑇∈𝐶, and 𝐿_𝜇∶𝐶→𝑅,𝐹∶𝑅×𝐶→𝑅 are given by𝐿_𝜇𝜙=𝜏₁₀𝐵𝜙+𝜇¹(0)𝜙₂(0)+𝜏₁₀𝐶𝜙+𝜇¹−𝜏^∗₂𝜏₁₀𝜙₁−𝜏^∗₂𝜏₁₀+𝜏₁₀𝐷𝜙+𝜇¹()𝜙−1₂(),𝐹=𝜏−1𝜇,𝜙₁₀+𝜇𝑓₁,𝑓₂^𝑇, (3.2)(3.3) respectively, where 𝜙(𝜃)=(𝜙₁(𝜃),𝜙₂(𝜃))^𝑇∈𝐶,𝑚𝐵=¹𝑚₂0𝑛₁𝑛,𝐶=00₂𝑛₃𝑚,𝐷=³0,𝑓00₁=−𝜙₁(0)𝜙₁(−1)+𝑙₁𝜙²₁(0)+𝑙₂𝜙²₂(0)+𝑙₃𝜙₁(0)𝜙₂(0)+𝑙₄𝜙³₁(0)+𝑙₅𝜙³₂(0)+𝑙₆𝜙²₁(0)𝜙₂(0)+𝑙₇𝜙₁(0)𝜙²₂(0)𝑓+h.o.t.,₂=𝑘₁𝜙₁−𝜏^∗₂𝜏₁₀𝜙₂(0)+𝑘₂𝜙₂(0)𝜙₂−𝜏^∗₂𝜏₁₀+𝑘₃𝜙²₁−𝜏^∗₂𝜏₁₀+𝑘₄𝜙₁−𝜏^∗₂𝜏₁₀𝜙₂−𝜏^∗₂𝜏₁₀+𝑒₁𝜙²₁−𝜏^∗₂𝜏₁₀𝜙₂(0)+𝑒₂𝜙₁−𝜏^∗₂𝜏₁₀𝜑₂−𝜏^∗₂𝜏₁₀𝜙₂(0)+𝑒₃𝜙³₁−𝜏^∗₂𝜏₁₀+𝑒₄𝜙₁−𝜏^∗₂𝜏₁₀𝜙²₂−𝜏^∗₂𝜏₁₀+𝑒₅𝜙³₂−𝜏^∗₂𝜏₁₀+h.o.t. (3.4) From the discussion in Section 2, we know that, if 𝜇=0, then system (3.1) undergoes a Hopf bifurcation at the positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂) and the associated characteristic equation of system (3.1) has a pair of simple imaginary roots ±𝜔^∗𝜏₁₀.

By the representation theorem, there is a matrix function with bounded variation components 𝜂(𝜃,𝜇),𝜃∈[−1,0] such that𝐿_𝜇𝜙=⁰₋₁𝑑𝜂𝜃,𝜇𝜙(𝜃),for𝜙∈𝐶. (3.5) In fact, we can choose𝜂=𝜏𝜃,𝜇¹⁰()𝜏+𝜇𝐵+𝐶+𝐷,𝜃=0,₁₀()−𝜏+𝜇𝐶+𝐷,𝜃∈^∗2𝜏₁₀,𝜏,0₁₀𝜏+𝜇𝐷,𝜃∈−1,−^∗2𝜏₁₀,0,𝜃=−1. (3.6)

For 𝜙∈𝐶([−1,0],𝑅²), define 𝐴𝜇(𝜃)𝜙=𝑑𝜙𝑑𝜃,−1≤𝜃<0,₀₋₁𝜙(𝑠)𝐹𝑑𝜂𝑠,𝜇,𝜃=0,𝑅𝜙=0,−1≤𝜃<0,𝜇,𝜙,𝜃=0. (3.7) Then, (3.1) is equivalent to the abstract differential equation𝑢̇_𝑡=𝐴𝜇𝑢_𝑡+𝑅𝜇𝑢_𝑡, (3.8) where 𝑢_𝑡(𝜃)=𝑢(𝑡+𝜃),𝜃∈[−1,0].

For 𝜓∈𝐶([0,1],(𝑅²)^∗), define𝐴^∗𝜓(𝑠)=−𝑑𝜓(𝑠)(],𝑑𝑠,𝑠∈0,1₀₋₁𝑑𝜂^𝑇(𝑡,0)𝜓(−𝑡),𝑠=0. (3.9) For 𝜙∈𝐶([−1,0],𝑅²) and 𝜓∈𝐶([0,1],(𝑅²)^∗), define the bilinear form𝜓,𝜙=𝜓(0)𝜙(0)−⁰₋₁^𝜃_𝜉=0𝜓^𝑇(𝜉−𝜃)𝑑𝜂(𝜃)𝜙(𝜉)𝑑𝜉, (3.10) where 𝜂(𝜃)=𝜂(𝜃,0), the 𝐴=𝐴(0) and 𝐴^∗ are adjoint operators. By the discussions in Section 2, we know that ±𝑖𝜔^∗𝜏₁₀ are eigenvalues of 𝐴(0), and they are also eigenvalues of 𝐴^∗ corresponding to 𝑖𝜔^∗𝜏₁₀ and −𝑖𝜔^∗𝜏₁₀, respectively. By direct computation, we can obtain𝑞(𝜃)=(1,𝛼)^𝑇𝑒^{𝑖𝜔∗𝜏10𝜃},𝑞^∗(𝑠)=𝑀1,𝛼^∗𝑒^{𝑖𝜔∗𝜏10𝑠}1,𝑀=𝐾, (3.11) where𝛼=𝑖𝜔^∗−𝑚₁−𝑚₃𝑒^{−𝑖𝜔∗𝜏10}𝑚₂,𝛼^∗=−𝑖𝜔^∗+𝑚₁+𝑚₃𝑒^{−𝑖𝜔∗𝜏10}𝑛₂𝑒^{−𝑖𝜔∗𝜏∗2/𝜏10},𝐾=1+𝛼𝛼^∗+𝜏₁₀𝑚₃𝑒^{𝑖𝜔∗𝜏10}+𝜏^∗₂𝜏₁₀𝑛₂𝛼𝑒^{𝑖𝜔∗(𝜏∗2/𝜏10)}+𝛼𝛼^∗𝑛₃𝜏^∗₂𝜏₁₀𝑒^{𝑖𝜔∗(𝜏∗2/𝜏10)}. (3.12) Furthermore, ⟨𝑞^∗(𝑠),𝑞(𝜃)⟩=1 and ⟨𝑞^∗(𝑠),𝑞(𝜃)⟩=0.

Next, we use the same notations as those in Hassard [7] and we first compute the coordinates to describe the center manifold 𝐶₀ at 𝜇=0. Let 𝑢_𝑡 be the solution of (3.1) when 𝜇=0.

Define𝑧(𝑡)=𝑞^∗,𝑢_𝑡,𝑊(𝑡,𝜃)=𝑢_𝑡(𝜃)−2Re𝑧(𝑡)𝑞(𝜃), (3.13) on the center manifold 𝐶₀, and we have𝑊(𝑡,𝜃)=𝑊(𝑧(𝑡),𝑧(𝑡),𝜃), (3.14) where𝑊(𝑧(𝑡),𝑧(𝑡),𝜃)=𝑊(𝑧,𝑧)=𝑊₂₀𝑧²2+𝑊₁₁𝑧𝑧+𝑊₀₂𝑧²2+⋯, (3.15) and 𝑧 and 𝑧 are local coordinates for center manifold 𝐶₀ in the direction of 𝑞^∗ and 𝑞^∗. Noting that 𝑊 is also real if 𝑢_𝑡 is real, we consider only real solutions. For solutions 𝑢_𝑡∈𝐶₀ of (3.1),𝑧̇(𝑡)=𝑖𝜔^∗𝜏₁₀𝑧+𝑞^∗(𝜃)𝐹(0,𝑊(𝑧,𝑧,𝜃))+2Re𝑧𝑞(𝜃)=^def𝑖𝜔^∗𝜏₁₀𝑧+𝑞^∗(0)𝐹₀. (3.16) That is𝑧̇(𝑡)=𝑖𝜔^∗𝜏₁₀𝑧+𝑔(𝑧,𝑧), (3.17) where𝑔(𝑧,𝑧)=𝑔₂₀𝑧²2+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2+⋯. (3.18) Hence, we have to obtain the expression of 𝑔(𝑧,𝑧) (see Appendix B). Then, it is easy to obtain the expression of 𝑔₂₀,𝑔₁₁,𝑔₀₂,𝑔₂₁ (see Appendix B).

For unknown 𝑊^(𝑖)₂₀(𝜃),𝑊^(𝑖)₁₁(𝜃),(𝑖=1,2) in 𝑔₂₁, we still need to compute them.

From (3.8), (3.13), we have𝑊^′=𝐴𝑊−2Re𝑞^∗(0)𝑓𝑞(𝜃),−1≤𝜃<0,𝐴𝑊−2Re𝑞^∗(0)𝑓𝑞(𝜃)+𝐹,𝜃=0,_def=𝐴𝑊+𝐻(𝑧,𝑧,𝜃), (3.19) where𝐻(𝑧,𝑧,𝜃)=𝐻₂₀(𝜃)𝑧²2+𝐻₁₁(𝜃)𝑧𝑧+𝐻₀₂(𝜃)𝑧²2+⋯. (3.20) Comparing the coefficients, we obtain𝐴−2𝑖𝜏₁₀𝜔^∗𝑊₂₀=−𝐻₂₀(𝜃),𝐴𝑊₁₁(𝜃)=−𝐻₁₁(𝜃), (3.21)(3.22) and we know that, for 𝜃∈[−1,0),𝐻(𝑧,𝑧,𝜃)=−𝑞^∗(0)𝑓₀𝑞(𝜃)−𝑞^∗(0)𝑓₀𝑞(𝜃)=−𝑔(𝑧,𝑧)𝑞(𝜃)−𝑔(𝑧,𝑧)𝑞(𝜃). (3.23) Comparing the coefficients of (3.23) with (3.20) gives that𝐻₂₀(𝜃)=−𝑔₂₀𝑞(𝜃)−𝑔₀₂𝑞(𝜃),𝐻₁₁(𝜃)=−𝑔₁₁𝑞(𝜃)−𝑔₁₁𝑞(𝜃). (3.24)(3.25) From (3.21), (3.24), and the definition of 𝐴, we get𝑊̇₂₀(𝜃)=2𝑖𝜔^∗𝜏₁₀𝑊₂₀(𝜃)+𝑔₂₀𝑞(𝜃)+𝑔₀₂𝑞(𝜃). (3.26) Noting that 𝑞(𝜃)=𝑞(0)𝑒^{𝑖𝜔∗𝜏10𝜃}, we have𝑊₂₀(𝜃)=𝑖𝑔₂₀𝜔^∗𝜏₁₀𝑞(0)𝑒^{𝑖𝜔∗𝜏10𝜃}+𝑖𝑔₀₂3𝜔^∗𝜏₁₀𝑞(0)𝑒^{−𝑖𝜔∗𝜏10𝜃}+𝐸₁𝑒^{2𝑖𝜔∗𝜏10𝜃}, (3.27) where 𝐸₁=(𝐸₁⁽¹⁾,𝐸₁⁽²⁾)^𝑇∈𝑅² is a constant vector.

Similarly, from (3.22), (3.25), and the definition of 𝐴, we have𝑊̇₁₁(𝜃)=𝑔₁₁𝑞(𝜃)+𝑔₁₁𝑞(𝜃),𝑊₁₁(𝜃)=−𝑖𝑔₁₁𝜔^∗𝜏₁₀𝑞(0)𝑒^{𝑖𝜔∗𝜏10𝜃}+𝑖𝑔₁₁𝜔^∗𝜏₁₀𝑞(0)𝑒^{−𝑖𝜔∗𝜏10𝜃}+𝐸₂, (3.28)(3.29) where 𝐸₂=(𝐸₂⁽¹⁾,𝐸₂⁽²⁾)^𝑇∈𝑅² is a constant vector

In what follows, we shall seek appropriate 𝐸₁,𝐸₂ in (3.27), (3.29), respectively. It follows from the definition of 𝐴 and (3.24), (3.25) that𝑑𝜂(⁰₋₁𝜃)𝑊₂₀(𝜃)=2𝑖𝜔^∗𝜏₁₀𝑊₂₀(0)−𝐻₂₀(0),₀₋₁𝑑𝜂(𝜃)𝑊₁₁(𝜃)=−𝐻₁₁(0), (3.30)(3.31) where 𝜂(𝜃)=𝜂(0,𝜃).

From (3.21), we have𝐻₂₀(0)=−𝑔₂₀𝑞(0)−𝑔₀₂𝑞(0)+2𝜏₁₀𝐻₁,𝐻₂^𝑇,𝐻₁₁(0)=−𝑔₁₁𝑞(0)−𝑔₁₁(0)𝑞(0)+2𝜏₁₀𝑃₁,𝑃₂^𝑇, (3.32)(3.33) where𝐻₁=−𝑒^{−𝑖𝜔∗𝜏10}+𝑙₁+𝑙₂𝛼²+𝑙₃𝐻𝛼,₂=𝑘₁𝛼+𝑘₂𝛼²𝑒^{−𝑖𝜔∗𝜏∗2}+𝑘₃+𝑘₄𝛼𝑒^{−2𝑖𝜔∗𝜏∗2},𝑃₁1=−2𝑒^{−𝑖𝜔∗𝜏10}+𝑒^{𝑖𝜔∗𝜏10}+𝑙₁+𝑙₂|𝛼|²{𝛼},𝑃+2Re₂=𝑘₁Re𝛼𝑒^{𝑖𝜔∗𝜏∗2}+𝑘₂|𝛼|²𝑒^{−𝑖𝜔∗𝜏∗2}+𝑒^{𝑖𝜔∗𝜏∗2}+𝑘₃+𝑘₄Re{𝛼}. (3.34) Noting that𝑖𝜔^∗𝜏₁₀𝐼−⁰₋₁𝑒^{𝑖𝜔∗𝜏10𝜃}𝑑𝜂(𝜃)𝑞(0)=0,−𝑖𝜔^∗𝜏₁₀𝐼−⁰₋₁𝑒^{−𝑖𝜔∗𝜏10𝜃}𝑑𝜂(𝜃)𝑞(0)=0, (3.35) and substituting (3.27) and (3.32) into (3.30), we have2𝑖𝜔^∗𝜏₁₀𝐼−⁰₋₁𝑒^{2𝑖𝜔∗𝜏10𝜃}𝑑𝜂(𝜃)𝐸₁=2𝜏₁₀𝐻₁,𝐻₂^𝑇. (3.36) That is2𝑖𝜔^∗−𝑚₃𝑒^{−2𝑖𝜔∗𝜏10}−𝑚₂−𝑛₂𝑒^{−2𝑖𝜔∗𝜏∗2}2𝑖𝜔^∗−𝑛₁𝑒^{−2𝑖𝜔∗𝜏10}−𝑛₃𝑒^{−2𝑖𝜔∗𝜏∗2}𝐸₁=2𝐻₁,𝐻₂^𝑇. (3.37) It follows that𝐸₁⁽¹⁾=Δ₁₁Δ₁,𝐸₁⁽²⁾=Δ₁₂Δ₁, (3.38) whereΔ₁=2𝑖𝜔^∗−𝑚₃𝑒^{−2𝑖𝜔∗𝜏10}2𝑖𝜔^∗−𝑛₁𝑒^{−2𝑖𝜔∗𝜏10}−𝑛₃𝑒^{−2𝑖𝜔∗𝜏∗2}−𝑚₂𝑛₂𝑒^{−2𝑖𝜔∗𝜏∗2},Δ₁₁=2𝐻₁2𝑖𝜔^∗−𝑛₁𝑒^{−2𝑖𝜔∗𝜏10}−𝑛₃𝑒^{−2𝑖𝜔∗𝜏∗2}+2𝐻₂𝑚₂,Δ₁₂=2𝐻₂2𝑖𝜔^∗−𝑚₃𝑒^{−2𝑖𝜔∗𝜏10}+2𝐻₁𝑛₂𝑒^{−2𝑖𝜔∗𝜏∗2}. (3.39) Similarly, substituting (3.28) and (3.33) into (3.31), we have𝑑𝜂(⁰₋₁𝜃)𝐸₂=2𝜏₁₀𝑃₁,𝑃₂^𝑇. (3.40) That is𝑚₁+𝑚₃𝑚₂𝑛₂𝑛₁+𝑛₃𝐸₂=2−𝑃₁,−𝑃₂^𝑇. (3.41) It follows that𝐸₂⁽¹⁾=Δ₂₁Δ₂,𝐸₂⁽²⁾=Δ₂₂Δ₂, (3.42) whereΔ₂=𝑚₁+𝑚₃𝑛₁+𝑛₃−𝑚₂𝑛₂,Δ₂₁=2𝑚₂𝑃₂−2𝑃₁𝑛₁+𝑛₃,Δ₂₂=2𝑛₂𝑃₁−2𝑃₂𝑚₁+𝑚₃. (3.43) From (3.27), (3.29), (3.38), (3.42), we can calculate 𝑔₂₁ and derive the following values:𝑐₁(0)=𝑖2𝜔^∗𝜏₁₀𝑔₂₀𝑔₁₁−2𝑔₁₁²−𝑔₀₂²3+𝑔₂₁2,𝜇₂𝑐=−Re₁(0)Re𝜆^′𝜏₁₀,𝛽₂=2Re𝑐₁(0),𝑇₂𝑐=−Im₁(0)+𝜇₂Im𝜆^′𝜏₁₀𝜔^∗𝜏₁₀. (3.44) These formulas give a description of the Hopf bifurcation periodic solutions of (3.1) at 𝜏=𝜏₁₀ on the center manifold. From the discussion above, we have the following result.

Theorem 3.1. The periodic solution is forward (backward) if 𝜇₂>0(𝜇₂<0); the bifurcating periodic solutions are orbitally asymptotically stable with asymptotical phase (unstable) if 𝛽₂<0(𝛽₂>0); the periods of the bifurcating periodic solutions increase (decrease) if 𝑇₂>0(𝑇₂<0).

4. Numerical Examples

In this section, we present some numerical results of system (1.1) to verify the analytical predictions obtained in the previous section. From Section 3, we may determine the direction of a Hopf bifurcation and the stability of the bifurcation periodic solutions. Let us consider the following system:𝑢̇₁(𝑡)=𝑢₁(𝑡)1−𝑢₁𝑡−𝜏₁−𝑢₁(𝑡)𝑢₂(𝑡)0.05+𝑢₁(𝑡)+0.6𝑢₂(𝑡),𝑢̇₂(𝑡)=𝑢₂(𝑡𝑢)−0.5+¹𝑡−𝜏₂0.05+𝑢₁𝑡−𝜏₂+0.6𝑢₂𝑡−𝜏₂, (4.1) which has a positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂)≈(0.3038,0.4230). When 𝜏₁=0, then we can easily obtain that (H1) and (H3) are satisfied. Take 𝑛=0, for example, by some computation by means of Matlab 7.0, we get 𝜔₀≈0.1342,𝜏₂₀≈1.52. From Lemma 2.2, we know that the transversal condition is satisfied. Thus, the positive equilibrium 𝐸₀≈(0.3038,0.4230) is asymptotically stable for 𝜏₂<𝜏₂₀≈1.52 and unstable for 𝜏₂>𝜏₂₀≈1.52 which is shown in Figure 1. When 𝜏₂=𝜏₂₀≈1.52, (4.1) undergoes a Hopf bifurcation at the positive equilibrium 𝐸₀≈(0.3038,0.4230), that is, a small amplitude periodic solution occurs around 𝐸₀≈(0.3038,0.4230) when 𝜏₁=0 and 𝜏₂ is close to 𝜏₂₀=1.52 which is shown in Figure 2.

(a)

(b)

(c)

(d)

Figure 1 

Trajectory portrait and phase portrait of system (4.1) with 𝜏1=0,𝜏2=1.5<𝜏20≈1.52. The positive equilibrium 𝐸0(0.3038,0.4230) is asymptotically stable. The initial value is (0.2,0.2).

(a)

(b)

(c)

(d)

Figure 2 

Trajectory portrait and phase portrait of system (4.1) with 𝜏1=0,𝜏2=1.6>𝜏20≈1.52. Hopf bifurcation occurs from the positive equilibrium 𝐸0(0.3038,0.4230). The initial value is (0.2,0.2).

Let 𝜏₂=1.5∈(0,1.52) and choose 𝜏₁ as a parameter. We have 𝜏₁₀≈1.2930. Then, the positive equilibrium is asymptotically when 𝜏₁∈[0,𝜏₁₀). The Hopf bifurcation value of (4.1) is 𝜏₁₀≈1.2930. By the algorithm derived in Section 3, we can obtain𝜆^′𝜏₁₀≈0.5018−7.2021𝑖,𝑐₁(0)𝜇≈−2.0231−3.3225𝑖,₂≈0.2646,𝛽₂≈−4.2342,𝑇₂≈8.3701. (4.2) Furthermore, it follows that 𝜇₂>0 and 𝛽₂<0. Thus, the positive equilibrium 𝐸₀≈(0.3038,0.4230) is stable when 𝜏₁<𝜏₁₀ as is illustrated by the computer simulations (see Figure 3). When 𝜏₁ passes through the critical value 𝜏₁₀, the positive equilibrium 𝐸₀≈(0.3038,0.4230) loses its stability and a Hopf bifurcation occurs, that is, a family of periodic solutions bifurcations from the positive equilibrium 𝐸₀≈(0.3038,0.4230). Since 𝜇₂>0 and 𝛽₂<0, the direction of the Hopf bifurcation is 𝜏₁>𝜏₁₀, and these bifurcating periodic solutions from 𝐸₀≈(0.3038,0.4230) at 𝜏₁₀ are stable, which are depicted in Figure 4.

(a)

(b)

(c)

(d)

Figure 3 

Trajectory portrait and phase portrait of system (4.1) with 𝜏2=0.5,𝜏1=1.2<𝜏10≈1.2930. The positive equilibrium 𝐸0(0.3038,0.4230) is asymptotically stable. The initial value is (0.2,0.2).

(a)

(b)

(c)

(d)

Figure 4 

Trajectory portrait and phase portrait of system (4.1) with 𝜏2=0.5,𝜏1=1.35>𝜏10≈1.2930. Hopf bifurcation occurs from the positive equilibrium 𝐸∗=(0.32,0.96). The initial value is (0.2,0.2).

5. Conclusions

In this paper, we have investigated local stability of the positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂) and local Hopf bifurcation of a Beddington-DeAngelis functional response predator-prey model with two delays. We have showed that, if one of the conditions (a), (b), (c), and (d) holds and (H1),(H3), and (H5) are satisfied, and 𝜏₂∈[0,𝜏₂₀), then the positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂) is asymptotically stable when 𝜏₁∈(0,𝜏₁₀), as the delay 𝜏₁ increases, the positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂) loses its stability and a sequence of Hopf bifurcations occur at the positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂), that is, a family of periodic orbits bifurcates from the the positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂). At last, the direction of Hopf bifurcation and the stability of the bifurcating periodic orbits are discussed by applying the normal form theory and the center manifold theorem. A numerical example verifying our theoretical results is carried out. In addition, we must point out that, although Gakkhar et al. [4] have also investigated the existence of Hopf bifurcation for system (1.1) with respect to positive equilibrium 𝐸₀(𝑢^∗₁,𝑢^∗₂), it is assumed that 𝜏₁=𝜏₂. For 𝜏₁≠𝜏₂, only numerical simulations are carried out to discuss the existence of Hopf bifurcation. In this paper, under the case 𝜏₁≠𝜏₂, we investigate the existence of Hopf bifurcation quantitatively. Our work generalizes the known results of Sunita Gakkhar et al. [4]. Similarly, we can investigate the Hopf bifurcation of system (1.2) by choosing the delay 𝜏₂ as bifurcation parameter. We will further investigate the topic elsewhere in the near future.

Appendices

A.

The expressions of 𝑚_𝑖,𝑛_𝑖(𝑖=1,2,3), and 𝐹_𝑖(𝑖=1,2) are as follows: 𝑚₁=1−𝑢^∗₁−𝑎₁𝑎+𝑢^∗₁+𝑏𝑢^∗₂𝑢^∗₂−𝑢^∗₁𝑢^∗₂𝑎+𝑢^∗₁+𝑏𝑢^∗₂,𝑚₂𝑎=−¹𝑎+𝑢^∗₁+𝑏𝑢^∗₂𝑢^∗₁−𝑢^∗₁𝑢^∗₂𝑎+𝑢^∗₁+𝑏𝑢^∗₂,𝑚₃=−𝑢^∗₁,𝑛₁=−𝑑𝑑₁+𝑑₁𝑢^∗₁𝑎+𝑢^∗₁+𝑏𝑢^∗₂,𝑛₂=𝑑₁𝑢^∗₁𝑢^∗₂𝑎+𝑢^∗₁+𝑏𝑢^∗₂,𝑛₃𝑑=−¹𝑏𝑢^∗₁𝑢^∗₂𝑎+𝑢^∗₁+𝑏𝑢^∗₂²,𝐹₁=−𝑢₁(𝑡)𝑢₁𝑡−𝜏₁+𝑙₁𝑢²₁(𝑡)+𝑙₂𝑢²₂(𝑡)+𝑙₃𝑢₁(𝑡)𝑢₂(𝑡)+𝑙₄𝑢³₁(𝑡)+𝑙₅𝑢³₂(𝑡)+𝑙₆𝑢²₁(𝑡)𝑢₂(𝑡)+𝑙₇𝑢₁(𝑡)𝑢²₂(𝑡)𝐹+h.o.t.,₂=𝑘₁𝑢₁𝑡−𝜏₂𝑢₂(𝑡)+𝑘₂𝑢₂(𝑡)𝑢₂𝑡−𝜏₂+𝑘₃𝑢²₁𝑡−𝜏₂+𝑘₄𝑢₁𝑡−𝜏₂𝑢₂𝑡−𝜏₂+𝑒₁𝑢²₁𝑡−𝜏₂𝑢₂(𝑡)+𝑒₂𝑢₁𝑡−𝜏₂𝑢₂𝑡−𝜏₂𝑢₂(𝑡)+𝑒₃𝑢³₁𝑡−𝜏₂+𝑒₄𝑢₁𝑡−𝜏₂𝑢²₂𝑡−𝜏₂+𝑒₅𝑢³₂𝑡−𝜏₂+h.o.t., (A.1) where 𝑙₁𝑎=−¹𝑎+𝑢^∗₁+𝑏𝑢^∗₂−𝑢^∗₂𝑎+𝑢^∗₁+𝑏𝑢^∗₂+𝑢^∗₁𝑢^∗₂𝑎+𝑢^∗₁+𝑏𝑢^∗₂²,𝑙₂𝑎=−¹𝑎+𝑢^∗₁+𝑏𝑢^∗₂−𝑏𝑢^∗₁𝑎+𝑢^∗₁+𝑏𝑢^∗₂+𝑏²𝑢^∗₁𝑢^∗₂𝑎+𝑢^∗₁+𝑏𝑢^∗₂²,𝑙₃𝑎=−¹𝑎+𝑢^∗₁+𝑏𝑢^∗₂𝑢1−^∗1+𝑏𝑢^∗₂𝑎+𝑢^∗₁+𝑏𝑢^∗₂+2𝑏𝑢^∗₁𝑢^∗₂𝑎+𝑢^∗₁+𝑏𝑢^∗₂²,𝑙₄𝑎=−¹𝑎+𝑢^∗₁+𝑏𝑢^∗₂𝑢^∗₂𝑎+𝑢^∗₁+𝑏u^∗₂²−𝑢^∗₁𝑢^∗₂𝑎+𝑢^∗₁+𝑏𝑢^∗₂³,𝑙₅𝑎=−¹𝑎+𝑢^∗₁+𝑏𝑢^∗₂𝑏²𝑢^∗₁𝑎+𝑢^∗₁+𝑏𝑢^∗₂²−𝑏³𝑎+𝑢^∗₁+𝑏𝑢^∗₂³,𝑙₆𝑎=−¹𝑎+𝑢^∗₁+𝑏𝑢^∗₂−1𝑎+𝑢^∗₁+𝑏𝑢^∗₂+𝑢^∗₁+2𝑏𝑢^∗₂𝑎+𝑢^∗₁+𝑏𝑢^∗₂²−3𝑏𝑢^∗₁𝑢^∗₂𝑎+𝑢^∗₁+𝑏𝑢^∗₂³,𝑙₇𝑎=−¹𝑎+𝑢^∗₁+𝑏𝑢^∗₂−𝑏𝑎+𝑢^∗₁+𝑏𝑢^∗₂+𝑏²𝑢^∗₂+2𝑏u^∗₁𝑎+𝑢^∗₁+𝑏𝑢^∗₂²−3𝑏²𝑢^∗₁𝑢^∗₂𝑎+𝑢^∗₁+𝑏𝑢^∗₂³,𝑘₁=𝑑₁1𝑎+𝑢^∗₁+𝑏𝑢^∗₂−𝑢^∗₁𝑎+𝑢^∗₁+𝑏𝑢^∗₂²,𝑘₂𝑑=−¹𝑏𝑎+𝑢^∗₁+𝑏𝑢^∗₂²,𝑘₃=𝑑₁𝑢^∗₂−1𝑎+𝑢^∗₁+𝑏𝑢^∗₂²+𝑢^∗₁𝑎+𝑢^∗₁+𝑏𝑢^∗₂³,𝑘₄=𝑑₁𝑢^∗₂−𝑏𝑎+𝑢^∗₁+𝑏𝑢^∗₂²+2𝑏𝑢^∗₁𝑎+𝑢^∗₁+𝑏𝑢^∗₂³,𝑒₁=𝑑₁𝑢^∗₂1𝑎+𝑢^∗₁+𝑏𝑢^∗₂³−𝑢^∗₁𝑎+𝑢^∗₁+𝑏𝑢^∗₂⁴,𝑒₂=𝑑₁−𝑏𝑎+𝑢^∗₁+𝑏𝑢^∗₂²+2𝑏𝑢^∗₁𝑎+𝑢^∗₁+𝑏𝑢^∗₂³,𝑒₃=𝑑₁𝑢^∗₂1𝑎+𝑢^∗₁+𝑏𝑢^∗₂³−𝑢^∗₁𝑎+𝑢^∗₁+𝑏𝑢^∗₂⁴,𝑒₄=𝑑₁𝑢^∗₂𝑏²𝑎+𝑢^∗₁+𝑏𝑢^∗₂³,𝑒₅𝑑=−¹𝑢^∗₁𝑢^∗₂𝑏³𝑎+𝑢^∗₁+𝑏𝑢^∗₂⁴. (A.2)

B.

The expressions of 𝑔(𝑧,𝑧),𝑔₂₀,𝑔₁₁,𝑔₀₂, and 𝑔₂₁ are as follows: 𝑔(𝑧,𝑧)=𝑞^∗(0)𝐹₀(𝑧,𝑧)=𝐹0,𝑢_𝑡=𝑀𝜏₁₀−𝑒^{−𝑖𝜔∗𝜏10}+𝑙₁+𝑙₂𝛼²+𝑙₃𝛼+𝛼^∗×𝑘₁𝛼𝑒^{−𝑖𝜔∗𝜏∗2}+𝑘₂𝛼²𝑒^{−𝑖𝜔∗𝜏∗2}+𝑘₃𝑒^{−2𝑖𝜔∗𝜏∗2}+𝑘₄𝛼𝑒^{−2𝑖𝜔∗𝜏∗2}𝑧²+𝑀𝜏₁₀−𝑒^{𝑖𝜔∗𝜏10}−𝑒^{−𝑖𝜔∗𝜏10}+2𝑙₁+2𝑙₂|𝛼|²+2𝑙₃Re{𝛼}+𝛼^∗×2𝑘₁Re𝛼𝑒^{𝑖𝜔∗𝜏∗2}+𝑘₂|𝛼|²𝑒^{𝑖𝜔∗𝜏∗2}+𝑒^{−𝑖𝜔∗𝜏∗2}+2𝑘₃+2𝑘₄Re{𝛼}𝑧𝑧+𝑀𝜏₁₀−𝑒^{𝑖𝜔∗𝜏10}+𝑙₁+𝑙₂𝛼²+𝑙₃𝛼+𝛼^∗×𝑘₁𝛼𝑒^{𝑖𝜔∗𝜏∗2}+𝑘₂|𝛼|²𝑒^{𝑖𝜔∗𝜏∗2}+𝑘₃𝑒^{2𝑖𝜔∗𝜏∗2}+𝑘₄𝛼𝑒^{𝑖𝜔∗𝜏∗2}𝑧²+𝑀𝜏₁₀−12𝑊⁽¹⁾₂₀(0)𝑒^{𝑖𝜔∗𝜏10}+12𝑊⁽¹⁾₂₀(0)+𝑊⁽¹⁾₁₁(0)𝑒^{−𝑖𝜔∗𝜏10}+𝑊⁽¹⁾₁₁(0)+𝑙₁𝑊⁽¹⁾₂₀(0)+2𝑊⁽¹⁾₁₁(0)+𝑙₂𝛼𝑊⁽²⁾₂₀(0)+2𝛼𝑊⁽²⁾₁₁(0)+𝑙₃𝑊⁽²⁾₁₁(0)+12𝑊⁽²⁾₂₀(0)+12𝛼𝑊⁽¹⁾₂₀(0)+𝛼𝑊⁽¹⁾₁₁(0)+3𝑙₄+3𝑙₅𝛼²𝛼+𝑙₆(2𝛼+𝛼)+𝑙₇2|𝛼|²+𝛼²+𝛼^∗×𝑘₁𝑊⁽²⁾₁₁(0)𝑒^{−𝑖𝜔∗𝜏∗2}+12𝛼𝑊⁽¹⁾₂₀−𝜏^∗₂𝜏₁₀+12𝑊⁽²⁾₂₀(0)𝑒^{𝑖𝜔∗𝜏∗2}+𝑊⁽²⁾₁₁(0)𝑒^{−𝑖𝜔∗𝜏∗2}+𝑘₂𝛼𝑊⁽²⁾₁₁−𝜏^∗₂𝜏₁₀+12𝛼𝑊⁽²⁾₂₀−𝜏^∗₂𝜏₁₀+12𝛼𝑊⁽²⁾₂₀(0)𝑒^{𝑖𝜔∗𝜏∗2}+𝛼𝑊⁽²⁾₁₁(0)𝑒^{−𝑖𝜔∗𝜏∗2}+𝑘₃2𝑊⁽¹⁾₁₁−𝜏^∗₂𝜏₁₀𝑒^{−𝑖𝜔∗𝜏∗2}+𝑊⁽¹⁾₂₀−𝜏^∗₂𝜏₁₀𝑒^{𝑖𝜔∗𝜏∗2}+𝑘₄𝑊⁽²⁾₁₁−𝜏^∗₂𝜏₁₀𝑒^{−𝑖𝜔∗𝜏∗2}+12𝑊⁽²⁾₂₀−𝜏^∗₂𝜏₁₀𝑒^{𝑖𝜔∗𝜏∗2}+12𝛼𝑊⁽¹⁾₂₀−𝜏^∗₂𝜏₁₀𝑒^{𝑖𝜔∗𝜏∗2}+𝛼𝑊⁽¹⁾₁₁−𝜏^∗₂𝜏₁₀𝑒^{−𝑖𝜔∗𝜏∗2}+𝑒₁2𝛼+𝛼𝑒^{−2𝑖𝜔∗𝜏∗2}+𝑒₂𝛼²+|𝛼|²+|𝛼|²𝑒^{−2𝑖𝜔∗𝜏∗2}+3𝑒₃𝑒^{−𝑖𝜔∗𝜏∗2}+𝑒₄𝛼²𝑒^{−𝑖𝜔∗𝜏∗2}+2|𝛼|²𝑒^{−𝑖𝜔∗𝜏∗2}+𝑒₅𝛼²𝑒^{−𝑖𝜔∗𝜏∗2}+2|𝛼|²𝛼𝑒^{−𝑖𝜔∗𝜏∗2}𝑧²𝑔𝑧+h.o.t.,₂₀=2𝑀𝜏₁₀−𝑒^{−𝑖𝜔∗𝜏10}+𝑙₁+𝑙₂𝛼²+𝑙₃𝛼+𝛼^∗×𝑘₁𝛼𝑒^{−𝑖𝜔∗𝜏∗2}+𝑘₂𝛼²𝑒^{−𝑖𝜔∗𝜏∗2}+𝑘₃𝑒^{−2𝑖𝜔∗𝜏∗2}+𝑘₄𝛼𝑒^{−2𝑖𝜔∗𝜏∗2},𝑔₁₁=𝑀𝜏₁₀−𝑒^{𝑖𝜔∗𝜏10}−𝑒^{−𝑖𝜔∗𝜏10}+2𝑙₁+2𝑙₂|𝛼|²+2𝑙₃{𝛼}+Re𝛼^∗2𝑘₁Re𝛼𝑒^{𝑖𝜔∗𝜏∗2}+𝑘₂|𝛼|²𝑒^{𝑖𝜔∗𝜏∗2}+𝑒^{−𝑖𝜔∗𝜏∗2}+2𝑘₃+2𝑘₄{𝛼,𝑔Re}₀₂=2𝑀𝜏₁₀−𝑒^{𝑖𝜔∗𝜏10}+𝑙₁+𝑙₂𝛼²+𝑙₃𝛼+𝛼^∗𝑘₁𝛼𝑒^{𝑖𝜔∗𝜏∗2}+𝑘₂|𝛼|²𝑒^{𝑖𝜔∗𝜏∗2}+𝑘₃𝑒^{2𝑖𝜔∗𝜏∗2}+𝑘₄𝛼𝑒^{𝑖𝜔∗𝜏∗2},𝑔₂₁=2𝑀𝜏₁₀−12𝑊⁽¹⁾₂₀(0)𝑒^{𝑖𝜔∗𝜏10}+12𝑊⁽¹⁾₂₀(0)+𝑊⁽¹⁾₁₁(0)𝑒^{−𝑖𝜔∗𝜏10}+𝑊⁽¹⁾₁₁(0)+𝑙₁𝑊⁽¹⁾₂₀(0)+2𝑊⁽¹⁾₁₁(0)+𝑙₂𝛼𝑊⁽²⁾₂₀(0)+2𝛼𝑊⁽²⁾₁₁(0)+𝑙₃𝑊⁽²⁾₁₁(0)+12𝑊⁽²⁾₂₀(0)+12𝛼𝑊⁽¹⁾₂₀(0)+𝛼𝑊⁽¹⁾₁₁(0)+3𝑙₄+3𝑙₅𝛼²𝛼+𝑙₆(2𝛼+𝛼)+𝑙₇2|𝛼|²+𝛼²+𝛼^∗𝑘₁𝑊⁽²⁾₁₁(0)𝑒^{−𝑖𝜔∗𝜏∗2}+12𝛼𝑊⁽¹⁾₂₀−𝜏^∗₂𝜏₁₀+12𝑊⁽²⁾₂₀(0)𝑒^{𝑖𝜔∗𝜏∗2}+𝑊⁽²⁾₁₁(0)𝑒^{−𝑖𝜔∗𝜏∗2}+𝑘₂𝛼𝑊⁽²⁾₁₁−𝜏^∗₂𝜏₁₀+12𝛼𝑊⁽²⁾₂₀−𝜏^∗₂𝜏₁₀+12𝛼𝑊⁽²⁾₂₀(0)𝑒^{𝑖𝜔∗𝜏∗2}+𝛼𝑊⁽²⁾₁₁(0)𝑒^{−𝑖𝜔∗𝜏∗2}+𝑘₃2𝑊⁽¹⁾₁₁−𝜏^∗₂𝜏₁₀𝑒^{−𝑖𝜔∗𝜏∗2}+𝑊⁽¹⁾₂₀−𝜏^∗₂𝜏₁₀𝑒^{𝑖𝜔∗𝜏∗2}+𝑘₄𝑊⁽²⁾₁₁−𝜏^∗₂𝜏₁₀𝑒^{−𝑖𝜔∗𝜏∗2}+12𝑊⁽²⁾₂₀−𝜏^∗₂𝜏₁₀𝑒^{𝑖𝜔∗𝜏∗2}+12𝛼𝑊⁽¹⁾₂₀−𝜏^∗₂𝜏₁₀𝑒^{𝑖𝜔∗𝜏∗2}+𝛼𝑊⁽¹⁾₁₁−𝜏^∗₂𝜏₁₀𝑒^{−𝑖𝜔∗𝜏∗2}+𝑒₁2𝛼+𝛼𝑒^{−2𝑖𝜔∗𝜏∗2}+𝑒₂𝛼²+|𝛼|²+|𝛼|²𝑒^{−2𝑖𝜔∗𝜏∗2}+3𝑒₃𝑒^{−𝑖𝜔∗𝜏∗2}+𝑒₄𝛼²𝑒^{−𝑖𝜔∗𝜏∗2}+2|𝛼|²𝑒^{−𝑖𝜔∗𝜏∗2}+𝑒₅𝛼²𝑒^{−𝑖𝜔∗𝜏∗2}+2|𝛼|²𝛼𝑒^{−𝑖𝜔∗𝜏∗2}. (B.1)

Acknowledgment

This work is supported by National Natural Science Foundation of China (no. 10771215), Doctoral Foundation of Guizhou College of Finance and Economics (2010), the soft Science and Technology Program of Guizhou Province (no. 2011LKC2030), and Nature and Technology Foundation of Guizhou Province (J[2012]2100)