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

Authors: Youming Lei; Yanyan Wang
Journal: Shock and Vibration
 (2017)

---

## Abstract

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.


---

## Body

## 1. 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 order-expanded 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 been used to analyze the stochastic systems [16, 17]. For example, the Chebyshev polynomial approximation [18–20], one form of orthogonal polynomial approximation, is applied to explore dynamical behaviors in the integer-order systems with bounded random parameters, for example, period-doubling 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–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:(1)Dcqx=y,Dcqy=x-x3-αy+fcos⁡wt,where Dcq is the Caputo derivative, x, y are state variables, α is a random parameter, and fcos⁡wt 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 Section2 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.


## 2. 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(2)Dcqxt=Jn-qdndtnxt,where n≔q is just the value of q rounded up to the nearest integer and Jβ is the Riemann-Liouville integral operator defined by(3)Jβxt=1Γβ∫0tt-τβ-1xτdτ,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 predictor-corrector approach for the following fractional differential equation equipped with suitable initial conditions(4)Dcqxt=ft,xt,0≤t≤T,xk0=x0k,k=0,1,…,q-1,which is equivalent to the Volterra integral equation(5)xt=∑k=0q-1x0ktkk!+1Γq∫0tt-τq-1fτ,xτdτ.For simplicity, it is assumed that we are working on a uniform grid tn=nh, (n=0,1,…,N) and h=T/N. The integral in (5) can be approximate to (6)∫0tn+1tn+1-τq-1gτdτ≈∫0tn+1tn+1-τq-1g~n+1τdτ,where g~n+1τ is the piecewise linear interpolant for gτ with nodes and knots chosen at the tjj=0,1,…,n+1, noting that fτ,xτ in (5) is written as gτ in (6). Applying the trapezoidal quadrature formula and the product rectangle rule to the right hand side of (6) respectively, we can have(7)∫0tn+1tn+1-τq-1g~n+1τdτ=hqqq+1∑j=0n+1aj,n+1gtj,∫0tn+1tn+1-τq-1g~n+1τdτ=∑j=0nbj,n+1gtj,where (8)aj,n+1=nq+1-n-qn+1q,j=0,n-j+2q+1+n-jq+1-2n-j+1q+1,1≤j≤n,1,j=n+1,bj,n+1=hqqn-j+1q-n-jq,0≤j≤n.Then the predictor and corrector formulae solving (4) are given by(9)xhPtn+1=∑k=0q-1x0ktn+1kk!+1Γq∑j=0nbj,n+1ftj,xhtj,xhtn+1=∑k=0q-1x0ktn+1kk!+hqΓq+2ftn+1,xhPtn+1+hqΓq+2∑j=0naj,n+1ftj,xhtj,where xhPtn+1 is the predicted value and xhtn+1 is the corrected value. In [29], it is shown that the error of this mathematical approach is estimated as(10)ε=maxj=0,1,…,N⁡xtj-xhtj=Ohp,where(11)p=min1+q,2.


### 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 Figure1 and can be expressed as follows [17]: (12)pu=2π1-u2u≤1,0u>1.Moreover, the Wigner semicircle distribution is orthogonal to Chebyshev polynomials of the second kind presented as [17](13)Hnu=∑k=0n/2-1kn-k!k!n-2k!2un-2k.Then, we have(14)H0u=1,H1u=2u,H2u=4u2-1,H3u=8u3-4u,H4u=16u4-12u2+1,⋮The recurrent formula for the second kind of Chebyshev polynomials is expressed in the form(15)uHnu=12Hn-1u+Hn+1u.The orthogonality of the second kind of Chebyshev polynomials is represented as(16)∫-112π1-u2HiuHjudu=1i=j,0i≠j.Since the weighting function in (16) is just the same as the Wigner semicircle PDF pu in (12), the left-hand side of (16) can be viewed as the expectation of the product HiuHju. Based on the orthogonality of Chebyshev polynomials, any measurable function fu∈L2 can be represented in the following form:(17)fu=∑i=0∞ciHiu,where(18)ci=∫-11pufuHiudu.Similarly, the results for measurable functions of several mutually independent random variables can also be obtained [17, 30].Figure 1
The Wigner semicircle PDFp(u) curve for random variable u.


