domain: math nlin

text:
Chunxiao Yang
email: chunxiao.yang@ls2n.fr1

Ina Taralova
email: ina.taralova@ls2n.fr

Jean Jacques Loiseau
email: jean-jacques.loiseau@ls2n.fr

Fractional Chaotic system solutions and their impact on chaotic behaviour

Keywords: Fractional calculus, Numerical solution, Fractional Chaotic system, Non-linear dynamics

This paper is devoted to the analysis of calculation methods for solving fractional chaotic systems and the impact of these different approaches on the behavior of the fractional chaotic system. Two widely used time domain fractional differential equations solving approaches are discussed, the fractional ABM correctorpredictor method based on Caputo fractional derivative definition, and the long memory calculation approach based on Grunwald fractional derivative. These numerical solutions calculation methods are employed to depict the phase portrait of a class of commensurate fractional chaotic systems. The Lyapunov exponent and bifurcation diagrams of the systems over various fractional orders and parameters are illustrated to detect the impact on the dynamics of the chaotic system applying different calculation approaches.

Introduction

Chaos is a random-like behavior exhibited by many nonlinear dynamic systems. The very first proponent on this topic can be dated back to 1880 while the three body problem was studied [START_REF] Poincaré | Sur le problème des trois corps et les équations de la dynamique[END_REF]. Eighty years latter, when Edward Lorenz worked on weather prediction, the so-called 'Lorenz attractor' was found [START_REF] Lorenz | The Predictability of Hydrodynamic Flow[END_REF]. By giving it a description and a poetic name of 'butterfly effect', the gate of mathematical and scientific world of Chaos was opened. Since then, many researchers have tried to uncover the deterministic laws behind the apparently random states of disorder of different chaotic systems. One of the characteristics of chaotic system is that it is very sensitive to the initial conditions as described by butterfly effect. This sensitivity can be measured by Lyapunov Exponent(LE) which calculates the rate of exponential divergence of trajectories starting from two close initial conditions. This characteristics also contributes to the application of chaotic systems in many domains of science and engineering, such as biology [START_REF] Liz | Chaos in discrete structured population models[END_REF], economics [START_REF] Kyrtsou | Evidence for chaotic dependence between US inflation and commodity prices[END_REF], finance [START_REF] Fernando | Applying the theory of chaos and a complex model of health to establish relations among financial indicators[END_REF], cryptography [START_REF] Qiao | Efficient Pseudo-chaotic Number Generator for Cryptographic Applications[END_REF] [START_REF] Babaei | A novel text and image encryption method based on chaos theory and DNA computing[END_REF] and etc. In the meantime, fractional calculus is considered as the generalization of classical integer-order integration and differentiation operators to real, or complex orders [START_REF] Petráš | Fractional-order nonlinear systems[END_REF]. Many mathematicians have discussed the fractional calculus since 1695 by introducing different mathematical characterisations (definitions) for fractional derivative and integration. In many cases, these characterisations are equivalent if the initial conditions are ad hoc [START_REF] Diethelm | The Analysis of Fractional Differential EquationsThe Analysis of Fractional Differential Equations[END_REF], and the most well known three are Rienmann-Liouville(RL), Grunwald-Letnikov(GL) and Caputo characterisations. The analysis and discussion of fractional calculus remained purely in the domain of mathematics for centuries. It was not until 1980s that the application of fractional calculus in the domain of science and engineering has started to be studied and explored. Due to the memory effect possessed by fractional calculus, it is considered to be suitable to model many real-life systems. After years of research, the fractional differential equations have now been used in diverse disciplines like physics, biology, and economics, etc. [START_REF] Mainardi | Fractional Calculus and Waves Linear Viscoelasticity: An Introduction to Mathematical Models[END_REF] [START_REF] Tarasov | Macroeconomic models with long dynamic memory: Fractional calculus approach[END_REF]. The fractional chaotic system also attracts a lot of attention. The difficulties for this research owes to the intricate geometric interpretation of fractional derivatives [START_REF] Li | A Novel Image Encryption Algorithm Based on a Fractional-Order Hyperchaotic System and DNA Computing[END_REF] and the fact that there exist, as mentioned above, different definitions for fractional derivatives. One basically considers continuous systems, and uses numerical methods to approximate the solution. In the case of a fractional system, the discrete approximating system may inherit the chaotic behaviour of the initial continuous system, but this relationship is somehow complex. What adds to the intricacy is that the chaotic behavior of the approximating can be different for different numerical methods employed to solve the fractional differential equations [START_REF] Yang | Design of a Fractional Pseudo-Chaotic Random Number Generator[END_REF]. Therefore, the understanding of the impact on the chaoticity of the system applying one or another numerical calculation approaches is of great importance, in order to choose the most appropriate one for a given application. In the following, two numerical calculation methods under GL and Caputo characeterizations for fractional differential equations are recalled. Then, we employ both methods to obtain the states of two fractional chaotic systems extended from classical integer order chaotic system. The impact on the chaoticity of the systems applying the two approaches has been analyzed in terms of LE and from the aspect of bifurcation diagram and time responses.

