Parametric influences of a modified fractional Van der Pol oscillator under Gaussian colored noise

Abstract

We investigate the bistable dynamic behaviors of a fractional and modified Van der Pol system subjected to Gaussian colored noise, especially focusing on the control of bifurcation parameters. Firstly, according to the minimal mean square error concept and the generalized harmonic balance technique, the fractional derivative term can be expressed as a linear combination of the damping force and restoring force. Consequently, the initial system can be simplified as an equivalent integer-order Van der Pol system. And the Probability Density Function (PDF) of the system’s stationary amplitude can be obtained by utilizing the stochastic averaging method. Subsequently, the critical condition for stochastic P-bifurcation of system amplitude is determined by using the singularity theory. Finally, the qualitative analysis of the stationary PDF curves of amplitude is discussed in each region divided by the transition set curves. The agreement between the analytical solutions and the numerical outcomes obtained through Monte Carlo simulation confirms the theoretical analysis presented in this paper. Additionally, the results attained in this paper have a potential guidance to control the system’s response by optimizing the parametric design of the fractional-order controller.

Introduction

Fractional calculus is an extension of calculus with non-integer orders, and it can capture the memory properties of viscoelastic substances compared with integer-order derivatives. As the expression of fractional derivative involves convolution, which enables to represent the memory and cumulative effects over time. Hence, the fractional derivative has been proved to be a better and more precise mathematical instrument for characterizing memory properties^1–4 and has been widely applied in various domains including anomalous diffusion, non-Newtonian fluid mechanics, soft matter physics, viscoelastic mechanics, and diverse reaction processes.^5–7 Therefore, considering the above characters and the ubiquitous ambient noise in engineering, it is crucial and meaningful to investigate the impacts of fractional-order parameters and noise excitation on the dynamic properties of the stochastic systems.

In general, additive noise comes from the internal fluctuations of the system, while multiplicative noise comes from the fluctuations of the external environment of the system, and in most cases, multiplicative and additive noise are interrelated and may have a great influence on system dynamics. Cao et al. ⁸ studied the bistable system under the action of interrelated noise and gave the corresponding uniform colored noise approximation. Li et al. ⁹ studied the non-equilibrium phase transition under correlated noise and found that the mutual correlation between multiplicative and additive noise plays an important role in the phase transition of the system. Mei et al.^10–12 studied the nonlinear dynamic system under the action of cross-correlation noise, and found that the correlation between noises could induce resonance and suppression phenomena in the curve of average first crossing time. Jin et al.^13,14 studied the stochastic resonance of linear systems under the combined action of dependent multiplicative colored noise and periodic modulation noise, and found the real stochastic resonance and the traditional stochastic resonance phenomenon. Jin et al. ¹⁵ studied the random bifurcation of bistable Duffing–Van der Pol systems under the action of correlated multiplicative and additive Gaussian white noises, and calculated the steady-state probability density function (PDF) curve of system amplitude by using the stochastic averaging method, then discussed the influence of noise intensity, noise cross-correlation coefficient, and system parameters on the stochastic bifurcation of the system.

There are random excitations in natural sciences, engineering sciences, and social sciences, such as strong winds, road or track irregularities, ground movements caused by strong earthquakes, etc. The response of a system under random excitations is fundamentally different from that under deterministic excitations, which can only be described by probability or statistical methods. Therefore, in addition to studying the characteristics of dynamic systems under deterministic excitations, it is necessary to pay attention to the characteristics of dynamic systems under random excitations, such as the vibration of engineering structures in earthquakes,^16,17 the vibration of ships and offshore platforms caused by wind and waves,^18,19 the vibration of aircraft caused by jet noise and high-altitude turbulence,^20,21 and the vibration of high-rise buildings and large bridges caused by wind and earthquakes.^22–24 Therefore, in order to accurately grasp the dynamics of the system under random excitation, it is necessary to study the theory of stochastic dynamics.

And in fact, all practical systems exhibit nonlinear characteristics to some extent. In these nonlinear systems, random factors exist widely, such as high-altitude turbulence encountered by wings, ²⁵ excitation from tracks encountered by trains, ²⁶ and excitation from wind waves, ocean currents, and internal waves in ocean structures.^27,28 Under the action of random excitation, the system response will switch between different steady states, and the steady-state probability density of system response will have multiple peaks.

