Abstract
A modified stochastic averaging method on a DuffingRayleigh oscillator with strongly nonlinearity subject to Gaussian colored noise excitation was proposed. The socalled He’s energy balance method was applied to obtain the averaged frequency of the conservative system. Subsequently, the stochastic averaging method of strong nonlinearity was used. The modified method can offer more concise approximate expressions of the drift and diffusion coefficients without weakening the accuracy of predicting the responses of the systems too much. The stationary responses of probability density of amplitudes, together with joint probability density of displacement and velocity are studied to verify the presented approach. The reliability of the systems was also investigated to offer further support. Digital simulations were carried out and the output of that are coincide with the theoretical approximations well.
1. Introduction
In recent half century, stochastic vibration has drawn accelerating interests. An effective and powerful method dealing with stochastic vibrations is the stochastic averaging method. The stochastic averaging method initialized by Landau and Stratonovich [1] and Khasminskii [2] is also called as the method of quasiconservative averaging. Later, Zhu and Lin [3, 4] rederived an approach named the stochastic averaging of energy envelope by taking into account of the conservative part of the WongZakai correction terms of the system. This prevalent approach has been used and developed by many researchers [512]. The responses of the singledegreeof freedom (SDOF) or multidegreeoffreedom (MDOF) systems subject to nonwhite wideband random excitation [68] or combined harmonic and white noise [10] have all been studied in detail. After that, Zhu [11, 12] extended the stochastic averaging of energy envelope method into Hamilton systems.
Xu and Chung [13] have proposed a method to deal with noisefree strongly nonlinear Oscillators based on the socalled generalized harmonic functions. This method has been extended by Zhu and Huang [14, 15] to handle strongly nonlinear oscillators with lightly linear and (or) nonlinear damping subject to weakly external and (or) parametric excitations of wideband random processes. The key technique of the Zhu and Huang’s method is using expressions of the displacement $x$ and velocity $y$ as: $x=A\mathrm{c}\mathrm{o}\mathrm{s}\theta \left(t\right)$, $y=\dot{x}=A\nu (A,\theta )\mathrm{s}\mathrm{i}\mathrm{n}\theta \left(t\right)$, where $\nu (A,\theta )$ is the instantaneous frequency and $A$ is the amplitude. The frequency $\nu (A,\theta )$ can be averaged to $\omega \left(A\right)$ by expending it into a string of Fourier series and integrating them with respect to $\theta $ form 0 to $2\pi $. A typical DuffingVanderpol oscillator has been studied by this method [15], the stability of the stationary response and the stochastic Hopfbifurcation has been investigated. The precision of this method has been verified by digital simulations. However, this method is less than flawless in view of the fact that it yields very long expressions of drift coefficient and diffusion coefficient which are quite time consuming during calculations. Thus, there remains a need to shorten the expressions of the drift coefficient and diffusion coefficient without weakening the accuracy too much. Furthermore, this method is mainly applied on white noise or nonwhite wideband random excitations, we try to extend this method to oscillators under colored noise excitations.
In this paper, a DuffingRayleigh oscillator under colored noise excitation was studied as a representative example. A socalled energy balance method calculating the averaged frequency of conservative systems proposed by He [16] has been used to obtain concise expression of the averaged $\omega \left(A\right)$. The feature of the method is constructing the Hamilton of the system and keeping the input and the output of it in balance. This method has been extended by many researches into nonconservative systems [1721] and the accuracy has been proved to be good enough. In the following sections, The Zhu and Huang’s method was modified by the He’s energy balance method. The stationary probability density function (PDF) of amplitude and the joint probability of the displacement and velocity were given simultaneously based on the shortened frequency $\omega \left(A\right)$. The reliability and the first passage failure were also studied to offer further support. Digital simulations were carried out and the results were consistent with the theoretical approximations.
2. Stochastic averaging method
A DuffingRayleigh oscillator with light damping and colored noise excitation was presented as follows. The linear and nonlinear damping coefficients are assumed to be of ${\epsilon}^{2}$ order, $\mu \to {\epsilon}^{2}\mu $, $\gamma \to {\epsilon}^{2}\gamma $, the excitation $\eta \left(t\right)$ is assumed to be of ${\epsilon}^{1}$ order:
As well known, the colored noise $\eta \left(t\right)$ can be obtained by a white noise passing through a linear filter:
where $\lambda $ denotes the correlation time and $\xi \left(t\right)$ denotes a Gaussian white noise with 0 mean and intensity of 2$D$. The correlation function of $\eta \left(t\right)$ is:
The conservative system of Eq. (1) is:
The Hamilton and the potential energy are given as:
where $H$ is the total energy of the oscillator while $G\left(x\right)$ is the potential energy.
Assuming all the trajectories in domain $U$ of phase plan ($x,\dot{x}$) are periodic surrounding (0, 0). The approximate periodic solution of Eq. (1) can be written as:
where:
in which $A$ is the amplitude of displacement, which is related to $H$ as follows:
and $\nu (A,t)$ is the instantaneous frequency of the oscillation. For the convenience of calculating, the instantaneous frequency $\nu (A,t)$ can be transformed to be the averaged frequency $\omega \left(A\right)$ by the He’s energy balance method [16].
For $\theta \left(t\right)=0$, we get $x\left(0\right)=A$, $\dot{x}\left(0\right)=0$, substituting it into Eq. (7) yields the Hamilton ${H}_{\theta =0}=\frac{1}{2}{\omega}^{2}{A}^{2}+\frac{1}{4}\alpha {A}^{4}$. Then, replacing $\nu (A,t)$ by $\omega \left(A\right)$ in Eq. (8) and substituting it into Eq. (7) yields the Hamilton ${H}_{\theta =\omega t}$. we obtain the following residual equation:
Since Eq. (7) is only an approximation to the exact solution, $R$ cannot be made zero everywhere. Collocation at $\omega t=\pi /4$, we get:
By letting $R\left(t\right)={H}_{\theta =\frac{\pi}{4}}{H}_{\theta =0}$ equals to zero, we can get the averaged frequency $\omega \left(A\right)$:
For Eq. (1) we can use the following approximate relation:
where $\phi \left(t\right)$ is the phase.
With the transformations accomplished, Eq. (1) becomes:
where:
${m}_{2}=\frac{1}{g\left(A\right)}[\mu A\omega \left(A\right)\mathrm{s}\mathrm{i}\mathrm{n}\theta +\gamma \omega \left(A{)}^{3}{A}^{3}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \right]\omega \left(A\right)\mathrm{c}\mathrm{o}\mathrm{s}\theta ,$
${b}_{1}=\frac{A}{g\left(A\right)}\omega \left(A\right)\mathrm{s}\mathrm{i}\mathrm{n}\theta ,{b}_{2}=\frac{1}{g\left(A\right)}\omega \left(A\right)\mathrm{c}\mathrm{o}\mathrm{s}\theta .$
Eq. (14) is of standard form. $A\left(t\right)$ converges weakly to a diffusive Markov process even in a semiinfinite time interval. The Ito equation of the limiting diffusion process is of the form:
where the drift coefficient $m\left(A\right)$ and the diffusion coefficient $\sigma \left(A\right)$ are as follows:
${\sigma}^{2}\left(A\right)=\u27e8{\int}_{\mathrm{\infty}}^{+\mathrm{\infty}}\left({\left.{b}_{1i}\right}_{t}{\left.{b}_{1j}\right}_{t+\tau}\right){R}_{ij}\left(\tau \right)d\tau \u27e9.$
$W\left(t\right)$ is the standard unit Brown's motion, $\u27e8\cdot \u27e9$ is averaging with respect to time, ${R}_{ij}\left(\tau \right)$ is the correlation function. For white noise, it is well known that:
we can obtain the following drift and diffusion coefficients of stochastic equations for amplitude $A\left(t\right)$:
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+\frac{D{\lambda}_{}^{2}\left({A}^{2}\alpha {\omega}_{0}^{2}+2{\omega}_{0}^{4}\right)}{A{\left({A}^{2}\alpha +{\omega}_{0}^{2}\right)}^{3}\left(3{A}^{2}\alpha +4{\lambda}_{}^{2}+4{\omega}_{0}^{2}\right)},$
${\sigma}^{2}\left(A\right)=\frac{{D}_{1}{\lambda}_{1}^{2}\left(3{A}^{2}\alpha +4{\omega}_{0}^{2}\right)}{{\left({A}^{2}\alpha +{\omega}_{0}^{2}\right)}^{2}\left(3{A}^{2}\alpha +4{\lambda}_{1}^{2}+4{\omega}_{0}^{2}\right)}.$
3. Stationary probability density function
The averaged FokkerPlankKolmogorov (FPK) equation associated with Ito Eq. (15) is of the form:
where $p(A,t)$ is the transition probability density of displacement amplitude $A$ nontrivial stationary solution exists due to nonvanishing of external excitation. Under assumption of zero probability flow at the two boundaries $A=0$ and $A=\mathrm{\infty}$, the stationary solution of FPK Eq. (19) for system Eq. (15) is of the form:
where constant $N$ is the normalization constant.
Substituting Eq. (18) into Eq. (20) yields:
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\bullet \mathrm{e}\mathrm{x}\mathrm{p}\left[\frac{27{A}^{10}{\alpha}^{3}\gamma}{640D{\lambda}_{}^{2}}+{A}^{8}\left(\frac{9{\alpha}^{2}\gamma}{128D}\frac{99{\alpha}^{2}\gamma {\omega}_{0}^{2}}{512D{\lambda}_{}^{2}}\right)+{A}^{6}\left(\frac{{\alpha}^{2}\mu}{8D{\lambda}_{}^{2}}\frac{7\alpha \gamma {\omega}_{0}^{2}}{32D}\frac{5\alpha \gamma {\omega}_{0}^{4}}{16D{\lambda}_{}^{2}}\right)\right.$
$\left.\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+{A}^{2}\left(\frac{\mu {\omega}_{0}^{2}}{2D}\frac{\mu {\omega}_{0}^{4}}{2D{\lambda}^{2}}\right)+{A}^{4}\left(\frac{\alpha \mu}{4D}\frac{7\alpha \mu {\omega}_{0}^{2}}{16D{\lambda}_{1}^{2}}\frac{3\gamma {\omega}_{0}^{4}}{16D}\frac{3\gamma {\omega}_{0}^{6}}{16D{\lambda}_{}^{2}}\right)\right].$
The stationary probability density of total energy $H$ or energy envelope can be obtained as follows:
where:
is the inverse function of $G\left(A\right)=H$.
The stationary probability density of displacement $x$ and velocity $y$ can be further obtained from $p\left(H\right)$ as follows:
where $T\left(H\right)$ can be obtained from $T\left(A\right)=2\pi /\omega \left(A\right)$ by the transformation $A={G}^{1}\left(H\right)$.
Fig. 1Stationary probability density function of amplitude A (solid lines: stochastic averaging Eq. (21); Dashed lines: method in reference [15]; dots: digital simulation)
a)
b)
Comparisons of theoretical analysis and digital simulation were displayed in Fig. 1 and Fig. 2. The damping coefficients were set as $\mu =$ 0.1, $\gamma =$ 0.1, the nature frequency was set as $\omega =$ 1, and the nonlinear stiffness coefficient was set as $\alpha =\text{1}$. The densities of white noise were set as $D=$0.5 and $D=$1.0, and the correlation time was set as $\lambda =\text{2}$ and $\lambda =$ 10, respectively. Totally 1 000 sets of white noise were imposed to the system Eq. (1). Each set noise includes 20 000 numbers. The numerical solutions of the oscillator Eq. (1) could be obtained by an order2 stochastic RungeKutta algorithm [22]. The last 10 000 dots were kept as the steady state responses. The stationary probability density function (PDF) of the amplitude with different $D$ and same $\lambda $ was shown in Fig. 1(a). As we know, large noise density leads to large horizontal coordinate $A$ of the peak of PDF and low peak value $p\left(A\right)$. Also, the stationary probability density function of the amplitude with different $\lambda $ and same $D$ was illustrated in Fig. 1(b). The solid lines were results obtained by the proposed method in this paper and the dashed lines were given by the method in reference [15]. We have to admit that the method in reference [15] offered better accuracy than the method we used in this paper. But, it was also obvious that the presented method can offer acceptable analytical approximate solutions in view of the fact the solid lines were quite against the dashed lines except for the peak positions. Furthermore, considering the reduction of the calculation and the error between the dashed line and the digital results was negligible, the presented method is valid and effective. It is obvious that larger time delay coefficient $\lambda $ results in wider diffusion range of amplitude.
The joint probability of the displacement $x$ and velocity $y$ is shown is Fig. 2, where Fig. 2(a) is the theoretical result given by Eq. (23) and the Fig. 2(b) is the digital simulation result, with parameters set as $D=$0.5 and $\lambda =\text{2}$. It is seen that the analytical results are fairly in agreement with those from digital simulation.
Fig. 2The joint probability of the displacement and velocity: a) theoretical result, b) digital simulation result
a)
b)
To give a intuitive impression of the system excited by colored noise, three samples of the phase trajectory were illustrated in the Fig. 3. It is obvious that both the large noise density and large correlation time lead to thickly cycles of trajectory spreading wider. This phenomenon can be seen clearly in Fig. 1. The peaks of the stationary probability density functions were drifted to the right side of the coordinate with noise density and large correlation time increased.
Fig. 3Samples of Phase trajectories: a) D= 0.5, λ= 2, b) D= 1, λ= 2, c) D= 0.5, λ= 10
a)
b)
c)
If the linear damping is increased ten times larger than the formal value $\mu =\text{0.1}$ i.e. $\mu =\text{1}$, and the other parameters are remained unchanged, the stationary response of amplitude is given as Fig. 4. The difference between the analytical result and the digital simulation cannot be neglected It is obvious that our method will lose much accuracy. Because this stochastic averaging method is only effective on weak damping systems.
Fig. 4Stationary probability density function of amplitude A (solid lines: stochastic averaging Eq. (21); Dasheddotted lines: digital simulation) with D= 0.5, λ= 2, μ= 1
4. Reliability function and the probability of first passage failure time
In this section, the reliability function and the probability of first passage failure time were studied to predict the responses of the system. And this procedure can offer further support to our previous study on the drift coefficient and the diffusion coefficient. The conditional reliability function, denoted by $R\left(\left.T\right{A}_{0}\right)$ is defined as the probability of amplitude $A$ could endure within the safe domain $\mathrm{\Omega}=(0,{A}_{c})$ with the time increasing form $t=0$ to $t=T$ given an initial amplitude ${A}_{0}$. Reliability function is mathematically defined as:
The conditional transition probability density is governed by the backward Kolmogorov (BK) equation with drift coefficient $m\left(A\right)$ and diffusion coefficient $\sigma \left(A\right)$ defined by Eq. (18). It can be shown that the conditional reliability function satisfies the following BK equation:
The initial condition and the boundary conditions associated with Eq. (25) are:
It is reasonable to assume that the firstpassage occurs once $A$ exceeds certain critical value $A={A}_{c}$ for the first time. Obviously, conditional damage probability can be defined as following:
The conditional probability density of first passage time can be obtained from the conditional reliability function as follows:
The critical value of the amplitude ${A}_{c}$ was set as 0.8, other parameters were set the same as those in the last section. The digital simulations were carried out, and the results were detailed in Fig. 5. The solid lines denote the theoretical solutions and the dots denote the digital simulations. Fig. 5(a) showed the reliability and Fig. 5(b) showed the probability density of the first passage failure with density $D=$0.2, $D=$0.4 and $\lambda =\text{0.5}$. It is clear that the reliability decreases with time and the larger the noise density $D$ is, the faster the reliability decreases. It is also obvious that larger noise density yields higher probability density of first passage failure. The influence of the time delay coefficient $\lambda $ on the reliability was illustrated in Fig. 5(c) and Fig. 5(d). With noise density fixed as $D=$0.5, time delay parameter $\lambda $ was chosen as $\lambda =$ 0.6 and $\lambda =$ 1.0 respectively. The Fig. 5(c) and Fig. 5(d) showed that larger time delay parameter $\lambda $ will result in faster decrease velocity of reliability and higher probability density of the first passage failure time. The numerical simulations are consistent with the theoretical prediction well enough, which means the drift coefficient and diffusion coefficient used in Eq. (25) are correct. Consequently, it is reasonable to believe that the stochastic averaging method we modified in this paper is accurate in calculating the drift coefficient and diffusion coefficient. The expression of the drift coefficient and diffusion coefficient is not very complex and the accuracy is acceptable.
Fig. 5The reliability function and the probability of fist passage failure time (solid lines: theoretical solutions; dots: digital simulations)
a)
b)
c)
d)
5. Conclusions
A modified stochastic averaging method dealing with colored noise excited oscillators is proposed in this paper, with the help of the He’s energy balance method the expression of the drift coefficient and diffusion coefficient are shorten significantly without weakening the accuracy too much. The amount of calculating has been reduced consequently. The DuffingRayleigh oscillator subject to colored noise excitation was taken as an example to verify the method. The probability of stationary responses of the amplitude and the joint probability of displacement and velocity are all studied. Furthermore, the reliability and the first passage failure are also studied. All the analytical approximations are consistent with the digital simulations.
References