## 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(2)Dcqxt=Jn-qdndtnxt,where n≔q is just the value of q rounded up to the nearest integer and Jβ is the Riemann-Liouville integral operator defined by(3)Jβxt=1Γβ∫0tt-τβ-1xτdτ,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 predictor-corrector approach for the following fractional differential equation equipped with suitable initial conditions(4)Dcqxt=ft,xt,0≤t≤T,xk0=x0k,k=0,1,…,q-1,which is equivalent to the Volterra integral equation(5)xt=∑k=0q-1x0ktkk!+1Γq∫0tt-τq-1fτ,xτdτ.For simplicity, it is assumed that we are working on a uniform grid tn=nh, (n=0,1,…,N) and h=T/N. The integral in (5) can be approximate to (6)∫0tn+1tn+1-τq-1gτdτ≈∫0tn+1tn+1-τq-1g~n+1τdτ,where g~n+1τ is the piecewise linear interpolant for gτ with nodes and knots chosen at the tjj=0,1,…,n+1, noting that fτ,xτ in (5) is written as gτ in (6). Applying the trapezoidal quadrature formula and the product rectangle rule to the right hand side of (6) respectively, we can have(7)∫0tn+1tn+1-τq-1g~n+1τdτ=hqqq+1∑j=0n+1aj,n+1gtj,∫0tn+1tn+1-τq-1g~n+1τdτ=∑j=0nbj,n+1gtj,where (8)aj,n+1=nq+1-n-qn+1q,j=0,n-j+2q+1+n-jq+1-2n-j+1q+1,1≤j≤n,1,j=n+1,bj,n+1=hqqn-j+1q-n-jq,0≤j≤n.Then the predictor and corrector formulae solving (4) are given by(9)xhPtn+1=∑k=0q-1x0ktn+1kk!+1Γq∑j=0nbj,n+1ftj,xhtj,xhtn+1=∑k=0q-1x0ktn+1kk!+hqΓq+2ftn+1,xhPtn+1+hqΓq+2∑j=0naj,n+1ftj,xhtj,where xhPtn+1 is the predicted value and xhtn+1 is the corrected value. In [29], it is shown that the error of this mathematical approach is estimated as(10)ε=maxj=0,1,…,N⁡xtj-xhtj=Ohp,where(11)p=min1+q,2.


## 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 Figure1 and can be expressed as follows [17]: (12)pu=2π1-u2u≤1,0u>1.Moreover, the Wigner semicircle distribution is orthogonal to Chebyshev polynomials of the second kind presented as [17](13)Hnu=∑k=0n/2-1kn-k!k!n-2k!2un-2k.Then, we have(14)H0u=1,H1u=2u,H2u=4u2-1,H3u=8u3-4u,H4u=16u4-12u2+1,⋮The recurrent formula for the second kind of Chebyshev polynomials is expressed in the form(15)uHnu=12Hn-1u+Hn+1u.The orthogonality of the second kind of Chebyshev polynomials is represented as(16)∫-112π1-u2HiuHjudu=1i=j,0i≠j.Since the weighting function in (16) is just the same as the Wigner semicircle PDF pu in (12), the left-hand side of (16) can be viewed as the expectation of the product HiuHju. Based on the orthogonality of Chebyshev polynomials, any measurable function fu∈L2 can be represented in the following form:(17)fu=∑i=0∞ciHiu,where(18)ci=∫-11pufuHiudu.Similarly, the results for measurable functions of several mutually independent random variables can also be obtained [17, 30].Figure 1
The Wigner semicircle PDFp(u) curve for random variable u.


