Zaid M Odibat
email: z.odibat@gmail.com

Prince Abdullah Bin Ghazi

Adaptive feedback control and synchronization of non-identical chaotic fractional order systems

Keywords: Chaos synchronization, Fractional order system, Caputo fractional derivative, Stability, Feedback control

Adaptive feedback control and synchronization of non-identical chaotic fractional order systems

Zaid M.

Introduction

As a kind of characteristic of nonlinear systems, chaos is a bounded unstable dynamic behavior that exhibits sensitive dependence on initial conditions and includes infinite unstable periodic motions. Chaos has been investigated and studied in mathematical and physical communities in the last few decades because of its grate applications in many fields such as secure communication, data encryption, flow dynamics and biomedical engineering [START_REF] Chen | Chaos Control: Theory and Applications[END_REF]. The research efforts have devoted to the chaos control and chaos synchronization problems in many dynamical systems [START_REF] Yamada | Stability theory of synchroized motion in coupled-oscillator systems[END_REF][START_REF] Pecora | Synchronization in chaotic systems[END_REF][START_REF] Ott | Controling chaos[END_REF][START_REF] Aziz-Alaoui | Synchronization of chaos[END_REF].

On the other hand, the development of models based on fractional order differential systems has recently gained popularity in the investigation of dynamical systems. Fractional derivatives provide an excellent instrument for the description of memory and hereditary properties of various materials and processes. The main reason for using the integer-order models was the absence of solution methods for fractional differential equations. The advantages or the real objects of the fractional order systems [START_REF] Petráš | A note on the fractional-order Chua's system[END_REF] are that we have more degrees of freedom in the model and that a "memory" is included in the model.

Recently, studying fractional order systems has become an active research area. The chaotic dynamics of fractional order systems began to attract much attention in recent years. It has been shown that fractional order systems, as generalizations of many well-known systems, can also behave chaotically, such as the fractional Duffing system [START_REF] Ge | Chaos in a fractional order modified Duffing system[END_REF], the fractional Chua system [START_REF] Petráš | A note on the fractional-order Chua's system[END_REF][START_REF] Hartley | Chaos in a fractional order Chua's system[END_REF], the fractional Rössler system [START_REF] Li | Chaos and hyperchaos in fractional order Rössler equations[END_REF], the fractional Chen system [START_REF] Li | Chaos in Chen's system with a fractional order[END_REF][START_REF] Li | Chaos in the fractional order Chen system and its control[END_REF][START_REF] Lu | A note on the fractional-order Chen system[END_REF], the fractional Lorenz system [START_REF] Grigorenko | Chaotic dynamics of the fractional Lorenz system[END_REF], the fractional Arneodo's system [START_REF] Arneodo | Asymptotic chaos[END_REF][START_REF] Lu | Chaotic dynamics and synchronization of fractional-order Arneodos systems[END_REF], the fractional Lü system [START_REF] Deng | Chaos synchronization of the fractional Lü system[END_REF], the fractional Newton-Leipnik [START_REF] Sheu | Chaos in the Newton-Leipnik system with fractional order[END_REF] and the fractional Chen-Lee system [START_REF] Tam | Parametric study of the fractional-order Chen-Lee system[END_REF]. In [START_REF] Hartley | Chaos in a fractional order Chua's system[END_REF][START_REF] Li | Chaos and hyperchaos in fractional order Rössler equations[END_REF][START_REF] Li | Chaos in Chen's system with a fractional order[END_REF][START_REF] Li | Chaos in the fractional order Chen system and its control[END_REF] it has been shown that some fractional order systems can produce chaotic attractors with order less than 3.

Meanwhile, chaotic dynamics of fractional order systems began to attract much attention in recent years.

A challenging problem is the control and synchronization of chaotic systems. Recent studies show that chaotic fractional order systems can also be synchronized [START_REF] Li | Synchronization in fractional-order differential systems[END_REF]- [START_REF] Zhu | Chaos and synchronization of the fractional-order Chua's system[END_REF]. In many literatures, synchronization among fractional order systems is only investigated through numerical simulations that are based on the stability criteria of linear fractional order systems, such as the work presented in [START_REF] Sheu | Chaos in a new system with fractional order[END_REF]- [START_REF] Zhu | Chaos synchronization of the fractional-order Chen's system[END_REF], or based on Laplace transform theory, such as the work presented in [START_REF] Wang | Extending synchronization scheme to chaotic fractional-order Chen systems[END_REF]- [START_REF] Zhu | Chaos and synchronization of the fractional-order Chua's system[END_REF].

In the present paper, we study the synchronization of fractional order chaotic systems with different fractional orders via nonlinear control. By taking the fractional version of Chen system (which belongs to the double scroll attractor family) and Rössler system (which belongs to the one scroll attractor family) as examples, we show that fractional order systems with different fractional orders can be synchronized. Based on stability results of fractional order systems, using the drive-response concept, an adaptive feedback control is constructed to achieve synchronization between two 3D fractional order systems that have nonidentical orders. The synchronization controllers are investigated theoretically and then, numerical simulations are presented to verify the theoretical analysis.

Basic concepts

There are several definitions of a fractional derivative of order α > 0 [START_REF] Oldham | The Fractional Calculus[END_REF]- [START_REF] Caputo | Linear models of dissipation whose Q is almost frequency independent, Part II[END_REF]. The two most commonly used are Riemann-Liouville and Caputo definitions. Each definition uses Riemann-Liouville fractional integration and derivatives of whole order. The difference between the two definitions is in the order of evaluation. Riemann-Liouville fractional integral operator of order α ≥ 0 of the function f (t) is defined as,

J α f (t) = 1 Γ(α) t 0 (t -τ ) α-1 f (τ )dτ, t > 0. (1)
Some properties of the operator J α can be found, for example, in [START_REF] Podlubny | Fractional Differential Equations[END_REF][START_REF] Gorenflo | Fractional calculus: Integral and differential equations of fractional order[END_REF]. We recall only the following, for μ ≥ -1, α, β ≥ 0 and γ > -1, we have,

J α J β f (t) = J α+β f (t), J α t γ = Γ(γ + 1) Γ(α + γ + 1) t α+γ .
The Laplace transform of Riemann-Liouville fractional integral satisfies,

L{J α f (t)} = s -α L{f (t)}. (2)
In this study, Caputo definition is used and the fractional derivative of f (t) is defined as,

D α f (t) = J m-α D m f (t), D α f (t) = 1 Γ(m -α) t 0 f (m) (τ ) (t -τ ) α-m+1 dτ, (3) for m -1 < α ≤ m, m ∈ IN, t > 0.
Caputo's definition has the advantage of dealing properly with initial value problems in which the initial conditions are given in terms of the field variables with their integer order which is the case in most physical processes. Fortunately, the Laplace transform of Caputo fractional derivative satisfies,

L{D α f (t)} = s α L{f (t)} - m-1 k=0 f (k) (0 + )s α-1-k , (4)
where m -1 < α ≤ m. The Laplace transform of Caputo fractional derivative requires the knowledge of the (bounded) initial values of the function and of its integer derivatives of order k = 1, 2, • • • , m -1.

Stability analysis of fractional systems

Stability analysis of fractional order systems, which is of main interest in control theory, has been thoroughly investigated where necessary and sufficient conditions have been derived [START_REF] Matignon | Stability results of fractional differential equations with applications to control processing[END_REF]- [START_REF] Tavazoei | A note on the stability of fractional order systems[END_REF] (see also references therein). In this section, we recall the main stability results. For this object, we consider the following n dimensional fractional order system,

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ d α 1 x1 dt α 1 = f1(x1, x2, • • • , xn), d α 2 x2 dt α 2 = f2(x1, x2, • • • , xn), . . . d αn xn dt αn = fn(x1, x2, • • • , xn), ( 5
)
where α i is a rational number between 0 and 1 and

d α i dt α i is the Caputo fractional derivative of order α i , for i = 1, 2, • • • , n. Assume that α i = k i /m i , (k i , m i ) = 1, k i , m i ∈ IN , for i = 1, 2, • • • , n.
Let m be the least common multiple of the denominators m i 's of α i 's.

