# Dynamic Properties of a Forest Fire Model

Authors: Na Min; Hongyang Zhang; Zhilin Qu
Journal: Abstract and Applied Analysis
 (2012)

---

## 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.


---

## Body

## 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:
(1.1)ut=d1uxx+u+f(u),x∈R,t≥0,
where u=u(x,t) can be seen as the area of the burned forest. uxx is a diffusion term of u in space, and d1 is the diffusion coefficient. f(u) 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)vt=d2vxx-cv+g(v).
In this equation, v=v(x,t) is the area where the fire has been put out. vxx is a diffusion term of v in space, and d2 is the diffusion coefficient. c is the resurgence probability of v. g(v) 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,u and v influence each other. Thus, f and g must be functions of u,v. We define g(u,v) by referring to the combustion model [13]:
(1.3)g(u,v)=uvb(u+1),b>0.Sinceg(u,v) has opposite effect on the fire area (or u), we can also define f(u,v) by taking into account KPP model [8]:
(1.4)f(u,v)=-au2-uvb(u+1).Then we get a new model:(1.5)ut-d1uxx=u-au2-uvb(u+1),vt-d2vxx=-cv+uvb(u+1),ux(0,t)=vx(0,t)=0,ux(lπ,t)=vx(lπ,t)=0,u(x,0)≥0,v(x,0)≥0,x∈(0,lπ).
Define
(1.6)H=L2[0,π]×L2[0,π]={(fg):f,g∈L2[0,π]},
and an inner product is given by
(1.7)〈(f1g1),(f2g2)〉=〈f1,f2〉L2+〈g1,g2〉L2=2π∫0π‍f1¯f2dx+2π∫0π‍g1¯g2dx,
where (f1,g1)T∈H,(f2,g2)T∈H. From the standpoint of biology, we are only interested in the dynamics of model (1.5) in the region:
(1.8)R+2={(u,v)∣u>0,v>0}.


## 2. Stability Analysis

