Hopf Bifurcations of a Stochastic Fractional-Order Van der Pol System

Abstract

The Hopf bifurcation of a fractional-order Van der Pol (VDP for short) system with a random parameter is investigated. Firstly, the Chebyshev polynomial approximation is applied to study the stochastic fractional-order system. Based on the method, the stochastic system is reduced to the equivalent deterministic one, and then the responses of the stochastic system can be obtained by numerical methods. Then, according to the existence conditions of Hopf bifurcation, the critical parameter value of the bifurcation is obtained by theoretical analysis. Then, numerical simulations are carried out to verify the theoretical results.

1. Introduction

Fractional calculus is a topic of more than 300 years old. The idea of fractional calculus has been known since the regular calculus. However, its development is very rapid just in recent several decades due to its application in physics, engineering, secure communications, and so on. Meanwhile, it has been found that many fractional-order nonlinear systems can demonstrate chaotic behavior, such as fractional-order Chua circuit [1], fractional-order Lorenz system [2], and fractional-order Chen system [3]. These examples and many other similar samples perfectly clarify the importance of consideration and analysis of dynamical systems with fractional-order models [4–6].

The VDP oscillator represents a nonlinear system with an interesting behavior that exhibits naturally in several applications. It has been used for study and design of many models including biological phenomena, such as the heartbeat, neurons, acoustic models, and radiation of mobile phones, and as a model of electrical oscillators. The deterministic VDP has very rich dynamical behaviors, such as saddle-node bifurcation, symmetry-breaking bifurcation, period-doubling bifurcation, Hopf bifurcation, and chaos. Recently, many researchers studied the bifurcation and periodic solutions of the VDP system in detail [7–11]. As research on fractional differential equations goes up, the VDP model with fractional order was proposed and investigated in [12, 13]. It was found that when the order is less than 1, the fractional-order VDP system also has rich dynamics similar to the corresponding integer-order one.

In the real world, uncertainty due to measuring, materials, manufacturing, and assembling is inevitable. These uncertainties usually can be modeled as random parameters with certain statistics. A stochastic system means a system with random parameters of given statistics, which is a mathematical model close to the real world. There are several methods to deal with the dynamical problems of stochastic systems. The first one is the Monte Carlo method [14], which is simple and universal but usually involves a great deal of computational effort. The second one is the stochastic finite element method [15], which involves the least computation but is usually restricted to systems with random variables of small perturbations only. And the third one is orthogonal polynomials [16–20], which is based on the expansion theory of orthogonal polynomials, without the limitation of small perturbations. There are a lot of references about the analysis of stochastic systems via orthogonal polynomials. For example, in [21, 22], the Chebyshev polynomial approximation was applied to analyze the bifurcation and chaos of stochastic Duffing system and stochastic Duffing-Van der Pol system. The period-doubling bifurcation for double-well stochastic Duffing system was investigated via Lagurre polynomial approximation in [23]. Recently, the synchronization for a stochastic fractional-order system with a random parameter via Laguerre polynomial was analyzed in [24]. However, the literature on the investigation of a stochastic fractional-order VDP system is few.

Inspired by the above discussion, in this paper, Hopf bifurcations of a fractional-order VDP system with a random parameter are studied. Firstly, the Chebyshev polynomial approximation method is applied to investigate the stochastic fractional-order system. Based on this method, the stochastic system is reduced to the equivalent deterministic one. Then, according to the existence conditions of Hopf bifurcation, the critical parameter value of the bifurcation is obtained by theoretical analysis. Meanwhile, numerical simulations are performed to verify the theoretical results. Besides, the Hopf bifurcation with the variation of derivative order is also observed by numerical computations.

The paper is organized as follows. In Section 2, the definitions for the fractional calculus and numerical algorithms are given. In Section 3, the Chebyshev polynomials are introduced as the basis of orthogonal polynomial approximation. In Section 4, the transformation of stochastic fractional-order VDP system into its equivalent deterministic one by Chebyshev polynomial approximation is shown. Section 5 is devoted to studying the Hopf bifurcation of the stochastic system. Finally, we summarize the results in Section 6.