First, if the system (5) is a linear system, that is

[f 1 (x), f 2 (x), • • • , f n (x)] T = [a ij ] n
i,j=1 x = Ax, where x ∈ IR n , then we have the following results:

• If α = α 1 = α 2 = • • • = α n ,
then the fractional order system (5) is asymptotically stable iff |arg(spec(A))| > απ/2. In this case the components of the state decay towards 0 like t -α [START_REF] Matignon | Stability results of fractional differential equations with applications to control processing[END_REF].

• If α i 's are rational numbers between 0 and 1, then the system ( 5) is asymptotically stable if all roots λ of the equation

det diag(λ mα1 , λ mα2 , • • • , λ mαn ) -A = 0 satisfy |arg(λ)| > γπ/2, where γ = 1/m [39].
Second, if function f i has second continuous partial derivatives in a ball centered at an equilibrium point

x * = (x 1 , x 2 , • • • , x n ), that is f i (x 1 , x 2 , • • • , x n ) = 0, for i = 1, 2, • • • , n
, then we have the following results:

• If α = α 1 = α 2 = • • • = α n , then the equilib- rium point x * of system (5) is asymptotically stable if |arg(spec(J| x * ))| > απ/2,
where the matrix J is the Jacobian matrix of the system (5) that is defined as J = [ ∂fi ∂xj ] n i,j=1 [START_REF] El-Sayed | Equilibrium points, stability and numerical solutions of fractionalorder predator-prey and rabies models[END_REF].

• If α i 's are rational numbers between 0 and 1, then the equilibrium point x * of system ( 5) is asymptotically stable if all roots λ of the equa-

tion det diag(λ mα1 , λ mα2 , • • • , λ mαn ) -J| x * = 0 satisfy |arg(λ)| > γπ/2, where γ = 1/m [41].
Fractional order systems are, at least, as stable as their integer order counterpart, because systems with memory are typically more stable than their memoryless counterpart [START_REF] El-Sayed | Equilibrium points, stability and numerical solutions of fractionalorder predator-prey and rabies models[END_REF]. The previous stability results play an important role in studying the existence of chaotic attractors and the synchronization of fractional order systems.

fractional Chen system

Chen and Ueta [START_REF] Chen | Yet another chaotic attractor[END_REF] introduced, in 1999, the Chen system which is similar but not topologically equivalent to Lorenz system. It is a chaotic system with a double scroll attractor. The fractional version of Chen system reads as,

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ d α 1 x dt α 1 = a(y -x), d α 2 y dt α 2 = (c -a)x -xz + cy, d α 3 z dt α 3 = xy -bz, ( 6
)
where 0 < α 1 , α 2 , α 3 ≤ 1. Integer order Chen system displays chaotic attractors, for example, when (a, b, c) = [START_REF] Hilfer | Applications of Fractional Calculus in Physics[END_REF][START_REF] Pecora | Synchronization in chaotic systems[END_REF][START_REF] Wang | Extending synchronization scheme to chaotic fractional-order Chen systems[END_REF]. Simulations are performed to obtain chaotic behavior of fractional order Chen system for different fractional orders α, when (a, b, c) = [START_REF] Hilfer | Applications of Fractional Calculus in Physics[END_REF][START_REF] Pecora | Synchronization in chaotic systems[END_REF][START_REF] Wang | Extending synchronization scheme to chaotic fractional-order Chen systems[END_REF]. For example, chaotic attractors are found in [START_REF] Yan | On chaos synchronization of fractional differential equations[END_REF] when α = (0.95, 0.95, 0.95). In [START_REF] Li | Chaos in the fractional order Chen system and its control[END_REF] chaotic behaviors are found when α = (0.9, 0.9, 0.9). Moreover, in [START_REF] Wang | Extending synchronization scheme to chaotic fractional-order Chen systems[END_REF] and [START_REF] Tavazoei | Chaotic attractors in incommensurate fractional order systems[END_REF], it is found that for the parameters α = (0.985, 0.99, 0.98) and α = (0.8, 1, 0.9), respectively, fractional order Chen system can display chaotic attractors.

The Jacobian matrix of system ( 6), evaluated at the equilibrium point x * = (x * , y * , z * ), is

J C = ⎛ ⎜ ⎝ -35 35 0 -7 -z * 28 -x * y * x * -3 ⎞ ⎟ ⎠ (8)
Now, according to stability results and the results presented in [START_REF] Tavazoei | Chaotic attractors in incommensurate fractional order systems[END_REF], a necessary condition for fractional Chen system to exhibit chaotic attractors similar to its integer order counterpart is:

min i {|arg(λ i )|} ≤ γπ/2, ( 9
)
where γ = 1/m and λ i 's are the roots of the equation det diag(λ mα1 , λ mα2 , λ mα3 ) -J C | x * = 0, for every equilibrium point x * . Otherwise, one of these equilibrium points becomes asymptotically stable and attracts the nearby trajectories. (

) 10
According to [START_REF] Li | Chaos and hyperchaos in fractional order Rössler equations[END_REF], the equilibrium points P 2 and P 3 are asymptotically stable if,

α < 2 π tan -1 14.88464 4.21398 = 0.82436. ( 11
)
Therefore, system (6) has the necessary condition α > 0.82436 for exhibiting double scroll chaotic attractor.

For α = (1, 0.9, 0.9) and according to the last two equilibrium points P 2 and P 3 , the equation [START_REF] Lu | A note on the fractional-order Chen system[END_REF] and so we get,

det diag(λ mα1 , λ mα2 , λ mα3 ) -J C | x * = 0 becomes, λ 28 -25λ 19 + 35λ 18 -21λ 10 + 105λ 9 + 4410 = 0,
min i {|arg(λ i )|} = 0.12458 < π/20. ( 13
)
Hence, system (6) satisfies the necessary condition for exhibiting a double scroll attractor when α = (1, 0.9, 0.9). Using similar analysis we can confirm that fractional Chen system (6) satisfies the necessary condition (9) for exhibiting a double scroll attractor when α = (0.95, 0.95, 0.95), α = (0.9, 0.9, 0.9), α = (0.985, 0.99, 0.98) and α = (0.8, 1, 0.9), in case of (a, b, c) = [START_REF] Hilfer | Applications of Fractional Calculus in Physics[END_REF][START_REF] Pecora | Synchronization in chaotic systems[END_REF][START_REF] Wang | Extending synchronization scheme to chaotic fractional-order Chen systems[END_REF].

In order to observe synchronization behavior in two coupled chaotic fractional order Chen systems with different fractional orders, we build the master and the slave fractional order Chen systems as,

M : ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ d α 1 xm dt α 1 = a(ym -xm), d α 2 ym dt α 2 = (c -a)xm -xmzm + cym, d α 3 zm dt α 3 = xmym -bzm, ( 14
)
S : ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ d β 1 xs dt β 1 = a(ys -xs) + u1(t), d β 2 ys dt β 2 = (c -a)xs -xszs + cys + u2(t), d β 3 zs dt β 3 = xsys -bzs + u3(t), ( 15
)
where α i ≥ β i , for i = 1, 2, 3, subscripts m and s stand for the master system and slave system, respectively, and u(t) = (u 1 (t), u 2 (t), u 3 (t)) is the nonlinear controller to be designed. Our aim is to determine the controller u(t) for the global synchronization of the non-identical fractional order Chen systems ( 14) and [START_REF] Lu | Chaotic dynamics and synchronization of fractional-order Arneodos systems[END_REF]. For this purpose, we define the synchronization error as,

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ e1(t) = xs(t) -xm(t), e2(t) = ys(t) -ym(t), e3(t) = zs(t) -zm(t). ( 16
)
Now, we design the controller u(t) = (u 1 (t), u 2 (t), u 3 (t)) as follows,

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ u1 = a(J α 1 -β 1 -1)[ys -xs] + k1J α 1 -β 1 [xs -xm], u2 = xszs -J α 2 -β 2 [xmzm] + (J α 2 -β 2 -1) [(c -a)xs + cys] + k2J α 2 -β 2 [ys -ym], u3 = J α 3 -β 3 [xmym] -xsys -b(J α 3 -β 3 -1) [zs] + k3J α 3 -β 3 [zs -zm], ( 17
)
where J γ is the Riemann-Liouville fractional integral operator of order γ, defined in (1), k = (k 1 , k 2 , k 3 ) is the coupling matrix. The constants k 1 , k 2 and k 3 will be determined such that e(t) -→ 0 as t -→ +∞.

