On a fractional-order delay Mackey-Glass equation

Abstract

In this paper, a fractional-order Mackey-Glass equation with constant delay is considered. The local stability of the fixed points is analyzed. Moreover, a discretization process is applied to convert the fractional-order delay equation to its discrete analog. A numerical simulation including Lyapunov exponent, phase diagrams, bifurcation, and chaos is carried out using Matlab to ensure theoretical results and to reveal more complex dynamics of the equation after discretization.

El-Sayed et al. Advances in Dierence Equations (2016) 2016:137 DOI 10.1186/s13662-016-0863-x

R E S E A R C H Open Access

http://crossmark.crossref.org/dialog/?doi=10.1186/s13662-016-0863-x&domain=pdf

Web End = On a fractional-order delay Mackey-Glass equation

Ahmed MA El-Sayed1*, Sanaa M Salman2 and Naemaa A Elabd3

*Correspondence: mailto:[email protected]

Web End [email protected]

1Department of Mathematics, Faculty of Science, Alexandria University, Alexandria, EgyptFull list of author information is available at the end of the article

Abstract

In this paper, a fractional-order Mackey-Glass equation with constant delay is considered. The local stability of the xed points is analyzed. Moreover,a discretization process is applied to convert the fractional-order delay equation to its discrete analog. A numerical simulation including Lyapunov exponent, phase diagrams, bifurcation, and chaos is carried out using Matlab to ensure theoretical results and to reveal more complex dynamics of the equation after discretization.

Keywords: fractional-order delay Mackey-Glass equations; xed points; local stability; discretization; Lyapunov exponent; bifurcation; chaos

1 Introduction

Delay dierential equations (DDEs) arise in the mathematical description of systems whose time evolution depends explicitly on a past state of the system, as for example in the case of delayed feedback. Neural systems [], respiration regulation [], agricultural commodity markets [], nonlinear optics, and neutrophil populations in the blood [] are but a few systems in which delayed feedback leads naturally to a description in terms of a delay dierential equation. We will restrict our attention to systems modeled by evolutionary delay equations that can be expressed in the form

x (t) = f x(t), x(t ) , x(t) n, t . (.)

Here the state of the system at time t is x(t), whose rate of change depends explicitly, via the function f , on the past state x(t ) where is a xed time delay. More general delay equations might be considered: multiple time delays, variable time delays, continuously distributed delays, and higher derivatives all arise in applications and lead to more complicated evolution equations. Nevertheless, equations of the form (.) constitute a suciently broad class of systems to be of practical importance, and they will provide adequate fodder for the types of problems we wish to consider.

DDEs arise in many areas of mathematical modeling: for example, population dynamics (taking into account the gestation times), infectious diseases (accounting for the incubation periods), physiological and pharmaceutical kinetics (modeling, for example, the bodys reaction to CO, etc. in circulating blood), chemical kinetics (such as mixing reactants), the navigational control of ships and aircraft (with, respectively, large and short lags), and more general control problems (see for example []).

2016 El-Sayed et al. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

El-Sayed et al. Advances in Dierence Equations (2016) 2016:137 Page 2 of 11

On the other hand, fractional calculus is a generalization of classical dierentiation and integration to arbitrary (non-integer) order []. Many mathematicians and applied researchers have tried to model real processes using fractional calculus []. In recent years dierential equations with fractional-order have attracted many researchers because of their applications in many areas of science and engineering. Analytical and numerical techniques have been implemented to study such equations. The fractional calculus has allowed the operations of integration and dierentiation to be applied for any fractional order [].

We recall the basic denitions (Caputo) and properties of fractional-order dierentiation and integration.

Denition The fractional integral of order

R+ of the function f (t), t > is dened

by

Iaf (t) =

t

(t s)

() f (s) ds,

and the fractional derivative of order (n , n) of f (t), t > is dened by

Daf (t) = InDnf (t), D = d

dt .

In addition, the following results are the main features in fractional calculus. Let ,

R+, (, ): Ia : L L, and if f (x) L, then Ia Iaf (x) = I+af (x). limn Iaf (x) = Inaf (x) uniformly on [a, b], n = , , , . . . , where Iaf (x) =

xa f (s) ds.

lim Iaf (x) = f (x) weakly. If f (x) is absolutely continuous on [a, b], then lim Daf (x) = df(x)dx.