2. Fundamentals of Fractional Derivative

2.1. Definition

Fractional calculus is a generalization of integration and differentiation to a noninteger-order integrodifferential operator _𝑎𝐷^𝑞_𝑡 which is defined by _𝑎𝐷^𝑞_𝑡=𝑑^𝑞𝑑𝑡^𝑞𝑅𝑞𝑞>01𝑅=0_𝑡𝑎(𝑑𝜏)^−𝑞𝑅𝑞<0, (1) where 𝑞 is the derivative order which can be a complex number and 𝑅(𝑞) is the real part of 𝑞. The numbers 𝑎 and 𝑡 are the limits of the operator. There are many definitions for fractional derivatives. Three most frequently used ones are the Grunwald-Letnikov definition, the Riemann-Liouville, and the Caputo definitions.

The Grunwald-Letnikov definition (GL) derivative with fractional-order 𝑞 is described by _𝑎^GL𝐷^𝑞_𝑡𝑓(𝑡)=lim_ℎ→01ℎ^{𝑞[(𝑡−𝑎)/ℎ]}_𝑗=0(−1)^𝑗𝑞𝑗𝑓𝑡−𝑗ℎ, (2) where the symbol [•] means the integer part.

The Riemann-Liouville (RL) definition of fractional derivatives is given by _𝑎^RL𝐷^𝑞_𝑡𝑓(𝑡)=𝑑^𝑛𝑑𝑡^𝑛1Γ𝑛−𝑞^𝑡_𝑎𝑓(𝜏)(𝑡−𝜏)^{𝑞−𝑛+1}𝑑𝜏,𝑛−1<𝑞<𝑛, (3) where Γ(•) is the gamma function.

The Caputo (𝐶) fractional derivative is defined as follows: ^𝐶_𝑎𝐷^𝑞_𝑡𝑓(𝑡)=1Γ𝑛−𝑞^𝑡_𝑎(𝑡−𝜏)^{𝑛−𝑞−1}𝑓^(𝑛)(𝜏)𝑑𝜏,𝑛−1<𝑞<𝑛. (4)

It is well known that the initial conditions for the fractional differential equations with Caputo derivatives take the same form as for the integer-order ones, which is very suitable for practical problems [25]. Therefore, we will use the Caputo definition for the fractional derivatives in this paper.

2.2. Numerical Algorithms

Obtaining numerical solutions of a fractional differential equation is not easily compared with that for ordinary differential equations. There are two approximation methods which can frequently be used to numerical computations on chaos and bifurcations with fractional differential equations. One is an improved version of Adams-Bashforth-Moulton algorithm based on the predictor-correctors scheme [26], which is a time-domain approach. The other is a method, known as frequency domain approximation [27], based on numerical analysis of fractional-order systems in the frequency domain.

Simulations for fractional-order systems using the time domain methods are complicated and due to long memory characteristics of these systems require a very long simulation time but on the other hand, it is more accurate [28]. Therefore, we employ the improved predictor-corrector algorithm for fractional-order differential equations in this paper.