In recent years, numerous researchers have investigated the dynamic characteristics of nonlinear systems with multistable states under various noise excitations and obtained productive outcomes. Baleanu et al. ²⁹ studied the fractional Lagrangian of Pais–Uhlenbeck oscillator, and obtained Euler–Lagrangian equation numerically based on the Grünwald–Letnikov approach. Rath et al. ³⁰ investigated the characteristic features of Yao–Cheng non-linear oscillator by utilizing both the analytical and numerical approaches and indicated that the solution of the non-linear oscillator showed the oscillatory behavior. Baleanu et al. ³¹ studied the dynamical behavior of the motion for a simple pendulum with a mass decreasing exponentially in time, derived the corresponding non-integer Euler–Lagrange equation, and then simulated the approximate results with respect to different non-integer orders. A number of studies^32–38 discussed the stochastic behaviors of the Van der Pol–Duffing oscillators driven by Gaussian white noise, Lévy noise, colored noise, and the combination of harmonic and random noise, respectively, the researchers derived the expression for the stationary PDF of the system’s amplitude and explored the stochastic P-bifurcation behaviors of the systems. Chen et al. ³⁹ utilized the radial basis function neural networks (RBFNN) method to analyze the transient response of randomly excited non-smooth vibro-impact systems, and expressed the solution of the system response as a serious of Gaussian activation functions with time-dependent weights. Luan and Huang ⁴⁰ discussed random stationary responses of nonlinear oscillators including Preisach hysteresis by the Recurrent Neural Networks method, and solved the stationary probability density from the FPK equation as well as various-order response statistics. Liu et al. ⁴¹ presented a novel method to identify the high dimensional stochastic dynamical systems driven by a non-Gaussian $\alpha$ -stable Lévy noise and illustrated three examples to demonstrate the feasibility. Xu et al. ⁴² explored the complicated dynamical behaviors of a conceptual airfoil excited by an external harmonic force and an extreme random load, and found that the extreme random load can lead to stochastic transition and stochastic resonance behaviors. He et al. ⁴³ investigated a nonlinear Helmholtz–Duffing oscillator with small amplitude based on the harmonic balance method, and showed that the calculation process could become more convenient by introducing the special functions. Sun and Zhang ⁴⁴ discussed the necessary conditions of system resonance by analyzing the displacement response of a multi-degree-of-freedom dynamic system, and found that its displacement response should present the corresponding modal shape while the natural frequency of the system is equal to the excitation frequency. Xu et al. ⁴⁵ have taken a class of generalized Duffing-type system as an example to reveal the three mechanisms for the occurrence of extreme events, and the results provided a guidance for further prediction as well as the avoidance of extreme events. Zhang et al. ⁴⁶ investigated the stochastic dynamics of an asymmetric tri-stable hybrid vibration energy harvester driven by colored noise, and analyzed the influences of the asymmetric parameter as well as the colored noise intensity on the steady state response and stochastic resonance of the system, the results showed that the energy harvesting performance can be effectively enhanced by reducing asymmetry and correlation time. Mohammed et al. ⁴⁷ considered the system of stochastic equations of the ion sound and Langmuir waves with multiplicative noise in the Itô sense, and obtained the new exact stochastic solutions by applying three distinct methods such as: He’s emi-inverse method, the Riccati–Bernoulli sub-ODE method and the sine-cosine method. Alshammari et al. ⁴⁸ investigated the stochastic-fractional Broer–Kaup equations perturbed by the multiplicative Wiener, and acquired the rational, hyperbolic, and elliptic stochastic solutions for SBKEs by utilizing the Jacobi elliptic function method. Al-Askar et al. ⁴⁹ addressed the stochastic Davey–Stewartson equations forced by multiplicative noise, and found the new rational, elliptic, hyperbolic, and trigonometric functions by using the mapping method. Albosaily et al. ⁵⁰ discussed the exact solutions to the stochastic fractional-space Allen–Cahn equation by utilizing the tanh–coth method, and then studied the effect of stochastic term on the solutions of the equation. Al-Askar et al. ⁵¹ considered the stochastic fractional-space Burgers’ equation forced by multiplicative noise, and then, generalized some previously results where this equation was not studied before in the Itô sense, further, discussed the effect of the stochastic term on the stability of the solutions by utilizing graphical representations.

Phenomena such as large deformation and large deflection often occur in engineering problems, in these situations where high computational accuracy is required, then high-order restoring force models and high-order nonlinear damping are needed to describe the effects of nonlinear stiffness and damping of the system,^52–56 leading to the occurrence of multiple steady-state phenomena. Many biological systems, such as the firing activity of biological neurons,^57,58 gene networks, ⁵⁹ and cellular protein transcription,^60,61 have also observed the phenomenon of multiple steady states. From the current research on stochastic systems, it can be seen that the deterministic part of most systems is multistable, and under the action of random excitation, the system response will switch between different steady states, then the steady-state probability density of the stochastic system response will exhibit multiple peaks. However, the mechanisms behind these complex phenomena are still unclear, which hinders a deeper understanding of such dynamical systems.

The Van der Pol equation was first discovered and established by Van der Pol during his research on nonlinear oscillation circuits of transistors. As a typical nonlinear term, Hartlen and Curri ⁶² proposed a control equation that used Van der Pol nonlinear oscillator as lift coefficients to represent wake effects and further studied the effect of cylinder motion on vortex induced lift based on the research of Bishop and Hasson. ⁶³ The model was also applied to investigate the vortex induced vibration responses of cables, etc., and predict the multiple steady-state phenomenon of these systems.

Therefore, as a typical representative of nonlinear multiple steady-state systems, it is necessary to study the dynamic behavior and generation mechanism of complex phenomenon of Van der Pol system under random excitation, which can more clearly determine the steady-state response behavior of stochastic systems under the influence of different factors, facilitating to accurately analyze the influence of system parameters and providing theoretical guidance and potential engineering application value for practical problems.

