Stochastic P-bifurcation in a generalized Van der Pol oscillator with fractional delayed feedback excited by combined Gaussian white noise excitations

Abstract

This paper investigates the stochastic P-bifurcation behavior of bistability in a generalized Van der Pol system with fractional time-delay feedback under additive and multiplicative Gaussian white noise excitations. First, using the minimal mean square error principle, the fractional derivative is found to be equivalent to a linear combination of damping and restoring forces, and the original system is transformed into an equivalent integer-order system. Second, the stationary probability density function of system amplitude is obtained by stochastic averaging, and according to singularity theory, the critical parameters for stochastic P-bifurcation of the system are found. Finally, the different types of stationary probability density function curves of the system amplitude are qualitatively analyzed by choosing the corresponding parameters in each region divided by the transition set curves. The consistency between the analytical solutions and the Monte Carlo simulation results verifies the theoretical analysis in this paper. The method used in this paper can directly guide the design of the fractional-order controller to adjust the response of the system.

Introduction

Fractional calculus is a generalization of integer-order calculus, which has a history of more than 300 years. Due to the limitation of the definition of integer-order derivative, the classical integer operators cannot express memory properties and do not have sufficient parameters to handle the different shapes of the hysteresis loops describing the behaviors of viscoelastic materials and structures. While fractional derivatives contain convolution, which can describe a memory effect and express a cumulative effect over time, hence they are more suitable to describe memory characteristics^1,2 and have become a powerful mathematical tool to study fields such as anomalous diffusion, non-Newtonian fluid mechanics, viscoelastic mechanics, and soft matter physics. Compared with integer-order calculus, the fractional derivative can describe various reaction processes more accurately with fewer parameters,^3,4 thus, it is necessary and significant to investigate the fractional differential equations on the typical mechanical properties and the influences of fractional-order parameters on the system.

Recently, many scholars have studied the dynamic behavior of nonlinear multi-stable systems under different noise excitations and achieved fruitful results. Yildirim and Gülkanat ⁵ obtained an approximate solution of a fractional Zakharov–Kuznetsov equation by using the homotopy perturbation method and revealed the effectiveness of the method through several examples. He and Ji ⁶ suggested the two-scale fractional calculus, for a large scale we have, for example, continuum mechanics, while on a small scale, saying a molecular scale, water becomes discontinuous, which can be considered as a fractal medium. Wang et al. ⁷ established the governing equations for the nonlinear transverse vibration of an axially moving viscoelastic beam with finite deformation using the Hamiltonian principle and produced nanoscale crimped fibers using stuffer box crimping and bubble electrospinning. He and Ji ⁸ and He ⁹ proposed a simple approach to the nonlinear oscillators and obtained the Taylor series solution for the well-known Lane–Emden equation. Li et al. ¹⁰ solved a paradox in an electrochemical sensor by a fractal modification of the surface coverage model and elucidated a simple solution process to the fractal model. He et al. ¹¹ pointed out the so-called enhanced variational iteration method for a nonlinear equation arising in electrospinning and vibration-electrospinning process is the standard variational iteration method, and an effective algorithm using the variational iteration algorithm-II is suggested for Bratu-like equation arising in electrospinning. Anjum and He ¹² suggested an easier approach by the Laplace transform to determining the Lagrange multiplier, making the method accessible to researchers facing various nonlinear problems, and adopted a nonlinear oscillator as an example to elucidate the identification process and the solution process. Roul ¹³ presented the analytical and numerical solutions of a degenerate parabolic equation with time-fractional derivatives arising in the spatial diffusion of biological populations by using the homotopy-perturbation method. Wu and He ¹⁴ elucidated that the homotopy perturbation method is valid for nonlinear oscillators with negative linear terms, and conditions for the periodic solutions can be easily obtained. Cai et al. ¹⁵ applied the fractal derivative to modeling viscoelastic behavior and obtained the creep modulus and relaxation compliance for the proposed fractal Maxwell and Kelvin models by utilizing the methodology of scaling transformation. Wang and Wang ¹⁶ modified the reduced differential transform method for obtaining the approximate analytical solutions of the fractional heat transfer equations. Wang and Liu ¹⁷ established a modified reduced differential transform method and a new iterative Elzaki transform method, and then applied them to obtain the analytical solutions of the time-fractional Navier–Stokes equations. In addition, the authors investigated the Van der Pol–Duffing oscillators under the colored noise, combined harmonic and random excitations, respectively; moreover, the stochastic P-bifurcation behaviors of the noise oscillators were discussed by analyzing changes in the stationary probability density function (PDF) of the systems.^18,19 Wu and Hao ²⁰ investigated the stochastic P-bifurcation of tri-stability in a generalized Duffing–Van der Pol oscillator system excited by multiplicative colored noise, obtained an analytical expression of the stationary PDF of system amplitude, and analyzed the influences of noise intensity and system parameters on stochastic P-bifurcation of the system. Chen and Zhu ²¹ studied the response of a Duffing system with fractional damping under combined white noise and harmonic excitations, and showed that variation in the order of fractional derivative can cause stochastic P-bifurcation of the system. Li et al. ²² investigated the stochastic P-bifurcation behavior of a bistable Van der Pol–Duffing oscillator with fractional derivative excited by additive and multiplicative Gaussian colored noise excitations and found that changes in the linear damping coefficient, the order of fractional derivative, and the noise intensity can each lead to stochastic P-bifurcation of the system. Liu et al. ²³ investigated the stochastic stability of a Duffing oscillator with fractional derivative damping under combined harmonic and Poisson white noise parametric excitations, and analyzed the asymptotic Lyapunov stability with probability one of the original system by using the largest Lyapunov exponent. Shen et al. ²⁴ investigated the primary resonance of Duffing oscillator with two kinds of fractional-order derivatives analytically and analyzed the effects of two kinds of fractional-order derivatives on dynamic behaviors of the system. Yao and Wang ²⁵ put forward a new approximate analytical method for a system of fractional differential equations and gave some examples to verify the correctness of the method proposed.

