Period-Doubling Bifurcation of Stochastic Fractional-Order Duffing System via Chebyshev Polynomial Approximation

Fractional-order calculus is more competent than integer-order one when modeling systems with properties of nonlocality and memory effect. And many real world problems related to uncertainties can be modeled with stochastic fractional-order systems with random parameters. Therefore, it is necessary to analyze the dynamical behaviors in those systems concerning both memory and uncertainties.The period-doubling bifurcation of stochastic fractional-order Duffing (SFOD for short) system with a bounded random parameter subject to harmonic excitation is studied in this paper. Firstly, Chebyshev polynomial approximation in conjunction with the predictor-corrector approach is used to numerically solve the SFOD system that can be reduced to the equivalent deterministic system. Then, the global and local analysis of period-doubling bifurcation are presented, respectively. It is shown that both the fractional-order and the intensity of the random parameter can be taken as bifurcation parameters, which are peculiar to the stochastic fractional-order system, comparing with the stochastic integer-order system or the deterministic fractional-order system.Moreover, the Chebyshev polynomial approximation is proved to be an effective approach for studying the period-doubling bifurcation of the SFOD system.


Introduction
Fractional calculus [1] had been proposed more than 300 years along with classical calculus, but really came to life over the last few decades.A particular feature is that some fractional models [2,3], which can describe the real object more accurately than classical "integer" counterpart, have been developed by scientists for its nonlocal advantage of fractional derivative and are widely used in many fields.For example, the Scott-Blair model and the Kelvin-Voigt model provide a powerful instrument for describing the properties of viscoelastic materials such as gels and rubbers [4,5].Fractional Fokker-Planck equations can be found effective for describing the inherent abnormal-exponential or heavy tail decay processes which can occur in biomechanics, chemistry, and acoustics [6,7].The GARCH model proposed by Bollerslev and the EGARCH model presented by Nelson have also been demonstrated to be able to model long memory in stock market volatility [8].But in fact, due to the random errors incurred by measuring, manufacturing, and assembling, the parameters for the fractional systems can be modeled as random parameters with certain distribution, which are generally bounded in engineering.Consequently, the fractional-order system with bounded random parameters is more appropriate to model the real phenomenon.
There have been several mathematical methods used to obtain the response of the system with random parameters numerically, for example, Monte Carlo method [9], stochastic perturbation method [10,11], and orthogonal polynomial approximation method proposed by Spanos and Ghanem [12] and Jensen and Iwan [13] and developed by Li et al. [14,15] who presented a new scheme for the estimation of the mean values of parameters in multiple degree of freedom structural systems and a dynamic condensation algorithm for the orderexpanded equation.The first method is universal but requires a very long simulation time.The second method based on the assumption of only small stochastic perturbation is involved with least computation.The last method is more effective than the first two ones for its advantages of the computation speed and the random perturbation assumption and has 2 Shock and Vibration been used to analyze the stochastic systems [16,17].For example, the Chebyshev polynomial approximation [18][19][20], one form of orthogonal polynomial approximation, is applied to explore dynamical behaviors in the integer-order systems with bounded random parameters, for example, perioddoubling bifurcation, symmetry-breaking bifurcation, Hopf bifurcation, and chaos.Although the approach for dealing with stochastic integer-order systems is well developed, there is no general method to solve stochastic fractional-order systems effectively.In this paper the Chebyshev polynomial approximation is extended to the stochastic fractional system.
The stochastic fractional-order Duffing (SFOD) system which means the fractional-order Duffing (FOD for short) system with a bounded random parameter is studied.The FOD system, the foundation of complex fractional systems, is a generalization of Duffing system, based on multiple ways, for example, fractionalizing the Duffing system in a state space [21][22][23] and fractionalizing the Duffing system by considering a fractional damping term [24,25].The SFOD system we consider is inspired from fractionalizing the Duffing system in a state space with negative linear stiffness, damping, and periodic excitation, which can be represented as follows: where    is the Caputo derivative, ,  are state variables,  is a random parameter, and  cos() is a harmonic excitation.Many researchers have been devoted to study the dynamical behaviors of the FOD system and conclude that the fractional version has different behaviors compared to the integer-order system.In [21,22], using numerical simulations it has been shown that the FOD system and the modified FOD system (nonautonomous systems) with order less than 2 can still behave in a chaotic manner, whereas for the integer nonautonomous system the minimum order must be two.In [23] the parametric range for undamped oscillations in FOD oscillator was found using the stability theorem for fractional system where the system order is 1.7, and numerical simulation results justify the analysis.While all the above is restricted to the deterministic FOD system without a random parameter, the one with a random parameter has not been fully exploited.
Inspired by the discussion above, we focus on analyzing the period-doubling bifurcation of the SFOD system using Chebyshev polynomial approximation.The paper is organized as follows.In Section 2 basic preliminaries including the definition of fractional derivative, numerical algorithms for fractional differential equations and Chebyshev polynomial are introduced.Section 3 is devoted to derive the equivalent deterministic system of the SFOD system using Chebyshev polynomial approximation.In Section 4 the period-doubling bifurcation in SFOD system where both the fractional-order and the intensity of the random parameter can be taken as bifurcation parameters is explored, and the Chebyshev polynomial approximation is proved to be effective.The last section is devoted to a brief conclusion.

