In the past twenty years, many scientists paid their attention on the fractional-order chaotic dynamical systems (see Genesio-Tesi system [1], Rabinovich system [2], and Lü system [3] et al.). They presented chaotic attractors indeed occur in the fractional-order model with order less than 3. Sheu and Chen [4] found that the lowest order of the fractional-order Newton-Leipnik system is 2.82. In 2004, Li and Peng [5] discovered the rich dynamical behavior displayed in the fractional-order Chen system such as the fixed points, limit cycles, periodic motions, and chaotic motions.

The original Shimizu-Morioka model [6] is described by the following ordinary differential equation: { x ˙ = P ( x , y , z ) = y y ˙ = Q ( x , y , z ) = x − β ˜ y − x z z ˙ = R ( x , y , z ) = − α ˜ z + x 2 (1.1) where ( x , y , z ) ∈ ℝ 3 are the state variables and α ˜ , β ˜ are positive real parameters. This model has been proposed as a simplified and an alternative model for studying the dynamics of the well-known Lorenz system [7] for large Rayleigh numbers (Ra), in which the complex behavior of the trajectories has been discovered by means of computer simulation. As in the Lorenz model, the Shimizu-Morioka model is invariant, with respect to the substitution ( x , y , z ) → ( − x , − y , z ) . The model received much attention due to its stability to describe bifurcation of the associated Lorenz-like strange attractors [8], for example, taking α ˜ = 0.45 and β ˜ = 0.75 (Figure 1).

Intrigued by the above interesting work, many researchers [9, 10] focused their study on the dynamical behavior analysis of the Shimizu-Morioka model. In particular, articles [11, 12] use feedback control laws and the delay feedback control method [13, 14] to study the local and global stabilization and bifurcation of the Shimizu-Morioka chaotic model.

If the dynamical system of Eq. 1.1 follows ∂ P ∂ x + ∂ Q ∂ y + ∂ R ∂ z = − ( α ˜ + β ˜ ) < 0 (1.2) then the system is known to be a dissipative one.

In 1992, British scholar Rucklideg studied two-dimensional convection problems of solute gradients and magnetic fields, and introduced the following chaotic system [15]: { x ⌣ ˙ = − a x ⌣ + b y ⌣ − y ⌣ z ⌣ y ⌣ ˙ = x ⌣ z ⌣ ˙ = − z ⌣ + y ⌣ 2 (1.3) where ( x ⌣ , y ⌣ , z ⌣ ) ∈ ℝ 3 are the state variables and a , b are the positive real parameters. When α ˜ ≠ 0 , by transformation: x = α ˜ 3 2 y ⌣ , y = α ˜ 5 2 x ⌣ , z = α ˜ 2 z ⌣ , t = α ˜ − 1 t ⌣ (1.4) transformed system Eq. 1.1 to { x ⌣ ˙ = − β ˜ α − 1 x ⌣ + α ˜ − 2 y ⌣ − y ⌣ z ⌣ y ⌣ ˙ = x ⌣ z ⌣ ˙ = − z ⌣ + y ⌣ 2 (1.5)

Therefore, system Eq. 1.5 and system Eq. 1.3 are equivalent when a = − β ˜ α ˜ − 1 , b = 1 / α ˜ 2 , and α ˜ ≠ 0 .

In this article, Section 2 provides a brief review of the fractional-order operator and discretization fractional-order Shimizu-Morioka model using numerical algorithm. In Section 3, the complex dynamical behaviors of the Shimizu-Morioka model with a fractional order are studied numerically in four cases. Finally, conclusions are given in Section 4.

In this section, we first give out the fractional-order differential operator and the Shimizu-Morioka model with a fractional order. Furthermore, we use the predictor–correctors scheme to discrete the fractional-order Shimizu-Morioka model. Last, we discuss the necessary condition for the existence of chaotic attractors.

There are several definitions of the fractional differential and integral operator, including Grünwald-Letnikov operator, Riemann-Liouville operator, and Caputo operator [16–18]. In this study, we use the following Caputo-type fractional derivative [19]. D ∗ α y ( x ) = J m − α y m ( x ) ,   α > 0 (2.6) where m = [ α ] is the first integer which is not less than α, y m is the ordinary m-order derivative, J β is the β-order Riemann-Liouville integral operator defined by J β z ( x ) = 1 Γ ( β ) ∫ 0 x ( x − t ) β − 1 z ( t ) d t ,   β > 0 (2.7) where Γ ( β ) is the gamma function.

The classical Riemann-Liouville fractional derivative is defined by D α y ( x ) = d m d x m J m − α y m ( x ) (2.8) which requires the homogeneous initial conditions. The main reason why we chose the Caputo-type fractional derivative is that the inhomogeneous initial conditions are also permitted.

The integer-order Shimizu-Morioka model Eq. 1.1 has been extended to the fractional-order Shimizu-Morioka model, which could describe the memory and hereditary properties of the model better. The fractional-order Shimizu-Morioka model is described as follows—in which the standard derivative will be replaced by the fractional-order derivative. { d q 1 x d t q 1 = y d q 2 y d t q 2 = x − β ˜ y − x z d q 3 z d t q 3 = − α ˜ z + x 2 (2.9) where 0 < q 1 , q 2 , q 3 ≤ 1 and the order is denoted by q = ( q 1 , q 2 , q 3 ) .

As for model Eq. 2.9, we derive the predictor–correctors scheme which is the generation of Adamas-Bashforth-Moulton one [16, 20]. The following fractional-order differential equation D ∗ α y ( t ) = f [ t , y ( t ) ] , 0 ≤ t ≤ T y ( k ) ( 0 ) = y 0 ( k ) , k = 0,1 , ⋯ , m − 1 (2.10) is equivalent to the Volterra integral equation y ( t ) = ∑ k = 0 [ α ] − 1 t k k ! y 0 ( k ) + 1 Γ ( α ) ∫ 0 t ( t − s ) α − 1 f [ s , y ( s ) ] d s (2.11)

Set h = T / N , t n = n h , and n = 0,1 , … , N ∈ ℤ + , then (2) can be discretized as follows: y n ( t n + 1 ) = ∑ k = 0 [ α ] − 1 t n + 1 k k ! y 0 ( k ) + h α Γ ( α + 2 ) { f [ t n + 1 , y h p ( t n + 1 ) ] + ∑ j = 0 n a j , n + 1 f [ t j , y h ( t j ) ] } where y h p ( t n + 1 ) = ∑ k = 0 [ α ] − 1 t n + 1 k k ! y 0 ( k ) + 1 Γ ( α ) ∑ j = 0 n h α α [ ( n − j + 1 ) α − ( n − j ) α ] f [ t j , y h ( t j ) ] a j , n + 1 = { n α + 1 − ( n − α ) ( n + 1 ) α ,                   j = 0 ( n − j + 2 ) α + 1 + ( n − j ) α + 1 − 2 ( n − j + 1 ) α + 1 ,     1 ≤ j ≤ n 1 ,                                     j = n + 1 The error estimate is max j = 0,1 , … , N | y ( t j ) − y h ( t j ) | = O ( h p ) , p = min ( 2,1 + α ) .

Applying the above formula, system Eq. 2.9 can be discretized as follows: { x n + 1 = x 0 + h q 1 Γ ( q 1 + 2 ) ( y n + 1 p + ∑ j = 0 n α 1 , j , n + 1 y j ) y n + 1 = y 0 + h q 2 Γ ( q 2 + 2 ) [ x n + 1 p − β ˜ y n + 1 p − x n + 1 p z n + 1 p + ∑ j = 0 n α 2 , j , n + 1 ( x j − β ˜ y j − x j z j ) ] z n + 1 = z 0 + h q 3 Γ ( q 3 + 2 ) [ − α ˜ z n + 1 p + ( x n + 1 p ) 2 + ∑ j = 0 n α 3 , j , n + 1 ( − α ˜ z j + x j 2 ) ] (2.12) where { x n + 1 p = x 0 + 1 Γ ( q 1 ) ∑ j = 0 n h q 1 q 1 [ ( n − j + 1 ) q 1 − ( n − j ) q 1 ] y j y n + 1 p = y 0 + 1 Γ ( q 2 ) ∑ j = 0 n h q 2 q 2 [ ( n − j + 1 ) q 2 − ( n − j ) q 2 ] ( x j − β y j − x j z j ) z n + 1 p = z 0 + 1 Γ ( q 3 ) ∑ j = 0 n h q 3 q 3 [ ( n − j + 1 ) q 3 − ( n − j ) q 3 ] ( − α z j + x j 2 ) and α i , j , n + 1 = { n q i + 1 − ( n − q i ) ( n + 1 ) q i ,                   j = 0 ( n − j + 2 ) q i + 1 + ( n − j ) q i + 1 − 2 ( n − j + 1 ) q i + 1 ,   1 ≤ j ≤ n ,   i = 1,2,3 1 ,                                   j = n + 1

The fractional-order Shimizu-Morioka model of system Eq. 2.9 discretes to system Eq. 2.12.

Now, we discuss the necessary condition for the existence of chaotic attractors in the fractional-order Shimizu-Morioka model. Set d q 1 x / d t q 1 = 0 , d q 2 y / d t q 2 = 0 , d q 3 z / d t q 3 = 0 , we get the following equilibrium points of system Eq. 2.9. E 0 = ( 0,0,0 ) , E 1 = ( α ˜ , 0,1 )

The Jacobian matrices at the equilibrium points E 0 and E 1 are J ( E 0 ) = [ 0 1 0 1 − β ˜ 0 0 0 − α ˜ ] , J ( E 1 ) = [ 0 1 0 1 − β ˜ − α ˜ 0 0 − α ˜ ]

The eigenvalues at E 0 are λ 1 = − 0.45 , λ 2 = − 1.443 , and λ 3 = 0.693 , and the eigenvalues at E 1 are λ 1 = 0.1061 + 0.7912 i , λ 2 = − 1.412 , and λ 3 = 0.1061 + 0.7912 i . E 0 and E 1 are saddle points.

Suppose λ is the unstable eigenvalue of the saddle points, then the necessary condition for the fractional-order system Eq. 2.9 to remain chaotic is keeping the eigenvalue λ in the unstable region. By [21], if the eigenvalue λ is in the unstable region, then the following condition is satisfied. | arg   λ i | > q π 2 where | arg   λ i | denotes the argument of the eigenvalue λ. That is, q > 2 π | arg   λ i |

So, the necessary condition for the existence of chaotic attractors in the fractional-order system Eq. 2.9 is q > 2 π | arg   λ i | = 2 π arctan 0.7912 0.1061 = 0.9151 which implied that when q i ( i = 1,2,3 ) > 0.9151 , system Eq. 2.9 has chaos, and when q i ( i = 1,2,3 ) < 0.9151 , system Eq. 2.9 has no chaos.

In what follows, some numerical simulations of system Eq. 2.12 will be studied. We chose the parameters α ˜ = 0.45 , β ˜ = 0.75 , and the initial value ( x 0 , y 0 , z 0 ) = ( 1,1,2 ) . The phase portraits and time histories are used to research the dynamical behaviors of system Eq. 2.9. Four cases are considered as follows.

This article mainly discussed the dynamical behaviors of the fractional-order Shimizu-Morioka model. We find that chaos does exist in the fractional-order model with order less than 3. Future work that requires further consideration regarding this topic includes theoretical analysis of system Eq. 2.9, the largest Lyapunov exponent in the state space, the linear and nonlinear feedback controller, synchronization of this kind of system, and in-depth studies on chaos control for the fractional state.