# Chaos Prediction in Fractional Delayed Energy-Based Models of Capital Accumulation

Authors: Mohamed El-Borhamy; Tamer Medhat; Manal E. Ali
Journal: Complexity
 (2021)

---

## Abstract

This paper presents the nonlinear dynamic analysis of energy-based models arisen from the applied systems characterized by the energy transport in the presence of fractional order derivative and time delay. The studied model is the fractional version of Bianca-Ferrara-Dalgaard-Strulik (BFDS) model of economy which is viewed as a transport network for energy in which the law of motion of capital occurs. By considering the time delay as bifurcation parameter, a proof to investigate the existence of Hopf bifurcation and the phase lock solutions using the Poincare-Linstedt and the harmonic balance methods is given. At definite values of time delay, period-doubling bifurcations followed up by the consequences of chaotic states are detected. Simulation results assure that the BFDS model can generate new (hyper) chaotic attractors beyond half order derivatives through the effect of the time delay on that system. In accordance with the literatures related to the problem of chaos, the concluded results confirm the proposed theorem by El-Borhamy in which the time delay possesses the ability to change the dynamic state of nonlinear systems from regular to chaotic within the fractional order derivative domain.


---

## Body

## 1. Introduction

Nowadays, the mathematical modeling using the fractional calculus, leading to the fractional differential equations (FDEs) or systems, becomes a new and hot challenge topic due to their frequent appearance in various physical and engineering applications, cf. [1, 2]. Some of the semianalytical methods are developed for solving fractional problems such as the adomian decomposition method, the homotopy perturbation method, the homotopy analysis method, and the perturbation techniques, cf. [3–6]. Another powerful method which can also give an explicit form for the solution is the variational iteration method (VIM), cf. [7, 8]. Consequently, considerable attention has been given to the qualitative and numerical analysis of systems with fractional order of the physical interest, because it is difficult to find exact solutions, so that a necessary object for several investigations of any scientific study in physical engineering and economical disciplines in this stream are focused on the stability, bifurcation, and chaos analysis, cf. [9–12]. The fractional order systems are characterized by two important mathematical features. First is the existence of an unlimited memory due to the fractional order while the integer order systems have a limited memory. The importance of this characterization roughly generalizes the mathematical treatment of the studied systems. Secondly, the existence of regular oscillations in fractional order systems requires very restricted conditions on the used fractional operator or the sustained periodic motion is not easy to be detected as well, cf. [13, 14]. In contrast, integer order systems possess the existence of regular oscillations which are periodic solutions or isolated periodic solutions (limit cycles).Retarded functional differential equations form a class of equations characterized by the existence of the time delay such that their solutions are influenced by the history. Such retarded equations generate infinite-dimensional dynamical systems, and they can be classified into various types of mathematical problems with delay (there are point delays, distributed delays, state-dependent delays, and integrals within or taken over the delay), cf. [15, 16].It is generally accepted that the presence of delays in systems is a potent source of nonstationary phenomena such as periodic oscillations, instabilities, and chaos, cf. [17, 18]. This can manifest itself as the loss of stability of an otherwise stable-steady state if the delay exceeds a certain threshold related to the dominant time-scale of a system. However, also there exists an evidence for that the time delay enhanced the stability, and short delays were shown to stabilize otherwise unstable dynamical systems. Also, it has been suggested that delays can damp out oscillations; this was shown with models for coupled oscillators under the condition that the delays in mutual interactions exceed a threshold value, cf. [19]. Indeed, in fractional systems with a history, the existence of time delay can support the system to create a sustained periodic output with regular modes.The study of stability and routes to chaos has become crucial for delayed fractional systems which has been taken care by many researchers, cf. [20–22]. Particularly, the study of bifurcation and chaos is an interesting and frustrating phenomenon and it is studied and developed intensively since founding the first chaotic attractor by Lorenz in [23]. In general, any (hyper) chaotic system is a nonlinear deterministic system that displays a complex and unpredictable behavior. Several methods are used to measure the existence of chaotic attractors, for instance, calculation of Lyapunov exponents, bifurcation diagrams, phase plane diagrams, and Poincare’s maps. Somehow, Lyapunov exponents are considered an excellent tool to decide whether the attractor is (hyper) chaotic or not, cf. [24, 25]. Meaningly, if the system is chaotic, then at least one of the Lyapunov exponents is positive and in the case of a hyperchaotic system, two or more Lyapunov exponents are positive.The contribution of this work is centered on to generate a new chaotic attractor from delayed fractional model of capital accumulation for the considered BFDS model in [26]. The study will be based on the influence of the fractional order and the time delay parameters on the chaotic dynamic state. Basic properties of the new chaotic model will be analyzed by means of bifurcation diagrams, phase plane diagrams, and Lyapunov exponents.The present work is organized as follows:Section 2 presents the fractional version of BFDS model. Section 3 handles the linear stability analysis of the fractional model and the characteristic function. Section 4 introduces Hopf bifurcation analysis for the system. Section 5 presents existence of the limit cycle using the perturbation analysis of the system. Section 6 presents the influence of time delay on the chaotic state of the system and the fractional order derivative as well. In the last, the conclusion is given.


