Dynamic Properties of a Forest Fire Model

Abstract

The reaction-diffusion equations have been widely used in physics, chemistry, and other areas. Forest fire can also be described by such equations. We here propose a fighting forest fire model. By using the normal form approach theory and center manifold theory, we analyze the stability of the trivial solution and Hopf bifurcation of this model. Finally, we give the numerical simulations to illustrate the effectiveness of our results.

1. Introduction

The forest fire is an important issue in the world. It has brought us huge losses. It not only burns our forests but also destroys the local ecological environment. Many factors lead to forest fires. Several authors have studied them in depth [1–6]. Some important organizations, especially the USDA Forest Service, have also researched them in their themes [7].

Reaction-diffusion equations have been applied in forest fire model for several years. Some authors analyzed the dynamical behavior of the fire front propagations using hyperbolic reaction-diffusion equations [8]. Lots of articles related to percolation theory [9] and self-organized criticality [10] are trying to provide a different dynamical model for the spread of the fire. In this paper, the model describes the condition that people are putting out the fire when the fire is spreading. We analyze dynamic properties of the reaction-diffusion equations. Kolmogorov et al. proposed the famous KPP model [11] in the 1930s. From then on, it had been applied in various fields including forest fire: 𝑢_𝑡=𝑑₁𝑢_𝑥𝑥+𝑢+𝑓(𝑢),𝑥∈𝑅,𝑡≥0, (1.1) where 𝑢=𝑢(𝑥,𝑡) can be seen as the area of the burned forest. 𝑢_𝑥𝑥 is a diffusion term of 𝑢 in space, and 𝑑₁ is the diffusion coefficient. 𝑓(𝑢) is a nonlinear function. The equation can describe the speed of fire spreading. Zeldovich et al. gave the famous theory of combustion and explosions [12]. We can get inspiration from it:

The people will go to put out the fire as soon as they realize the forest fire. We can use a reaction-diffusion equation to describe it. 𝑣_𝑡=𝑑₂𝑣_𝑥𝑥−𝑐𝑣+𝑔(𝑣). (1.2) In this equation, 𝑣=𝑣(𝑥,𝑡) is the area where the fire has been put out. 𝑣_𝑥𝑥 is a diffusion term of 𝑣 in space, and 𝑑₂ is the diffusion coefficient. 𝑐 is the resurgence probability of 𝑣. 𝑔(𝑣) is a nonlinear function which represents the ability of people to put out the fire.

Now, let us consider the two reaction-diffusion equations together. As we know, 𝑢 and 𝑣 influence each other. Thus, 𝑓 and 𝑔 must be functions of 𝑢,𝑣. We define 𝑔(𝑢,𝑣) by referring to the combustion model [13]: 𝑔(𝑢,𝑣)=𝑢𝑣𝑏(𝑢+1),𝑏>0. (1.3)

Since 𝑔(𝑢,𝑣) has opposite effect on the fire area (or 𝑢), we can also define 𝑓(𝑢,𝑣) by taking into account KPP model [8]: 𝑓(𝑢,𝑣)=−𝑎𝑢²−𝑢𝑣𝑏().𝑢+1 (1.4)

Then we get a new model: 𝑢_𝑡−𝑑₁𝑢_𝑥𝑥=𝑢−𝑎𝑢²−𝑢𝑣𝑏(),𝑣𝑢+1_𝑡−𝑑₂𝑣_𝑥𝑥=−𝑐𝑣+𝑢𝑣𝑏(),𝑢𝑢+1_𝑥(0,𝑡)=𝑣_𝑥(0,𝑡)=0,𝑢_𝑥(𝑙𝜋,𝑡)=𝑣_𝑥()𝑢()()().𝑙𝜋,𝑡=0,𝑥,0≥0,𝑣𝑥,0≥0,𝑥∈0,𝑙𝜋 (1.5) Define 𝐻=𝐿²[0,𝜋]×𝐿²[]=𝑓𝑔0,𝜋∶𝑓,𝑔∈𝐿²[0,𝜋], (1.6) and an inner product is given by 𝑓₁𝑔₁,𝑓₂𝑔₂=𝑓₁,𝑓_2𝐿2+𝑔₁,𝑔_2𝐿2=2𝜋^𝜋₀𝑓₁𝑓₂2𝑑𝑥+𝜋^𝜋₀𝑔₁𝑔₂𝑑𝑥, (1.7) where (𝑓₁,𝑔₁)^𝑇∈𝐻,(𝑓₂,𝑔₂)^𝑇∈𝐻. From the standpoint of biology, we are only interested in the dynamics of model (1.5) in the region: 𝑅²₊={(𝑢,𝑣)∣𝑢>0,𝑣>0}. (1.8)