## 3. 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(19)α=α-+σu,where α¯ is the mean value of α, u is a bounded random variable defined on [-1,1] with the Wigner semicircle PDF, and σ is regarded as the intensity of the random variable u. Then the SFOD system with a random parameter can be written as(20)Dcqx=y,Dcqy=x-x3-α-+σuy+fcos⁡wt.The response of system (20) is stochastic due to the random variable u, and it can be viewed as the function with respect to time and random variable; namely, (21)x=xt,u,y=yt,u.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:(22)xt,u=∑i=0NxitHiu,yt,u=∑i=0NyitHiu,where Hiu represents the ith Chebyshev polynomial and N refers to the largest order of the polynomials.Substituting (22) into (20), we can get(23)Dcq∑i=0NxitHiu=∑i=0NyitHiu,Dcq∑i=0NyitHiu=∑i=0NxitHiu-∑i=0NxitHiu3-α-∑i=0NyitHiu-σu∑i=0NyitHiu+fcos⁡wt.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 ∑i=0NxitHiu3 in (23) can be represented by the following form:(24)∑i=0NxitHiu3=∑i=03NSitHiu,where Sit can be obtained by Maple (V16) and the term σu∑i=0NyitHiu in (23) can also be expressed by the following form:(25)σu∑i=0NyitHiu=σ2∑i=0NyitHi-1u+Hi+1u=σ2∑i=0Nyi-1t+yi+1tHiu+yNtHN+1u,where y-1 and yN+1 are supposed to be zero in the approximation.Then substituting (24) and (25) into (23), we have (26)Dcq∑i=0NxitHiu=∑i=0NyitHiu,Dcq∑i=0NyitHiu=∑i=0NxitHiu-∑i=03NSitHiu-α-∑i=0NyitHiu-σ2∑i=0Nyi-1t+yi+1tHiu+fcos⁡wt.Multiplying both sides of (26) by Hiu for i=0,1,…,N and then taking expectation with respect to the random variable u, we finally obtain the equivalent deterministic system of the SFOD system. Note that, if N→∞, ∑i=0NxitHiu and ∑i=0NyitHiu are strictly equivalent to the responses xt,u and yt,u, respectively. However, in fact we can only suppose that N is finite; that is, response (22) is just an approximate value. According to references [19, 20], we take N=4 and further the equivalent deterministic system can be expressed approximately as(27)Dcqx0t=y0t,Dcqy0t=x0t-S0t-α-y0t-σ2y1t+fcos⁡wt,Dcqx1t=y1t,Dcqy1t=x1t-S1t-α-y1t-σ2y0t+y2t,Dcqx2t=y2t,Dcqy2t=x2t-S2t-α-y2t-σ2y1t+y3t,Dcqx3t=y3t,Dcqy3t=x3t-S3t-α-y3t-σ2y2t+y4t,Dcqx4t=y4t,Dcqy4t=x4t-S4t-α-y4t-σ2y4t.The ensemble mean response of the SFOD system (EMR for short) can be expressed as (28)Ext,u=∑i=04xitEHiu=x0t,Eyt,u=∑i=04yitEHiu=y0t.Particularly if u≡0 or σ=0, the original stochastic system (20) is reduced to the mean-parameter system(29)Dcqx=y,Dcqy=x-x3-α-y+fcos⁡wt,of which the determinate responses xt and yt can be approximated as (30)xt≈xt,0=∑i=04xitHi0=x0t-x2t+x4t,yt≈yt,0=∑i=04yitHi0=y0t-y2t+y4t,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 xit and yit, i=0,1,2,3,4, and further EMR namely,Ext,uandEyt,u and SRM namely,xt,0andyt,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).


## 4. 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 q 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 q 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 (q or σ) to verify whether period-doubling bifurcation happens in the equivalent deterministic system (27) as the fractional-order q or the random intensity σ varies.