## 2. Fractional Delayed BFDS Model

By building on the contribution by Bianca-Ferrara and Dalgaard-Strulik in [26–28] for the economy model, this work will develop the mathematical modeling for an economy growth which is viewed as a transport network for energy by including the fractional order derivative. In general, the energy conservation principles of dynamical systems have been proposed for the modeling of the law of motion of certain grown phenomena in physical sciences, cf. [28–33], so that the law of motion of capital occurs with a fractional time rate of energy and exhibiting whenever the system behavior is dependent at least in part on its history. The principle of energy conservation implies(1)Et−τ=μkt−τ+νDtqkt+nkt,q∈0,1,where μ is the energy requirement to operate and maintain the generic capital good, ν is the energy costs to create a new capital good, kt denotes capital stock assuming that the population grows at a constant rate of n, and Dtqkt represents Riemann–Liouville fractional derivative of kt. Specifically, the generalization of the BFDS model by including the energy conservation equation with a fractional derivative and time delay takes care of the previous occurring dynamics.Due to Dalgaard and Strulik mathematical model in [28], the energy is derived as(2)Et=εkt−τa,where 0<a<1 is a real constant, proportional to the dimension and efficiency of the network,(3)a=DD+η,where D is the dimension of the network, η∈1,D depends on the efficiency of the network, and ε>0 is a real constant in the sense that it is independent of capital per worker. In specific words, the base ground of the scaling coefficient, a, should fall in 1/2,D/D+1 and for a three-dimensional network, the scaling exponent a should fall in the interval 1/2,3/4, depending on the efficiency of the system. But, in this work, the domain of the scaling coefficient a has been extended within the interval 0,1.Consequently, the law of motion for capital is described by the following general nonlinear fractional delay differential equation:(4)Dtotqkt+nkt=rεkt−τa−μkt−τb,r=1ν,b=1,for some initial functions kt=ϕt∀t∈−τ,0, where the ratio μ/ν captures the physical phenomenon of capital depreciation.The common definitions concerning the fractional derivative are usually being Gru¨nwald–Letnikov fractional derivative, Riemann–Liouville fractional derivative, and Caputo’s fractional derivative which are mostly used to express the fractionalized dynamical systems, cf. [1, 2]. For a function kt defined on an interval to,T, the Riemann−Liouville fractional integral of kt of order q>0 is defined by(5)IRLtoqkt=1Γq∫tott−ξq−1kξdξ,t≥to,and Riemann–Liouville fractional derivative of kt of order q>0 is defined by(6)DtotqRLkt=dℓdtℓItoℓ−qkt,where ℓ−1<q<ℓ(i.e., ℓ is the integer order) and Γq is the gamma function. The definition of Caputo’s fractional derivative of order q>0 is(7)DtoCtqkt=Itoℓ−qkℓt.The relation between Caputo and Riemann–Liouville fractional derivatives is illustrated by the following expression:(8)DtoRLtqkt=DtoCtqkt+∑j=1ℓ−1kjtoΓj−q+1t−toj−q.Moreover, Caputo’s fractional derivative ofkt=c (constant) is always 0, whereas in the case of zero value of to=0, Riemann–Liouville fractional derivative of kt=c is different from zero as(9)DRLtqc=cΓ1−qt−q.However, by settingto⟶−∞ in both definitions and requiring reasonable behavior of kt and its derivatives, we end up with the same formula for both of them as the following:(10)D−∞RLtqkt=D−∞Ctqkt=D−∞tqkt=1Γℓ−q∫−∞tkℓξt−ξq−ℓ+1dξ.The fractional derivatives of sine or cosine functions in the Riemann–Liouville sense have no periodic behavior and they are(11a)D0RLtqcosωt=ωt−qE2,1−q−ω2t2,(11b)D0RLtqsinωt=ωt1−qE2,2−q−ω2t2,where Eq1,q2t is the two-parameter Mittag-Leffler function,(12)Eq1,q2t=∑j=0∞tjΓq1k+q2,q1>0,q2>0.When the time tends to infinity, it is corresponding to Weyl’s definition of fractional derivative, and then equation (11) converges to the following periodic functions:(13a)D0tqcosωt=ωqcosωt−πq2,t⟶∞,(13b)D0tqsinωt=ωqsinωt−πq2,t⟶∞.Also, by settingto⟶−∞ in equation (6), it is corresponding to Liouville’s definition of fractional derivative, and it also gives the following periodic results:(14a)D−∞tqcosωt=ωqcosωt+πq2,(14b)D−∞tqsinωt=ωqsinωt+πq2,(14c)D−∞tqexpωt=ωqexpωt.