Applying Laplace transform to systems ( 14), ( 15) and ( 17), letting

X(s) = L((x(t)), Y (s) = L(y(t)), Z(s) = L(z(t)), E 1 (s) = L((e(t)), E 2 (s) = L((e 2 (t)), E 3 (s) = L((e 3 (t))
, and using formulas ( 2) and ( 4), we obtain,

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ s α 1 Xm(s) = s α 1 -1 xm(0) + a(Ym -Xm), s α 2 Ym(s) = s α 2 -1 ym(0) + (c -a)Xm -L{xmzm} + cYm, s α 3 Zm(s) = s α 3 -1 zm(0) + L{xmym} -bZm, ( 18
) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ s β 1 Xs(s) = s β 1 -1 xs(0) + a(Ys -Xs) + U1, s β 2 Ys(s) = s β 2 -1 ys(0) + (c -a)Xs -L{xszs} + cYs + U2, s β 3 Zs(s) = s β 3 -1 zs(0) + L{xsys} -bZs + U3, ( 19
) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ U1 = a(s β 1 -α 1 -1)(Ys -Xs) + k1s β 1 -α 1 (Xs -Xm), U2 = L{xszs} -s β 2 -α 2 L{xmzm} + (s β 2 -α 2 -1) ((c -a)Xs + cYs) + k2s β 2 -α 2 (Ys -Ym), U3 = s β 3 -α 3 L{xmym} -L{xsys} -b(s β 3 -α 3 -1) (Zs) + k3s β 3 -α 3 (Zs -Zm). ( 20
)
Multiplying the first equation in ( 19) by s α1-β1 , the second equation by s α2-β2 and the third equation by s α3-β3 then, from the definition of the error functions ( 16), we get,

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ s α 1 E1(s) = s α 1 -1 e1(0) + a(E2 -E1) + k1E1, s α 2 E2(s) = s α 2 -1 e2(0) + (c -a)E1 + cE2 + k2E2, s α 3 E3(s) = s α 3 -1 e3(0) -bE3 + k3E3. ( 21
)
Applying the inverse Laplace transform, using formula (4), to system (21), we obtain the following linear fractional order system,

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ d α 1 e1 dt α 1 = a(e2 -e1) + k1e1, d α 2 e2 dt α 2 = (c -a)e1 + ce2 + k2e2, d α 3 e3 dt α 3 = -be3 + k3e3. ( 22
)
In case of α i = β i , for i = 1, 2, 3, then the controller u(t), defined in [START_REF] Sheu | Chaos in the Newton-Leipnik system with fractional order[END_REF], reduces to

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ u1 = k1[xs -xm], u2 = xszs -xmzm + k2[ys -ym], u3 = xmym -xsys + k3[zs -zm]. ( 23
)
The Jacobian matrix for the error system ( 22) is

J = ⎛ ⎜ ⎝ -a + k1 a 0 c -a c+ k2 0 0 0 -b + k3 ⎞ ⎟ ⎠ , (24)
and so, the characteristic equation det diag(λ r1 , λ r2 , λ r3 ) -J = 0 can be written as,

((λ r 1 + a -k1)(λ r 2 -c -k2) + a(a -c))(λ r 3 + b -k3) = 0, ( 25
)
where r i = mα i , for i = 1, 2, 3. According to the stability results, the drive system ( 14) and the response system [START_REF] Lu | Chaotic dynamics and synchronization of fractional-order Arneodos systems[END_REF] are synchronized ( e(t) -→ 0 as t -→ +∞) if all roots λ of Eq. ( 25) satisfy |arg(λ)| > γπ/2, where γ = 1/m. Now, in case of α 1 = α 2 = α 3 = α and (a, b, c) = (35, 3, 28), then Eq. ( 25) reduces to, ((λ+35-k 1 )(λ-28-k 2 )+245)(λ+3-k 3 ) = 0, [START_REF] Li | The synchronization of three fractional differential systems[END_REF] and so, systems ( 14) and ( 15) are synchronized if k 1 , k 2 and k 3 satisfy the laws,

⎧ ⎨ ⎩ k1 + k2 -7 ∓ w < 0, i f A= w ∈ IR w k 1 +k 2 -7 < tan(απ/2), if A = iw ∈ iIR (27) k3 < 3, ( 28
)
where A = (7-k1 -k2) 2 27) and [START_REF] Wang | Extending synchronization scheme to chaotic fractional-order Chen systems[END_REF].

Therefore, under the controller,

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ u1 = 35(J 0.05 -1)[ys -xs] + k1J 0.05 [xs -xm], u2 = xszs -J 0.05 [xmzm] + (J 0.05 -1) [-7xs + 28ys] + k2J 0.05 [ys -ym], u3 = J 0.05 [xmym] -xsys -3(J 0.05 -1) [zs] + k3J 0.05 [zs -zm], ( 29
)
the drive system ( 14) and the response system (15) are synchronized. The error functions evolution, in this case, is shown in Fig. 1. From Fig. 1, it is obvious that the components of the error system [START_REF] Zhou | Chaos control and synchronization in a fractional neuron network system[END_REF] decay towards zero as t -→ +∞. So, we can numerically conclude that the designed controller can effectively control the chaotic fractional order Chen systems ( 14) and ( 15) with non-identical orders to achieve synchronization.

Example 2 If we take (a, b, c) = [START_REF] Hilfer | Applications of Fractional Calculus in Physics[END_REF][START_REF] Pecora | Synchronization in chaotic systems[END_REF][START_REF] Wang | Extending synchronization scheme to chaotic fractional-order Chen systems[END_REF], α = (1, 0.9, 0.9), β = (0.9, 0.9, 0.9), k 1 = 35, k 2 = -28 and k 3 = 2, then Eq. ( 25) reduces to,

(λ 19 + 245)(λ 9 + 1) = 0. ( 30
)
Simply, we can show that all roots of Eq. ( 30) lie in the region |arg(λ)| > π/20. Therefore, under the controller

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ u1 = 35(J 0.1 -1)[ys -xs] + k1J 0.1 [xs -xm], u2 = xszs -xmzm + k2(ys -ym), u3 = xmym -xsys + k3(zs -zm), ( 31
)
the drive system ( 14) and the response system [START_REF] Lu | Chaotic dynamics and synchronization of fractional-order Arneodos systems[END_REF] are synchronized. The error functions evolution, in this case, is shown in Fig. 2. It is clear, from Fig. 2, that the components of the error system [START_REF] Zhou | Chaos control and synchronization in a fractional neuron network system[END_REF] decay towards zero as t -→ +∞.

fractional Rössler system

Now, we consider Rössler system [START_REF] Rössler | An equation for continuous chaos[END_REF], which is a nonlinear system that can exhibit one scroll chaotic attractor. Its fractional version reads as,

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ d α1 x dt α1 = -(y + z), d α2 y dt α2 = x + ay, d α3 z dt α3 = z(x -c) + b, ( 32
)
where α = (α 1 , α 2 , α 3 ) is subject to 0 < α 1 , α 2 , α 3 ≤ 1. Integer Rössler system, when α = (1, 1, 1), is chaotic, for example, when (a, b, c) = (0.2, 0.2, 5). Simulations are performed to obtain chaotic behavior of fractional order Rössler system and the results demonstrate that chaos indeed exists with order less than 3. For example, chaotic attractors are found in [START_REF] Li | Chaos and hyperchaos in fractional order Rössler equations[END_REF] when α = (0.9, 0.9, 0.9) and (a, b, c) = (0.4, 0.2, 10). Also, in [START_REF] Peng | Synchronization of fractional order chaotic systems[END_REF], it is found that for the parameters α = (1, 1, 0.8) and (a, b, c) = (0.2, 0.2, 5) fractional order Rössler system can display chaotic behaviors.

For example, fractional Rössler system, when (a, b, c) = (0.2, 0.2, 5), has two equilibrium points, P1 : (0.00801, -0.04006, 0.04006), P2 : (4.99199, -24.95994, 24.95994). [START_REF] Oldham | The Fractional Calculus[END_REF] The Jacobian matrix of system [START_REF] Zhu | Chaos and synchronization of the fractional-order Chua's system[END_REF], evaluated at the equilibrium point x * = (x * , y * , z * ), is

J R = ⎛ ⎜ ⎝ 0 -1 -1 1 0.2 0 z * 0 x * -5 ⎞ ⎟ ⎠ (34)
Now, according to stability results and the results presented in [START_REF] Tavazoei | Chaotic attractors in incommensurate fractional order systems[END_REF], a necessary condition for fractional Rössler system to exhibit chaotic attractor similar to its integer order counterpart is:

min i {|arg(λ i )|} ≤ γπ/2, (35)
where γ = 1/m and λ i 's are the roots of the equation det diag(λ mα1 , λ mα2 , λ mα3 ) -J R | x * = 0, for every equilibrium point x * . Otherwise, one of these equilibrium points becomes asymptotically stable and attracts the nearby trajectories.

For α = (1, 1, 0.8) and according to the first equilibrium point P 1 , the equation det diag(λ mα1 , λ mα2 , λ mα3 ) -J R | x * = 0 becomes, λ 14 + 4.99199λ 10 -0.2λ 9 -0.95833λ 5 + λ 4 + 4.98397 = 0, [START_REF] Gorenflo | Fractional calculus: Integral and differential equations of fractional order[END_REF] and so we get,

min i {|arg(λ i )|} = 0.29486 < π/10. ( 37
)
Hence, system (32) satisfies the necessary condition for exhibiting a one scroll attractor when α = (1, 1, 0.8) and (a, b, c) = (0.2, 0.2, 5). For α = (0.9, 0.9, 0.9) and (a, b, c) = (0.4, 0.2, 10), then the equilibrium point (0.00801, -0.02002, 0.02002) of system (32) has the following eigenvalues λ 1 = -9.99001, λ 2,3 = 0.199008 ∓ 0.979690 i, and so, |arg(λ 2,3 )| = 1.370390 < (0.9)π/2. Therefore, system (32) satisfies the necessary condition for exhibiting one scroll attractor when α = (0.9, 0.9, 0.9) and (a, b, c) = (0.4, 0.2, 10).

In order to observe synchronization behavior in two coupled chaotic fractional order Rössler systems with different fractional orders, we build the master and the slave fractional order Rössler systems as,

M : ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ d α 1 xm dt α 1 = -(ym + zm), d α 2 ym dt α 2 = xm + aym, d α 3 zm dt α 3 = zm(xm -c) + b, (38)
S :

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ d β 1 xs dt β 1 = -(ys + zs) + u1(t), d β 2 ys dt β 2 = xs + ays + u2(t), d β 3 zs dt β 3 = zs(xs -c) + u3(t) + b, ( 39
)
where α i ≥ β i , for i = 1, 2, 3, subscripts m and s stand for the master system and slave system, respectively, and u(t) = (u 1 (t), u 2 (t), u 3 (t)) is the nonlinear controller to be designed. Our aim is to determine the controller u(t) for the global synchronization of non-identical fractional order Rössler systems [START_REF] Matignon | Stability results of fractional differential equations with applications to control processing[END_REF] and [START_REF] Deng | Stability analysis of linear fractional differential system with multiple time delays[END_REF]. For this purpose, we design the controller u(t) as follows,

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ u1 = (J α 1 -β 1 -1)[-xs -zs] + k1J α 1 -β 1 [xs -xm], u2 = (J α 2 -β 2 -1)[xs + ays] +k2J α 2 -β 2 [ys -ym], u3 = J α 3 -β 3 [xmzm] -xszs -c(J α 3 -β 3 -1)[zs] -b + b t α 3 -β 3 Γ(1+α 3 -β 3 ) + k3J α 3 -β 3 [zs -zm], ( 40
)
where J γ is the Riemann-Liouville fractional integral operator of order γ, defined in (1), k = (k 1 , k 2 , k 3 ) is the coupling matrix. Following the same analysis presented in the previous section, if we apply Laplace transform to systems [START_REF] Matignon | Stability results of fractional differential equations with applications to control processing[END_REF], ( 39) and ( 40), using formulas ( 2) and ( 4), we obtain,

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ s α 1 E1(s) = s α 1 -1 e1(0) -(E1 + E3) + k1E1, s α 2 E2(s) = s α 2 -1 e2(0) + E1 + aE2 + k2E2, s α 3 E3(s) = s α 3 -1 e3(0) -cE3 + k3E3. ( 41
)
Applying the inverse Laplace transform, using formula (4), to system (41), we obtain the following linear fractional order system,

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ d α 1 e1 dt α 1 = -(e1 + e3) + k1e1, d α 2 e2 dt α 2 = e1 + ae2 + k2e2, d α 3 e3 dt α 3 = -ce3 + k3e3. ( 42
)
In case of α i = β i , for i = 1, 2, 3, then the controller u(t), defined in [START_REF] El-Sayed | Equilibrium points, stability and numerical solutions of fractionalorder predator-prey and rabies models[END_REF], reduces to

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ u1 = k1[xs -xm], u2 = k2[ys -ym], u3 = xmzm -xszs + k3[zs -zm]. (43)
The Jacobian matrix for the error system ( 42) is

J = ⎛ ⎜ ⎝ -1 + k1 0 -1 1 a + k2 0 0 0 -c + k3 ⎞ ⎟ ⎠ , (44)
and so, the characteristic equation det diag(λ r1 , λ r2 , λ r3 ) -J = 0 can be written as,

(λ r 1 + 1 -k1)(λ r 2 -a -k2)(λ r 3 + c -k3) = 0, ( 45
)
where r i = mα i , for i = 1, 2, 3. According to the stability results, the drive system [START_REF] Matignon | Stability results of fractional differential equations with applications to control processing[END_REF] and the response system (39) are synchronized ( e(t) -→ 0 as t -→ +∞) if all roots λ of Eq. (45) satisfy |arg(λ)| > γπ/2, where γ = 1/m. Now, in case of α 1 = α 2 = α 3 = α, then Eq. (45) reduces to,

(λ + 1 -k 1 )(λ -a -k 2 )(λ + c -k 3 ) = 0, (46)
and so, systems [START_REF] Matignon | Stability results of fractional differential equations with applications to control processing[END_REF] and ( 39) are synchronized if k 1 , k 2 and k 3 satisfy the law,

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ k 1 < 1, k 2 < -a, k 3 < c. ( 47
)
Example 3 Taking (a, b, c) = (0.4, 0.2, 10), α = (1, 1, 1) and β = (0.9, 0.9, 0.9). If we select k 1 = 0.5, k 2 = -1 and k 3 = 8, then they satisfy the law (47).

Therefore, under the controller,

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ u1 = (J 0.1 -1)[-xs -zs] + k1J 0.1 [xs -xm], u2 = (J 0.1 -1)[xs + ays] + k2J 0.1 [ys -ym], u3 = J 0.1 [xmzm] -xszs -c(J 0.1 -1)[zs] -0.2 + 0.2 t 0.1 Γ(1.1) + k3J 0.1 [zs -zm], (48)
the drive system (38) and the response system (39) are synchronized. The error functions evolution, in this case, is shown in Fig. 3. From Fig. 3, it is obvious that the components of the error system (42) decay towards zero as t -→ +∞. So, we can numerically conclude that the designed controller can effectively control the chaotic fractional order Rössler systems ( 38) and ( 39) with non-identical orders to achieve synchronization. (50)

the drive system (38) and the response system (39) are synchronized. The error functions evolution, in this case, is shown in Fig. 4. It is clear, from Fig. 4, that the components of the error system (42) decay towards zero as t -→ +∞.

Conclusion

In the present paper, we study the control and synchronization problems of non-identical chaotic fractional order systems. We present theoretical results for drive-response synchronization between fractional order systems with different fractional orders, based on stability results of fractional order systems and Laplace transform theory, such as the fractional extension of Chen and Rössler systems. The designed adaptive nonlinear controller that applied to the response system affects the system dynamics to realize synchronization. The controller is designed such that the components of the error system decay towards zero as the time variable, t, tends to infinity. The numerical simulations show the effectiveness and the feasibility of the proposed scheme. Finally, the recent appearance of fractional order systems as models in many fields makes it necessary to investigate the issues of chaos, control and synchronization of such systems and we hope that this work is a step in this direction. 14) and [START_REF] Lu | Chaotic dynamics and synchronization of fractional-order Arneodos systems[END_REF], when (a, b, c) = [START_REF] Hilfer | Applications of Fractional Calculus in Physics[END_REF][START_REF] Pecora | Synchronization in chaotic systems[END_REF][START_REF] Wang | Extending synchronization scheme to chaotic fractional-order Chen systems[END_REF], α = (1, 0.9, 0.9), β = (0.9, 0.9, 0.9) and (k1, k2, k3) = (35, -28, 2). β = (0.9, 0.9, 0.9) and (k1, k2, k3) = (0.5, -1, 8).

  Fractional Chen system (6), when (a, b, c) =[START_REF] Hilfer | Applications of Fractional Calculus in Physics[END_REF][START_REF] Pecora | Synchronization in chaotic systems[END_REF][START_REF] Wang | Extending synchronization scheme to chaotic fractional-order Chen systems[END_REF], has three equilibrium points,

4 Now, if α 1 =

 41 α 2 = α 3 = α, then the eigenvalues of the equilibrium points are, ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ P1 : λ1 = -30.83587, λ2 = 23.83587, λ3 = -3, P2 : λ1 = -18.42796, λ2,3 = 4.21398 ± 14.88464 i, P3 : λ1 = -18.42796, λ2,3 = 4.21398 ± 14.88464 i.

- 4 ( 2 . 1

 421 245-(35-k1 )(28+k2)) 1/Example Taking (a, b, c) = (35, 3, 28), α = (1, 1, 1) and β = (0.95, 0.95, 0.95). If we select k 1 = 20, k 2 = -15 and k 3 = 2, then they satisfy the laws (

Example 4

 4 If we take (a, b, c) = (0.2, 0.2, 5), α = (1, 1, 0.8), β = (1, 1, 0.8), k 1 = 0, k 2 = -1.2 and k 3 = 3, then Eq. (45) reduces to, (λ 5 + 1)(λ 5 + 1)(λ 4 + 2) = 0. (49) Simply, we can show that all roots of Eq. (49) lie in the region |arg(λ)| > π/10. Therefore, under the controller ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ u1 = k1[xsxm], u2 = k2[ysym], u3 = xmzmxszs + k3[zszm],

e 1 , e 2 , e 3 Figure 1 :

 131 Figure 1: Synchronization of fractional Chen systems (14) and (15), when (a, b, c) = (35, 3, 28), α = (1, 1, 1), β = (0.95, 0.95, 0.95) and (k1, k2, k3) = (20, -15, 2).

e 1 , e 2 , e 3 Figure 2 :

 132 Figure 2: Synchronization of fractional Chen systems (14) and[START_REF] Lu | Chaotic dynamics and synchronization of fractional-order Arneodos systems[END_REF], when (a, b, c) =[START_REF] Hilfer | Applications of Fractional Calculus in Physics[END_REF][START_REF] Pecora | Synchronization in chaotic systems[END_REF][START_REF] Wang | Extending synchronization scheme to chaotic fractional-order Chen systems[END_REF], α = (1, 0.9, 0.9),

e 1 ,

 1 e 2 , e 3

Figure 3 :

 3 Figure 3: Synchronization of fractional Rössler systems (38) and (39), when (a, b, c) = (0.4, 0.2, 10), α = (1, 1, 1),

e 1 , e 2 , e 3 Figure 4 :

 134 Figure 4: Synchronization of fractional Rössler systems (38) and (39), when (a, b, c) = (0.2, 0.2, 5), α = (1, 1, 0.8), β = (1, 1, 0.8) and (k1, k2, k3) = (0, -1.2, 3).

NODY9609_source.tex; 27/10/2009; 17:30 p. 4

NODY9609_source.tex; 27/10/2009; 17:30 p. 6

NODY9609_source.tex; 27/10/2009; 17:30 p. 7


Marius-F Danca
email: danca@rist.ro

Chaotic behavior of a class of discontinuous dynamical systems of fractional-order

Keywords: Mathematics Subject Classification (2000) 34A36, 26A33, 65P20 Fractional derivative, discontinuous dynamical system, Filippov regularization, differential inclusion, numerical method

Chaotic behavior of a class of discontinuous dynamical systems of fractional-order

Marius-F.

Introduction

Dynamical systems (d.s.), discontinuous with respect to the state variable, occur in a large number of problems, especially from mechanics (dry friction with stick and slip modes, impacts, oscillating systems with viscous damping, elasto-plasticity, alternatively forced vibrations, braking processes with locking phases), elec-trical engineering (electrical circuits and networks with switches, power electronics), automatic and optimal control theory, convex optimizations, game theory, nonsmooth control systems synthesis, uncertain systems, walking machines, biological and physiological systems and everywhere non-smooth characteristics are used to represent switches (see e.g. [START_REF] Andronov | Theory of oscillators[END_REF][START_REF] Buhite | An ordinary differential equation from the theory of plasticity[END_REF][START_REF] Clarke | Optimization and non-smooth analysis[END_REF][START_REF] Deimling | Multivalued differential equations and dry friction problems[END_REF][START_REF] Schilling | An algorithm to solve boundary value problem for differential equations and applications in optimal control[END_REF][START_REF] Wiercigroch | Applied nonlinear dynamics and chaos of mechanical systems with discontinuities[END_REF] and some references therein). In other words, in the real word non-smoothness is a common phenomenon. The underlying mathematical models are generally described by a set of first-order differential equations with discontinuous components. Particularly, the discontinuity is due to the discontinuity of the state variable, of the associated vector field, of the Jacobian, or to higher-order discontinuity.

In practical examples, the discontinuity appears because of piecewise continuous or switch-like functions (e.g. signum, absolute value, Heaviside function -also known as the "unit step function"-maximum function, etc.).

On the other hand, there are many continuous d.s. which display fractional-order dynamics. These mathematical phenomena allow to describe a real object more accurately than do classical "integer" dynamics. One of the main reasons to use the integer-order models was the absence of methods for fractional differential equations.

The real objects or phenomena such as viscoelastic systems, dielectric polarization, electrochemistry, percolation, material science, polymer modeling, theory of ultra-slow processes, electromagnetic waves, evolution of complex systems, financial processes, secure communication, etc. are generally fractional (see [START_REF] Bagley | Fractional order state equations for the control of viscoelastically damped structures[END_REF][START_REF] Nakagava | Basic characteristics of a fractance device[END_REF][START_REF] Oustaloup | La Derivation Non Entiere: Theorie[END_REF][START_REF] Sun | Linear approximation of transfer function with a pole of fractional order[END_REF][START_REF] Podlubny | Fractional Differential Equations[END_REF][START_REF] Ichise | An analog simulation of noninteger order transfer functions for analysis of electrode process[END_REF][START_REF] Podlubny | Analogue realization of fractional-order controllers[END_REF][START_REF] Laskin | Fractional market dynamics[END_REF][START_REF] Kusnezov | Quantum Levy processes and fractional kinetics[END_REF]).

We consider in this paper a class of unified Filippov discontinuous d.s. of fractional-order and prove that the underlying Initial Value Problem (I.V.P.) admits solutions which can be numerically determined. We also prove numerically that in these systems of lower than third order chaos may appear (as it is known, in the case of integer order, according to the well known Poincar-Bendixon theorem, chaos appears only at systems of minimum order three).

These systems are modeled by the following autonomous I.V.P.

d q x(t) dt q = g(x(t))+A s(x(t)), x(0) = x 0 , t ∈ I = [0, ∞), (1)
with g : R n → R n , A = (b i,j ) n n real constant matrix, the vector function s : R n → R n being composed of signum functions (the most encountered case), i.e.

s(x) = ⎛ ⎜ ⎝ sgn(x 1 ) . . . sgn(x n ) ⎞ ⎟ ⎠ , ( 2
)
q the fractional order, 0 < q ≤ 1 (although there are also some applications for q > 1, as is the case in many applications).

The right-hand side depends on p ∈ R, the control parameter.

Hereafter, we impose the following assumption on the discontinuous right-hand functions (H1):The right-hand side is piecewise continuous (continuous on a finite number of m open domains

D i ⊂ R n , i = 1, 2, . . .

, m, and has finite (possibly different) limits from different boundary points, i.e. bounded discontinuities), g being continuous (linear function in the great majority of practical examples). The null set of the discontinuity (switch) points is

M = R n \ ∪ D i .
At every M point, the discontinuity surface is a hyperplane (plane for n = 3).

The considered class of systems is autonomous. Thus, the I.V.P. (1) can be written as follows

d q x dt q = g(x) + A s(x), x(0) = x 0 , t ∈ I. ( 3
)
To achieve our goal, we consider the I.V.P. (3) in a more familiar framework, such as differential inclusions and differential equations of fractional-order. The paper is organized as follows: In Section 2, the discontinuous I.V.P. (3) is transformed, via Filippov's regularization, into a set-valued I.V.P. and then into a continuous differential equation of fractional-order via Cellina's Theorem [START_REF] Aubin | Differential inclusions set-valued maps and viability theory[END_REF][START_REF] Aubin | Set-valued analysis[END_REF]. In Section 3 the numerical integration of differential equations of fractional-order, which will be used to integrate the obtained I.V.P. (3), is presented briefly. In Section 4 this algorithm is used to investigate the lowest-order q for which chaos vanishes in two representative chaotic systems.

Continuous approximation of I.V.P. (3)

One of the first questions of set-valued analysis is how to relate set-valued and single-valued functions, so as to avoid dealing with the set-valued functions. Because the considered examples in this paper are modeled by three-dimensional real I.V.P., the next results will be particularized in the Euclidean space R n , some of them are valid even in more general (such as Hilbert) spaces.

First, let us consider the IVP (3) in the common case q = 1 (i.e. discontinuous I.V.P. )

dx dt = g(x) + A s(x), x(0) = x 0 , t ∈ I. ( 4
)
Due to the discontinuity of the right-hand side, the I.V.P. may not have any solution (see [START_REF] Filippov | Differential Equations with Discontinuous Right-Hand Sides[END_REF] for the background on existence and uniqueness).

To overcome this situation, Filippov proposed the idea of restarting the I.V.P. as a set-valued I.V.P. via differential inclusion [START_REF] Filippov | Differential Equations with Discontinuous Right-Hand Sides[END_REF]

dx dt ∈ F (x), x(0) = x 0 , for a.a. t ∈ I, ( 5
)
where F : R n ⇒ R n is a convex set-valued vector function defined on the set of all subsets of R n (for the background on set-valued functions we refer e.g. to [START_REF] Aubin | Differential inclusions set-valued maps and viability theory[END_REF][START_REF] Aubin | Set-valued analysis[END_REF]).

Let us consider the following discontinuous I.V.P. enjoying the assumptions H1

dx dt = f (x), x(0) = x 0 , t ∈ I.
One of the simplest definitions for F is the following

F (x) = ε>0 μ(M)=0 conv (f {x ∈ R n : x < ε}\M ), (6)
F being the smallest closed and convex set containing all limit values of f , and μ the Lebesgue measure. On D i , where the function f is continuous, F (x) consists of one point, which coincides with the value of f at this point, F (x) = {f (x)}. Thus, on D i , the I.V.P. becomes a single-valued continuous problem dx/dt = f (x). At the discontinuity points, the set F (x) is a subset of R n given e.g. by [START_REF] Buhite | An ordinary differential equation from the theory of plasticity[END_REF].

Remark 1 In the practical examples, ε must be small enough, so that the motion of the physical system can be arbitrarily close to a certain solution of the differential inclusion. For example, the set-valued version of the usual sign function is defined by By applying the Filippov regularization to the I.V.P. ( 4), it becomes a set-valued I.V.P.

Sgn(x) = ⎧ ⎨ ⎩ {+1} x > 0, [-1, 1] x = 0, {-1} x < 0.
dx dt ∈ g(x) + A S(x), x(0) = x 0 , f or a.a. t ∈ I, ( 7
)
where S is the set-valued form of s (2), i.e. S(x) = (Sgn(x 1 ), Sgn(x 2 ), . . . , Sgn(x n ))

t Now, the obtained set-valued I.V.P. admits at least a generalized (Filippov) solution (for existence and uniqueness we refer to [START_REF] Filippov | Differential Equations with Discontinuous Right-Hand Sides[END_REF] and for numerical methods for differential inclusions [START_REF] Taubert | Converging multistep Methods for Initial Value Problems Involving Multivalued Maps[END_REF][START_REF] Dontchev | Difference Methods for Differentiale Inclusions[END_REF]).

Remark 2 It is easy to see that the set-valued function AS is upper semicontinuous with closed and convex values (see e.g. [START_REF] Filippov | Differential Equations with Discontinuous Right-Hand Sides[END_REF] p. 85 or [START_REF] Aubin | Differential inclusions set-valued maps and viability theory[END_REF] p. 101).

Next, following the Filippov procedure applied to the right-hand side of the I.V.P. (3), one obtains the following differential inclusion of fractional-order

d q x dt q ∈ g(x) + A S(x), x(0) = x 0 , for a.a. t ∈ I, ( 8
)
Definition 1 [START_REF] Aubin | Set-valued analysis[END_REF], [START_REF] Kastner-Maresch | Difference methods with selection strategies for differential inclusions[END_REF] selection of a set-valued function

F : R n ⇒ R n is a single-valued function f : R n -→ R n satisfying h(x) ∈ F (x), ∀x ∈ R n .
Convex set-valued functions admit approximate selections [START_REF] Aubin | Set-valued analysis[END_REF], i.e. single-valued functions h ε : R n -→ R n for which Graph(h ε ) ⊂ Graph(F ) + εB and h ε (x) belongs to the convex hull of the image of F. The constructive proof of Cellina's theorem (see e.g. [START_REF] Aubin | Differential inclusions set-valued maps and viability theory[END_REF]) is a great advantage because it allows the construction of h ε .

There are several strategies to choose selections for a general class of set-valued problem (see e.g. [START_REF] Kastner-Maresch | Difference methods with selection strategies for differential inclusions[END_REF]).

Following the preceding theorem, we can enounce the following main result Proposition 1 Let the I.V.P. ( 3) with g continuous. Then, the I.V.P. can be approximated by the following continuous I.V.P. of fractional-order

d q x dt q = f (x), x(0) = x 0 , t ∈ I, with f : R n -→ R n f (x) = g(x) + h ε (x) for x ∈ M , g(x) + As(x) for x / ∈ M,
where h ε is an approximate selection of the set-valued function AS.

Proof The proof is obvious if we consider that first, the I.V.P. ( 3) is transformed, via Filippov regularization, in the set-valued I.V.P. [START_REF] Deimling | Multivalued differential equations and dry friction problems[END_REF]. Next, Theorem 1 can be applied in the neighborhoods of the discontinuity points.

Remark 3 i) Obviously, ε is a vector and in applications we have chosen the same value for all its components;; ii) As proved in [START_REF] Danca | On a possible approximation of discontinuous dynamical systems[END_REF], this local approximation can be made even with smooth functions;

iii) The proof of the above theorem could also be done via maximal monotonicity of AS which ensures its closedness and convexity (see e.g. [START_REF] Aubin | Differential inclusions set-valued maps and viability theory[END_REF] p.141, Proposition 2).

