An Efficient Nonstandard Finite Difference Scheme for Chaotic Fractional-Order Chen System

In this paper, an efficient nonstandard finite difference scheme for the numerical solution of chaotic fractional-order Chen system is developed. In the new method, an appropriate nonlocal framework in conjunction with the Grünwald-Letnikov approximation are applied for the discretization of fractional differential system. By constructing the discretization with the nonstandard finite difference scheme, high resolution of the system can be obtained, and the numerical instabilities of the nonlinear fractional-order Chen chaotic system can be also addressed to some extent. In addition, a new fractional derivative of the Caputo type is employed in the context of fractional-order Chen system to further decrease the computational complexity in the long-term treatment of fractional model. Numerical simulations demonstrate the applicability, accuracy and efficiency of the developed method.


I. INTRODUCTION
During the last decades fractional calculus has become a powerful tool to describe the dynamics of complex systems in numerous branches of science and engineering [1]- [6]. Fractional-order differential equations, therefore found numerous applications in modelling fractional-order dynamics of systems from interdisciplinary fields, including visco-elasticity, nono-bio mechanics, control and robotics, quantitative finance, electromagnetic waves, biologic science, etc [7]- [12]. The most important advantage of using fractional-order differential equations in these and other applications is their nonlocal property. It is well known that the integer-order differential operator is local, while fractional-order differential operator is nonlocal. This means that the next state of a system depends not only upon its current state but also upon all of its historical states. This makes studying the fractional-order dynamic systems an active area of research.
Chaotic system represents basic dynamical model displaying, under specific choices of their entry parameters, The associate editor coordinating the review of this manuscript and approving it for publication was Mehmet Alper Uslu. the occurrence of a chaotic attractor [13]- [17]. The involvement of fractional derivatives to chaotic systems will bring further requirements on their qualitative analysis. It has been shown that some fractional-order systems, e.g., fractional-order Lorenz system, fractional-order Chua's system, fractional-order Rössler system, and fractional-order Chen system, etc, behave chaotically [18]- [22]. In recent years, a great deal of effort has been paid on developing suitable numerical methods for the solutions of fractionalorder chaotic systems as the chaotic behaviors of these systems are essential in various stabilization and synchronization issues [23]- [27]. Commonly used methods include Grünwald-Letnikov (GL) method, decomposition method, variational iteration scheme, finite difference scheme, and predictor-corrector method, etc [28]- [31]. However, these methods have some limitations, and low accuracy or even wrong results from these methods have been reported in the detection of chaos in the fractional-order nonlinear systems [32], [33]. For the study of fractional-order chaotic systems, decomposition methods do not work as one requires to find behavior of the solutions as time tends to infinity. Another limitation of decomposition methods is the tedious computations of Adomian polynomials [34].
Although difficulties which arise in the decomposition methods can be overcome by the variational iteration method, the term with the fractional derivative in this method is considered as a restricted variation, making the identification of the Lagrange multipliers very inaccurate [35]. The finite difference method has been considered as a standard tool for solving fractional-order systems. However, the solution from this method may become artificially chaotic for nonlinear equations, which yields spurious instabilities for discretization parameters. Although such spurious instabilities and oscillations can commonly be eliminated by applying small grid-sizes, the computational cost related to investigating the long-term treatment of a considered chaotic model is rather expensive [36]. Predictor-corrector method has become a generally accepted tool for solving fractional-order chaotic systems due to its high resolution. However, due to the wellknown nonlocal nature of fractional operators, the computational complexity of this method may grow prohibitively for systems of global dynamics, and as a result, the computational burden limits the application of this method in the prediction of long-terms behaviors of fractional-order chaotic systems [37]- [39]. Therefore, it is necessary to develop robust and stable numerical methods to eliminate instabilities of realistic parameters of the fractional-order chaotic system, and also to test its long-term accuracy and efficiency.
As an alternative to the standard finite difference scheme, nonstandard finite difference scheme has been developed to offset the weakness of some issues related to numerical instabilities which may be found in finite difference scheme. The main advantage of this scheme is that the essential qualitative features of the mathematical model are preserved independently of the chosen step-size. It has been observed that nonstandard finite difference schemes can more easily preserve certain properties and structures obeyed by the original systems, and have better convergence properties when compared with most of the existing methods [40], [41]. Recently, the nonstandard finite difference schemes for the discrete models have been applied in nonlinear dynamical systems such as the dynamics of HIV transmission, linear heat/diffusion transport problems, and a generalized Nagumo reaction-diffusion model, etc [41]- [43]. In particular, the nonstandard finite difference scheme has been successfully developed to investigate the dynamic treatments of fractional-order Bloch and Memristive systems very recently. Nevertheless, the construction of nonstandard finite difference schemes in the predictions of chaotic behavior of fractional-order nonlinear systems is rather limited.
The main contribution of the present study is the efficient construction of a nonstandard finite discretization scheme in conjunction with a new definition of fractional derivative for solving chaotic fractional-order Chen system. By constructing the discretization with the nonstandard finite difference scheme, high resolution of the chaotic fractional-order Chen system can be obtained. More importantly, spurious solutions and numerical instabilities that arise in finite difference method can be eliminated to great extent. On the other hand, in order to circumvent the difficulty of the computational complexity for long-term treatment of chaotic systems, a new definition of fractional derivative of Caputo type [44], instead of the extensively used Riemann-Liouville definition, is employed in the context of fractional-order Chen system. Since the length of the interval of the new derivative can be adaptively chosen depending on the specific problem, the computational burden for the solution of the chaotic system can be significantly decreased, and as a result, the developed method is more efficient than the most used predictor-corrector method in capturing chaotic behavior of the fractional-order Chen system. To the best knowledge of the authors, the construction of nonstandard finite difference scheme for discretization, as well as the replacement of the Riemann-Liouville derivative by a new fractional derivative, in the context of chaotic fractional-order Chen system have not been reported. This work made a new contribution on these regards.
The rest of paper is organized as follows. We begin by briefly recalling some basic concept of nonstandard finite difference scheme and the new definition of fractional derivative of Caputo type in Section II. The construction of nonstandard finite difference scheme to the Grunwald-Letnikov discretization process for solving chaotic fractional-order Chen system, under the new definition of fractional derivative, is presented in Section III. Followed by the numerical simulation of the chaotic fractional-order Chen system in Section IV, to demonstrate the accuracy and efficiency of the developed method. Comparisons of the developed method with the predictor-corrector method are made.