## 3. Linear Stability Analysis

### 3.1. The Linearized Form

Equilibria or the stationary points, of course, coincide with the corresponding points for the zero fractional rate and zero delay, i.e.,Dtqkt=0, at τ=0(corresponding to the classical case when q⟶1, Dtqkt⟶k′t=0). Hence, there exists a unique positive steady state k∗, satisfying the relation εk∗a−1=μ+nν. After setting the following translation kt=xt+k∗, one can obtain(15)Dtotqxt+k∗+nxt+k∗=rεxt−τ+k∗a−μxt−τ+k∗.Hence, Taylor’s expansion at the zero equilibrium reads(16)Dtotqxt+nxt=εaνk∗a−1−μνxt−τ+12νaa−1k∗a−2εx2t−τ+16νaa−1a−2k∗a−3εx3t−τ+⋯,and its linearized form is(17)Dtotqxt=−nxt+kp+nxt−τ,where kp=a−1μ/ν<0. Equation (17) can be reduced to(18)Dtotqxt=−nxt+k˜pxt−τ,where k˜p=kp+n.


### 3.2. Particular Solutions of the Linearized Form

The prototype linear rate equation without delay and fractional derivative for a rate constantkp will have the following decaying solution:(19)xt=xoekpt−to,ifq=1,τ=0,n=0,where xo is an arbitrary constant.In fact, on the other hand, the same equation with a constant delay and integer order derivative will have multiple solutions ranging from monotonic decay to damped oscillations, stable periodic oscillations, and undamped growing oscillations based on the range of time lag values. For instance, a particular solution for the abovementioned linear rate equation with constant delay and integer order derivative at the zero rate of population can be obtained as(20)xt=xosinkpτ1−coskpτcoskpt−to−xosinkpt−to,ifq=1,n=0.Some methods are developed to obtain the solution for the delay differential equations (DDEs) for a given set of initial values, where the functionxt is defined over the interval −τ,0 such as the method of steps and the extended one-step techniques, cf. [34, 35].To follow-up, the same equation, without a constant delay but with a fractional derivative in Caputo’s sense, will have the following decaying solution:(21)xt=xoEq,1kpt−toq,whenτ=0,n=0.The general solution of a linear fractional delay problemn=0, seek in Liouville’s sense, is(22)xt=xoekpt−to,whenτ=1−qkplnkp≥0.


### 3.3. The Characteristic Function and Eigenvalues

