# Stochastic Response of Dynamical Systems with Fractional Derivative Term under Wide-Band Excitation

Authors: Jun-Feng Zhao
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Transient solution of a fractional stochastic dynamical system under wide-band noise excitation is investigated. Generalized Harmonic Balance technique is firstly used to approximate restoring force of the given system as an amplitude-dependent form. In this way, stochastic averaging method then can be applied to transform the system into an Ito differential equation. Furthermore, the fractional derivative in the integral-differential form can be equivalent to a combination of periodic functions after the averaging procedure. As the following, Galerkin method therein is utilized to obtain the transient probability density functions by solving associated Fokker-Planck-Kolmogorov (FPK) equation. As an example, the Rayleigh oscillator is studied to illustrate the efficiency and accuracy of the proposed approaches. Numerical results show that exact stationary solution and transient solution derived from Galerkin method are in good agreement with those from Monte Carlo Simulation.


---

## Body

## 1. Introduction

Engineering structures such as high buildings and huge bridges often vibrate if they are subjected to random excitation, such as earthquake or strong wind. The vibration behavior of these engineering structures is usually characterized by a dynamical system under noise excitation in mathematical framework. Therefore, how to determine the solution or response, that is, vibrating displacement and velocity, of engineering structural system is an important issue in the field of structural dynamics and mechanical engineering.Generally, there are two different kinds of system response that researchers often considered. One is stationary response which can be obtained by setting the time tending to infinity and therein only displays the future feature of system response. The other is transient response which shows the dynamical evolution of system solution with respect to each time instant. Comparatively, transient response is more important due to its comprehensive expression of system evolution with the time than stationary response but also more difficult to get owing to time involved.It has been proved [1, 2] that only few one-order nonlinear systems or higher-dimensional linear systems have achieved the exact solutions by now; most generally structural systems, for example, the fractional systems we considered in the present paper, have to develop approximated techniques to get numerical solution. In the past decades, Calculus of Variations [3], perturbation method [4], Non-Gauss Closure method [5], and Galerkin method [6] have ever been applied to derive approximated transient solution for systems without any fractional derivative terms and without wide-band noise in them. As for those systems endowed with fractional derivative terms, new methods have to be developed and applied to explore dynamical behaviors of systems especially transient response. The difficulty associated with fractional derivative is mainly caused by its complicated integral-derivative mathematical form with time delay in the integration; this mathematical formula means the first task to deal with fractional derivative is to transform or replace it by other special functions. Therefore, how to deal with fractional derivative becomes a big challenge to those researchers who focus on studying fractional dynamical systems. According to my knowledge, only few works by now have been concerned with discussion on the solution of fractional dynamical system and most of them mainly pay attention to analyze stationary response or resonance response under Gaussian white noise.For example, Chen et al. [7] analyzed stationary response of a Duffing oscillator with hardening stiffness and fractional derivative subjected to Gaussian white noise. Hu et al. [8] investigated stationary response of a strongly nonlinear oscillator with fractional derivative damping under bounded noise excitation. Yang et al. [9] estimated stationary response of a nonlinear system with Caputo-type fractional derivative under Gaussian white noise. These three works are all achieved on the basis of stochastic averaging method [10]. On the other hand, Liu et al. [11] used multiple scales to investigate principal resonance response of SDOF system with small fractional derivative damping under narrow-band random parametric excitation. Xu et al. [12] utilized Lindstedt-Poincare and multiple scales method to investigate a stochastic viscoelastic system with fractional damping; they obtained resonance response and moment response. Matteo et al. [13] applied Wiener path integral to study nonstationary response of nonlinear oscillators with fractional derivative elements under Gaussian white-noise excitation. Li et al. [14] have ever attempted to study the solution for one kind of stochastic dynamical systems with fractional derivative by using equivalent linearization method in 2014. To sum up, all these works have ever given much insight into the cases of stationary response and Gaussian white-noise excitation for fractional dynamical systems. Nevertheless, as for transient response and more practically wide-band noise, due to their instantaneous feature with the time change and different forms of power spectrums, analysis of them is very rare and valuable.Different from these works, the present paper aims to extend Galerkin method often used in deterministic systems to study transient response for those systems under wide-band noise excitation in which damping or restoring force is of fractional order. By selecting the eigenfunctions of FPK equation with respect to the oscillator envelope as a set of basic functions, this method may provide substantially computational advantages in the implementation of the approximate procedure based on eigenfunctions’ orthogonal properties. Except that, we use an effective approximate method to replace complicated fractional derivative as a series of periodic functions.This paper is organized as follows. Harmonic Balance technique is applied to transform the original system into an equivalent linear one in Section2. After that, stochastic averaging method will reduce the envelope process to satisfy a diffusive differential equation in Section 3. In Section 4, Galerkin method is used to get transient response of system by way of solving associated Fokker-Planck equation. At the end of this paper, numerical results associated with Rayleigh oscillator derived from proposed approaches are drawn and compared with results from Monte Carlo Simulation to show the efficiency.