In order to get the approximate solution of a fractional-order chaotic system by the improved predictor-corrector algorithm, the following equation is considered: 𝑑^𝑞𝑥𝑑𝑡^𝑞()𝑥=𝑓𝑡,𝑥,0≤𝑡≤𝑇_(𝑘)(0)=𝑥₀^(𝑘),𝑘=0,1,…,𝑞−1, (5) where ⌈𝑞⌉ is just the value 𝑞 rounded up to the nearest integer and 𝑥^(𝑘) is the ordinary 𝑘th derivative of 𝑥. Formula (5) is equivalent to the Volterra integral equation 𝑥(𝑡)=^{⌈𝑞⌉−1}_𝑘=0𝑥₀^(𝑘)𝑡^𝑘+1𝑘!Γ𝑞^𝑡₀(𝑡−𝜏)^𝑞−1𝑓(𝜏,𝑥(𝜏))𝑑𝜏. (6) Now, for the sake of simplicity, we assume that we are working on a uniform grid {𝑡_𝑛=𝑛ℎ∶𝑛=0,1,…,𝑀} with some integer 𝑀 and set ℎ=𝑇/𝑀. Using the standard quadrature techniques for the integral in (6), denote 𝑔(𝜏)=𝑓(𝜏,𝑥(𝜏)); the integral is replaced by the trapezoidal quadrature formula at point 𝑡_𝑛+1^𝑡𝑛+1₀𝑡_𝑛+1−𝜏^𝑞−1𝑔(𝜏)𝑑𝜏≈^𝑡𝑛+1₀𝑡_𝑛+1−𝜏^𝑞−1𝑔_𝑛+1(𝜏)𝑑𝜏, (7) where 𝑔_𝑛+1 is the piecewise linear interpolant for 𝑔 with nodes and knots chosen at the 𝑡_𝑗, (𝑗=0,1,…,𝑛+1). After some elementary calculations, the right hand side of (7) gives ^𝑡𝑛+1₀𝑡_𝑛+1−𝜏^𝑞−1𝑔_𝑛+1(𝜏)ℎ𝑑𝜏=^𝑞𝑞𝑞+1^𝑛+1_𝑗=0𝛼_𝑗,𝑛+1𝑔𝑡_𝑗. (8) And if we use the product rectangle rule, the right hand of (7) can be written as ^𝑡𝑛+1₀𝑡_𝑛+1−𝜏^𝑞−1𝑔_𝑛+1(𝜏)𝑑𝜏=^𝑛+1_𝑗=0𝛽_𝑗,𝑛+1𝑔𝑡_𝑗, (9) where 𝛼_𝑗,𝑛+1=𝑛^𝑞+1−𝑛−𝑞(𝑛+1)^𝑞,𝑗=0𝑛−𝑗+2^𝑞+1+𝑛−𝑗^𝑞+1−2𝑛−𝑗+1^𝑞+1,𝛽1≤𝑗≤𝑛,1,𝑗=𝑛+1,_𝑗,𝑛+1=ℎ^𝑞𝑞𝑛−𝑗+1^𝑞−𝑛−𝑗^𝑞,1≤𝑗≤𝑛. (10) Then the predictor and corrector formulae for solving (6) are given, respectively, by 𝑥^𝑃_ℎ𝑡_𝑛+1=^{⌈𝑞⌉−1}_𝑘=0𝑥₀^(𝑘)𝑡^𝑘_𝑛+1+1𝑘!Γ𝑞^𝑛_𝑗=0𝛽_𝑗,𝑛+1𝑓𝑡_𝑗,𝑥_ℎ𝑡_𝑗,𝑥_ℎ𝑡_𝑛+1=^{⌈𝑞⌉−1}_𝑘=0𝑥₀^(𝑘)𝑡^𝑘_𝑛+1+ℎ𝑘!Γ^𝑞𝑓𝑡𝑞+2^𝑛+1,𝑥^𝑃_ℎ𝑡_𝑛+1+ℎ^𝑞Γ𝑞+2^𝑛_𝑗=0𝛼_𝑗,𝑛+1𝑓𝑡_𝑗,𝑥_ℎ𝑡_𝑗. (11) The approximation accuracy of scheme (11) is 𝑂(ℎ^{min{2,1+𝑞}}).

3. The Chebyshev Polynomials

In general, random parameters in some engineering structures are bounded in nature. So the probability density function (PDF for short) model for the bounded random variables is taken to be an arch-like PDF of the form [29]: 𝜌(𝑢)=2𝜋√1−𝑢²|𝑢|0|𝑢|≤1>1. (12)