Using the fundamental tools of fractional calculus, the characteristic function of the linear fractional time delay case is(23)Λλ=k˜pe−λτ−λq−n.It is easy to see that,(24)ifΛ0=k˜p−n<0,Λ1τ=k˜pe−1τq−n≥0,Then, one can find that there is at least oneλ<0, satisfying Λλ=0. Using such value of λ, the solution x˜t=xoe−λt becomes a positive solution with(25)xt>x˜ton−τ,0,andx0>xo,so that(26)xt>xoeλt,∀t≥0,as long as xt remains positive. But xt could become zero without first reaching x˜t, hence(27)xt>xoekpt∀t≥0,n=0,under the condition(28)kp≥1τqe,n=0.By considering equation (23),(29)Λλ=0⟶λq=−n+k˜pe−λτ.Then, using Lagrange’s expansion, letz be that root of(30)Fz=Fzo+ξϕz,which has the value Fz=Fzo when ξ=0. If Fz and ϕz are analytic functions inside and on a circular contour C containing z=zo, then the roots satisfy(31)Fz=Fzo+∑m=1∞ξmm!dm−1dsm−1F′sϕsms=zo.Hence, if we consider that(32)zo=−n1/q,Fλ=λq,ϕλ=e−λτ,then(33)Fλ=Fzo+∑m=1∞k˜pmm!dm−1dsm−1F′sϕsms=zo,so that one obtains(34)λq=zoq+q∑m=1∞k˜pmm!dm−1dsm−1sq−1e−msτs=zo.Then,(35)λqq=−n+k˜pzoq−1e−zoτ+12k˜p2q−1zoq−2−2τzoq−1e−2zoτ+⋯.Then, the eigenvalues can be represented by(36)λ=q−n+k˜pzoq−1e−zoτ+12k˜p2q−1zoq−2−2τzoq−1e−2zoτ+….1/q.Regarding the zero rate of population casen=0 in equation (29), one can get the following eigenvalues using the following expansion:(37)1kpλq=1−λτ1!+λ2τ22!−λ3τ33!+⋯.If0<τ≪1, we have the following linearized part of equation (37):(38)λ=1τ−1τkpλq.Then, the root of the characteristic function is approximated by(39)λ=1τ−qτqkp+2q2q−12!1τ2qkp2−3q3q−13q−23!1τ3qkp3+⋯.For the case ofτqkp=−1, the eigenvalues of equation (39) become(40)λ=1τ+q1!+2q2q−12!+3q3q−13q−23!+⋯.In equation (39), in the case of integer derivative q=1, it can be reduced to(41)λ=1τ+∑m=1∞−1m1τmqkpm.If we consider in equation (29) that q = 1, then the eigenvalue is given by(42)λ=−n−1τ∑m=1∞mm−1emnτ−τk˜pmm!,if the following condition holds:(43)k˜p<1τe−nτ+1.By considering thatλ=α+iω and to obtain the transition values α=0, we assume that λ=iωc, and one obtains(44)λ=ωce±π/2i=ωccosπ2+isin±π2,which is a root of equation (23) and ωc>0 is a real number.Then, one gets(45)ωcqcosqπ2+isin±qπ2+n−k˜pcosωcτc−isin±ωcτc=0.Separate the real and imaginary parts to obtain(46)ωcqcosqπ2+n−k˜pcosωcτc=0,ωcqsinqπ2−k˜psinωcτc=0.From squaring and adding the above two equations, one obtains(47)ωc2q+2nωcqcosπ2q−2k˜pωcqcosqπ2+ωcτc+k˜p2+n2=0.Thus, one can get the following constraint:(48)qπ2±−kp1/qτc≠2m+1π2,m∈ℤ,ifn=0,which means that the characteristic equation has no purely imaginary roots. Consequently, we can state the following theorem.Theorem 1.
Ifqπ/2±−kp1/qτc≠2m+1π/2, where m∈ℤ, then the zero solution of the linearized fractional time delay equation is globally asymptotically stable at zero rate of population.Figures1 and 2 illustrate the validity of Theorem 1 at certain value kp and q.Figure 1
Stability diagrams atτ = 1 when q = 0.25, 0.5, and 0.75 from left to right, respectively (n=m=0).Figure 2
Instability diagrams atτ = 1 when q = 0.25, 0.5, and 0.75 from left to right, respectively (n=0).


### 3.4. Local Stability Analysis of the Algebraic Retarded Equation

The algebraic retarded equation is obtained whenq⟶0 as a nondifferentiable version for equation (4) as follows:(49)kt=r1+nεkt−τa−μkt−τ.By considering that the equilibrium solution iske, it reads(50)ke=1+nν+με1/a−1,so that the equilibrium solution, represented by equation (50), is locally asymptotically stable, if all the roots λ satisfy the following condition:(51)λ<1,where(52)λ=rn+1a1+nν+a−1μ1/τ.


## 3.1. The Linearized Form

Equilibria or the stationary points, of course, coincide with the corresponding points for the zero fractional rate and zero delay, i.e.,Dtqkt=0, at τ=0(corresponding to the classical case when q⟶1, Dtqkt⟶k′t=0). Hence, there exists a unique positive steady state k∗, satisfying the relation εk∗a−1=μ+nν. After setting the following translation kt=xt+k∗, one can obtain(15)Dtotqxt+k∗+nxt+k∗=rεxt−τ+k∗a−μxt−τ+k∗.Hence, Taylor’s expansion at the zero equilibrium reads(16)Dtotqxt+nxt=εaνk∗a−1−μνxt−τ+12νaa−1k∗a−2εx2t−τ+16νaa−1a−2k∗a−3εx3t−τ+⋯,and its linearized form is(17)Dtotqxt=−nxt+kp+nxt−τ,where kp=a−1μ/ν<0. Equation (17) can be reduced to(18)Dtotqxt=−nxt+k˜pxt−τ,where k˜p=kp+n.


## 3.2. Particular Solutions of the Linearized Form