Preliminaries on fractional calculus and fractional systems

In this section, some preliminaries on fractional calculus and fractional systems are introduced to give a rough idea on the topic. The widely-accepted stability criteria for a commensurate fractional system is also illustrated.

Fractional calculus

As mentioned before, the fractional calculus studies the fractional derivative and integral which can be considered as the extension of classical integer order differentiation and integration to real or complex orders. In the long history of the study of fractional calculus, many mathematicians have contributed and introduced different characterisations(referred as 'definitions' in many paper) towards the topic. Here after, we give two well-known definitions Grünwald-Letnikov (GL) and Caputo definitions [14][15].

The fractional derivatives under GL characterisation can be writen as

a D α t f (t) = lim h→0 1 h α [ t-a h ] j=0 (-1) j α j f (t -jh) (1)
The term [ t-a h ] in equation ( 1) stands for the integer part of t-a h ; a and t are the bounds of the derivative operation a D α t for f (t); α represents the fractional derivative order. The term α j in (1) is defined in [START_REF] Lorenz | The Predictability of Hydrodynamic Flow[END_REF], where Γ (.) is the Eular Gamma function in the form of [START_REF] Liz | Chaos in discrete structured population models[END_REF].

α j = Γ (α + 1) Γ (j + 1)Γ (α -j + 1) (2)
Γ (α) = ∞ 0 t α-1 e t dt (3)
The caputo type fractional derivative holds the form as following,

c a D α t f (t) = 1 Γ (n -α) t a f (n) (τ ) (t -τ ) α-n+1 dτ, for n -1 < α < n (4)
where α denotes the fractional derivative order; a and t are the bounds for the operation; n is the smallest integer greater than α; Γ (.) is the Eular Gamma function in equation (3); and f (n) (t) is the n-th derivative of f (t). The Caputo type fractional derivative is often used for engineering application since the fractional differential equations with this type of derivative can provide the applied problem with an interpretive initial condition.

Fractional system

A fractional system is a dynamic system which can be modeled by fractional differential equations [START_REF] West | Physics of fractal operators[END_REF]. A general form of fractional system is as follows,

0 D αi t x i (t) = f i (x 1 (t) , x 2 (t) , ..., x n (t) , t) x i (0) = c i , i = 1, 2, ..., n. (5)
In ( 5), x i (0)(i = 1, 2, ...n) denotes the initial conditions for each component constituting the state vectors; α i (i = 1, 2, ...n) is the fractional derivative order for i-th differential equations consisting the system, and f i is a linear or nonlinear function.

The equilibrium points of system ( 5) can be obtained by solving equation f i (x) = 0(i = 1, 2, ..., n). If a commensurate system with α i = α, i = 1, 2, ..., n is considered, then, according to the stability theorem defined in [START_REF] Tavazoei | A necessary condition for double scroll attractor existence in fractional-order systems[END_REF], the equilibrium points are locally asymptotically stable if the eigenvalue of the Jacobian matrix of system (5) satisfies the following equation evaluated at equilibria.

|arg (eig (J))| = |arg (λ i )| > α π 2 , i = 1, 2, ..., n (6)
where J denotes the Jacobian matrix of (5), λ i (i = 1, 2, ...n) are its eigenvalues.

Numerical calculation methods for fractional differential equations

In this section, two numerical solutions calculation methods for fractional differential equations are introduced. The methods are based on Grünwald-Letnikov and Caputo fractional derivative characterisations.

Grünwald-Letnikov calculation method