For the dynamics of time-delay systems, Chen et al. ²⁶ proposed a stochastic averaging technique, which can be used to study randomly excited strongly nonlinear systems with a delayed feedback fractional-order proportional-derivative controller, and obtained stationary PDF of the system. Wen et al. ²⁷ studied the deterministic and autonomous Duffing systems with fractional time-delay coupled feedback, and found that fractional time-delay coupled feedback plays the roles of both velocity time-delay feedback and displacement time-delay feedback. Jiang and Wei ²⁸ considered a classical Van der Pol oscillator with general time-delay feedback and found that there are the Bogdanov–Takens bifurcation, triple-zero, and Hopf-zero singularities in the system by analyzing the distribution of the associated characteristic roots.

Because of the complexity of fractional derivatives, analyzing them is difficult, and the influences of system parameters on vibration characteristics are mostly studied numerically, which are usually limited to the qualitative analysis. It is difficult to find the critical condition of parametric influence, which affects the analysis and design of such systems, in part because the bistable stochastic P-bifurcation of the fractional delayed feedback system has not been reported. Accordingly, we take the nonlinear vibration of a generalized Van der Pol oscillator excited by both additive and multiplicative Gaussian white noise excitations simultaneously as an example, and obtain the critical parametric conditions for stochastic P-bifurcation using the singularity method. Furthermore, we compare the Monte Carlo simulation results with the analytical solutions obtained by stochastic averaging, their consistency verifies the theoretical analysis in this paper.

Derivation of the equivalent system

There are many definitions of fractional derivatives and the following definitions are mainly introduced.

The Caputo derivative of the function $x(t)$ defined on the interval $\left[a,b\right]$ is formulated as

$${}_a{}^{C}D{\hspace{0.17em}}^p[x(t)]=\frac{1}{\Gamma (m-p)}{\displaystyle {\int }_a^t\frac{x^{(m)}(u)}{{(t-u)}^{1+p-m}}}\text{d}u$$

(1)

where $p$ represents the order of fractional derivative ${}_a{}^{C}D{\hspace{0.17em}}^p[x(t)]$ , $m-1<p\le m,m\in N$ , $t\in [a,b]$ , $\Gamma (m)$ is the Euler Gamma function, and $x^{(m)}(t)$ is the m order derivative of $x(t)$ .

The Riemann–Liouville derivative of the function $x(t)$ defined on the interval $\left[a,b\right]$ is formulated as

$${}_{a}D{\hspace{0.17em}}^p[x(t)]=\frac{1}{\Gamma (m-p)}\frac{d^m}{d{t}^m}{\displaystyle {\int }_a^t\frac{x(u)}{{(t-u)}^{1+p-m}}}\text{d}u$$

(2)

where $p$ represents the order of fractional derivative ${}_{a}D{\hspace{0.17em}}^p[x(t)]$ , $m-1<p\le m,m\in N$ , $t\in [a,b]$ and $\Gamma (m)$ is the Euler Gamma function.

The fractal derivative of the function $x(t)$ is formulated as^29–31

