# Complex Behavior Analysis of a Fractional-Order Land Dynamical Model with Holling-II Type Land Reclamation Rate on Time Delay

Authors: Li Wu; Zhouhong Li; Yuan Zhang; Binggeng Xie
Journal: Discrete Dynamics in Nature and Society
 (2020)

---

## Abstract

In this paper, a fractional-order land model with Holling-II type transformation rate and time delay is investigated. First of all, the variable-order fractional derivative is defined in the Caputo type. Second, by applying time delay as the bifurcation parameter, some criteria to determine the stability and Hopf bifurcation of the model are presented. It turns out that the time delay can drive the model to be oscillatory, even when its steady state is stable. Finally, one numerical example is proposed to justify the validity of theoretical analysis. These results may provide insights to the development of a reasonable strategy to control land-use change.


---

## Body

## 1. Introduction

Since the 20th century, the global environment has changed at an unprecedented speed, and a series of major global environmental problems pose a serious threat to the survival and development of mankind. From the perspective of earth system science, the generation of global environmental problems can be recognized as the result of the interaction between the earth’s atmosphere, hydrosphere, biosphere, lithosphere, and human activities [1, 2]. Land use/land cover change (LUCC) is as a key link connecting the four circles of the earth system, being an ideal entry point for the study of natural and human processes, and has rapidly become the focus of global change research [3, 4]. LUCC dynamic change is to describe, evaluate, interpret, and predict the quantity, quality, spatial distribution, classification, rate, and process of land use/land cover change by accumulating spatiotemporal continuous land use/land cover data and using the mathematical model [5]. The driving forces of land use change are mainly divided into natural and social systems. Notably, hydrology, topography, and geomorphology [6, 7] are the main driving forces in the natural system, while economic development, population growth, and policies are the main driving forces in the social and economic systems [8]. Generally speaking, it is the interwoven factors in these two systems that cause land use changes [9], but compared with the natural factors (such as climate, soil, and topography), human factors (such as regional policy, economic development, and population growth) play a dominant role [10].Since the key of LUCC research is to understand the driving force and driving mechanism [11], a variety of system analysis and mathematical statistic methods have been well applied in the research of driving force. For example, the internal and human driving force of land use change is studied through the analysis of the multiple linear regression model [12]. The system dynamic model was used to quantitatively diagnose the contribution of each driving factor to land use change [13]. The CLUE comprehensive model was used to establish the interaction between land-use change and its influencing factors, and land change and its spatial distribution were simulated [14–16]. Therefore, the comprehensive mathematical model will be a development trend of LUCC research at present and in the future.In 1997, Dobson et al. considered the land model with four coupled differential equations [17]:(1)F˙t=−dPtFt+sUt,A˙t=dPtFt−aAt+bUt,U˙t=aAt−sUt−bUt,P˙t=rPt1−hAtPt,where F stands for an original area of pristine forest habitat, A denotes the agricultural land, unused land is U, and P is the population density with time t, for more detail, one can read [17] and the references cited therein. In [18], Chen et al. studied the dynamical analysis of the land model with Holling-II type land reclamation rate, and they found the positive equilibrium point, the basic reproductive rate R0, and, furthermore, some sufficient conditions for the global stability of the positive equilibrium and one of the boundary equilibrium.On the one hand, it is well known that time delays are unavoidable for population modelling. It is therefore important to consider the dynamics for the population models with time delays, and dynamic complex analysis is obviously one of the most important problems. For example, delay differential equations have been applied in the stability of population dynamic systems [19–25], impulsive effects, and control systems [26–28].On the other hand, there are two types of differential equations, namely, integer-order and fractional-order differential equations. Traditionally, the fractional-order differential equations have enjoyed a preference over integer-order differential equations because of the mathematical tractability of fractional systems. Over the past 30 decades, fractional-order differential equations have been applied in many fields, such as design and control of various ecological systems [29–32], secure communication [33–35], and system control [36–42].Recently, the existence of Hopf bifurcation of differential equations has been studied as an important qualitative behavior of integer-order differential equations [43]. Moreover, fractional calculus is merged into complicated, dynamical systems which extremely renovate the theory of the design and control performance for complex systems. It has been discovered that physical phenomena in nature can be depicted more accurately by fractional-order models in comparison with classical integer-order ones [44, 45]. Some scholars introduced fractional calculus into population models and constructed fractional population or epidemic models [46–50], fractional dynamic systems, and fractional neural networks [51–53]. However, to the best of our knowledge, there are few studies to investigate the existence of Hopf bifurcation to the fractional-order land model with time delay.Motived by the above ideas, in this paper, we will consider the following fractional-order land model with Holling-II type response and time delay:(2)DϕFt=−dFtPt−τ1+Pt−τ+sUt,DϕRt=dFtPt−τ1+Pt−τ−aRt+bUt,DϕUt=aRt−sUt−bUt,DϕPt=rPt1−hRtPt−τ,where ϕ∈0,1 is the fractional order, Ft,Rt,andUt represent the land of wood and grass, survival land, and unused land, respectively, and Pt is the density of the predator population at time t; Rt (land area of survival) of a period of time 1/a becomes unused land (area Ut), which in turn recovers through natural succession or ecological restoration to become a forest after a time interval 1/s. The unused land may also be restored to survival viable land after a time interval 1/b. Average clearing ability is described by constant d. r is the growth rate into the population. τ>0 is the time delay required for the gestation of the mature population.The main contributions can be summed up in three aspects:(1)
A novel fractional-order land model with Holling-II type land reclamation rate and time delay is formulated(2)
Two primary dynamical properties—stability and oscillation—of the delayed fractional-order land model are investigated(3)
The influences of the order on the Hopf bifurcation are obtainedThroughout this paper, we address the following assumptions.Assumption 1.
bd+ds−ahs>0.Assumption 2.
abd+ads−ha2s>0.
Suppose Assumptions1 and 2 hold, then system (2) has an unique positive equilibrium point E∗=F∗,R∗,U∗,P∗, which is described by(3)F∗=asah+bh+b+hs+sb+sad+as+bd+ds,R∗=bd+ds−ahsad+as+bd+ds,U∗=abd+ads−ha2sb+sad+as+bd+ds,P∗=bd+ds−ahshad+as+bd+ds.
The remainder of the current paper is organized as follows: In Section2, some definitions and lemmas of fractional calculus are recalled. In Section 3, the fractional land model with Holling-II type land reclamation rate and time delay is investigated by using time delay as the bifurcation parameter, and the conditions of Hopf bifurcation are presented. One numerical example is given to illustrate the effectiveness of our main results in Section 4. Finally, conclusions are drawn in the last section.