The explicit numerical approximation of q-th derivative under GL characterisation at the points kh, (h = 1, 2, ...) is expressed as follows [START_REF] Podlubny | Fractional Differential Equations[END_REF]

(k-Lm)/h D α t k f (t) ≈ h -α k j=0 (-1) j α j f (t k-j ). (7)
In expression [START_REF] Babaei | A novel text and image encryption method based on chaos theory and DNA computing[END_REF], L m is the memory length; t k = kh, where h is the calculation time step; the binomial coefficient (-1) j α j can be denoted as c (α) j (j = 0, 1, ...) which is expressed use the following expression [START_REF] Dorcak | Numerical Models for the Simulation of the Fractional-Order Control Systems[END_REF],

c (α) 0 = 1, c (α) j = 1 - 1 + α j c (α) j-1 . (8)
Thus, the general numerical solution of fractional differential equation described by equation( 9) can be expressed as given in [START_REF] Mainardi | Fractional Calculus and Waves Linear Viscoelasticity: An Introduction to Mathematical Models[END_REF].

a D α t y(t) = f (y(t), t) (9)
y(t k ) = f (y(t k ), t k )h α - k j=ν c (α) j y(t k-j )) (10)
The sum in [START_REF] Mainardi | Fractional Calculus and Waves Linear Viscoelasticity: An Introduction to Mathematical Models[END_REF] stands for the memory term. If a 'long memory effect' is considered, then the lower index ν = 1 for all k, otherwise ν = 1 for k < (L m /h) and ν = k -L m for k > (L m /h).

Fractional ABM corrector-predictor method

The fractional ABM corrector-predictor method is another widely used time domain numerical calculation method in the domain of engineering. It is a generalization of the classical Adams-Bashforth-Moulton integrator which is used for the numerical calculation of classical first order problem.

From the analytical point of view, the fractional differential equations under Caputo characterization with initial conditions y k (0) = y k 0 , k = 0, 1, 2...m -1 where m := α , is equivalent to Volterra integral equation expressed as follows,

y(x) = α -1 k=0 y (k) 0 x k k! + 1 Γ (α) x 0 (x -t) α-1 f (t, y(t))dt (11)
The algorithm is developed on a uniform grid {t n = nh : n = 0, 1, ...N }. The basic idea of the algorithm is to obtain the approximation of the latter point on the grid from the former point. Detailed formula derivation for the algorithm can be found in [START_REF] Diethelm | A predictor-corrector approach for the numerical solution of fractional differential equations[END_REF]. Here, we only give out the derived equations for the next states values in equation ( 12)- [START_REF] Caputo | Linear Models of Dissipation whose Q is almost Frequency Independent-II[END_REF].

y h (t n+1 ) = α -1 k=0 t k n+1 k! y (k) 0 + h α Γ (α + 2) f (t n+1 , y P h (t n+1 )) + h α Γ (α + 2) n j=0 a j,n+1 f (t j , y h (t j )), (12)
y P h (t n+1 ) = α -1 k=0 t k n+1 k! y (k) 0 + 1 Γ (α) n j=0 b j,n+1 f (t j , y h (t j )). ( 13
)
a j,n+1 =      n α+1 -(n -α)(n + 1) α , if j = 0, (n -j + 2) α+1 + (n -j) α+1 -2(n -j + 1) α+1 , if 1 ≤ j ≤ n, 1, if j = n + 1. ( 14
) b j,n+1 = h α α ((n + 1 -j) α -(n -j) α ) (15)
In the above equations, y h (t n+1 ) stands for the next state, y P h (t n+1 ) denotes the predictor value for the next state, a and b are coefficients.

Fractional chaotic Chen and Lu systems

Fractional chaotic Chen systems

The system equation for fractional Chen system can be expressed as following [START_REF] Lu | A note on the fractional-order Chen system[END_REF],

f c (x) =    D αc x 1 (t) = a c (x 2 (t) -x 1 (t)) D αc x 2 (t) = (c c -a c ) x 1 (t) -x 1 (t) x 3 (t) + c c x 2 (t) D αc x 3 (t) = x 1 (t) x 2 (t) -b c x 3 (t) (16)
In the equation, D αc denotes the fractional derivative with order α c , (a c , b c , c c ) are the parameters of the system. The system is an extension from integer order chaotic Chen system studied in [START_REF] Lu | A new chaotic attractor coined[END_REF]. The equilibria of the system can be obtained through the same way as its original integer order system, by setting the right hand sand system equation equal to zero f c (x * ) = 0 as given below,

   a c (x 2 (t) -x 1 (t)) = 0 (c c -a c ) x 1 (t) -x 1 (t) x 3 (t) + c c x 2 (t) = 0 x 1 (t) x 2 (t) -b c x 3 (t) = 0 (17)