### 4.1.q as Bifurcation Parameter

In this case, the fractional-orderq, 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 q, bifurcation diagram of the fractional-order q is shown in Figure 2 with the fractional-order q on the horizontal axis and the ensemble mean response (EMR) x0t (28) on the vertical one. As to the numerical simulation of bifurcation diagram, we take 50 forcing periods T=2π/w as transient value and 50 periods are represented. Also, Δt=0.01 and α-=0.5,f=8.3,w=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(31)x0=y0=1.0;x00=y00=1.0,xi0=yi0=0,i=1,2,3,4.Figure 2
Bifurcation diagram of the fractional-orderq with α-=0.5,f=8.3,w=1.0,σ=0.001.From the bifurcation diagram depicted in Figure2, we can see that the period of the ensemble mean response (EMR) varies from 1T to 2T and then to 4T and finally period-doubling bifurcations lead to chaos as the fractional-order q varies. Then, we choose some given values of the fractional-order q, as 0.984, 0.988 and 0.9895, to validate that the ensemble mean response (EMR) experiences periodic responses with 1T, 2T, and 4T.Whenq=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 1T, 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.Figure 3
Period 1T for q=0.984: (a) phase portrait for DRM-SRM; (b) time history diagram for DRM-SRM; (c) phase portrait for EMR-MCS; (d) time history diagram for EMR-MCS.
(a)
(b)
(c)
(d)Changingq from 0.984 to 0.988, we can see that there is a steady-state periodic EMR with period 2T 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.Figure 4
Period 2T for q=0.988: (a) phase portrait for DRM-SRM; (b) time history diagram for DRM-SRM; (c) phase portrait for EMR-MCS; (d) time history diagram for EMR-MCS.
(a)
(b)
(c)
(d)Further increasingq to 0.9895, the system has a steady-state periodic EMR with period 4T (see Figures 5(c) and 5(d)). Furthermore, we verify the effectiveness of Chebyshev polynomial approximation by contrasting the consistencies over DRM-SRM and EMR-MCS, respectively, as depicted in Figures 5(a), 5(b), 5(c), and 5(d).Figure 5
Period 4T for q=0.9895: (a) phase portrait for DRM-SRM; (b) time history diagram for DRM-SRM; (c) phase portrait for EMR-MCS; (d) time history diagram for EMR-MCS.
(a)
(b)
(c)
(d)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-orderq, and period-doubling bifurcation takes place as the fractional-order q varies from 0.984 to 0.9895, followed by a set of periodic responses with period 1T, 2T, and 4T. That is to say, the fractional-order q 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].


### 4.2.σ 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) x0t (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,f=8.3, w=1.0, and q=0.984 are taken as the system parameters, and the time step size and the initial conditions are the same as the above.Figure 6
Bifurcation diagram of the random intensityσ with α-=0.5, f=8.3, w=1.0, and q=0.984.From the bifurcation diagram depicted in Figure6, it can be observed that the period of the ensemble mean response (EMR) varies from 1T to 2T 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 1T and 2T.Forσ=0.01, in Figure 7 it can be observed that the system has a steady-state periodic EMR with period 1T, which is consistent with the bifurcation diagram in Figure 6. Moreover, the consistency over EMR-MCS (illustrated in Figures 7(a) and 7(b)) reveals the effectiveness of Chebyshev polynomial approximation.Figure 7
Period 1T for σ=0.01: (a) phase portrait for EMR-MCS; (b) time history diagram for EMR-MCS.
(a)
(b)Increasingσ to 0.06, we can see that there is a steady-state periodic EMR with period 2T 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) and 8(b).Figure 8
Period 2T for σ=0.06: (a) phase portrait for EMR-MCS; (b) time history diagram for EMR-MCS.
(a)
(b)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 1T and 2T. 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.


## 4.1.q as Bifurcation Parameter