Basic Preliminaries
2.1.Fractional Derivative.As noted in [26], there are several fractional derivative definitions, that is, Caputo derivative, Riemann-Liouville derivative, and Grünwald-Letnikov derivative.However, the Caputo derivative, making the initial conditions of fractional differential equations with the same form as for integer-order differential equations possible, is the most popular fractional derivative in practice.This is actually the reason for us to choose the Caputo derivative defined as where  fl ⌈⌉ is just the value of  rounded up to the nearest integer and   is the Riemann-Liouville integral operator defined by where  > 0 and  represents Gamma function.
After the description of the definition of fractional derivative, we now introduce numerical algorithms for fractional differential equations.Two methods are commonly used for numerical simulation of fractional-order systems.One is the time-domain method known as the predictor-corrector approach which is complicated but accurate [27]; the other is the frequency-domain approximation which provides a simple procedure for simulating the fractional-order systems numerically but has some limitations [28].Considering the accuracy of calculation, we employ the predictor-corrector approach.
Then we briefly recall the idea behind the predictorcorrector approach for the following fractional differential equation equipped with suitable initial conditions which is equivalent to the Volterra integral equation For simplicity, it is assumed that we are working on a uniform grid   = ℎ, ( = 0, 1, . . ., ) and ℎ = /.The integral in (5) can be approximate to where g+1 () is the piecewise linear interpolant for () with nodes and knots chosen at the   ( = 0, 1, . . .,  + 1), noting that (, ()) in ( 5) is written as () in ( 6).Applying the trapezoidal quadrature formula and the product rectangle rule to the right hand side of ( 6) respectively, we can have where Then the predictor and corrector formulae solving (4) are given by where   ℎ ( +1 ) is the predicted value and  ℎ ( +1 ) is the corrected value.In [29], it is shown that the error of this mathematical approach is estimated as where 2.2.Chebyshev Polynomial.The choice of the orthogonal polynomial basis functions depends on the probability density function (PDF) of the random variable.In fact, random parameters existing in engineering structures are generally bounded, so we choose the random variable with an Wigner semicircle distribution, which is shown in Figure 1 and can be expressed as follows [17]: Moreover, the Wigner semicircle distribution is orthogonal to Chebyshev polynomials of the second kind presented as [17]   () = Then, we have The recurrent formula for the second kind of Chebyshev polynomials is expressed in the form The orthogonality of the second kind of Chebyshev polynomials is represented as Since the weighting function in ( 16) is just the same as the Wigner semicircle PDF () in ( 12), the left-hand side of ( 16) can be viewed as the expectation of the product   ()  ().
Based on the orthogonality of Chebyshev polynomials, any measurable function () ∈  2 can be represented in the Shock and Vibration following form: where Similarly, the results for measurable functions of several mutually independent random variables can also be obtained [17,30].