## 2. Equivalent System

Consider the lightly damped oscillators governed by a second-order differential equation. The motion equation of such oscillators with fractional derivative under wide-band noise excitation may be written in the form(1)x¨t+2ζω0x˙t+βDαxt+ω02xt+εfxt,x˙t=Wt,where x is generalized displacement and the superscript dot represents differentiation with respect to time t. ζ, ω0, and β are all constants: they are critical damping rate, natural angular frequency, and coefficient of fractional derivative, respectively. It is assumed that ε is a small real number.There are three classical definitions on fractional derivative in the field of mathematics. The first one is in the sense of Grunwald-Letnikov, defined as the limit of difference form of integer-order calculus and often used in numerical calculation. The second is Riemann-Liouville definition, defined as derivative of a generalized integral with a time delay in it, often used in theoretic analysis in continuous systems. The third one is Caputo definition, defined as a generalized integral with the derivative of time delay in integration, often used in theoretic analysis in continuous systems too. Comparatively, Riemann-Liouville and Caputo definitions are all generalization from the Grunwald-Letnikov one, and they are nearly equivalent on some conditions. The essential difference is that hypersingularity is required on initial condition in the sense of Riemann-Liouville definition, but only weak singularity, that is, integer-order derivative with engineering and physical significance on initial condition, is required in the sense of Caputo definition. Consider that, in the present paper, we take the Caputo definition forDαx(t), and only 0<α<1 is considered.(2)Dαxt=1Γ1-α∫0tx˙t-τταdτ,0<α<1,where Γ() is the Gamma function. f(x,x˙) is a function with respect to generalized displacement and generalized velocity; W(t) is wide-band stationary and ergodic random noise with spectral density:(3)Sω=Dπ1ω2-ωk22+4ζk2ωk2ω2.ζk and ωk are all constants. The samples of wide-band noise can be produced by the following second-order linear filter equation, in which G(t) is the Gaussian white noise with zero mean and constant intensity D.(4)W¨t+2ζkωkW˙t+ωk2Wt=Gt.In order to develop an analytical procedure to estimate the probability density function (PDF) for the system response, it is needed to introduce the generalized Van der Pol transformation from the system response (x,x˙) to envelope response a(t). According to the assumption on parameters, the sample motion of system (1) is quasiperiodic; hence, we can take(5)xt=atcos⁡φtx˙t=-atωasin⁡φtin which φ(t)=ω(a)t+γ(t), a(t) is the envelope process and γ(t) is the phase process; they are all stochastic processes slowly varying with respect to time t; ω(a) represents effective angular frequency of system (1). It has been illustrated [15] that function f(x,x˙) partly plays a role of classical damping force and partly contributes to the restoring force. Therefore, it is feasible to replace it by the combination of linear restoring force and linear damping; that is,(6)fx,x˙=2ω0ζeax˙+ωe2ax,where ζe(a) and ωe2(a) can be computed by applying Harmonic Balance Technique [16]; that is,(7)ζea=-12ωaω0πa∫02πfacos⁡φ,-aωasin⁡φsin⁡φdφ,(8)ωe2a=1πa∫02πfacos⁡φ,-aωasin⁡φcos⁡φdφ.Obviously, the effective angular frequency is something with natural frequency of system, ωe, and fractional derivative. Thus, the original nonlinear system (1) can be reconstructed as a quasilinear system, which is(9)x¨+2ω0ζ+εζeax˙+βDαxt+ω02+εωe2ax=ξtwith initial conditions x(0)=0 and x˙(0)=0. Since the envelope process a(t) varies slowly with respect to time, therefore a(t-τ)≈a(t) for all small delay τ. At the same time, first-order Taylor is expanding φ(t-τ) at the point of t; then the time-delay term in fractional derivative can be approximated by(10)x˙t-τ=-at-τωasin⁡φt-τ≈-atωasin⁡φcos⁡ωaτ-cos⁡φsin⁡ωaτ.As a consequence, the integral in the expression of fractional derivative Dαx(t) can be rewritten as(11)∫0tx˙t-τταdτ=atωacos⁡φ∫0tsin⁡ωaτταdτ-atωasin⁡φ∫0tcosωaτταdτ.The following asymptotic integrals are introduced herein [17] in order to solve the above integral, which will be very helpful to simplify the whole calculation for the expression of fractional derivative.(12)∫0tcos⁡ωτταdτ=ωα-1Γ1-αsin⁡πα2+sin⁡ωtωtα+o1ωtα+1,∫0tsin⁡ωτταdτ=ωα-1Γ1-αcos⁡πα2+cos⁡ωtωtα+o1ωtα+1,where o() denotes infinitesimal of the higher order. Substituting (12) into (11) and ignoring the higher-order terms, the fractional derivative can be replaced by(13)∫0tx˙t-τταdτ=atcos⁡φωαaΓ1-αcos⁡πα2-cos⁡ωtωtα-atsin⁡φωαaΓ1-αsinπα2+sin⁡ωtωtα.Substituting (5) and (13) into (9), and after some mathematical procedure, the equation governing the evolution of a(t) can be obtained(14)a˙=-2aω0ζ+εζesin2⁡φ-sin⁡φωaξt+βaωα-1asin⁡φcos⁡φcos⁡πα2-sin⁡φcos⁡πα2-βasin⁡φωacos⁡φcos⁡ωt-sin⁡φsin⁡ωtΓ1-αtα.