The Mackey-Glass equation is a nonlinear time delay dierential equation, which was proposed as a model of hematopoiesis, given by

dx dt =

x + xc

x, (.)

where , c, , are real parameters, and x represents the value of the variable x at time (t ). Depending on the values of the parameters, this equation displays a range of periodic and chaotic dynamics.

In this work, we will show that considering a fractional-order derivative with delay in equation (.) will exhibit more complex and richer dynamics.

Consider the fractional-order delay Mackey-Glass equation given in the form

Dx(t) =

x(t ) + x(t )c x(t

), t (, T], (.)

with the initial condition

x() = x, (.)

where (, ], R+, and c > . In equation (.), we consider delay in the last term.

El-Sayed et al. Advances in Dierence Equations (2016) 2016:137 Page 3 of 11

2 Discretization process

In this part, we apply the discretization process represented in [, ], and [] for discretizing the delay fractional-order Mackey Glass equation with piecewise constant arguments given by

Dx(t) =

x([ tr ]r )

+ xc([ tr ]r )

x t

r

r , (.)

with initial condition (.).The steps of the discretization process are as follows. Let t [, r), then tr [, ). That is,

Dx(t) =

x() + xc() x(

), (.)

and the solution of (.) is given by

x(t) = x + I

x() + xc() x(

)

= x +

x() + xc() x(

)

t

(t s)

() ds

= x + t

( + )

x() + xc() x(

) .

Let t [r, r), then tr [, ). That is,

Dx(t) =

x(r )

+ xc(r ) x(r

), t [r, r), (.)

and the solution of (.) is given by

x(t) = x(r) + Ir

x(r )

+ xc(r ) x(r

)

= x(r) +

x(r )

+ xc(r ) x(r

)

t

(t s)

() ds

r

= x(r) + (t r)

( + )

x(r )

+ xc(r ) x(r

) .

Repeating the process we can easily get

x(t) = xn(nr) + (t nr)

( + )

xn(nr )

+ xcn(nr ) xn(nr

) , t nr, (n + )r .

Let t (n + )r, we obtain the discretization

xn+(n + ) = xn + (r)

( + )

xn(r )

+ xcn(r ) xn(r

) . (.)

It is worth to pay attention here that Eulers discretization method is an approximation for the derivative while the predictor-corrector method is an approximation for the integral.

El-Sayed et al. Advances in Dierence Equations (2016) 2016:137 Page 4 of 11

However, our proposed discretization method here is an approximation for the right-hand side of the system under consideration as is pretty clear from (.). Moreover, we have noticed that when , the discretization will be Eulers discretization.

In the following, we will discuss two cases of the delay: Case I: = r, and Case II: = r.

3 Case I: = r

In this case we have a second-order dierence equation given by

xn+ = xn + r

( + )

xn

+ xcn xn

. (.)

Existence and stability of xed points

To nd the xed points of system (.), we rst split it into two rst-order dierence equations as follows:

xn+ = xn + r

( + )

yn

+ ycn yn

,

(.)

yn+ = xn.

For all values of the parameter , system (.) has one xed point, namely, x = (, ). For > , we have an additional xed point, which is x = ( c , c ).

In order to study the local stability of these xed points, we need the moduli of the eigenvalues of the Jacobian matrix evaluated at each of the xed points []. The Jacobian matrix of system (.) evaluated at any xed point (x, y) is obtained by

J =

r

(+) (+y

c(c) (+yc) )

.

The eigenvalues associated to the Jacobian matrix are

, = . . + R(s ),

where

R = r

( + ), S =

+ yc( c)

( + yc) .

The xed points x, x of the system equation (.) are stable if |i| < , i = , . In order to study the qualitative behavior of the solution of system (.) we rely on the Jury criteria given generally by

. F := + T + D > ,. TC := T + D > ,. H := D > ,where the trace and determinant of the Jacobian matrix are given, respectively, by

T: Tr(J) = J + J,

D: Det(J) = JJ JJ.

El-Sayed et al. Advances in Dierence Equations (2016) 2016:137 Page 5 of 11

Proposition The xed point x is locally asymptotically stable if < ( + /R), and losses stability via a ip bifurcation when > and via a Neimark-Sacker bifurcation when > r

(+) r .

Proof The Jacobian matrix at the rst xed point x is obtained by

J =

r

(+) ( )

,

which has two eigenvalues,

, = . . + R( ).