The singularity of the equilibrium points can also be acquired through the classical method as given below, by evaluating the eigenvalue of the jacobian matrix of the system at equilibrium points.

det (λ c I-J c ) =   λ c + a c -a c 0 -c l + a l + x * 3 λ c -c c x * 1 -x * 2 -x * 1 λ c + b c   = 0 (18)
J c in equation ( 18) represents the Jacobian matrix of the system equation, I is the identity matrix, λ c denotes the eigenvalue, and (x * 1 , x * 2 , x * 3 ) stands for the equilibrium point. The singularity of the three equilibrium points of fractional Chen system for system parameters (a c , b c , c c ) = (35, 3, 28) can be obtained through above analytical expressions and are given in Table 1.

Fractional chaotic Lu system

The system equation for fractional chaotic Lu system extended from integer order Lu system can be described as follows [START_REF] Deng | Chaos synchronization of the fractional Lü system[END_REF],

f l (x) =    D α l x 1 (t) = a l (x 2 (t) -x 1 (t)) D α l x 2 (t) = -x 1 (t) x 3 (t) + c l x 2 (t) D α l x 3 (t) = x 1 (t) x 2 (t) -b l x 3 (t) (19)
where D α l denotes the fractional derivative with order α l , a l , b l , and c l are the parameters of the system. The equilibrium points of the system can be acquired calculating the solutions of the following system of equations,

   a l (x 2 (t) -x 1 (t)) = 0 -x 1 (t) x 3 (t) x 1 (t) + c l x 2 (t) = 0 x 1 (t) x 2 (t) -b l x 3 (t) = 0 (20)
The singularity of the equilibria can be obtained the same way as discussed previously for the fractional Chen system through the following identities,

det (λ l I-J l ) =   λ l + a l -a l 0 x * 3 λ l -c c x * 1 -x * 2 -x * 1 λ l + b l   = 0 (21)
where J l in equation ( 21) represents the Jacobian matrix of the fractional Lu system, λ l denotes the eigenvalue, and (x * 1 , x * 2 , x * 3 ) stands for the equilibrium point. When the parameters of the system is set to (a c , b c , c c ) = (36, 3, 20), three equilibirum points E * 1 = (0, 0, 0), E * 2 = (7.460, 7.460, 20) and E * 3 = (-7.460, -7.460, 20) can be obtained applying equation [START_REF] Lu | A note on the fractional-order Chen system[END_REF]. The singularity of the equilibria is also given in Table 1.

Solutions for the chaotic systems applying different approaches

In this section, the solutions for fractional Chen and Lu solution are obtained applying both GL method and fractional ABM corrector-predictor method dis-cussed in section 3. The impact of the two approaches on the chaotic behavior of the systems are also discussed.

Chaotic system applying GL method