Because of the amount of necessary calculus for fractional numerical integration, we shall use a simple form for h ε .

To be precise, let us consider one component of As in the I.V.P. (3), the scalar function s j (x j ) = a ij sgn(x j ). In a ε-neighborhood of x j , we have chosen in our numerical experiments for h ε the known sigmoid function, widely used in practical applications. Its scalar variant has the following form ( Fig. 1)

h ε (x) = 2 1 + e -x/ε -1.
To connect continuously h ε in some ε-neighborhood of x j , we shall use the following form

h ε (x j ) = 2a ij 1 + e -xj /ε -a ij . ( 9
)
The influence of ε is underlined in Fig. 2 d and Fig. 2 e, where the visible (for ε = 1E -2) corners, typical for discontinuous d.s., become smooth if the width of the neighborhood enlarges. Therefore, for physical reasons ε needs to be chosen relatively small.

Conclusion

The discontinuity hyperplanes of I.V.P. ( 3) can be replaced in their ε-neighborhoods with continuous surfaces, via approximate selections.

Numerical approach of discontinuous systems of fractional-order

In this section we present summarily a numerical scheme for continuous differential equations of fractional-order which will be used for our I.V.P.(3).

Numerical integration of differential equations of fractional-order

There are several definitions of the fractional differential operator.