II. PRELIMINARIES
This section gives a brief introduction to the nonstandard finite difference method and the new definition of fractional derivative of Caputo type that required for our subsequent development in this paper. Detailed information can be found, for example, in [40], [44].

A. THE NONSTANDARD FINITE DIFFERENCE SCHEME
The nonstandard finite difference schemes are discrete model of ordinary and partial differential equations that have been used to model important phenomena in science and engineering. The basic idea is to use the denominator function, instead of the step size in conventional numerical schemes, to make the numerical solution as close to the exact solution as possible at the mesh points. According to [40], the nonstandard finite difference scheme is constructed based on a set of rules, which include (1) the orders of the discrete derivatives should be equal to those of the corresponding derivatives of differential equations, (2) the denominator functions for the discrete derivatives must, in general, be expressed in terms of more complicated functions of the step-sizes than those conventionally used, (3) nonlinear terms should be approximated in a nonlocal way, and etc. VOLUME 8, 2020 For a given differential equation with the form where f (y(t)) is the nonlinear term of the differential equation. By applying the nonstandard finite difference scheme, the discrete form of Eq.(1) becomes where φ(h) is the denominator function of the step size h, h = t/n → 0, and φ(h) should satisfy the following condition which is in agreement with Rule 2 mentioned above Examples of function φ(h) that satisfy Eq.(3) are h, sin(h), or e h − 1, etc. Moreover, in order to satisfy Rule 3, the right-hand side of Eq.(1) should be replaced by the following nonlocal discrete approximation [45] where a is the coefficients of the nonlinear term f (y(t)) in Eq.(1).