According to the Jury criteria [, ], where T = , D = r

() (+) , the rst condition is always satised, while the second and third may be violated. That is, the xed point x loses stability via a ip bifurcation when > , and via a Neimark-Sacker bifurcation when >

r (+)r .

Proposition The xed point x of system (.) is stable if < cR

+cR , and it loses stability

via a pitchfork bifurcation if > cR

+cR , via a ip bifurcation if > , and via a Neimark-Sacker bifurcation if < cR

cR .

Proof Calculating the Jacobian matrix at the second xed point x of system (.) we obtain

J =

r

c() (+)

,

which has two eigenvalues,

, = . .

+ Rc(

) ,

where the trace and determinant of J(x) are given, respectively, by

T = , D = Rc

.

According to the Jury criteria, the three conditions may be all violated. That is, x loses stability via a pitchfork bifurcation if > cR

+cR , via a ip bifurcation if > , and via a Neimark-Sacker bifurcation if < cR

cR .

4 Case II: = 2r

In this section, we take the delay to be = r in equation (.). Applying the discretization process we end up with a system of third-order dierence equations given by

xn+ = xn + r

( + )

xn

+ xcn xn

. (.)

El-Sayed et al. Advances in Dierence Equations (2016) 2016:137 Page 6 of 11

To study the xed points of system (.) we rst split it into three rst-order dierence equations as follows:

xn+ = xn + r

( + )

yn+ = xn, (.)

zn+ = yn.

In the following, we study the local stability of the xed points of the system (.).

Existence and stability of xed points

System (.) has the following xed points: For all parameter values, there is only one xed point x x = (, , ). For > , there is an additional xed point x x = ( c , c , c ).

By considering a Jacobian matrix for one of these xed points and calculating their eigenvalues, we can investigate the stability of each xed point based on the roots of the system characteristic equation. The Jacobian matrix is given by

J =

(+) , S = +z

( + ) r .

The remaining conditions give the following inequalities: a < > R , |b| > b ( (R R)) > (R R), c > |c| ( (R R)) (R R) > (R R) + (R R).

Linearizing the system (.) about x x yields the following characteristic equation:

P() = cR(

) . (.)

yn

+ zcn zn

,

R(S )

,

c(c) (+zc) .

Linearizing the system (.) about x x yields the following characteristic equation:

P() = R( ). (.)

Let

a = , a = , a = , a = R( ).

From the Jury test, if P() > , P() < , and a < , |b| > b, c > |c|, where b = a, b = a aa, b = a aa, c = b b, and c = bb bb, then the roots of P() satisfy < and thus x is asymptotically stable.

The rst condition gives < , while the second condition gives

> r

where R = r

El-Sayed et al. Advances in Dierence Equations (2016) 2016:137 Page 7 of 11

Let

From the Jury test, if P() > , P() < , and a < , |b| > b, c > |c|, where b = a, b = a aa, b = a aa, c = b b, and c = bb bb, then the roots of P() satisfy

< and thus x is asymptotically stable.

We are going to check these conditions at x x:

p() > > ,

p() < <

cr cr (+) ,

|b| > b ( + cR ( )) > cR().

Thus, any condition may be violated resulting in instability of x x.

5 Numerical simulation

In this section, a numerical simulation is carried out with the aid of Matlab to illustrate our theoretical results and to reveal the more complex dynamics of equation (.) in the two cases = r and = r. In all numerical simulations, we take c = , and r = .. First of all, let us consider system (.). Indeed, if one is interested in determining whether a dynamical system is chaotic or not, often just a few of the largest Lyapunov characteristic exponents (LCEs) may provide the answer. This actually is so because a positive LCE is a good indicator for chaos. Since for non-chaotic systems all LCEs are non-positive, the presence of a positive LCE has often been used to help determine if a system is chaotic or not. In this paper, we compute the LCEs via the Householder QR-based methods described in []. For system (.), we get when = ., c = , and r = ., LCE = ., and LCE = . as shown in Figure . We vary the parameter and x the other parameters, c, r, and . Bifurcation diagrams of system (.) are also shown in Figure for dierent values of the fractional-order parameter . If we consider = ., it is shown that the xed point x(, ) is stable if < , and at > it losses stability via a ip bifurcation. Afterwards, a stable periodic solution of period appears, then the periodic solution of period becomes unstable, and a periodic solution of period appears and chaos happens. Figure shows the dierent phase plane for system (.) for = .. For = ., Figure (a) shows an invariant closed curve bifurcating from x(, ), while for = ., Figure (b) shows a chaotic attractor. Now we vary the parameter from to and x to plot the bifurcation diagram for system (.) as a function of as shown in Figure .