Figure 1 shows the arch-like PDF curve for bounded random variable 𝑢. As the orthogonal polynomial basis for this kind of PDF, the second kind of Chebyshev polynomials is the only choice. The general expansion for the second kind of Chebyshev polynomials can be described as follows: 𝐻_𝑛(𝑢)=^[𝑛/2]_𝑘=0(−1)^𝑘(𝑛−𝑘)!()!()𝑘!𝑛−2𝑘2𝑢^𝑛−2𝑘. (13) Then, we have 𝐻₀(𝑢)𝐻=1,₁(𝑢)𝐻=2𝑢,₂(𝑢)=4𝑢²𝐻−1,₃(𝑢)=8𝑢³𝐻−4𝑢,₄(𝑢)=16𝑢⁴−12𝑢²+1.⋯. (14) The recurrent formula for the Chebyshev polynomials of the second kind is 𝑢𝐻_𝑛(𝑢)=12𝐻_𝑛−1(𝑢)+𝐻_𝑛+1(𝑢). (15) The orthogonality of Chebyshev polynomial of the second kind can be expressed as ¹₋₁2𝜋√1−𝑢²𝐻_𝑖(𝑢)𝐻_𝑗(𝑢)10.𝑑𝑢=𝑖=𝑗𝑖≠𝑗 (16) Equation (16) represents a weighted orthogonal relationship. The weighting function is just the same as that for the PDF of the random parameter which obeys arch distribution, 𝜌(𝑢) in (12), and the left-hand side of (16) can be regarded as the expectation of the product 𝐻_𝑖(𝑢)𝐻_𝑗(𝑢). Owing to the orthogonality of Chebyshev polynomial, any measurable function 𝑓(𝑢)∈𝐿² can be expanded as 𝑓(𝑢)=^∞_𝑖=0𝑐_𝑖𝐻_𝑖(𝑢), (17) where 𝑐_𝑖=∫¹₋₁𝜌(𝑢)𝑓(𝑢)𝐻_𝑖(𝑢)𝑑𝑢. This expansion is called an orthogonal decomposition of random function 𝑓(𝑢), which is the theoretical base of orthogonal decomposition methods.

Figure 1 

The arch-like PDF curve for random variable 𝑢.

4. The Stochastic Fractional-Order VDP System

In this paper, we consider the fractional-order VDP system with a random parameter. The system can be described by the following differential equations: 𝐷^𝑞𝐷𝑥=𝑦_𝑞𝑦=𝜇𝑦−𝑎𝑥²𝑦+𝑏𝑥, (18) where 𝑞 is the fractional order and 𝑥, 𝑦 are state variables. 𝑎, 𝑏 are deterministic parameters and 𝜇 is a random parameter, and it can be expressed as 𝜇=𝜇+𝛿𝑢, (19) where 𝜇 and 𝛿/2 are the mean value and standard deviation of 𝜇, respectively, and 𝛿 is regarded as the intensity of 𝜇. 𝑢 is a random variable defined on [−1,1].