The prototype linear rate equation without delay and fractional derivative for a rate constantkp will have the following decaying solution:(19)xt=xoekpt−to,ifq=1,τ=0,n=0,where xo is an arbitrary constant.In fact, on the other hand, the same equation with a constant delay and integer order derivative will have multiple solutions ranging from monotonic decay to damped oscillations, stable periodic oscillations, and undamped growing oscillations based on the range of time lag values. For instance, a particular solution for the abovementioned linear rate equation with constant delay and integer order derivative at the zero rate of population can be obtained as(20)xt=xosinkpτ1−coskpτcoskpt−to−xosinkpt−to,ifq=1,n=0.Some methods are developed to obtain the solution for the delay differential equations (DDEs) for a given set of initial values, where the functionxt is defined over the interval −τ,0 such as the method of steps and the extended one-step techniques, cf. [34, 35].To follow-up, the same equation, without a constant delay but with a fractional derivative in Caputo’s sense, will have the following decaying solution:(21)xt=xoEq,1kpt−toq,whenτ=0,n=0.The general solution of a linear fractional delay problemn=0, seek in Liouville’s sense, is(22)xt=xoekpt−to,whenτ=1−qkplnkp≥0.


## 3.3. The Characteristic Function and Eigenvalues

Using the fundamental tools of fractional calculus, the characteristic function of the linear fractional time delay case is(23)Λλ=k˜pe−λτ−λq−n.It is easy to see that,(24)ifΛ0=k˜p−n<0,Λ1τ=k˜pe−1τq−n≥0,Then, one can find that there is at least oneλ<0, satisfying Λλ=0. Using such value of λ, the solution x˜t=xoe−λt becomes a positive solution with(25)xt>x˜ton−τ,0,andx0>xo,so that(26)xt>xoeλt,∀t≥0,as long as xt remains positive. But xt could become zero without first reaching x˜t, hence(27)xt>xoekpt∀t≥0,n=0,under the condition(28)kp≥1τqe,n=0.By considering equation (23),(29)Λλ=0⟶λq=−n+k˜pe−λτ.Then, using Lagrange’s expansion, letz be that root of(30)Fz=Fzo+ξϕz,which has the value Fz=Fzo when ξ=0. If Fz and ϕz are analytic functions inside and on a circular contour C containing z=zo, then the roots satisfy(31)Fz=Fzo+∑m=1∞ξmm!dm−1dsm−1F′sϕsms=zo.Hence, if we consider that(32)zo=−n1/q,Fλ=λq,ϕλ=e−λτ,then(33)Fλ=Fzo+∑m=1∞k˜pmm!dm−1dsm−1F′sϕsms=zo,so that one obtains(34)λq=zoq+q∑m=1∞k˜pmm!dm−1dsm−1sq−1e−msτs=zo.Then,(35)λqq=−n+k˜pzoq−1e−zoτ+12k˜p2q−1zoq−2−2τzoq−1e−2zoτ+⋯.Then, the eigenvalues can be represented by(36)λ=q−n+k˜pzoq−1e−zoτ+12k˜p2q−1zoq−2−2τzoq−1e−2zoτ+….1/q.Regarding the zero rate of population casen=0 in equation (29), one can get the following eigenvalues using the following expansion:(37)1kpλq=1−λτ1!+λ2τ22!−λ3τ33!+⋯.If0<τ≪1, we have the following linearized part of equation (37):(38)λ=1τ−1τkpλq.Then, the root of the characteristic function is approximated by(39)λ=1τ−qτqkp+2q2q−12!1τ2qkp2−3q3q−13q−23!1τ3qkp3+⋯.For the case ofτqkp=−1, the eigenvalues of equation (39) become(40)λ=1τ+q1!+2q2q−12!+3q3q−13q−23!+⋯.In equation (39), in the case of integer derivative q=1, it can be reduced to(41)λ=1τ+∑m=1∞−1m1τmqkpm.If we consider in equation (29) that q = 1, then the eigenvalue is given by(42)λ=−n−1τ∑m=1∞mm−1emnτ−τk˜pmm!,if the following condition holds:(43)k˜p<1τe−nτ+1.By considering thatλ=α+iω and to obtain the transition values α=0, we assume that λ=iωc, and one obtains(44)λ=ωce±π/2i=ωccosπ2+isin±π2,which is a root of equation (23) and ωc>0 is a real number.Then, one gets(45)ωcqcosqπ2+isin±qπ2+n−k˜pcosωcτc−isin±ωcτc=0.Separate the real and imaginary parts to obtain(46)ωcqcosqπ2+n−k˜pcosωcτc=0,ωcqsinqπ2−k˜psinωcτc=0.From squaring and adding the above two equations, one obtains(47)ωc2q+2nωcqcosπ2q−2k˜pωcqcosqπ2+ωcτc+k˜p2+n2=0.Thus, one can get the following constraint:(48)qπ2±−kp1/qτc≠2m+1π2,m∈ℤ,ifn=0,which means that the characteristic equation has no purely imaginary roots. Consequently, we can state the following theorem.Theorem 1.
Ifqπ/2±−kp1/qτc≠2m+1π/2, where m∈ℤ, then the zero solution of the linearized fractional time delay equation is globally asymptotically stable at zero rate of population.Figures1 and 2 illustrate the validity of Theorem 1 at certain value kp and q.Figure 1
Stability diagrams atτ = 1 when q = 0.25, 0.5, and 0.75 from left to right, respectively (n=m=0).Figure 2
Instability diagrams atτ = 1 when q = 0.25, 0.5, and 0.75 from left to right, respectively (n=0).


