- Email: [email protected]
Noise-induced chaos in a quadratically nonlinear oscillator
Chaos, Solitons and Fractals 30 (2006) 920–929 www.elsevier.com/locate/chaos
Noise-induced chaos in a quadratically nonlinear oscillator Chunbiao Gan Department of Mechanics, CMEE, Zhejiang University, Hangzhou 310027, China Accepted 30 August 2005
Abstract The present paper focuses on the noise-induced chaos in a quadratically nonlinear oscillator. Simple zero points of the stochastic Melnikov integral theoretically mean the necessary rising of noise-induced chaotic response in the system based on the stochastic Melnikov method. To quantify the noise-induced chaos, the boundary of the systemÕs safe basin is firstly studied and it is shown to be incursively fractal when chaos arises. Three cases are considered in simulating the safe basin of the system, i.e., the system is excited only by the harmonic excitation, by both the harmonic and the Gaussian white noise excitations, and only by the Gaussian white noise excitation. Secondly, the leading Lyapunov exponent by RosensteinÕs algorithm is shown to quantify the chaotic nature of the sample time series of the system. The results show that the boundary of the safe basin can also be fractal even if the system is excited only by the external Gaussian white noise. Most importantly, the almost-harmonic, the noise-induced chaotic and the thoroughly random responses can be found in the system. Ó 2005 Elsevier Ltd. All rights reserved.
1. Introduction Since 1990s, stochastic Melnikov method has been applied to study the effect of noise on homoclinic or heteroclinic bifurcation and noise-induced chaos [1–6]. It is a consequence of (1) WigginsÕ [7] result that systems capable of chaotic behavior under harmonic excitation can also behave chaotically if the excitation is quasi-periodic, and (2) the application of the stochastic Melnikov approach to systems excited by physically realizable stochastic processes, which can be approximated arbitrarily closely by random quasi-periodic sums [3]. For both deterministic and stochastic systems the Melnikov conditions are necessary but not sufficient. Results yielded by the Melnikov approach are therefore weak, though less so far transient than for steady state motions. Generally, it is not easy to quantify noise-induced chaos. One knows that an important goal of studying dynamical systems is to determine their global structures, and one of these global structures is the boundary of basin. Study on safe basin in oscillatory system with potential well can be found in [8–14]. The non-stationary effect on safe basin and the analysis of noise-induced chaos in a forced softening Duffing oscillator are presented in [15,16]. Because of the coexistence of period and chaotic attractors, the basin boundaries of attractors are usually fractal and naturally incursive, in the sense that they are related to homoclinic or heteroclinic intersections of stable and unstable manifolds of the saddle points in the system, which usually means that chaos arises in such system. E-mail addresses: [email protected], [email protected] 0960-0779/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2005.08.157
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929
921
The most striking feature of chaos is the unpredictability of the future, which is called sensitive dependence on initial conditions. It can be quantified by the Lyapunov exponent, which characterizes the average exponential rate of divergence or convergence of nearby orbits in phase space of the system and quantifies the strength of chaos. Till now, most authors employed the algorithm presented by Wolf et al. [17] to compute the Lyapunov exponent even for a stochastic dynamic system. However, as pointed out by Kantz and Schreiber [18], the algorithm by Wolf is not very robust and one can easily obtain the wrong results. WolfÕs algorithm does not allow one to test for the presence of exponential divergence, but just assumes its existence and thus yields a finite exponent for stochastic data also, where the true exponent is infinite. While WolfÕs algorithm only uses a delay reconstruction of phase space, there is another class of algorithms which also involves the approximation of the underlying deterministic dynamics. References are Sano and Sawada [19], Eckmann et al. [20] and Rosenstein et al. [21]. They show that, though the output data are generally affected by noise, the influence from the noise can be minimized by using the averaging statistics. This paper studies the effect of the Gaussian white noise on the boundary of safe basin and the noise-induced chaotic response in a quadratically nonlinear oscillator. In Section 2, the stochastic Melnikov integral is analytically presented based on the dynamical theory and the nature of a linear time-invariant filter. Section 3 presents the variation of safe basin numerically when one adjusts the amplitude of the harmonic excitation in the deterministic case, or the strength of the Gaussian white noise excitation in the stochastic case of the system, from which one can see that fractal boundary appears even if the oscillator is excited only by the Gaussian white noise excitation. In Section 4, noise-induced chaotic responses are simulated and testified by employing the algorithm for the leading Lyapunov exponent from Rosenstein et al. [21]. Section 5 gives the summary and discussion.
2. Stochastic Melnikov integral and its single zero points The quadratically nonlinear oscillator considered in the present paper can be written as follows: x_ ¼ y y_ ¼ x þ x2 þ eðly þ f1 cos xt þ f2 nðtÞÞ;
ð1Þ
where 0 < e 1, l and f1 are the damping coefficient and the amplitude of the harmonic excitation respectively, and n(t) is taken as the Gaussian white noise, f2 is the strength of this stochastic excitation. The spectrum density for the Gaussian white noise is assumed to be S0 with the strength D(=2pS0). For system (1) the preliminary assumptions for the Melnikov approach are those listed by Wiggins [7] as well as the constraint for the random excitation: n(t) is uniformly continuous and bounded in ensemble, i.e., for any d1 > 0, there exists d2 > 0, such that jn(t1) n(t2)j < d1 when jt1 t2j < d2 for any times t1 and t2 and any idealization of n(t). One can easily check that system (1) satisfies the above assumptions, and the hyperbolic fixed points of system (1) will not disappear under the harmonic and the random excitations, so the Melnikov approach can also be employed to analyze the stochastic system (1). From [3,6,7], the Melnikov integral for system (1) can be expressed as Mðt0 ; l; f1 ; f2 Þ ¼ M d ðt0 ; l; f1 Þ þ M r ðt0 ; f2 Þ;
ð2Þ
where Md(t0; l, f1) is the deterministic part of M(t0; l, f1, f2), i.e., the mean value of the Melnikov integral, and Mr(t0; f2) is the corresponding random part of M(t0; l, f1, f2) due to the Gaussian white noise excitation. By the direct integration and the residue theory, one can obtain that 6 12pf1 x2 epx cos xt0 . M d ðt0 ; l; f1 Þ ¼ l þ 5 ð1 e2px Þ
ð3Þ
The random portion Mr(t0; f2) due to the Gaussian white noise, instead of directly integrated, can be estimated by considering the convolution integral as a filtering process, i.e., a stationary process n(t) passing through a linear time-invariant filter given by the homoclinic orbit [3,6]. Because of the linear time-invariant filtering, the random process Mr(t0; f2) is still stationary and of zero mean (E[Mr(t0; f2)] = 0). Its corresponding spectra may be calculated by directly multiplying the original noise spectrum S nðtþt0 Þ and the transfer function of the linear time-invariant filter. For the Gaussian white noise excitation, one has Z þ1 jF ðXÞj2 S nðtþt0 Þ ðXÞdX; ð4Þ r2M r ¼ 1
R þ1 where F ðXÞ ¼ f2 1 y 0 ðtÞeiXt dt ¼ 12pf2 X2 epX i=ð1 e2pX Þ and S nðtþt0 Þ ðXÞ S 0 is the spectrum density of the Gaussian white noise process. Integrating Eq. (4) numerically yields
922
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929
r2M r ¼ 0:764p2 S 0 f22 .
ð5Þ
When f2 5 0 in system (1), the Melnikov integral measures the random distance between the stable and the unstable manifolds. The threshold value for the rising of the chaotic response will depend on the property of the random excitation process, and may deviate from the mean value. From energy point of view [6], the necessary condition is given by 4p2 f12 x4 e2px =ð1 e2px Þ2 þ 0:764p2 S 0 f22 P
1 2 l. 25
ð6Þ
Inequality (6) tells one that the chaotic response may be induced by the Gaussian white noise excitation. Moreover, it deserves to note that inequality (6) is not sufficient, and other measures are needed for the identification of the systemÕs chaotic response. In the following two sections, the incursively fractal boundary of safe basin of the system and the positive leading exponent are presented respectively. In the present paper, we all set e = 0.1, x = 1.0 and the stepsize h = 0.01 if not pointed out specifically.
3. Incursively fractal boundary of safe basin The decrease of safe basinÕs area is often called basin erosion. Since those intersections of stable and unstable manifolds map one to another, the manifolds present a convoluted structure that extends through a wide region of the phase space. The transient orbits near such a fractal basin boundary will be as convoluted as the boundary itself. So they are expected to cross the critical line very often, increasing dramatically the number of unsafe initial conditions. In view of this point, the erosion of the safe basin is usually related to fractal basin boundary of attractors. The safe basin without erosion is drawn within the region G surrounded by the coincident stable and unstable manifolds of the conservative case of system (1) with l = f1 = f2 = 0 (see Fig. 1), which is defined by sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ( ffi ffi) 1 2 1 3 1 2 1 3 ð7Þ G ¼ ðx; yÞ 2 Rj 0:5 6 x 6 1; 2 H x þ x 6 y 6 2 H x þ x 2 3 2 3 for generating 400 grid lines in x-direction and an increasing series of grid lines from 1 to 200 when 0.5 6 x 6 0.25 and a decreasing series from 199 to 1 when 0.25 < x 6 1.0 in y-direction. Each grid point is set as an initial condition to perform the simulation of system (1). When the Hamiltonian value of any phase point (x(t), y(t)) of the systemÕs trajectory is larger than H within 100 thousand steps, this motion initiating from the corresponding initial point is considered to be unsafe and this initial point is discarded.
Fig. 1. Safe basin without erosion in the deterministic case of system (1), where l = f1 = f2 = 0.0.
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929
923
For the stochastic case of the system, 16 sample responses are simulated with the same initial condition (picked up from the grid points) and 16 different realizations of the random processes due to Monte-Carlo method [22]. Only the initial point being safe for all the 16 sample responses is recorded. 3.1. Erosion of safe basin and effect of the Gaussian white noise In Fig. 1, the safe basin is a closed and apparently smooth curve when f1 = 0. Following the increase of the driving amplitude f1, this orbit will undergo a bifurcation cascade to chaos, and the boundary of the safe basin shows to be fractal (see Fig. 2, in which H is set to be 1/6). In general, the area of the safe basin decreases following the decrease of the Hamiltonian H. To learn the effect of the noise excitation on the erosion of safe basin, Fig. 3 presents a sequence of safe basins when the Gaussian white noise is added to the system, i.e., f2 5 0. The incursively fractal fingers are also observed, which means that chaotic responses still exist in the stochastic case of system (1). From Fig. 3, one knows that the Gaussian
Fig. 2. The boundary of safe basin in the deterministic case of system (1), where l = 1.0, f2 = 0.
Fig. 3. The effect of the Gaussian white noise on the boundary of safe basin in system (1), where l = 1.0, f1 = 1.0, D = 1.0.
924
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929
white noise excitation can aggravate the erosion of the systemÕs safe basin. In addition, fractal boundary can also appear when the Hamiltonian H is taken as other positive value, which is not given here. From Fig. 4, one knows that the Gaussian white noise excitation can play a dispersive role to the systemÕs dynamics, the response of the system can still remain in the safe basin within 100 thousand steps even if the system without damping item is driven by a periodical force. In this case, incursively fractal boundary of the safe basin also means that noise excitation can induce chaos, which can be further shown in the following studies. 3.2. Noise-induced fractal boundary of safe basin From Figs. 3 and 4, the response initiating from the potential well may not escape from the well within 100 thousand steps when the system is excited by both the harmonic and the Gaussian white noise excitation, and the boundary of the safe basin can show to be incursively fractal. According to inequality (6) given in Section 2, the Melnikov condition can be satisfied when the strength of the random excitation is large enough, and chaotic response may arise in the system under the single Gaussian white noise excitation. Hence, the boundary of the safe basin may appear to be incursively fractal also. Fig. 5 testifies this point, and the area of the safe basin decreases following the increase of the strength of the Gaussian white noise excitation.
Fig. 4. The dispersive role of the Gaussion white noise on the systemÕs dynamics, where l = 0.0, f1 = 0.1, D = 1.0.
Fig. 5. The boundary of safe basin when the system is excited only by the Gaussian white noise, where l = 0.0, f1 = 0.0, D = 1.0.
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929
925
4. Leading Lyapunov exponent In Section 3, it is interesting to find that the boundary of the safe basin of the system can be incursively fractal when the Gaussian white noise excitation is added to the system. From the theory of dynamical systems, chaotic response must arise in the system. In this section, RosensteinÕs approach [21] is used to compute the responseÕs leading Lyapunov exponent to identify the noise-induced chaos according to the Melnikov condition presented in Section 2. The basic algorithm for the computation of the responseÕs leading Lyapunov exponent k1 is listed in [21], the following figures (Figs. 6–9) display some numerical results for the time series of system (1). Fig. 6 shows the simulation results for period-3 and chaotic motions in the deterministic case of system (1) (f2 = 0), where the embedding dimension m is chosen to be 3 due to TakenÕs theorem, i.e., m > 2n (n is the degree-of-freedom of the system). In Fig. 6(a) and (b), the systemÕs response is harmonic (Fig. 6(a)) and its Leading Lyapunov exponents is approximately equal to zero by the least-square fit (Fig. 6(b)). While f1 = 1.56, an obvious slope appears in the curve shown in Fig. 6(d), which means a positive leading exponent for this chaotic motion (Fig. 6(c)). After a short transition, there is a long linear region that is used to extract the leading Lyapunov exponent, and the curve will saturate at long times since the system is bounded in phase space and the average divergence cannot exceed the ‘‘length’’ of the attractor. These data are not shown in the figures. It is important to notice that the scaling region flattens with increasing embedding dimension m, which means m cannot be set to too large, i.e., one cannot distinguish low-dimensional systems from high-dimensional ones. Therefore, in the results shown in Fig. 7(b) and (d), Fig. 8(b) and (d) and Fig. 9(b) and (d), we all set m = 3. Moreover, the accurate values of the leading Lyapunov exponents are not accurately calculated, the main goal of this paper is just to show nonchaotic and chaotic responses in stochastic system (1). -2.50
1.0
Displacement (m)
-2.55
0.5
0.0
-2.60
-2.65
-0.5 -2.70
200
300
400
500
600
700
800
0.0
0.5
1.0
1.5
2.0
2.5
3.0
Time (s)
Time (s)
(b) Leading exponent for the response in (a)
(a) Period–3 response (f1 = 1.545) 1.0
0.5
Displacement (m)
-1.4
0.0
-1.5
-1.6
-1.7 -0.5 200
300
400
500
600
700
800
0.0
0.5
1.0
1.5
2.0
2.5
Time (s)
Time (s)
(a) Chaotic response (f1 = 1.56)
(d) Leading exponent for the response in (c)
Fig. 6. Harmonic and chaotic responses in system (1), where l = 1.0, f2 = 0.0, x = 1.0.
3.0
926
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929 -2.05
0.8
0.6
Displacement (m)
-2.10 0.4
0.2
0.0
-2.15
-2.20
-0.2
-0.4
-2.25 200
300
400
500
600
700
0.0
800
0.5
1.0
1.5
2.0
2.5
3.0
Time (s)
Time (s)
(b) Leading exponent for the response in (a)
(a) Almost–harmonic response (f1 = 1.5,, D = 0.01) 1.0
-0.8
0.5
-0.6
0.0
-1.0
-1.2
-1.4
-0.5 200
300
400
500
600
700
800
0.0
0.5
1.0
1.5
2.0
2.5
Time (s)
Time (s)
(c) Chaotic response (f1 = 1.5, D = 1.0)
(d) Leading exponent for the response in (c)
3.0
Fig. 7. Almost-harmonic and chaotic responses in system (1), where l = 1.0, f2 = 0.1, x = 1.0.
From the simulation result shown in Fig. 7(c) and (d), the harmonic response in the deterministic case will change to be chaotic when the Gaussian white noise excitation is added to the system. This can be learnt from the fact that obvious slope appears in the curve of the mean divergence of the response. Also, one can easily check that the Melnikov necessary condition (6) is satisfied for the parametersÕ values taken in Fig. 7. Here, the response in Fig. 7(a) is called almost-harmonic because its leading Lyapunov exponent equals approximately zero. To answer if chaos can arise when system (1) is excited only by the Gaussian white noise, i.e., f1 = 0, f2 5 0, we first check the Melnikov necessary condition (6) given in Section 2. In this case, only the damping and the random items appear in the Melnikov integral. This means that when the variance of the random part is enough large (see Inequality (6)), the stable and the unstable manifolds of system (1) will intersect with each other very often but randomly, and chaos may arise in this case. Fig. 8 displays the results when l = 0 in the system. For weak noise excitation, regular nature dominates the systemÕs response with flat outset in the curve, so it is not chaotic and called almost-harmonic due to its zero leading Lyapunov exponent (see Fig. 8(a) and (b)). Following the increase of the strength of the noise excitation, obvious slope appears in the curve and positive leading exponent can be obtained for this response (see Fig. 8(c) and (d)). Fig. 8 shows that the response may be chaotic or not even if Inequality (6) is satisfied. It is wonderful to learn from the simulation results shown in Fig. 9, initial ‘‘jumps’’ from a small separation occur near t = 0 which means k1 = 1 (Fig. 9(b) and (d)). If one increases the strength of the Gaussian white noise excitation, bounded solutions cannot be obtained including the chaotic ones. This means that the Gaussian white noise excitation cannot induce chaos and thoroughly random nature (k1 = 1) dominates the systemÕs responses when l = 5. As for other positive values of l, similar simulation process can be performed which is not given here.
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929
927
-4.65
0.005
Displacement (m)
0.010
0.000
-0.005
-4.70
-4.75
-0.010 -4.80
200
300
400
500
600
700
800
0.0
0.5
1.0
1.5
Time (s)
2.0
2.5
3.0
Time (s)
(b) Leading exponent for the response in (a)
(a) Almost–harmonic response (f2 = 0.01)
0.5
Displacement (m)
0.0
0.0
-0.5
-1.0
-1.5
-0.5 200
300
400
500
600
Time (s)
(c) Chaotic response (f2 = 0.2)
700
800
0.0
0.5
1.0
1.5
2.0
2.5
3.0
Time (s)
(d) Leading exponent for the response in (c)
Fig. 8. Almost-harmonic and chaotic responses in system (1), where l = 0.0, f1 = 0.0, D = 1.0.
5. Summary and discussion The present paper demonstrates the effect of the Gaussian white noise on the boundary of the safe basin and the noise-induced chaotic response in the quadratically nonlinear oscillator. By employing the stochastic Melnikov method, the necessary condition for the rising of chaos is obtained, from which one can know that the almost-harmonic, the chaotic and the thoroughly random responses can exist in the system even though the stochastic Melnikov condition is satisfied. The algorithm for leading Lyapunov exponent by Rosenstein et al. is used for the identification of noiseinduced chaos in the stochastic oscillatory system. Due to the particularity of the system, the effect of noise on the boundary of the safe basin is discussed. From the numerical results presented in Section 3, the erosion can be aggravated when the driving amplitude of the harmonic excitation, or the strength of the Gaussian white noise excitation is increased. The boundary of the safe basin can also become, by and by, fractal following the increase of the strength of the excitation when the system without damping force is excited only by the Gaussian white noise. Noise can induce chaos. When the Gaussian white noise excitation is imposed on the system, the closed loop by the stable and the unstable manifolds in the corresponding conservative system can also be broken down, but will intersect randomly (periodically in the harmonic excitation case in view of one order) with each other. When the system is excited only by the Gaussian white noise and the damping coefficient takes positive value given in Section 4, all the sample responses are not chaotic but move thoroughly randomly by the simulations of their leading Lyapunov exponents (all equal to infinity).
928
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929 0.06
-1.8
Displacement (m)
0.04
0.02
0.00
-0.02
-0.04
-2.0
-2.2
-2.4 200
300
400
500
600
700
0.0
800
0.5
1.0
1.5
2.0
2.5
3.0
Time (s)
Time (s)
(a) Thoroughly random response (f2 = 0.1, D = 1.0)
(b) Leading exponent for the response in (a) -0.4
0.6
Displacement (m)
0.4
0.2
0.0
-0.6
-0.8
-1.0
-0.2 -1.2
-0.4 200
300
400
500
600
700
800
Time (s)
(c) Thoroughly random response (f2 = 1.0, D = 1.0)
0.0
0.5
1.0
1.5
2.0
2.5
3.0
Time (s)
(d) Leading exponent for the response in (c)
Fig. 9. Thoroughly random responses in system (1) and their corresponding leading exponents, where l = 5.0, f1 = 0.0.
Acknowledgement The work was supported by the financial support from the National Natural Science Foundation of China under Grant No. 10302025.
References [1] Bulsara AR, Schieve WC, Jacobs EW. Homoclinic chaos in systems perturbed by weak Langevin noise. Phys Rev A 1990;41(2):668–81. [2] Frey M, Simiu E. Equivalence between motions with noise-induced jumps and chaos with Smale horseshoes. In: Proceedings of the 9th Engineering Mechanics Conference. ASCE; 1992. p. 660–3. [3] Frey M, Simiu E. Noise-induced chaos and phase space flux. Physica D 1993;63:321–40. [4] Xie WC. Effect of noise on chaotic motion of buckled column under periodic excitation. Nonlinear Stochastic Dyn 1994;AMD Vol. 192/DE Vol. 73:215–25. [5] Simiu E, Frey M. Melnikov processes and noise-induced exits from a well. J Eng Mech 1996;122(3):263–70. [6] Lin H, Yim SCS. Analysis of a nonlinear system exhibiting chaotic, noisy chaotic, and random behaviors. ASME J Appl Mech 1996;63:509–16. [7] Wiggins S. Global bifurcations and chaos: analytical methods. New York: Springer; 1988. [8] Soliman MS, Thompson JMT. Integrity measures quantifying the erosion of smooth and fractal basins of attraction. J Sound Vib 1989;35:453–75. [9] McDonald SW, Grebogi C, Ott E, Yorke JA. Fractal basin boundaries. Physica D 1985;17:125–53. [10] Moon FC, Li GX. Fractal basin boundaries and homoclinic orbits for periodic motion in a two-well potential. Phys Rev Lett 1985;55:1439–42.
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929
929
[11] Soliman MS. Fractal erosion of basins of attraction in coupled nonlinear systems. J Sound Vib 1995;182:727–40. [12] Senjanovic I, Parunov J, Cipric G. Safety analysis of ship rolling in rough sea. Chaos, Solitons & Fractals 1997;4:659–80. [13] Freitas MST, Viana RL, Grebogi C. Erosion of the safe basin for the transversal oscillations of a suspension bridge. Chaos, Solitons & Fractals 2003;18:829–41. [14] Xu J, Lu QS, Huang KL. Controlling erosion of safe basin in nonlinear parametrically excited systems. ACTA Mech Sinica 1996;12:281–8. [15] Gan CB, Lu QS, Huang KL. Non-stationary effects on safe basins of a softening Duffing oscillator. ACTA Mech Solida Sinica 1998;11(3):253–60. [16] Gan C. Noise-Induced chaos and basin erosion in softening Duffing oscillator. Chaos, Solitons & Fractals 2005;25:1069–81. [17] Wolf A, Swift JR, Swinney HL, Vastano JA. Determining Lyapunov exponents from a time series. Physica D 1985;16:285–317. [18] Kantz H, Schreiber T. Nonlinear time series analysis. Cambridge, England: Cambridge University Press; 1997. [19] Sano M, Sawada Y. Measurement of the Lyapunov spectrum from a chaotic time series. Phys Rev Lett 1985;55:1082–5. [20] Eckmann JP, Kamphorts SO, Rulle D, Cliliberto S. Lyapunov exponents from a time series. Phys Rev A 1986;34:4971–9. [21] Rosenstein MT, Collins JJ, Luca CJ. A practical method for calculating leading Lyapunov exponents from small data sets. Physica D 1993;65:117–34. [22] Shinozuka M. Digital simulation of random processes and its applications. J Sound Vib 1972;25(1):111–28.
Noise-induced chaos in a quadratically nonlinear oscillator Chunbiao Gan Department of Mechanics, CMEE, Zhejiang University, Hangzhou 310027, China Accepted 30 August 2005
Abstract The present paper focuses on the noise-induced chaos in a quadratically nonlinear oscillator. Simple zero points of the stochastic Melnikov integral theoretically mean the necessary rising of noise-induced chaotic response in the system based on the stochastic Melnikov method. To quantify the noise-induced chaos, the boundary of the systemÕs safe basin is firstly studied and it is shown to be incursively fractal when chaos arises. Three cases are considered in simulating the safe basin of the system, i.e., the system is excited only by the harmonic excitation, by both the harmonic and the Gaussian white noise excitations, and only by the Gaussian white noise excitation. Secondly, the leading Lyapunov exponent by RosensteinÕs algorithm is shown to quantify the chaotic nature of the sample time series of the system. The results show that the boundary of the safe basin can also be fractal even if the system is excited only by the external Gaussian white noise. Most importantly, the almost-harmonic, the noise-induced chaotic and the thoroughly random responses can be found in the system. Ó 2005 Elsevier Ltd. All rights reserved.
1. Introduction Since 1990s, stochastic Melnikov method has been applied to study the effect of noise on homoclinic or heteroclinic bifurcation and noise-induced chaos [1–6]. It is a consequence of (1) WigginsÕ [7] result that systems capable of chaotic behavior under harmonic excitation can also behave chaotically if the excitation is quasi-periodic, and (2) the application of the stochastic Melnikov approach to systems excited by physically realizable stochastic processes, which can be approximated arbitrarily closely by random quasi-periodic sums [3]. For both deterministic and stochastic systems the Melnikov conditions are necessary but not sufficient. Results yielded by the Melnikov approach are therefore weak, though less so far transient than for steady state motions. Generally, it is not easy to quantify noise-induced chaos. One knows that an important goal of studying dynamical systems is to determine their global structures, and one of these global structures is the boundary of basin. Study on safe basin in oscillatory system with potential well can be found in [8–14]. The non-stationary effect on safe basin and the analysis of noise-induced chaos in a forced softening Duffing oscillator are presented in [15,16]. Because of the coexistence of period and chaotic attractors, the basin boundaries of attractors are usually fractal and naturally incursive, in the sense that they are related to homoclinic or heteroclinic intersections of stable and unstable manifolds of the saddle points in the system, which usually means that chaos arises in such system. E-mail addresses: [email protected], [email protected] 0960-0779/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2005.08.157
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929
921
The most striking feature of chaos is the unpredictability of the future, which is called sensitive dependence on initial conditions. It can be quantified by the Lyapunov exponent, which characterizes the average exponential rate of divergence or convergence of nearby orbits in phase space of the system and quantifies the strength of chaos. Till now, most authors employed the algorithm presented by Wolf et al. [17] to compute the Lyapunov exponent even for a stochastic dynamic system. However, as pointed out by Kantz and Schreiber [18], the algorithm by Wolf is not very robust and one can easily obtain the wrong results. WolfÕs algorithm does not allow one to test for the presence of exponential divergence, but just assumes its existence and thus yields a finite exponent for stochastic data also, where the true exponent is infinite. While WolfÕs algorithm only uses a delay reconstruction of phase space, there is another class of algorithms which also involves the approximation of the underlying deterministic dynamics. References are Sano and Sawada [19], Eckmann et al. [20] and Rosenstein et al. [21]. They show that, though the output data are generally affected by noise, the influence from the noise can be minimized by using the averaging statistics. This paper studies the effect of the Gaussian white noise on the boundary of safe basin and the noise-induced chaotic response in a quadratically nonlinear oscillator. In Section 2, the stochastic Melnikov integral is analytically presented based on the dynamical theory and the nature of a linear time-invariant filter. Section 3 presents the variation of safe basin numerically when one adjusts the amplitude of the harmonic excitation in the deterministic case, or the strength of the Gaussian white noise excitation in the stochastic case of the system, from which one can see that fractal boundary appears even if the oscillator is excited only by the Gaussian white noise excitation. In Section 4, noise-induced chaotic responses are simulated and testified by employing the algorithm for the leading Lyapunov exponent from Rosenstein et al. [21]. Section 5 gives the summary and discussion.
2. Stochastic Melnikov integral and its single zero points The quadratically nonlinear oscillator considered in the present paper can be written as follows: x_ ¼ y y_ ¼ x þ x2 þ eðly þ f1 cos xt þ f2 nðtÞÞ;
ð1Þ
where 0 < e 1, l and f1 are the damping coefficient and the amplitude of the harmonic excitation respectively, and n(t) is taken as the Gaussian white noise, f2 is the strength of this stochastic excitation. The spectrum density for the Gaussian white noise is assumed to be S0 with the strength D(=2pS0). For system (1) the preliminary assumptions for the Melnikov approach are those listed by Wiggins [7] as well as the constraint for the random excitation: n(t) is uniformly continuous and bounded in ensemble, i.e., for any d1 > 0, there exists d2 > 0, such that jn(t1) n(t2)j < d1 when jt1 t2j < d2 for any times t1 and t2 and any idealization of n(t). One can easily check that system (1) satisfies the above assumptions, and the hyperbolic fixed points of system (1) will not disappear under the harmonic and the random excitations, so the Melnikov approach can also be employed to analyze the stochastic system (1). From [3,6,7], the Melnikov integral for system (1) can be expressed as Mðt0 ; l; f1 ; f2 Þ ¼ M d ðt0 ; l; f1 Þ þ M r ðt0 ; f2 Þ;
ð2Þ
where Md(t0; l, f1) is the deterministic part of M(t0; l, f1, f2), i.e., the mean value of the Melnikov integral, and Mr(t0; f2) is the corresponding random part of M(t0; l, f1, f2) due to the Gaussian white noise excitation. By the direct integration and the residue theory, one can obtain that 6 12pf1 x2 epx cos xt0 . M d ðt0 ; l; f1 Þ ¼ l þ 5 ð1 e2px Þ
ð3Þ
The random portion Mr(t0; f2) due to the Gaussian white noise, instead of directly integrated, can be estimated by considering the convolution integral as a filtering process, i.e., a stationary process n(t) passing through a linear time-invariant filter given by the homoclinic orbit [3,6]. Because of the linear time-invariant filtering, the random process Mr(t0; f2) is still stationary and of zero mean (E[Mr(t0; f2)] = 0). Its corresponding spectra may be calculated by directly multiplying the original noise spectrum S nðtþt0 Þ and the transfer function of the linear time-invariant filter. For the Gaussian white noise excitation, one has Z þ1 jF ðXÞj2 S nðtþt0 Þ ðXÞdX; ð4Þ r2M r ¼ 1
R þ1 where F ðXÞ ¼ f2 1 y 0 ðtÞeiXt dt ¼ 12pf2 X2 epX i=ð1 e2pX Þ and S nðtþt0 Þ ðXÞ S 0 is the spectrum density of the Gaussian white noise process. Integrating Eq. (4) numerically yields
922
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929
r2M r ¼ 0:764p2 S 0 f22 .
ð5Þ
When f2 5 0 in system (1), the Melnikov integral measures the random distance between the stable and the unstable manifolds. The threshold value for the rising of the chaotic response will depend on the property of the random excitation process, and may deviate from the mean value. From energy point of view [6], the necessary condition is given by 4p2 f12 x4 e2px =ð1 e2px Þ2 þ 0:764p2 S 0 f22 P
1 2 l. 25
ð6Þ
Inequality (6) tells one that the chaotic response may be induced by the Gaussian white noise excitation. Moreover, it deserves to note that inequality (6) is not sufficient, and other measures are needed for the identification of the systemÕs chaotic response. In the following two sections, the incursively fractal boundary of safe basin of the system and the positive leading exponent are presented respectively. In the present paper, we all set e = 0.1, x = 1.0 and the stepsize h = 0.01 if not pointed out specifically.
3. Incursively fractal boundary of safe basin The decrease of safe basinÕs area is often called basin erosion. Since those intersections of stable and unstable manifolds map one to another, the manifolds present a convoluted structure that extends through a wide region of the phase space. The transient orbits near such a fractal basin boundary will be as convoluted as the boundary itself. So they are expected to cross the critical line very often, increasing dramatically the number of unsafe initial conditions. In view of this point, the erosion of the safe basin is usually related to fractal basin boundary of attractors. The safe basin without erosion is drawn within the region G surrounded by the coincident stable and unstable manifolds of the conservative case of system (1) with l = f1 = f2 = 0 (see Fig. 1), which is defined by sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ( ffi ffi) 1 2 1 3 1 2 1 3 ð7Þ G ¼ ðx; yÞ 2 Rj 0:5 6 x 6 1; 2 H x þ x 6 y 6 2 H x þ x 2 3 2 3 for generating 400 grid lines in x-direction and an increasing series of grid lines from 1 to 200 when 0.5 6 x 6 0.25 and a decreasing series from 199 to 1 when 0.25 < x 6 1.0 in y-direction. Each grid point is set as an initial condition to perform the simulation of system (1). When the Hamiltonian value of any phase point (x(t), y(t)) of the systemÕs trajectory is larger than H within 100 thousand steps, this motion initiating from the corresponding initial point is considered to be unsafe and this initial point is discarded.
Fig. 1. Safe basin without erosion in the deterministic case of system (1), where l = f1 = f2 = 0.0.
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929
923
For the stochastic case of the system, 16 sample responses are simulated with the same initial condition (picked up from the grid points) and 16 different realizations of the random processes due to Monte-Carlo method [22]. Only the initial point being safe for all the 16 sample responses is recorded. 3.1. Erosion of safe basin and effect of the Gaussian white noise In Fig. 1, the safe basin is a closed and apparently smooth curve when f1 = 0. Following the increase of the driving amplitude f1, this orbit will undergo a bifurcation cascade to chaos, and the boundary of the safe basin shows to be fractal (see Fig. 2, in which H is set to be 1/6). In general, the area of the safe basin decreases following the decrease of the Hamiltonian H. To learn the effect of the noise excitation on the erosion of safe basin, Fig. 3 presents a sequence of safe basins when the Gaussian white noise is added to the system, i.e., f2 5 0. The incursively fractal fingers are also observed, which means that chaotic responses still exist in the stochastic case of system (1). From Fig. 3, one knows that the Gaussian
Fig. 2. The boundary of safe basin in the deterministic case of system (1), where l = 1.0, f2 = 0.
Fig. 3. The effect of the Gaussian white noise on the boundary of safe basin in system (1), where l = 1.0, f1 = 1.0, D = 1.0.
924
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929
white noise excitation can aggravate the erosion of the systemÕs safe basin. In addition, fractal boundary can also appear when the Hamiltonian H is taken as other positive value, which is not given here. From Fig. 4, one knows that the Gaussian white noise excitation can play a dispersive role to the systemÕs dynamics, the response of the system can still remain in the safe basin within 100 thousand steps even if the system without damping item is driven by a periodical force. In this case, incursively fractal boundary of the safe basin also means that noise excitation can induce chaos, which can be further shown in the following studies. 3.2. Noise-induced fractal boundary of safe basin From Figs. 3 and 4, the response initiating from the potential well may not escape from the well within 100 thousand steps when the system is excited by both the harmonic and the Gaussian white noise excitation, and the boundary of the safe basin can show to be incursively fractal. According to inequality (6) given in Section 2, the Melnikov condition can be satisfied when the strength of the random excitation is large enough, and chaotic response may arise in the system under the single Gaussian white noise excitation. Hence, the boundary of the safe basin may appear to be incursively fractal also. Fig. 5 testifies this point, and the area of the safe basin decreases following the increase of the strength of the Gaussian white noise excitation.
Fig. 4. The dispersive role of the Gaussion white noise on the systemÕs dynamics, where l = 0.0, f1 = 0.1, D = 1.0.
Fig. 5. The boundary of safe basin when the system is excited only by the Gaussian white noise, where l = 0.0, f1 = 0.0, D = 1.0.
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929
925
4. Leading Lyapunov exponent In Section 3, it is interesting to find that the boundary of the safe basin of the system can be incursively fractal when the Gaussian white noise excitation is added to the system. From the theory of dynamical systems, chaotic response must arise in the system. In this section, RosensteinÕs approach [21] is used to compute the responseÕs leading Lyapunov exponent to identify the noise-induced chaos according to the Melnikov condition presented in Section 2. The basic algorithm for the computation of the responseÕs leading Lyapunov exponent k1 is listed in [21], the following figures (Figs. 6–9) display some numerical results for the time series of system (1). Fig. 6 shows the simulation results for period-3 and chaotic motions in the deterministic case of system (1) (f2 = 0), where the embedding dimension m is chosen to be 3 due to TakenÕs theorem, i.e., m > 2n (n is the degree-of-freedom of the system). In Fig. 6(a) and (b), the systemÕs response is harmonic (Fig. 6(a)) and its Leading Lyapunov exponents is approximately equal to zero by the least-square fit (Fig. 6(b)). While f1 = 1.56, an obvious slope appears in the curve shown in Fig. 6(d), which means a positive leading exponent for this chaotic motion (Fig. 6(c)). After a short transition, there is a long linear region that is used to extract the leading Lyapunov exponent, and the curve will saturate at long times since the system is bounded in phase space and the average divergence cannot exceed the ‘‘length’’ of the attractor. These data are not shown in the figures. It is important to notice that the scaling region flattens with increasing embedding dimension m, which means m cannot be set to too large, i.e., one cannot distinguish low-dimensional systems from high-dimensional ones. Therefore, in the results shown in Fig. 7(b) and (d), Fig. 8(b) and (d) and Fig. 9(b) and (d), we all set m = 3. Moreover, the accurate values of the leading Lyapunov exponents are not accurately calculated, the main goal of this paper is just to show nonchaotic and chaotic responses in stochastic system (1). -2.50
1.0
Displacement (m)
-2.55
0.5
0.0
-2.60
-2.65
-0.5 -2.70
200
300
400
500
600
700
800
0.0
0.5
1.0
1.5
2.0
2.5
3.0
Time (s)
Time (s)
(b) Leading exponent for the response in (a)
(a) Period–3 response (f1 = 1.545) 1.0
0.5
Displacement (m)
-1.4
0.0
-1.5
-1.6
-1.7 -0.5 200
300
400
500
600
700
800
0.0
0.5
1.0
1.5
2.0
2.5
Time (s)
Time (s)
(a) Chaotic response (f1 = 1.56)
(d) Leading exponent for the response in (c)
Fig. 6. Harmonic and chaotic responses in system (1), where l = 1.0, f2 = 0.0, x = 1.0.
3.0
926
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929 -2.05
0.8
0.6
Displacement (m)
-2.10 0.4
0.2
0.0
-2.15
-2.20
-0.2
-0.4
-2.25 200
300
400
500
600
700
0.0
800
0.5
1.0
1.5
2.0
2.5
3.0
Time (s)
Time (s)
(b) Leading exponent for the response in (a)
(a) Almost–harmonic response (f1 = 1.5,, D = 0.01) 1.0
-0.8
0.5
-0.6
0.0
-1.0
-1.2
-1.4
-0.5 200
300
400
500
600
700
800
0.0
0.5
1.0
1.5
2.0
2.5
Time (s)
Time (s)
(c) Chaotic response (f1 = 1.5, D = 1.0)
(d) Leading exponent for the response in (c)
3.0
Fig. 7. Almost-harmonic and chaotic responses in system (1), where l = 1.0, f2 = 0.1, x = 1.0.
From the simulation result shown in Fig. 7(c) and (d), the harmonic response in the deterministic case will change to be chaotic when the Gaussian white noise excitation is added to the system. This can be learnt from the fact that obvious slope appears in the curve of the mean divergence of the response. Also, one can easily check that the Melnikov necessary condition (6) is satisfied for the parametersÕ values taken in Fig. 7. Here, the response in Fig. 7(a) is called almost-harmonic because its leading Lyapunov exponent equals approximately zero. To answer if chaos can arise when system (1) is excited only by the Gaussian white noise, i.e., f1 = 0, f2 5 0, we first check the Melnikov necessary condition (6) given in Section 2. In this case, only the damping and the random items appear in the Melnikov integral. This means that when the variance of the random part is enough large (see Inequality (6)), the stable and the unstable manifolds of system (1) will intersect with each other very often but randomly, and chaos may arise in this case. Fig. 8 displays the results when l = 0 in the system. For weak noise excitation, regular nature dominates the systemÕs response with flat outset in the curve, so it is not chaotic and called almost-harmonic due to its zero leading Lyapunov exponent (see Fig. 8(a) and (b)). Following the increase of the strength of the noise excitation, obvious slope appears in the curve and positive leading exponent can be obtained for this response (see Fig. 8(c) and (d)). Fig. 8 shows that the response may be chaotic or not even if Inequality (6) is satisfied. It is wonderful to learn from the simulation results shown in Fig. 9, initial ‘‘jumps’’ from a small separation occur near t = 0 which means k1 = 1 (Fig. 9(b) and (d)). If one increases the strength of the Gaussian white noise excitation, bounded solutions cannot be obtained including the chaotic ones. This means that the Gaussian white noise excitation cannot induce chaos and thoroughly random nature (k1 = 1) dominates the systemÕs responses when l = 5. As for other positive values of l, similar simulation process can be performed which is not given here.
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929
927
-4.65
0.005
Displacement (m)
0.010
0.000
-0.005
-4.70
-4.75
-0.010 -4.80
200
300
400
500
600
700
800
0.0
0.5
1.0
1.5
Time (s)
2.0
2.5
3.0
Time (s)
(b) Leading exponent for the response in (a)
(a) Almost–harmonic response (f2 = 0.01)
0.5
Displacement (m)
0.0
0.0
-0.5
-1.0
-1.5
-0.5 200
300
400
500
600
Time (s)
(c) Chaotic response (f2 = 0.2)
700
800
0.0
0.5
1.0
1.5
2.0
2.5
3.0
Time (s)
(d) Leading exponent for the response in (c)
Fig. 8. Almost-harmonic and chaotic responses in system (1), where l = 0.0, f1 = 0.0, D = 1.0.
5. Summary and discussion The present paper demonstrates the effect of the Gaussian white noise on the boundary of the safe basin and the noise-induced chaotic response in the quadratically nonlinear oscillator. By employing the stochastic Melnikov method, the necessary condition for the rising of chaos is obtained, from which one can know that the almost-harmonic, the chaotic and the thoroughly random responses can exist in the system even though the stochastic Melnikov condition is satisfied. The algorithm for leading Lyapunov exponent by Rosenstein et al. is used for the identification of noiseinduced chaos in the stochastic oscillatory system. Due to the particularity of the system, the effect of noise on the boundary of the safe basin is discussed. From the numerical results presented in Section 3, the erosion can be aggravated when the driving amplitude of the harmonic excitation, or the strength of the Gaussian white noise excitation is increased. The boundary of the safe basin can also become, by and by, fractal following the increase of the strength of the excitation when the system without damping force is excited only by the Gaussian white noise. Noise can induce chaos. When the Gaussian white noise excitation is imposed on the system, the closed loop by the stable and the unstable manifolds in the corresponding conservative system can also be broken down, but will intersect randomly (periodically in the harmonic excitation case in view of one order) with each other. When the system is excited only by the Gaussian white noise and the damping coefficient takes positive value given in Section 4, all the sample responses are not chaotic but move thoroughly randomly by the simulations of their leading Lyapunov exponents (all equal to infinity).
928
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929 0.06
-1.8
Displacement (m)
0.04
0.02
0.00
-0.02
-0.04
-2.0
-2.2
-2.4 200
300
400
500
600
700
0.0
800
0.5
1.0
1.5
2.0
2.5
3.0
Time (s)
Time (s)
(a) Thoroughly random response (f2 = 0.1, D = 1.0)
(b) Leading exponent for the response in (a) -0.4
0.6
Displacement (m)
0.4
0.2
0.0
-0.6
-0.8
-1.0
-0.2 -1.2
-0.4 200
300
400
500
600
700
800
Time (s)
(c) Thoroughly random response (f2 = 1.0, D = 1.0)
0.0
0.5
1.0
1.5
2.0
2.5
3.0
Time (s)
(d) Leading exponent for the response in (c)
Fig. 9. Thoroughly random responses in system (1) and their corresponding leading exponents, where l = 5.0, f1 = 0.0.
Acknowledgement The work was supported by the financial support from the National Natural Science Foundation of China under Grant No. 10302025.
References [1] Bulsara AR, Schieve WC, Jacobs EW. Homoclinic chaos in systems perturbed by weak Langevin noise. Phys Rev A 1990;41(2):668–81. [2] Frey M, Simiu E. Equivalence between motions with noise-induced jumps and chaos with Smale horseshoes. In: Proceedings of the 9th Engineering Mechanics Conference. ASCE; 1992. p. 660–3. [3] Frey M, Simiu E. Noise-induced chaos and phase space flux. Physica D 1993;63:321–40. [4] Xie WC. Effect of noise on chaotic motion of buckled column under periodic excitation. Nonlinear Stochastic Dyn 1994;AMD Vol. 192/DE Vol. 73:215–25. [5] Simiu E, Frey M. Melnikov processes and noise-induced exits from a well. J Eng Mech 1996;122(3):263–70. [6] Lin H, Yim SCS. Analysis of a nonlinear system exhibiting chaotic, noisy chaotic, and random behaviors. ASME J Appl Mech 1996;63:509–16. [7] Wiggins S. Global bifurcations and chaos: analytical methods. New York: Springer; 1988. [8] Soliman MS, Thompson JMT. Integrity measures quantifying the erosion of smooth and fractal basins of attraction. J Sound Vib 1989;35:453–75. [9] McDonald SW, Grebogi C, Ott E, Yorke JA. Fractal basin boundaries. Physica D 1985;17:125–53. [10] Moon FC, Li GX. Fractal basin boundaries and homoclinic orbits for periodic motion in a two-well potential. Phys Rev Lett 1985;55:1439–42.
C. Gan / Chaos, Solitons and Fractals 30 (2006) 920–929
929
[11] Soliman MS. Fractal erosion of basins of attraction in coupled nonlinear systems. J Sound Vib 1995;182:727–40. [12] Senjanovic I, Parunov J, Cipric G. Safety analysis of ship rolling in rough sea. Chaos, Solitons & Fractals 1997;4:659–80. [13] Freitas MST, Viana RL, Grebogi C. Erosion of the safe basin for the transversal oscillations of a suspension bridge. Chaos, Solitons & Fractals 2003;18:829–41. [14] Xu J, Lu QS, Huang KL. Controlling erosion of safe basin in nonlinear parametrically excited systems. ACTA Mech Sinica 1996;12:281–8. [15] Gan CB, Lu QS, Huang KL. Non-stationary effects on safe basins of a softening Duffing oscillator. ACTA Mech Solida Sinica 1998;11(3):253–60. [16] Gan C. Noise-Induced chaos and basin erosion in softening Duffing oscillator. Chaos, Solitons & Fractals 2005;25:1069–81. [17] Wolf A, Swift JR, Swinney HL, Vastano JA. Determining Lyapunov exponents from a time series. Physica D 1985;16:285–317. [18] Kantz H, Schreiber T. Nonlinear time series analysis. Cambridge, England: Cambridge University Press; 1997. [19] Sano M, Sawada Y. Measurement of the Lyapunov spectrum from a chaotic time series. Phys Rev Lett 1985;55:1082–5. [20] Eckmann JP, Kamphorts SO, Rulle D, Cliliberto S. Lyapunov exponents from a time series. Phys Rev A 1986;34:4971–9. [21] Rosenstein MT, Collins JJ, Luca CJ. A practical method for calculating leading Lyapunov exponents from small data sets. Physica D 1993;65:117–34. [22] Shinozuka M. Digital simulation of random processes and its applications. J Sound Vib 1972;25(1):111–28.