Stochastic P-bifurcation in a tri-stable Van der Pol system with fractional derivative under Gaussian white noise

1. Introduction

Fractional calculus is a generalization of integer-order calculus that has been used for over 300 years to extend the order of calculus operation from traditional integer order to the case of non-integer order. Integer-order derivatives cannot express the memory property of viscoelastic substances. The definition of fractional derivatives contain convolution, which expresses a memory effect, showing a cumulative effect over time. Compared with traditional integer-order calculus, fractional calculus is a more suitable mathematical tool for describing memory characteristics [1-4], and it has become a powerful mathematical tool to study fields such as anomalous diffusion, non-Newtonian fluid mechanics, viscoelastic mechanics, and soft matter physics. Since the fractional derivative can describe various reaction processes more accurately, many problems can be better described by fractional differential equations [5-11]. The classical integer operators do not have sufficient parameters to handle the different shapes of the hysteresis loops describing the behaviors of viscoelastic materials and structures, while the fractional derivative can accurately describe the constitutive relation of viscoelastic materials with fewer parameters, so the studies of fractional differential equations on the typical mechanical properties and the influences of fractional order parameters on the system are very necessary and significant.

In recent years, many scholars have studied the dynamic behavior of nonlinear multi-stable systems under different noise excitations and achieved fruitful results. For integer-order systems, Zhu et al. put forward a magnetic shape memory alloy microgripper and studied its nonlinear dynamic characteristics under stochastic perturbation [12]. Guo et al. investigated the asymptotic behaviors of the solution for the Hasegawa-Mima equation driven by external noise, and the existence of the random attractor for the corresponding random dynamical system is proved [13]. The authors studied the stationary PDF of amplitude of Duffing-Van der Pol oscillator under levy noise, color noise, combined harmonic and random noise in Refs. 14-18 respectively, and discussed the stochastic P-bifurcation of the noise oscillator by analyzing the changes of stationary PDF of the systems, as well as obtained the analytical results of the bimodal stationary PDF of the systems, it showed that both the system parameters and the noise intensity can induce stochastic P-bifurcation of the systems. Hao and Wu studied the stochastic P-bifurcation of a tri-stable Van der Pol-Duffing oscillator excited by Gaussian white noise, multiplicative colored noise, as well as the combined additive and multiplicative Gaussian white noise, and obtained the analytical expressions of the stationary PDF of the system amplitude, then analyzed the influences of system parameters and noise intensity on the stochastic P-bifurcation of the system [19-21]. Deng et al. proposed a new motor bearing fault diagnosis method based on integrating fuzzy entropy and support vector machine and comparing fault feature extraction by using empirical wavelet transform and Hilbert transform with the theoretical calculation [22]. Deng et al. improved particle swarm optimization and propose a novel intelligent diagnosis method aiming at the problem that the most existing fault diagnosis methods could not effectively recognize the early faults in the rotating machinery, the empirical mode decomposition and fuzzy information entropy [23]. Zhao et al. proposed a new fault feature extraction method based on integrating ensemble empirical mode decomposition, mode selection and multi-scale fuzzy entropy to accurately diagnose fault, the method is used to decompose the vibration signal into a series of intrinsic mode functions with a different physical significance [24]. Deng et al. proposed a genetic and ant colony adaptive collaborative optimization algorithm to overcome the deficiencies of weak local search ability in genetic algorithms and slow global convergence speed in ant colony optimization algorithm in solving complex optimization problems [25]. Deng et al. proposed an efficient multi-objective optimization model of gate assignment problem and studied an improved adaptive particle swarm optimization algorithm based on making full use of the advantages of Alpha-stable distribution and dynamic fractional calculus [26]. Zhao et al. proposed an enhanced empirical wavelet transform based on maximum-minimum length curve method to realize fault diagnosis of motor bearings and verified the effectiveness of the method for fault diagnosis [27]. Deng et al. proposed an improved chaotic ant colony optimization algorithm based on adaptive multi-variant strategies, improving the search efficiency, search accuracy and robustness of the algorithm [28]. Deng et al. presented a novel two-stage hybrid swarm intelligence optimization algorithm and a novel parallel hybrid intelligence optimization algorithm that combine the evolution ideas of the genetic algorithms, particle swarm optimization and ant colony optimization based on the compensation for solving the traveling salesman problem and optimizing the structure and parameters of the fuzzy neural network [29, 30]. Zhao et al. proposed a damage degree identification method based on high-order difference mathematical morphology gradient spectrum entropy to solve the problem that fault signal of rolling bearings are weak and difficult to be quantitatively measured and testified the effectiveness of the proposed method using the vibration signal of motors under no-load and load states [31]. Deng et al. proposed an improved self-adaptive differential evolution algorithm with multiple strategies algorithm using a different search strategy and a parallel evolution mechanism, improved the robustness of the algorithm [32].