## 3.4. Local Stability Analysis of the Algebraic Retarded Equation

The algebraic retarded equation is obtained whenq⟶0 as a nondifferentiable version for equation (4) as follows:(49)kt=r1+nεkt−τa−μkt−τ.By considering that the equilibrium solution iske, it reads(50)ke=1+nν+με1/a−1,so that the equilibrium solution, represented by equation (50), is locally asymptotically stable, if all the roots λ satisfy the following condition:(51)λ<1,where(52)λ=rn+1a1+nν+a−1μ1/τ.


## 4. Hopf Bifurcation Analysis

Now, we are going to derive an expression for the critical time delay that separates stable and unstable regions bifurcated fromk∗, the unique positive equilibrium for the mathematical model. Moreover, we are going to prove a theorem which characterizes the nature of the equilibrium point k∗.Letλj=αjτ+iωjτ in general, so at αjτc=0 and ωjτc=ωc by substituting in the characteristic equation, then(53)ωc2q+2ncosπq2ωcq+n2−kp2=0.Then, one obtains(54a)ωcq=−ncosπq2±kp2−n2sin2π2q,(54b)τc=1ωcarctanπj−ωpqsinπq/2ωpqcosπq/2+n,j=0,1,2,….Thus, atn=0, one easily obtains(55)ωc=±kp1/qandτc=qωcπq2+2πj,j=0,1,2,….Therefore, to characterize the existence of the transition curves, we state the following theorem.Theorem 2.
Letk∗ be the unique positive equilibrium point for the mathematical model. Then, there exists a positive number τc such that the equilibrium k∗ is asymptotically stable for τ∈0,τc and unstable for τ>τc.Proof.
Letλj=αjτ+iωjτ denote a root of equation (23) near τ=τc, satisfying αjτc=0 and ωjτ=ωc. Differentiating the characteristic equation (equation (23)) with respect to τ, obtaining the change of sign of the real part α=ℛλ, one obtains(56)signℛdλτdτ−1|τ=τc=signqωc2nωcqcosπ/2q+ωc2qn2+2nωcqcosπ/2q+ωc2q,
Atn=0, this implies that(57)signℛdλτdτ−1|τ=τc=signqωc2,and at q=1, this gives(58)signℛdλτdτ−1|τ=τc=sign1n2+ωc2.
Thus, from equation (56), if we have that dαjτ/dτ|τ=τc>0, this implies that all the roots crossing the imaginary axis at iωc from left to right as τ increase and this result yields in the loss of stability. Then, the conclusion holds.


## 5. Construction of Periodic Solution