## 2. Basic Tools of Fractional Calculus

In this section, we recall some definitions and lemmas of fractional calculus, which can be used in the proofs of the main results in Section3.There are several definitions of fractional derivatives. The Gru¨nwald–Letnikov definition, the Riemann–Liouville definition, and the Caputo definition are usually used to deal with fractional-order systems. Since the Caputo derivative only requires the initial conditions which are based on the integer-order derivative and represents well-understood features of physical situation, it is more applicable to real-world problems. Hence, the Caputo fractional-order derivative is adopted in this paper.Definition 1.
(see [44]). The fractional integral of order α for a function gt is defined as(4)Iαgt=1Γα∫t0tt−sα−1gsds,where t≥t0, α>0, Γ⋅ is the gamma function, and Γs=∫0∞ts−1e−tdt.Definition 2.
(see [44]). The Caputo fractional derivative of order α for a function gt∈Cnt0,∞,R is defined by(5)Dαgt=1Γn−α∫t0tgnst−sα−n+1ds,where t≥t0 and n is a positive integer such that n−1≤α<n.
Moreover, when0<α<1,(6)Dαgt=1Γ1−α∫t0tg′st−sαds.Lemma 1.
(see [41]). The following autonomous system(7)Dαy=Jy,y0=y0,where 0<α<1, y∈Rn, and A∈Rn×n, is asymptotically stable if and only if argλi>απ/2i=1,2,…,n. In this case, each component of the states decays towards 0 like t−α. Also, this system is stable if and only if argλi≥απ/2, and those critical eigenvalues that satisfy argλi=απ/2 have geometric multiplicity one.