B. THE NEW DEFINITION OF FRACTIONAL DERIVATIVE
It is known that most of the existing fractional derivatives, ranging form Riemann-Liouville to Caputo definitions, are global in time, i.e., the fractional operators are defined over the entire past history of the analysed process. This means that the interval of integration may be very large with the increase of time t, and as a result, the computational complexity of these fractional derivatives and the resulting fractional differential equations may be prohibitive for systems of global dynamics. In order to circumvent this difficulty, a new fractional derivative of the Caputo type was developed very recently [44]. Let f a differentiable function, then the new derivative of fractional order p is defined as: with the kernel K p (t − s) given as where f (n) denotes the common n-order derivative, and the changeable integral interval l(t, p) is the function of time t and fractional order p, depending on the specific problem at hand.
Obviously, the integral interval in the new derivative can be adaptively chosen so that the computational burden can be significantly decreased. In this sense, the new derivative is no longer global as the existing definitions. In addition, although the new derivative is defined in the framework of the Caputo type, it is naturally associated with the Reimann-Liouville definition. In this way, the well established Grünwald-Letnikov approach for numerically solving the Reimann-Liouville differential differential equations can be readily modified for approximating solutions of fractional differential equations involving new derivatives [44].

III. CONSTRUCTION OF NONSTANDARD FINITE DIFFERENCE SCHEME FOR CHAOTIC FRACTIONAL-ORDER CHEN SYSTEM
In this section, we develop a straightforward method based on the nonstandard finite difference scheme and the Grünwald-Letnikov approximation to obtain numerical solutions of the fractional-order Chen system. In order to further improve the efficiency of the developed method, the new fractional derivative that described in Section II-B is employed in the context of chaotic fractional-order Chen system.

A. THE FRACTIONAL-ORDER CHEN SYSTEM WITH RIEMANN-LIOUVILLE DERIVATIVE
The Chen system, which was introduced in [46], is a chaotic system with a double scroll attractor. The discovery of Chen system reflected the existence of the Lorenz chaotic attractor via a linear state-feedback control used in a chaotification process in a stable Lorenz system. The fractional-order Chen system, which has been derived by replacing standard derivatives with fractional derivative, is described by where (x, y, z) are the state variables, (a, b, c) are positive real numbers, and the symbol d p i dt p i = D p i t (i = 1, 2, 3) stands for the Riemann-Liouville derivative operator of a real order 0 < p 1 , p 2 , p 3 < 1. When p 1 = p 2 = p 3 = 1, system in Eq.(7) reduces to the Chen system, which is chaotic when (a, b, c) = (35, 3, 28).
The Grünwald-Letnikov discritization is firstly used to approximate the Riemann-Liouville fractional derivative based on the finite differences on an uniform grid in [0, t] with step size h as where the Grünwald-Letnikov coefficients ω p j is computed by the following recursive formula Inspired by the nonstandard finite difference scheme, i.e., by replacing the step size h in Eq.(8) with the denominator function φ(h), the Riemann-Liouville fractional derivative is further approximated as Thus, by virtue of Eq.(10), the fractional-order Chen system in Eq. (7) can be discritisized as where φ 1 (h), φ 2 (h), and φ 3 (h) are the denominator functions. Following the method in [36], the denominator functions with the following form are adopted Further, according to the rules in nonstandard finite difference scheme, i.e., nonlinear terms should be approximated in a nonlocal way [45], the nonlinear terms in Eq.(7) are approximated by the following substitutions: a)x, x n+1 z n → xz, (1 + c)y n − y n+1 → cy, 2x n+1 y n − x n y n+1 → xy, By doing some algebraic manipulation to Eq.(11) and Eq. (13), the nonstandard finite difference scheme on the chaotic fractional-order Chen system is finally constructed asBy doing some algebraic manipulation to Eq.(11) and Eq. (13), the nonstandard finite difference scheme on the chaotic fractional-order Chen system is finally constructed as, (14) shown at the bottom of this page. It can be seen that, similar as the finite difference scheme, the developed nonstandard finite difference scheme for the discretization of fractional-order Chen system is straightforward, and can be readily generalized to a wide class of existing fractional-order nonlinear systems, ranging from fractional Lorenz system to fractional-order Chua's system. More importantly, by the replacement of the step size with a more complicated denominator function, as well as a nonlocal approximation of the nonlinear terms in the fractional-order differential equations, spurious solutions and numerical instabilities encountered in the existing finite difference scheme can be eliminated.