Next, we turn to the second case, when = r. Figure shows the bifurcation diagram for system (.) as a function of . If = ., the gure shows that the xed point x x becomes unstable when > as it loses stability via a ip bifurcation. Then the stable period- orbit appears at = ., which in turn loses stability; then chaos appears. Finally, Figure shows the phase plane for system (.) for = .. Figure (a) shows a stable xed point x x for = ., Figure (b) shows a double scroll for = ., and Figure (c),(d) shows chaotic attractors for = , and = ., respectively.

6 Conclusion

In this paper, the dynamic behavior of a fractional-order delay Mackey-Glass equation is investigated after applying a discretization process to it. We have considered two dierent

a < <

a = , a = , a = , a = cR(

) .

cr cr (+) ,

El-Sayed et al. Advances in Dierence Equations (2016) 2016:137 Page 8 of 11

(a) (b)

(c) (d)

(e) (f)

Figure 1 Bifurcation diagram and Lyapunov exponent for system (3.2) with different values of the fractional-order parameter , c = 6 and r = 0.5.

cases for the delay , the rst is when = r, and the second is when = r, where r is the discretization parameter. Stability of the xed points and local bifurcations of xed points of the discretized systems in the two cases was are analyzed. A numerical simulation was carried out to ensure our theoretical analysis and to reveal the more complex dynamics of the system.

El-Sayed et al. Advances in Dierence Equations (2016) 2016:137 Page 9 of 11

(a) (b)

Figure 2 Phase plane for system (3.2) with = 0.95, c = 6 and r = 0.5.

(a) (b)

Figure 3 Bifurcation diagrams for system (3.2) as a function of .

(a) (b)

Figure 4 Bifurcation diagram and chaos for system (4.2) as a function of .

El-Sayed et al. Advances in Dierence Equations (2016) 2016:137 Page 10 of 11

(a) (b)

(c) (d)

Figure 5 Phase plane for system (3.1) with = 0.95, c = 6, and r = 0.5.

Competing interests

The authors declare that they have no competing interests.

Authors contributions

The authors declare that the study was realized in collaboration with the same responsibility. All authors read and approved the nal manuscript.

Author details

1Department of Mathematics, Faculty of Science, Alexandria University, Alexandria, Egypt. 2Faculty of Education, Alexandria University, Alexandria, Egypt. 3Faculty of Science, Alexandria University, Alexandria, Egypt.

Acknowledgements

The authors would like to thank the referees of this manuscript for their valuable comments and suggestions.

Received: 10 February 2016 Accepted: 11 May 2016

References

1. an der Heiden, U: Delays in physiological systems. J. Math. Biol. 8, 345-364 (1979)2. Glass, L, Mackey, MC: Pathological conditions resulting from instabilities in physiological control systems. Ann. N.Y. Acad. Sci. 316, 214-235 (1979)

3. Mackey, MC: Commodity price uctuations: price dependent delays and non- linearities as explanatory factors.J. Econ. Theory 48(2), 497-509 (1989)4. Baleanu, D, Magin, RL, Bhalekar, S, Daftardar-Gejji, V: Chaos in the fractional order nonlinear Bloch equation with delay. Commun. Nonlinear Sci. Numer. Simul. 25(1-3), 41-49 (2015)

5. Wu, G-C, Baleanu, D: Discrete chaos in fractional delayed logistic maps. Nonlinear Dyn. 80(4), 1697-1703 (2015)6. Jarad, F, Abdeljawad, T, Baleanu, D: Higher order fractional variational optimal control problems with delayed arguments. Appl. Math. Comput. 218(18), 9234-9240 (2012)

7. Kilbas, AA, Srivastava, HM, Trujillo, JJ: Theory and Applications of Fractional Dierential Equations. Elsevier, Amsterdam (2006)

8. Podlubny, I: Fractional Dierential Equations. Mathematics in Science and Engineering, vol. 198. Academic Press, San Diego (1999)

El-Sayed et al. Advances in Dierence Equations (2016) 2016:137 Page 11 of 11