The Caputo fractional derivative of order q > 0, where q ∈ (m -1, m] for some m ∈ N for a differential function x : [0, t] -→ R, is the most common tool, widely used in real applications, allowing for interpretation in a physically meaningful way (see e.g. [START_REF] Diethelm | Predictor-Corrector Approach for the Numerical Solution of Fractional Differential Equations[END_REF], [START_REF] Podlubny | Fractional Differential Equations[END_REF] and [START_REF] Diethelm | Analysis of Fractional Differential Equations[END_REF])

D q x = 1 Γ (m -q) t 0 x (m) (s) (t -s) q+1-m ds,
where Γ is the known Gamma function.

To be precise, let us consider in the sequel the continuous autonomous fractional I.V.P.

D q x = f (x), 0 ≤ t ≤ T, x (k) (0) = x (k) 0 , k = 0, 1, . . . , m -1. ( 10
)
Under the continuity condition of the function f , there is at least a solution (see e.g. [START_REF] Diethelm | Predictor-Corrector Approach for the Numerical Solution of Fractional Differential Equations[END_REF] where the most general case of nonautonmous problems can be considered). Because D q has an m-dimensional kernel (m being just the value rounded up to the nearest integer i.e. m = q ), m initial conditions need to be specified. Certainly, for the usual case 0 < q ≤ 1 we have to specify just one condition.