## 3. Main Results

In this section, by applying the previous analytic technique, we shall investigate the stability and bifurcation of system (2) with time delay by taking time delay as the bifurcation parameter. The conditions of delay-induced bifurcation will be derived.From reference [18], the dynamics of system (2) crucially depend on the basic reproduction number R0, which is given by(8)R0=bd+dsash.Letx1t=Ft−F∗,x2t=Rt−R∗,x3t=Ut−U∗,andx4t=Pt−P∗, and then we have(9)Dϕx1t=−dx1t+F∗x4t−τ+P∗1+x4t−τ+P∗+sx3t+U∗,Dϕx2t=dx1t+F∗x4t−τ+P∗1+x4t−τ+P∗−ax2t+R∗+bx3t+U∗,Dϕx3t=ax2t+R∗−sx3t+U∗−bx3t+U∗,Dϕx4t=rx4t+P∗1−hx2t+R∗x4t−τ+P∗.Taking advantage of the Taylor expansion formula, the linearized system of system (9) at the zero equilibrium is(10)Dϕx1t=a11x1t+a13x3t+b14x4t−τ,Dϕx2t=a21x1t+a22x2t+a23x3t+b24x4t−τ,Dϕx3t=a32x2t+a33x3t,Dϕx4t=a42x2t+b44x4t−τ,where a11=−dP∗/1+P∗,a13=s,b14=−dF∗/1+P∗2,a21=dP∗/1+P∗,a22=−a,a23=b,b24=dF∗/1+P∗2,a32=a,a33=−b+s,a42=hrP∗2/R∗2,a44=r−hrP∗/R∗,andb44=−hrP∗/R∗.The associated characteristic equation of system (10) can be obtained as(11)detsϕ−a110−a13−b14e−sτ−a21sϕ−a22−a23−b24e−sτ0−a32sϕ−a3300−a420sϕ−a44−b44e−sτ=0,which leads to(12)P1s+P2se−sτ=0,where(13)P1s=−a11s2ϕ−a22s2ϕ−a33s2ϕ+a11a22sϕ−a23a32sϕ+a11a33sϕ+a22a33sϕ−a12a21a32+a11a23a32−a11a22a33+s3ϕsϕ−a44,P2s=a21a12a32b44−a42b14sϕ−a33−sϕ−a11a42b24sϕ−a33+b44−a33sϕ+a22a33−sϕ−a23a32+s2ϕ.LetP1s=A1+iB1 and P2s=A2+iB2, and from equation (12), we have(14)A1+iB1+A2+iB2e−sτ=0.Assume thats=iw=wcosπ/2+isinπ/2 is a root of equation (12), ω>0. Substituting s into equation (12) and separating the real and imaginary parts, one can have(15)A2coswτ+B2sinwτ=−A1,B2coswτ−A2sinwτ=−B1.Applying equation (15), direct calculation yields(16)coswτ=−A1A2+B1B2A22+B22=ℱw,sinwτ=−−A2B1+A1B2A22+B22=Gw,where(17)A1=−a23a32+a22a33+a22+a33a44+a11a22+a33+a44w2ϕcosπϕ−a11+a22+a33+a44w3ϕcos3πϕ2+a12a21a32+a22a33−a23a32a44+a11−a23a32+a22a33+a22+a33a44wϕ−cosπϕ2+a12a21a32+a11a22a33−a23a32a44+w4ϕcos2πϕ,B1=a21a12a32b44−a42b14wϕcosπϕ2−a33+a42b24−a11+a33wϕ−cosπϕ2+a11a33+w2ϕcosπϕ+b44a11+a22+a33w2ϕcosπϕ+a11a22−a23a32+a11+a22a33wϕ−cosπϕ2+a11a22a33−a23a32−w3ϕcos3πϕ2,A2=a21a12a32b44−a42b14wϕcosπϕ2−a33+a42b24−a11+a33wϕ−cosπϕ2+a11a33+w2ϕcosπϕ+b44a11+a22+a33w2ϕcosπϕ+a11a22−a23a32+a11+a22a33wϕ−cosπϕ2+a11a22a33−a23a32−w3ϕcos3πϕ2,B2=−wϕsinπϕ2a42b24−a11−a33+2wϕcosπϕ2+b44−2a11+a22+a33×wϕcosπϕ2+a11a22−a23a32+a11a33+a22a33+w2ϕ2cosπϕ+1+a21a42b14.From (16), we obtain(18)τk=1warccosℱw+2kπ,k=0,1,2,….Define(19)τ0∗=minτk,k=0,1,2,…,where τk is defined by equation (18).To derive the condition of the occurrence for Hopf bifurcation, we further give the following assumption.Assumption 3.
(20)W1V1+W2V2V12+V22≠0,where(21)W1=w0A2sinw0τ0−B2cosw0τ0,W2=w0A2cosw0τ0+B2sinw0τ0,V1=A1′+A2′−τ0A2cosw0τ0+B2′−τ0B2sinw0τ0,V2=B1′+B2′−τ0B2cosw0τ0A2′−τ0A2sinw0τ0.Lemma 2.
Letsτ=ντ+iwτ be the root of equation (18) near τ=τj satisfying ντj=0, wτj=w0, then the following transversality condition holds(22)Redsdτw=w0,τ=τ0≠0.Proof.
By using implicit function theorem and differentiating (12) with respect to τ, we have(23)P1′sdsdτ+P2′se−sτdsdτ+P2se−sτ−τdsdτ−s=0,where Pi′s are the derivatives of Pisi=1,2.
Thus, we obtain(24)dsdτ=WsVs,where(25)Ws=sP2se−sτ,Vs=P1′s+P2′se−sτ−τP2se−sτ.
LetPiR and PiI be the real and imaginary parts of Pisi=1,2, respectively. We further suppose that W1 and W2 are the real and imaginary parts of Ws, respectively, and V1, V2 are the real and imaginary parts of Vs, respectively.
Hence,(26)Redsdττ=τ0,w=w0=W1V1+W2V2V12+V22.
Therefore, based on Assumption3, the transversality condition is satisfied. This completes the proof of Lemma 2.
Denote(27)β1=α1,β2=α11α3α2,β3=α110α3α2α10α4α3,β4=α3β3,where(28)α1=−a11−a22−a33−a44−b44,α2=a22b44−a42b24+a11b44+a33b44+a11a22+a33a22+a44a22−a23a32+a11a33+a11a44+a33a44,α3=a23a32b44−a21a42b14+a11a42b24+a33a42b24−a11a22b44−a11a33b44−a22a33b44−a13a21a32+a11a23a32+a23a44a32−a11a22a33−a11a22a44−a11a33a44−a22a33a44,α4=a21a33a42b14−a11a33a42b24+a13a21a32a44−a11a23a32a44+a11a22a33a44+a13a21a32b44−a11a23a32b44+a11a22a33b44.Lemma 3.
If Assumptions1 and 2 hold, then the positive equilibrium point F∗,R∗,U∗,P∗ of the delayed fractional-order land model (2) is asymptotically stable when τ=0.Proof.
Whenτ=0, system (2) changes into(29)λ4+β1λ3+β2λ2+β3λ+β4=0.
If the conditions ofβi>0i=1,2,3,4 hold and, furthermore, if Assumptions 1 and 2 are satisfied, it is easy to check from the Routh–Hurwitz criterion that the four eigenvalues of the characteristic equation (29) have negative real parts. Hence, the positive equilibrium point F∗,R∗,U∗,P∗ of the delayed fractional-order land model (2) is asymptotically stable when τ=0.
From the above analysis, together with Lemma3, we will give the following theorem.Theorem 1.
For the model (2):(i)
Under Assumptions1–3, the equilibrium point E∗=F∗,R∗,U∗,P∗ is global asymptotically stable for τ∈0,+∞(ii)
Under Assumptions1–3:(a)
The equilibrium pointE∗=F∗,R∗,U∗,P∗ is locally asymptotically stable for τ∈0,τ0(b)
Model (2) undergoes a Hopf bifurcation at the E∗=F∗,R∗,U∗,P∗ when τ=τ0, i.e., it has a branch of periodic solutions bifurcating from the equilibrium near τ=τ0