$$\frac{\mathit{dx}\text{(}t\text{)}}{d{t}^p}=\underset{t^{\prime }-t}{\text{lim}}\frac{x(t)-x(t^{\prime })}{t^p-t^{\prime }{\hspace{0.17em}}^p}$$

(3)

where $p$ represents the fractal in time and the order of the fractal derivative.

The Riemann–Liouville derivative and Caputo derivative are most commonly used, the initial condition corresponding to Riemann–Liouville derivative has no physical meaning; however, the initial condition of the system described by Caputo derivative has both clear physical meaning and forms the same as in the integer-order differential equation. So in this paper, we adopt the Caputo derivative as defined in equation (1).

For a given physical system, because the moment when the oscillator begins to vibrate is always $t=0$ , and the Caputo derivative is often used in the following form

$${}_0{}^{C}D{\hspace{0.17em}}^p[x(t)]=\frac{1}{\Gamma (m-p)}{\displaystyle {\int }_0^t\frac{x^{(m)}(u)}{{(t-u)}^{1+p-m}}}\text{d}u$$

(4)

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

In this paper, we study the generalized Van der Pol oscillator system with fractional-order time-delay coupled feedback excited by Gaussian white noise excitations

$$\ddot{X}(t)-(-\epsilon +{\alpha }_1{X}^2(t)-{\alpha }_2{X}^4(t)+{\alpha }_3{X}^6(t))\dot{X}(t)+w^{2}X(t)+{}_0{}^{C}D{\hspace{0.17em}}^p[X(t-\tau )]={\xi }_1(t)+X(t){\xi }_2(t)$$

(5)

where $\epsilon$ represents the linear damping coefficient; ${\alpha }_1$ , ${\alpha }_2$ , and ${\alpha }_3$ represent the nonlinear damping coefficients of the system; $w$ is the natural frequency; and $\tau$ is the time-delay introduced in the system. ${}_0{}^{C}D{\hspace{0.17em}}^p[X(t-\tau )]$ is the $p$ ( $0\le p\le 1$ ) order Caputo derivative of $X(t-\tau )$ with respect to $t$ , which is defined by equation (4). ${\xi }_k(t)(k=1,2)$ are two independent Gaussian white noises, which satisfy

$$E[{\xi }_k(t)]=0,E[{\xi }_k(t){\xi }_k(t-s)]=2D_k\delta (s)(k=1,2)$$

(6)

where $D_k$ $(k=1,2)$ denote the intensities of Gaussian white noises ${\xi }_k(t)$ , respectively, and $\delta (s)$ is the Dirac function.

The fractional derivative has the contributions of damping force and restoring force,^32–34 hence, we introduce the equivalent system as follows

$$\ddot{X}(t)-(-\epsilon +{\alpha }_1{X}^2(t)-{\alpha }_2{X}^4(t)+{\alpha }_3{X}^6(t)+C(p,\tau ))\dot{X}(t)+(K(p,\tau )+w^2)X(t)={\xi }_1(t)+X(t){\xi }_2(t)$$

(7)

where $C(p,\tau )$ and $K(p,\tau \text{)}$ are coefficients of the equivalent damping and restoring forces of fractional derivative ${}_0{}^{C}D{\hspace{0.17em}}^p[X(t-\tau )]$ , respectively.

Applying the equivalent methods mentioned in Li et al. ²² and Yang et al. ³⁴ the concrete forms of $C(p,\tau )$ and $K(p,\tau )$ are as follows

$$C(p,\tau )=-w^{p-1}\sin (p\pi /2-w\tau ),K(p,\tau )=w^p\text{cos}(p\pi /2-w\tau )$$

(8)

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

$$\ddot{X}(t)-\gamma \dot{X}(t)+w_0^{2}X(t)={\xi }_1(t)+X(t){\xi }_2(t)$$

(9)

where