It is well known that the stochastic function space, which is constituted with the responses of nonlinear dynamic systems with random parameters, has proved to be a Hilbert space with respect to convergence in the mean square. Therefore, the fractional-order system (18) with random parameters can be reduced into the equivalent deterministic system by using of the orthogonal polynomial expansion. It follows from the orthogonal polynomial approximation that the responses of system (18) can be expressed approximately by the following series: 𝑥(𝑡,𝑢)=^𝑁_𝑖=0𝑥_𝑖(𝑡)𝐻_𝑖(𝑢)𝑦(𝑡,𝑢)=^𝑁_𝑖=0𝑦_𝑖(𝑡)𝐻_𝑖(𝑢), (20) where 𝐻_𝑖(𝑢) represents the 𝑖th Chebyshev polynomial and 𝑁 is the largest order of the polynomials. Then substituting (20) into the system (18), we can get 𝐷^𝑞𝑁_𝑖=0𝑥_𝑖(𝑡)𝐻_𝑖(𝑢)=−^𝑁_𝑖=0𝑦_𝑖(𝑡)𝐻_𝑖(𝑢)𝐷^𝑞𝑁_𝑖=0𝑦_𝑖(𝑡)𝐻_𝑖(𝑢)=𝜇^𝑁_𝑖=0𝑦_𝑖(𝑡)𝐻_𝑖(𝑢)+𝛿𝑢^𝑁_𝑖=0𝑦_𝑖(𝑡)𝐻_𝑖(𝑢)−𝑎^𝑁_𝑖=0𝑥_𝑖(𝑡)𝐻_𝑖(𝑢)^2𝑁_𝑖=0𝑦_𝑖(𝑡)𝐻_𝑖(𝑢)+𝑏^𝑁_𝑖=0𝑥_𝑖(𝑡)𝐻_𝑖(𝑢), (21) and with the aid of the recurrent formulas of Chebyshev polynomial, all the nonlinear terms in the system can be further reduced into linear combination of 𝐿_𝑖(𝑢) as follows: ^𝑁_𝑖=0𝑥_𝑖(𝑡)𝐻_𝑖(𝑢)^2𝑁_𝑖=0𝑦_𝑖(𝑡)𝐻_𝑖(𝑢)²=^3𝑁_𝑖=0𝑆_𝑖(𝑡)𝐻_𝑖(𝑢), (22) where 𝑆_𝑖(𝑡) can be obtained with the help of the computer algebraic systems, such as Maple or MATLAB. Meanwhile, the term 𝛿𝑢(∑^𝑁_𝑖=0𝑦_𝑖(𝑡)𝐻_𝑖(𝑢)) can be reduced to 𝑢𝛿^𝑁_𝑖=0𝑦_𝑖(𝑡)𝐻_𝑖(𝑢)=𝛿^𝑁_𝑖=0𝑦_𝑖(𝑡)𝐻_𝑖−1(𝑢)+𝐻_𝑖+1(𝑢)2. (23) Substituting (22) and (23) into (21), then we can have 𝐷^𝑞𝑁_𝑖=0𝑥_𝑖(𝑡)𝐻_𝑖(𝑢)=−^𝑁_𝑖=0𝑦_𝑖(𝑡)𝐻_𝑖(𝑢)𝐷^𝑞𝑁_𝑖=0𝑦_𝑖(𝑡)𝐻_𝑖(𝑢)=𝜇^𝑁_𝑖=0𝑦_𝑖(𝑡)𝐻_𝑖(𝑢)+𝛿^𝑁_𝑖=0𝑦_𝑖(𝑡)𝐻_𝑖−1(𝑢)+𝐻_𝑖+1(𝑢)2−𝑎^𝑁_𝑖=0𝑆_𝑖(𝑡)𝐻_𝑖(𝑢)+𝑏^𝑁_𝑖=0𝑥_𝑖(𝑡)𝐻_𝑖(𝑢). (24) Multiply both sides of (24) by 𝐻_𝑖(𝑢), 𝑖=0,1,2,…,𝑁 in sequence and take expectation with respect to 𝑢. Then we can finally obtain the equivalent deterministic system. Remember that if 𝑁→∞ in (24), then its left is strictly equivalent to the response of stochastic system (18). Otherwise, (24) is just approximately valid with a minimal residual error. Due to the requirement of computational precision, we take 𝑁=4 in the following numerical analysis. Therefore, we get the equivalent deterministic system approximately as 𝐷^𝑞𝑥₀=−𝑦₀𝐷^𝑞𝑦₀=−𝑎𝑆₀+𝑏𝑥₀+𝜇𝑦₀+12𝛿𝑦₁𝐷^𝑞𝑥₁=−𝑦₁𝐷^𝑞𝑦₁=−𝑎𝑆₁+𝑏𝑥₁+𝜇𝑦₁+12𝛿𝑦₀+𝑦₂𝐷^𝑞𝑥₂=−𝑦₂𝐷^𝑞𝑦₂=−𝑎𝑆₂+𝑏𝑥₂+𝜇𝑦₂+12𝛿𝑦₁+𝑦₃𝐷^𝑞𝑥₃=−𝑦₃𝐷^𝑞𝑦₃=−𝑎𝑆₃+𝑏𝑥₃+𝜇𝑦₃+12𝛿𝑦₂+𝑦₄𝐷^𝑞𝑥₄=−𝑦₄𝐷^𝑞𝑦₄=−𝑎𝑆₄+𝑏𝑥₄+𝜇𝑦₄+12𝛿𝑦₃. (25) We can get the numerical solutions 𝑥_𝑖, 𝑦_𝑖, (𝑖=0,1,2,3,4) of system (25) by the improved predictor-corrector algorithm. Therefore, the approximate random responses of the stochastic system (18) can be expressed as 𝑥(𝑡,𝑢)≈⁴_𝑖=0𝑥_𝑖(𝑡)𝐻_𝑖(𝑢)𝑦(𝑡,𝑢)≈⁴_𝑖=0𝑦_𝑖(𝑡)𝐻_𝑖(𝑢). (26) The ensemble mean responses (EMRs for short) of the stochastic fractional-order system are described as 𝐸[𝑥(𝑡,𝑢)]≈⁴_𝑖=0𝑥_𝑖(𝑡)𝐸𝐻_𝑖(𝑢)=𝑥₀(𝑡)𝐸𝑦(𝑡,𝑢)≈⁴_𝑖=0𝑦_𝑖(𝑡)𝐸𝐻_𝑖(𝑢)=𝑦₀(𝑡). (27) Especially for 𝑢=0, namely, 𝜇=𝜇, this sample system is significant for reference. It is usually called a mean parameter system, of which the responses (SRMs for short) may be approximated as 𝑥(𝑡,0)≈⁴_𝑖=0𝑥_𝑖(𝑡)𝐻_𝑖(0)=𝑥₀(𝑡)−𝑥₂(𝑡)+𝑥₄(𝑡)𝑦(𝑡,0)≈⁴_𝑖=0𝑦_𝑖(𝑡)𝐻_𝑖(0)=𝑦₀(𝑡)−𝑦₂(𝑡)+𝑦₄(𝑡). (28)