With the numerical solution of fractional differential equation calculated under GL method derived as in [START_REF] Mainardi | Fractional Calculus and Waves Linear Viscoelasticity: An Introduction to Mathematical Models[END_REF], the calculation for the states of fractional Chen system and fractional Lu system (expression ( 16) and ( 19)) can be expressed by the following identities ( 22) and ( 23), respectively. To be mentioned is that in our work, the 'long memory effect' is adopted applying GL method which means that the number ν in equaitons ( 22) and ( 23) is equal to 1. The time step h in above equations is set to a fixed value 0.001. We plotted the phase portraits of the two systems with fractional orders α c = 0.9 and α l = 0.95 in Fig. 1a and1b, respectively. The parameters and initial conditions for Chen system are (35, 3, 28) and (-9, -5, 14). Those of Lu system are chosen to be (36, 3, 20) and (0.2, 0.5, 0.3).

               x 1 (n) = (a c (x 2 (n) -x 1 (n -1)))h αc - n j=ν c (αc) j x 1 (n -j) x 2 (n) = ((c c -a c )x 1 (n) -x 1 (n)x 3 (n -1) + c c x 3 (n -1))h αc - n j=ν c (αc) j x 2 (n -j) x 3 (n) = (x 1 (n)x 2 (n) -b c x 3 (n -1))h αc - n j=ν c (αc) j x 3 (n -j) (22)                x 1 (n) = (a l (x 2 (n -1) -x 1 (n -1)))h α l - n j=ν c (α l ) j x 1 (n -j) x 2 (n) = (-x 1 (n)x 3 (n -1) + c l x 2 (n -1))h α l - n j=ν c (α l ) j x 2 (n -j) x 3 (n) = (x 1 (n)x 2 (n) -b l x 3 (n -1))h α l - n j=ν c (α l ) j x 3 (n -j) (23

Chaotic systems applying ABM corrector-predictor approach

Based on the fractional ABM corrector-predictor numerical calculation approach for the solution of fractional differential equations given in equations ( 12)-( 15), the states of fractional Chen system applying ABM predictor corrector approach can be expressed as follows,

X c (n + 1) = X c (0) + h αc Γ (α c + 2) f c (X P c (n + 1)) + h αc Γ (α c + 2) n j=0 a j,n+1 f c (X c (j)) (24)
X P c (n + 1) = X c (0) + 1 Γ (α c ) n j=0 b 1 j,n+1 f c (X c (j)) a j,n+1 =    n αc+1 -(n -α c )(n + 1) αc , ifj = 0, (n -j + 2) αc+1 + (n -j) αc+1 -2(n -j + 1) αc+1 , if1 ≤ j ≤ n, 1, ifj = n + 1. b j,l+1 = h αc α c ((n + 1 -j) αc -(n -j) αc ) (25)
In the above expressions, X c (n + 1), X c (n) and X P c (n + 1) are state vectors composed of all the state components x 1 , x 2 , and x 3 ; α c is the fractional order between (0, 1); f c stands for the Chen system equations. The formula for the calculation of the states of fractional Lu system can be obtained by substituting the state vectors, fractional order and system equations in equations ( 24)-(25) with X l , α l and f l where 0 < α l < 1. The phase portraits of the two systems acquired employing corrector-predictor approach are given in Fig. 1c and1d, respectively. The fractional orders, parameters and initial conditions are the same as those for the GL method.

Impact on system chaoticity with chosen methods

For the work in this section, we used the same parameters and initial conditions for the two systems as adopted in the previous section, which are (a c , b c , c c ) = (35, 3, 28), X c (0) = (-9, -5, 14); (a l , b l , c l ) = (36, 3, 20), X l (0) = (0.2, 0.5, 0.3), respectively. The time step h is set to 0.005. The MATLAB code [START_REF] Garrappa | Predictor-corrector PECE method for fractional differential equations[END_REF] for ABM corrector-predictor method and [START_REF] Danca | Matlab Code for Lyapunov Exponents of Fractional-Order Systems[END_REF] is employed for the following simulation and the calculation of LE. According to the stability criteria introduced by the inequality (6) , the reference [START_REF] Tavazoei | A necessary condition for double scroll attractor existence in fractional-order systems[END_REF] states that for a fractional system D α x = f (x) to remain chaotic, a necessary condition is keeping the eigenvalues λ in the unstable region, which gives the following equation for the fractional derivative order α.

α > 2 π tan -1 |Im(λ)| Re(λ) (26)
(a) GL method Chen system (αc = 0.9) (b) GL method Lu system (α l = 0.95) (c) GL method Chen system (αc = 0.9 (d) ABM method Lu system (α l = 0.95 where λ denotes the eigenvalues of the Jacobian matrix of the system, α is the commensurate fractional order. Therefore, for the given parameter values, the fractional chaotic Chen system should have a fractional order α c greater than or equal to 0.8244. In Fig. 2, we plot the phase portrait of fractional Chen system at boundary fractional values 0.82 and 0.83 applying both GL and ABM corrector-predictor methods. The time response of the last 2000 states obtained through both methods are also given. The states calculated by GL method is in red and ABM corrector-predictor in blue. It is not difficult to observe from Fig. 2a and 2c that with order 0.82 there are only one red point in the figure, which indicates that the states stays at the same fixed point applying GL method. Whereas for the applied ABM method(blue dots), they appear to have a shape of the attractors. When the system order is equal to 0.83, both methods display the shapes with attractors. This indicates that when applying GL calculation method with long memory effect, the system's dynamic behavior is in accordance with the stability criteria given by equation (26). While the ABM calculation method applied in this paper provides the system with a smaller derivative order for the system to be chaotic. The time response figures given by Fig. 2b and 2d confirms the founding. The It is clear that for derivative order 0.82, the red attractors stays at the same value for the three state vector components x 1 , x 2 and x 3 , while the blue curves appear to be oscillating. We also give the Lyapunov exponent and bifurcation diagrams over different fractional orders of fractional Chen and Lu systems in Fig. 3. For each fractional derivative orders, 10 4 states were generated and the LEs were calculated throughout the iterations. The LE spectrum curves in 3a and 3b are obtained by combining LE values of the last iteration for every evaluated orders. The plots show that only x 1 component possesses LE value greater than 0 applying both methods. It can be observed that applying ABM corrector-predictor approach, for the fractional Chen system, the LE of x 1 greater than 0 appears before order 0.53, whereas for GL method, the LE exceeds 0 after fractional order of 0.8. The LEs for fractional Lu system calculated using both methods show the similar results, with ABM method having a smaller chaotic fractional derivative value. This is in accordance with our previous findings concerning the phase portrait and time response which draws to the conclusion that GL method give a more accurate approximation of original fractional system.

Apart from this, from the y-coordinates of the bifurcation diagram where the system is non-chaotic, it can be observed that the solution obtained using ABM method stays at the equilibrium point as obtained through analytical analysis.

(a) Chen system GL method (b) Chen system ABM method (c) Lu system GL method (d) Lu system ABM method The LEs results and bifurcation diagram over different parameters of fractional Lu system are also given in Fig. 4 to illustrate the dynamics possesses by the system. We set the system fractional order fixed to 0.9. It can be observed that applying different numerical calculation methods, the system dynamic is quite different. It is worth mentioning that the results for different parameters are conducted by changing one parameter at a time and fixing the other two unchanged.

Conclusion

In this paper, we recalled two numerical solutions calculation methods for fractional differential equations adopting Grünward-Leinikov and Caputo charac- terization of fractional derivative, respectively. Two fractional chaotic systems, fractional Chen system and fractional Lu system are discussed and their discretized states were calculated employing both methods. The results show that compared to the adopted ABM corrector-predictor method, the GL approach with long memory effect provide the original fractional system with a better approximation in coherence with the analytical studies. At the contrary, employing ABM method, the approximation accuracy appears to be deteriorate. However, in terms of chaoticity, it has a greater chaotic range for fractional derivatives.

Fig. 1 :

 1 Fig. 1: Phase Portrait of fractional Chen and Lu systems characterized by GL and ABM method

( a )Fig. 2 :

 a2 Fig. 2: Phase Portrait and time response of Chen system at boundary fractional order values

Fig. 3 :

 3 Fig. 3: LE and bifurcation results for Chen and Lu systems over different fractional derivatives employing different methods (a c , b c , c c ) = (35, 3, 28), (a l , b l , c l ) = (36, 3, 20)

( a )

 a Lu system GL method a l LE results (b) Lu system ABM method a l LE results (c) Lu system GL method b l LE results (d) Lu system ABM method b l LE results (e) Lu system GL method c l LE results (f) Lu system ABM method C l LE results

Fig. 4 :

 4 Fig. 4: LE and bifurcation results for Chen and Lu systems over different fractional parameters employing different methods

Table 1 :

 1 ) Fractional Chen and Lu systems' equilibria and their singularity

	System	Equilibrium	λ1	Eigenvalue λ2, λ3	singularity
	Fractional	(0, 0, 0)	-30.8359	23.8359,-3	Saddle
	Chen	(-7.9379,-7.9379,21) -18.4280 4.2140±14.8846i Saddle Focus
	system	(7.9379,7.9379,21) -18.4280 4.2140±14.8846i Saddle Focus
	Fractional	(0, 0, 0)	-36		20, -3	Saddle
	Lü	(-7.460, -7.460, 20) -22.6516 1.8258 ± 13.6887i Saddle Focus
	system	(7.460, 7.460, 20) -18.4280 1.8258 ± 13.6887i Saddle Focus