In this section, we investigate the existence of periodic solution without periodic input; otherwise, it is generated by the existence of time delay. Necessary conditions for generating periodic solutions, and the stability of periodic solution of equation (4) at τc are given by using the method based on Poincare-Lindstedt method and harmonic balance method, cf. [5, 36, 37].Following the procedures in [17], we start by considering Taylor’s expansion of equation (4) up to the third order at the zero equilibrium using Liouville’s sense of fractional derivative as obtained by equation (16). Then, we begin by introducing a small parameter ε via the scaling(59)xt,ε=∑i=0∞εi+1xit.Then, for values of delayτ which lie close to τc, expansion series is introduced as follows:(60)τε=τc+ετ1+ε2τ2+⋯.Next, we stretch time by replacing the independent variablet by s, as follows:(61)s=ωt,anddqdtq=ωqdqdsq.After substitution, we get(62)ωqDsqxs=−nxs+bxs−ωτ+a3x2s−ωτ+b4x3s−ωτ+⋯,where b=εa/νk∗a−1−μ/ν,a3=1/2νaa−1k∗a−2ε,b4=1/6νaa−1a−2k∗a−3ε. Expansions of ω,  ωq,  xs , and xs−ωτ in a power series in ε,(63)ω=ωc+εω1+ε2ω2+⋯.By recognizing thatω1=0 and τ1=0, then we obtain(64)xs−ωτ=xs−ωcτc−ε2ωcτ2+τcω2x′s−ωcτc.Finally, substituting and collecting terms, we obtain(65)Orderofε1:ωcqDsqx0s+nx0s−bx0s−ωcτc=0,Orderofε2:ωcqDsqx1s+nx1s−bx1s−ωcτc=a3x02s−ωcτc,Orderofε3:ωcqDsqx2s+nx2s−bx2s−ωcτc=b4x03s−ωcτc+2a3x0s−ωcτcx1s−ωcτc−bωcτ2+τcω2x0′s−ωcτc−qω2Dsqx0s.We take the solution of thex0s equation as follows:(66)x0s=A0coss.Then, we obtain an expression forx1s:(67)x1s=A1sin2s+B1cos2s+C1.Next, step substitute the expressions ofx0 and x1 into x2 equation after long simplifications; the coefficients of the resonant terms, sins and coss, are equated to zero. This results in Ao2 as follows:(68)Ao2=2−n+bbωcτ2+ω2τcsinωcτc−qω2cosπ/2q3/2a+bb4cosωcτc−2a32cosωcτc,or,Ao2=2a+bbωcτ2+ω2τccosωcτc+qω2sinπ/2q2a3sinωcτc−3/2a+bb4sinωcτc.By including the zero rate of populationn=0, one can obtain(69)A02=MNτ2,where(70)M=8kpkp222qωc2q+ωc4q−kp2q+1ωc3qcosqπ2kpωcsinωcτckpτccosωcτc+qωcq−sinqπ2−kpωccosωcτckpτcsinωcτc+qωcq−cosqπ2,N=−2kp2q−1ωc3qBsinqπ2kpτccosωcτc+qωcq−1sinqπ2−4Bsinqπ2+kpτcsinωcτc+qωpq−1cosqπ24Bcosqπ2+2qωc3qkpBcosqπ2−Bωc4qkpτccosωcτc+qωcq−1sinqπ24Bcosωcτc+kpτcsinωcτc+qωcq−1cosqπ24Bsinωcτc−Bkp222qωc2q+ωc4q−kp2q+1ωc3qcosqπ2kpτccosωcτc+qωcq−1sinqπ28Bcosωcτc+kpτcsinωcτc+qωcq−1cosqπ28Bsinωcτc+6kpCkp222qωc2q+ωc4q−kp2q+1ωc3qcosqπ2kpτccosωcτc+qωcq−1sinqπ2cosωcτc+kpτcsinωcτc+qωcq−1cosqπ2sinωcτc,where(71)B=aa−1μk∗−12ν,C=aa−1a−2μk∗−26ν.Thus,A0 is real so A02>0, this means that τ2 should have the same sign as M/N. If q=1, this leads to the classical model with(72)A02=80k∗2ωc6π−4a2+10πaτ2.In order forA02 to be positive, τ2 is required to be more than zero, which means that the limit cycle is born out of unstable fixed point. Since the stability of the limit cycle should be the opposite of the stability of fixed point from which it is born, our conclusion is that the limit cycle is stable and we have supercritical Hopf bifurcation, so that the following theorem is concluded.Theorem 3.
Forq∈0,1, a stable limit cycle is born out of unstable fixed point with radius A0.Numerically, these results are verified by Figures3–6 to confirm this conclusion at the values a=0.001, ε=1, μ=1, and ν=2. The economy oscillator settles down to a steady oscillation with uniform periods at each fractional order derivative. This is clearly appeared in the phase plane for t∈900,1000 by avoiding the transients which have completely disappeared at the scale shown.Figure 3
The birth of Hopf-bifurcated limit cycle from the fixed pointk∗ for the studied case when q=0.25.Figure 4
The birth of Hopf-bifurcated limit cycle from the fixed pointk∗ for the studied case when q=0.5.Figure 5
The birth of Hopf-bifurcated limit cycle from the fixed pointk∗ for the studied case when q=0.75.Figure 6
The birth of Hopf-bifurcated limit cycle from the fixed pointk∗ for the studied case when q=1.


## 6. Routes to Chaos