The fractional differential operators are computationally time expensive as compared to their integerorder counterparts. Some numerical methods for [START_REF] Wiercigroch | Applied nonlinear dynamics and chaos of mechanical systems with discontinuities[END_REF], essentially ad hoc methods for very specific types of differential equations are presented e.g. in [START_REF] Podlubny | Fractional Differential Equations[END_REF], [START_REF] Diethelm | An algorithm for the numerical solution of differential equations of fractional order[END_REF][START_REF] Lubich | Discretized fractional calculus[END_REF][START_REF] Shokooh | A comparison of numerical methods applied to a fractional model of damping materials[END_REF].

In agreement with the standard mathematical theory ( [START_REF] Samko | Fractional Integrals and Derivatives: Theory and Applications[END_REF] Sect. 42), the initial conditions for the I.V.P. [START_REF] Wiercigroch | Applied nonlinear dynamics and chaos of mechanical systems with discontinuities[END_REF] should have the form

d q-k dt q-k x (t) | 0+ = b k , k = 1, . . . , m.
However, because in practical applications, these values are frequently not available, and it may not even be clear what their physical meaning is, one can specify the initial conditions in the classical form as they are commonly used in initial value problems with integerorder equations [START_REF] Diethelm | Analysis of Fractional Differential Equations[END_REF].