And as the complexity of the fractional derivative, it is only possible to qualitatively analyze the parametric influences of the fractional system, and can’t obtain the critical values of the parametric effects. However, the practical importance of the critical parameter conditions can’t be ignored when analyzing and designing the fractional order systems. Furthermore, since the multiplicative parameter noise excitation is more complicated than additive external noise excitation in the exploration of system behaviors, we investigate the nonlinear vibration of fractional-order stochastic system with the fractional derivative and noise excitation in this paper. Comparing with the Van der Pol systems ⁶⁴ which considered the effect of additional and multiplicative Gaussian white noises simultaneously, and the bistable Van der Pol systems ⁶⁵ which discussed the influences of additional and multiplicative Gaussian white noises, respectively, in this paper, we investigate a modified and bistable Van der Pol system with fractional element which is excited by additive and multiplicative Gaussian colored noises as an dynamic model. The more complex Gaussian colored noise which can be generated through a first-order low pass filter to the Gaussian white noise, and this means that it should has a correlation time $\tau$ in the noise rather than the Gaussian white noise. Based on the singularity method and the help of stochastic averaging method, the critical parameter conditions for the stochastic P-bifurcation of the system are acquired. And then, all the types of the system’s stationary PDF in each sub-region divided by the transition set curves are analyzed. Finally, we compare the numerical results derived from Monte Carlo simulation with the analytical solutions acquired through the stochastic averaging technique, and the consistency between them confirms the validity of the theoretical analysis process in this paper.

Derivation of the equivalent system

In fact, the Riemann–Liouville derivative and Caputo derivative are most commonly used. However, the initial conditions of Riemann–Liouville derivative hasn’t explicit physical meaning, whereas the initial condition of Caputo derivative not only possesses a distinct physical meaning but also has the identical form as the differential equation of integer order. Hence, we use the Caputo fractional derivative in this study as

