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

#### 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 implieswhere is the energy requirement to operate and maintain the generic capital good, is the energy costs to create a new capital good, denotes capital stock assuming that the population grows at a constant rate of , and represents Riemann–Liouville fractional derivative of . 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 aswhere is a real constant, proportional to the dimension and efficiency of the network,where is the dimension of the network, depends on the efficiency of the network, and 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, , should fall in and for a three-dimensional network, the scaling exponent should fall in the interval , depending on the efficiency of the system. But, in this work, the domain of the scaling coefficient has been extended within the interval .

Consequently, the law of motion for capital is described by the following general nonlinear fractional delay differential equation:for some initial functions , where the ratio captures the physical phenomenon of capital depreciation.

The common definitions concerning the fractional derivative are usually being Grnwald–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 defined on an interval , the RiemannLiouville fractional integral of of order is defined byand Riemann–Liouville fractional derivative of of order is defined bywhere (i.e., is the integer order) and is the gamma function. The definition of Caputo’s fractional derivative of order is

The relation between Caputo and Riemann–Liouville fractional derivatives is illustrated by the following expression:

Moreover, Caputo’s fractional derivative of (constant) is always 0, whereas in the case of zero value of , Riemann–Liouville fractional derivative of is different from zero as

However, by setting in both definitions and requiring reasonable behavior of and its derivatives, we end up with the same formula for both of them as the following:

The fractional derivatives of sine or cosine functions in the Riemann–Liouville sense have no periodic behavior and they arewhere is the two-parameter Mittag-Leffler function,

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:

Also, by setting in equation (6), it is corresponding to Liouville’s definition of fractional derivative, and it also gives the following periodic results:

#### 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., , at (corresponding to the classical case when , ). Hence, there exists a unique positive steady state , satisfying the relation . After setting the following translation , one can obtain

Hence, Taylor’s expansion at the zero equilibrium readsand its linearized form iswhere . Equation (17) can be reduced towhere .

##### 3.2. Particular Solutions of the Linearized Form

The prototype linear rate equation without delay and fractional derivative for a rate constant will have the following decaying solution:where 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

Some methods are developed to obtain the solution for the delay differential equations (DDEs) for a given set of initial values, where the function is defined over the interval 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:

The general solution of a linear fractional delay problem , seek in Liouville’s sense, is

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

It is easy to see that,

Then, one can find that there is at least one , satisfying . Using such value of , the solution becomes a positive solution withso thatas long as remains positive. But could become zero without first reaching , henceunder the condition

By considering equation (23),

Then, using Lagrange’s expansion, let be that root ofwhich has the value when . If and are analytic functions inside and on a circular contour containing , then the roots satisfy

Hence, if we consider thatthenso that one obtains

Then,

Then, the eigenvalues can be represented by

Regarding the zero rate of population case in equation (29), one can get the following eigenvalues using the following expansion:

If , we have the following linearized part of equation (37):

Then, the root of the characteristic function is approximated by

For the case of , the eigenvalues of equation (39) become

In equation (39), in the case of integer derivative , it can be reduced to

If we consider in equation (29) that *q* = 1, then the eigenvalue is given byif the following condition holds:

By considering that and to obtain the transition values , we assume that , and one obtainswhich is a root of equation (23) and is a real number.

Then, one gets

Separate the real and imaginary parts to obtain

From squaring and adding the above two equations, one obtains

Thus, one can get the following constraint:which means that the characteristic equation has no purely imaginary roots. Consequently, we can state the following theorem.

Theorem 1. *If , where , then the zero solution of the linearized fractional time delay equation is globally asymptotically stable at zero rate of population.*

Figures 1 and 2 illustrate the validity of Theorem 1 at certain value and *q*.

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

The algebraic retarded equation is obtained when as a nondifferentiable version for equation (4) as follows:

By considering that the equilibrium solution is , it readsso that the equilibrium solution, represented by equation (50), is locally asymptotically stable, if all the roots satisfy the following condition:where

#### 4. Hopf Bifurcation Analysis

Now, we are going to derive an expression for the critical time delay that separates stable and unstable regions bifurcated from , the unique positive equilibrium for the mathematical model. Moreover, we are going to prove a theorem which characterizes the nature of the equilibrium point .

Let in general, so at and by substituting in the characteristic equation, then

Then, one obtains

Thus, at , one easily obtains

Therefore, to characterize the existence of the transition curves, we state the following theorem.

Theorem 2. *Let be the unique positive equilibrium point for the mathematical model. Then, there exists a positive number such that the equilibrium is asymptotically stable for and unstable for .*

*Proof. *Let denote a root of equation (23) near , satisfying and . Differentiating the characteristic equation (equation (23)) with respect to , obtaining the change of sign of the real part , one obtainsAt , this implies thatand at , this givesThus, from equation (56), if we have that , this implies that all the roots crossing the imaginary axis at 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 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

Then, for values of delay which lie close to , expansion series is introduced as follows:

Next, we stretch time by replacing the independent variable by , as follows:

After substitution, we getwhere Expansions of , , , and in a power series in ,

By recognizing that and , then we obtain

Finally, substituting and collecting terms, we obtain

We take the solution of the equation as follows:

Then, we obtain an expression for :

Next, step substitute the expressions of and into equation after long simplifications; the coefficients of the resonant terms, and , are equated to zero. This results in as follows:

By including the zero rate of population , one can obtainwherewhere

Thus, is real so , this means that should have the same sign as . If , this leads to the classical model with

In order for to be positive, 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. *For , a stable limit cycle is born out of unstable fixed point with radius .*

Numerically, these results are verified by Figures 3–6 to confirm this conclusion at the values , , , and . 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 by avoiding the transients which have completely disappeared at the scale shown.

#### 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 solution and the time delay under the step size 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 , the system has period-doubling behavior according to the fixed parameters, , , , and , for .

The figures started with , , and , 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:

Figures 8–10 at , , and 1, 2.5, and 5, respectively. Figures 11–13 at , , and 1, 2.5, and 5, respectively. Figures 14–16 at , , and 1, 1.5, and 2, respectively. Figures 17–19 at , , and 1, 1.3, and 1.45, respectively.

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 parameter . Moreover, this model with 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.

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

#### Data Availability

The data used to support the findings of this study are available from the authors upon request.

#### Disclosure

This study was performed as part of the employment of the authors: University of Tanta and University of Kafrelsheikh, Egypt.

#### Conflicts of Interest

The authors declare no conflicts of interest.

#### Authors’ Contributions

All authors have read and agreed to the published version of the manuscript.