## 3. Procedure of Stochastic Averaging

Since the envelopea(t) varies slowly with respect to time, according to the Stratonovich-Khasminskii limit theory [18], the slow-varying process a(t) will converge to one-dimensional diffusive Markov equation in the sense of Gaussian white-noise excitation after deterministic averaging and stochastic averaging, called stochastic averaging method. Stochastic averaging method was initially proposed by Stratonovich in solving the response of stochastic dynamical systems; its basic idea is to approximate the system response by a Fokker-Planck-Kolmogorov equation, which is derived from a diffusive Markov equation with transient probability density governed by a partial differential equation; more relative investigations and mathematical foundation can be found in [18, 19].Suppose the system response is ergodic in the state space; then, stochastic averaging method leads the envelope processa(t) to the following first-order stochastic differential equation based on averaging operator t=limT→∞(1/T)∫0Tdt:(15)a˙=-aω0ζ+εζe-aβωα-1a2sinπα2+πSω2aω2a+πSω1/2ωaηt,where η(t) is a zero-mean and delta correlated random process and S(ω) is power spectral function defined in (3).Consequently, the Fokker-Planck-Kolmogorov (FPK) equation associated with transient solution of fractional dynamical system (1) can be casted as(16)∂tpa,t=-∂aMapa,t+12∂a2Napa,t,where the symbol p(a,t) is the probability density function of envelope process a(t). Furthermore, we assume that the system is initially at rest; this condition can be expressed as p(a,0)=δ(a), where δ(a) is Dirac delta function. Drift function M(a) and diffusion function N(a) in FPK (16) are, respectively, of the following form:(17)Ma=-aω0ζ+εζe-aβωα-1a2sinπα2+πSω2aω2a(18)Na=πSωωa.