$\{\begin{array}{l}{{}_a{}^{C}D}^p\left[x\left(t\right)\right]=\frac{1}{\Gamma \left(m-p\right)}{{\displaystyle \int }}_a^t\frac{x^{\left(m\right)}\left(u\right)}{{\left(t-u\right)}^{1+p-m}}du,m<p<m+1\\ {{}_a{}^{C}D}^p\left[x\left(t\right)\right]=x^{\left(m\right)}\left(t\right),p=m\end{array}$

(1)

where $m\in N$ , $t\in \left[a,b\right]$ , and $x^{\left(m\right)}\left(t\right)$ is the $m$ -th order derivative of $x\left(t\right)$ , and $\Gamma \left(m\right)$ is the Gamma function that fulfills $\Gamma \left(m+1\right)=m\Gamma \left(m\right)$ .

For a given physical system, the starting moment of an oscillator is $t=0$ and the Caputo derivative is generally denoted by

${{}_0{}^{C}D}^p\left[x\left(t\right)\right]=\frac{1}{\Gamma \left(m-p\right)}{{\displaystyle \int }}_0^t\frac{x^{\left(m\right)}\left(u\right)}{{\left(t-u\right)}^{1+p-m}}du$

(2)

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

In this paper, we investigate the modified and generalized Van der Pol system with fractional damping, which is driven by the Gaussian colored noise excitation

$\{\begin{array}{l}\ddot{x}\left(t\right)-u\left(-\epsilon +{\alpha }_1{x}^2\left(t\right)-{\alpha }_2{x}^4\left(t\right)\right)·{{}_0{}^{C}D}^p\left[x\left(t\right)\right]+w^{2}x\left(t\right)=\left[\alpha +\left(1-\alpha \right)x\left(t\right)\right]\zeta \left(t\right)\\ \dot{\zeta }\left(t\right)+\frac{1}{\tau }\zeta \left(t\right)=\frac{1}{\tau }\xi \left(t\right)\end{array}$

(3)

where $\mu$ satisfies $0<\mu <1$ , $x\left(t\right)$ is the displacement of the system, $\epsilon$ represents the coefficient of linear damping, ${\alpha }_1$ and ${\alpha }_2$ stand for the nonlinear damping coefficients of the system, respectively, and $w$ is the

intrinsic frequency of the system. ${{}_0{}^{C}D}^p\left[x\left(t\right)\right]$ is the $p$ -order ( $0<p<1$ ) Caputo fractional derivative of $x\left(t\right)$ , which is defined by the equation (2). And $\xi \left(t\right)$ is the Gaussian white noise excitation, which satisfies

$\{\begin{array}{l}E\left[\xi \left(t_1\right)\xi \left(t_2\right)\right]=2D\delta \left(t_1-t_2\right)\\ E\left[\xi \left(t\right)\right]=0\end{array}$

(4)

where $D$ represents the strength of the white noise, $\alpha \left(0\le \alpha \le 1\right)$ is the weight between the additive and multiplicative noise components ( $\alpha =1$ corresponding to the purely additive Gaussian colored noise and $\alpha =0$ corresponding to the purely multiplicative Gaussian colored noise).

In the meantime, the Gaussian colored noise can be acquired by applying a first-order low pass filter to the Gaussian white noise, and as demonstrated in equation (3), which has zero mean and the auto-correlation function as

$\{\begin{array}{l}E\left[\zeta \left(t_1\right)\zeta \left(t_2\right)\right]=\frac{D}{\tau }\mathit{exp}\left[-\frac{\left|t_1-t_2\right|}{\tau }\right]\\ S\left(w\right)=\frac{D}{\pi \left(1+{\tau }^2{w}^2\right)}\end{array}$

(5)

where the parameters $\tau$ and $D$ denote the correlation time and the strength of colored noise, respectively.

In order to introduce the equivalent system, we consider the fractional derivative as the combination of both damping force and restoring force ⁶⁶ as

$\ddot{x}\left(t\right)-u\left(-\epsilon +{\alpha }_1{x}^2\left(t\right)-{\alpha }_2{x}^4\left(t\right)\right)·\left[C\left(p,w\right)\dot{x}\left(t\right)+K\left(p,w\right)x\left(t\right)\right]+w^{2}x\left(t\right)=\left[\alpha +\left(1-\alpha \right)x\left(t\right)\right]\zeta \left(t\right)$

(6)

where $C\left(p,w\right)$ and $K\left(p,w\right)$ represents the coefficients of the equivalent damping and restoring forces about the fractional derivative ${{}_0{}^{C}D}^p\left[x\left(t\right)\right]$ , respectively.

And applying the minimal mean square error method, ⁶⁶ we can get the ultimate forms of $C\left(p,w\right)$ and $K\left(p,w\right)$

$\{\begin{array}{l}K\left(p,w\right)=w^p\hspace{0.17em}\mathit{cos}\left(p\pi /2\right)\\ C\left(p,w\right)=w^{p-1}\hspace{0.17em}\mathit{sin}\left(p\pi /2\right)\end{array}$

(7)

Hence, the isovalent Van der Pol oscillator with equation (3) can be expressed as

$\ddot{x}\left(t\right)+{w_0}^{2}x\left(t\right)-u{\alpha }_1{w}^p\hspace{0.17em}\mathit{cos}\left(p\pi /2\right)x^3\left(t\right)+u{\alpha }_2{w}^p\hspace{0.17em}\mathit{cos}\left(p\pi /2\right)x^5\left(t\right)-w^{p-1}\hspace{0.17em}\mathit{sin}\left(p\pi /2\right)·\gamma ·\dot{x}\left(t\right)=[\alpha +(1-\alpha )x\left(t\right)]\zeta \left(t\right)$

(8)

where

$\{\begin{array}{l}\gamma =u\left(-\epsilon +{\alpha }_1{x}^2\left(t\right)-{\alpha }_2{x}^4\left(t\right)\right)\\ {w_0}^2=w^2+u\epsilon w^p\hspace{0.17em}\mathit{cos}\left(p\pi /2\right)\end{array}$

(9)

Stationary PDF of the system amplitude

Linearizing the cubic and quintic stiffness terms and taking the undetermined damping and stiffness coefficients as functions of the system amplitude, then the tantamount expression of system (8) can be rewritten as ⁶⁷

$\ddot{x}\left(t\right)+{w_0}^{2}x\left(t\right)+C\left(a\right)\dot{x}\left(t\right)+K\left(a\right)x-u{w}^{p-1}\hspace{0.17em}\mathit{sin}\left(p\pi /2\right)\left(-\epsilon +{\alpha }_1{x}^2\left(t\right)-{\alpha }_2{x}^4\left(t\right)\right)·\dot{x}\left(t\right)=\left[\alpha +\left(1-\alpha \right)x\left(t\right)\right]\zeta \left(t\right)$

(10)

In order to ascertain the coefficients $C\left(a\right)$ and $K\left(a\right)$ in equation (10), the error between system (8) and system (10) is determined by

$e=-u{\alpha }_1{w}^p\hspace{0.17em}\mathit{cos}\left(p\pi /2\right)x^3\left(t\right)+u{\alpha }_2{w}^p\hspace{0.17em}\mathit{cos}\left(p\pi /2\right)x^5\left(t\right)-C\left(a\right)\dot{x}\left(t\right)-K\left(a\right)x\left(t\right)$

(11)

Supposing that the system (10) exists the steady-state solution as the form

$\{\begin{array}{l}x\left(t\right)=a\left(t\right)\mathit{cos}\hspace{0.17em}\phi \left(t\right)\\ \phi \left(t\right)=w_{0}t+\theta \left(t\right)\end{array}$

(12)

where ${w_0}^2=w^2+u\epsilon w^p\hspace{0.17em}\mathit{cos}\left(p\pi /2\right)$ , utilizing the generalized harmonic balance technique and minimizing the error (11) in the mean square sense, and the undetermined coefficients $C\left(a\right)$ and $K\left(a\right)$ can be obtained as follows ⁶⁷

$\{\begin{array}{l}C\left(a\right)=0\\ K\left(a\right)=u{w}^p\hspace{0.17em}\mathit{cos}\left(p\pi /2\right)\frac{\left(5{\alpha }_2{a}^2-6{\alpha }_1\right)a^2}{8}\end{array}$

(13)

Replacing equation (13) into equation (10), the equivalent system can be given by

$\ddot{x}\left(t\right)+{\Omega }^{2}x\left(t\right)-u{w}^{p-1}\hspace{0.17em}\mathit{sin}\left(p\pi /2\right)\left(-\epsilon +{\alpha }_1{x}^2\left(t\right)-{\alpha }_2{x}^4\left(t\right)\right)·\dot{x}\left(t\right)=\left[\alpha +\left(1-\alpha \right)x\left(t\right)\right]\zeta \left(t\right)$

(14)

where ${\Omega }^2=w^2+u\epsilon w^p\hspace{0.17em}\mathit{cos}\left(p\pi /2\right)+u{w}^p\hspace{0.17em}\mathit{cos}\left(p\pi /2\right)\frac{\left(5{\alpha }_2{a}^2-6{\alpha }_1\right)a^2}{8}$ .

Supposing that the system (14) possesses the solution with periodic form, we make the transform as follows ⁶⁸

$\{\begin{array}{l}X=x\left(t\right)=a\left(t\right)\mathit{cos}\hspace{0.17em}\Phi \left(t\right)\\ Y=\dot{x}\left(t\right)=-a\left(t\right)\Omega \hspace{0.17em}\mathit{sin}\hspace{0.17em}\Phi \left(t\right)\\ \Phi \left(t\right)=\Omega t+\theta \left(t\right)\end{array}$

(15)

where $\Omega$ is natural frequency of the equivalent system (14), $a\left(t\right)$ and $\theta \left(t\right)$ represent the amplitude and phase processes of the system’s response, respectively, and they are both random processes.

Substituting equation (15) into equation (14), we attain that

$\{\begin{array}{l}\frac{da}{dt}=F_{11}\left(a,\theta \right)+G_{11}\left(a,\theta \right)\zeta \left(t\right)\\ \frac{d\theta }{dt}=F_{21}\left(a,\theta \right)+G_{21}\left(a,\theta \right)\zeta \left(t\right)\end{array}$

(16)

in which

$\{\begin{array}{l}{F}_{11}\left(a,\theta \right)=u{w}^{p-1}\hspace{0.17em}\mathit{sin}\left(p\pi /2\right)a\hspace{0.17em}{\mathit{sin}}^2\hspace{0.17em}\Phi \left(-\epsilon +{\alpha }_1{a}^2\hspace{0.17em}{\mathit{cos}}^2\hspace{0.17em}\Phi -{\alpha }_2{a}^4\hspace{0.17em}{\mathit{cos}}^4\hspace{0.17em}\Phi \right)\\ F_{21}\left(a,\theta \right)=u{w}^{p-1}\hspace{0.17em}\mathit{sin}\left(p\pi /2\right)\mathit{sin}\hspace{0.17em}\Phi \hspace{0.17em}\mathit{cos}\hspace{0.17em}\Phi \left(-\epsilon +{\alpha }_1{a}^2\hspace{0.17em}{\mathit{cos}}^2\hspace{0.17em}\Phi -{\alpha }_2{a}^4\hspace{0.17em}{\mathit{cos}}^4\hspace{0.17em}\Phi \right)\\ G_{11}=-\alpha \frac{\mathit{sin}\hspace{0.17em}\Phi }{\Omega }-\left(1-\alpha \right)\frac{a\hspace{0.17em}\mathit{cos}\hspace{0.17em}\Phi \hspace{0.17em}\mathit{sin}\hspace{0.17em}\Phi }{\Omega }\\ G_{21}=-\alpha \frac{\mathit{cos}\hspace{0.17em}\Phi }{a\Omega }-\left(1-\alpha \right)\frac{{\mathit{cos}}^2\hspace{0.17em}\Phi }{\Omega }\end{array}$

(17)

Equation (16) can be regarded as the Stratonovich stochastic differential equation, and by adding the Wong–Zakai correction term, the corresponding Itô stochastic differential equation can be acquired as

$\{\begin{array}{l}da=\left[F_{11}\left(a,\theta \right)+F_{12}\left(a,\theta \right)\right]dt+G_{11}\left(a,\theta \right)dB\left(t\right)\\ d\theta =\left[F_{21}\left(a,\theta \right)+F_{22}\left(a,\theta \right)\right]dt+G_{21}\left(a,\theta \right)dB\left(t\right)\end{array}$

(18)

where $B\left(t\right)$ is the standardized Wiener process and

$\{\begin{array}{l}{F}_{12}\left(a,\theta \right)=\frac{1}{2}{{\displaystyle \int }}_{-\infty }^{\infty }\left[\frac{\partial G_{11}\left(a,\theta \right)}{\partial a}{G}_{11}\left(a,\theta ,t+h\right)+\frac{\partial G_{11}\left(a,\theta \right)}{\partial \theta }{G}_{21}\left(a,\theta ,t+h\right)\right]R\left(h\right)dh=\frac{\alpha D\hspace{0.17em}{\mathit{cos}}^2\left(\Phi \right)}{a{\Omega }^2\left(1+{\Omega }^2{\tau }^2\right)}+\left(1-\alpha \right)\frac{a\left({\mathit{sin}}^2\hspace{0.17em}2\Phi +4\hspace{0.17em}\mathit{cos}\hspace{0.17em}2\Phi \hspace{0.17em}{\mathit{cos}}^2\hspace{0.17em}\Phi \right)D}{4{\Omega }^2\left(1+{\Omega }^2{\tau }^2\right)}\\ F_{22}\left(a,\theta \right)=\frac{1}{2}{{\displaystyle \int }}_{-\infty }^{\infty }\left[\frac{\partial G_{21}\left(a,\theta \right)}{\partial a}{G}_{11}\left(a,\theta ,t+h\right)+\frac{\partial G_{21}\left(a,\theta \right)}{\partial \theta }{G}_{21}\left(a,\theta ,t+h\right)\right]R\left(h\right)dh=-\frac{\alpha D\hspace{0.17em}\mathit{sin}\left(2\Phi \right)}{a^2{\Omega }^2\left(1+{\Omega }^2{\tau }^2\right)}-\left(1-\alpha \right)\frac{D\hspace{0.17em}\mathit{sin}\left(2\Phi \right){\mathit{cos}}^2\hspace{0.17em}\Phi }{{\Omega }^2\left(1+{\Omega }^2{\tau }^2\right)}\end{array}$

(19)

Based on the stochastic averaging process and averaging equation (18) over $\Phi$ , ⁶⁹ we derive the following averaged Itô equation

$\{\begin{array}{l}da=m_1\left(a\right)dt+{\sigma }_{11}\left(a\right)dB\left(t\right)\\ d\theta =m_2\left(a\right)dt+{\sigma }_{21}\left(a\right)dB\left(t\right)\end{array}$

(20)

where

$\{\begin{array}{l}{m}_1\left(a\right)=u{w}^{p-1}\hspace{0.17em}\mathit{sin}\left(p\pi /2\right)\left(-\frac{\epsilon }{2}a+\frac{1}{8}{\alpha }_1{a}^3-\frac{1}{16}{\alpha }_2{a}^5\right)+\frac{\alpha D}{2a{\Omega }^2\left(1+{\Omega }^2{\tau }^2\right)}+\left(1-\alpha \right)\frac{3aD}{8{\Omega }^2\left(1+{\Omega }^2{\tau }^2\right)}\\ {{\sigma }_{11}}^2\left(a\right)=\frac{{\alpha }^{2}D}{{\Omega }^2\left(1+{\Omega }^2{\tau }^2\right)}+{\left(1-\alpha \right)}^2\frac{a^{2}D}{4{\Omega }^2\left(1+{\Omega }^2{\tau }^2\right)}\\ m_2\left(a\right)=0\\ {{\sigma }_{21}}^2\left(a\right)=\frac{{\alpha }^{2}D}{a^2{\Omega }^2\left(1+{\Omega }^2{\tau }^2\right)}+{\left(1-\alpha \right)}^2\frac{3D}{4{\Omega }^2\left(1+{\Omega }^2{\tau }^2\right)}\end{array}$

(21)

Equations (20) and (21) show that $da$ doesn’t depend on $\theta$ , and the averaged Itô equation of $a\left(t\right)$ is independent of $\theta \left(t\right)$ . Therefore, the random process $a\left(t\right)$ is a one-dimensional diffusion process.

Hence, the FPK equation of $a\left(t\right)$ can be expressed as its reduced form

$0=-\frac{\partial }{\partial a}\left[m_1\left(a\right)p\left(a\right)\right]+\frac{1}{2}\frac{{\partial }^2}{\partial a^2}\left[{{\sigma }_{11}}^2\left(a\right)p\left(a\right)\right]$

(22)

And the boundary conditions satisfy

$\{\begin{array}{l}p\left(a\right)=c,c\in \left(-\infty ,+\infty \right)asa=0\mathsf{\hspace{0.33em}}\\ p\left(a\right)\to 0,\partial \overline{p}/\partial a\to 0asa\to \infty \end{array}$

(23)

According to the boundary conditions shown in equation (23), the steady-state PDF of system amplitude $a$ can be derived as

$p\left(a\right)=\frac{C}{{{\sigma }_{11}}^2\left(a\right)}\mathit{exp}\left[{{\displaystyle \int }}_0^a\frac{2m_1\left(u\right)}{{{\sigma }_{11}}^2\left(u\right)}du\right]$

(24)

where C is the constant after normalization and obeys

$C={\left[{{\displaystyle \int }}_0^{\infty }\left(\frac{1}{{{\sigma }_{11}}^2\left(a\right)}\mathit{exp}\left[{{\displaystyle \int }}_0^a\frac{2m_1\left(u\right)}{{{\sigma }_{11}}^2\left(u\right)}du\right]\right)da\right]}^{-1}$

(25)

Substituting equation (21) into equation (24), we can get the detailed expression of stationary PDF of the system amplitude $a$ as

$p\left(a\right)=\frac{4C{\mathrm{\Delta }}_1}{D{\mathrm{\Delta }}_2}{a}^{\frac{4\alpha }{{\mathrm{\Delta }}_2}}\hspace{0.17em}\mathit{exp}\left(-\frac{{\mathrm{\Delta }}_1{\mathrm{\Delta }}_3+18D\left(\alpha -1\right)a^2}{12D{\mathrm{\Delta }}_2}\right)$

(26)

where $C$ is the normalized constant and

$\{\begin{array}{l}{\mathrm{\Delta }}_1={\Omega }^2\left({1+\Omega }^2{\tau }^2\right)\\ {\mathrm{\Delta }}_2=\left(a^2+4\right){\alpha }^2-2a^2\alpha +a^2\\ {\mathrm{\Delta }}_3=u{w}^{p-1}\hspace{0.17em}\sin \left(p\pi /2\right)\left(24\epsilon a^2-3{\alpha }_1{a}^4+{\alpha }_2{a}^6\right)\\ {\Omega }^2=w^2+u\epsilon w^p\hspace{0.17em}\mathit{cos}\left(p\pi /2\right)+u{w}^p\hspace{0.17em}\mathit{cos}\left(p\pi /2\right)\frac{\left(5{\alpha }_2{a}^2-6{\alpha }_1\right)a^2}{8}\end{array}$

(27)

Stochastic P-bifurcation of the system

Stochastic P-bifurcation refers to the variety in the quantity of peaks observed in the stationary PDF curve. In this section, we utilize singularity theory to discuss the effects of parameters on the stochastic P-bifurcation behaviors of the system, and further determine the critical parameter conditions.

For the sake of convenience, $p\left(a\right)$ can be expressed by

$p\left(a\right)=\frac{4C}{D}R\left(a,D,\epsilon ,w,p,\tau ,{\alpha }_1,{\alpha }_2\right)\mathit{exp}\left[Q\left(a,D,\epsilon ,w,p,\tau ,{\alpha }_1,{\alpha }_2\right)\right]$

(28)

in which

$\{\begin{array}{l}R(a,D,\epsilon ,w,p,\tau ,{\alpha }_1,{\alpha }_2)=\frac{{\Omega }^2\left(1+{\Omega }^2{\tau }^2\right)}{\left(a^2+4\right){\alpha }^2-2a^2\alpha +a^2}{a}^{\frac{4\alpha }{\left(a^2+4\right){\alpha }^2-2a^2\alpha +a^2}}\\ Q(a,D,\epsilon ,w,p,\tau ,{\alpha }_1,{\alpha }_2)=\frac{-\left[{\Omega }^2\left(1+{\Omega }^2{\tau }^2\right)u{w}^{p-1}\hspace{0.17em}\mathit{sin}\left(p\pi /2\right)·\left(24\epsilon a^2-3{\alpha }_1{a}^4+{\alpha }_2{a}^6\right)\right]}{12D{\mathrm{\Delta }}_2}\\ {\Omega }^2=w^2+u\epsilon w^p\hspace{0.17em}\mathit{cos}\left(p\pi /2\right)+u{w}^p\hspace{0.17em}\mathit{cos}\left(p\pi /2\right)\frac{\left(5{\alpha }_2{a}^2-6{\alpha }_1\right)a^2}{8}\end{array}$

(29)

According to the singularity theory, ⁷⁰ the stationary PDF of system amplitude should satisfy the two conditions as follows

$\{\begin{array}{l}\frac{\partial p\left(a\right)}{\partial a}=0\\ \frac{{\partial }^{2}p\left(a\right)}{\partial a^2}=0\end{array}$

(30)

And substituting equation (28) into equation (30), we acquire

$H=\left\{R^{\prime }+R{Q}^{\prime }=0,R^{″}+2R^{\prime }{Q}^{\prime }+R{Q}^{″}+R{Q}^{\prime 2}=0\right\}$

(31)

where $H$ is the critical condition for the varieties about the number of the PDF curve’s peaks.

In the following section, we investigate the influences of $\mathit{p}$ , $\mathit{\tau }$ , and $\mathit{D}$ on the system. Since the relationship of the three-dimensional surface is not easy to describe and display, here we only give the two-dimensional section of the transition set, and we investigate the influences of the correlated time $\mathit{\tau }$ in the colored noise and the noise intensity D below. Without loss of generality, we take the value of the fractional derivative’s order $\mathit{p}=\mathbf{0.8}$ and the other parameters are taken as ${\alpha }_1=1.51$ , ${\alpha }_2=2.85$ , $w=1$ , $\epsilon =0.1$ , $u=1$ , and $\alpha =1$ , which means that the system is driven by purely additive Gaussian colored noise.

The influence of unfolding parameters $\mathit{\tau }$ and $\mathit{D}$ for stochastic P-bifurcation

Based on equations (29) and (31), we can derive the transition set for stochastic P-bifurcation of the system, and the unfolding parameters $\tau$ and $D$ are shown in Figure 1 as below

Figure 1.

Transition set curves (taking $\tau$ and $D$ as the unfolding parameters).

According to the singularity theory, the topological structure of the stationary PDF curves remains qualitatively consistent among various points in the identical region. By selecting one point ( $\tau ,D$ ) in each region, we can acquire all the different types of the stationary PDF curves which are qualitatively diverse, and the unfolding parametric plane about $\tau$ and $D$ is divided into three sub-areas by the transition set curve. To ensure convenience, every section in Figure 1 is assigned to one number.

Initially, we investigate the stationary PDF of amplitude $p\left(a\right)$ for one point ( $\tau ,D$ ) in each of the three sub-regions depicted in Figure 1. Subsequently, we contrast the theoretical solutions with the numerical results acquired through Monte Carlo simulation ⁷¹ of the original fractional-order system (3), and the corresponding outcomes are displayed in Figure 2.

Figure 2.

PDF $p\left(a\right)$ of the system in three sub-regions of Figure 1(taking $\tau$ and $D$ as the unfolding parameters). (a) Parameter ( $\tau ,D$ ) in area 1 of Figure 1. (b) Parameter ( $\tau ,D$ ) in area 2 of Figure 1. (c) Parameter ( $\tau ,D$ ) in area 3 of Figure 1.

According to Figure 1, the unfolding parametric area where the PDF arises bimodal is surrounded by two curves. When the unfolding parameters ( $\tau ,D$ ) is considered as $\tau =0.6,D=0.001$ in area 1, the peak of the PDF $p\left(a\right)$ is near the origin and the system has a stable equilibrium as shown in Figure 2(a); when the unfolding parameters ( $\tau ,D$ ) is considered as $\tau =0.5,D=0.005$ in area 2, the PDF $p\left(a\right)$ of the system has two peaks, and there exists a stable limit cycle whose corresponding amplitude $a$ is distant from the origin and the probability near the origin is also non-zero, consequently, the limit cycle and equilibrium coexist in the system simultaneously, as shown in Figure 2(b); when the unfolding parameters ( $\tau ,D$ ) is considered as $\tau =0.2,D=0.009$ in area 3, the PDF $p\left(a\right)$ only has one peak which is far away from the origin, and the system presented as a stable limit cycle at this time, as shown in Figure 2(c).

The influence of unfolding parameters $\mathit{p}$ and $\mathit{D}$ for stochastic P-bifurcation

Without losing generality, we take the value of correlation time $\mathit{\tau }=\mathbf{0.5}$ in the colored noise and the other parameters are taken as ${\alpha }_1=1.51$ , ${\alpha }_2=2.85$ , $w=1$ , $\epsilon =0.1$ , $u=0.5$ , and $\alpha =1$ . And according to equations (29) and (31), we can also obtain the transition set for stochastic P-bifurcation of the system, and the unfolding parameters $p$ and $D$ are shown in Figure 3 as below.

Figure 3.

Transition set curves (taking $p$ and $D$ as the unfolding parameters).

We analyze the stationary PDF $p\left(a\right)$ of system amplitude for one point ( $p,D$ ) in each of the three sub-regions displayed in Figure 3. And then, we compare the theoretical solutions with the numerical results obtained through Monte Carlo simulation ³⁸ of the original system (3), and the corresponding outcomes are derived in Figure 4.

Figure 4.

PDF $p\left(a\right)$ of the system in three sub-regions of Figure 3(taking $p$ and $D$ as the unfolding parameters). (a) Parameter ( $p,D$ ) in area 1 of Figure 3 (b) Parameter ( $p,D$ ) in area 2 of Figure 3. (c) Parameter ( $p,D$ ) in area 3 of Figure 3.

From Figure 3, we can see that the unfolding parametric area where the PDF occurs bimodal is also surrounded by the two transition set curves. And when the unfolding parameter ( $p,D$ ) is taken as $p=0.3,D=0.0004$ in area 1, the PDF $p\left(a\right)$ has a stable equilibrium as shown in Figure 4(a); when the unfolding parameter ( $p,D$ ) is taken as $p=0.4,D=0.002$ in area 2, the PDF $p\left(a\right)$ has a stable limit cycle far away from the origin and the probability is not zero near the origin, there are both the limit cycle and equilibrium in the system simultaneously as shown in Figure 4(b); when the unfolding parameter ( $p,D$ ) is taken as $p=0.8,D=0.006$ in area 3, the PDF $p\left(a\right)$ appears in the form of a stable limit cycle at the moment as shown in Figure 4(c).

In brief, the stationary PDF in any two adjacent areas in Figures 1 and 3 exhibits distinct qualitative property. Regardless of the values of the unfolding parameters, the system will undergo stochastic P-bifurcation behavior as long as the unfolding parameters pass through any line in Figures 1 and 3. Hence, the transition set curves are just right the critical parametric conditions of the stochastic P-bifurcation about the system. The analytic solutions presented in Figures 2 and 4 are well in accordance with the results obtained from Monte Carlo simulation of the original system (3), which provides further validation to the theoretical analysis and demonstrates the applicability of the methods presented in this paper for analyzing the stochastic P-bifurcation behavior of fractional-order systems.

In comparison to the controllers with integer order, the fractional-order controllers exhibit superior dynamic performances and robustness. ⁷¹ Over the last few years, a considerable number of controllers with fractional-order derivative have been proposed.^72–75 Through the analysis mentioned above, we have identified the critical parameters that the stochastic P-bifurcation will take place in system (3), which can enable the system to convert between monostable and bistable states by choosing the unfolding parameters appropriately. And this should provide the theoretical guidance for analyzing and designing the fractional-order controllers.

Conclusion

This paper investigated the stochastic P-bifurcation of a modified fractional-order Van der Pol system driven by Gaussian colored noise excitation, and the following conclusions can be drawn:

(1) The results indicate that each of the fractional derivative’s order $p$ , the correlation time $\tau$ and the noise intensity $D$ can induce stochastic P-bifurcation of the system, and the number of peaks in the stationary PDF curves of the system amplitude can be switched from one to two by adjusting the unfolding parameters.

(2) The parametric region where the system presents bistable behaviors is surrounded by the two transition set curves, and the transition set curves obtained are just right the critical parametric conditions of the stochastic P-bifurcation about the system.

(3) The system’s response can be retained at the monostability or a small vibration amplitude near the equilibrium by selecting the appropriate unfolding parameters, which can provide the theoretical guidance in the design of such systems and prevent damage and instability caused by nonlinear jumps or large amplitude vibrations of the system.

(4) It also demonstrates that the theoretical approach utilized in this study is viable to explore the stochastic P-bifurcation phenomena of nonlinear oscillators with fractional-order derivative element.

However, the system studied in the article is the single degree of freedom system, and the complexity as well as the abstraction of state space increase the difficulty to analyze the high-dimensional dynamic system. The investigation of two degree of freedom systems or even higher dimensional and coupled systems should be the next research focus in future.