Bifurcation analysis of fractional duffing system based on improved short memory principle method

. In this paper, the dynamic characteristics of fractional Duffing system are analyzed and studied by using the improved short memory principle method. This method has small amount of calculation and high precision, and can effectively improve the problem of large amount of calculation caused by the memory of fractional order. The influence of frequency change on the dynamic performance of the fractional Duffing system is studied using nonlinear dynamic analysis methods, such as Phase Portrait, Poincare Map and Bifurcation Diagram. Moreover, the dynamic behaviour of the fractional Duffing system when the fractional order and excitation amplitude changes are investigated. The analysis shows that when the excitation frequency changes from 0.43 to 1.22, the bifurcation diagram contains four periodic and three chaotic motion regions. Periodic motion windows are found in the three chaotic motion regions. It is confirmed that the frequency and amplitude of the external excitation and the fractional order of damping have a greater impact on system dynamics. Thus, attention shall be paid to the design and analysis of system dynamics.


Introduction
Fractional calculus was proposed almost at the same time as integer calculus. Its development was relatively slow for nearly 300 years due to the lack of a clear physical meaning of fractional calculus, and it was mainly studied in the field of pure mathematics. In 1974, Oldham and Spanier [1] co-authored the first monograph on fractional calculus titled 'Applied Fractional Calculus', which paved its way towards the application of fractional calculus. In the 1970s, Mandelbrot [2] pointed out that many fractional dimensions exist in nature. Since then, the theory of fractional calculus developed rapidly, and numerous applications of fractional calculus emerged [3]- [9]. The application of fractional calculus in mechanical engineering gradually increased in recent years mainly because many physical objects had fractional characteristics, such as viscoelasticity, damping, mechanical friction and impact [10]- [17].
Many popular fractional systems, such as the Lorenz system [5], the Rossler system [13], Chua's circuit [6] and the Duffing system [9], were used for research. The Duffing system was applied in many fields, including fluid-induced vibration, large-amplitude oscillation of centrifugal speed control systems and mathematical modeling [9], [18], [19]. The Duffing system that considers fractional differential damping is equivalent to introducing polymer damping into the classic Duffing system. Compared with the traditional integer-order damping model, the fractional-order damping model can more accurately describe the dynamic behavior of the system [20]- [22]. In Reference [23], the nonlinear dynamic behavior of the fractional-order Duffing system was studied using the Adams-Bashforth-Moulton predictor-corrector method. The Adomian decomposition method was utilized to find the numerical solution to the combined fractional differential equation of the Duffing equation and the Van der Pol equation in Reference [24]. The Melnikov method transforms a fractional-order system into an equivalent integer-order system, which can effectively detect chaos in the fractional-order Duffing system [25]. Based on the definition formula with an amount of calculations, reference [26] studies the bifurcation and vibration resonance of fractional double-damping Duffing time delay system driven by external excitation signal with two wildly different frequencies. In Ref. [27], phase diagrams, bifurcation diagrams and Lyapunov exponents are used to determine the presence of chaos over a wide range of fractional orders, based on the Caputo fractional difference form of Duffing maps.
In this study, the improved short memory principle method is used to study the dynamic characteristics of the fractional Duffing system. Firstly, the Grünwald-Letnikov definition of fractional calculus is adopted. Secondly, the classical short memory principle is improved, and the numerical solution of the system is obtained. Lastly, the influence of excitation frequency on the dynamic performance of the fractional Duffing system is studied using nonlinear dynamic analysis methods, such as phase portrait, Poincare map and bifurcation diagram. Moreover, the dynamic behavior of the fractional Duffing system when the external excitation amplitude and fractional order change are investigated.

Fractional differential definition
The three most frequently used definitions for fractional derivatives are Grünwald-Letnikov, Riemman-Liouville and Caputo definitions. The Grünwald-Letnikov definition is used in this work, as it is expressed as: where is the fractional order, ℎ is the time step and [•] means rounding. If the time step ℎ is small enough, the above-mentioned formula can be written as: where is the coefficient of binomial 1 − : Eq. (3) is changed to the recursive form: Its first item = 1.