Dynamic trajectories can be drawn to check the existence of aperiodic or chaotic motion in BFDS model system. More or less most of the analytical methods cannot provide enough information to specify the perfect existence of chaotic motion. Therefore, methods such as bifurcation diagrams, Lyapunov exponent diagrams, and Poincare’s maps or phase trajectory diagrams can identify the behavior of the system motion through given parameters. This study will consider the most sensitive parameter produced chaos, which is the time delay, cf [17]. In most cases, studying the effect of the time delay to create chaos in delayed systems is indispensable and more important due to its inherence with them.In our analysis, started by using the bifurcation diagrams, we obtain valuable insights into system’s nonlinear analysis, since with a range of such parameter variation. Some what, the system dynamic is mostly to be observed under the influence of that parameter. Bifurcation diagrams are drawn under the variation of time delay with constant interval. The bifurcation can be easily detected by graph the relationship between the solutionkt and the time delay under the step size Δτ=0.001 for the shown figures in Figure 7. It can be vividly observed that its existence affects the characteristic of the dynamic trajectories. Roughly speaking, at different values of fractional orders q, the system has period-doubling behavior according to the fixed parameters, n=1, ε=15, μ=0.95, and ν=0.29, for a=0.125.Figure 7
Bifurcation diagrams ata = 0.125 and q = 0.5, 0.75, and 1, respectively.The figures started withq=0.5, q=0.75, and q=1, respectively. It is obvious from each figure that, with further increasing of the time delay, the model response gradually started by a stable response up to a certain threshold and after that it enters a chaotic zone by the route of period-doubling bifurcation supported by the value of fractional derivative. This indicates that it might be predicted a strange attractor, representing the chaotic and irregular motion in the model as shown in the figures:Figures8–10 at a=0.125, q=0.25, and τ=1, 2.5, and 5, respectively. Figures 11–13 at a=0.125, q=0.5, and τ=1, 2.5, and 5, respectively. Figures 14–16 at a=0.125, q=0.75, and τ=1, 1.5, and 2, respectively. Figures 17–19 at a=0.125, q=1, and τ=1, 1.3, and 1.45, respectively.Figure 8
Solutions and phase trajectories ata = 0.125, q = 0.25, and τ = 1.Figure 9
Solutions and phase trajectories ata = 0.125, q = 0.25, and τ = 2.5.Figure 10
Solutions and phase trajectories ata = 0.125, q = 0.25, and τ = 5.Figure 11
Solutions and phase trajectories ata = 0.125, q = 0.5, and τ = 1.Figure 12
Solutions and phase trajectories ata = 0.125, q = 0.5, and τ = 2.5.Figure 13
Solutions and phase trajectories ata = 0.125, q = 0.5, and τ = 5.Figure 14
Solutions and phase trajectories ata = 0.125, q = 0.75, and τ = 1.Figure 15
Solutions and phase trajectories ata = 0.125, q = 0.75, and τ = 1.5.Figure 16
Solutions and phase trajectories ata = 0.125, q = 0.75, and τ = 2.Figure 17
Solutions and phase trajectories ata = 0.125, q = 1, and τ = 1.Figure 18
Solutions and phase trajectories ata = 0.125, q = 1, and τ = 1.3.Figure 19
Solutions and phase trajectories ata = 0.125, q = 1, and τ = 1.45.Indeed, these results assure that the prediction of the chaos is not far away to be obtained within some specific values of the model parameters taking into the account that the branch points have to be changed with further increasing of the parametera. Moreover, this model with q∈0,1 can generate hyperchaotic behavior due to the two positive values of Lyapunov exponents as shown in Figure 20, so that the obtained figures indicate the generation of new strange attractor for fractional order equation with delay based on the three measurements (bifurcation diagram, phase trajectories, and Lyapunov exponents). These results confirm the results in [17] by the responsibility of the time delay to create chaotic motion in delayed system with the support of the fractional order derivative.Figure 20
Lyapunov exponents ata=0.125 and τ=15 for q=0.5, q=0.75, and q=1, respectively.


## 7. Conclusion

This work has studied the dynamic change of a single nonlinear equation represented by the BFDS model with delay supported by the fractional order derivative within the interval (0, 1]. The dynamic of the system fluctuates due to the time delay parameter to generate stable periodic, aperiodic, or chaotic behavior which has been evaluated based on three measurements: bifurcation diagrams, phase diagrams, and Lyapunov exponents. Simulation results showed that the BFDS model can generate period-doubling bifurcation, leading to various new chaotic attractors based on the time delay parameter and the fractional order derivative. Comparatively, the obtained figures indicate that the generation of new hyperchaotic strange attractor is mostly predictable for every value of the fractional derivative of the model. Consequently, based on the results inferred from the BFDS model, it has been confirmed that the time delay is able to change the dynamic state of nonlinear systems from regular to chaotic which is in agreement with El-Borhamy’s theorem in [17] within the fractional order derivative domain q∈0,1.



---

*Source: 8751963-2021-03-09.xml*