9. Miller, KS, Ross, B: An Introduction to the Fractional Calculus and Fractional Dierential Equation. Wiley, New York (1993)

10. Yang, XJ, Srivastava, HM, Cattani, C: Local fractional homotopy perturbation method for solving fractal partial dierential equations arising in mathematical physics. Rom. Rep. Phys. 67(3), 752-761 (2015)

11. Yang, XJ, Srivastava, HM: An asymptotic perturbation solution for a linear oscillator of free damped vibrations in fractal medium described by local fractional derivatives. Commun. Nonlinear Sci. Numer. Simul. 29(1), 499-504 (2015)

12. Yang, XJ, Baleanu, D, Srivastava, HM: Local fractional similarity solution for the diusion equation dened on Cantor sets. Appl. Math. Lett. 47, 54-60 (2015)

13. Wang, J, Ye, Y, Pan, X, Gao, X, Zhuang, C: Fractional zero-phase ltering based on the Riemann-Liouville integral. Signal Process. 98(5), 150-157 (2014)

14. Wang, J, Ye, Y, Gao, Y, Gao, X, Qian, S: Fractional compound integral with application to ECG signal denoising. Circuits Syst. Signal Process. 34, 1915-1930 (2015)

15. Wang, J, Ye, Y, Pan, X, Gao, X: Parallel-type fractional zero-phase ltering for ECG signal denoising. Biomed. Signal Process. Control 18, 36-41 (2015)

16. Wang, J, Ye, Y, Gao, X: Fractional 90 phase-shift ltering based on the double-sided Grunwald-Letnikov dierintegrator. IET Signal Process. 9(4), 328-334 (2015)

17. Bhrawy, AH, Alhamed, YA, Baleanu, D, Al-Zahrani, AA: New spectral techniques for systems of fractional dierential equations using fractional-order generalized Laguerre orthogonal functions. Fract. Calc. Appl. Anal. 17, 1137-1157 (2014)

18. Hafez, RM, Ezz-Eldien, SS, Bhrawy, AH, Ahmed, EA, Baleanu, D: A Jacobi Gauss-Lobatto and Gauss-Radau collocation algorithm for solving fractional Fokker-Planck equations. Nonlinear Dyn. 82(3), 1431-1440 (2015)

19. Bhrawy, AH, Hafez, RM, Alzahrani, E, Baleanu, D, Alzahrani, AA: Generalized Laguerre-Gauss-Radau scheme for the rst order hyperbolic equations in a semi-innite domain. Rom. J. Phys. 60, 918-934 (2015)

20. Bhrawy, AH, Abdelkawy, MA, Alzahrani, AA, Baleanu, D, Alzahrani, E: A Chebyshev-Laguerre Gauss-Radau collocation scheme for solving time fractional sub-diusion equation on a semi-innite domain. Proc. Rom. Acad., Ser. A 16, 490-498 (2015)

21. Doha, EH, Bhrawy, AH, Baleanu, D, Ezz-Eldien, SS, Hafez, RM: An ecient numerical scheme based on the shifted orthonormal Jacobi polynomials for solving fractional optimal control problems. Adv. Dier. Equ. 2015, 15 (2015)

22. El Raheem, ZF, Salman, SM: On a discretization process of fractional-order logistic dierential equation. J. Egypt. Math. Soc. 22, 407-412 (2014)

23. El-Sayed, AMA, El-Raheem, ZF, Salman, SM: Discretization of forced Dung system with fractional-order damping. Adv. Dier. Equ. 2014, 66 (2014)

24. El-Sayed, AMA, Salman, SM: On a discretization process of fractional-order Riccati dierential equation. J. Fract. Calc. Appl. 4(2), 251-259 (2013)

25. Elaidy, SN: An Introduction to Dierence Equations, 3rd edn. Undergraduate Texts in Mathematics. Springer, New York (2005)

26. Palis, J, Takens, F: Hyperbolicity and Sensitive Chaotic Dynamics and Homoclinic Bifurcation. Cambridge University Press, Cambridge (1993)

27. Puu, T: Attractors, Bifurcation and Chaos: Nonlinear Phenomena in Economics. Springer, Berlin (2000)28. Udwadia, FE, von Bremen, H: A note on the computation of the largest p-Lyapunov characteristic exponents for nonlinear dynamical systems. J. Appl. Math. Comput. 114, 205-214 (2000)