† Corresponding author. E-mail:

Project supported by the National Natural Science Foundation of China (Grant Nos. 11472212, 11532011, and 11502201).

This paper aims to investigate the stochastic response of the van der Pol (VDP) oscillator with two kinds of fractional derivatives under Gaussian white noise excitation. First, the fractional VDP oscillator is replaced by an equivalent VDP oscillator without fractional derivative terms by using the generalized harmonic balance technique. Then, the stochastic averaging method is applied to the equivalent VDP oscillator to obtain the analytical solution. Finally, the analytical solutions are validated by numerical results from the Monte Carlo simulation of the original fractional VDP oscillator. The numerical results not only demonstrate the accuracy of the proposed approach but also show that the fractional order, the fractional coefficient and the intensity of Gaussian white noise play important roles in the responses of the fractional VDP oscillator. An interesting phenomenon we found is that the effects of the fractional order of two kinds of fractional derivative items on the fractional stochastic systems are totally contrary.

Fractional calculus^{[1–5]} can be traced back to the Leibniz’s letter to L’Hospital in 1695. Unfortunately, for quite a long time, fractional calculus remained in an abstract mathematical field and developed very slowly. In the past few decades, fractional calculus has attracted considerable attention, partly due to the fact that the fractional-order mathematical models can be used to make a more accurate description than the integer-order models. Thus, fractional calculus can be considered as an old and yet novel branch of mathematics.^{[6,7]}

Because of the extensive existence of the random factor,^{[8,9]} it is of practical significance to study the stochastic fractional systems. In order to obtain the approximate analytical solution of fractional systems, a lot of approaches have been used by many authors. The stochastic averaging method has been used by many researchers to obtain the analytical solutions. Specifically, Huang and Jin^{[10]} investigated the stationary response and stability of a single-degree-of-freedom (SDOF) strongly nonlinear stochastic system with light damping modeled by a fractional derivative. Chen *et al*.^{[11]} and Chen and Zhu^{[12]} studied the stationary responses. Chen and Zhu^{[13]} investigated the stochastic jump and bifurcation. Chen *et al*.^{[14]} studied first passage failure; Chen *et al*.^{[15,16]} fractional control of stochastic fractional systems. Yang *et al*.^{[17,18]} studied the stationary response of a nonlinear system with Caputo-type fractional derivative under Gaussian white noise excitation by the stochastic averaging method. Zhang *et al*.^{[19]} discussed the response of a Duffing–Rayleigh system with fractional derivative under Gaussian white noise excitation. There are several other methods to deal with stochastic systems besides the stochastic averaging method. Spanos and Zeldin^{[20]} put forward a frequency-domain method for systems with fractional derivatives. Agrawal^{[21]} proposed an analytical approach for stochastic dynamic systems endowed with a fractional derivative by the eigenvector expansion method and Laplace transforms. Xu *et al*.^{[22,23]} studied the response of the stochastic oscillator with fractional derivative by using a new technique based on the Lindstedt–Poincaré (LP) method and multiple-scale method. Liu *et al*.^{[24]} investigated the principal resonance responses of SDOF systems with small fractional derivative damping subjected to narrow-band random parametric excitation by the multiple scale method.

However, the systems mentioned in the above literature only contain one kind of fractional derivative, little work is related to the study of deterministic or stochastic systems containing two kinds of fractional derivative to the best of the authors’ knowledge. Huang *et al*.^{[25]} discussed the statistical behaviors of the response of stochastic systems endowed with two kinds of fractional derivatives by using the Laplace transform in conjunction with the weighted generalized Mittag–Leffler function. Shen *et al*.^{[26]} investigated the primary resonance of the Duffing oscillator with two kinds of fractional-order derivatives. Liao^{[27]} studied the Duffing oscillator with two kinds of fractional-order derivative by using the improved version of the constrained optimization harmonic balance method.

Inspired by the above discussion, in this paper we are to investigate the stationary response of the VDP oscillator with two kinds of fractional derivatives excited by Gaussian white noise.

The rest of this paper is organized as follows. The fractional VDP oscillator is reduced to the equivalent stochastic VDP oscillator without fractional derivative terms by using the generalized harmonic function technique in Section 2. The approximate analytical solutions of the equivalent VDP oscillator are obtained by the stochastic averaging method in Section 3. The predictor–corrector algorithm for the solution of the initial value problem with fractional derivative is briefly introduced in Section 4. Numerical simulations are performed to verify the theoretical results in Section 5. Finally, the conclusions are given in Section 6.

Consider a fractional VDP oscillator^{[28–30]} with two kinds of fractional derivatives subjected to Gaussian white noise. The motion of the system is governed by the following equation:

*t*;

*β*

_{1},

*β*

_{2}, and

*ω*are constant coefficients;

*W*(

*t*) is a stationary Gaussian white noise with correlation function

*E*[

*W*(

*t*)

*W*(

*t*+

*τ*)] = 2

*Dδ*(

*τ*). There are many definitions of the fractional derivative. In the present paper, the Caputo-type fractional derivative is adopted as follows:

*n*− 1 <

*α*<

*n*,

*n*∈

*N*, and Γ(·) is a Gamma function.

According to Refs. [11], [14], and [26], the term associated with the fractional derivative not only serves as the classical damping force, but also contributes to the restoring force. Based on the generalized harmonic function technique,^{[11,14]} this term can be replaced by a linear restoring force and a linear damping force:

In order to obtain the values of *C*_{1} (*α*_{1}), *C*_{2} (*α*_{2}), *K*_{1} (*α*_{1}), and *K*_{2} (*α*_{2}), we first introduce two important formulae:

*C*

_{1}(

*α*

_{1}).

*s*=

*t*−

*u*and d

*s*= −d

*u*, one can obtain

*C*

_{11}and the second part as

*C*

_{12}, then

*C*

_{111}and the second part as

*C*

_{112}, then

Second, we present the detailed derivation procedure for *K*_{1} (*α*_{1}).

*s*=

*t*−

*u*and d

*s*= − d

*u*, one can obtain

*K*

_{11}and the second part as

*K*

_{12}, then

*K*

_{121}and the second part as

*K*

_{122}, then

*C*

_{2}(

*α*

_{2}),

*K*

_{2}(

*α*

_{2}):

In this section, the stochastic averaging method is applied to the VDP oscillator (

According to Ref. [20], we can assume that the solution of system (

*A*,

*Φ*, and

*Θ*are random processes. One can obtain the following equations for the amplitude

*A*and the phase angle

*Θ*:

*A*(

*t*) and

*Φ*(

*t*) converge weakly into a two-dimensional diffusion Markov process as

*ε*→ 0, in a time interval 0 ≤

*t*≤

*T*, where

*T*∼

*O*(

*ε*

^{−1}). Equation (

*A*(

*t*) is independent of

*Θ*(

*t*), the limiting process

*A*(

*t*) is a one-dimensional diffusion process governed by

*p*=

*c*in which

*c*∈ (−∞, +∞) at

*A*= 0 and

*p*→ 0,

*∂ p*/

*∂ A*→ 0 as

*A*→ ∞. Thus, based on these boundary conditions, the stationary solution of Eq. (

*C*is the normalization constant.

The stationary probability density of the system Hamiltonian

In the paper, the following algorithm^{[31]} is used to deal with the numerical solution of differential equations of fractional order, equipped with suitable initial conditions, i.e.,

*a*

_{i,n+1}is in the following form:

*b*

_{i,n+1}represents the weight of the predictor and is expressed as

*O*(

*h*

^{min(2,1+α)}). The pseudo-code is available in the appendix of Ref. [31].

In order to demonstrate the efficiency of the proposed method, some numerical results are obtained for the fractional VDP oscillator (*α*_{1} and *α*_{1}, the fractional coefficients *β*_{1} and *β*_{2}, and the intensity of Gaussian white noise on the fractional VDP oscillator.

In this section, the time step is 0.01 and the total computation time is 10^{4} s in the implementation of the Monte Carlo simulation. The results from the analytical solutions and the Monte Carlo simulations are shown in the following figures by solid line (—) and star (_{*}), respectively.

*α*

_{1}∈ (0,1)

We examine the effect of the order of Caputo-type fractional derivative *α*_{1} on the fractional VDP oscillator (*p*(*A*), displacement *p*(*x*), and velocity *p*(*y*) from an analytical solution for different values of the fractional derivative order are shown in Figs. *α*_{1}, the greater the probability of the oscillator with large amplitude (displacement or velocity) is.

*α*

_{2}∈ (1,2)

We examine the effect of the order of Caputo-type fractional derivative *α*_{2} on the fractional VDP oscillator (*α*_{2}, the smaller the probability of the oscillator with large amplitude (displacement or velocity) is, which is totally contrary to the conclusion in Subsection 5.1.

*β*

_{1}

We investigate the effect of the coefficient of Caputo-type fractional derivative *β*_{1} on the fractional VDP oscillator (*β*_{1}, the larger the probability of the oscillator with large amplitude (displacement or velocity) is.

*β*

_{2}

Now we come to examine the effect of the coefficient of Caputo-type fractional derivative *β*_{2} on the fractional VDP oscillator (*β*_{2}, the larger the probability of the oscillator with large amplitude (displacement or velocity) is, which is the same as the conclusion in Subsection 5.3.

*D*

Finally, we consider the effect of the intensity of Gaussian white noise *D* on the fractional VDP oscillator (*D*, the smaller the probability of the oscillator with large amplitude (displacement or velocity) is.

In this paper, we investigate the stochastic response of VDP oscillator with two kinds of fractional derivatives under Gaussian white noise excitation. The equivalent VDP oscillator without fractional derivative terms for the fractional VDP oscillator is obtained by using the generalized harmonic balance technique. Then, the analytical solution of the equivalent VDP oscillator can be derived by the stochastic averaging method. Finally, the analytical results are validated by those from the Monte Carlo simulation of the original fractional VDP oscillator.

The numerical results not only demonstrate the accuracy of the proposed approach but also show that the fractional orders *α*_{1} and *α*_{2}, the fractional coefficients *β*_{1} and *β*_{2} and the intensity *D* of Gaussian white noise play important roles in the responses of the fractional VDP oscillator.

More specifically, the smaller the fractional derivative order *α*_{1}, the greater the probability of the oscillator with large amplitude (displacement or velocity) is; the effect of fractional order *α*_{2} is contrary to that of fractional order *α*_{1}, i.e., the smaller the fractional derivative order *α*_{2}, the smaller the probability of the oscillator with large amplitude (displacement or velocity) is. The smaller the fractional coefficient *β*_{1} or *β*_{2}, the larger the probability of the oscillator with large amplitude (displacement or velocity) is. The smaller the intensity *D* of Gaussian white noise, the smaller the probability of the oscillator with large amplitude (displacement or velocity) is.

It is noted that the approach is only suitable for the fractional VDP oscillator with two kinds of fractional derivatives subjected to weakly Gaussian white noise and light damping.

**Reference**