5. Hopf Bifurcations

A Hopf bifurcation occurs when a periodic solution or limit cycle, surrounding an equilibrium point, arises or goes away as a parameter varies. In this section, the Hopf bifurcation for the stochastic fractional-order VDP system will be studied with theoretical and numerical computations. Obviously, both the stochastic VDP and equivalent deterministic systems (18) and (25) only have one equilibrium point, namely, the origin (0,0). The Jacobian matrix 𝐽 for the linearized system evaluated at the origin is 𝑏𝐽=0−1000000001𝜇021𝛿000000000−10000002𝛿0𝑏1𝜇021𝛿000000000−100000002𝛿𝑏1𝜇021𝛿000000000−100000002𝛿𝑏1𝜇02𝛿1000000000−100000002𝛿𝑏𝜇. (29) The corresponding characteristic equation is 𝑓(𝜆)=𝑎₀𝜆¹⁰+𝑎₁𝜆⁹+𝑎₂𝜆⁸+𝑎₃𝜆⁷+𝑎₄𝜆⁶+𝑎₅𝜆⁵+𝑎₆𝜆⁴+𝑎₇𝜆³+𝑎₈𝜆²+𝑎₉𝜆+𝑎₁₀, (30) where 𝑎₀,𝑎₁,…,𝑎₁₀ are the coefficients of (30) and can be obtained by computing via Maple or MATLAB. According to the bifurcation theory of nonlinear dynamical systems and taking 𝑤 as a bifurcation parameter [22], we can obtain the existence conditions of Hopf bifurcation for the system. (1)The eigenvalues of the linearized system about the fixed point are 𝛼(𝑤)±𝑖𝛽(𝑤).(2)The real and image parts of the eigenvalues satisfy 𝛼(𝑤_𝑐)=0, 𝛽(𝑤_𝑐)>0, 𝛼^′(𝑤_𝑐)≠0, when 𝑤=𝑤_𝑐.(3)The real parts of other eigenvalues are not equal to zeros.