B. THE FRACTIONAL-ORDER CHEN SYSTEM WITH NEW DERIVATIVE
In order to improve the computational efficiency in the solution of fractional differential system, the new fractional derivative with changeable integral interval that described in Section II-B, instead of the classical Riemann-Liouville derivative, is employed in the context of fractional-order system. The resulting chaotic fractional-order Chen system has the form where symbol D p i (i = 1, 2, 3) stands for the new fractional derivative operator that has been presented in Section II-B. By following the similar procedures that have been described in Section III-A, the Grünwald-Letnikov approximationbased nonstandard finite difference scheme on the new system given in Eq. (15) can be derived as, (16) shown at the bottom of this page, where l(t, p) denotes the maximum  integer which is less than or equal to l(t, p). Obviously, it is easy to determine the changeable integral length l(t, p), also known as memory length, such that the fractional-order Chen system defined in Eq.(15) coincides with the existing one. In this case, the nonstandard finite difference schemes presented in Eq. (16) reduce to that in Eq. (14). We note that the changeable integral length l(t, p), has to be determined in prior in the implementation of the developed nonstandard finite difference scheme, although in some cases, a flexible choice of l(t, p) can be used to meet the requirement of real world process, depending on the complexity of the problem. In this study, we develop a short memory principle-based method to determine l(t, p) such that the computational complexity for solving Eq.(16) can be decreased to a reasonable level. Since the short memory principle naturally associate the new fractional derivative with the Grünwald-Letnikov method, the set of Grünwald-Letnikov discretization coefficients can be used to determine the changeable integral length l(t, p) as where ω p j is the Grünwald-Letnikov coefficients as given in Eq. (9). With this choice of memory length l(t, p), we define a memory intensity measure, L f , as L f = l(t, p)/t, and it is straightforward to show that 0 ≤ L f ≤ 1. Figure 1 describes the relation between the memory intensity measure, L f , and the value of fractional order with different number of discrete points n. Obviously, the memory length l(t, p) should be adopted dependent on the value of memory intensity L f , i.e., a small value of L f leads to a relatively small memory length, and thereby a smaller amount of computational cost is required in the computation of the resulting fractional differential equation. It is noted that, with the integer-order case, i.e., p = 1, the memory intensity L f becomes zero, illustrating that the no memory effect exists in the system. In order to further investigate the effect of the choice of memory length l(t, p) on the accuracy of the numerical results, the resulting truncation error caused by the memory length l(t, p) is further determined as where M = sup τ ∈(0,t−l(t,p)) |f (y(t))|. In practice, Eq. (18) can be used to measure the truncated error that caused by the local approximation in the solution of fractional-order Chen system. With the truncation error in Eq. (18), one can further derive the required memory length depending on the target accuracy. As a result, the employment of the new derivative supplies more degree of flexibility to implement time-dependent analysis to meet the requirement of a realworld complex process. Due to the changeable nature of the new fractional operator, the computational complexity related to investigating the long-term chaotic behavior of fractionalorder Chen system can be decreased to great extent.

C. STABILITY OF THE DEVELOPED METHOD
In this section, stability properties of the developed nonstandard finite difference schemes in Eq. (14) and Eq. (16) are studied to determine the stability bound of the proposed method. We firstly rewrite Eq. (14) in an explicit form as      x n+1 = A 1 x n + B 1 y n y n+1 = A 2 y n + B 2 x n+1 + C 2 x n+1 z n z n+1 = A 3 z n + B 3 x n+1 y n + C 3 x n y n+1 (19) where, (20) as shown at the bottom of this page. We then consider the state variable x. Assume that x 0 is the only term that has an error with the perturbed valuex 0 = x 0 + δ x 0 . Then, it is need to determine the effect of this error on the next solution x 1 . The perturbed valuex 0 produces a perturbation ofx 1 of the formx 1 = x 1 + δ x 1 . According to Eq. (19), we obtain (20) 98414 VOLUME 8, 2020 Therefore we have δ x 1 = A 1 δ x 0 which means that the error is amplified by the factor A 1 when the nonstandard finite difference scheme advanced by one step. After k steps, we have For the method to be stable, it is necessary to prove that |A 1 | ≤ 1 for all step size h sufficiently small. From Eq. (20), it is not difficult to find that condition |A 1 | ≤ 1 always holds for the chaotic fractional-order Chen system.
For the rest of the state variables y and z, we still assume that y 0 and z 0 are the only term that has an error with the perturbed valueȳ 0 = y 0 + δ y 0 , andz 0 = z 0 + δ z 0 , respectively. Then, from Eq. (19), we derive the following relation after one step, (22) as shown at the bottom of this page.
Following the similar procedure as formulated above, the following condition should be satisfied to keep the method stable for all step size h sufficiently small. We note that, for fractional-order Chen system, it is not difficult to verify that the conditions |A 1 | ≤ 1, |C 2 | ≤ 1, and |C 3 | ≤ 1 in Eq.(23) always hold. Thus, stability condition of the developed nonstandard finite difference scheme can be simplified as |A 2 | ≤ 1 and |B 2 | ≤ 1. Furthermore, stability bound of the step size h can be determined from solving the following inequality as where ϕ 2 (h) is the denominator function as shown in Eq. (12). We emphasize that although stability condition in Eq.(23) is derived based on the nonstandard finite difference scheme in Eq. (14), it is still works for that in Eq.(16) because only Grünwald-Letnikov coefficients ω 11 , ω 21 , and ω 31 are used in the above derivation. Therefore, stability bound of the step size h given in Eq.(24) is valid for both developed nonstandard finite difference schemes in Eq. (14) and Eq.(16).