Improved short memory principle method
The classical short memory principle truncates the memory time. In Eq. (2), the upper limit of summation becomes ℎ ⁄ , where is the fixed memory time length, and the part exceeding is discarded. The discarded part often brings errors that cannot be ignored. In particular, when fractional order approaches zero, the binomial coefficient changes slowly.
The improved short memory principle method changes the truncation of the memory time in the classic short memory principle to the number of terms of the binomial coefficient. Then, the finite binomial coefficient is repeatedly applied to gradually enlarging time steps for error reduction. The expression of a time step being enlarged for the first time is: The definition of Grünwald-Letnikov is used in this study, and the formulas are no longer marked. In Eq. (5), is the number of binomial truncated terms and is the magnification of the step length. Similarly, the expression of the second enlargement of step size is: In the same manner, other expressions of step size can be derived.

Basis for step enlargement
Step enlargement is based on the initial time step. After the step size is enlarged, although the interval for taking the function value increases, the binomial coefficient of the function value also increases, and the magnification is almost the same as the step size magnification. For example, Eq. (5) can be rewritten as: The second term at the right end of the formula above is the result of the step size being enlarged for the first time. The binomial coefficient before the function value − ℎ is ⁄ . If the step size is not enlarged, the binomial coefficient in front of the function value shall be . Table 1 shows the ratio of the corresponding binomial coefficients when = 100, = 0.5 and are 5, 10 and 20, respectively. is the ratio of binomial coefficient ⁄ after the first step is enlarged to binomial coefficient . is the ratio of binomial coefficient ⁄ after the second step is enlarged to binomial coefficient . The ratio is almost equal to step magnification and , and the ratio increases with the increase in step magnification. After this increase, the first binomial coefficient increases; in particular, when = 20, it is larger than the number of initial steps in the enlarged step. A large magnification causes a large error. In addition, whether the adopted function value is representative in the interval, that is, whether the adopted function value is close to the average value of the function value in the interval, which is also a question to be considered. The error in the long-term calculation is limited because expanding the step size to derive the function value is a dynamic ergodic process. In summary, it is required to try to intercept as many binomial coefficients as possible if the calculation conditions permit, and a small step size magnification shall be selected to obtain a highly accurate numerical solution.

Fractional duffing system
The form of the classic integer-order Duffing system is: where is the first-order differential operator, is the second-order differential operator, is the mass, is the damping coefficient, is the linear stiffness coefficient, is the nonlinear stiffness coefficient and and are the amplitude and frequency of external excitation, respectively. If the damping term is changed to the fractional order, it becomes a fractional-order Duffing system: where is the fractional order of the damping of the fractional Duffing system. In this study, the improved short memory principle is used to solve the fractional Duffing equation.
Fractional order = 1.2 and excitation amplitude = 0.6 are adopted to analyze the influence of excitation frequency on the system. The calculation parameters are as follows: the basic step size is ℎ = 0.001, the number of truncated items is = 10,000, and the step length magnification is = 5. Excitation frequency is changed from 0.43 to 1.23 in the experiment to analyze the influence of excitation frequency on the fractional Duffing system, and the change step is 0.001. A total of 250 cycles is calculated for each frequency point, and the first 150 cycles are discarded. The bifurcation diagram of the system obtained with as the control parameter is shown in Fig. 1. respectively. The phase portrait is a universal limit cycle; correspondingly, the Poincaré map exhibits an isolated point. Notably, 0.488 < ≤ 0.518 is the first chaotic motion region of the system. The phase portrait and Poincaré map when = 0.509 are shown in Figs. 2(b1) and 2(b2), respectively. On the Poincaré map, many discrete points are concentrated in two areas. With the increase in excitation frequency , the system enters the period-3 motion region, and its frequency range is 0.519 ≤ ≤ 0.577. The phase portrait and Poincaré map when = 0.55 are shown in Figs. 2(c1) and 2(c2), respectively. The Poincaré map appears as three isolated points, and 0.578 ≤ ≤ 0.657 is the second chaotic region of the system. The phase portrait when = 0.6 is shown as Fig. 2(d1), and scattered points can be seen on the corresponding Poincaré map in Fig. 2(d2). When ≥ 0.658, the system enters a stable period-5 motion region, and its interval is 0.658 ≤ ≤ 0.848. The phase portrait and Poincaré map when = 0.8 are shown in Figs. 2(e1) and 2(e2), respectively. The Poincaré map has five scattered points, and 0.849 ≤ ≤ 0.995 is the third chaotic motion region. Fig. 2(f1) shows the phase portrait when = 0.88, and the corresponding Poincaré map is given in Fig. 2(f2). As further increases, the system enters periodic motion from chaotic motion. Fig. 1 indicates that when ≥ 1.17, the system moves stably in period 1. A period-doubling bifurcation transition region exists between the periodic motion region and the third chaotic motion region. The phase portrait and Poincaré map of the transition region when = 1.02 are shown in Figs. 2(g1) and 2(g2), respectively.
Notably, periodic motion windows are present in the three chaotic motion regions. The bifurcation diagram of the periodic motion window in the first chaotic motion region is shown in Fig. 3(a). The chaotic motion transits from period-doubling motion to period-2 motion. The phase portrait and Poincaré map when = 0.506 in the period window are shown in Figs. 3(b) and 3(c), respectively. The bifurcation diagram of the periodic motion window in the second chaotic motion region is shown in Fig. 4(a). In the period window, = 0.635 is period-4 motion, and its phase portrait and Poincaré map are shown in Figs Notably, an interesting change occurs in the periodic window of the second chaotic motion region. Although the system moves in period 4 of the window, the bifurcation diagram has eight dotted lines instead of four solid lines. Although the system exhibits period-4 motion in the period window, the bifurcation diagram has eight dotted lines instead of four solid lines. This observation shows that the motion of the system jumps with the change in frequency in the period window, but it is still period-4 motion. For example, the phase portraits when = 0.6352 and = 0.6353 are shown in Figs. 6(a) and 6(b), respectively. The phase portraits of two adjacent frequency points are rotated by exactly 180 degrees. When the frequency is further refined, this jumping property remains. Figs. 6(c) and 6(d) present phase portraits with of 0.63280 and 0.63281, respectively.