Firstly, we consider the location stability [14, 15] and the number of the equilibria of model (1.5) in R+2. We can also study autowave solutions [16] of the model. The interior equilibrium point is a root of the following equation:
(2.1)u-au2-uvb(u+1)=0,-cv+uvb(u+1)=0.
It is obvious that (2.1) has an only real solution Y0=(u0,v0), where
(2.2)u0=bc1-bc,v0=ab(1a-u0)(1+u0),
and b<1/c(a+1).Now, we analyze the asymptotic stability of(u0,v0) by Lyapunov function.Lemma 2.1.
For the model (1.5),
(1) ifa≥1,  (u0,v0) is global asymptotic stability.
(2) ifa<1 and (1-a)/c≤b≤1/(ac+c),(u0,v0) also has global asymptotic stability.Proof.
Defining(2.3)ω(u,v)=∫0lπ‍∫u0u‍r/b(r+1)-cr/(r+1)drdx+1b∫0lπ‍∫v0v‍s-v0sdsdx,
we can get
(2.4)∂ω∂t=1b∫0lπ‍(h(u)-h(u0))(p(u)-p(u0))dx+Y(t),
where(2.5)h(u)=uu+1,p(u)=(1-au)(u+1),Y(t)=-d1c∫0lπ‍h′(u)h2(u)ux2dx+d2v0b∫0lπ‍vx2v2dx.
In what follows, we split it into two cases to prove. If a≥1, for all u>0,p'(u)<0, since h'(u)=1/(u+1)2>0, we can get
(2.6)(h(u)-h(u0))(p(u)-p(u0))≤0.
If a<1 and (1-a)/c≤b≤1/(ac+c) (equal to v0≤b), we can still get (2.6).  That is to say
(2.7)wt(u,v)<0.
We prove the conclusion.Because of the conclusion of Lemma2.1, we always assume a<1 and 0<b<(1-a)/(2ac-c). Introducing perturbations u*=u-u0,v*=v-v0, and replace (u*,v*) with (u,v), for which model (1.5) yields
(2.8)ut-d1uxx=u+u0-a(u+u0)2-(u+u0)(v+v0)b(u+u0+1),vt-d2vxx=-c(v+v0)+(u+u0)(v+v0)b(u+u0+1),ux(0,t)=vx(0,t)=0,ux(lπ,t)=vx(lπ,t)=0,u(x,0)≥0,v(x,0)≥0,x∈(0,lπ).
Now we can get the linearized system of parametric model (2.8) at (0,0),
(2.9)(utvt)=(Δ~+L(b))(uv),
where
(2.10)Δ~=(0d1∂2∂x2d2∂2∂x20),L(b)=(u0(1-a-au0)1+u0-c1-au01+u00).
The eigenvalues of Δ~ are as follows:
(2.11){-d2n2l2,-d1n2l2}n=0+∞,
and the corresponding eigenvectors as follows:
(2.12){βn1,βnk2}n=0+∞,
where
(2.13)βn1=(0cos(nlx)),βn2=(cos(nlx)0).
Define  for all y∈H(2.14)y=∑k=1n‍YkT(βk1βk2),Yk=(〈y,βk1〉〈y,βk2〉).
It is easy to get, λ∈(Δ~+L(b)), if and only if the equation
(2.15)∑k=1n‍YkT(Eλ-L(b)-(-d1k2l200-d2k2l2))(βk1βk2)=0
is held.We obtain(2.16)|Eλ-L(b)-(-d1n2l200-d2n2l2)|=0.
Rewrite it as
(2.17)λ2-Tn(b)λ+Dn(b)=0,
where
(2.18)Tn(b)=bc(1-a-2abc1-bc)-(d1+d2)n2l2,Dn(b)=c(1-bc-abc)-d2bc(1-a-2abc1-bc)n2l2+d1d2n4l4.
From (2.18), when (1-a)/(ac+c)<b<(1-a)/c is held, we can get Tn(b)<0,Dn(b)>0. So the system's eigenvalues have negative real part, and (u0,v0) has local asymptotic stability. Then, we can conclude that the system has Hopf bifurcation [14] in b∈(0,(1-a)/(ac+c)).Define(2.19)B={b0∣Tn(b0)=0,Dn(b0)>0,Tj(b0)≠0,Dj(b0)≠0,∀j≠n}.b0∈B∪(0,(1-a)/(ac+c)), α(b)±iω(b) are characteristic roots of Δ~+L(b), where
(2.20)α(b)=12(A(b)-d1+d2l2n2),ω(b)=Dn(b)-α2(b),A(b)=bc(1-a-2abc1-bc).
Now we compute transversality condition:
(2.21)α′(b0)=12a(1-b0c)2(1a-1-4b0c1-b0c-(b0c1-b0c)2)×{>0,0<b0<1c(1-2a1+a),<0,1c(1-2a1+a)<b0<1-aac+c.
Now we consider A(0)=A(b0B)=0 and A(b) is positive in (0,b0B). So we can get the maximum value of A(b) (defined asA(b*)):Define(2.22)ln=nd1+d2A(b*),n∈ℕ,A(b*)=a(1a+1-2)2,
for all l∈(ln,ln+1],0≤j≤n,bj,-B and bj,+B are two roots of the equation
(2.23)A(b)=d1+d2l2j2.
It is easy to get
(2.24)0<b1,-B<⋯<bn,-B<1c(1-2a1+a)<bn,+B<⋯<b1,+B<b0B.
Then we give the condition of Dn(bj,±B)≠0 especially Dn(bj,±B)>0.As we know(2.25)Dn(b)≥ac-d2A(b*)n2l2+d1d2n4l4.
Then Dn(b)>0 is held if and only if
(2.26)d1d2>0,(d2A(b*))2-4d1d2ac>0.Theorem 2.2.
Assume thatd1,d2,c>0,0<a<1, and the equation is held:
(2.27)d2A2(b*)-4d1ac>0,
where
(2.28)ln=nd1+d2A(b*),n∈ℕ,A(b*)=a(1a+1-2)2,
then for all l∈(ln,ln+1], existing b=bj,±B or b=b0B; there are Hopf bifurcations at the real solution of model (1.5).
Furthermore(2.29)0<b1,-B<⋯<bn,-B<1c(1-2a1+a)<bn,+B<⋯<b1,+B<b0B.


## 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 b0 or bj,±(j=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 transformationu*=u-u0,v*=v-v0, and replacing (u*,v*) with (u,v), the parametric system (1.5) is equivalent to the following functional differential equation (FDE) system:
(3.1)∂U∂t=(Δ~+L(b0))+F(b0,U),
where
(3.2)U=(u,v)T,F(b0,U)=(f-fuu-fvvg-guu-gvv).
The adjoint operator of Ln(b) is defined as
(3.3)L*(b0)=(d1∂2∂x2+fu(b0)fv(b0)gu(b0)d2∂2∂x2+gv(b0)).
It is easy to get
(3.4)〈u,L(b0)v〉=〈L*(b0)u,v〉.
From the discussions in Section 2, define q*=(an*,bn*)Tcos(n/l)x. We have
(3.5)L*(b0)q*=-iωq*,〈q*,q〉=1,〈q*,q¯〉=0.
Decompose X as X=Xc⊕Xs, where Xc={zq+zq¯∣z∈C} and Xs={u∈x∣〈q*,u〉=0}. For all (u,v)∈X, existing z∈ℂ and ω=(ω1,ω2)∈Xs,  we can obtain
(3.6)(uv)=zq+zq¯+(ω1ω2).
Rewrite (3.1) as
(3.7)z˙=iω0z+〈q*,F0*〉,ω˙=L(b0)ω+H(z,z¯,ω),
where
(3.8)F0*=zq+zq¯+ω,H(z,z¯,ω)=F0*-〈q*,F0*〉q-〈q*¯,F0*〉q¯.
Using the same notations as in [11],
(3.9)F0*(U)=12Q(U,U)+16C(U,U,U)+O|U4|,
where U=(u,v) and Q,C are symmetrical multilinear functions. We can compute
(3.10)Q(q,q)=(An1An2)cos2nlx,Q(q,q¯)=(Bn1Bn2)cos2nlx,C(q,q,q¯)=(Cn1Cn2)cos3nlx,
where
(3.11)An1=fuuan2+2fuvanbn+fvvbn2,An2=guuan2+2guvanbn+gvvbn2,Bn1=fuu|an|2+fuv(anbn¯+an¯bn)+fvv|bn|2,Bn2=guu|an|2+guv(anbn¯+an¯bn)+gvv|bn|2,Cn1=fuuu|an|2an+fuuv(2|an|2bn+an2bn¯)+fuvv(2bn2an+bn2an¯)+fvvv|bn|2bn,Cn2=guuu|an|2an+guuv(2|an|2bn+an2bn¯)+guvv(2bn2an+bn2an¯)+gvvv|bn|2bn.
Define
(3.12)H(z,z¯,ω)=12H20z2+H11zz¯+12H02z¯2+o(|z||ω|),
where
(3.13)H20=Q(q,q)-〈q*,Q(q,q)〉q-〈q*¯,Q(q,q)〉q¯,H11=Q(q,q¯)-〈q*,Q(q,q¯)〉q-〈q*¯,Q(q,q¯)〉q¯.
On the center manifold, we have
(3.14)ω=12ω20z2+ω11zz¯+12ω02z¯2+o(|z|3).
We can obtain
(3.15)w20=[2iω0I-Ln(b0)]-1H20,ω11=[Ln(b0)]-1H11.
Comparing (3.9) and (3.13), we can get
(3.16)H20={Q(q,q)-(An1An2)cos2nlx-(An1An2)(12cos2n2x+12),n∈ℕ*,(A01A02)-〈q*,Q(q,q)〉(a0b0)-〈q*,Q(q,q)〉(a0¯b0¯),n=0.
Similarly
(3.17)ω11={-12[L(bo)]-1(Cn1Cn2)(cos2n2+1),n∈ℕ*,-[L(b0)]-1[(c01c02)-〈q*,Q(q,q¯)〉(a01b02)-〈q*¯,Q(q,q¯)〉(a0¯b0¯)],n=0.
Then on the center manifold rewrite (dU) as
(3.18)dzdt=iω0z+〈q*,F0*〉=iω0z+∑2≤i+j≤3‍giji!j!ziz¯j+o(|z|4),
where
(3.19)g20=〈q*,Q(q,q)〉=4cω0a2-a(1-a)2ω0-2ca2(3a-1)i(1-a2)ω0,g11=〈q*,Q(q,q¯)〉=a(a-q)ω0-2ca2i(1+a)ω0,g02=〈q*,Q(q¯,q¯)〉=a(1-a)2+2cω0a2-4ca2i(1-a2)ω0,g21=2〈q*,Q(w11,q)〉+〈q*,Q(w20,q¯)〉+〈q*,C(q,q,q¯)〉=-12a3(1-a)ω0-8ca3ω0+4ca3(3-5a)i(1-a)(1+a)2ω0.
Using conclusions in [14] we can get
(3.20)C1(b)=g20g11(3α(b)+iω(b))2(α2(b)+ω2(b))+|g11|2α(b)+iω(b)+|g02|22(α(b)+3iω(b))+g212,
then
(3.21)C1(b0)=i2ω0(g20g11-2|g11|2-13|g02|2)+g212,μ2=-Re{c1(0)}Re{b′(τn)},β2=2Re{c1(0)},T2=2πω0(1+τ2s2)+o(s4),
where
(3.22)τ2=-1ω0[Im(c1(b0))-Re(c1(b0))α′(b0)ω′(b0)].Now we give a conclusion.Conclusion.
(1) The sign ofμ2 determines the direction of Hopf bifurcation. When μ2>0, the Hopf bifurcation is supercritical; when μ2<0, the Hopf bifurcation is subcritical.
(2)β2 determines the stability of bifurcated periodic solutions. When β2<0, the periodic solutions are stable; when β2>0, the periodic solutions are unstable.
(3)T2 determines the period of bifurcated periodic solutions. When T2>0, the period increases; when T2<0, the period decreases.


## 4. Example

In this section, we use a numerical simulation to illustrate the analytical results we obtained in previous sections.Letx∈(0,lπ),d1=1,d2=3,c=4,a=0.0588. The system (1.5) is
(4.1)ut-uxx=u-0.0588u2-uvb(u+1),vt-3vxx=-4v+uvb(u+1),ux(0,t)=vx(0,t)=0,ux(lπ,t)=vx(lπ,t)=0,u(x,0)≥0,v(x,0)≥0,x∈(0,lπ).
Now we determine the direction of a Hopf bifurcation with b∈B 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. b0B=0.2222, b*=0.1667, A(b*)=0.4706, ln=2.9155n, and (2.27) is held. When n=2, l=3.1162, (ln,ln+1]=(2.9155,5.8302], l∈(ln,ln+1]. We can get B=(0.1253,0.1944,0.2222). The only positive equilibrium point of (4.1) is Y0=(4b/(1-4b),b(0.0588-u0)(1+u0)/0.0588). When b=0.1253, we can compute (4.2)Rec1(0.1253)=0.2856>0,μ2=-0.3423<0,T2=1.2346>0.
The positive equilibrium point of (4.1) is unstable and the Hopf bifurcation is supercritical. The positive equilibrium point Y0 of system (4.1) is locally asymptotically stable when b=0.1 as is illustrated by computer simulations in Figure 1. And periodic solutions occur from Y0 when b=0.1253 as is illustrated by computer simulations in Figure 2. When b=1.3, we can easily show that the positive equilibrium point Y0 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 b.Whenb=0.1, the positive equilibrium point Y0 is asymptotically stable.
(a)
(b)Whenb=0.1253, periodic solutions occur from Y0.
(a)
(b)Whenb=1.3, the positive equilibrium point Y0 is unstable.
(a)
(b)



---

*Source: 948480-2012-10-23.xml*