IV. NUMERICAL SIMULATIONS
In this section, numerical solutions of the chaotic fractionalorder Chen system with the Riemann-Liouville derivative and the new derivative are presented based on the developed nonstandard finite difference schemes in Eq. (14) and Eq.(16), respectively. We fix the two parameters a and b at 35 and 3, and set equal p to the fractional orders p 1 , p 2 , and p 3 in the system. Also, the initial values are set as (x, y, z) = (1, 2, 3), if not specified.

A. LOCAL STABILITY OF FRACTIONAL-ORDER CHEN SYSTEM
It is well known that the fractional-order Chen system in Eq.(7) has three equilibria      a)b, (2c − a)b, 2c − a), if (2c − a)b > 0, and that the Jacobian matrix of system evaluated at equilibrium point (x * , y * , z * ) isif (2c − a)b > 0, and that the Jacobian matrix of system evaluated at equilibrium point (x * , y * , z * ) is Also, according to the stability theorem of linear fractional systems, the necessary condition to fractional-order systems have chaos is [47] π/2M − min where M is the lowest common multiple of the denominators of p i , and λ i are the roots of the following Equation Based on Eq. (27), it is not difficult to show that the fractionalorder Chen system in Eq. (7) is asymptotically stable with no chaos when the parameters are set as (a, b, c) = (35, 3, 23) and p i = 0.9, i = 1, 2, 3, since the instability measure of the system in this case would be, (29) as shown at the bottom of this page, where λ i are the roots of For the case (a, b, c) = (35,3,24) and p i = 0.9, i = 1, 2, 3, Eq.(28) transforms to be and in this case, the relative value for the instability measure is calculated as, (32) as shown at the bottom of this page, illustrating the existence of chaotic attractor for the given parameters in the fractional-order Chen system. In order to illustrate the performance of the developed nonstandard finite difference schemes in Eq. (14) and Eq.(16) (denoted by new method-1 and new method-2, respectively), the numerical results are compared with those from the predictor-corrector method (denoted by PC method) and the finite difference method (denoted by FD method). Figure 2 represents the phase portrait for chaotic solutions under the case (a, b, c) = (35,3,24) and p i = 0.9, i = 1, 2, 3. It can be seen that the chaotic attractors of the fractional-order Chen system are successfully captured by the developed nonstandard finite difference schemes (e.g., both new method-1 and new method-2). The chaotic attractors obtained from different methods are similar illustrating the developed method can achieve comparable accuracy with the well-recognized predictor-corrector method. Phase plots shown in Figure 3 demonstrate that the fractionalorder Chen system is asymptotically stable and converges to the same equilibrium points under the case (a, b, c) = (35,3,23) and p i = 0.9, i = 1, 2, 3 with different methods. Again, it is found that the equilibrium points that obtained from both developed nonstandard finite difference schemes and the predictor-corrector method are almost the same, verifying the high resolution of the developed method.
In order to further examine the accuracy and stability of the developed method, more numerical results about the  , 2c − a), which the fractional-order Chen system converge to, are listed in Table 1 and Table 2. We firstly considered the cases for (a, b, c) = (35,3,23) and (35,3,20) with different values of fractional order, varying from 0.1 to 0.9, respectively. Note that the exact equilibrium points that correspond to these two cases are computed as E c = (±5.7446, ±5.7446, 11) and E c = (±3.8730, ±3.8730, 11), respectively. From Table 1, it can be seen that, with a large fractional order, i.e., p ≥ 0.7, results from different methods are in good accordance with the exact one, illustrating that the dynamical behaviors obtained from the developed nonstandard finite difference schemes are consistent with those from the predictor-corrector method. With the decreasing of the fractional order, i.e., p ≤ 0.4, results from standard finite difference method clearly diverge to the exact one, while  the developed nonstandard finite difference schemes and the predictor-corrector method can still achieve satisfactory accuracy.
We next examined the effect of step size on the stability of the method. According to Eq.(24), the upper bound of step size is computed as 0.01. That is, when h ≤ 0.01, the developed nonstandard finite difference scheme is numerically stable. Table 2 lists the numerical solutions of the equilibrium points of the system with fixed fractional order p = 0.9, and different values of step size, varying from 0.0005 to 0.01. It is obvious that the results from both developed methods and predictor-corrector method are in good accordance with the exact solutions for all cases, while the results from the finite difference method diverge from the exact solution with the increasing of the step size. Figure 4 geometrically describes the relations between the solution of state variables x, y, z and the step size h with different methods. It can be seen that the numerical instability that caused by the finite difference method can be eliminated by the developed nonstandard finite difference schemes to great extent. In order to further verify the stability condition, relation between the solution of state variables x and the step size h from new method-1 is shown in Figure 5. It is found that if the step size allowed by the stability criterion is applied, the solution will converge otherwise oscillation behaviour of the solution will be obtained. Similar results have been also found for other state variables y and z, and also for new method-2, these results are not given in this paper for simplicity.
In order to illustrate efficiency of the developed method, we further compare the CPU time in the solution of fractionalorder Chen system with different values of step size. Table 3 lists the CPU time required by the new method-1, finite difference method and predictor-corrector method with a relatively large step size h = 0.01. It is seen that the CPU time required by the finite difference method and the new method-1, is about 94s, which is much less than that from the predictor-corrector method, 32054.6s, at the fixed time t = 500s, illustrating that the developed nonstandard finite difference scheme is much more efficient than the predictorcorrector method to achieve the comparable accuracy. With a smaller step size h = 0.005, the CPU time required by new method-1, new method-2 and finite difference method are further compared in Table 4. Since the predictor-corrector method is more time-consuming with a smaller step size, we did not list the CPU time from this method. Again, CPU time from new method-1 and finite difference method is very close. It is also found that the memory length l(t, p) in new method-2 is about 30% less than the time interval in new method-1. As a result, more than 30% of the computational cost from the new method-1 can be reduced by employing the changeable memory length l(t, p) in new method-2, demonstrating the advantage of the adoption of new derivative in the developed numerical scheme for the fractional-order Chen system.

V. CONCLUSION
In this paper, a nonstandard finite discretization scheme in conjunction with a new fractional derivative has been successfully developed for the numerical solution of chaotic fractional-order Chen system. Based on the nonstandard finite difference scheme and the Grünwald-Letnikov approximation, the developed method is straightforward to implement and can achieve high resolution of dynamic behavior of system. The stability of the developed scheme is investigated, and it has proved that the nonstandard finite difference scheme is conditionally stable. If the step size allowed by the stability criterion is applied, the solution will converge otherwise oscillation behaviour will be obtained. In addition, a new definition of fractional derivative is employed in the context of fractional-order Chen system to decrease the computational complexity in the long-term treatment of fractional-order Chen system. Numerical simulations show that the developed nonstandard finite discretization method can achieve comparable accuracy with the well-recognized predictor-correction method, while the computational burden of the developed method is much less than that of the predictor-correction method. Additional 30% of computational time can be further reduced by employing the new fractional derivative instead of the Riemann-Liouville one. We point out that the developed method may be extended to some other linear or nonlinear chaotic systems of fractional order.