In this case, the fractional-orderq, 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 q, bifurcation diagram of the fractional-order q is shown in Figure 2 with the fractional-order q on the horizontal axis and the ensemble mean response (EMR) x0t (28) on the vertical one. As to the numerical simulation of bifurcation diagram, we take 50 forcing periods T=2π/w as transient value and 50 periods are represented. Also, Δt=0.01 and α-=0.5,f=8.3,w=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(31)x0=y0=1.0;x00=y00=1.0,xi0=yi0=0,i=1,2,3,4.Figure 2
Bifurcation diagram of the fractional-orderq with α-=0.5,f=8.3,w=1.0,σ=0.001.From the bifurcation diagram depicted in Figure2, we can see that the period of the ensemble mean response (EMR) varies from 1T to 2T and then to 4T and finally period-doubling bifurcations lead to chaos as the fractional-order q varies. Then, we choose some given values of the fractional-order q, as 0.984, 0.988 and 0.9895, to validate that the ensemble mean response (EMR) experiences periodic responses with 1T, 2T, and 4T.Whenq=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 1T, 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.Figure 3
Period 1T for q=0.984: (a) phase portrait for DRM-SRM; (b) time history diagram for DRM-SRM; (c) phase portrait for EMR-MCS; (d) time history diagram for EMR-MCS.
(a)
(b)
(c)
(d)Changingq from 0.984 to 0.988, we can see that there is a steady-state periodic EMR with period 2T 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.Figure 4
Period 2T for q=0.988: (a) phase portrait for DRM-SRM; (b) time history diagram for DRM-SRM; (c) phase portrait for EMR-MCS; (d) time history diagram for EMR-MCS.
(a)
(b)
(c)
(d)Further increasingq to 0.9895, the system has a steady-state periodic EMR with period 4T (see Figures 5(c) and 5(d)). Furthermore, we verify the effectiveness of Chebyshev polynomial approximation by contrasting the consistencies over DRM-SRM and EMR-MCS, respectively, as depicted in Figures 5(a), 5(b), 5(c), and 5(d).Figure 5
Period 4T for q=0.9895: (a) phase portrait for DRM-SRM; (b) time history diagram for DRM-SRM; (c) phase portrait for EMR-MCS; (d) time history diagram for EMR-MCS.
(a)
(b)
(c)
(d)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-orderq, and period-doubling bifurcation takes place as the fractional-order q varies from 0.984 to 0.9895, followed by a set of periodic responses with period 1T, 2T, and 4T. That is to say, the fractional-order q 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].


## 4.2.σ 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) x0t (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,f=8.3, w=1.0, and q=0.984 are taken as the system parameters, and the time step size and the initial conditions are the same as the above.Figure 6
Bifurcation diagram of the random intensityσ with α-=0.5, f=8.3, w=1.0, and q=0.984.From the bifurcation diagram depicted in Figure6, it can be observed that the period of the ensemble mean response (EMR) varies from 1T to 2T 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 1T and 2T.Forσ=0.01, in Figure 7 it can be observed that the system has a steady-state periodic EMR with period 1T, which is consistent with the bifurcation diagram in Figure 6. Moreover, the consistency over EMR-MCS (illustrated in Figures 7(a) and 7(b)) reveals the effectiveness of Chebyshev polynomial approximation.Figure 7
Period 1T for σ=0.01: (a) phase portrait for EMR-MCS; (b) time history diagram for EMR-MCS.
(a)
(b)Increasingσ to 0.06, we can see that there is a steady-state periodic EMR with period 2T 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) and 8(b).Figure 8
Period 2T for σ=0.06: (a) phase portrait for EMR-MCS; (b) time history diagram for EMR-MCS.
(a)
(b)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 1T and 2T. 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.


## 5. 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 period-doubling 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, symmetry-breaking bifurcation, and chaos.



---

*Source: 4162363-2017-04-16.xml*