## 4. Solution of the FPK Equation

Searching for solution of the FPK equation is always a hard work to get the response of the stochastic dynamical systems, and many efforts have been denoted to solving this problem. In some special cases such as first-order nonlinear systems with Gaussian excitations or second- and higher-order linear systems with constant coefficients, exact solution of FPK equation can be achieved. However, in most general cases, we have to explore approximate solutions for governing FPK equation of high-order and nonlinear systems. In this section, we will discuss the stationary solution in a closed form and transient solution in an approximate way for envelope responsea(t).


### 4.1. Stationary Solution

By reviewing the partial differential equation (16), stationary solution will be obtained if time tends to infinity; that means the stationary probability density function of envelope process a(t), which we denote by ps(a), can be found if ∂p(a,t)/∂t=0 in (16). Based on Spanos et al. work in [15], the stationary PDF ps(a) satisfies an ordinary differential equation, which is(19)-ddaMaps+12ddaNadpsda+psddaNa=0.The exact solution of this differential equation may be determined in a closed form [18, 19] as(20)psa=cNaexp⁡2∫MaNada,where c is normalization constant. Correspondingly, the exact stationary response of original system (1) associated with generalized displacement and velocity might derive from conversion back.(21)px,x˙=psaa=x2+x˙2/ω2a.


### 4.2. Transient Solution

In this section, the Galerkin method will be developed to estimate the transient solution of the fractional dynamical system (1) on the basis of PFK equation. It has been mentioned that Galerkin method is an important numerical analysis on the basis of variation of parameters. One can characterize the solution of a differential equation by a finite set of basis functions with some constraints on the function space. In this way, original differential equation can be converted to be a set of solvable high-dimensional linear equations in a weak sense. Following Galerkin method, the approximate transient solution denoted by p~(a,t) can be expressed by the following:(22)p~a,t=pLa,t+∑n=0NyntAna,where pL(a,t) is transient solution to the linear system when parameter ε=0, which is governed by(23)pLa,t=aσ21-exp⁡-2ζω0texp⁡-a22σ21-exp⁡-2ζω0t,in which σ2 is the stationary variance of x(t) for ε=0; the eigenvalues λn and eigenfunctions An(a) can be obtained from partial differential equation (16); they are given by(24)λn=2ζω0n,n∈Z+Ana=aσ2exp⁡-a22σ2Lna22σ2,n∈Z+,where Ln is the Laguerre polynomial of order n. It has been proved [18] that the eigenfunction An(a) fulfills the orthogonality condition on the basis of the properties of Laguerre polynomials and satisfies the following two recurrence formulae:(25)An,Am=∫0∞AnaAmaA0a=δnm,n,m∈Z+,where δnm is the Kronecker delta symbol. Except that, it can be proved that eigenfunction An(a) satisfies the following iterative formula and differential relationship:(26)n+1An+1a=2n+1-a22σ2Ana-nAn-1aaAn′a=2n+1An+1a-2n+1Ana.In (22), yn(t) is a series of unknown deterministic functions with respect to t, and N is an appropriate order, which will be decided by the approximation accuracy between the exact solution and approximated solution that we expected.Substitutingp~(a,t) into (16) and at the same time replacing the exact transient solution p(a,t) by the approximate one, that is, p~(a,t), consequently, the residual error will occur. We denote the residual error by Δr, which should be the function with regard to a(t) and yn(t); hence,(27)Δra,Yt=∂tp~a,t-∂tpa,t,where Y(t)=[y0(t),y1(t),…,yN(t)]T is a vector function composed of N+1 unknown functions. Obviously, Δr[a,Y(t)] should be zero if p~(a,t) is exactly equivalent to p(a,t) in the strong sense. However, p~(a,t) is usually not zero because p~(a,t) is only assumed to be an approximation of p(a,t), but not a definitely exact solution. Therefore, according to Galerkin technique, the N+1 unknown function yn(t) in (27) can be estimated by imposing the condition that the projection of the residual error Δr[a,Y(t)] on an appropriate set of independent weighing functions is zero. Taking orthogonality condition of An(a) into account, multiplying (27) at two sides byAj(a)/A0(a), and then integrating from zero to infinity and utilizing the orthogonality relation (25), we have(28)∫0∞AjaA0aΔra,Ytda=0,j∈Z+.This is a set of linear first-order ordinary differential equations containing N+1 unknown functions yn(t); we can simplify (28) by considering the orthogonality of eigenfunctions An(a). In this way, a reduced equation is derived, which can be expressed as a vector form(29)Y˙t=CYt+Ft,where C is a matrix of order (N+1)×(N+1) and F(t) is an N+1 dimensional vector; their elements can be obtained by solving (28). At this point, how to find the transient solution for the fractional system (1) turns out to be how to get unknown function yn(t) from solving this linear differential equation (29). Specifically, we can define the following operator for eigenfunction:(30)Im,n,k=∫0∞akAmaAnaA0admdn,n,k=0,1,2,…Then, it can be proved that the above integrals satisfy a recursive formula governed as follows by using the properties of Laguerre polynomials:(31)Im,n,k+2=2σ22n+1Im,n,k-nIm,n-1,k-n+1Im,n+1,k.With the assistance of these recursive formulas, substituting back yn(t) into (22), then the transient solution of envelope response can be completely obtained.