Landau P. S., Stratonovich R. L. Theory of stochastic transitions of various systems between different states. Proceedings of Moscow University, Vol. 3, 1962, p. 3345, (in Russian).

Khasminskii R. Z. On the behavior of a conservative system with small friction and small random noise. Journal of Applied Mathematics and Mechanics, Vol. 28, 1964, p. 11261130, (in Russian).

Zhu W. Q. Stochastic averaging of the energy envelope of nearly Lyapunov systems. Proceedings of the IUTAM Symposium on Random Vibrations and Reliability, Akademie, Berlin, 1983, p. 347357.

Zhu W. Q., Lin Y. K. Stochastic averaging of energy envelope. Journal of Engineering Mechanics, Vol. 117, Issue 8, 1991, p. 18901905.

RedHorse J. R., Spanos P. D. A generalization to stochastic averaging in random vibration. International Journal of NonLinear Mechanics, Vol. 27, Issue 1, 1992, p. 85101.

Cai G. Q. Random vibration of nonlinear systems under nonwhite excitations. Journal of Engineering Mechanics, Vol. 121, Issue 5, 1995, p. 637639.

Dementberg M., Cai G. Q., Lin Y. K. Application of quasiconservative averaging to a nonlinear system under nonwhite excitation. International Journal of NonLinear Mechanics, Vol. 30, Issue 5, 1995, p. 677685.