$$\{\begin{array}{l}\gamma =-\epsilon +{\alpha }_1{X}^2-{\alpha }_2{X}^4+{\alpha }_3{X}^6-w^{p-1}\sin (p\pi /2-w\tau )\\ w_0^2=w^2+w^p\cos (p\pi /2-w\tau )\end{array}$$

(10)

The stationary PDF of system amplitude

In our first example, we examine the system in equation (9), with linear and nonlinear damping coefficients $\epsilon =0.2$ , ${\alpha }_1=1.51$ , ${\alpha }_2=2.85$ , ${\alpha }_3=-1.693$ , natural frequency $w=1$ , and time-delay $\tau =0.5$ . For convenience in discussing parametric influence, the bifurcation diagram of system amplitude with the variation of fractional order $p$ is shown in Figure 1 when $D_1=D_2\text{=}0$ .

Figure 1.

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

As can be seen from Figure 1, the Hopf bifurcation point is marked by $p_H=0.189$ , and the Fold bifurcation point is marked by $p_F=0.245$ ; there are two attractors when $p$ changes in [0.189, 0.245]: equilibrium and limit cycle, and the corresponding result is shown in Figure 2.

Figure 2.

Phase diagram of the deterministic system (at $p$ =0.22). Equilibrium and limit cycle.

For a system with linear and nonlinear damping coefficients $\epsilon =-0.5$ , ${\alpha }_1=1.51$ , ${\alpha }_2=2.85$ , ${\alpha }_3=-1.693$ , natural frequency $w=1$ , and the fractional order $p=0.5$ , the bifurcation diagram of system amplitude with the variation of time-delay $\tau$ is shown in Figure 3 when $D_1=D_2\text{=}0$ .

Figure 3.

Bifurcation diagram of the deterministic system (with variation in $\tau$ ).

As can be seen from Figure 3, the Fold bifurcation point is marked by ${\tau }_F=0.159$ , and the Hopf bifurcation point is marked by ${\tau }_H=0.262$ ; it also shows that there are two attractors where $0.159\le p\le 0.262$ : limit cycle and equilibrium, and the corresponding result is shown in Figure 4.

Figure 4.

Phase diagram of the deterministic system (at $\tau$ =0.2). Equilibrium and limit cycle.

Assuming that the solution of system (9) has the periodic form, and we introduce the following transformation ³⁵

$$\{\begin{array}{l}X=A(t)\cos \Phi (t)\\ \dot{X}=-A(t)w_0\sin \Phi (t)\\ \Phi (t)=w_{0}t+\theta (t)\end{array}$$

(11)

where $w_0$ is natural frequency of the equivalent system (9), $A(t)$ and $\theta (t)$ represent the amplitude and phase processes of system response, respectively, and they are both random processes.

Substituting equation (11) into equation (9), we can obtain

$$\{\begin{array}{c}\frac{\mathit{dA}}{\mathit{dt}}=F_{11}(A,\theta )+G_{11}(A,\theta ){\xi }_1(t)+G_{12}(A,\theta ){\xi }_2(t)\\ \frac{d\theta }{\mathit{dt}}=F_{21}(A,\theta )+G_{21}(A,\theta ){\xi }_1(t)+G_{22}(A,\theta ){\xi }_2(t)\end{array}$$

(12)

in which

$$\{\begin{array}{l}{F}_{11}(A,\theta )=A{\sin }^2\Phi (-\epsilon +{\alpha }_1{A}^2{\cos }^2\Phi -{\alpha }_2{A}^4{\cos }^4\Phi +{\alpha }_3{A}^6{\cos }^6\Phi -w^{p-1}\sin \left(\frac{p\pi }{2}-w\tau \right))\\ F_{21}(A,\theta )=\sin \Phi \cos \Phi (-\epsilon +{\alpha }_1{A}^2{\cos }^2\Phi -{\alpha }_2{A}^4{\cos }^4\Phi +{\alpha }_3{A}^6{\cos }^6\Phi -w^{p-1}\sin \left(\frac{p\pi }{2}-w\tau \right))\\ G_{11}=-\frac{\sin \Phi }{w_0},G_{12}=-\frac{A\sin \Phi \cos \Phi }{w_0}\\ G_{21}=-\frac{\cos \Phi }{A{w}_0},G_{22}=-\frac{{\cos }^2\Phi }{w_0}\end{array}$$

(13)

Equation (12) can be treated as the Stratonovich differential equation, and by adding the relevant correction term, we can transform it into the corresponding Itô stochastic differential equation as follows

$$\{\begin{array}{c}\mathit{dA}=\left[F_{11}(A,\theta )+F_{12}(A,\theta )\right]dt+\sqrt{2D_k}{G}_{1k}(A,\theta )d{B}_k(t)\\ \begin{array}{l}d\theta =\left[F_{21}(A,\theta )+F_{22}(A,\theta )\right]dt+\sqrt{2D_k}{G}_{2k}(A,\theta )d{B}_k(t)\\ k=1,2\end{array}\end{array}$$

(14)

where $B_k(t)$ are independent and normalized Wiener processes, in addition

$$\{\begin{array}{c}{F}_{12}(A,\theta )=D_1\frac{\partial G_{11}}{\partial A}{G}_{11}+D_1\frac{\partial G_{11}}{\partial \theta }{G}_{21}+D_2\frac{\partial G_{12}}{\partial A}{G}_{12}+D_2\frac{\partial G_{12}}{\partial \theta }{G}_{22}\\ F_{22}(A,\theta )=D_1\frac{\partial G_{21}}{\partial A}{G}_{11}+D_1\frac{\partial G_{21}}{\partial \theta }{G}_{21}+D_2\frac{\partial G_{22}}{\partial A}{G}_{12}+D_2\frac{\partial G_{22}}{\partial \theta }{G}_{22}\end{array}$$

(15)

By stochastic averaging of averaging equation (14) over $\Phi$ , ³⁶ we can obtain the following averaged Itô equations

$$\{\begin{array}{c}\mathit{dA}=m_1(A)\mathit{dt}+{\sigma }_{11}(A)d{B}_1(t)\text{+}{\sigma }_{12}(A)d{B}_2(t)\\ d\theta =m_2(A)\mathit{dt}+{\sigma }_{21}(A)d{B}_1(t)\text{+}{\sigma }_{22}(A)d{B}_2(t)\end{array}$$

(16)

where $B_1(t)$ and $B_2(t)$ are two unit Wiener processes that are independent of each other and

$$\{\begin{array}{l}{m}_1(A)=-\frac{1}{2}(w^{p-1}\sin (p\pi /2-w\tau )+\epsilon )A+\frac{1}{8}{\alpha }_1{A}^3-\frac{1}{16}{\alpha }_2{A}^5+\frac{5}{128}{\alpha }_3{A}^7+\frac{D_1}{2A{w}_0^2}+\frac{3D_{2}A}{8w_0^2}\\ {\sigma }_{11}{\hspace{0.17em}}^2(A)=\frac{D_1}{w_0^2},{\sigma }_{12}{\hspace{0.17em}}^2(A)\text{=}\frac{D_2{A}^2}{4w_0^2}\\ m_2(A)=0\\ {\sigma }_{21}{\hspace{0.17em}}^2(A)=\frac{D_1}{A^2{w}_0^2},{\sigma }_{22}{\hspace{0.17em}}^2(A)\text{=}\frac{3D_2}{4w_0^2}\end{array}$$

(17)

Equations (16) and (17) show that the averaged Itô equation of $A(t)$ is independent of $\theta (t)$ , and the process $A(t)$ is actually a one-dimensional diffusion process. Then the reduced Fokker–Planck–Kolmogorov equation of $A(t)$ can be written as follows

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

(18)

The boundary conditions satisfy

$$\{\begin{array}{l}p(A)=c,c\in (-\infty ,+\infty )\mathit{as}A=0\hfill \\ p(A)\to 0,\partial p(A)/\partial A\to 0\mathit{as}A\to \infty \hfill \end{array}$$

(19)

Based on the boundary conditions (19), the stationary PDF of system amplitude $A$ can be obtained as

$$p(A)=\frac{C}{{\sigma }_{11}{\hspace{0.17em}}^2(A)+{\sigma }_{12}{\hspace{0.17em}}^2(A)}\exp \left[{\displaystyle {\int }_0^A\frac{2m_1(u)}{({\sigma }_{11}{\hspace{0.17em}}^2(u)+{\sigma }_{12}{\hspace{0.17em}}^2(u))}\text{d}u}\right]$$

(20)

where C is the normalization constant that satisfies

$$C={\left[{\displaystyle {\int }_0^{\infty }\left(\frac{1}{{\sigma }_{11}{\hspace{0.17em}}^2(A)+{\sigma }_{12}{\hspace{0.17em}}^2(A)}\exp \left[{\displaystyle {\int }_0^A\frac{2m_1(u)}{({\sigma }_{11}{\hspace{0.17em}}^2(u)+{\sigma }_{12}{\hspace{0.17em}}^2(u))}\text{d}u}\right]\right)}\text{d}A\right]}^{-1}$$

(21)

Substituting equation (17) into equation (20), the explicit expression of stationary PDF of system amplitude $A$ can be described as

$$p(A)=4CA{w}_0^2{(4D_1+A^2{D}_2)}^{-\frac{{\Delta }_1}{D^4{\hspace{0.17em}}_2}}\exp \left(\frac{{\Delta }_2}{96D_2^3}\right)$$

(22)

where

$$\{\begin{array}{l}{\Delta }_1=2w_0^2[(\epsilon +\sin (p\pi /2-w\tau )w^{p-1})D_2^3+{\alpha }_1{D}_1{D}_1^2+2{\alpha }_2{D}_1^2{D}_2+5{\alpha }_3{D}_1^3]\\ {\Delta }_2=A^2{w}_0^2(48{\alpha }_1{D}_2^2+96{\alpha }_2{D}_1{D}_2+240{\alpha }_3{D}_1^2-12{\alpha }_2{D}_2^2{A}^2-30{\alpha }_3{D}_1{D}_2{A}^2+5{\alpha }_3{D}_2^2{A}^4)\end{array}$$

(23)

Stochastic P-bifurcation of system amplitude

Stochastic P-bifurcation means the changes in the number of peaks in the PDF curve. To obtain the critical parametric conditions for stochastic P-bifurcation, we analyze the parametric influences on stochastic P-bifurcation of the system by using singularity theory in this section.

For convenience, $p(A)$ can be expressed as follows

$$p(A)=4CR(A,D,\epsilon ,w,p,{\alpha }_1,{\alpha }_2,{\alpha }_3)\mathit{exp}[Q(A,D,\epsilon ,w,p,{\alpha }_1,{\alpha }_2,{\alpha }_3)]$$

(24)

where

$$\{\begin{array}{l}\begin{array}{l}R(A,D_1,D_2,\epsilon ,w,p,{\alpha }_1,{\alpha }_2,{\alpha }_3)=A{w}^2{\hspace{0.17em}}_0{(4D_1+D_2{A}^2)}^{-\frac{2w_0^2}{D_2^4}}[(\epsilon +\sin (p\pi /2-w\tau )w^{p-1})D_2^3\hfill \\ +{\alpha }_1{D}_1{D}_1^2+2{\alpha }_2{D}_1^2{D}_2+5{\alpha }_3{D}_1^3]\hfill \end{array}\hfill \\ \begin{array}{l}Q(A,D_1,D_2,\epsilon ,w,p,{\alpha }_1,{\alpha }_2,{\alpha }_3)=\frac{A^2{w}_0^2}{96D_2^3}(48{\alpha }_1{D}_2^2+96{\alpha }_2{D}_1{D}_2+240{\alpha }_3{D}_1^2-12{\alpha }_2{D}_2^2{A}^2\hfill \\ -30{\alpha }_3{D}_1{D}_2{A}^2+5{\alpha }_3{D}_2^2{A}^4)\hfill \end{array}\end{array}$$

(25)

Based on singularity theory, ³⁷ the stationary PDF of system amplitude needs to meet the following two conditions

$$\frac{\partial p(A)}{\partial A}=0,\frac{{\partial }^{2}p(A)}{\partial A^2}=0$$

(26)

Substituting equation (24) into equation (26), we can obtain the following condition^20,22

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

(27)

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

The influence of the fractional order $p$ 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 to represent the influences of the fractional order $p$ and the noise intensities $D_1$ and $D{}_2\hspace{0.17em}$ below.

According to the distribution of deterministic attractors in Figure 1, without loss of generality, we take the different values of the fractional order $p$ (other parameters are the same as those in Figure 1) in the monostable and bistable intervals separately and calculate the corresponding transition sets according to equations (25) and (27). As the transition set is the empty set when $p$ is taken in the interval [0, 0.189], we show the transition sets only when $p$ is taken as (a) $p=0.24$ , in the interval [0.190, 0.245]; and (b) $p=0.27$ , in the interval [0.246, 0.3], and as shown in Figure 5(a) and (b), respectively.

Figure 5.

Transition sets under the different values of fractional order $p$ (taking $D_2$ and $D_1$ as unfolding parameters). (a) The fractional order $p=0.24$ and (b) the fractional order $p=0.27$ .

From the above analysis and the deterministic bifurcation diagram shown in Figure 1, it can be concluded that when the fractional order $p$ is taken in the interval [0, 0.189], the deterministic system corresponding to system (5) has only one attractor: limit cycle, and the transition set of the system under additive and multiplicative noise excitations together is the empty set, which implies that it cannot make the system appear bistable by adding the noise excitations. When $p$ is taken in the bistable interval [0.190,0.245] of the deterministic system (5), without losing generality, taking $p=0.24$ as an example, the transition set curves of the system under two noise excitations together are shown in Figure 5(a). Because the deterministic system is bistable at this time, the small noise intensities $D_1$ and $D_2$ can arouse the bistable characteristic of the stochastic system, and with the increasing of $D_1$ and $D_2$ , the system becomes monostable again. When $p$ is taken in the interval [0.246, 0.3], taking $p=0.27$ as an example, the deterministic system corresponding to system (5) also has an attractor: equilibrium, and the transition set curves under two noise excitations together are shown in Figure 5(b). It is noteworthy that the system response exhibits a bistable region, and the bimodal region of the stationary PDF curve of system amplitude is surrounded by an approximate triangle. At this moment, although there is only one equilibrium in the deterministic system, the appropriate sizes for $D_1$ and $D_2$ can excite the bistability of the stochastic system.

Based on singularity theory, types of the stationary PDF curves of system amplitude at different points ( $D_2,D_1$ ) in the same region are qualitatively identical. By taking one point ( $D_2,D_1$ ) in each region, we can obtain all varieties of the stationary PDF curves that are qualitatively different. The unfolding parametric plane $D_2-D_1$ is divided into two sub-regions by the transition set curve; for the sake of convenience, each region in Figure 5 is marked with a number.

Without loss of generality, we analyze the stationary PDF of amplitude $p\text{(}A\text{)}$ and the joint PDF $p\text{(}X\text{,}\dot{X}\text{)}$ only for one point ( $D_2,D_1$ ) in each of the two sub-regions in Figure 5(b), and then compare the analytical solutions with the Monte Carlo simulation results from original system (5) using the numerical method for fractional derivative, ³⁸ and show the corresponding results in Figures 6 and 7.

Figure 6.

PDFs for $p\text{(}A\text{)}$ in different sub-regions in Figure 5(b) (taking $D_2$ and $D_1$ as unfolding parameters). (a) Parameter ( $D_2,D_1$ ) in region 1 in Figure 5(b) and (b) parameter ( $D_2,D_1$ ) in region 2 in Figure 5(b).

Figure 7.

Joint PDFs for $p\text{(}X\text{,}\dot{X}\text{)}$ in different sub-regions in Figure 5(b) (taking $D_2$ and $D_1$ as unfolding parameters). (a) Parameter ( $D_2,D_1$ ) in region 1 in Figure 5(b) and (b) parameter ( $D_2,D_1$ ) in region 2 in Figure 5(b).

As shown in Figure 5(b), the parametric region where the PDF occurs bimodal is surrounded by an approximately triangular region. And when we take the parameter ( $D_2,D_1$ ) in region 1, the PDFs for $p\text{(}A\text{)}$ and $p\text{(}X\text{,}\dot{X}\text{)}$ each have a stable equilibrium as shown in Figures 6(a) and 7(a); in region 2, the PDFs for $p\text{(}A\text{)}$ and $p\text{(}X\text{,}\dot{X}\text{)}$ each have a stable limit cycle, and the probability near the origin is not zero at the moment, there are both the equilibrium and limit cycle in the system simultaneously, as shown in Figures 6(b) and 7(b).

The influence with variation of time-delay $\tau$ on the system

Similar to the above, according to the distributions of deterministic attractors in Figure 3, without loss of generality, we take the different values of time-delay $\tau$ (other parameters are the same as those in Figure 3) in the monostable and bistable intervals separately and calculate the corresponding transition sets according to equations (25) and (27), as the transition set is the empty set when $\tau$ is taken in the interval [0.263, 0.5], we show the transition sets only when the time-delay $\tau$ is taken as (a) $\tau =0.1$ in the interval [0, 0.159]; (b) $\tau =0.17$ in the interval [0.160, 0.262] respectively, and as shown in Figure 8.

Figure 8.

Transition sets under the different values of time-delay $\tau$ (taking $D_2$ and $D_1$ as unfolding parameters). (a) The time-delay $\tau =0.1$ and (b) the time-delay $\tau =0.17$ .

From the deterministic bifurcation diagram shown in Figure 3 and the transition sets shown in Figure 8, it can be seen that when the time-delay $\tau$ is taken in the interval [0, 0.159], without losing generality, taking $\tau =0.1$ as an example, the deterministic system corresponding to system (5) has only one attractor: equilibrium, and the transition set curves under two noise excitations together are shown in Figure 8(a). It is also noteworthy that the system exhibits a bistable region, and the bimodal region of the stationary PDF of system amplitude is also surrounded by an approximate triangle. At this moment, although there is only one equilibrium in the deterministic system, the appropriate sizes for $D_1$ and $D_2$ can arouse the bistable characteristic of the stochastic system. When $\tau$ is taken in the bistable interval [0.160, 0.262], taking $\tau =0.17$ as an example, the transition set curves of the system under two noise excitations together are shown in Figure 8(b). As the deterministic system is bistable at this time, the small noise intensities $D_1$ and $D_2$ can arouse the bistable characteristic of the stochastic system, and with the increasing of $D_1$ and $D_2$ , the system becomes monostable again. When $\tau$ is taken in [0.263, 0.5], the deterministic system corresponding to system (5) has also an attractor: limit cycle, and the transition set of the system is the empty set, which also implies that it cannot make the system become bistable by adding the noise excitations at this moment.

From Figure 8, the unfolding parametric plane $D_2-D_1$ is also divided into two sub-regions by the transition set curves, and without loss of generality, we analyze the stationary PDF of amplitude $p\text{(}A\text{)}$ and the joint PDF $p\text{(}X\text{,}\dot{X}\text{)}$ only for one point ( $D_2,D_1$ ) in each of the two sub-regions in Figure 8(a), then compare the analytical solutions with the Monte Carlo simulation results from original system (5) using the numerical method, ³⁸ and show the corresponding results in Figures 9 and 10.

Figure 9.

PDFs for $p\text{(}A\text{)}$ in different sub-regions in Figure 8(a) (taking $D_2$ and $D_1$ as unfolding parameters). (a) Parameter ( $D_2,D_1$ ) in region 1 in Figure 8(a) and (b) parameter ( $D_2,D_1$ ) in region 2 in Figure 8(a).

Figure 10.

Joint PDFs for $p\text{(}X\text{,}\dot{X}\text{)}$ in different sub-regions in Figure 8(a) (taking $D_2$ and $D_1$ as unfolding parameters). (a) Parameter ( $D_2,D_1$ ) in region 1 in Figure 8(a) and (b) parameter ( $D_2,D_1$ ) in region 2 in Figure 8(a).

As can be seen from Figure 8(a), the parametric region where the PDF occurs bimodal is also surrounded by an approximately triangular region. And when we take the parameter ( $D_2,D_1$ ) in region 1, the PDFs for $p\text{(}A\text{)}$ and $p\text{(}X\text{,}\dot{X}\text{)}$ each have a stable equilibrium as shown in Figures 9(a) and 10(a); in region 2, the PDFs for $p\text{(}A\text{)}$ and $p\text{(}X\text{,}\dot{X}\text{)}$ each have a stable limit cycle; moreover, the probability near the origin is not zero, which implies the limit cycle coexists with the equilibrium in the system simultaneously, and as shown in Figures 9(b) and 10(b).

According to all of the above, the results show that the stationary PDFs for $p\text{(}A\text{)}$ and $p\text{(}X\text{,}\dot{X}\text{)}$ in any two adjacent regions in Figures 5 and 8 are qualitatively different. And it also indicates that no matter the exact values of the unfolding parameters cross any curve in these figures, the system will occur stochastic P-bifurcation behavior. Thus, the transition set curves obtained are just the critical parametric conditions for stochastic P-bifurcation of the system. The analytic solutions shown in Figures 6 and 9 are well consistent with the Monte Carlo simulation results from original system (5), thus further verifying the theoretical analysis in this paper.

Compared with integer-order controllers, the fractional-order controllers have better dynamic performances and robustness, and recently, various fractional-order controllers have been developed.^39,40 And we obtained the critical conditions when the system (5) will exhibit stochastic P-bifurcation through the above analysis, which can make the system switch between monostable and bistable states by selecting the corresponding unfolding parameters; this can provide theoretical guidance for the design of fractional-order controllers.

Conclusion

In this paper, we studied the bistable stochastic P-bifurcation of a generalized Van der Pol system with fractional time-delay feedback excited by additive and multiplicative Gaussian white noise excitations simultaneously and discussed the influences of parameters $p$ , $\tau$ , $D_1$ , and $D_2$ on the system. Based on the minimal mean square error principle, the original system was transformed into an equivalent integer-order system, and we obtained the stationary PDF of system amplitude using the stochastic averaging method. Then, the critical parametric conditions for stochastic P-bifurcation of the system were obtained based on singularity theory, according to which we can maintain the system response at a small amplitude near the equilibrium or monostability by selecting the corresponding unfolding parameters, avoiding the instability and damage caused by the large amplitude vibration or nonlinear jump phenomenon of such systems and providing theoretical guidance for system design in practical engineering. The consistency between the Monte Carlo simulation results and the analytical solutions can also verify the theoretical analysis. It shows that the fractional order $p$ , time-delay $\tau$ , and noise intensities $D_1$ and $D_2$ can each arouse stochastic P-bifurcation of the system, and the number of peaks in the stationary PDF curves of system amplitude can be controlled from two to one by selecting the corresponding unfolding parameters. It also illustrates that the method used in this paper is feasible to analyze the stochastic P-bifurcation behaviors of nonlinear oscillators with fractional derivative.