2. Stability Analysis

Firstly, we consider the location stability [14, 15] and the number of the equilibria of model (1.5) in 𝑅²₊. We can also study autowave solutions [16] of the model. The interior equilibrium point is a root of the following equation: 𝑢−𝑎𝑢²−𝑢𝑣𝑏(𝑢+1)=0,−𝑐𝑣+𝑢𝑣𝑏(𝑢+1)=0. (2.1) It is obvious that (2.1) has an only real solution 𝑌₀=(𝑢₀,𝑣₀), where 𝑢₀=𝑏𝑐1−𝑏𝑐,𝑣₀1=𝑎𝑏𝑎−𝑢₀1+𝑢₀, (2.2) and 𝑏<1/𝑐(𝑎+1).

Now, we analyze the asymptotic stability of (𝑢₀,𝑣₀) by Lyapunov function.

Lemma 2.1. For the model (1.5),
(1) if 𝑎≥1,  (𝑢₀,𝑣₀) is global asymptotic stability.
(2) if 𝑎<1 and (1−𝑎)/𝑐≤𝑏≤1/(𝑎𝑐+𝑐),(𝑢₀,𝑣₀) also has global asymptotic stability.

Proof. Defining 𝜔(𝑢,𝑣)=₀^𝑙𝜋𝑢_𝑢0𝑟/𝑏(𝑟+1)−𝑐()1𝑟/𝑟+1𝑑𝑟𝑑𝑥+𝑏₀^𝑙𝜋𝑣_𝑣0𝑠−𝑣₀𝑠𝑑𝑠𝑑𝑥, (2.3) we can get 𝜕𝜔=1𝜕𝑡𝑏₀^𝑙𝜋ℎ(𝑢)−ℎ𝑢₀𝑝(𝑢)−𝑝𝑢₀𝑑𝑥+𝑌(𝑡), (2.4) whereℎ(𝑢)=𝑢(𝑢)=(),𝑌(𝑡)𝑢+1,𝑝1−𝑎𝑢)(𝑢+1=−𝑑₁𝑐₀^𝑙𝜋ℎ^′(𝑢)ℎ²(𝑢)𝑢²_𝑥𝑑𝑑𝑥+²𝑣₀𝑏₀^𝑙𝜋𝑣²_𝑥𝑣²𝑑𝑥. (2.5) In what follows, we split it into two cases to prove. If 𝑎≥1, for all 𝑢>0,𝑝^′(𝑢)<0, since ℎ^′(𝑢)=1/(𝑢+1)²>0, we can get ℎ(𝑢)−ℎ𝑢₀𝑝(𝑢)−𝑝𝑢₀≤0. (2.6) If 𝑎<1 and (1−𝑎)/𝑐≤𝑏≤1/(𝑎𝑐+𝑐) (equal to 𝑣₀≤𝑏), we can still get (2.6). That is to say 𝑤_𝑡(𝑢,𝑣)<0. (2.7) We prove the conclusion.

Because of the conclusion of Lemma 2.1, we always assume 𝑎<1 and 0<𝑏<(1−𝑎)/(2𝑎𝑐−𝑐). Introducing perturbations 𝑢^∗=𝑢−𝑢₀,𝑣^∗=𝑣−𝑣₀, and replace (𝑢^∗,𝑣^∗) with (𝑢,𝑣), for which model (1.5) yields 𝑢_𝑡−𝑑₁𝑢_𝑥𝑥=𝑢+𝑢₀−𝑎𝑢+𝑢₀²−𝑢+𝑢₀𝑣+𝑣₀𝑏𝑢+𝑢₀,𝑣+1_𝑡−𝑑₂𝑣_𝑥𝑥=−𝑐𝑣+𝑣₀+𝑢+𝑢₀𝑣+𝑣₀𝑏𝑢+𝑢₀,𝑢+1_𝑥(0,𝑡)=𝑣_𝑥(0,𝑡)=0,𝑢_𝑥(𝑙𝜋,𝑡)=𝑣_𝑥()𝑢()()().𝑙𝜋,𝑡=0,𝑥,0≥0,𝑣𝑥,0≥0,𝑥∈0,𝑙𝜋 (2.8) Now we can get the linearized system of parametric model (2.8) at (0,0), 𝑢_𝑡𝑣_𝑡=(𝑏)𝑢𝑣Δ+𝐿, (2.9) where Δ=0𝑑₁𝜕²𝜕𝑥²𝑑₂𝜕²𝜕𝑥²0(𝑏)=𝑢,𝐿⁰1−𝑎−𝑎𝑢₀1+𝑢₀−𝑐1−𝑎𝑢₀1+𝑢₀0. (2.10) The eigenvalues of Δ are as follows: −𝑑₂𝑛²𝑙²,−𝑑₁𝑛²𝑙^2+∞_𝑛=0, (2.11) and the corresponding eigenvectors as follows: 𝛽¹_𝑛,𝛽²_𝑛𝑘^+∞_𝑛=0, (2.12) where 𝛽¹_𝑛=0𝑛cos𝑙𝑥,𝛽²_𝑛=𝑛cos𝑙𝑥0. (2.13) Define for all 𝑦∈𝐻𝑦=^𝑛_𝑘=1𝑌^𝑇_𝑘𝛽¹_𝑘𝛽²_𝑘,𝑌_𝑘=⟨𝑦,𝛽¹_𝑘⟩⟨𝑦,𝛽²_𝑘⟩. (2.14) It is easy to get, 𝜆∈(Δ+𝐿(𝑏)), if and only if the equation ^𝑛_𝑘=1𝑌^𝑇_𝑘𝐸𝜆−𝐿(𝑏)−−𝑑₁𝑘²𝑙²00−𝑑₂𝑘²𝑙²𝛽¹_𝑘𝛽²_𝑘=0 (2.15) is held.

We obtain 𝐸𝜆−𝐿(𝑏)−−𝑑₁𝑛²𝑙²00−𝑑₂𝑛²𝑙²=0. (2.16) Rewrite it as 𝜆²−𝑇_𝑛(𝑏)𝜆+𝐷_𝑛(𝑏)=0, (2.17) where 𝑇_𝑛(𝑏)=𝑏𝑐1−𝑎−2𝑎𝑏𝑐−𝑑1−𝑏𝑐¹+𝑑₂𝑛²𝑙²,𝐷_𝑛(𝑏)=𝑐(1−𝑏𝑐−𝑎𝑏𝑐)−𝑑₂𝑏𝑐1−𝑎−2𝑎𝑏𝑐𝑛1−𝑏𝑐𝑙²²+𝑑₁𝑑₂𝑛⁴𝑙⁴. (2.18) From (2.18), when (1−𝑎)/(𝑎𝑐+𝑐)<𝑏<(1−𝑎)/𝑐 is held, we can get 𝑇_𝑛(𝑏)<0,𝐷_𝑛(𝑏)>0. So the system's eigenvalues have negative real part, and (𝑢₀,𝑣₀) has local asymptotic stability. Then, we can conclude that the system has Hopf bifurcation [14] in 𝑏∈(0,(1−𝑎)/(𝑎𝑐+𝑐)).

Define 𝐵=𝑏₀∣𝑇_𝑛𝑏₀=0,𝐷_𝑛𝑏₀>0,𝑇_𝑗𝑏₀≠0,𝐷_𝑗𝑏₀≠0,∀𝑗≠𝑛. (2.19)𝑏₀∈𝐵∪(0,(1−𝑎)/(𝑎𝑐+𝑐)), 𝛼(𝑏)±𝑖𝜔(𝑏) are characteristic roots of Δ+𝐿(𝑏), where 𝛼(𝑏)=12𝐴(𝑏)−𝑑₁+𝑑₂𝑙²𝑛²,𝜔(𝑏)=𝐷_𝑛(𝑏)−𝛼²(𝑏),𝐴(𝑏)=𝑏𝑐1−𝑎−2𝑎𝑏𝑐.1−𝑏𝑐 (2.20) Now we compute transversality condition: 𝛼^′𝑏₀=12𝑎1−𝑏₀𝑐²1𝑎−1−4𝑏₀𝑐1−𝑏₀𝑐−𝑏₀𝑐1−𝑏₀𝑐²×>0,0<𝑏₀<1𝑐1−2𝑎,11+𝑎<0,𝑐1−2𝑎1+𝑎<𝑏₀<1−𝑎.𝑎𝑐+𝑐 (2.21) Now we consider 𝐴(0)=𝐴(𝑏^𝐵₀)=0 and 𝐴(𝑏) is positive in (0,𝑏^𝐵₀). So we can get the maximum value of 𝐴(𝑏) (defined as 𝐴(𝑏^∗)):

Define 𝑙_𝑛=𝑛𝑑₁+𝑑₂𝐴𝑏^∗,𝑛∈ℕ,𝐴𝑏^∗=𝑎1𝑎+1−√2², (2.22) for all 𝑙∈(𝑙_𝑛,𝑙_𝑛+1],0≤𝑗≤𝑛,𝑏^𝐵_𝑗,− and 𝑏^𝐵_𝑗,+ are two roots of the equation 𝐴(𝑏)=𝑑₁+𝑑₂𝑙²𝑗². (2.23) It is easy to get 0<𝑏^𝐵_1,−<⋯<𝑏^𝐵_𝑛,−<1𝑐1−2𝑎1+𝑎<𝑏^𝐵_𝑛,+<⋯<𝑏^𝐵_1,+<𝑏^𝐵₀. (2.24) Then we give the condition of 𝐷_𝑛(𝑏^𝐵_𝑗,±)≠0 especially 𝐷_𝑛(𝑏^𝐵_𝑗,±)>0.

As we know 𝐷_𝑛(𝑏)≥𝑎𝑐−𝑑₂𝐴𝑏^∗𝑛²𝑙²+𝑑₁𝑑₂𝑛⁴𝑙⁴. (2.25) Then 𝐷_𝑛(𝑏)>0 is held if and only if 𝑑₁𝑑₂𝑑>0,₂𝐴𝑏^∗2−4𝑑₁𝑑₂𝑎𝑐>0. (2.26)

Theorem 2.2. Assume that 𝑑₁,𝑑₂,𝑐>0,0<𝑎<1, and the equation is held: 𝑑₂𝐴²𝑏^∗−4𝑑₁𝑎𝑐>0, (2.27) where 𝑙_𝑛=𝑛𝑑₁+𝑑₂𝐴𝑏^∗,𝑛∈ℕ,𝐴𝑏^∗=𝑎1𝑎+1−√2², (2.28) then for all 𝑙∈(𝑙_𝑛,𝑙_𝑛+1], existing 𝑏=𝑏^𝐵_𝑗,± or 𝑏=𝑏^𝐵₀; there are Hopf bifurcations at the real solution of model (1.5).
Furthermore 0<𝑏^𝐵_1,−<⋯<𝑏^𝐵_𝑛,−<1𝑐1−2𝑎1+𝑎<𝑏^𝐵_𝑛,+<⋯<𝑏^𝐵_1,+<𝑏^𝐵₀. (2.29)

3. Hopf Bifurcation

In the above section, we have already obtained the conditions which ensure that model (2.8) undergoes the Hopf bifurcation at the critical values 𝑏₀ or 𝑏_𝑗,±(𝑗=1,⋯). In the following part, we will study the direction and stability of the Hopf bifurcation based on the normal form approach theory and center manifold theory introduced by Hassard at al. [14].

Firstly, by the transformation 𝑢^∗=𝑢−𝑢₀,𝑣^∗=𝑣−𝑣₀, and replacing (𝑢^∗,𝑣^∗) with (𝑢,𝑣), the parametric system (1.5) is equivalent to the following functional differential equation (FDE) system: 𝜕𝑈=𝑏𝜕𝑡Δ+𝐿₀+𝐹𝑏₀,𝑈, (3.1) where 𝑈=(𝑢,𝑣)^𝑇,𝐹𝑏₀,𝑈=𝑓−𝑓_𝑢𝑢−𝑓_𝑣𝑣𝑔−𝑔_𝑢𝑢−𝑔_𝑣𝑣. (3.2) The adjoint operator of 𝐿_𝑛(𝑏) is defined as 𝐿^∗𝑏₀=𝑑₁𝜕²𝜕𝑥²+𝑓_𝑢𝑏₀𝑓_𝑣𝑏₀𝑔_𝑢𝑏₀𝑑₂𝜕²𝜕𝑥²+𝑔_𝑣𝑏₀. (3.3) It is easy to get 𝑢,𝐿𝑏₀𝑣=𝐿^∗𝑏₀𝑢,𝑣. (3.4) From the discussions in Section 2, define 𝑞^∗=(𝑎^∗_𝑛,𝑏^∗_𝑛)^𝑇cos(𝑛/𝑙)𝑥. We have 𝐿^∗𝑏₀𝑞^∗=−𝑖𝜔𝑞^∗,𝑞^∗,𝑞=1,𝑞^∗,𝑞=0. (3.5) Decompose 𝑋 as 𝑋=𝑋^𝑐⊕𝑋^𝑠, where 𝑋^𝑐={𝑧𝑞+𝑧𝑞∣𝑧∈𝐶} and 𝑋^𝑠={𝑢∈𝑥∣⟨𝑞^∗,𝑢⟩=0}. For all (𝑢,𝑣)∈𝑋, existing 𝑧∈ℂ and 𝜔=(𝜔₁,𝜔₂)∈𝑋^𝑠,  we can obtain 𝑢𝑣=𝑧𝑞+𝜔𝑧𝑞+¹𝜔₂. (3.6) Rewrite (3.1) as 𝑧̇=𝑖𝜔₀𝑧+𝑞^∗,𝐹^∗₀,𝜔̇=𝐿𝑏₀𝜔+𝐻(𝑧,𝑧,𝜔), (3.7) where 𝐹^∗₀=𝑧𝑞+𝑧𝑞+𝜔,𝐻(𝑧,𝑧,𝜔)=𝐹^∗₀−𝑞^∗,𝐹^∗₀𝑞−𝑞^∗,𝐹^∗₀𝑞. (3.8) Using the same notations as in [11], 𝐹^∗₀(𝑈)=12𝑄()+1𝑈,𝑈6𝐶(𝑈,𝑈,𝑈)+𝑂𝑈⁴, (3.9) where 𝑈=(𝑢,𝑣) and 𝑄,𝐶 are symmetrical multilinear functions. We can compute 𝑄=𝐴𝑞,𝑞^1𝑛𝐴²_𝑛cos²𝑛𝑙𝑥,𝑄𝑞,𝑞=𝐵¹_𝑛𝐵²_𝑛cos²𝑛𝑙𝑥,𝐶𝑞,𝑞,𝑞=𝐶¹_𝑛𝐶²_𝑛cos³𝑛𝑙𝑥, (3.10) where 𝐴¹_𝑛=𝑓_𝑢𝑢𝑎²_𝑛+2𝑓_𝑢𝑣𝑎_𝑛𝑏_𝑛+𝑓_𝑣𝑣𝑏²_𝑛,𝐴²_𝑛=𝑔_𝑢𝑢𝑎²_𝑛+2𝑔_𝑢𝑣𝑎_𝑛𝑏_𝑛+𝑔_𝑣𝑣𝑏²_𝑛,𝐵¹_𝑛=𝑓_𝑢𝑢𝑎_𝑛²+𝑓_𝑢𝑣𝑎_𝑛𝑏_𝑛+𝑎_𝑛𝑏_𝑛+𝑓_𝑣𝑣𝑏_𝑛²,𝐵²_𝑛=𝑔_𝑢𝑢𝑎_𝑛²+𝑔_𝑢𝑣𝑎_𝑛𝑏_𝑛+𝑎_𝑛𝑏_𝑛+𝑔_𝑣𝑣𝑏_𝑛²,𝐶¹_𝑛=𝑓_𝑢𝑢𝑢𝑎_𝑛²𝑎_𝑛+𝑓_𝑢𝑢𝑣2𝑎_𝑛²𝑏_𝑛+𝑎²_𝑛𝑏_𝑛+𝑓_𝑢𝑣𝑣2𝑏²_𝑛𝑎_𝑛+𝑏²_𝑛𝑎_𝑛+𝑓_𝑣𝑣𝑣𝑏_𝑛²𝑏_𝑛,𝐶²_𝑛=𝑔_𝑢𝑢𝑢𝑎_𝑛²𝑎_𝑛+𝑔_𝑢𝑢𝑣2𝑎_𝑛²𝑏_𝑛+𝑎²_𝑛𝑏_𝑛+𝑔_𝑢𝑣𝑣2𝑏²_𝑛𝑎_𝑛+𝑏²_𝑛𝑎_𝑛+𝑔_𝑣𝑣𝑣𝑏_𝑛²𝑏_𝑛. (3.11) Define 𝐻(𝑧,)=1𝑧,𝜔2𝐻₂₀𝑧²+𝐻₁₁𝑧1𝑧+2𝐻₀₂𝑧²+𝑜(|𝑧||𝜔|), (3.12) where 𝐻₂₀=𝑄𝑞,𝑞−𝑞^∗,𝑄𝑞,𝑞𝑞−𝑞^∗,𝑄𝑞,𝑞𝐻𝑞,₁₁=𝑄𝑞,𝑞−𝑞^∗,𝑄𝑞,𝑞𝑞−𝑞^∗,𝑄𝑞,𝑞𝑞. (3.13) On the center manifold, we have 1𝜔=2𝜔₂₀𝑧²+𝜔₁₁𝑧1𝑧+2𝜔₀₂𝑧²+𝑜|𝑧|³. (3.14) We can obtain 𝑤₂₀=2𝑖𝜔₀𝐼−𝐿_𝑛𝑏₀⁻¹𝐻₂₀,𝜔₁₁=𝐿_𝑛𝑏₀⁻¹𝐻₁₁. (3.15) Comparing (3.9) and (3.13), we can get 𝐻₂₀=𝑄−𝐴𝑞,𝑞^1𝑛𝐴²_𝑛cos²𝑛𝑙𝐴𝑥−^1𝑛𝐴²_𝑛12cos²𝑛21𝑥+2,𝑛∈ℕ^∗,𝐴¹₀𝐴²₀−𝑞^∗𝑎,𝑄𝑞,𝑞⁰𝑏₀−𝑞^∗,𝑄𝑞,𝑞𝑎₀𝑏₀,𝑛=0. (3.16) Similarly 𝜔₁₁=−12𝐿𝑏_𝑜⁻¹𝐶¹_𝑛𝐶²_𝑛cos²𝑛2+1,𝑛∈ℕ^∗,−𝐿𝑏₀⁻¹𝑐¹₀𝑐²₀−𝑞^∗,𝑄𝑞,𝑞𝑎¹₀𝑏²₀−𝑞^∗,𝑄𝑞,𝑞𝑎₀𝑏₀,𝑛=0. (3.17) Then on the center manifold rewrite (𝑑𝑈) as 𝑑𝑧𝑑𝑡=𝑖𝜔₀𝑧+𝑞^∗,𝐹^∗₀=𝑖𝜔₀𝑧+_{2≤𝑖+𝑗≤3}𝑔_𝑖𝑗𝑧𝑖!𝑗!^𝑖𝑧^𝑗+𝑜|𝑧|⁴, (3.18) where 𝑔₂₀=𝑞^∗,𝑄𝑞,𝑞=4𝑐𝜔₀𝑎²−𝑎(1−𝑎)²𝜔₀−2𝑐𝑎²(3𝑎−1)𝑖1−𝑎²𝜔₀,𝑔₁₁=𝑞^∗,𝑄𝑞,𝑞=𝑎𝑎−𝑞𝜔₀−2𝑐𝑎²𝑖(1+𝑎)𝜔₀,𝑔₀₂=𝑞^∗,𝑄𝑞,𝑞=𝑎(1−𝑎)²+2𝑐𝜔₀𝑎²−4𝑐𝑎²𝑖1−𝑎²𝜔₀,𝑔₂₁=2𝑞^∗,𝑄𝑤₁₁,𝑞+𝑞^∗,𝑄𝑤₂₀,𝑞+𝑞^∗,𝐶𝑞,𝑞,𝑞=−12𝑎³(1−𝑎)𝜔₀−8𝑐𝑎³𝜔₀+4𝑐𝑎^{3(3−5𝑎)𝑖}(1−𝑎)(1+𝑎)²𝜔₀. (3.19) Using conclusions in [14] we can get 𝐶₁(𝑏)=𝑔₂₀𝑔₁₁(3𝛼(𝑏)+𝑖𝜔(𝑏))2𝛼²(𝑏)+𝜔²(𝑏)+𝑔₁₁²𝛼(𝑏)(𝑏)+𝑔+𝑖𝜔2(𝛼(𝑏)⁰²²(𝑏+𝑔+3𝑖𝜔))2²¹, (3.20) then 𝐶₁𝑏₀=𝑖2𝜔₀𝑔₂₀𝑔₁₁−2𝑔₁₁²−13𝑔₀₂²+𝑔₂₁2,𝜇₂𝑐=−Re₁(0)Re𝑏^′𝜏_𝑛,𝛽₂=2Re𝑐₁(0,𝑇₂)=2𝜋𝜔₀1+𝜏₂𝑠²+𝑜𝑠⁴, (3.21) where 𝜏₂1=−𝜔₀Im𝑐₁𝑏₀−Re𝑐₁𝑏₀𝛼^′𝑏₀𝜔^′𝑏₀. (3.22)

Now we give a conclusion.

Conclusion. (1) The sign of 𝜇₂ determines the direction of Hopf bifurcation. When 𝜇₂>0, the Hopf bifurcation is supercritical; when 𝜇₂<0, the Hopf bifurcation is subcritical.
(2) 𝛽₂ determines the stability of bifurcated periodic solutions. When 𝛽₂<0, the periodic solutions are stable; when 𝛽₂>0, the periodic solutions are unstable.
(3) 𝑇₂ determines the period of bifurcated periodic solutions. When 𝑇₂>0, the period increases; when 𝑇₂<0, the period decreases.

4. Example

In this section, we use a numerical simulation to illustrate the analytical results we obtained in previous sections.

Let 𝑥∈(0,𝑙𝜋),𝑑₁=1,𝑑₂=3,𝑐=4,𝑎=0.0588. The system (1.5) is 𝑢_𝑡−𝑢_𝑥𝑥=𝑢−0.0588𝑢²−𝑢𝑣𝑏(),𝑣𝑢+1_𝑡−3𝑣_𝑥𝑥=−4𝑣+𝑢𝑣𝑏(),𝑢𝑢+1_𝑥(0,𝑡)=𝑣_𝑥(0,𝑡)=0,𝑢_𝑥(𝑙𝜋,𝑡)=𝑣_𝑥()𝑢()()().𝑙𝜋,𝑡=0,𝑥,0≥0,𝑣𝑥,0≥0,𝑥∈0,𝑙𝜋 (4.1) Now we determine the direction of a Hopf bifurcation with 𝑏∈𝐵 and the other properties of bifurcating periodic solutions based on the theory of Hassard et al. [14], as discussed before. By means of software MATLAB 7.0, we can get some figures to illustrate the effectiveness of our results. 𝑏^𝐵₀=0.2222, 𝑏^∗=0.1667, 𝐴(𝑏^∗)=0.4706, 𝑙_𝑛=2.9155𝑛, and (2.27) is held. When 𝑛=2, 𝑙=3.1162, (𝑙_𝑛,𝑙_𝑛+1]=(2.9155,5.8302], 𝑙∈(𝑙_𝑛,𝑙_𝑛+1]. We can get 𝐵=(0.1253,0.1944,0.2222). The only positive equilibrium point of (4.1) is 𝑌₀=(4𝑏/(1−4𝑏),𝑏(0.0588−𝑢₀)(1+𝑢₀)/0.0588). When 𝑏=0.1253, we can compute Re𝑐₁(0.1253)=0.2856>0,𝜇₂=−0.3423<0,𝑇₂=1.2346>0. (4.2) The positive equilibrium point of (4.1) is unstable and the Hopf bifurcation is supercritical. The positive equilibrium point 𝑌₀ of system (4.1) is locally asymptotically stable when 𝑏=0.1 as is illustrated by computer simulations in Figure 1. And periodic solutions occur from 𝑌₀ when 𝑏=0.1253 as is illustrated by computer simulations in Figure 2. When 𝑏=1.3, we can easily show that the positive equilibrium point 𝑌₀ is unstable as is illustrated in Figure 3. From the above results, we can conclude that the stability properties of the system could switch with parameter 𝑏.

(a)

(b)

Figure 1 

When 𝑏=0.1, the positive equilibrium point 𝑌0 is asymptotically stable.

(a)

(b)

Figure 2 

When 𝑏=0.1253, periodic solutions occur from 𝑌0.

(a)

(b)

Figure 3 

When 𝑏=1.3, the positive equilibrium point 𝑌0 is unstable.

Acknowledgments

This research was supported by the National Natural Science Foundations of Heilongjiang province (C200941). The authors are grateful to the anonymous reviewers for their valuable comments and suggestions which essentially improve the quality of the paper.