Chebyshev Polynomial Approximation for SFOD System
In this section we are devoted to derive the equivalent deterministic system of the SFOD system using Chebyshev polynomial approximation and further obtain the system responses by the predictor-corrector approach as above proposed.Consider the SFOD system (1) subject to harmonic excitation, in which we suppose that  is a random parameter represented by where  is the mean value of ,  is a bounded random variable defined on [−1, 1] with the Wigner semicircle PDF, and  is regarded as the intensity of the random variable .Then the SFOD system with a random parameter can be written as The response of system (20) is stochastic due to the random variable , and it can be viewed as the function with respect to time and random variable; namely, Due to the fact that the stochastic function space constituted with the responses of nonlinear dynamic systems with random parameters is a Hilbert space with respect to convergence in the mean square [30], response (21) can be expressed approximately by the following series: where   () represents the th Chebyshev polynomial and  refers to the largest order of the polynomials.
Substituting (22) into (20), we can get In fact, the product of two Chebyshev polynomials can be reduced into a linear combination of individual Chebyshev polynomial according to the recurrent formula (15), meaning that the term (∑  =0   ()  ()) 3 in ( 23) can be represented by the following form: where   () can be obtained by Maple (V16) and the term  ∑  =0   ()  () in ( 23) can also be expressed by the following form: where  −1 and  +1 are supposed to be zero in the approximation.
Then substituting ( 24) and ( 25) into ( 23), we have Multiplying both sides of ( 26) by   () for ( = 0, 1, . . ., ) and then taking expectation with respect to the random variable , we finally obtain the equivalent deterministic system of the SFOD system.Note that, if  → ∞, ∑  =0   ()  () and ∑  =0   ()  () are strictly equivalent to the responses (, ) and (, ), respectively.However, in fact we can only suppose that  is finite; that is, response ( 22) is just an approximate value.According to references [19,20], we take  = 4 and further the equivalent deterministic system can be expressed approximately as The ensemble mean response of the SFOD system (EMR for short) can be expressed as Particularly if  ≡ 0 or  = 0, the original stochastic system (20) is reduced to the mean-parameter system of which the determinate responses () and () can be approximated as known as the sample response of the mean-parameter system (SRM for short).
Applying the predictor-corrector approach to system (27), we can obtain the responses   () and   (), ( = 0, 1, 2, 3, 4), and further EMR (namely, [(, )] and [(, )]) and SRM (namely, (, 0) and (, 0)) can be worked out, respectively.Moreover, for illustrating the effectiveness of the Chebyshev polynomial approximation in the next section, it is necessary to obtain the determinate response of the mean-parameter system (29) (DRM for short) and the response of the original stochastic system (20) using Monte Carlo simulation method (MCS for short).

Period-Doubling Bifurcation in SFOD System
In this section the dynamical behaviors of the equivalent deterministic system ( 27) is analyzed by bifurcation diagram, phase portrait, and time history diagram.The main idea of studying dynamical behaviors is shown below.First, we choose the fractional-order  and the random intensity  as bifurcation parameters, which are peculiar to the stochastic fractional-order system, comparing with the stochastic integer-order system or the deterministic fractional-order system.Then, bifurcation diagrams of the fractional-order  and the random intensity , illustrating the variation of dynamical behaviors from the global perspective, are, respectively, shown.Based on the global analysis of dynamical behaviors in the equivalent deterministic system (27), we finally use phase portraits and time history diagrams for some given values ( or ) to verify whether period-doubling bifurcation happens in the equivalent deterministic system (27) as the fractional-order  or the random intensity  varies.