## 4.1. Stationary Solution

By reviewing the partial differential equation (16), stationary solution will be obtained if time tends to infinity; that means the stationary probability density function of envelope process a(t), which we denote by ps(a), can be found if ∂p(a,t)/∂t=0 in (16). Based on Spanos et al. work in [15], the stationary PDF ps(a) satisfies an ordinary differential equation, which is(19)-ddaMaps+12ddaNadpsda+psddaNa=0.The exact solution of this differential equation may be determined in a closed form [18, 19] as(20)psa=cNaexp⁡2∫MaNada,where c is normalization constant. Correspondingly, the exact stationary response of original system (1) associated with generalized displacement and velocity might derive from conversion back.(21)px,x˙=psaa=x2+x˙2/ω2a.


## 4.2. Transient Solution

In this section, the Galerkin method will be developed to estimate the transient solution of the fractional dynamical system (1) on the basis of PFK equation. It has been mentioned that Galerkin method is an important numerical analysis on the basis of variation of parameters. One can characterize the solution of a differential equation by a finite set of basis functions with some constraints on the function space. In this way, original differential equation can be converted to be a set of solvable high-dimensional linear equations in a weak sense. Following Galerkin method, the approximate transient solution denoted by p~(a,t) can be expressed by the following:(22)p~a,t=pLa,t+∑n=0NyntAna,where pL(a,t) is transient solution to the linear system when parameter ε=0, which is governed by(23)pLa,t=aσ21-exp⁡-2ζω0texp⁡-a22σ21-exp⁡-2ζω0t,in which σ2 is the stationary variance of x(t) for ε=0; the eigenvalues λn and eigenfunctions An(a) can be obtained from partial differential equation (16); they are given by(24)λn=2ζω0n,n∈Z+Ana=aσ2exp⁡-a22σ2Lna22σ2,n∈Z+,where Ln is the Laguerre polynomial of order n. It has been proved [18] that the eigenfunction An(a) fulfills the orthogonality condition on the basis of the properties of Laguerre polynomials and satisfies the following two recurrence formulae:(25)An,Am=∫0∞AnaAmaA0a=δnm,n,m∈Z+,where δnm is the Kronecker delta symbol. Except that, it can be proved that eigenfunction An(a) satisfies the following iterative formula and differential relationship:(26)n+1An+1a=2n+1-a22σ2Ana-nAn-1aaAn′a=2n+1An+1a-2n+1Ana.In (22), yn(t) is a series of unknown deterministic functions with respect to t, and N is an appropriate order, which will be decided by the approximation accuracy between the exact solution and approximated solution that we expected.Substitutingp~(a,t) into (16) and at the same time replacing the exact transient solution p(a,t) by the approximate one, that is, p~(a,t), consequently, the residual error will occur. We denote the residual error by Δr, which should be the function with regard to a(t) and yn(t); hence,(27)Δra,Yt=∂tp~a,t-∂tpa,t,where Y(t)=[y0(t),y1(t),…,yN(t)]T is a vector function composed of N+1 unknown functions. Obviously, Δr[a,Y(t)] should be zero if p~(a,t) is exactly equivalent to p(a,t) in the strong sense. However, p~(a,t) is usually not zero because p~(a,t) is only assumed to be an approximation of p(a,t), but not a definitely exact solution. Therefore, according to Galerkin technique, the N+1 unknown function yn(t) in (27) can be estimated by imposing the condition that the projection of the residual error Δr[a,Y(t)] on an appropriate set of independent weighing functions is zero. Taking orthogonality condition of An(a) into account, multiplying (27) at two sides byAj(a)/A0(a), and then integrating from zero to infinity and utilizing the orthogonality relation (25), we have(28)∫0∞AjaA0aΔra,Ytda=0,j∈Z+.This is a set of linear first-order ordinary differential equations containing N+1 unknown functions yn(t); we can simplify (28) by considering the orthogonality of eigenfunctions An(a). In this way, a reduced equation is derived, which can be expressed as a vector form(29)Y˙t=CYt+Ft,where C is a matrix of order (N+1)×(N+1) and F(t) is an N+1 dimensional vector; their elements can be obtained by solving (28). At this point, how to find the transient solution for the fractional system (1) turns out to be how to get unknown function yn(t) from solving this linear differential equation (29). Specifically, we can define the following operator for eigenfunction:(30)Im,n,k=∫0∞akAmaAnaA0admdn,n,k=0,1,2,…Then, it can be proved that the above integrals satisfy a recursive formula governed as follows by using the properties of Laguerre polynomials:(31)Im,n,k+2=2σ22n+1Im,n,k-nIm,n-1,k-n+1Im,n+1,k.With the assistance of these recursive formulas, substituting back yn(t) into (22), then the transient solution of envelope response can be completely obtained.