The numerical scheme to integrate [START_REF] Wiercigroch | Applied nonlinear dynamics and chaos of mechanical systems with discontinuities[END_REF], which is used in this paper, is a predictor-corrector algorithm (called hereafter DFF ), proposed by Diethlem, Ford and Freed in [START_REF] Diethelm | Predictor-Corrector Approach for the Numerical Solution of Fractional Differential Equations[END_REF] together with error estimates and numerical examples. The scheme is a generalization of the classical multistep method Adams-Bashforth-Moulton and uses the Caputo derivative.

Being of practical significance and helpful for solving a broad class of problems, the method has been constructed and analyzed for the fully general set of equations without any special assumptions, and is easy to implement on a computer Let us next assume that we are working on a uniform grid {t n = nh : n = 0, 1, ..., N } with some integer N with the step size h := T /N and the case of q ∈ (0, 1].

The predictor formula for the value at the point t n+1 , x p , is the fractional variant of the Adams-Bashforth method which for our particular case m = 1, has the expression

x p (t n+1 ) = x 0 + 1 Γ (q) n j=0 b j,n+1 f (x (t j )) ,
while the corrector formula (the fractional variant of the one-step implicit Adams Moulton method)

x (t n+1 ) = x 0 + h q Γ (q + 2) f (x p (t n+1 )) + h q Γ (q + 2) n j=0 a j,n+1 f (x (t j )) ,
where Γ is the Gamma function, and a and b are the corrector and predictor weights respectively given by the following formula

a j,n+1 = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ n q+1 -(n -q) (n + 1) q if j = 0, (n -j + 2) q+1 + (n -j) q+1 -2 (n -j + 1) q+1 if 1 ≤ j ≤ n, 1 i f j = n + 1, (11)
Table 1 Pseudo-code of DFF numerical scheme

x i [0] = x i0 [0], i = 1, 2, .
. . , Ne {initial conditions} for j := 1 to N do begin for i := 1 to Ne do {predictor}

p i := x i0 [0] + h q Γ (q + 1) j-1 k=0 b [j -k] f i (x 1 [k] , . . . , x Ne [k])
for i := 1 to Ne do {corrector}

x i [j] := x i 0 [0] + h q Γ (q + 2) f i (p 1 , . . . , p Ne ) + ((j -1) q+1 -(j -1 -q) j q )f (x 1 [0] , . . . , x Ne [0]) + j-1 k=1 a [j -k] f (x 1 [k] , . . . , x Ne [k]) end and b j,n+1 = h q q ((n + 1 -j) q -(n -j) q ) . ( 12
)
After a few modifications, designed to enhance the efficiency [START_REF] Diethelm | Predictor-Corrector Approach for the Numerical Solution of Fractional Differential Equations[END_REF], DFF scheme has the form presented as pseudo-code in Table 1 for the case m = 1 and a set of N e autonomous differential equations of fractionalorder q with the right-hand side (f 1 (x 1 , x 2 , . . . , x Ne ) , . . . , f Ne (x 1 , x 2 , . . . , x Ne )) t where:

x i0 , i = 1, 2, ..., N e are the initial conditions;

x[•] array of N + 1 real numbers containing the approximate solutions: x k ≈ x(kh), x being the exact solution;

p the predicted values; a, b, arrays of 3×(N +1) real numbers which, in this algorithm, are determined with the following variant [START_REF] Diethelm | Predictor-Corrector Approach for the Numerical Solution of Fractional Differential Equations[END_REF].