𝑞 as Bifurcation Parameter.
In this case, the fractionalorder , which is peculiar to the fractional-order system comparing with the integer-order system, is taken as bifurcation parameter.In order to systematically describe the dynamical behaviors resulting from the variation of fractional-order , bifurcation diagram of the fractional-order  is shown in Figure 2 with the fractional-order  on the horizontal axis and the ensemble mean response (EMR)  0 () (28) on the vertical one.As to the numerical simulation of bifurcation diagram, we take 50 forcing periods  = 2/ as transient value and 50 periods are represented.Also, Δ = 0.01 and  = 0.5,  = 8.   = 1.0,  = 0.001 are taken as the time step size and the system parameters, respectively.In addition, the initial conditions related to the original stochastic system (20) and the equivalent deterministic system are choose, respectively, as From the bifurcation diagram depicted in Figure 2, we can see that the period of the ensemble mean response (EMR) varies from 1 to 2 and then to 4 and finally period-doubling bifurcations lead to chaos as the fractionalorder  varies.Then, we choose some given values of the fractional-order , as 0.984, 0.988 and 0.9895, to validate that the ensemble mean response (EMR) experiences periodic responses with 1, 2, and 4.
When  = 0.984, from the phase portrait and time history diagram shown in Figures 3(c) and 3(d), it can be observed that the system has a steady-state periodic EMR with period 1, which is consistent with the bifurcation diagram in Figure 2.Moreover, the consistencies over DRM-SRM (illustrated in Figures 3(a) and 3(b)) and EMR-MCS (depicted in Figures 3(c) and 3(d)) reveal the effectiveness of Chebyshev polynomial approximation.
Changing  from 0.984 to 0.988, we can see that there is a steady-state periodic EMR with period 2 indicated in Figures 4(c) and 4(d), which also is observed from the bifurcation diagram in Figure 2. Also, in the Figures 4(a) and 4(b) it can be clearly seen that DRM and SRM are close to each other and so are EMR and MCS (see Figures 4(c) and 4(d)), which confirms the validity of the Chebyshev polynomial approximation.
Further increasing  to 0.9895, the system has a steadystate periodic EMR with period 4 (see Figures 5(c As to the analysis above, we can conclude that Chebyshev polynomial approximation is an effective method for studying the period-doubling bifurcation resulting from the variation of the fractional-order , and period-doubling bifurcation takes place as the fractional-order  varies from 0.984 to 0.9895, followed by a set of periodic responses with period 1, 2, and 4.That is to say, the fractionalorder  can lead to period-doubling bifurcation in stochastic fractional-order system, which is a new phenomenon comparing with the fact that we generally take the amplitude and the frequency of harmonic excitation as the period-doubling bifurcation parameters in the classical integer-order system with a random parameter [20,31].

𝜎 as Bifurcation Parameter.
In this case, the random intensity , which is peculiar to the stochastic system with random parameters comparing with the deterministic system, is taken as bifurcation parameter.In order to present the dynamical behaviors resulting from the variation of the random intensity  in the equivalent deterministic system (27) from the global perspective, the bifurcation diagram of the random intensity  is shown in Figure 6 with the random intensity  on the horizontal axis and the ensemble mean response (EMR)  0 () (28) on the vertical one.As to the numerical simulation of the bifurcation diagram, we also take 50 forcing periods as transient value and 50 periods are represented.In addition,  = 0.5,  = 8.3,  = 1.0, and  = 0.984 are taken as the system parameters, and the time step size and the initial conditions are the same as the above.
From the bifurcation diagram depicted in Figure 6, it can be observed that the period of the ensemble mean response (EMR) varies from 1 to 2 and finally period-doubling bifurcations lead to chaos when the random intensity  increases.Then, we choose some given values of the random intensity , as 0.01 and 0.06, to validate that the ensemble mean response (EMR) experiences periodic responses with 1 and 2.
For  = 0.01, in Figure 7 it can be observed that the system has a steady-state periodic EMR with period 1, which is consistent with the bifurcation diagram in Figure 6.Increasing  to 0.06, we can see that there is a steadystate periodic EMR with period 2 illustrated in Figure 8. Furthermore, we can demonstrate the validity of the Chebyshev polynomial approximation based on the fact that EMR and MCS are nearly close to each other as depicted in Figures 8(a Based on the discussion above, we can conclude that Chebyshev polynomial approximation is an effective method for studying the period-doubling bifurcation resulting from the variation of the random intensity , and period-doubling bifurcation happens when we range  from 0.01 to 0.06, followed by periodic responses with period 1 and 2.Therefore, the random intensity  can be taken as the bifurcation parameter, which is peculiar to the stochastic system with random parameters comparing to the deterministic system.

Conclusion
The period-doubling bifurcation of the SFOD system subject to harmonic excitation is explored.Chebyshev polynomial approximation is used to derive the equivalent deterministic system of the SFOD system and further obtain the system responses which are useful for studying the period-doubling bifurcation of the SFOD system by the predictor-corrector approach.The result reveals that both the fractional-order and the intensity of the random parameter can be taken as bifurcation parameters, which is peculiar to the stochastic fractional-order system, comparing with the stochastic integer-order system or the deterministic fractional-order system.Moreover, the Chebyshev polynomial approximation is proved to be an effective approach for studying the perioddoubling bifurcation of the SFOD system.Besides, maybe the proposed Chebyshev polynomial approximation can be further used to study other nonlinear dynamical behaviors in the SFOD system, for example, Hopf bifurcation, symmetrybreaking bifurcation, and chaos.+ 3 3 () 2  4 () +  0 () 3 +  2 () 3   +  4 ()