International Journal of Innovative Computing, Information and Control Volume 4, Number 1, February 2008

c ICIC International °2007 ISSN 1349-4198 pp. 0–0

CONVERGENCE PROPERTY OF NOISE-INDUCED SYNCHRONIZATION Katsutoshi Yoshida Department of Mechanical Systems Engineering Utsunomiya University 7-1-2 Yoto, Utsunomiya-shi, Tochigi 321-8585, Japan [email protected]

Yusuke Nishizawa ShinMaywa Industries, Ltd. 1-1 Shinmeiwa-cho, Takarazuka-shi, Hyogo 665-8550, Japan

Abstract. This paper describes how to characterize convergence properties of the synchronization system consisting of two identical nonlinear dynamical systems linked by a common noisy input only. We consider two types of input: the harmonic and random (HR) forcing and the narrow-band random (NR) forcing. Statistical linearization approach enables us to characterize the convergence properties while Liapunov exponents, which is a well-known identifier to provide a necessary condition of occurrence of the synchronization, fail to characterize them. The result shows that it is possible to detect the condition of slow convergence as multi-valued solutions of moment equations of the target system. Keywords: Noise, Synchronization, Stability, Statistical Equivalent Approach.

1. Introduction. The noise-induced synchronization of a dynamical system with its copies, as illustrated in Fig.1, can easily be found in nonlinear systems such as the discrete maps[1, 2], the Lorenz system[1], the Duffing oscillator[3], the single mode CO2 laser[4], and the uncoupled neurons[5]. One of the most important results of them is that the perfect synchronization may arise under some suitable conditions[4, 3, 1]. Moreover, the perfect synchronization exhibits significant degree of robustness against mismatches among the copies such as the parameters mismatch[3] and the independent random fluctuations of the copies[1, 6]. Furthermore, regarding the response of the discrete maps as a Markov process to derive the transition law of it, we have analytically shown that the perfect synchronization can be considered as an absorbing barrier of the Markov process[7]. It is also shown that our analysis is in good agreement with the Monte Carlo simulations and the experiments on the pair of multivibrators linked by a common noisy input only[7]. In these studies, however, little investigation have been done on transient behavior before reaching the perfect synchronization. In engineering applications, too slow convergence of synchronization errors would be equivalent to failing to converge, even if the synchronization is achieved mathematically. In this paper, we first investigate the transient probability densities of the synchronization system, showing that the slow convergence is related to the slow diffusions caused by multimodal probability densities. We then take a statistical linearization approach to 1

2

KATSUTOSHI YOSHIDA AND YUSUKE NISHIZAWA

Synchronized Random Outputs Vibration System Copy 1

Noise Generator

Copy 2 ・・・

・・・

Copy n Figure 1. Sketch of the noise-induced synchronization. detect this multimodal feature in analytical ways. It is shown that the slow convergence can be detected by multi-valued solutions of the moment equations.

2. Noise-induced Synchronization. We summarize the existence of the noise-induced synchronization and the necessary condition for it, and then, show numerically that very slow convergence of synchronization errors may arise under certain parameter conditions. 2.1. Existence of the Synchronization. In the discrete time case, we have already probed mathematically the existence of the perfect synchronization arising in the pair of logistic maps linked by a common noisy input only[7]. On the other hand, in the continuous time case, Maritan and Banavar[8] have reported that a pair of d-dimensional Langevin equations: x˙ = F (x) + η(t),

y˙ = F (y) + η(t).

(1)

where hηi (t)ηj (t0 )i = kδij (t − t0 ) (i, j = 1, · · · , d), can produce the perfect synchronization of the variables x and y in the following manner[8]. The joint probability density P (x, y, t) satisfies the Fokker-Plank-Kolmogorov (FPK) equation: ∑ dP =− dt i=1 d

(

∂ ∂ (Fi (x)P ) + (Fi (y)P ) ∂xi ∂yi

) (2)

+ k(∇2x + ∇2y + ∇x · ∇y )P where ∇x = [∂/∂x1 , · · · , ∂/∂xd ] and the initial condition is P (x, y, 0) = δ(x − x0 )δ(y − y 0 ). The equation can be rewritten in the form: 2d ∑ dP (r, t) =− ∂i Ji dt j=1

(3)

CONVERGENCE PROPERTY OF NOISE-INDUCED SYNCHRONIZATION