For the fractional systems, Shu et al. considered a stochastic coupled fractional Ginzburg-Landau equation under additive noise and proved the existence of random attractor for the random dynamical system [33]. Sandev et al. studied the generalized fractional Langevin equations in the presence of a harmonic potential and investigated the asymptotic behavior of the generalized stochastic oscillator as well as considered the case of high viscous damping [34]. Chen et al. studied the response of a Duffing oscillator with fractional derivative damping under combined harmonic and white noise excitations and showed that variations of the order of fractional derivative can cause stochastic P-bifurcation in the system [35]. Huang et al. discussed the response and the stationary PDF of a single degree of freedom nonlinear system under Gaussian white noise excitation [36]. Yang et al. studied the stochastic bifurcations in the nonlinear vibroimpact system with Caputo fractional derivative under Gaussian white noise excitation by using the Zhuravlev non-smooth transformation and stochastic averaging [37]. Li et al. studied the bistable stochastic bifurcation of the generalized Duffing-Van der Pol system with fractional derivative under additive and multiplicative colored noise excitations respectively, and found that the linear damping coefficient, the order of fractional derivative and the noise intensity can all lead to stochastic P-bifurcation in the system [38].

Due to the complexity of fractional derivatives, analyzing them is difficult, the vibration characteristics of the parameters can only be qualitatively analyzed, and critical conditions of parametric influences cannot be found. These problems affect the analysis and design of such systems, in part because the stochastic P-bifurcation of tri-stability for the fractional order system has not been reported. Accordingly, we take the nonlinear vibration of a generalized Van der Pol oscillator under white noise excitation as an example, giving the concrete form of the governing equation of the nonlinear fractional oscillator, then obtain the transition set curves of the fractional order system and critical parametric conditions for stochastic P-bifurcation of the system by the singularity method, then analyze the types of stationary PDF curves of the system in each region in the parameter plane. We compare numerical results from a Monte Carlo simulation with the analytical results obtained in this paper; it can be seen that the numerical solutions are in good agreement with the analytical solutions, verifying our theoretical analysis.

2. Equivalent Van der Pol system

Among the many definitions of fractional derivatives, the Riemann-Liouville derivative and Caputo derivative are most commonly used. However, while the initial conditions corresponding to the Riemann-Liouville derivative have no physical meaning, the initial conditions of systems described by the Caputo derivative have both clear physical meanings and form the same as the initial conditions for differential equations of integer order. Therefore, we adopt the Caputo-type fractional derivative as follows:

where $m-1<p\le m$, $m\in N$, $t\in [a,b]$, $\mathrm{\Gamma }\left(m\right)$ is the Euler Gamma function, and $x^{\left(m\right)}\left(t\right)$ is the m order derivative of $x\left(t\right)$.

For a given physical system, because the initial moment of the oscillator is $t=$ 0, the following form of the Caputo derivative is often used:

where $m-1<p\le m$, $m\in N$.

In this paper, we studied the generalized Van der Pol oscillator system with fractional derivative under Gaussian white noise excitation:

where $\epsilon$ represents the linear damping coefficient ${\alpha }_1$, ${\alpha }_2$, ${\alpha }_3$ and ${\alpha }_4$ represent the nonlinear damping coefficients of the system, and $w$ is the natural frequency of the system. ${}_0{}^C\mathrm{}{D}^p\left[x\right(t\left)\right]$ is the $p$ ($0\le p\le 1$) order Caputo derivative of $x\left(t\right)$, defined by Eq. (2). $\xi \left(t\right)$ is the Gaussian white noise, which satisfies:

where $\delta \left(\tau \right)$ is the Dirac function.

The fractional derivative term has contributions from both damping and restoring forces [38-42], hence, we introduce the following equivalent system:

