Stochastic Bifurcations of a Fractional-Order Vibro-Impact System Driven by Additive and Multiplicative Gaussian White Noises

Stochastic fractional-order systems or stochastic vibro-impact systems can present rich dynamical behaviors, and lots of studies dealing with stochastic fractional-order systems or stochastic vibro-impact systems are available now, while the discussion on the stochastic systems with both vibro-impact factors and fractional derivative element is rare. ­is paper is concerned with the stochastic bifurcation of a fractional-order vibro-impact system driven by additive and multiplicative Gaussian white noises. Firstly, we can remove the discontinuity of the original system with the help of nonsmooth transformation and obtain the equivalent stochastic system.­en, we adopt the stochastic averagingmethod to get the approximately analytical solutions. At last, an example is discussed in detail to assess the reliability of the developed approach. We also nd that the coecient of restitution factor, fractional derivative coecient, and fractional derivative order can induce the stochastic bifurcation.


Introduction
As the fractional-order models can more accurately describe the complex systems than the integer-order models, the investigation on the fractional-order systems attracts more and more attention.Fractional calculus enables us to understand the inherent complexity of the real world [1,2] by a new mathematical tool.Many excellent books [3,4] and articles [5][6][7][8][9][10][11][12][13] about fractional calculus are available.Stochastic perturbations are ubiquitous in the real world, so it is necessary to study the dynamical behaviors of the fractionalorder stochastic systems.A lot of methods have been put forward to study the fractional-order stochastic systems, such as the stochastic averaging method [14][15][16][17], multiple scales method [18][19][20], Wiener path integral technique [21], and statistical linearization-based technique [22].Some recent articles on this topic are as follows.Yang et al. [23] investigated the aperiodic stochastic resonance in a bistable fractional-order system induced by the fractional order and the noise intensity.Li et al. [24] estimated the reliability of stochastic dynamical systems under random excitations with a fractional-order proportional integral derivative controller.Denoël [25] carried out the multiple timescale spectral analysis of a noisy single-degree-of-freedom system with a fractional derivative constitutive term.Wang et al. [26] studied the global dynamics of fractional-order systems with the help of the short memory principle and generalized cell mapping method.Di Matteo et al. [27] developed a Galerkin scheme-based approach to determine the survival probability and rst-passage probability of a hysteretic system endowed with fractional derivative elements under Gaussian white noise.Li et al. [28] considered the bifurcation control of a Van der Pol oscillator using the fractional-order PID controller.However, these studies focused on the dynamical behaviors of smooth systems instead of nonsmooth systems.
Vibro-impact systems [29] as the typical nonsmooth systems can present rich dynamical behaviors because of its strong nonlinearity, so a great deal of attention has been devoted to this topic.e stochastic average method was adopted by many authors such as Huang et al. [30], Feng et al. [31,32], and Namachchivaya and Park [33].By comparing the stability domains of P-bifurcation and Dbifurcation, Kumar et al. [34] concluded that these bifurcations need not occur in same regimes.Wang et al. [35] put forward a new procedure based on the generalized cell mapping (GCM) method to explore the stochastic response of vibro-impact systems numerically.Based on Zhuravlev-Ivanov transformation and the iterative method of weighted residue, Chen et al. [36] proposed a new method to obtain the closed-form stationary solution of the vibroimpact system under Gaussian white noise excitation.Although much attention was devoted to the study of vibroimpact systems, little work focused on the investigation of vibro-impact systems with the fractional derivative damping under random excitation.So, in this paper, we will explore the response of a fractional-order vibro-impact oscillator driven by additive and multiplicative Gaussian white noises.
is paper is organized as follows.In Section 2, nonsmooth coordinate transformations are adopted to simplify the fractional-order vibro-impact oscillator.In Section 3, the detailed process to get the analytical solutions is presented.In Section 4.1, an example of fractional-order vibro-impact systems driven by additive and multiplicative Gaussian white noises is discussed in detail to assess the reliability of the developed approach.In Section 4.2, the stochastic bifurcations induced by the system parameters are exhibited.
e conclusions are given in Section 5.

System Description and Its Simplification
e motion equations of a fractional-order vibro-impact system under additive and multiplicative random excitations are of the following form: where ε > 0 is a small constant; β and ω 0 are constant coefficients; 0 < r ≤ 1 is the coefficient of restitution factor; _ x − and _ x + stand for the instant velocities just before and after collision, respectively; ξ 1 (t) and ξ 2 (t) are Gaussian white noises whose statistical properties are of the following form: D α x refers to the fractional derivative element in the Riemann-Liouville sense: where α is the fractional derivative order.
According to Refs.[32,37], the following equation is obtained: According to Refs.[32,37,39], we have , and we can get the following equivalent oscillator without impact term:

Stochastic Averaging Method
As ε > 0 is a small constant, β is a constant coefficient, so the fractional derivative term εβD α x is also small.e lightly damped oscillator ( 6) is subjected to weak random excitations; according to the stochastic averaging method [14], we can assume the solution of equation ( 6) as where A(t), Φ(t), and Ψ(t) are random processes.Substituting equation ( 8) into equation ( 7) and according to Ref. [14], we can obtain the equations for the amplitude A(t) and the phase angle Ψ(t): where (10) en, the averaged Itô equation for the limited process en, the most important step is the calculation of the first term of equation ( 12), i.e., Substituting equation (3) into equation ( 14), we have e Fourier series of the absolute value of the cosine function is where B n � (4/π)((− 1) n /(1 − 4n 2 )).
In order to smooth the solution, substituting equation ( 16) into equation (15), we have According to Refs.[40,41] and equations ( 15) and ( 17), equation ( 14) can be simplified as e other parts of equation (10) and the averaged diffusion coefficient σ 2 (A) can be obtained through mathematical calculation.