Cai G. Q., Lin Y. K., Xu W. Strongly Nonlinear System under NonWhite Random Excitations. Book Series, Stochastic Structural Dynamics, Balkema, Rotterdam, 1999, p. 1116.

Krenk S., Roberts J. B. Local similarity in nonlinear random vibration. Journal of Applied Mechanics, Vol. 66, Issue 1, 1999, p. 225235.

Huang Z. L., Zhu W. Q. Stochastic averaging of quasiintegrable Hamiltonian systems under combined harmonic and white noise excitations. International Journal of NonLinear Mechanics, Vol. 39, Issue 9, 2004, p. 14211434.

Zhu W. Q., Yang Y. Q. Stochastic averaging of quasinonintegrableHamiltonian systems. Journal of Applied Mechanics, Vol. 64, Issue 1, 1997, p. 157164.

Zhu W. Q., Huang Z. L., Yang Y. Q. Stochastic averaging of quasiintegrableHamiltonian systems. Journal of Applied Mechanics, Vol. 64, Issue 4, 1997, p. 975984.

Xu Z., Chung Y. K. Averaging method using generalized harmonic functions for strongly nonlinear oscillators. Journal of Sound and Vibration, Vol. 174, Issue 4, 1994, p. 563576.

Huang Z. L., Zhu W. Q. Stochastic averaging of strongly nonlinear oscillators under combined harmonic and whitenoise excitations. Journal of Sound and Vibration, Vol. 238, Issue 2, 2000, p. 233256.