## 4. Numerical Example and Application

In this section, we give one example to demonstrate the effectiveness of the proposed approach. The simulation results are based on theAdams–Bashforth–Moulton predictor-corrector scheme [54], and step length Δt=0.01.Consider the following fractional-order land dynamical model with time delay:(30)DϕFt=−dFtPt−τ1+Pt−τ+sUt,DϕRt=dFtPt−τ1+Pt−τ−aRt+bUt,DϕUt=aRt−sUt−bUt,DϕPt=rPt1−hRtPt−τ,in this case, taking ϕ=0.9, and according to Yuxi city land and population data of Yunan Province, China, let a=0.34,b=0.012,r=0.0048,s=0.004,h=0.05,andd=0.08. By computing, we obtain the positive equilibrium point as 0.075617,0.021169,0.44984,0.423379. Moreover, we have R0=bd+ds/ash=18.8235>1, ω0=0.002358, and τ0=728.403. Furthermore, the initial values are selected as F0,R0,U0,P0=0.07,0.02,0.4,0.75, and it is easy to verify that the conditions of (ii) in Theorem 1 are satisfied. Figures 1–8 depict the equilibrium point 0.075617,0.021169,0.44984,0.423379 which is locally asymptotically stable when τ=700<τ0, while Figures 9–16 display the equilibrium point 0.075617,0.021169,0.44984,0.423379 that is unstable, and Hopf bifurcation occurs when τ=τ0 (see Supplementary Materials (available here)).Figure 1
Time series solution of system (30) with ϕ=0.90, τ=700<τ0=728.403.Figure 2
Time series solution of system (30) with ϕ=0.90, τ=700<τ0=728.403.Figure 3
Time series solution of system (30) with ϕ=0.90, τ=700<τ0=728.403.Figure 4
Time series solution of system (30) with ϕ=0.90, τ=700<τ0=728.403.Figure 5
Portrait diagram of land model (30) with ϕ=0.90, τ=700<τ0=728.403.Figure 6
Portrait diagram of land model (30) with ϕ=0.90, τ=700<τ0=728.403.Figure 7
Portrait diagram of land model (30) with ϕ=0.90, τ=700<τ0=728.403.Figure 8
Portrait diagram of land model (30) with ϕ=0.90, τ=700<τ0=728.403.Figure 9
Time series solution of system (30) with ϕ=0.90, τ=740>τ0=728.403.Figure 10
Time series solution of system (30) with ϕ=0.90, τ=740>τ0=728.403.Figure 11
Time series solution of system (30) with ϕ=0.90, τ=740>τ0=728.403.Figure 12
Time series solution of system (30) with ϕ=0.90, τ=740>τ0=728.403.Figure 13
Portrait diagram of land model (30) with ϕ=0.90, τ=740>τ0=728.403.Figure 14
Portrait diagram of land model (30) with ϕ=0.90, τ=740>τ0=728.403.Figure 15
Portrait diagram of land model (30) with ϕ=0.90, τ=740>τ0=728.403.Figure 16
Portrait diagram of land model (30) with ϕ=0.90, τ=740>τ0=728.403.


## 5. Conclusion

The problem of bifurcation for a fractional land model with time delay has been studied in this paper. Firstly, by comprehensively analyzing the characteristic equation, the stability criterion and the conditions of existence of Hopf bifurcation are established. It was shown that the presence of delay heavily affects the stability of the model. Finally, according to China, Yunan Province, Yuxi city land data, some simulation examples have been exploited to illustrate the applicability and usefulness of the developed theoretical result.



---

*Source: 1053283-2020-07-25.xml*