Complexity
where C is a normalization constant.According to Ref. [14], the joint stationary probability density function of the displacement y and velocity _ y is as follows: e stationary PDF of the variables x 1 and x 2 can be obtained as e marginal stationary probability density functions p(x 1 ) and p(x 2 ) can be achieved as

Example
e motion equations we consider are expressed as where β, c 4 , c 2 , c 0 , and ω 0 are small constant coefficients.After introducing the nonsmooth coordinate transformations, we have e averaged drift coefficients and diffusion coefficients in equation (11) are According to equations ( 20) and ( 22)-( 24), we can obtain p(A), p(x 1 ), p(x 2 ), and p(x 1 , x 2 ).
(30) e relative error is Comparing equation ( 29) with equation ( 30), we can conclude that keeping more items indeed can improve the accuracy.From equation (31), when ω 0 � 1 and α � 0.5, the relative error is only 0.3343%.So, it is reasonable to keep the first 20 terms when dealing with the series.

Effectiveness of the Method.
In this section, the accuracy of the proposed method will be verified by comparison with the Monte Carlo simulation results.e solid lines are the analytical numerical results, while the discrete dots are the numerical results.We can obtain the analytical solution by substituting m(A) and σ 2 (A) into equation (20).By using the fourth-order Runge-Kutta algorithm, we can obtain numerical results from the original equation (25).
In Figure 1, a comparison between the numerical results and the analytical results is represented.
e system parameter values are listed in Table 1.A very good agreement can be found.So, the effectiveness of the proposed method is acceptable.
In order to further assess the effectiveness of the developed method, another comparison between the numerical results and the analytical results is carried out as shown in Figure 2. e system parameter values are listed in Table 2.A very good match between the numerical and the analytical results indicates that the developed procedure is effective.

Bifurcation Analysis.
In Section 4.1, we demonstrated the effectiveness of the proposed method.In this section, we turn our attention to the stochastic P-bifurcation induced by system parameters.As the investigation on the stochastic Pbifurcation enabled us to have a more clear understanding on the dynamical behavior of the system, especially on the long-run probability distributions, we will conduct the bifurcation analysis in this section.In this paper, the stochastic P-bifurcation or phenomenological bifurcation takes place when the structure of stationary probability density function has qualitative changes as parameters are varied.
First, we investigate the influence of the coefficient of restitution factor r on the stochastic bifurcation.e system parameter values are listed in Table 3. Figure 3 depicts the joint probability density functions for different r.It can be concluded that increasing the coefficient of restitution factor  Complexity r from 0.984 to 0.989 gives rise to the occurrence of stochastic P-bifurcation.Specifically, when r � 0.984, the joint probability density function has one peak, while when r � 0.989, the joint probability density function presents a crater-like structure.e qualitative transformation of the probability density function indicates the occurrence of stochastic P-bifurcation.In order to better understand the progress of the stochastic P-bifurcation, the corresponding section graphs of probability density functions when x 1 � 0.05 are presented in Figure 4.

Complexity
Second, we explore the influence of the fractional derivative coefficient β on the stochastic bifurcations.
e system parameter values are listed in Table 4. Figure 5 depicts the joint probability density functions for different β. Figure 6 presents the corresponding section graphs of joint probability density functions when x 1 � 0.05.It can be observed that decreasing fractional derivative coefficient β from 0.009 to 0.003 leads to the occurrence of stochastic Pbifurcation.ird, we discuss the influence of the fractional derivative order α on the stochastic bifurcations.e system parameter values are listed in Table 5. Figure 7 depicts the joint probability density functions for different α. Figure 8 presents the corresponding section graphs of joint probability density functions when x 1 � 0.05.According to similar analysis, it can be observed that decreasing fractional derivative order α from 0.6 to 0.3 contributes to the occurrence of stochastic P-bifurcation.

Conclusions
We carried out the investigation on the stochastic bifurcation of a fractional-order vibro-impact system under additive and multiplicative Gaussian white noise excitations.
ere are two challenges to study the stationary response of the fractional vibro-impact systems under Gaussian white    noises.e first one is how to deal with the discontinuity of the original system.e second one is how to get the explicit expression of the averaged drift coefficient when we utilize the stochastic averaging method.ese two challenges have been solved in this paper by the nonsmooth transformations and stochastic averaging method.An example is discussed in detail to assess the reliability of the developed approach.e results showed that the proposed method has a satisfactory accuracy.We also found that the coefficient of restitution factor, fractional derivative coefficient, and fractional derivative order can be treated as bifurcation parameters.

Figure 2 :
Figure 2: A comparison between the analytical solutions (solid blue lines) and the numerical results (discrete blue dots).

4 Complexity Table 1 :
Parameter values used in simulation.

Table 2 :
Parameter values used in simulation.

Table 3 :
Parameter values used in simulation.

Table 5 :
Parameter values used in simulation.