for k := 1 to N do a[k] := k q -(k -1) q b[k] := (k + 1) q+1 -2k k+1 + (k -1) q+1 end
The Gamma function Γ can be approximated via the so-called Lanczos approximation [START_REF] Press | Numerical Recipes in C, The Art of Scientific Computing[END_REF] using the formula (http://www.rskey.org/gamma.htm)

Γ (z) = 6 i=0 P i z i 6 i=0 (z + i) (z + 5.5) z+0.5 e -(z+5.5) , ( 13
)
for a complex variable z for which Re(z) > 0. The coefficients P are drawn in Table 2 Remark 4 i) Another way to deal with fractional derivatives is the use of the frequency domain approximation, based on the approximation of the fractional operators in the frequency domain. This technique, proposed by researchers on automatic control, is suitable for Simulink (Matlab) (see e.g. [START_REF] Charef | Fractal system as represented by singularity function[END_REF][START_REF] Petráš | Chaos in the fractional-order Volta's system: modeling and simulation[END_REF][START_REF] Ahmad | Chaos in fractional-order autonomous nonlinear systems[END_REF]); ii) As known, an attractor requires a relatively large integration time interval (for example for some of our experiments one needs N = 2E4÷2.5E4), and an array of this size is hard to implement in any compilers. The best solution we found for our examples was to save x into a file. Proof The I.V.P. (3) can be transformed into a setvalued I.V.P., via Filippovs regularization (Section 2) which can be continuously approximated (Proposition 1). The obtained single-valued I.V.P. may be numerically integrated via DFF (Section 3.1).

One of the most unpleasant aspects of the DFF scheme is that it requires each step, the entire backward inte-Fig. 2 Phase plots for fractional discontinuous Sprott system (14). a),b) Chaotic trajectories for q = 1 and q = 0.9. c) Attractive fixed point for q = 0.89. d) A chaotic trajectory for q = 0.95 and ε = 0.01 e) A chaotic trajectory for q = 0.95 and ε = 0.1.

gration history. In other words, each value x k is a determined function of all previous values x 0 , x 1 , . . . , x k-1

or, for a large k, this implies supplementary computational cost. Thus, unlike the classical derivatives of integer order, the fractional derivatives operators are not local (i.e. we cannot determine D q x(t) by using only a few values of x in a neighborhood of t).

Remark 5 Lipschitz continuity of g is necessary only to ensure the uniqueness either of the I.V.P. [START_REF] Clarke | Optimization and non-smooth analysis[END_REF] or (3) if the study of convergence properties and errors of some numerical method is required.

Error analysis is a difficult task, especially because of the numerical integration of fractional equations (see a detailed study in [START_REF] Diethelm | Analysis of Fractional Differential Equations[END_REF]). Therefore, in our study we have chosen empirically h so that, for the common case q = 1, the obtained attractor fits as well possible the known attractor shape in the phase space R 3 . Also, the optimal choice of the step length must ensure the maximum accuracy in the approximate solution at minimum computational cost.

In order to avoid the specific phenomena for discontinuous d.s. in the discontinuity surfaces (e.g. sliding modes [START_REF] Wiercigroch | Applied nonlinear dynamics and chaos of mechanical systems with discontinuities[END_REF]), the initial conditions need to be chosen in D i . As can be seen from the obtained images, the trajectories present some "corners" typical of discontinuous d.s.

Examples

Let us first introduce the following notion Notation Let us denote by d = 3q the order of the d.s. modeled by I.V.P. 3 and by d * the lowest value for which chaos persists.

While in the case of the continuous d.s. of integerorder chaos can appear only in nonlinear systems with order minimum 3, in nonlinear systems of fractionalorder, this is not the case. For example, in [START_REF] Hartley | Chaos in a fractional order Chua's system[END_REF] has been shown that a Chua circuit of order d as low as 2.7 can produce a chaotic attractor; in [START_REF] Arena | Chaos in a fractional order Duffing system[END_REF] it is shown that a sinusoidally nonautonomous driven Duffing system with d less than 2 can still behave in a chaotic manner; in [START_REF] Wu | Chaos in the fractional-order Lorenz system[END_REF] it is proved that the lowest order d for fractional Lorenz system is d * = 2.97 and in [START_REF] Lu | Chaotic dynamics of the fractional-order Lü system and its synchronization[END_REF] it is proved that the lowest-order chaotic system out of all the found chaotic systems reported in specialized literature, appears in the fractional Lü system.

In this section we explore the chaos disappearance, by searching d * with the algorithm proposed before, in two examples, of known discontinuous three-dimensional d.s. In the simulations, we only vary the fractional derivative order q while the control parameter p is fixed so the underlying d.s. of integer order behaves chaotically. Simulations are performed by lowering q from q = 1 (the integer case) until chaos behavior disappears.

All computations in this paper were performed in double precision arithmetic, N was chosen large enough (of order 10E4) to provide accurate results, the integration step size h of 1E -3 and ε = 1E -3 of order.

Remark 6 In all the practical examples we have found in the scientific texts (either two-dimensional, typical for dry friction problems, or three-dimensional, found generally in many areas of science) there is at most one signum function in each equation of I.V.P. This helps to facilitate the numerical integration (as well known, it is time-consuming to compute the exponential functions).

Sprott system

The first example considered in this paper was proposed by Sprott in [START_REF] Ahmad | Chaos in fractional-order autonomous nonlinear systems[END_REF] and [START_REF] Sprott | A New Class of Chaotic Circuit[END_REF] and it is a representative one for the considered class of d.s.

d q x 1 dt q = x 2 , ( 14
) d q x 2 dt q = x 3 , d q x 3 dt q = -x 1 -x 2 -px 3 + sgn(x 1 ). whith g(x) = ⎛ ⎝ x 2 x 3 -x 1 -x 2 -px 3 ⎞ ⎠ , A = ⎛ ⎝ 0 0 0 0 0 0 1 0 0 ⎞ ⎠ .
In this case, the chaos corresponding to p = 0.5 (see Fig. 2 a for q = 1 and Fig. 2 b for q = 0.9) vanishes at q = 0.89 when d * = 2.57 and when an attractive fixed point is obtained (Fig. 2 c).

Chen system

The next example is the fractional variant of the discontinuous Chen system [START_REF] Aziz-Alaoui | Asymptotic Analysis of a New Piecewise-Linear Chaotic System[END_REF] belonging to a more general class of discontinuous fractional systems

d q x 1 dt q = -p(x 2 -x 1 ), (15)
d q x 2 dt q = sgn(x 1 )(7px 3 ) + 0.7x 2 , d q x 3 dt q = sgn(x 2 )x 1 -0.168x 3 .

For q = 1 the attractor (chaotic for p = 1.18) is presented in Fig. 3 a. In this chase chaos seems to be more "persistent" (see Fig. 3 b,c) since it exists until q = 0.75 (see Fig. 3 d) when a stable limit cycle appears. Thus, d * = 2.25.

Conclusion

In this paper, a continuous approximation, via Filippov regularization, of a class of fractionalorder discontinuous I.V.P. has been considered. The obtained I.V.P. of fractional-order was integrated by using the algorithm proposed by Diethlem, Ford and Freed in [START_REF] Diethelm | Predictor-Corrector Approach for the Numerical Solution of Fractional Differential Equations[END_REF]. Discontinuous fractional-order Sprott and Chen systems have been considered. We have found that chaos still exists in these systems if the order is lower than three. In the cases of continuous fractional systems, the chaotic behavior is affected when the fractional-order decreases. However, in our opinion, choosing the optimal value for q (and implicitly for d), so that the mathematical model (3) can better illustrate the behavior and properties of a physical system, is much more important than finding the lowest d * (a simple reason could be that for d = d * chaos seems to vanish, in possible contradiction with the real system). 3 Phase plots for fractional discontinuous Chen system (15). a)-c) Chaotic trajectories for q = 1, q = 0.9 and q = 0.8 respectively. d) Stable limit cycle for q = 0.75.

Fig. 1

 1 Fig. 1 Graph of the sigmoid function for different ε values.

3. 2 Theorem 2

 22 Numerical integration of I.V.P.[START_REF] Aubin | Set-valued analysis[END_REF] Now, using the results in Section 2 and Subsection 3.1, we can enounce the main result of this paper Let I.V.P. (3) with g continuous. Then the problem can be numerically integrated.

Fig.

  Fig.3Phase plots for fractional discontinuous Chen system (15). a)-c) Chaotic trajectories for q = 1, q = 0.9 and q = 0.8 respectively. d) Stable limit cycle for q = 0.75.

Table 2 P

 2 coefficients for Γ function[START_REF] Oustaloup | La Derivation Non Entiere: Theorie[END_REF]

	P 0 = 75122.6331530
	P 1 = 80916.6278952
	P 2 = 36308.2951477
	P 3 = 8687.24529705
	P 4 = 1168.92649479
	P 5 = 83.8676043424
	P 6 = 2.50662827511

NODY9612_source.tex; 28/10/2009; 12:02 p. 7