For the high order of equivalent system (25), so it is hard to get the threshold value of bifurcation parameter directly. Therefore, we can get the threshold via the following theorem.

Theorem 1. Let 𝑓(𝜆)=𝑎₀𝜆^𝑛+𝑎₁𝜆^𝑛−1+⋯+𝑎_𝑛=0 be the characteristic equation evaluated at the equilibrium point, Δ_𝑛 the 𝑛-dimensional Routh-Hurwitz determinant, and 𝑎_𝑖, (𝑖=0,1,…,𝑛) the coefficient of the characteristic equation. If there exists a threshold of bifurcation parameter 𝑤_𝑐 which conforms to the following terms, then the first two existence conditions of Hopf bifurcation are satisfied:(1)Δ_𝑛−1(𝑤_𝑐)=0, Δ_𝑛−2(𝑤_𝑐)≠0, Δ_𝑛−3(𝑤_𝑐)≠0;(2)𝑎_𝑖(𝑤_𝑐)>0, (0,1,…,𝑛);(3)Δ^′_𝑛−1(𝑤_𝑐)≠0.

According to Theorem 1, the following parameter relations which can make the equality Δ₉=0 are derived: 𝛿𝜇=0,𝜇=±2𝛿,𝜇=±4√,𝜇=±3𝛿2,√𝜇=±3𝛿4𝛿,𝜇=4±√3𝛿4𝛿,𝜇=−4±√3𝛿4. (31) Then, (31) is substituted into the 8-dimensional and 7-dimensional Routh-Hurwitz determinants Δ₈, Δ₇, and only the following parameters can make Δ₈≠0, Δ₇≠0 satisfied: 𝛿𝜇=±2√𝜇=±3𝛿2. (32) Meanwhile, the relations which can make all the coefficients satisfied 𝑎_𝑖>0 (𝑖=0,1,…,10) are √𝜇=±3𝛿2. (33) We substitute the relations (33) into Δ^′₉(𝜇) and find that Δ^′₉(𝜇) is not equal to 0; therefore, according to Theorem 1, the first two existence conditions of Hopf bifurcation are satisfied. When (33) is substituted into the characteristic equation, all the eigenvalues of it can be obtained: 𝜆_1,2=±√𝜆𝑏𝑖,_3,4=𝛿4−√3𝛿4±144𝛿²−2√3𝛿²𝜆−16𝑏_5,6√=−3𝛿4±143𝛿²𝜆−16𝑏_7,8𝛿=−4−√3𝛿4±144𝛿²−2√3𝛿²𝜆−16𝑏_9,10√=−3𝛿2±123𝛿²−4𝑏. (34) Obviously, when the intensity of random parameter 𝛿≠0, the real parts of all the other eigenvalues are less than 0 except 𝜆₁ and 𝜆₂. The characteristic equation does not have zero solution for 𝑏≠0. Therefore, all the existence conditions of Hopf bifurcation are satisfied. When 𝜇_𝑐=−√3𝛿/2, the Hopf bifurcation occurs around the equilibrium (0,0).

The parameters are taken as 𝑎=1, 𝑏=1, and 𝛿=0.1, and the derivative order is 𝑞=0.995. And initial conditions are [𝑥₀(0),…,𝑥₄(0)]=(0.2,0.2,0.2,0.2)^𝑇 and [𝑦₀(0),…,𝑦₄(0)]=(−0.2,−0.2,−0.2,−0.2)^𝑇. According to the above analysis, we can get that the critical parameter value of the bifurcation is 𝜇_𝑐=−√3𝛿/2≈−0.0866, which implies that if 𝜇<𝜇_𝑐, the system will converge to the equilibrium; otherwise, it will converge to a limit cycle. The parameter values are chosen as 𝜇=−0.09 and 𝑢=−0.01, respectively. Then we obtain the time trajectories and phase diagrams for the EMR and SRM of the stochastic system by simulation computations; see Figures 2 and 3. It is clearly seen that the system converges to the fixed point (0,0) for 𝜇=−0.09 and to a limit cycle for 𝑢=−0.01, which implies that the Hopf bifurcation occurs in the stochastic system (18). Meanwhile, the SRM and EMR of phase plots coincide well enough. It means that the Chebyshev polynomial approximation method is effective for the stochastic fractional-order system (18).