where $c\left(a\right)$, $w^2\left(a\right)$ are the coefficients of equivalent damping and restoring forces of the fractional derivative ${}_0{}^C\mathrm{}{D}^{p}x$, respectively.

The difference between systems Eqs. (3) and (5) is:

The necessary conditions for minimal mean square error are [43]:

Substituting Eq. (6) into Eq. (7) yields:

Assume that the solution of system Eq. (3) has the form:

where $\phi \left(t\right)=wt+\theta \left(t\right)$.

Based on Eq. (9), we can obtain:

Substituting Eqs. (9), (10) into Eq. (8) yields:

To simplify Eqs. (11), (12) further, asymptotic integrals are introduced as follows:

Substituting Eq. (13) into Eqs. (11), (12), and averaging them across$\phi$produces the ultimate forms of $c\left(a\right)$ and $w^2\left(a\right)$ as follows:

Therefore, the equivalent Van der Pol oscillator associated with system (5) can be expressed as follows:

where:

3. The stationary PDF of amplitude

In our first example, we examine the system in Eq. (15), with linear and nonlinear damping coefficients $\epsilon =$ –0.1, ${\alpha }_1=$ 1.51, ${\alpha }_2=$ 2.85, ${\alpha }_3=$ 1.693, ${\alpha }_4=$ 0.312, natural frequency $w=$ 1. For convenience in discussing parametric influence, the bifurcation diagram of the limit cycle amplitude along with variation of the fractional order $p$ is shown in Fig. 1 when $D=$ 0.

It can be seen that there are 2 attractors in the system where 0.0637 $\le p\le$ 0.1377 (equilibrium and limit cycle, shown in Fig. 2(a)); there are 3 attractors where 0.1378 $\le p\le$ 0.1406 (equilibrium, small and large limit cycles, as shown in Fig. 2(b)); there are also 2 attractors where 0.1407 $\le p\le$ 0.1522 (equilibrium and small limit cycle, as shown in Fig. 2(c)).

In order to obtain the stationary PDF of the amplitude of system Eq. (15), we introduce the following transformation [44]:

where $a\left(t\right)$,$\theta \left(t\right)$ represent the amplitude process and the phase process of the system response, respectively. These are both random processes.

Fig. 1Bifurcation diagram of the deterministic system (with variation in p)

Fig. 2Phase diagrams of the deterministic system (for 3 different values of p)

a) Equilibrium and large limit cycle

b) Tri-stable state

c) Equilibrium and small limit cycle

Substituting Eq. (17) into Eq. (15), we can obtain:

in which:

Eq. (18) can be regarded as a Stratonovich stochastic differential equation. By adding the corresponding Wong-Zakai correction term, we transform it into the following Itô stochastic differential equation:

where $B\left(t\right)$ is a standard Wiener process and:

By the stochastic averaging method [45], averaging Eq. (20) regarding $\mathrm{\Phi }$, the following averaged Itô equation can be obtained as follows:

where:

Substituting Eq. (19) into Eq. (23), we obtain the explicit expression for the averaged drift and diffusion coefficients as follows:

Eqs. (22)-(24) show that the averaged Itô equation for $a\left(t\right)$ is independent of $\theta \left(t\right)$, the random process $a\left(t\right)$ is a one-dimensional diffusion process. The corresponding FPK equation associated with $a\left(t\right)$ can be wrote as:

The boundary conditions are:

Thus, based on these boundary conditions Eq. (26), we can obtain the stationary PDF of amplitude:

where $C$ is a normalized constant, which satisfies:

Substituting Eq. (24) into Eq. (27), we obtain the explicit expression for the stationary PDF of amplitude $a$:

in which:

4. Stochastic P-bifurcation

Stochastic P-bifurcation refers to the changes in the number of peaks in the PDF curve. To obtain the critical parametric conditions for stochastic P-bifurcation, we analyze the influence of parameters on stochastic P-bifurcation of the system using singularity theory.

For the sake of convenience, we can write $p\left(a\right)$ as follows:

in which:

According to singularity theory [46], the stationary PDF needs to satisfy the following two conditions:

Substituting Eq. (31) into Eq. (33), we can obtain the following condition [20, 38]:

where $H$ represents the condition for the changes of the number of peaks in the PDF curve.

Substituting Eq. (32) into Eq. (34), we can get the critical parametric conditions for stochastic P-bifurcation of the system with respect to fractional order$p$and noise intensity $D$ as follows:

