- Email: [email protected]
Study of the Duffing–Rayleigh oscillator subject to harmonic and stochastic excitations by path integration
Applied Mathematics and Computation 172 (2006) 1212–1224
www.elsevier.com/locate/amc
Study of the Duffing–Rayleigh oscillator subject to harmonic and stochastic excitations by path integration W.X. Xie a, W. Xu
a,*
, L. Cai
b
a
b
Department of Applied Mathematics, Northwestern Polytechnical University, Xi’an, Shaanxi 710072, PR China College of Astronautics, Northwestern Polytechnical University, Xi’an, Shaanxi 710072, PR China
Abstract In this paper, the Duffing–Rayleigh oscillator subject to harmonic and stochastic excitations is investigated via path integration based on the Gauss–Legendre integration formula. The method can successfully capture the steady state periodic solution of probability density function. This path integration method, using the periodicity of the coefficient of associated Fokker–Planck–Kolmogorov equation, is extended to deal with the averaged stationary probability density, and is efficient to computation. Meanwhile, the changes of probability density caused by the intensities of harmonic and stochastic excitations, are discussed in three cases through the instantaneous probability density and the averaged stationary probability density. 2005 Elsevier Inc. All rights reserved. Keywords: Path integration; Duffing–Rayleigh oscillator; Harmonic and stochastic excitations; Probability density
*
Corresponding author. E-mail address: [email protected] (W. Xu).
0096-3003/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.03.018
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
1213
1. Introduction So far, it is still hard to solve analytically the Fokker–Planck–Kolmogorov (FPK) equations for non-linear stochastic dynamical systems, which govern the transition probability density function (PDF) of the responses. Exact solutions are known just for some special systems (see [1–3] etc.). For this reason, several approximate and numerical techniques have been developed (see e.g. [4–7]). For the dynamical system subject to harmonic and stochastic excitations, one of the drift coefficients of the corresponding FPK equation is periodic. So there is no stationary solution and may exist steady state solution when t ! 1 [7]. Some approximate solutions for such systems using stochastic averaging method are given in Von Wagner and Wedig [4] and Huang et al. [7]. In fact, the method presented by Von Wagner and Wedig [4] can only pursue the approximate averaged stationary PDF of Duffing–Rayleigh oscillator subject to harmonic and white noise excitations by orthogonal function expansions, e.g. especially adjusted polynomials and Fourier series. The weighted functions of the new polynomials are obtained by the application of the stochastic averaging method. However, the calculation of double crater-like density for Duffing system is hard to perform. Huang et al. [7] get the averaging FPK equation by stochastic averaging method, and resort to the numerical path integration at last. Generally speaking, the approximate methods above have no way to capture the periodic changes of the instantaneous PDF for the response before averaging. But the numerical path integration method based on Gauss– Legendre integration formula proposed by Yu et al. [8,9], is a powerful tool to capture the evolution of the instantaneous PDF. Especially, Yu finds the evolution of the PDF to the periodic steady state PDF of Duffing system [9]. In this paper, a complicated system, the Duffing–Rayleigh oscillator subject to harmonic and stochastic excitations is of interest. The system has abundant dynamical behaviors, such as multi-value of response and hopping between stable responses. These phenomena were observed first in the Duffing–Rayleigh oscillator with multiple scales and second-order Gaussian closure method [10]. Taking advantages of this method developed by Yu and Lin [9], this paper not only obtains the periodic steady state PDF and the noise-induced state changes of PDF, but also can calculate the averaged stationary PDF by extending the path integration based on Gauss–Legendre integration formula, which is helpful to interpret the multi-value of response.
2. Duffing–Rayleigh oscillator subject to harmonic and stochastic excitations Consider the non-linear oscillator subject to harmonic and white noise excitations of the form
1214
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
b 2 X€ t þ c 1 þ 2 X_ t X_ t þ x21 ð1 þ eX 2t ÞX t ¼ x1 r1 sin xt þ x1 r2 nt ; x1
ð1Þ
where x1 is the natural circle frequency, c is the linear damping coefficient, b and e are the non-linearity parameters, r1 and x are intensity and frequency of the harmonic excitation respectively, and r2 is the intensity of the additive white noise nt. The corresponding FPK equation can be derived as oQ o o 1 o2 ¼ ½m1 ðx1 ; x2 ; tÞQ ½m2 ðx1 ; x2 ; tÞQ þ ½bðx1 ; x2 ; tÞQ; ot ox1 ox2 2 ox22
ð2Þ
where x1 ¼ x; x2 ¼ x_ ; Q ¼ qðx1 ; x2 ; tx0 ; x0 ; t0 Þ; 1
2
m1 ðx1 ; x2 ; tÞ ¼ x2 ;
m2 ðx1 ; x2 ; tÞ ¼ c 1 þ xb2 x22 x2 x21 ð1 þ ex21 Þx1 þ x1 r1 sin xt; 1
and
2
bðx1 ; x2 ; tÞ ¼ ðx1 r2 Þ .
Eq. (2) governs the transition PDF of the response for the dynamical system (1). The equations for the first and second order moments on the basis of Gaussian closure method are given as 8 m_ 10 > > > > > > > m_ 01 > > > > > > > > > > > > > m_ 20 > > > < m_ 11 > > > > > > > > > > > > > > m_ 02 > > > > > > > > > :
¼ m01 cb ¼ 2 ð3m01 m02 2m301 Þ x21 eð3m10 m20 2m310 Þ cm01 x1 x21 m10 þ x1 r1 sin xt ¼ 2m11 ¼ m02 þ
cb ð2m301 m10 3m02 m11 Þ þ x21 eð2m410 3m220 Þ x21
cm11 x21 m20 þ x1 r1 m10 sin xt cb ¼ 2 ð4m401 6m202 Þ þ x21 eð4m310 m01 6m20 m11 Þ 2cm02 x1 2x21 m11 þ 2x1 r1 m01 sin xt þ x21 r22
j
where mij ¼ E½X it X_ t .
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
1215
3. The solutions of path integration Here we compute the instantaneous PDF in three cases for system (1) using path integration [9]. Next we extend the method proposed in [9] to compute the averaged stationary PDF, similarly using the Gauss–Legendre formula. Firstly, introducing the definition of averaged PDF [11] Z 0 1 t þT P ðx1 ; x2 Þ ¼ pðx1 ; x2 ; tÞ dt; ð3Þ T t0 where T = 2p/x. Then by Gauss–Legendre formula the numerical computation of the averaged PDF is easy to be achieved. Integral (3) can be discretized into the following composite Gauss–Legendre quadrature form P ðx1 ; x2 Þ ¼
Lk K 1 X dk X ckl pðx1 ; x2 ; tkl Þ; T k¼1 2 l¼1
ð30 Þ
where K is the number of subintervals in [t 0 , t 0 + T], Lk is the number of quadrature points in sub-interval k, dk is the length of subinterval k, each tkl is the position of a Gauss quadrature point in time space different from [9], and ckl is the corresponding weight. The method to approximate p(x1, x2, tkl) refers to Yu and Lin [9]. When t 0 in Eq. (3) is large enough, the averaged PDF over every forcing period T is an asymptotically stationary PDF. 3.1. Case 1 b = 0.0 That is X€ t þ cX_ t þ x21 ð1 þ eX 2t ÞX t ¼ x1 r1 sin xt þ x1 r2 nt .
ð4Þ
Eq. (4) is a Duffing oscillator subject to sinusoidal and white noise excitations. The presence of the non-linearity causes multi-valued regions where more than one mean-square value of the response is possible. When c = 0.2, e = 0.3, x1 = 1.0, r1 = 1.0 and r2 = 0.0, we consider the deterministic case of Eq. (4) via method of harmonic balance [12]. Then the variation of the response amplitude vs. the frequency x of harmonic excitation is shown in Fig. 1. In Fig. 1, the curve between two vertical dashed lines shows that the amplitude has three solutions. The smallest and biggest of them are stable (depicted by solid lines). The last one is unstable (depicted by dashed line), i.e. it is infeasible solution. In practice, which stable solutions will prevail depends on the initial condition. Now let x = 1.6, and show the influence of stochastic excitation to Duffing oscillator. The steady state response of system changes from a limit cycle to a
1216
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
Fig. 1. Amplitude vs. frequency x for the Duffing oscillator (4) at c = 0.2, e = 0.3, x1 = 1.0, r1 = 1.0, r2 = 0.0.
diffused limit cycle with variety of the stochastic excitation [13]. From aspect of the PDF of response, it is expected that the steady state PDF of the response will have two peaks. Each peak corresponds to a stable solution of the response. In case of Duffing oscillators under harmonic excitation, there exist two stable periodic solutions for a specific range of excitation frequencies. If there is additionally a white noise excitation, with a certain probability, the response process can jump between the two stable solutions. The phenomenon can be captured by the numerical path integration, not only for the instantaneous PDF but also the averaged stationary PDF. The method of computation for the former has been carried out for Duffing oscillator. However, the numerical method for the latter is a new performance. The results obtained by path integration are displayed in Figs. 2–4. After 200T, PDF of Duffing oscillator (4) has reached steady state and has periodically changing pattern, shown in Fig. 2. For this particular set of parameters of system and excitations, the two dominant peaks are not equilibratory. One of them is great and concentrated, and the other is small and dispersed. They are the results of the disturbance by white noise and harmonic excitations. This analysis is in good agreement with the results in Yu and Lin [9]. The double crater-like shape of the averaged stationary PDF is successfully computed by Eq. (3 0 ) shown in Fig. 3. It is easy to find that there indeed exist two stable solutions of the response in a driven period. The averaging jump
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
1217
Fig. 2. Surface plots of PDF of the Duffing oscillator (4) at t = 1/4, 2/4, 3/4, 4/4 of a period when r2 = 0.4.
Fig. 3. The averaged stationary PDF of the Duffing oscillator (4) at 1T.
between the two peak happens every a period. With time pasted, the result becomes as Fig. 4. Next let r2 = 0.1. Fig. 5 clearly indicates when the intensity of stochastic excitation decreases, two peaks of the PDF separate from each other. But the periodicity of PDF still persists in Fig. 6. The harmonic excitation is the main deterministic factor to the shape of PDF. For the smaller white noise
1218
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
Fig. 4. The averaged stationary PDF of the Duffing oscillator (4) at 7T.
Fig. 5. Surface plots of PDF of the Duffing oscillator (4) at t = 5T when r2 = 0.1.
excitation, the probability of the transition between two peaks largely decreases. In Figs. 2–6, we select time step Dt = T/4, and set 100 · 100 Gaussian quadrature points in 50 · 50 subinterval within domain [5, 5] · [5, 5]. The initial distribution is given by ( ) 2 2 ð0Þ ð0Þ 1 ðx l Þ ð_ x l Þ 1 2 pðxð0Þ ; x_ ð0Þ Þ ¼ pffiffiffiffiffiffiffiffi exp ; ð5Þ 2p s1 s2 2s1 2s2 where l1 = 2.0, l2 = 1.8, s1 = s2 = 1.0.
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
1219
Fig. 6. Surface plots of PDF of the Duffing oscillator (4) at t = 1/4, 2/4, 3/4, 4/4 of a period when r2 = 0.1.
3.2. Case 2 e = 0.0 The system (1) can be written as b _2 _ € X t þ c 1 þ 2 X t X t þ x21 X t ¼ x1 r1 sin xt þ x1 r2 nt . x1
ð6Þ
In fact, system (6) is Rayleigh oscillator subject to harmonic and stochastic excitations. When r2 = 0.0, consider the deterministic case of system (6) via the method of multiple scales [12]. Let q = a2/4, a is amplitude, r = (xx1)/c, and r1 = ck, 0 < c 1. When k 2 ¼ 23, from Fig. 7 about amplitude q vs. frequency r, there are three values of amplitude for the special range of r. The solid line depicts the solution of response. The dash line separates the solution into two parts. However, different from Case 1, one part lying up the dash line is stable, and another is unstable. That is to say, just the largest vibration can realize in physics, and the remaining two are unpractical. For the steady state PDF of the response in a period, it is not expected to split, but periodic, just as shown in Fig. 8.
1220
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
Fig. 7.pffiffi The amplitude q vs. frequency r for the Rayleigh oscillator (6). c = 0.1, b ¼ 13, x1 = 1.0, r1 ¼ 306 and r2 = 0.0.
The averaged stationary PDF also indicates that there is only one state of response with large probability in average meaning. It is ready to verify by path integration, shown in Fig. 9. Thepffiffi applied parameters where are c = 0.1, b ¼ 13, x1 = 1.0, x = 1.016, r1 ¼ 306, r2 = 0.4. To apply path integration, then the computational parameters in Figs. 8 and 9 are set as 80 · 80 Gaussian quadrature points in 40 · 40 subinterval within domain [4, 4] · [4, 4]. The initial distribution is given by Eq. (5), here l1 = 0.0, l2 = 0.0, s1 = s2 = 0.1, and time step Dt = T/4. 3.3. Case 3 b 5 0.0, e 5 0.0 In Fig. 10, there is no range of multi-value of the amplitude when the frequency varies. For Case 3, whatever the initial value is given, after long time the phase portrait is always the same. So the instantaneous PDF is of single peak and periodic by numerical path integration, described in Fig. 11. Likewise, the crater-like averaged stationary PDF in a driven forcing period is easy to be calculated using the method of path integration proposed in this section, shown in Fig. 12. The result agrees well with that of Von Wagner and Wedig [4]. The shape of the averaged stationary PDF is a result of the disturbance by strong harmonic excitation and weak white noise.
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
1221
Fig. 8. Contour plots of PDF of the Rayleigh oscillator (6) at t = 1/4, 2/4, 3/4, 4/4 of a period.
Fig. 9. The averaged stationary PDF of the Rayleigh oscillator (6).
The applied parameters of Figs. 11 and 12 are c = 0.5, b = 0.3, e = 0.3, x1 = 1.0, x = 2.0, r1 = 1.0, r2 = 0.1 and 60 · 60 Gaussian quadrature points
1222
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
Fig. 10. The amplitude vs. frequency x for the Duffing–Rayleigh oscillator (1) at c = 0.5, b = 0.3, e = 0.3, x1 = 1.0, r1 = 1.0 and r2 = 0.0.
Fig. 11. Surface plots of PDF of the Duffing–Rayleigh oscillator (1) at t = 1/4, 2/4, 3/4, 4/4 of a period.
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
1223
Fig. 12. The averaged stationary PDF of the Duffing–Rayleigh oscillator (1).
in 30 · 30 subinterval within domain [1, 1] · [1, 1]. The initial distribution is given by Eq. (5), here l1 = 0.2, l2 = 0.1, s1 = s2 = 0.1 and time step is Dt = T/4.
4. Conclusions The application of path integration based on Gauss–Legendre formula to the system subject to harmonic and stochastic excitations ratifies the validity of this method to capture behavior of non-linear system. Up to now the analytical solution of FPK equation is a heavy challenging problem, and is only feasible for some special classes of non-linear systems. Recently, researchers are devoted to more powerful numerical methods. The results of this paper show that the method of path integration can capture the steady state periodic solution of PDF for the Duffing–Rayleigh oscillator. Path integration method depicts the facts of multi-value and periodic jump phenomenon in a forcing period indeed. On the other hand, we propose a numerical path integral method for calculating the averaged stationary PDF similarly via Gauss–Legendre formula. This method takes advantage of the periodic coefficient of FPK equation of dynamical system, and is efficient for calculation indeed. The double crater-like PDF of Duffing system (4) and the crater-like PDF of Duffing–Rayleigh system (1) are successfully calculated by this method. At the same time, from the shape of the averaged PDF, we can find and illustrate the facts of multi-value and jump phenomenon for the Duffing–Rayleigh oscillator.
1224
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
Finally, the noise-induced state change of PDF is also shown. Harder the PDF splits, weaker the intensity of stochastic excitation is in some special parameters of the Duffing–Rayleigh oscillator.
Acknowledgements The work is supported by the National Natural Science Foundation of China (Grant No. 10472091 and 10332030) and Natural Science Foundation of Shaanxi Province (Grant No. 2003A03). All calculations were performed at the Center of High Performance Computing of Northwestern Polytechnical University.
References [1] T.K. Caughey, F. Ma, The exact steady-state solution of a class of nonlinear stochastic systems, International Journal of Non-Linear Mechanics 17 (1982) 137–142. [2] M.F. Dimentberg, An exact solution to a certain non-linear random vibration problem, International Journal of Non-Linear Mechanics 17 (1982) 231–236. [3] G.Q. Cai, Y.K. Lin, Exact and approximate solutions for randomly excited MDOF nonlinear systems, International Journal of Non-Linear Mechanics 31 (1996) 647–655. [4] U. Von Wagner, W.V. Wedig, On The calculation of stationary solutions of multi-dimensional Fokker–Planck equations by orthogonal functions, Nonlinear Dynamics 21 (2000) 289–306. [5] M.Di. Paola, A. Sofi, Approximate solution of the Fokker–Planck–Kolmogorov equation, Probabilistic Engineering Mechanics 17 (2002) 369–384. [6] A. Naess, V. Moe, Efficient path integration methods for nonlinear dynamic systems, Probabilistic Engineering Mechanic 15 (2000) 221–231. [7] Z.L. Huang, W.Q. Zhu, Y. Suzuki, Stochastic averaging of strongly non-linear oscillators under combined harmonic and white noise excitations, Journal of Sound and Vibration 238 (2000) 233–256. [8] J.S. Yu, G.Q. Cai, Y.K. Lin, A new path integration procedure based on Gauss–Legendre scheme, International Journal of Non-Linear Mechanics 32 (1997) 759–768. [9] J.S. Yu, Y.K. Lin, Numerical path integration of a nonlinear oscillator subject to both sinusoidal and white noise excitations, International Journal of Non-Linear Mechanics 39 (2004) 1493–1500. [10] A.H. Nayfeh, S.J. Serhan, Response statistics of nonlinear systems to combined deterministic and random excitations, International Journal of Non-Linear Mechanics 25 (1990) 493–515. [11] A. Naess, Chaos and nonlinear stochastic dynamics, Probabilistic Engineering Mechanics 15 (2000) 37–47. [12] J.Q. Zhou, Y.Y. Zhu, Nonlinear Vibrations, XiÕan Jiaotong University Press, 1998 (in Chinese). [13] H.W. Rong, W. Xu, G. Meng, T. Fang, Response of nonlinear oscillator to combined harmonic and random excitation, Chinese Journal of Applied Mechanics 18 (2001) 32–36 (in Chinese).
www.elsevier.com/locate/amc
Study of the Duffing–Rayleigh oscillator subject to harmonic and stochastic excitations by path integration W.X. Xie a, W. Xu
a,*
, L. Cai
b
a
b
Department of Applied Mathematics, Northwestern Polytechnical University, Xi’an, Shaanxi 710072, PR China College of Astronautics, Northwestern Polytechnical University, Xi’an, Shaanxi 710072, PR China
Abstract In this paper, the Duffing–Rayleigh oscillator subject to harmonic and stochastic excitations is investigated via path integration based on the Gauss–Legendre integration formula. The method can successfully capture the steady state periodic solution of probability density function. This path integration method, using the periodicity of the coefficient of associated Fokker–Planck–Kolmogorov equation, is extended to deal with the averaged stationary probability density, and is efficient to computation. Meanwhile, the changes of probability density caused by the intensities of harmonic and stochastic excitations, are discussed in three cases through the instantaneous probability density and the averaged stationary probability density. 2005 Elsevier Inc. All rights reserved. Keywords: Path integration; Duffing–Rayleigh oscillator; Harmonic and stochastic excitations; Probability density
*
Corresponding author. E-mail address: [email protected] (W. Xu).
0096-3003/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.03.018
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
1213
1. Introduction So far, it is still hard to solve analytically the Fokker–Planck–Kolmogorov (FPK) equations for non-linear stochastic dynamical systems, which govern the transition probability density function (PDF) of the responses. Exact solutions are known just for some special systems (see [1–3] etc.). For this reason, several approximate and numerical techniques have been developed (see e.g. [4–7]). For the dynamical system subject to harmonic and stochastic excitations, one of the drift coefficients of the corresponding FPK equation is periodic. So there is no stationary solution and may exist steady state solution when t ! 1 [7]. Some approximate solutions for such systems using stochastic averaging method are given in Von Wagner and Wedig [4] and Huang et al. [7]. In fact, the method presented by Von Wagner and Wedig [4] can only pursue the approximate averaged stationary PDF of Duffing–Rayleigh oscillator subject to harmonic and white noise excitations by orthogonal function expansions, e.g. especially adjusted polynomials and Fourier series. The weighted functions of the new polynomials are obtained by the application of the stochastic averaging method. However, the calculation of double crater-like density for Duffing system is hard to perform. Huang et al. [7] get the averaging FPK equation by stochastic averaging method, and resort to the numerical path integration at last. Generally speaking, the approximate methods above have no way to capture the periodic changes of the instantaneous PDF for the response before averaging. But the numerical path integration method based on Gauss– Legendre integration formula proposed by Yu et al. [8,9], is a powerful tool to capture the evolution of the instantaneous PDF. Especially, Yu finds the evolution of the PDF to the periodic steady state PDF of Duffing system [9]. In this paper, a complicated system, the Duffing–Rayleigh oscillator subject to harmonic and stochastic excitations is of interest. The system has abundant dynamical behaviors, such as multi-value of response and hopping between stable responses. These phenomena were observed first in the Duffing–Rayleigh oscillator with multiple scales and second-order Gaussian closure method [10]. Taking advantages of this method developed by Yu and Lin [9], this paper not only obtains the periodic steady state PDF and the noise-induced state changes of PDF, but also can calculate the averaged stationary PDF by extending the path integration based on Gauss–Legendre integration formula, which is helpful to interpret the multi-value of response.
2. Duffing–Rayleigh oscillator subject to harmonic and stochastic excitations Consider the non-linear oscillator subject to harmonic and white noise excitations of the form
1214
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
b 2 X€ t þ c 1 þ 2 X_ t X_ t þ x21 ð1 þ eX 2t ÞX t ¼ x1 r1 sin xt þ x1 r2 nt ; x1
ð1Þ
where x1 is the natural circle frequency, c is the linear damping coefficient, b and e are the non-linearity parameters, r1 and x are intensity and frequency of the harmonic excitation respectively, and r2 is the intensity of the additive white noise nt. The corresponding FPK equation can be derived as oQ o o 1 o2 ¼ ½m1 ðx1 ; x2 ; tÞQ ½m2 ðx1 ; x2 ; tÞQ þ ½bðx1 ; x2 ; tÞQ; ot ox1 ox2 2 ox22
ð2Þ
where x1 ¼ x; x2 ¼ x_ ; Q ¼ qðx1 ; x2 ; tx0 ; x0 ; t0 Þ; 1
2
m1 ðx1 ; x2 ; tÞ ¼ x2 ;
m2 ðx1 ; x2 ; tÞ ¼ c 1 þ xb2 x22 x2 x21 ð1 þ ex21 Þx1 þ x1 r1 sin xt; 1
and
2
bðx1 ; x2 ; tÞ ¼ ðx1 r2 Þ .
Eq. (2) governs the transition PDF of the response for the dynamical system (1). The equations for the first and second order moments on the basis of Gaussian closure method are given as 8 m_ 10 > > > > > > > m_ 01 > > > > > > > > > > > > > m_ 20 > > > < m_ 11 > > > > > > > > > > > > > > m_ 02 > > > > > > > > > :
¼ m01 cb ¼ 2 ð3m01 m02 2m301 Þ x21 eð3m10 m20 2m310 Þ cm01 x1 x21 m10 þ x1 r1 sin xt ¼ 2m11 ¼ m02 þ
cb ð2m301 m10 3m02 m11 Þ þ x21 eð2m410 3m220 Þ x21
cm11 x21 m20 þ x1 r1 m10 sin xt cb ¼ 2 ð4m401 6m202 Þ þ x21 eð4m310 m01 6m20 m11 Þ 2cm02 x1 2x21 m11 þ 2x1 r1 m01 sin xt þ x21 r22
j
where mij ¼ E½X it X_ t .
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
1215
3. The solutions of path integration Here we compute the instantaneous PDF in three cases for system (1) using path integration [9]. Next we extend the method proposed in [9] to compute the averaged stationary PDF, similarly using the Gauss–Legendre formula. Firstly, introducing the definition of averaged PDF [11] Z 0 1 t þT P ðx1 ; x2 Þ ¼ pðx1 ; x2 ; tÞ dt; ð3Þ T t0 where T = 2p/x. Then by Gauss–Legendre formula the numerical computation of the averaged PDF is easy to be achieved. Integral (3) can be discretized into the following composite Gauss–Legendre quadrature form P ðx1 ; x2 Þ ¼
Lk K 1 X dk X ckl pðx1 ; x2 ; tkl Þ; T k¼1 2 l¼1
ð30 Þ
where K is the number of subintervals in [t 0 , t 0 + T], Lk is the number of quadrature points in sub-interval k, dk is the length of subinterval k, each tkl is the position of a Gauss quadrature point in time space different from [9], and ckl is the corresponding weight. The method to approximate p(x1, x2, tkl) refers to Yu and Lin [9]. When t 0 in Eq. (3) is large enough, the averaged PDF over every forcing period T is an asymptotically stationary PDF. 3.1. Case 1 b = 0.0 That is X€ t þ cX_ t þ x21 ð1 þ eX 2t ÞX t ¼ x1 r1 sin xt þ x1 r2 nt .
ð4Þ
Eq. (4) is a Duffing oscillator subject to sinusoidal and white noise excitations. The presence of the non-linearity causes multi-valued regions where more than one mean-square value of the response is possible. When c = 0.2, e = 0.3, x1 = 1.0, r1 = 1.0 and r2 = 0.0, we consider the deterministic case of Eq. (4) via method of harmonic balance [12]. Then the variation of the response amplitude vs. the frequency x of harmonic excitation is shown in Fig. 1. In Fig. 1, the curve between two vertical dashed lines shows that the amplitude has three solutions. The smallest and biggest of them are stable (depicted by solid lines). The last one is unstable (depicted by dashed line), i.e. it is infeasible solution. In practice, which stable solutions will prevail depends on the initial condition. Now let x = 1.6, and show the influence of stochastic excitation to Duffing oscillator. The steady state response of system changes from a limit cycle to a
1216
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
Fig. 1. Amplitude vs. frequency x for the Duffing oscillator (4) at c = 0.2, e = 0.3, x1 = 1.0, r1 = 1.0, r2 = 0.0.
diffused limit cycle with variety of the stochastic excitation [13]. From aspect of the PDF of response, it is expected that the steady state PDF of the response will have two peaks. Each peak corresponds to a stable solution of the response. In case of Duffing oscillators under harmonic excitation, there exist two stable periodic solutions for a specific range of excitation frequencies. If there is additionally a white noise excitation, with a certain probability, the response process can jump between the two stable solutions. The phenomenon can be captured by the numerical path integration, not only for the instantaneous PDF but also the averaged stationary PDF. The method of computation for the former has been carried out for Duffing oscillator. However, the numerical method for the latter is a new performance. The results obtained by path integration are displayed in Figs. 2–4. After 200T, PDF of Duffing oscillator (4) has reached steady state and has periodically changing pattern, shown in Fig. 2. For this particular set of parameters of system and excitations, the two dominant peaks are not equilibratory. One of them is great and concentrated, and the other is small and dispersed. They are the results of the disturbance by white noise and harmonic excitations. This analysis is in good agreement with the results in Yu and Lin [9]. The double crater-like shape of the averaged stationary PDF is successfully computed by Eq. (3 0 ) shown in Fig. 3. It is easy to find that there indeed exist two stable solutions of the response in a driven period. The averaging jump
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
1217
Fig. 2. Surface plots of PDF of the Duffing oscillator (4) at t = 1/4, 2/4, 3/4, 4/4 of a period when r2 = 0.4.
Fig. 3. The averaged stationary PDF of the Duffing oscillator (4) at 1T.
between the two peak happens every a period. With time pasted, the result becomes as Fig. 4. Next let r2 = 0.1. Fig. 5 clearly indicates when the intensity of stochastic excitation decreases, two peaks of the PDF separate from each other. But the periodicity of PDF still persists in Fig. 6. The harmonic excitation is the main deterministic factor to the shape of PDF. For the smaller white noise
1218
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
Fig. 4. The averaged stationary PDF of the Duffing oscillator (4) at 7T.
Fig. 5. Surface plots of PDF of the Duffing oscillator (4) at t = 5T when r2 = 0.1.
excitation, the probability of the transition between two peaks largely decreases. In Figs. 2–6, we select time step Dt = T/4, and set 100 · 100 Gaussian quadrature points in 50 · 50 subinterval within domain [5, 5] · [5, 5]. The initial distribution is given by ( ) 2 2 ð0Þ ð0Þ 1 ðx l Þ ð_ x l Þ 1 2 pðxð0Þ ; x_ ð0Þ Þ ¼ pffiffiffiffiffiffiffiffi exp ; ð5Þ 2p s1 s2 2s1 2s2 where l1 = 2.0, l2 = 1.8, s1 = s2 = 1.0.
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
1219
Fig. 6. Surface plots of PDF of the Duffing oscillator (4) at t = 1/4, 2/4, 3/4, 4/4 of a period when r2 = 0.1.
3.2. Case 2 e = 0.0 The system (1) can be written as b _2 _ € X t þ c 1 þ 2 X t X t þ x21 X t ¼ x1 r1 sin xt þ x1 r2 nt . x1
ð6Þ
In fact, system (6) is Rayleigh oscillator subject to harmonic and stochastic excitations. When r2 = 0.0, consider the deterministic case of system (6) via the method of multiple scales [12]. Let q = a2/4, a is amplitude, r = (xx1)/c, and r1 = ck, 0 < c 1. When k 2 ¼ 23, from Fig. 7 about amplitude q vs. frequency r, there are three values of amplitude for the special range of r. The solid line depicts the solution of response. The dash line separates the solution into two parts. However, different from Case 1, one part lying up the dash line is stable, and another is unstable. That is to say, just the largest vibration can realize in physics, and the remaining two are unpractical. For the steady state PDF of the response in a period, it is not expected to split, but periodic, just as shown in Fig. 8.
1220
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
Fig. 7.pffiffi The amplitude q vs. frequency r for the Rayleigh oscillator (6). c = 0.1, b ¼ 13, x1 = 1.0, r1 ¼ 306 and r2 = 0.0.
The averaged stationary PDF also indicates that there is only one state of response with large probability in average meaning. It is ready to verify by path integration, shown in Fig. 9. Thepffiffi applied parameters where are c = 0.1, b ¼ 13, x1 = 1.0, x = 1.016, r1 ¼ 306, r2 = 0.4. To apply path integration, then the computational parameters in Figs. 8 and 9 are set as 80 · 80 Gaussian quadrature points in 40 · 40 subinterval within domain [4, 4] · [4, 4]. The initial distribution is given by Eq. (5), here l1 = 0.0, l2 = 0.0, s1 = s2 = 0.1, and time step Dt = T/4. 3.3. Case 3 b 5 0.0, e 5 0.0 In Fig. 10, there is no range of multi-value of the amplitude when the frequency varies. For Case 3, whatever the initial value is given, after long time the phase portrait is always the same. So the instantaneous PDF is of single peak and periodic by numerical path integration, described in Fig. 11. Likewise, the crater-like averaged stationary PDF in a driven forcing period is easy to be calculated using the method of path integration proposed in this section, shown in Fig. 12. The result agrees well with that of Von Wagner and Wedig [4]. The shape of the averaged stationary PDF is a result of the disturbance by strong harmonic excitation and weak white noise.
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
1221
Fig. 8. Contour plots of PDF of the Rayleigh oscillator (6) at t = 1/4, 2/4, 3/4, 4/4 of a period.
Fig. 9. The averaged stationary PDF of the Rayleigh oscillator (6).
The applied parameters of Figs. 11 and 12 are c = 0.5, b = 0.3, e = 0.3, x1 = 1.0, x = 2.0, r1 = 1.0, r2 = 0.1 and 60 · 60 Gaussian quadrature points
1222
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
Fig. 10. The amplitude vs. frequency x for the Duffing–Rayleigh oscillator (1) at c = 0.5, b = 0.3, e = 0.3, x1 = 1.0, r1 = 1.0 and r2 = 0.0.
Fig. 11. Surface plots of PDF of the Duffing–Rayleigh oscillator (1) at t = 1/4, 2/4, 3/4, 4/4 of a period.
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
1223
Fig. 12. The averaged stationary PDF of the Duffing–Rayleigh oscillator (1).
in 30 · 30 subinterval within domain [1, 1] · [1, 1]. The initial distribution is given by Eq. (5), here l1 = 0.2, l2 = 0.1, s1 = s2 = 0.1 and time step is Dt = T/4.
4. Conclusions The application of path integration based on Gauss–Legendre formula to the system subject to harmonic and stochastic excitations ratifies the validity of this method to capture behavior of non-linear system. Up to now the analytical solution of FPK equation is a heavy challenging problem, and is only feasible for some special classes of non-linear systems. Recently, researchers are devoted to more powerful numerical methods. The results of this paper show that the method of path integration can capture the steady state periodic solution of PDF for the Duffing–Rayleigh oscillator. Path integration method depicts the facts of multi-value and periodic jump phenomenon in a forcing period indeed. On the other hand, we propose a numerical path integral method for calculating the averaged stationary PDF similarly via Gauss–Legendre formula. This method takes advantage of the periodic coefficient of FPK equation of dynamical system, and is efficient for calculation indeed. The double crater-like PDF of Duffing system (4) and the crater-like PDF of Duffing–Rayleigh system (1) are successfully calculated by this method. At the same time, from the shape of the averaged PDF, we can find and illustrate the facts of multi-value and jump phenomenon for the Duffing–Rayleigh oscillator.
1224
W.X. Xie et al. / Appl. Math. Comput. 172 (2006) 1212–1224
Finally, the noise-induced state change of PDF is also shown. Harder the PDF splits, weaker the intensity of stochastic excitation is in some special parameters of the Duffing–Rayleigh oscillator.
Acknowledgements The work is supported by the National Natural Science Foundation of China (Grant No. 10472091 and 10332030) and Natural Science Foundation of Shaanxi Province (Grant No. 2003A03). All calculations were performed at the Center of High Performance Computing of Northwestern Polytechnical University.
References [1] T.K. Caughey, F. Ma, The exact steady-state solution of a class of nonlinear stochastic systems, International Journal of Non-Linear Mechanics 17 (1982) 137–142. [2] M.F. Dimentberg, An exact solution to a certain non-linear random vibration problem, International Journal of Non-Linear Mechanics 17 (1982) 231–236. [3] G.Q. Cai, Y.K. Lin, Exact and approximate solutions for randomly excited MDOF nonlinear systems, International Journal of Non-Linear Mechanics 31 (1996) 647–655. [4] U. Von Wagner, W.V. Wedig, On The calculation of stationary solutions of multi-dimensional Fokker–Planck equations by orthogonal functions, Nonlinear Dynamics 21 (2000) 289–306. [5] M.Di. Paola, A. Sofi, Approximate solution of the Fokker–Planck–Kolmogorov equation, Probabilistic Engineering Mechanics 17 (2002) 369–384. [6] A. Naess, V. Moe, Efficient path integration methods for nonlinear dynamic systems, Probabilistic Engineering Mechanic 15 (2000) 221–231. [7] Z.L. Huang, W.Q. Zhu, Y. Suzuki, Stochastic averaging of strongly non-linear oscillators under combined harmonic and white noise excitations, Journal of Sound and Vibration 238 (2000) 233–256. [8] J.S. Yu, G.Q. Cai, Y.K. Lin, A new path integration procedure based on Gauss–Legendre scheme, International Journal of Non-Linear Mechanics 32 (1997) 759–768. [9] J.S. Yu, Y.K. Lin, Numerical path integration of a nonlinear oscillator subject to both sinusoidal and white noise excitations, International Journal of Non-Linear Mechanics 39 (2004) 1493–1500. [10] A.H. Nayfeh, S.J. Serhan, Response statistics of nonlinear systems to combined deterministic and random excitations, International Journal of Non-Linear Mechanics 25 (1990) 493–515. [11] A. Naess, Chaos and nonlinear stochastic dynamics, Probabilistic Engineering Mechanics 15 (2000) 37–47. [12] J.Q. Zhou, Y.Y. Zhu, Nonlinear Vibrations, XiÕan Jiaotong University Press, 1998 (in Chinese). [13] H.W. Rong, W. Xu, G. Meng, T. Fang, Response of nonlinear oscillator to combined harmonic and random excitation, Chinese Journal of Applied Mechanics 18 (2001) 32–36 (in Chinese).