(a) The time trajectory of EMR

(b) The phase diagram of EMR in 𝑥0−𝑦0 plane

(c) The time trajectory of SRM

(d) The phase diagram of SRM

Figure 2 

The time trajectories and phase diagrams for 𝑢=−0.09.

(a) The time trajectory of EMR

(b) The phase diagram of EMR in 𝑥0−𝑦0 plane

(c) The time trajectory of SRM

(d) The phase diagram of SRM

Figure 3 

The time trajectories and phase diagrams for 𝑢=−0.01.

When the others parameters are fixed and 𝜇=−0.02, we also found that Hopf bifurcation occurs with the variation of the derivative order. The stochastic system converges to the equilibrium when 𝑞≤0.98 and to a limit cycle for 𝑞≥0.99, which means that the Hopf bifurcation occurs when the derivative order 𝑞=0.99; see Figure 4, from which it is obtained that the derivative order can also be taken as a bifurcation parameter for fractional-order system.

(a) 𝑞=0.98

(b) 𝑞=0.99

Figure 4 

The phase diagrams for different values of order 𝑞.

The simulation time is taken as 𝑡=800, one integration step of time is taken as Δ𝑡=0.001, derivative order 𝑞=0.97, and divide the time interval [0,800] into 8 equal parts. The time trajectory and phase diagram for 𝑡=800 are shown in Figure 5, from which it appears that the 𝑥₀ converges to a fixed point as 𝑡≥300. However, the time trajectory is zoomed in when 𝑡 belongs to the interval [300,400]; see Figure 6(a), from which we can see that the trajectory gradually decreases but does not converge to a point directly. Secondly, the time trajectory is zoomed in when 𝑡∈[400,500]; see Figure 6(b). Meanwhile, we also obtain the trajectories for 𝑡∈[500,600], 𝑡∈[600,700], and 𝑡∈[700,800] which are shown in Figures 6(c)–6(e), respectively. Comparing these figures, it is seen that the amplitude of 𝑥₀ decreases as the time 𝑡 increases. Based on these observations, we can suggest that the EMRs of the stochastic system (18) have an extremely long transient process while converging to the equilibrium as 𝑡→∞.

(a)

(b)

Figure 5 

The time trajectory and phase diagram of EMR when 𝑡∈[0,800]. (a) The time trajectory. (b) The phase diagram in 𝑥0−𝑦0 plane.

(a) 𝑡∈[300,400]

(b) 𝑡∈[400,500]

(c) 𝑡∈[500,600]

(d) 𝑡∈[600,700]

(e) 𝑡∈[700,800]

Figure 6 

The time trajectories of EMR 𝑥0 for different time intervals.

6. Conclusions

In this paper, the Hopf bifurcation of a fractional-order VDP system with a random parameter is analyzed. An equivalent deterministic system is obtained with the approximation principle of Chebyshev orthogonal polynomials. Based on bifurcation theory of nonlinear dynamic systems, the critical parameter value of the bifurcation is obtained by theoretical analysis. Meanwhile, numerical simulations are performed to verify the theoretical results. Besides, the derivative order can also be taken as a bifurcation parameter, and the Hopf bifurcation is observed with the variation of the order.

It is worth mentioning that the orthogonal polynomial approximation is an effective method for integer-order stochastic systems, but the universality of the method for fractional-order ones still remains an open question. Due to the absences of analytical methods, the threshold value of bifurcation when the derivative order is taken as a bifurcation parameter cannot be obtained by theoretical analysis but only the numerical experiments.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This work is supported by the National Natural Science Foundation of China (NSFC) under Grant nos. 11332008, 11172223, and 11172224.