3

with r = [x1 , · · · , xd , y1 , · · · , yd ], Ji = fi P − k

2d ∑

Γij ∂j P,

j=1

{ 1 |i − j| = 0 or d, Γij = 0 otherwise, ∂ = [∂/∂x1 , · · · , ∂/∂xd , ∂/∂y1 , · · · , ∂/∂yd ], f = [F1 (x), · · · , Fd (x), F1 (y), · · · , Fd (y)]. Then, the stationary solution of the FPK equation Ps (x, y) = limt→∞ P (x, y, t) implies the relations: F (x)Ps (x, y) − k(∇x + ∇y )Ps (x, y) = 0, F (y)Ps (x, y) − k(∇x + ∇y )Ps (x, y) = 0.

(4)

Assuming the potential force F (x) = −∇x V (x), the probability densities of the position of the particle at infinite time will be Peq (x) ∝ e−βV (x) provided the system is ergodic. The unique solution of (4) normalizable can be obtained as Ps (x, y) ∝ δ(x − y) exp(βV (x)). From the above result, it can be concluded that a large class of the nonlinear vibration systems of the form (1) can produce the perfect synchronization. 2.2. Necessary Condition for the Synchronization. The necessary, but not sufficient condition of the existence of the synchronization is the negative sign of Liapunov exponents associated with the response of the system[3]. Let us consider the synchronization error ∆x = y − x of the pair (1). Substituting it into the second equation (x + ∆x)˙ = F (x) + Dx F (x)∆x + O(∆x)2 , provided ∆x ¿ 1, we obtain the first variation equation: ˙ = Dx F (x)∆x ∆x

(5)

which enables us to evaluate stability of the synchronization error ∆x, that is, ∆x is stable if the sign of Liapunov exponents from Eq.(5) is negative. 2.3. Mean Convergence Time. The theory of random dynamical systems[9] provides another approach to characterize the noise-induced perfect synchronization. Let (Ω, F, P ) be Wiener space and θt ω(·) := ω(t+·)−ω(t). A measurable function ϕ : R×Ω×Rn → Rn , (t, ω, x) 7→ ϕ(t, ω, x) is called a random dynamical system (RDS) over θ if it satisfies the cocycle property: ϕ(t + s, ω) = ϕ(t, θs ω) ◦ ϕ(s, ω). The transition operator generating stochastic processes produced by the system (1) yields a practical example of the RDS. In our application, the invariant quantity µω of transition operator ϕ which satisfies ϕ(t, ω)µω = µθt ω plays an important role. That is, the perfect synchronization occurs if the invariant µω coincides with the Dirac’s measure δω . This means that all the initial points subjected to the same sample path of excitation ω ∈ Ω converge to a simple point, in other words, they become perfectly synchronized. We numerically calculate the support

4

KATSUTOSHI YOSHIDA AND YUSUKE NISHIZAWA

of the time-varying function µω (t) s.t. limt→∞ µω (t) = µω by 1 ∑ = δ j , Xkj (n) := ϕ(tn , ωk )Xkj (0) m j=1 Xk (n) m

µω (t) ≈

µn,m ωk

(6)

where n is a discrete time t = n∆t, {Xkm (n)}M m=1 is the k-th set of M points composing numerical probability densities, and k represents the k-th sample path of the excitation η(t). To evaluates the convergence time of synchronization errors, we introduce a new quantity, { ¯ } K ¯ 1 ∑ j i hT i := min n ¯¯ sup |Xk (n) − Xk (n)| < ² (7) K k=1 1≤i,j≤M which we call a mean convergence time, where ² ¿ 1 is a criteria of the convergence. 3. Phenomenological Approach. 3.1. Slow Convergence of Synchronization Errors. We focus on a pair of identical piecewise linear systems: x¨ + cx˙ + kG(x; µ) = Q + u(t), y¨ + cy˙ + kG(y; µ) = Q + u(t).

(8)

where Q is a preload, u(t) is a random input, and the function G is a piecewise linear function defined by G(x; µ) = (x + µ) − |x + µ| + (x − µ) + |x − µ| describing a linear spring with a dead zone of the width µ.   x + µ (x ≤ −µ) G(x; µ) = 0 (−µ < x < µ)  x − µ (x ≥ µ). The synchronization system (8) produces the perfect synchronization whose errors decay deterministically. In the present work, we consider two types of the input u(t), the harmonic and random (HR) forcing: u(t) = h(t) := P cos(ωt) + sw(t),

(9)

and the narrow-band random (NR) forcing: u(t) = ν(t),

ν¨ + 2ζωn ν˙ + ωn2 ν = sw(t),

(10)

where w(t) is the standard Gaussian white noise. We choose c = 0.04, k = 1.0, Q = 0.3, µ = 0.7, ζ = 0.02 in both cases, and P = 0.2, s = 0.02 for the HR forcing, and ζ = 0.02, s = 1.414 × 10−2 ωn for the NR forcing. Fig.2 shows sample paths of the synchronization errors subjected to the HR forcing for ω = 0.81 (upper) and 1.07 (lower) respectively. There is significant difference in convergence times between the two conditions. The convergence times are t = 18973 and t = 195 respectively. This means that one spends time to converge more than 97 times as much as the other. Similar results are obtained for the NR forcing as shown in Fig.3. In this case, the convergence time is t = 2322 for ωn = 0.78 and t = 124 for ωn = 0.83.

x(t)-y(t)

CONVERGENCE PROPERTY OF NOISE-INDUCED SYNCHRONIZATION 5 4 3 2 1 0 -1 -2 -3 -4 -5 0

10000

20000

30000

40000

t 4 3

x(t)-y(t)

2 1 0 -1 -2 -3 -4 0

40

80

120

160

200

t Figure 2. Sample paths of synchronization errors as to the harmonic and random (HR) inputs for ω = 0.81 (upper) and 1.07 (lower) respectively. 4 3

x(t)-y(t)

2 1 0 -1 -2 -3 -4 0

1000

2000

3000

4000

300

400

t 4 3

x(t)-y(t)

2 1 0 -1 -2 -3 -4 0

100

200

t Figure 3. Sample paths of synchronization errors as to the narrow-band random (NR) inputs for ωn = 0.78 (upper) and 0.83 (lower) respectively.

5

6

KATSUTOSHI YOSHIDA AND YUSUKE NISHIZAWA

0 B

-0.02

A

-0.04

λ

-0.06 -0.08 -0.1 -0.12 -0.14 -0.16 -0.18 0.7

0.8

0.9

ω

1

1.1

1.2

0 B

-0.02

A

λ

-0.04 -0.06 -0.08 -0.1 -0.12 0.7

0.75

0.8

0.85

ωn

0.9

0.95

1

Figure 4. The largest Liapunov exponents of sample paths of the single system with the harmonic and random (HR) forcing (upper) and the narrow-band random (NR) forcing (lower) respectively. 3.2. Liapunov Exponents. As discussed in section 2.2, presence of the synchronization can be characterized by the largest Liapunov exponents being negative[3]. We investigate whether the exponents are also effective for evaluating the extent of convergence of the synchronization errors. Fig.4 shows the largest Liapunov exponent λ of the single system: x¨ + cx˙ + kG(x; µ) = Q + u(t)

(11)

which is one of the pair equations (8) subjected to the HR and NR forcing respectively. The parameter condition is the same as the sample processes shown in Fig.2. In case of the HR forcing (upper), we obtain λ ≈ −0.101 at ω = 0.81 (denoted as A) and λ ≈ −0.102 at ω = 1.07. This means that the both conditions produce nearly the same value of the largest Liapunov exponents while their converging times significantly differ from each other as discussed in section 3.1. Similar results are obtained for the NR forcing

CONVERGENCE PROPERTY OF NOISE-INDUCED SYNCHRONIZATION

7

case where we obtain nearly the same value λ ≈ −0.03 at both ω = 0.78 and ω = 0.83 (denoted as A and B respectively) . It is clear from these results that the largest Liapunov exponents do not relate to the convergence time. One explanation is that they are long time averages applicable to evaluating steady states, in other words, they can characterize the presence of the synchronization only. 3.3. Probability Densities. To take another approach, we focus on the transient probability densities of the single system (11), which satisfies FPK equations, for the HR forcing, ) ∂p ∂ ( s2 ∂ 2 p ∂ x p − = − − cx − k G(x + x; µ) + Q + P cos ωt p, (12) 2 2 0 ∂t 2 ∂x22 ∂x1 ∂x2 and for the NR forcing, ) ∂p s2 ∂ 2 p ∂ ∂ ( = − x p − − cx − k G(x + x ; µ) + Q p 2 2 0 1 ∂t 2 ∂ν22 ∂x1 ∂x2 ) ∂ ∂ ( − ν2 p − − 2ζωn ν2 − ωn2 ν1 p ∂ν1 ∂ν2

(13)

where x1 := x, x2 := x, ˙ ν1 := ν, ν2 := ν. ˙ To obtain direct solutions of these equations, we employ Monte-Carlo techniques. Of course, the stationary solutions of (12) and (13) are mathematically unique so that they do not have multi-valued solutions. It follows that the stationary probability densities of (11) do not depend on initial conditions. Since the same dynamical systems subjected to the same excitations can produce different outputs only if they start from distinct initial conditions, it is impossible to take anti-synchronized, or isolated, initial conditions with a unique stationary density of the form exp[−V (x)] > 0. It also seems that this is the main reason of why the Liapunov exponents, which is an average of stationary densities, fail to characterize the converging process. Therefore, we conclude that what is required to solve the problem is to investigate transient isolations of the probability densities. In order to provide physical insight of that transient isolations of the probability densities, Fig.5 shows the solutions of (12) at t = 60 × 2π/ω to demonstrate dependency on the forcing frequency ω and on the initial conditions x0 := (x1 (0), x2 (0)). All the densities are obtained from cumulative frequencies over the 10000 samples of numerical solutions of (11) starting from (x1 (0), x2 (0)) = (−0.5, 0) and (2, 2), denoted as (a) and (b) respectively. More mathematically, the initial functions are taken as p(x1 , x2 , 0) = δ(x1 + 0.5)δ(x2 ) and δ(x1 − 2)δ(x2 − 2) in the following calculations. It is obviously seen that the transient densities significantly depend on the initial conditions within the region ω ∈ [0.78, 0.9]. In case of ω = 0.76, both initial conditions (a) and (b) cause quite similar densities so that dependency on the initial conditions is hardly found in this case. As ω increases to 0.78, the outer ring found at ω = 0.76 is being replaced with the smaller ring for (a) while the previous state maintains for (b). Further increase of ω cause second change of the densities at ω = 0.92 where the outer ring of (b) at ω = 0.9 is being replaced with the smaller ring of (b) at ω = 0.92 while the densities of (a) maintains between ω = 0.9 and 0.92.

8

KATSUTOSHI YOSHIDA AND YUSUKE NISHIZAWA

ω = 0.76

(a) p(x1,x2,t)

(b)

p(x1,x2,t)

0.04

0.05

0

0 -6 -4 -2 0 x1

2

4

6

6 24 0 x2 -2 -6 -4

-6 -4 -2 0 x1

ω = 0.78

(a) p(x1,x2,t)

2

4

6

6 24 0 x2 -2 -6 -4

6

46 2 -2 0 x2 -6 -4

6

46 2 -2 0 x2 -4 -6

6

6 24 0 x2 -2 -6 -4

(b)

p(x1,x2,t)

0.08

0.05

0

0 -6 -4 -2 0 x1

2

4

6

46 2 -2 0 x2 -6 -4

-6 -4 -2 0 x1

ω = 0.90

(a) p(x1,x2,t)

2

4

(b)

p(x1,x2,t)

0.15

0.015

0

0 -6 -4 -2 0 x1

2

4

6

46 2 -2 0 x2 -4 -6

-6 -4 -2 0 x1

ω = 0.92

(a) p(x1,x2,t)

2

4

(b)

p(x1,x2,t)

0.15

0.1

0

0 -6 -4 -2 0 x1

2

4

6

6 24 0 x2 -2 -6 -4

-6 -4 -2 0 x1

2

4

Figure 5. Transient probability densities of the FPK equation of the HR forcing (12) at t = 60 × 2π/ω starting from (a) x0 = (−0.5, 0), and (b) x0 = (2, 2).

CONVERGENCE PROPERTY OF NOISE-INDUCED SYNCHRONIZATION

9

This dependency on initial conditions, in other words, transient isolations of the densities can be understood as stochastic counter parts of hysteric jump phenomena, although in contrast to the deterministic case, these isolations are transient ones and become diffused as time progressed because of the uniqueness of solutions of FPK equations. In this transient sense, we may say that the FPK equation (12) has multi-valued solutions within the hysteric region ω ∈ [0.78, 0.90]. The hysteric region may indicate the slow convergence of synchronization errors. Because, as was shown in Fig.2, the sample path at ω = 0.81 having very long anti-synchronized transient is apparently contained by the hysteric region ω ∈ [0.78, 0.90]. On the other hand, outside of this hysteric region, rapid convergence of the sample path presenting is observed at ω = 1.07. Therefore, we may conclude that it will be possible to detect the slow convergence of synchronization errors if this transiently happened hysteric region can be detected in a certain manner. Although the Monte-Carlo technique provides the robustest way to achieve this, it will consume a significant amount of computing resources. 4. Statistical Linearization Approach. In order to detect the hysteric region arising in the transient probability densities analytically, we introduce a standard statistical linearization technique. It will be shown that statistically linearized equations of moments also exhibit hysteric regions which precisely correspond to those of the transient densities. 4.1. Moment Equations. From FPK equations (12) and (13), we have the moment equations of (11) of the form:  m˙ 1 = m2 ,      m˙ 2 = −cm2 − khGi + Q + P cos(ωt), (14) σ˙11 = 2σ12 ,    σ˙12 = σ22 − cσ12 − kh(x1 − m1 )Gi,    σ˙22 = −2cσ22 − 2kh(x2 − m2 )Gi + s2 , and  m˙ 1 = m2 , m˙ 2 = −cm2 − khGi + Q + m3 ,      m˙ 3 = m4 , m˙ 4 = −2ζωn m4 − ωn2 m3 ,      σ˙11 = 2σ12 , σ˙12 = σ22 − cσ12 − kh(x1 − m1 )Gi + σ13 ,     σ˙13 = σ23 + σ14 ,    σ˙ = σ − 2ζω σ − ω 2 σ , 14 24 n 14 n 13  σ˙22 = −2cσ22 − 2kh(x2 − m2 )Gi + 2σ23 ,     σ˙23 = −cσ23 − kh(ν1 − m3 )Gi + σ33 + σ24 ,      σ˙24 = −cσ24 − kh(ν2 − m4 )Gi + σ34 − 2ζωn σ24 − ωn2 σ23 ,      σ˙33 = 2σ34 , σ˙34 = σ44 − 2ζωn σ34 − ωn2 σ33 ,    σ˙44 = −4ζωn σ44 − 2ωn2 σ34 + s2

(15)

respectively, where G = G(x1 ; µ) and mi = hXi i, σij = hXi Xj i, (X1 , X2 , X3 , X4 ) = (x1 , x2 , ν1 , ν2 ). The nonlinear terms are linearized by statistical equivalent techniques[10]

KATSUTOSHI YOSHIDA AND YUSUKE NISHIZAWA



30

B

A

20 10 0

4000 3000 2000



10

1000 0 0.7

0.8

0.9

ω

1

1.1

1.2

Figure 6. Variance of response hx21 i and mean convergence time hT i of the pair (8) subjected to the HR forcing. 2.2 A

1.8 1.4 1

1600 1200 800





B

400 0 0.7 0.75 0.8 0.85 0.9 0.95

ωn

1

Figure 7. Variance of response hx21 i and mean convergence time hT i of the pair (8) subjected to the NR forcing. as follows.

∫ hG(x; µ)i ≈

∞

G(x; µ)¯ p(x)dx, ∫ ∞ h(x − m)G(x; µ)i ≈ (x − m)G(x; µ)¯ p(x)dx

(16)

−∞

√ where p¯(x) := exp[−(x − m)2 /2σ]/ 2πσ.

(17)

−∞

4.2. Characterization of the Slow Convergence. Fig.6 shows the variance hx21 i = σ22 and the mean convergence time hT i for the HR forcing, plotted as a function of ω. The

CONVERGENCE PROPERTY OF NOISE-INDUCED SYNCHRONIZATION

11

variance is calculated as a numerical solution of the moment equations (14). The solid and broken lines indicate the forward and backward sweeps of the frequency ω respectively. The mean convergence time hT i is estimated by substituting numerical solutions of (8) starting from 5 × 5 uniform Cartesian grids on the region (x1 , x2 ) ∈ [−10, 10] × [−10, 10] into Eq.(7) where M = 5 × 5 = 25, K = 100, ² = 10−5 . The plots of hT i are saturated at the given value hT imax = 3000. It is clearly seen in Fig.6 that the mean convergence time hT i rapidly increase within the hysteric jumps of the variance of response hx21 i. This means that the hysteric jumps of variance in the statistical equivalent sense act as an identifier to indicate the slow convergence of the synchronization errors. Fig.7 shows the result for the NR forcing. Clearly, similar results are obtained where the mean convergence time takes the maximal value near the hysteric jumps of the variance. Again, the hysteric jumps of variance indicate the slow convergence. From these result, it can be concluded that the multi-valued solutions of variance arising in the statistically linearized moment equations will be an effective index to detect the slow convergence of synchronization errors analytically. Our statistically equivalent approach will save computing resources compered with direct simulations in terms of the MonteCarlo techniques.

5. Conclusion. We have investigated how to characterize the convergence of the synchronization errors of the synchronization system which consists of the pair of the piecewise linear systems subjected to the common random excitation. We first have demonstrated numerically that although the occurrence of the synchronization is mathematically guaranteed, the convergence speed of the synchronization errors significantly depends on the parameter conditions, and the speed can not be characterized by the largest Liapunov exponents which is one of the most common tools for synchronization phenomena. To characterize this dependency, we have examine the transient probability densities of the system synchronized and found that the slow convergence depends on the slow diffusions caused by the multimodal probability densities and that this multimodal feature can be detected as the multi-valued solutions of the statistically linearized moment equations. It follows that it is possible to detect the slow convergence if one can detect the multimodal feature of the probability densities. We then have considered the analytical method to detect the multimodal feature. To achieve this, we have proposed the statistical equivalent approach using the multi-valued solutions of variance calculated from the statistically linearized moment equations. The result shows that the multi-valued solutions indicates the slow convergence of synchronization errors both for the harmonic and random inputs and for the narrow-band random inputs. The above result leads to the conclusion that the multi-valued solutions of variance arising in the statistically linearized moment equations will be an effective index to detect the slow convergence of synchronization errors analytically. Although our statistically equivalent approach will save computing resources compered with the direct Monte-Carlo simulations, it contains unknown elements such as the statistical linearization errors. In view of this, the main problem in the near future would be investigating how the multiple solutions in the statistical equivalent sense relate the

12

KATSUTOSHI YOSHIDA AND YUSUKE NISHIZAWA

transient probability densities, which is one of the well-known problems which has not been fully understood yet. REFERENCES [1] R. Toral, C.R. Mirasso, E. Hern´ andez-Garc´ıa, and O. Piro: Analytical and numerical studies of noise-induced synchronization of chaotic systems, Chaos, 11–3, 665/673 (2001). [2] H. Suetani, T. Horita, and S. Mizutani: Noise-induced enhancement of fluctuation and spurious synchronization in uncoupled type-I intermittent chaotic systems, Phys. Rev. E, 69, 016219 (2004). [3] A. Stefa´ nski and T. Kapitaniak: Synchronization of mechanical systems driven by chaotic or random excitation, J. Sound and Vibration, 260, 565/576 (2003). [4] C.S. Zhou, J. Kurths, E. Allaria, S. Boccaletti, R. Meucci, and F.T. Arecchi: Constructive effects of noise in homoclinic chaotic systems, Phys. Rev. E, 67, 066220 (2003). [5] A.B. Neiman and D.F. Russell: Synchronization of noise-induced bursts in noncoupled sensory neurons, Phys. Rev. Lett., 88–13, 138103 (2002). [6] K. Yoshida and K. Sato: Noise-induced synchronization without coupling, Trans. JSME Series C (in Japanese), 70–696, 2228/2234 , Aug. 2004. [7] K. Yoshida, K. Sato, and A. Sugamaga: Noise-induced synchronization of uncoupled nonlinear systems, J. Sound and Vibration, 290, 34/47 , Feb. 2006. [8] A. Maritan and J.R. Banavar: Chaos, noise, and synchronization, Phys. Rev. Lett., 72, 1451/1454 (1994). [9] L. Arnold, G. Bleckert, and K.R. Schenk-Hoppe: The stochastic Brusselator: Parametric noise destroys Hopf bifurcation, Stochastic Dynamics, 71/92 , Springer-Verlag (1999). [10] Y.K. Lin and G.Q. Cai: Probabilistic Structural Dynamics, ch. 7, 281/304 , McGraw-Hill, (1979).