## 5. Examples and Results

As an example, we consider a Rayleigh oscillator subjected to wide-band noise excitation, in which(32)fx,x˙=c1x˙+c2x˙3,where c1 and c2 are all constants. Substituting this function into system (1), then system (1) can be represented as(33)x¨t+2ζω0x˙t+βDαxt+ω02xt+εc1x˙+c2x˙3=Wt.Based on formulas (7) and (8), it can be verified that(34)ζea=4c1+3c2ω2aa28ω0,ωea=0.If so, the exact stationary solution of the envelope a(t) is determined according to (20), which is governed by(35)psa=cω0πSωexp⁡-ω02ζω0+εc1/2+β/2ω0α-1sin⁡πα/2πSωa2-3ω04εc2a416πSω+ln⁡a.Relatively, the differential equation governing unknown functions yn(t) satisfies the following differential equations according to (29):(36)y˙jt+λj-2jc1ε2+βω0α-22sin⁡πα2yjt-c1ε2+βω0α-22sin⁡πα2∫0∞AjaA0apLda+∫0∞aAjaA0a∂apLda+∑n=0Nynt2n+2I0,j,n+1-3c2εω0283∫0∞a2AjaA0apLda+∫0∞a3AjaA0a∂apLda-3c2εω024∑n=0Nynt1-nI2,j,n+1+nI2,j,n+1=0,j=0,1,…,N.Figure1 shows the exact solutions of the system envelope according to the mathematical formula (35) and Monte Carlo Simulation. The solid lines represent the exact solution, and the hollow circles represent the simulated solution. In this example, we take parameters as ζ=0.04,ω0=2,ε=0.45,β=0.045,c1=-0.1,c2=0.001,D=0.81,ωk=2.1,ζk=0.05,andα=0.2. Observing that two results are in good agreement with each other, therefore, this figure illustrates the efficiency of stochastic averaging method and Galerkin method.Figure 1
Stationary solution of the envelope response obtained from exact expression and Monte Carlo Simulation.Joint probability density function with respect to generalized displacement and velocity is also plotted in Figure2. It is seen that the joint probability density function is an upside-down mountain shape with one peak, which shows the entire evolution principle of system solution. Except that, Figures 3 and 4 show the stationary displacement and velocity response of original system by calculating the marginal probability density; the good agreements are also achieved between exact solutions and MCS.Figure 2
Joint stationary PDF response of generalized displacement and velocity.Figure 3
Stationary generalized displacement solution of the fractional system.Figure 4
Stationary generalized velocity solution of the fractional system.Figure5 shows the solutions to the differential equation (36) for N=7, where eight unknown functions yn(t), n=0,1,2,…,7, are displayed with respect to time t, where parametric values are the same as those in Figure 1. It can be seen that y1(t) and y2(t) decrease from zero monotonously as time increases but others tend to be constant very fast, which means these two functions are more important than the others, because they nearly dominate the whole approximated series in (29) as large value of N is considered.Figure 5
The time evolution of unknown functionsyn(t) as time increase.Figure6 displays the transient solution p~(a,t) at different time instants by way of proposed Galerkin method and Monte Carlo Simulation, where we take N=15. Observe that the maximum of the transient solutions goes down as the increase of time, and the peak points move toward the right with the increase of envelope values. Particularly, the transient solution tends to be stationary very fast once time reaches 150 s and more. The comparisons between analytical solution and those derived from Monte Carlo Simulation have demonstrated that the methods we used in this paper are in good agreement at each time instant. In fact, Monte Carlo Simulation is one of the powerful tools to find statistical properties of stochastic dynamical systems on the basis of probability theory. Generally, the first step is to generate amounts of samples according to the distribution of noise excitation and then solve the differential equation by Rugge-Kutta method from one initial point. Repeat the procedure above to get a great large amount of samples; at last, we can get the probability density of envelope response at any time based on the law of large numbers.Figure 6
The transient solution of the system at different time instant along with the envelope.At the given time instantt=150, Figure 7 illustrates that the solutions we get from the proposed Galerkin method also coincide with both the exact solution and MCS; this agreement exactly proves the efficiency of our proposed methods. Moreover, the effect caused by the terms of truncated series is plotted and shown in Figure 8. Note that 5 and higher terms in truncated series are enough to get an expected result of approximate transient PDF solutions.Figure 7
The stationary solution of the fractional system obtained by three different approaches.Figure 8
The influence caused by the orderN in approximated series.


## 6. Conclusion

The transient solution of a dynamical system with fractional derivative term is discussed from the point of view of analytical and numerical points. We firstly apply generalized Harmonic Balance technique to transform the system to be a quasilinear one. Based on the framework of stochastic averaging method, the fractional derivative damping is averaged to be a function associated with envelope. In this way, the system envelope can be described by a Markov diffusive equation with an FPK equation followed. Galerkin method is used to get time-evolutionary probability density function for envelope response. Numerical results concerning Rayleigh oscillator including stationary solution and transient solution have shown the efficiency of the proposed approaches and agreement with Monte Carlo Simulation.



---

*Source: 9638523-2016-11-23.xml*