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

and Applied Analysis 3 the integral is replaced by the trapezoidal quadrature formula at point t n+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 fractionalorder 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][5][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][8][9][10][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][17][18][19][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 fractionalorder 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.

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 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 where the symbol [•] means the integer part.
The Riemann-Liouville (RL) definition of fractional derivatives is given by where Γ(•) is the gamma function.
The Caputo () fractional derivative is defined as follows: 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.

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 fractionalorder chaotic system by the improved predictor-corrector algorithm, the following equation is considered: 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 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 where g+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 And if we use the product rectangle rule, the right hand of ( 7) can be written as where Then the predictor and corrector formulae for solving (6) are given, respectively, by The approximation accuracy of scheme (11) is (ℎ min{2,1+} ).

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]: 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: Then, we have The recurrent formula for the Chebyshev polynomials of the second kind is The orthogonality of Chebyshev polynomial of the second kind can be expressed as 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 () ∈  2 can be expanded as where   = ∫ 1 −1 ()()  ().This expansion is called an orthogonal decomposition of random function (), which is the theoretical base of orthogonal decomposition methods.

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: where  is the fractional order and ,  are state variables.,  are deterministic parameters and  is a random parameter, and it can be expressed as 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: where   () represents the th Chebyshev polynomial and  is the largest order of the polynomials.Then substituting (20) into the system (18), we can get 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: (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 Substituting ( 22) and ( 23) into ( 21), then we can have 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 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 The ensemble mean responses (EMRs for short) of the stochastic fractional-order system are described as 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

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 The corresponding characteristic equation is where  0 ,  1 , . . .,  10 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 () ± ().
(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 () =  0   +  1  −1 + ⋅ ⋅ ⋅ +   = 0 be the characteristic equation evaluated at the equilibrium point, Δ  the n-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: (2)   (  ) > 0, (0, 1, . . ., ); (3) Δ  −1 (  ) ̸ = 0.According to Theorem 1, the following parameter relations which can make the equality Δ 9 = 0 are derived: Then, (31) is substituted into the 8-dimensional and 7dimensional Routh-Hurwitz determinants Δ 8 , Δ 7 , and only the following parameters can make Δ 8 ̸ = 0, Δ 7 ̸ = 0 satisfied: (32) Meanwhile, the relations which can make all the coefficients satisfied   > 0 ( = 0, 1, . . ., 10) are We substitute the relations (33) into Δ  9 () and find that Δ  9 () 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: Obviously, when the intensity of random parameter  ̸ = 0, the real parts of all the other eigenvalues are less than 0 except  1 and  2 .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), . . .,  4 (0)] = (0.2, 0.2, 0.2, 0.2)  and [ 0 (0), . . .,  4 (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.09and  = −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.09and 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).
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.
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  0 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  0 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  → ∞.

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.

Figure 1 :
Figure 1: The arch-like PDF curve for random variable .

Figure 4 :Figure 5 :
Figure 4: The phase diagrams for different values of order .