Fractional order duffing system with fractional order or excitation amplitude fluctuations
This paper also examines the dynamic characteristics of the system when the fractional order changes. The excitation amplitude is = 1, the excitation frequency is = 1, the basic step length is ℎ = 0.001, the basic truncation term is = 10,000, and the step length magnification is = 5. its memory weakens as the fractional order increases. Thus, the number of truncated items is set to decrease by 1000 items with each increase in by 0.2. The bifurcation diagram is shown in Fig. 7(a). In the process of changing fractional order from 0.2 to 2, three periodic and three chaotic motion regions appear. The transition is from the first chaotic zone to a stable period-1 motion region after period-doubling bifurcation. As continues to increase, it enters the second chaotic motion region through period-doubling bifurcation. Fig. 7(b) presents the bifurcation diagram of the system when excitation amplitude fluctuates from 0.1 to 2.1. The fractional order is = 0.8, the excitation frequency is = 1, the basic step length is ℎ = 0.001, the basic truncation item number is = 10,000, and the step length magnification is = 5. The figure shows that three periodic and two chaotic motion regions appear during the fluctuation of excitation amplitude. The first is a stable period-1 motion; the first chaotic motion region is entered through period-doubling bifurcation, followed by entrance to the period-3 motion region. A part of the period-6 movement is contained in the period-3 motion region, and the second chaotic motion region is entered from period-3 motion. As the excitation amplitude continues to increase, the system transitions from the second chaotic region to the period-1 motion region through the period-doubling bifurcation. a) b) Fig. 7. a) Bifurcation diagram of versus , 0.2 < < 2 and b) bifurcation diagram of versus , 0.1 < < 2.1

Conclusions
An improved short memory principle method is introduced to solve the fractional Duffing system numerically. And based on high computational efficiency, some previously unobserved phenomena are obtained. The influence of excitation frequency on the dynamic performance of the system is studied using nonlinear dynamic analysis methods, such as phase portrait, Poincaré map and bifurcation diagram. The bifurcation diagram has four periodic and three chaotic motion regions. The third chaotic motion region exhibits an obvious transition from period-doubling motion to stable period-1 motion. In addition, periodic motion windows are found in the three chaotic motion regions. The periodic motion of the periodic window in the second chaotic motion region has a jumping property. The influence of excitation amplitude and fractional order on the dynamic system performance is also studied. The results show that a change in the system parameters of the fractional Duffing equation affects the dynamic system performance. Thus, such a change shall be considered during the dynamic system design and analysis.