A stochastic averaging method on the strongly nonlinear Duffing-Rayleigh oscillator under Gaussian colored noise excitation

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 quasi-conservative averaging. Later, Zhu and Lin [3, 4] re-derived an approach named the stochastic averaging of energy envelope by taking into account of the conservative part of the Wong-Zakai correction terms of the system. This prevalent approach has been used and developed by many researchers [5-12]. The responses of the single-degree-of freedom (SDOF) or multi-degree-of-freedom (MDOF) systems subject to non-white wide-band random excitation [6-8] 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 noise-free strongly nonlinear Oscillators based on the so-called generalized harmonic functions. This method has been extended by Zhu and Huang [14, 15] to handle strongly non-linear oscillators with lightly linear and (or) non-linear damping subject to weakly external and (or) parametric excitations of wide-band 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 Duffing-Vanderpol oscillator has been studied by this method [15], the stability of the stationary response and the stochastic Hopf-bifurcation 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 non-white wide-band random excitations, we try to extend this method to oscillators under colored noise excitations.

In this paper, a Duffing-Rayleigh oscillator under colored noise excitation was studied as a representative example. A so-called 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 non-conservative systems [17-21] 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 Duffing-Rayleigh 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:

Eq. (14) is of standard form. $A\left(t\right)$ converges weakly to a diffusive Markov process even in a semi-infinite 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:

$W\left(t\right)$ is the standard unit Brown's motion, $⟨\cdot ⟩$ 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)$:

3. Stationary probability density function

The averaged Fokker-Plank-Kolmogorov (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 non-vanishing 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:

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 order-2 stochastic Runge-Kutta 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); Dashed-dotted 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 first-passage 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 Duffing-Rayleigh 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.