where amplitude $a$ satisfies:

with parameters $\epsilon =$ –0.1, ${\alpha }_1=$ 1.51, ${\alpha }_2=$ 2.85, ${\alpha }_3=$ 1.693, ${\alpha }_4=$ 0.312, $w=$ 1, according to Eqs. (35), (36), we obtain the transition set for stochastic P-bifurcation of the system with unfolding parameters $p$ and $D$ shown in Fig. 3.

Fig. 3Transition set curves of system Eq. (3)

According to singularity theory, the topological structures of the stationary PDF curves of different points ($p$, $D$) in the same region are qualitatively the same. By taking a point ($p$, $D$) in each region, we obtain all of the varieties of the stationary PDF curves that are qualitatively different. The unfolding parameter $p$-$D$ plane is divided into five sub-regions by the transition set curve; for convenience, each region in Fig. 3 is marked with a number.

We analyze the stationary PDF curves for a point ($p,D$) from each of the five sub-regions of Fig. 3, and show the corresponding results in Fig. 4.

Fig. 4Stationary PDF curves for the amplitude of system (3)

a) Parameter ($p$, $D$) in region 1 of Fig. 3

b) Parameter ($p$, $D$) in region 2 of Fig. 3

c) Parameter ($p$, $D$) in region 3 of Fig. 3

d) Parameter ($p$, $D$) in region 4 of Fig. 3

e) Parameter ($p$, $D$) in region 5 of Fig. 3

As can be seen from Fig. 3, the parameter region where the PDF curve appears multimodal is surrounded by two approximately triangular regions. In particular, the region 4 where the two triangular regions are coincident forms a trimodal region of the stationary PDF curve for the system. When the parameter ($p$, $D$) is taken in the region 1, the PDF curve has a distinct peak far away from the origin, as shown in Fig. 4(a); in region 2, the PDF curve has two distinct peaks farther away from the origin, there are both small and large limit cycles in the system, as shown in Fig. 4(b); in region 3, the PDF curve still has a distinct peak farther away from the origin, but the probability is obviously not zero near the origin, there are both the equilibrium and large limit cycle in the system simultaneously, as shown in Fig. 4(c); in region 4, the PDF curve has three peaks, it shows that the equilibrium coexists with the small and large limit cycles in the system which is tri-stable, as shown in Fig. 4(d); in the region 5, the amplitude corresponding to the peak deviating from the origin of the PDF curve is smaller than the corresponding amplitude in Fig. 4(c), there are both the equilibrium and small limit cycle in the system simultaneously, as shown in Fig. 4(e).

The results show that the stationary PDF curves of the system amplitudes in any two adjacent regions in Fig. 3 are qualitatively different. Regardless of the exact values of the unfolding parameters, if they cross any line in the figure, the system will demonstrate stochastic P-bifurcation behavior. Therefore, the transition set curves are just the critical parametric conditions for the stochastic P-bifurcation of the system. The analytic results shown in Fig. 4 are in good agreement with the numerical results obtained by Monte Carlo simulation, further verifying the theoretical analysis and showing that it is feasible to use the methods in this paper to analyze the stochastic P-bifurcation behavior of fractional order system.

5. Conclusions

In this paper, we studied the tri-stable stochastic P-bifurcation of a generalized Van der Pol system with fractional derivative under Gaussian white noise excitation. According to the minimal mean square error principle, we transformed the original system into an equivalent integer-order system and obtained the stationary PDF of the system amplitude using stochastic averaging. In addition, we obtained the critical parametric conditions for stochastic P-bifurcation of the system using singularity theory, based on this, the system response can be maintained at the small amplitude near the equilibrium by selecting the appropriate unfolding parameters, providing theoretical guidance for system design in practical engineering and avoiding the instability and damage caused by the large amplitude vibration or nonlinear jump phenomenon of the system. The numerical results by Monte Carlo simulation of the original system also verify the theoretical results obtained in this paper. We conclude that the fractional order $p$ and noise intensity $D$ can each cause the stochastic P-bifurcation behavior of the system, and the number of peaks in the stationary PDF curve of system amplitude can change from one to three by selecting the appropriate unfolding parameters ($p$, $D$), it also shows that the method used in this paper is feasible to analyze the stochastic P-bifurcation behaviors of the nonlinear oscillator with fractional derivative.