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.
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 Section 2. 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 formwhere is generalized displacement and the superscript dot represents differentiation with respect to time . , , 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 for , and only is considered.where is the Gamma function. is a function with respect to generalized displacement and generalized velocity; is wide-band stationary and ergodic random noise with spectral density: and are all constants. The samples of wide-band noise can be produced by the following second-order linear filter equation, in which is the Gaussian white noise with zero mean and constant intensity .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 to envelope response . According to the assumption on parameters, the sample motion of system (1) is quasiperiodic; hence, we can takein which , is the envelope process and is the phase process; they are all stochastic processes slowly varying with respect to time ; represents effective angular frequency of system (1). It has been illustrated [15] that function 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,where and can be computed by applying Harmonic Balance Technique [16]; that is,Obviously, the effective angular frequency is something with natural frequency of system, , and fractional derivative. Thus, the original nonlinear system (1) can be reconstructed as a quasilinear system, which iswith initial conditions and . Since the envelope process varies slowly with respect to time, therefore for all small delay . At the same time, first-order Taylor is expanding at the point of ; then the time-delay term in fractional derivative can be approximated byAs a consequence, the integral in the expression of fractional derivative can be rewritten asThe 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.where denotes infinitesimal of the higher order. Substituting (12) into (11) and ignoring the higher-order terms, the fractional derivative can be replaced bySubstituting (5) and (13) into (9), and after some mathematical procedure, the equation governing the evolution of can be obtained
3. Procedure of Stochastic Averaging
Since the envelope varies slowly with respect to time, according to the Stratonovich-Khasminskii limit theory [18], the slow-varying process 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 process to the following first-order stochastic differential equation based on averaging operator :where is a zero-mean and delta correlated random process and 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 aswhere the symbol is the probability density function of envelope process . Furthermore, we assume that the system is initially at rest; this condition can be expressed as , where is Dirac delta function. Drift function and diffusion function in FPK (16) are, respectively, of the following form:
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 response .
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 , which we denote by , can be found if in (16). Based on Spanos et al. work in [15], the stationary PDF satisfies an ordinary differential equation, which is
The exact solution of this differential equation may be determined in a closed form [18, 19] aswhere is normalization constant. Correspondingly, the exact stationary response of original system (1) associated with generalized displacement and velocity might derive from conversion back.
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 can be expressed by the following:where is transient solution to the linear system when parameter , which is governed byin which is the stationary variance of for ; the eigenvalues and eigenfunctions can be obtained from partial differential equation (16); they are given bywhere is the Laguerre polynomial of order . It has been proved [18] that the eigenfunction fulfills the orthogonality condition on the basis of the properties of Laguerre polynomials and satisfies the following two recurrence formulae:where is the Kronecker delta symbol. Except that, it can be proved that eigenfunction satisfies the following iterative formula and differential relationship:
In (22), is a series of unknown deterministic functions with respect to , and is an appropriate order, which will be decided by the approximation accuracy between the exact solution and approximated solution that we expected.
Substituting into (16) and at the same time replacing the exact transient solution by the approximate one, that is, , consequently, the residual error will occur. We denote the residual error by , which should be the function with regard to and ; hence,where is a vector function composed of unknown functions. Obviously, should be zero if is exactly equivalent to in the strong sense. However, is usually not zero because is only assumed to be an approximation of , but not a definitely exact solution. Therefore, according to Galerkin technique, the unknown function in (27) can be estimated by imposing the condition that the projection of the residual error on an appropriate set of independent weighing functions is zero. Taking orthogonality condition of into account, multiplying (27) at two sides by, and then integrating from zero to infinity and utilizing the orthogonality relation (25), we haveThis is a set of linear first-order ordinary differential equations containing unknown functions ; we can simplify (28) by considering the orthogonality of eigenfunctions . In this way, a reduced equation is derived, which can be expressed as a vector formwhere is a matrix of order and is an 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 from solving this linear differential equation (29). Specifically, we can define the following operator for eigenfunction:Then, it can be proved that the above integrals satisfy a recursive formula governed as follows by using the properties of Laguerre polynomials:With the assistance of these recursive formulas, substituting back 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 whichwhere and are all constants. Substituting this function into system (1), then system (1) can be represented asBased on formulas (7) and (8), it can be verified thatIf so, the exact stationary solution of the envelope is determined according to (20), which is governed byRelatively, the differential equation governing unknown functions satisfies the following differential equations according to (29):
Figure 1 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 Observing that two results are in good agreement with each other, therefore, this figure illustrates the efficiency of stochastic averaging method and Galerkin method.
Joint probability density function with respect to generalized displacement and velocity is also plotted in Figure 2. 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 5 shows the solutions to the differential equation (36) for , where eight unknown functions , , are displayed with respect to time , where parametric values are the same as those in Figure 1. It can be seen that and 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 is considered.
Figure 6 displays the transient solution at different time instants by way of proposed Galerkin method and Monte Carlo Simulation, where we take . 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.
At the given time instant , 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.
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.
Competing Interests
The author declares that they have no competing interests.
Acknowledgments
The work reported in this paper is supported by Natural Science Basic Research Plan in Shaanxi Province of China (no. 2015JM1028), the National Natural Science Foundation of China (nos. 11302157, 11202155), and National Natural Science Foundation of Fujian Province (no. 2014J014).