Zhu W. Q., Huang Z. L., Suzuki Y. Response and stability of strongly nonlinear oscillators under wideband random excitation. International Journal of NonLinear Mechanics, Vol. 36, Issue 8, 2001, p. 12351250.

He J. H. Preliminary report on the energy balance for nonlinear oscillations. Mechanics Research Communications, Vol. 29, Issues 23, 2002, p. 107111.

Davodi A. G., Ganji D. D., Azami R., et al. Application of improved amplitudefrequency formulation to nonlinear differential equation of motion equations. Modern Physics Letters. Vol. 23, Issue 28, 2009, p. 34273436.

Ganji S. S., Ganji D. D., Karimpour S. He’s energy balance and He’s variation methods for nonlinear oscillations in engineering. International Journal of Modern Physics B, Vol. 23, Issue 3, 2009, p. 461471.

Mehdipour I., Ganji D. D., Mozaffari M. Applications of the energy balance method to nonlinear vibrating equations. Current Applied Physics, Vol. 10, Issue 1, 2010, p. 104112.

HuiLi Zhang. Periodic solutions for some strongly nonlinear oscillations by He’s energy balance method. Computers and Mathematics with Applications, Vol. 58, Issues 1112, 2009, p. 24802485.

Vounesian Davood, Askari Hassan, Saadatnia Zia, et al. Frequency analysis of strongly nonlinear generalized Duffing oscillators using He’s frequencyamplitude formulation and He’s energy balance method. Computers and Mathematics with Applications, Vol. 59, Issue 9, 2010, p. 32223228.

Honeycutt Rebecca L. Stochastic RungeKutta algorithms. Part 2: colored noise. Physical Review A, Vol. 45, Issue 2, 1992, p. 604610.
About this article
The authors gratefully acknowledge the support of the Natural Science Foundation of China (NSFC) through Grant Nos. 11402186 and 11302144 and the Tianjin Research Program of Application Foundation and Advanced Technology through Grant Nos. 14JCQNJC05600, 13JCYBJC17900, 14JCQNJC05300.