A probabilistic linearization method for non-linear systems subjected to additive and multiplicative excitations

International Journal of Non-Linear Mechanics 41 (2006) 1191 – 1205 www.elsevier.com/locate/nlm

A probabilistic linearization method for non-linear systems subjected to additive and multiplicative excitations S. Lacquaniti, G. Ricciardi ∗ Dipartimento di Ingegneria Civile (DIC), University of Messina, Italy Received 8 August 2006; accepted 5 December 2006

Abstract A general method to obtain approximate solutions for the random response of non-linear systems subjected to both additive and multiplicative Gaussian white noises is presented. Starting from the concept of linearization, the proposed method of “Probabilistic Linearization” (PL) is based on the replacement of the Fokker–Planck equation of the original non-linear system with an equivalent one relative to a linear system subjected to additive excitation only. By means of the general scheme of the weighted residuals, the unknown coefficients of the equivalent system are determined. Assuming a Gaussian probability density function of the response process and by choosing the weighting functions in a suitable way, the equivalence of the proposed method, called “Gaussian Probabilistic Linearization” (GPL), with the “Gaussian Stochastic Linearization” (GSL) applied to the coefficients of the Itô differential rule is evidenced. In addition, the generalization of the proposed method, called “Generalized Gaussian Probabilistic Linearization” (GGPL), is presented. Numerical applications show as, varying the choice of the weighting functions, it is possible to obtain different linearizations, with a variable degree of accuracy. For the two examples considered, different suitable combinations of the weighting functions lead to different equivalent linear systems, all characterized by the exact solution in terms of variance. 䉷 2007 Elsevier Ltd. All rights reserved. Keywords: Random excitations; Non-linear systems; Probabilistic linearization; Gaussian probabilistic linearization; Fokker–Planck equation

1. Introduction For non-linear systems subjected to both additive and multiplicative random excitations, the probability density function of the response is known only in very special cases [1–3]. In all the cases in which the exact solution is not available one needs to adopt approximate methods. A method to derive approximate solutions, frequently used, that gives a partial probabilistic characterization of the response is given by the moment equation approach. Since for a nonlinear system the moment differential equations constitute an infinite hierarchy, closure methods are needed to express the higher order moments in terms of lower ones. In this context, the Gaussian Closure (GC) method is the simplest and most commonly used method; it is based on the hypothesis of

∗ Corresponding author. Tel.: +39 090 3977163; fax: +39 090 3977480.

E-mail address: [email protected] (G. Ricciardi). 0020-7462/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijnonlinmec.2006.12.002

Gaussianity of the response process and it allows one to approximate the second order moments of the response. Another common method to derive approximate is based on the concept of linearization. In particular the Stochastic Linearization (SL) is the most used and versatile approximate technique in the study of non-linear random problems. There is a vast literature devoted to this subject; it is based on the replacement of the original non-linear system by an equivalent linear one, whose coefficients are obtained by minimizing the difference in a least mean square sense [4–7]. This method leads to different solutions on the basis of the particular criteria of linearization adopted. An exhaustive review of these criteria can be found in [8,9]. In this paper, after introducing some general concepts on stochastic dynamics [10,11], the application of the SL method to non-linear systems will be illustrated, showing that this approach depends on the presence or not of multiplicative excitations. In particular, in the case of multiplicative excitations, two procedures that lead to different solutions are present in the

1192

S. Lacquaniti, G. Ricciardi / International Journal of Non-Linear Mechanics 41 (2006) 1191 – 1205

literature: the linearization on the equation of motion proposed by Chang and Young [12] and two different linearizations on the Itô differential equations proposed by Brukner and Lin [13] and by Kottalam et al. [14] and Wu [15], almost simultaneously. It is important to note that the linearization method proposed by Kottalam et al. and afterwards by Wu is equivalent to a linearization on the Itô differential rule [16] directly, and consequently only in this case the SL method coincides with the GC method. In the second part of this paper, a new approach of SL is proposed. It is based on the linearization applied on Fokker–Planck equation and can be seen as a general method of approximate solution valid for non-linear systems under additive and multiplicative excitations. The proposed method, called “Probabilistic Linearization” (PL), is based on the replacement of the Fokker–Planck equation of the original non-linear system with an equivalent one relative to a linear system subjected to external excitation only, which is closed in some sense to the original system [17]. The equivalent system is determined by means of the general scheme of weighted residuals, by imposing that the error to be orthogonal to a set of linearly independent weighting functions opportunely chosen [18]. In this way, the Gaussian solution of the equivalent system may be considered as an approximate solution of the original non-linear one. Satisfying in an approximate way the Fokker–Planck equation of both the original non-linear system and the equivalent linear one, the proposed PL technique leads to a residual error expressed in terms of the unknown coefficients of the equivalent linear system, which are obtained by applying the weighted residuals technique. In this way, the equivalent coefficients depend on the probability density function and on the weighting functions adopted. As it will be evidenced, by assuming a Gaussian probability density function of the response and by choosing in a suitable way the weighting functions, the proposed method called “Gaussian Probabilistic Linearization” (GPL) coincides with the “Gaussian Closure” (GC) and the “Gaussian Stochastic Linearization” (GSL) applied to the coefficients of the Itô differential rule. More accurate solutions than GPL can be achieved by choosing judiciously the weighting functions; in these cases the generalization of the proposed method is here called “Generalized Gaussian Probabilistic Linearization” (GGPL). Unlike the classical stochastic linearization procedure, the application of the GGPL method to non-linear oscillators will evidence that as a unique (universal) linearization, valid in general dose not exist, but it is possible to perform different linearizations by choosing the weighting functions differently for the specific cases [19,20]. In particular, the numerical examples will illustrate as a suitable choice of the weighting functions leads to exact solution in terms of the variance of the response process; in this case the GGPL method is a “true” linearization of the non-linear system [21]. 2. Some preliminary concepts on stochastic dynamics Let us consider a non-linear system whose dynamic behavior, described by the n-vector of state variables X(t), is governed

by the following system of first-order differential equations:  d gik (X)Wk (t) Xi (t) = fi (X) + dt m

(i = 1, 2, . . . , n),

(1)

k=1

where Xi (t) are the components of the response state vector X(t), while Wk (t) with k = 1, 2, . . . , m are m random excitations. Functions fi (X) and gik (X) of the state vector X(t) are generally non-linear. In the following we assume that the random excitations are zero-mean Gaussian white noise processes, fully characterized by the second order cross-correlation functions E[Wk (t1 )Wl (t2 )] = Qkl (t1 − t2 ), (·) being the Dirac’s delta function, Qkl = 2Skl the constant strengths of the white noise processes, Skl the cross-power spectral density of Wk (t) and Wl (t). If gik (X) are constant, i.e. gik (X) = gik , then the system is said to be excited by additive excitations only; on the contrary, if the functions gik (X) depend on the response process X(t) the system is said to be excited by multiplicative excitations. Eq. (1) can be written in the form of Itô stochastic differential equations as follows: dXi (t) = mi (X) dt +

m 

gik (X) dk (t),

(2)

k=1

where k (t) (k = 1, 2, . . . , m) are m Wiener processes whose increments are characterized by cross-correlation functions E[dk (t1 ) dl (t2 )] = Qkl (t1 − t2 ) dt1 dt2 , while mi (X) are the drift coefficients related to the coefficients of the equations of motion (1) by the following expression: mi (X) = fi (X) + zi (X).

(3)

In particular, the second terms on the right-hand side of expression (3), given by 1  j Qkl gj l (X) gik (X) 2 jXj n

zi (X) =

m

m

(4)

j =1 k=1 l=1

are the well-known Wong–Zakai (or Stratonovich) correction terms [22]. From Eq. (3) it is easy to note that the Wong–Zakai correction terms zi (X) vanish in the case of purely additive excitations only. By introducing the standard Wiener processes Bj (t) (with j = 1, 2, . . . , n), whose increments dBj (t) are characterized by cross-correlation functions E[dBj (t1 ) dBl (t2 )] = j l (t1 − t2 ) dt1 dt2 , Eq. (2) can be rewritten in the following equivalent Itô form: dXi (t) = mi (X) dt +

n 

ij (X) dBj (t),

(5)

j =1

where ij (X) are the diffusion coefficients related to the coefficients of the equations of motion (1) by the following expressions: 2ij (X) =

m m   k=1 l=1

Qkl gik (X)gj l (X).

(6)

S. Lacquaniti, G. Ricciardi / International Journal of Non-Linear Mechanics 41 (2006) 1191 – 1205

Exact solutions for non-linear systems subjected to both additive and multiplicative random excitations are difficult to obtain. In many cases, therefore, it is necessary to adopt approximate solutions. In the case in which the response of the system can be represented by a Markov process, a frequently used approximate solution is given by the moment equation approach. This approach is based on the use of the Itô differential rule, given by: ⎧⎡ ⎫ ⎤ n ⎨ n ⎬   j(X) ⎣mi (X) dt + d(X) = ij (X) dBj (t)⎦ ⎩ jXi ⎭ +

n n 1 

2

2ij (X)

i=1 j =1

j2 (X) jXi jXj

 dt ,

(7)

where (X) is a differentiable function of the response vector process X(t). By applying the stochastic average to both members of Eq. (7) and dividing by dt, the following differential equation governing the evolution of the average of the function (X) are obtained:   n  j j(X) E mi (X) E[(X)] = jXi jt i=1

  n n 1  j2 (X) 2 + . E ij (X) 2 jXi jXj

(8)

i=1 j =1

Setting (X)=X1k1 X2k2 . . . Xnkn , the differential equation (8) rule all the moments of order k = k1 + k2 + · · · + kn of the response state vector X(t). Unfortunately, for non-linear systems, these equations constitute an infinite hierarchy; in fact the moment differential equations up to a given order contain moments of higher order and, consequently, this system is not solvable. Closure methods are needed to solve the moment differential equations, based on some hypothesis on the response process in order to express the higher order moments in terms of lower ones. The Gaussian Closure (GC) method in conjunction with the moment equations approach is the simplest and most commonly used method; it is based on the hypothesis of Gaussianity of the response process and it allows to approximate the second order moments of the response. In general, the moment equation approach, in conjunction to a closure technique, gives a partial probabilistic characterization of the response process only in terms of statistical moments until a fixed order. A complete probabilistic characterization in terms of the joint probability density function of the response process X(t) could be obtained by solving the following second-order partial differential equation, called the Fokker–Planck–Kolmogorov (FPK) equation [10,11]:  n jp NL (x, t)  j NL − [ai (x)p NL (x, t)] jt jxi i=1

+

where p NL (x, t) is the joint probability density function of the NL (x) are the firstresponse vector X(t), while aiNL (x) and bij and the second-order derivative moments, respectively, which are related to the drift and diffusion coefficients as follows: aiNL (x) = [mi (X)]X=x ⎡ ⎤ m  m n   1 j = ⎣fi (X) + Qkl gj l (X) gik (X)⎦ 2 jxj j =1 k=1 l=1

 n n 1  j2 NL [bij (x)p NL (x, t)] = 0 2 jxi jxj i=1 j =1

(9)

,

X=x

(10a)

j =1

i=1

1193

 NL bij (x) = [2ij (X)]X=x =

m  m 

 ,

Qkl gik (X)gj l (X)

k=1 l=1

X=x

(10b) where the superscript NL means that the functions are relative to the original non-linear system. Unfortunately, complete analytical solutions of this equation are known for first-order systems only [23,24]. For higher order systems, analytical solutions are known in stationary conditions, and generally for non-linear oscillators with additive excitations only [25]. The problem of obtaining exact solutions becomes much more difficult if multiplicative excitations are also present. In this case, by using the concept of “detailed balance” [26], a general procedure to obtain the exact stationary solutions has been developed [27,28]. Nevertheless, it has been shown that, when both additive and multiplicative excitations are present, certain relations between the system parameters and spectral densities of the excitations are implied in order to find exact solutions [3]. In general, such relations do not hold in practical cases; therefore, it is necessary to use approximate approaches. 3. The SL method Among various approximate methods, the well-known Stochastic Linearization (SL) procedure is the most used in solving non-linear random vibration problems. As an illustration of the method, let us consider a singledegree-of-freedom (SDOF) oscillator with non-linear stiffness and damping under additive and multiplicative excitations, whose dynamic behavior is governed by the following secondorder differential equation: ˙ = X¨ + H (X, X)

m 

˙ j (t). gj (X, X)W

(11)

j =1

The basic idea of the SL approach is to replace the original non-linear system by a linear one in such a way that the difference between the two systems is minimized in some statistical sense. This basic idea has been developed by many researchers, obtaining various versions of this method [8,9]. Since the application of this approach depends on the presence or not of multiplicative-type excitation, in the next two subsections it seems suitable to consider separately the case of

1194

S. Lacquaniti, G. Ricciardi / International Journal of Non-Linear Mechanics 41 (2006) 1191 – 1205

purely additive excitations and the case of multiplicative-type excitations. 3.1. Purely additive excitation For non-linear oscillators with random additive excitations only, the equation of motion of the system is ˙ = X¨ + H (X, X)

m 

gj Wj (t),

(12)

j =1

where gj are constant coefficients. Following the classical approach of SL method, the non-linear equation of motion (12) is replaced by the following linearized one: X¨ + eq X˙ + keq X =

m 

In the case of systems with multiplicative-type excitations, the applications of the SL method are based on various procedures that lead to different results, as shown in the next subsections. 3.2.1. Linearization on the equation of motion Chang and Young [12] apply the SL procedure in the classical form, as in the case of non-linear systems with additive excitations only. Following this approach, Eq. (11) is replaced by the “linearized” one: H ˙ X¨ + H eq X + keq X =

m 

j (jeq X˙ + keq X)Wj (t).

(15)

j =1

gj Wj (t).

(13)

j =1

˙ is replaced by a linear Thus, the non-linear function H (X, X) function that depends on the unknown equivalent coefficients keq and eq . They are found by minimizing the differences between the two systems represented by Eqs. (12) and (13) in least-square sense. This is equivalent to impose the conditions jE[e2 ]/j = 0, with  = keq , eq , where e = H (X1 , X2 ) − eq X2 − keq X1 is the difference between Eqs. (12) and (13), ˙ and X1 = X(t) and X2 = X(t). This procedure yields the wellknown expressions of the equivalent coefficients: E[H (X1 , X2 )X1 ] , keq = 2X1

3.2. Multiplicative-type excitation

E[H (X1 , X2 )X2 ] eq = , 2X2

˙ is reIn this way, not only the non-linear function H (X, X) ˙ placed by a linear one, but the m non-linear functions gj (X, X), that define the multiplicative character of the excitations, are H, linearized too. The 2m+2 unknown equivalent coefficients keq j

j H eq , keq and eq (j = 1, 2, . . . , m) are determined by minimizing the mean square of the following m + 1 errors committed passing from Eq. (10) to Eq. (15): H eH (X1 , X2 ) = H (X1 , X2 ) − H eq X2 − keq X1 , j ej (X1 , X2 ) = gj (X1 , X2 ) − jeq X2 − keq X1

(j = 1, 2, . . . , m). (14)

where 2X1 =E[X12 ] and 2X2 =E[X22 ] are the response variances. It must be stressed that Eqs. (14) cannot be evaluated, since the joint probability density function of the response processes is unknown. Various versions of SL method can be obtained on the basis of the joint probability density function used in order to calculate the unknown statistics in Eqs. (14). A consistent choice in accordance with the Gaussian solution of the equivalent linear oscillator with additive excitations (13) is that one that uses the Gaussian probability density function. On the basis of this choice, the classical SL method is wellknown as the Gaussian Stochastic Linearization (GSL). Unfortunately the GSL method gives accurate results for weakly nonlinear systems only. This drawback is due to the inadequacy of the Gaussian assumption to represent the non-Gaussian characteristic of the response for systems that exhibit strong nonlinear behavior. In these cases, a best result can be obtained by adopting a non-Gaussian joint probability density function, that allows obtaining a non-Gaussian expression of the equivalent coefficients of the linearized system; this alternative approach is well-known as non-Gaussian Stochastic Linearization (NGSL) [29–32]. It is important to point out that, as before illustrated, in the case of additive excitations only, the classical SL procedure leads to the same results of the GC if it is directly applied on the coefficients of the non-linear equation of motion [7,33].

(16a)

(16b)

In the stationary state, one find for the equivalent coefficients H , H expressions similar to that given by Eq. (14), obtained keq eq in the case of purely additive excitations, while the 2m unknown j equivalent coefficients keq and jeq are given as follows: j = keq

E[gj (X1 , X2 )X1 ] , 2X1

jeq =

E[gj (X1 , X2 )X2 ] . 2X2

(17)

A possible choice in order to calculate the averages at the righthand side of these equations could be to adopt the probability density function relative to the “quasi linear” oscillator (15); unfortunately, the analytical expression of this function is unknown. Alternatively, it is possible to assume a Gaussian probability density function but the solution achieved by means of this approach is different from that obtained by applying the GC method. 3.2.2. Linearization on the Itô-type differential equations An alternative version of SL method in the case of multiplicative excitations could consist in linearizing the following Itô equations associated to the equation of motion of the nonlinear oscillator (10): dX1 = X2 dt, dX2 = m(X1 , X2 ) dt +

(18a) m  j =1

gj (X1 , X2 ) dj (t),

(18b)

S. Lacquaniti, G. Ricciardi / International Journal of Non-Linear Mechanics 41 (2006) 1191 – 1205

where

Therefore, the method consists in the linearization of the drift coefficient m(X1 , X2 ) and the replacement of the diffusion term (X1 , X2 ) by the square root of a second degree polynomial. In particular, by minimizing the mean square of the error e1 , given by Eq. (21a), with respect to keq and eq , and of the following one:

m(X1 , X2 ) = − H (X1 , X2 ) 1  j Qj l gl (X1 , X2 ) gj (X1 , X2 ) 2 jX2 m

+

m

j =1 l=1

(19) is the drift coefficient which contains the Wong–Zakai correction terms. Following this approach, Eqs. (18) is replaced by the following linearized ones: dX1 = X2 dt,

(20a)

dX2 = (keq X1 + eq X2 ) dt +

m 

bj dj (t),

(20b)

j =1

where bj are m additional equivalent coefficients. Minimizing the mean square of the following errors committed passing from Eq. (18) to Eq. (20): e1 (X1 , X2 ) = m(X1 , X2 ) − keq X1 − eq X2 ,

(21a)

e2,j (X1 , X2 ) = gj (X1 , X2 ) − bj

(21b)

(j = 1, 2, . . . , m)

one obtain the following expressions for the 2 + m equivalent coefficients: keq =

E[m(X1 , X2 )X1 ] , 2X1

eq =

(22)

dX1 = X2 dt,

(23a)

dX2 = m(X1 , X2 ) dt + (X1 , X2 ) dB(t),

(23b)

⎡

(X1 , X2 ) = ⎣

m m  

⎤1/2 j l (X1 , X2 )⎦

j =1 l=1

⎡ =⎣

m m  

⎤1/2 Qj l gj (X1 , X2 )gl (X1 , X2 )⎦

(24)

j =1 l=1

is the diffusion term. Eqs. (23) are replaced by the following “linearized” ones: dX1 = X2 dt,

(25a)

dX2 = (keq X1 + eq X2 ) dt + (c11 X12 + c12 X1 X2 + c22 X22

+ d 1 X1 + d 2 X2 )

1/2

dB(t).

− (2eq + c12 )X1 X2 − (2keq + c22 )X22 − d 1 X1 − d 2 X 2

(26)

with respect to dr and crs (r, s =1, 2), the unknown coefficients of the “linearized” oscillator represented by Eqs. (25) are determined. Although the explicit expressions of the obtained equivalent coefficients involve statistical moments of order higher than two, this approach leads to differential equations of the first and second order moments that are identical for the original and the linearized systems. Nevertheless also, this approach leads to different results with respect to GC. Kottalam et al. [14] and Wu [15], almost simultaneously, have proposed a modified version of the SL procedure that leads to the same results of the GC approach. Besides the error e1 committed on the drift coefficients, they introduce a new expression of error measured on each component of the diffusion terms j l (X1 , X2 ), given by the following expression: (j, l = 1, 2, . . . , m)

It is worth noting that this approach leads to different results to that determined by the GC. A different version of the approach based on the linearization of the Itô equations has been proposed by Bruckner and Lin [13]. Following this approach, the Itô equations (18) are transformed in a form involving standard Wiener process B(t), as follows:

where

e2 (X1 , X2 ) = 2X2 m(X1 , X2 ) + 2 (X1 , X2 ) − c11 X12

e¯2,j l (X1 , X2 ) = Qj l gj (X1 , X2 )gl (X1 , X2 ) − Qj l b¯j l

E[m(X1 , X2 )X2 ] , 2X2

bj = E[gj (X1 , X2 )].

1195

(27)

where b¯j l = bj bl . In this way, by minimizing the mean square error with respect to the unknown coefficients b¯j l , we find b¯j l = E[gj (X1 , X2 )gl (X1 , X2 )]

(j, l = 1, 2, . . . , m)

(28)

It is important to note that Eq. (27) represents an expression of the error committed passing from each diffusion term j l relative to the non-linear oscillator, whose Itô equations are given by Eqs. (18), to each diffusion coefficient Qj l b¯j l relative to the equivalent linear oscillator subjected to additive excitations only, whose Itô equations are given by Eqs. (20). Although they describe the method as applying SL to the coefficients of the Itô equations, it may be more appropriate to say that they apply the SL procedure to the coefficients of the Itô differential rule, directly. Following this approach proposed by Falsone [16], let (X1 , X2 ) be a differentiable function of the response processes, the Itô differential rule corresponding to the non-linear oscillator (11): ⎛ j(X1 , X2 ) j(X1 , X2 ) X2 dt + d(X1 , X2 ) = ⎝ jX1 jX2 ⎡ ⎤⎞ m  × ⎣m(X1 , X2 ) dt + gj (X1 , X2 ) dj (t)⎦⎠ j =1

⎞ m m 1 ⎝ j2 (X1 , X2 )   2 + j l (X1 , X2 )⎠ dt 2 jX22 ⎛

j =1 l=1

(25b)

(29)

1196

S. Lacquaniti, G. Ricciardi / International Journal of Non-Linear Mechanics 41 (2006) 1191 – 1205

is replaced by the following linearized one: ⎛ j(X1 , X2 ) j(X1 , X2 ) d(X1 , X2 ) = ⎝ X2 dt + jX1 jX2 ⎡ × ⎣(keq X1 + eq X2 ) dt +

m 

⎤⎞ bj dj (t)⎦⎠

NL (x) and p L (x) be the joint probability Let, moreover, pex density functions of the response process X(t); in particular, these functions are the stationary solutions of the following reduced FPK equations corresponding to the original non-linear system and to the equivalent linear one, respectively: ex LNL FP [p (x)] =

j =1

⎛

i=1

⎞

1 j2 (X1 , X2 )   ¯ ⎠ + ⎝ bij dt. 2 jX22 m

m

 n n 1  j2 NL ex − [b (x)p (x)] = 0, 2 jxi jxj ij

(30)

i=1 j =1

j =1 l=1

Comparing Eq. (29) with Eq. (30), we find the same errors represented by Eq. (21a) on the drift coefficients and by Eq. (27) on each component of the diffusion coefficient; by minimizing the expression (21a) with respect to keq and eq , and the expressions (27) with respect to the unknowns b¯j l , we obtain the expressions (22a) and (28) for the equivalent coefficients. By following this procedure, in the case of non-linear systems subjected to multiplicative excitations, it is possible to affirm that SL and GC are two equivalent approaches even if SL is directly applied to the coefficients of the Itô differential rule. As was evidenced in the previous sections, in the case of multiplicative-type excitations, the SL method shows different versions depending on the linearization procedure that can be performed on the equation of motion, or on the Itô-type differential equations, or, finally, on the coefficients of the Itô differential rule. In the next section a new approach of SL based on the FPK equation will be proposed, which can be seen like a general method of approximate solution for non-linear systems subjected to both additive and multiplicative random excitations. 4. A PL method Unlike the classical SL method, which minimizes the mean value of an expression of error between the non-linear and linear systems measured on the differential equations of motion, the proposed method minimizes an expression of residual error measured on the associated FPK equations. L To illustrate this approach, let LNL FP (x) and LFP (x) be the forward FPK differential operators corresponding to the nonlinear system represented by Eq. (1), and the equivalent linear one under additive Gaussian white noise excitations only, whose equation of motion is given by the following expression:   d W¯ ik (t), Xi = ij Xj + dt n

m

j =1

k=1

 n  j NL ex [a (x)p (x)] jxi i

(31)

where ij , (for i, j =1, 2, . . . , n) are n2 unknown equivalent coefficients, while W¯ ik (t), (for i =1, 2, . . . , n and k =1, 2, . . . , m) are m × n “equivalent” Gaussian white noise processes with j ¯ ij (t1 − t2 ). cross-correlation functions: E[W¯ ik (t1 )W¯ l (t2 )] = Q kl

(32) LLFP [p L (x)] =

 n  j L L [a (x)p (x)] jxi i i=1

 n n 1  j2 L L − [b (x)p (x)] = 0. 2 jxi jxj ij i=1 j =1

(33) NL (x) are given by Eqs. (10a) In Eq. (32), aiNL (x) and bij L (x) and (10b), respectively, while in Eq. (33) aiL (x) and bij are the first and second derivative moments corresponding to the equivalent linear system (31), given by the following expressions:

aiL (x) = aiL (x, i ) =

n 

ij xj ,

j =1 L L ¯ ij ) = bij (x) = bij (x, Q

m m  

¯ ij Q kl

(34a,b)

k=1 l=1

where (for fixed i = 1, 2, . . . , n), ij are the components of the ¯ ij are the n-vector i , while (for fixed i, j = 1, 2, . . . , n), Q kl ¯ ij of order m. Since Eq. (33) components of the square matrix Q corresponds to the equivalent linear system (31) subjected to ¯ ij ) is a Gaussian additive excitations only, pL (x) = p G (x, i , Q probability density function, which depends on the unknown ¯ ij , parameters collected in the vectors i and in the matrices Q which appear in the FPK equation (33). The proposed Probabilistic Linearization (PL) method se¯ ij and, consequently, the lects the vectors i and the matrices Q equivalent linear system (31), such that the corresponding FPK equation, represented by Eq. (33), is the closest in some sense to the FPK equation (32) of the original non-linear system [17]. In particular, the proposed approach finds the equivalent linear ¯ ij ) best system, whose probability density function p G (x, i , Q minimizes an expression of residual error measured on the FPK equations. In order to illustrate the procedure, we consider a joint probability density function p(x) ˜ that satisfies in approximate way the FPK equation of both the non-linear system and the equivalent linear one; in fact, the substitution of p(x) ˜ into Eqs. (32)

S. Lacquaniti, G. Ricciardi / International Journal of Non-Linear Mechanics 41 (2006) 1191 – 1205

and (33) gives the following two residual errors:  n  j NL (NL) (x) = LNL [ p(x)] ˜ = [a (x) p(x)] ˜ FP jxi i i=1

−

 n n 1  j2 NL [bij (x)p(x)] ˜  = 0, 2 jxi jxj

(35)

i=1 j =1



(L)

¯ = LLFP [p(x)] (x, , Q) ˜ =

 n  j L i [a (x,  )p(x)] ˜ jxi i

 n n 1  j2 L ¯ ij )p(x)] [bij (x, Q ˜ jxi jxj 2 i=1 j =1

= 0,

(36)

where  = [1 2 . . . n ] is the matrix collecting the vectors ¯ is a block matrix, whose (i, j ) block is given by the i and Q ¯ ij . matrix Q ¯ (as a funcThe two residual errors (NL) (x) and (L) (x, , Q) tion of x), represent indicators of departure of p(x) ˜ from the NL (x) and p L (x). Subtracttwo probability density functions pex ing Eq. (36) from Eq. (35), we find the following expression of the global residual error that will be used afterwards: ¯ = (NL) (x) − (L) (x, , Q) ¯ = L[p(x)], (x, , Q) ˜

(37)

where  n  j i L(·) = [ai (x,  )(·)] jxi i=1

 n n 1  j2 ij ¯ − [bij (x, Q )(·)] 2 jxi jxj

(38)

i=1 j =1

ai (x, i ) = aiNL (x) − aiL (x, i )

(39)

NL L ¯ ij ) = bij ¯ ij ) bij (x, Q (x) − bij (x, Q

(40)

represent the errors in terms of the first and second derivative moments, respectively. The global residual error given by Eq. (37) can be written in the following form: n  i=1

1   (2) ¯ ij ), ij (x, Q 2 n

(1) i (x, i ) −

n

(41)

i=1 j =1

where (1) i (x, i ) =

j [ai (x, i )p(x)] ˜  = 0, jxi

(2) ¯ ij ) = ij (x, Q

In order to select the equivalent linear system, the method of weighted residuals will be used here. This is an approximate approach that consists in a general minimization scheme in which the equivalent linear system is chosen such that the following integral:  ¯ ¯ dx = 0

wε (, Q) = wε (x)(x, , Q) Rn

(ε = 1, 2, . . . , R)

(44)

vanishes for a number R of selected weighting functions wε (x). Since the global residual error (41) can be splitted in two parts, the procedure of minimization will be performed on each of these parts, separately. In this way, the vectors i of the un¯ ij of the unknown equivalent coefficients and the matrices Q known strengths of the equivalent Gaussian white noises are obtained by imposing that every component is orthogonal with respect to a set of linearly independent weighting functions, suitably selected:  (1) (1)

wr ,i (i ) = wr (x)i (x, i ) dx = 0 Rn

is a particular differential operator obtained subtracting the FPK operators for the non-linear system and the equivalent linear ¯ ij ), given by: one. In Eq. (38) ai (x, i ) and bij (x, Q

¯ = (x, , Q)

are n + n2 components that contribute to the global residual error. From expressions (42) and (43) it is possible to note that (1) the first component i is a function of the vector i only, (2) whereas the second component ij is a function of the matrix ¯ ij only. Q In the next section, a procedure to minimize each residual error component, and consequently the global one, will be illustrated. 5. Minimization of the components of the global residual error

i=1

−

1197

j2 ¯ ij )p(x)] [bij (x, Q ˜ = 0 jxi jxj

(r = 1, 2, . . . , n),  (2) (2) ¯ ij ) = ¯ ij ) dx = 0

wuv ,ij (Q wuv (x)ij (x, Q

(45a)

(u, v = 1, 2, . . . , m).

(45b)

Rn

In Eqs. (45a,b), the subscript w of indicates that depends on the particular choice of the weighting functions. For the present purpose, we assume that the weighting functions wr (x) and wuv (x) are differentiable with respect to xi . Upon integration by parts on every xi and using the expressions (42) and (43), Eqs. (45a,b) will be written as follows:  jwr (x) (1)

wr ,i (i ) = ai (x, i ) p(x) ˜ dx = 0 jxi Rn (r = 1, 2, . . . , n)

(42) (2) ¯ ij ) =

wuv ,ij (Q

(43)



(46a)

¯ ij ) j wuv (x) p(x) dx = 0 bij (x, Q jxi jxj Rn

(u, v = 1, 2, . . . , m)

2

(46b)

1198

S. Lacquaniti, G. Ricciardi / International Journal of Non-Linear Mechanics 41 (2006) 1191 – 1205

By selecting a set of functions wr (x) and wuv (x), one obtains a corresponding set of constraints of the type of Eqs. (46a,b). In particular, for fixed i = 1, 2, . . . , n and for j = 1, 2, . . . , n the n equivalent coefficients ij collected in the vector i can be determined from expressions (46a) only, whereas, for fixed i, j = 1, 2, . . . , n and for k, l = 1, 2, . . . , m, the remaining m2 ¯ ij , collected in the matrix Q ¯ ij , can be unknown strengths Q kl determined in similar way from expression (46b). In fact, with reference to the equivalent coefficients ij , by substituting Eqs. (9a), (34a) and (39) in Eq. (46a), the latter becomes: ⎞ ⎤ ⎡⎛ n  jw (X) r (1) ⎦=0

wr ,i (i ) = E˜ ⎣⎝mi (X) − ij Xj ⎠ jXi

¯ ij ). By assuming m2 represents each component of ij (x, Q ¯ ij are weighting functions wkl (x), the equivalent strengths Q kl determined by imposing the following condition:  (2),ij (2),ij ¯ ij ¯ ij ) dx = 0

wkl (Qkl ) = wkl (x)kl (x, Q kl

(i = 1, 2, . . . , n),

The previous expression represent an alternative form for Eq. (50); in particular, unlike Eq. (50), by means Eq. (53) it is possible to determine the equivalent strengths by adopting ¯ ij , a weighting function wkl (X) to obtain each unknown Q kl directly. In these conditions, when indexes i, j = 1, 2, . . . , n ¯ ij and k, l = 1, 2, . . . , m are fixed, each equivalent strengths Q kl assumes the following explicit form:   j2 w (X) E˜ gik (X)gj l (X) jXi kljXj (i, j = 1, 2, . . . , n) ¯ ij = Qkl Q   2 kl j w (X) (k, l = 1, 2, . . . , m) E˜ jXi kljXj

j =1

(47)

˜ denotes the ensemble average with respect to the where E[·]  ˜ = probability density function p(x), ˜ i.e. E[·] ˜ dx. In Rn (·)p(x) i order to find the n equivalent coefficients j it is necessary to select n weighting functions wr (x). When the index i = 1, 2, . . . , n is fixed, and for r = 1, 2, . . . , n, Eq. (47) becomes a linear system of n equations with the n equivalent coefficients ij as unkonwn. This linear system can be written in the following matrix form: Ci  i = v i

(i = 1, 2, . . . , n)

(48)

where Ci is a square matrix of order n and vi is a n-vector, whose elements Cji r and vri are given by the following expressions, respectively:     jwr (X) jwr (X) i i ˜ ˜ Cj r = E Xj , vr = E mi (X) . (49) jXi jXi ¯ ij equivalent strengths, by Moreover, with reference to the Q kl substituting Eqs. (10b), (34b) and (40) into the expression (46b), for i, j = 1, 2, . . . , n fixed we find:   m m 2 w (X)  j uv ij (2),ij ¯ ij ¯ )

wuv (Q )=E˜ (Qkl gik (X)gj l (X)−Q kl jXi jXj k=1 l=1

=0

(u, v = 1, 2, . . . , m).

(50)

The previous expression represent a system of m2 equations, whose solution furnishs the unknown m2 equivalent strengths ¯ ij . Q kl In order to simplify the previous procedure, for fixed i, j = 1, 2, . . . , n, in what follows the second component of the resid(2) ual global error ij , given by Eq. (43), will be split into the following form: (2) ¯ ij ) = ij (x, Q

m m  

(2),ij

kl

¯ ij )  = 0, (x, Q kl

(51)

k=1 l=1

(2)

Rn

(k, l = 1, 2, . . . , m)

(53)

By following a similar procedure to that one used to obtain Eq. (50), Eq. (53) becomes:   2 ij j wkl (X) (2),ij ¯ ij ˜ ¯

wkl (Qkl ) = E (Qkl gik (X)gj l (X) − Qkl ) =0 jXi jXj (k, l = 1, 2, . . . , m)

(54)

(55) The linear system of Eqs. (48) and expressions (55) represent the generalized form of the unknown vectors i and of the ¯ ij of the FPK equation (33) corresponding to the matrices Q equivalent linear system. The proposed PL method allows one to obtain different approximate solutions in terms of the unknown equivalent vectors ¯ ij , by selecting different weighti and equivalent matrices Q ing functions wr (X) and wkl (X). Moreover, the expressions of the equivalent coefficients and equivalent strengths depend on the particular choice of the probability density function used in order to calculate the unknown averages in Eqs. (48) and (55). In this way, the proposed PL method can be considered as a more general approach for generating approximate solutions. In particular in the next sections some remarkable cases of the proposed approach with reference to an SDOF non-linear oscillator will be illustrated. 6. The PL method for a SDOF non-linear oscillator Let us consider a SDOF non-linear oscillator subjected to both additive and multiplicative Gaussian white noise excitations, whose equation of motion is given as follows: ˙ = X¨ + H (X, X)

m 

˙ k (t). gk (X, X)W

(56)

k=1

where (2),ij

kl

¯ ) = (Qkl gik (x)gj l (x) − Q ¯ ) (x, Q kl kl ij

ij

(52)

In order to illustrate the proposed approach, Eq. (56) will be written as a system of two first-order differential

S. Lacquaniti, G. Ricciardi / International Journal of Non-Linear Mechanics 41 (2006) 1191 – 1205

equations:

L ¯ 22 (Q ) = b22

X˙ 1 = X2 ,

(57a)

X˙ 2 = −H (X1 , X2 ) +

m 

(57b)

gk (X1 , X2 )Wk (t).

k=1

Let p NL (x1 , x2 ) be the exact stationary probability density function of the response state vector process X(t)=[X1 (t) X2 (t)]T . The reduced FPK equation associated with Eqs. (57a,b) assumes the following expression: j NL j NL [a1 (x1 , x2 )p NL (x1 , x2 )] + [a (x1 , x2 )p NL (x1 , x2 )] jx1 jx2 2 −

1 j2 NL [b (x1 , x2 )p NL (x1 , x2 )] = 0, 2 jx22 22

(58)

where the first and the second-order derivative moments are related to the equation of motion by the following expressions:

1  Qkl gl (x1 , x2 ) 2 m

m

k=1 l=1

× NL b22 (x1 , x2 ) =

j gk (x1 , x2 ) jx2

m m  

(59a,b)

(60)

Qkl gk (x1 , x2 )gl (x1 , x2 )

k=1 l=1 NL = bNL = bNL = 0. By adopting the proposed procebeing b11 12 21 dure of PL, the reduced FPK equation (58) is replaced by the following one:

k=1 l=1

¯ 2+ ¯ 1 − X X˙ 2 = −kX

m 

W¯ k (t),

(64)

k=1 l=1

(65)

while v2 is a vector of two components given by the following expressions (for i = 1, 2)   jwi (X1 , X2 ) vi2 = E˜ m(X1 , X2 ) jX2   m  m  j 1 = E˜ −H (X1 , X2 ) + Qkl gl gk 2 jX2

(61)

(62b)

=

(63a,b)

(67)

x2 =X2

k¯ = 21

a1L (x1 , x2 ) = x2 ,

k=1 l=1

By using expression (59a,b) and by solving the linear system of two equations (65) with unknowns 21 and 22 , we find the following explicit expressions for the equivalent coefficients k¯ ¯ and :

(62a)

where k¯ = 21 and ¯ = 22 are the unknown equivalent coefficients, while W¯ k (t) are m equivalent Gaussian white noises with ¯ kl (t1 − t2 ). In cross-correlation functions E[W¯ k (t1 )W¯ l (t2 )] = Q Eq. (61) the following expressions:

jwi (X1 , X2 ) jX2



where m(X1 , X2 ) = m2 (X1 , X2 ) = a2NL (x1 , x2 )| x1 =X1 .

k=1

a2L (x1 , x2 ) = a2L (x1 , x2 , 2 ) = − 21 x1 − 22 x2 ,

¯ kl Q

C2 2 = v2 ,

×

corresponding to an equivalent linear oscillator subjected to additive excitations only, whose equations of motion have the following expressions: X˙ 1 = X2 ,

m  m 

are the first and the second-order derivative moments relative L = bL = bL = 0. to the equivalent linear oscillator, being b11 12 21 L NL In this case, since a1 (x1 , x2 ) = a1 (x1 , x2 ) = x2 , we find immediately that 11 = 0 and 12 = 1, i.e. 1 = [0 1]T . Therefore, the proposed procedure of probabilistic linearization allows to determinate the remaining equivalent coefficients 21 and 22 , ¯ T and the m2 uncollected in the vector 2 = [ 21 22 ]T = [k¯ ] ¯ kl collected in the square matrix known equivalent strengths Q ¯ 22 of order m. Q By selecting two weighting functions w1 (X1 , X2 ) and w2 (X1 , X2 ), the equivalent coefficients 21 and 22 will be determined by solving the system of equations (48), which in this case assumes the following simplified form:

j L j L [a1 (x1 , x2 )p L (x1 , x2 )] + [a (x1 , x2 )p L (x1 , x2 )] jx1 jx2 2 1 L j2 L − b22 [p (x1 , x2 )] = 0 2 jx22

¯ 22 Q kl =

where C2 is a square matrix of order 2, whose elements are given by:   jwj (X1 , X2 ) 2 ˜ (ij = 1, 2) (66) Cij = E Xi jX2

a1NL (x1 , x2 ) = x2 a2NL (x1 , x2 ) = − H (x1 , x2 ) +

m  m 

1199

  ˜ ˜ 2 w  ]E[m(X ˜ ˜ 2 w  ]E[m(X E[X 1 , X2 )w1 ] − E[X 1 , X2 )w2 ] 2 1 , ˜ 1 w  ]E[X ˜ 2 w  ] − E[X ˜ 2 w  ]E[X ˜ 1 w ] E[X 1 2 1 2 (68a)

¯ = 22 =

  ˜ ˜ 1 w  ]E[m(X ˜ ˜ 1 w  ]E[m(X E[X 1 , X2 )w2 ] − E[X 1 , X2 )w1 ] 1 2 ˜ 1 w  ]E[X ˜ 2 w  ] − E[X ˜ 2 w  ]E[X ˜ 1 w ] E[X 1 2 1 2 (68b)

where wi (X1 , X2 ) = jwi (X1 , X2 )/jX2 for i = 1, 2.

1200

S. Lacquaniti, G. Ricciardi / International Journal of Non-Linear Mechanics 41 (2006) 1191 – 1205

¯ kl , let In order to find the m2 unknown equivalent strengths Q wkl (X1 , X2 ) (for k, l =1, 2, . . . , m) be m2 weighting functions, it is enough to apply expression (55), directly; in this case it assumes the following simplified form: ¯ kl = kl Qkl Q   j2 wkl (X1 ,X2 ) ˜ E gk (X1 , X2 )gl (X1 , X2 ) jX22   kl = j2 wkl (X1 ,X2 ) E˜ jX 2

(69a,b)

2

It is important to note that, on the basis of the expressions ¯ kl depend on: (69a,b), the equivalent coefficients 21 , 22 and Q (a) the probability density function used in order to calculate the averages appearing in the expressions (68a,b) and (69a,b); (b) on the two weighting functions w1 (X1 , X2 ) and w2 (X1 , X2 ); (c) on the m2 functions wkl (X1 , X2 ) chosen. In this way, the proposed PL can be considered as a general procedure of SL made on the FPK equation of the original non-linear oscillator.

the Itô differential rule. In this way, in the case of multiplicativetype excitations, the results obtained with this proposed approach, named Gaussian Probabilistic Linearization (GPL) coincide with those obtained by means of the Gaussian Probabilisitc Linearization (GSL) (on the Itô differential rule) and GC approaches. 7.2. Generalized Gaussian Probabilistic Linearization (GGPL) The propose GPL method can be generalized. In fact, in the case in which we assume as weighting functions two functions w1 (X1 , X2 ) and w2 (X1 , X2 ) different to those given by Eq. (70), by adopting the Gaussian probability density function ¯ 22 ), the expressions (68a,b) and (69a,b) of the pG (x1 , x2 ; 2G , Q G equivalent coefficients become: k¯ = k¯G =

EG [X2 w2 ]EG [m(X1 , X2 )w1 ] − EG [X2 w1 ]EG [m(X1 , X2 )w2 ] EG [X1 w1 ]EG [X2 w2 ] − EG [X2 w1 ]EG [X1 w2 ]

7. Some particular cases of the proposed PL method

(72a)

In this section will be illustrated some special cases of the PL method based on different choices of the weighting functions used in the expressions of the equivalent coefficients and on the probability density function used in order to calculate the averages in the expressions (68a,b) and (69a,b).

¯ = ¯ G =

EG [X1 w1 ]EG [m(X1 , X2 )w2 ]−EG [X1 w2 ]EG [m(X1 , X2 )w1 ] EG [X1 w1 ]EG [X2 w2 ] − EG [X2 w1 ]EG [X1 w2 ]

By adopting the following weighting functions: w2 (X1 , X2 ) = w22 (X1 , X2 ) = X22

,

(72b) ¯ kl = Q ¯ kl,G = ¯ kl,G Qkl Q   j2 w (X ,X ) EG gk (X1 , X2 )gl (X1 , X2 ) kljX21 2 2   Qkl . = 2 j wkl (X1 ,X2 ) EG jX 2

7.1. Gaussian Probabilistic Linearization (GPL)

w1 (X1 , X2 ) = X1 X2 ,

,

(72c)

2

(70) and by evaluating the averages by means of the joint Gaus¯ 22 ), which sian probability density function pG (x1 , x2 ; 2G , Q G 2 ¯ 22 of depend on the unknown equivalent coefficients G and Q G the equivalent linear oscillator, we obtain the following Gaussian expressions of the equivalent coefficients: EG [m(X1 , X2 )X1 ] k¯ = kG = , EG [X12 ] EG [m(X1 , X2 )X2 ] ¯ = G = EG [X22 ]

(71a,b)

¯ kl = Qkl,G = kl,G Qkl Q = EG [gk (X1 , X2 )gl (X1 , X2 )]Qkl (71c)  ¯ 22 ) dx1 dx2 . Eqs. (71a–c) where EG [·]= R2(·)pG (x1 , x2 ; 2G , Q G have been obtained by taking into account that the response processes X1 (t) and X2 (t) are independent and with zero mean. It is worth noting that the expressions of equivalent coefficients, here achieved, coincide with the expressions obtained by Falsone [16], who apply the approach of SL on the coefficients of

In general, a more accurate solution in terms of the equivalent coefficients can be obtained by a suitable choice of the weighting functions. In this way, on the basis of an appropriate form of the weighting functions and by using the Gaussian probability density function, the proposed approach named Generalized Gaussian Probabilistic Linearization (GGPL) allows to obtain an improved expression of the equivalent coefficients and consequently, a better approximation of the FPK equation of the original non-linear oscillator. 8. Numerical applications In this section, in order to show the results obtained by means the GGPL method, two non-linear oscillators under both additive and multiplicative excitations are illustrated. In both cases, on the basis of the particular choice of the weighting functions, the proposed method can lead to approximate solutions in terms of variance of the response process that appear more accurate than the GPL (or GSL). Moreover, from these examples and for specific set of parameters, it will be observed that certain combinations of the weighting functions can lead to exact variance solutions. In this way, by following the terminology of

S. Lacquaniti, G. Ricciardi / International Journal of Non-Linear Mechanics 41 (2006) 1191 – 1205

Kozin [21], the proposed method can be interpreted as a “true” linearization of the non-linear system.

1 (2,2)

0.6 (2,4) exact + 10 % error

˙ 2 (t) + W3 (t), = XW 1 (t) + XW

(73)

where and  are non-negative constants, W1 (t), W2 (t) and W3 (t) are uncorrelated Gaussian white noises with zero mean and strengths Q11 , Q22 and Q33 , respectively. In particular, this example is performed by assuming that the multiplicative excitations have the same strengths Q11 =Q22 =Q. Let X1 =X ˙ the drift and diffusion coefficients are given by and X2 = X, the following expressions: mNL 1 (X1 , X2 ) = X2 2 3 1 mNL 2 (X1 , X2 ) = − X2 − (X1 X2 + X2 ) − X1 + 2 QX 2

(74a,b)

- 10 % error

0.2 0

4

6

the following expressions:     jw jw EG X12 X2 jX12 + EG X23 jX12 1   − Q, ¯ = +  jw1 2 EG X2 jX2 

j2 w X12 jX211

EG



EG

10

 Q,

¯ 22 = Q

EG

j2 w X22 jX222



EG

(79a)



2

j2 w22 jX22

(78)



2

j2 w11 jX22



 Q,

¯ 33 = Q33 . Q (75)

8

Fig. 1. Pairs of values of Q and  for which the GGPL method with = 2 and for different values of  leads to a percentage error within ±10% for the variance of the response process.

¯ 11 = Q

j [x2 p ex (x1 , x2 )] jx1   j 1 2 3 − x2 − (x1 x2 + x2 ) − x1 + Qx 2 + 2 jx2  1 j2 × pex (x1 , x2 ) − [((x12 + x22 )Q 2 jx22

2



(74c)

In the stationary state the reduced FPK equation is:

+ Q33 )p (x1 , x2 )] = 0.

(2,3)

0.4

X¨ + X˙ + (X 2 X˙ + X˙ 3 ) + X

ex

Q

As a first example, let us consider the non-linear oscillator subjected to both additive and multiplicative excitations whose dynamic behavior is governed by the following differential equation:

+ X22 )Q + Q33

(2,2.5)

0.8

8.1. A van der Pol-type oscillator

222 (X1 , X2 ) = (X12

1201

(79b)

(79c)

For this particular case, the exact probability density function p ex (x1 , x2 ) of the response process is known [2] By applying the proposed GGPL method, Eq. (75) is replaced by the following equation:

On the basis of the particular choice of the weighting functions w1 , w11 and w22 , and since the response is assumed Gaussian, the proposed method allows one to obtain the variance of ¯ the response process in terms of the equivalent coefficients , ¯ 11 , Q ¯ 22 , Q ¯ 33 by means of the following expression: Q

j j ¯ 2 − x1 )p L (x1 , x2 )] [x2 p L (x1 , x2 )] + [(−x jx1 jx2

2X1 = 2X2 =

2 1 ¯ ¯ 22 + Q ¯ 33 ) j [p L (x1 , x2 )] = 0, − (Q 11 + Q 2 jx22

(76)

where p L (x1 , x2 ) is the probability density function of the response process of the following equivalent linear oscillator subjected to additive excitations only: X¨ + ¯ X˙ + X = W¯ 1 (t) + W¯ 2 (t) + W¯ 3 (t).

(77)

¯ Q ¯ 22 , Q ¯ 33 can be deter¯ 11 , Q The equivalent coefficients , mined by using the expressions (72b,c). By taking into account that the response process is with zero mean and the displacement and velocity are two independent processes, they assume

¯ 22 + Q ¯ 33 ¯ 11 + Q Q . ¯ 2

(80)

In Figs. 1 and 2, for assigned = Q33 = 1.0, the continuous lines represent the values of Q and  for which the proposed method leads to the exact variance of the response, for different weighting functions. The following notation has been adopted ( , ), that represent the order of the power of the following weighting functions used: w1 (X1 , X2 ) = |X2 | w11 (X1 , X2 ) = w22 (X1 , X2 ) = |X2 |

(81)

For =  = 2 the solutions obtained with the proposed method, indicated with the notation (2, 2), are coincident with the GPL

1202

S. Lacquaniti, G. Ricciardi / International Journal of Non-Linear Mechanics 41 (2006) 1191 – 1205 1

2.2 =1 =2 =3 =4 =5 GPL

(3,2)

(4,2)

1.8

0.8

2 X2 /X,ex

Q

(2.5,2)

0.6

1.4

1

(2,2)

0.4

exact +10 % error -10 % error

0.6

0.2 0

0.2

0.4

0.6

0.8

0.2

1

0



Fig. 2. Pairs of values of Q and  for which the GGPL method with  = 2 and for different values of leads to a percentage error within ±10% for the variance of the response process.

1.6

1.5

2 

2.5

3

1.2

1

0.8

3.5

4

 = 1.0  = 0.8  = 0.6  = 0.4  = 0.2  =0 GPL

1.1

2 X2 /X,ex

2 X2 /X,ex

1.2

1

Fig. 4. Normalized variance of the response process for a van der Pol oscillator with parameters: = 3.65,  = 1.4, Q = 12 and Q33 = 20. Comparison of the solution obtained by GPL and GGPL method versus  for different values of .

=1 =2 =3 =4 =5 GPL

1.4

0.5

1

0.9

0.6 0

1

2

3

4

5

6

7



0.8 0

0.2

0.4

0.6

0.8

1

1.2

1.4



Fig. 3. Normalized variance of the response process for a van der Pol oscillator with parameters: = 1.3,  = 4.0, Q = 2 and Q33 = 5. Comparison of the solution obtained by GPL and GGPL method versus  for different values of .

(or GSL) method. Moreover, in the same figures, the shading regions define the values of Q and  for which the proposed method leads to an error in term of the response variance contained within ±10%. By adopting the same weighting functions given in Eqs. (81), in Figs. 3 and 4 the variances of the response normalized to exact value 2X /2X,ex are depicted, for different values of and . In particular, Fig. 3 is relative to a van der Pol oscillator with parameters = 1.3,  = 4.0, Q = 2 and Q33 = 5, for which the GPL method underestimates the exact variance; in this case for fixed 1  5, the proposed method leads to the exact value by adopting a order of power   . As an example, for = 1, the exact variance is obtained for   1; while for = 4, the exact variance is obtained for   5. Fig. 4 is relative to the case in which = 3.65,  = 1.4, Q = 12 and Q33 = 20, for which the GPL overestimates the exact variance; in this

Fig. 5. Normalized variance of the response process for a van der Pol oscillator with parameters: = 1.3,  = 4.0, Q = 2 and Q33 = 5. Comparison of the solutions obtained by GPL and GGPL method versus for different values of .

second case, for fixed 1 5 the proposed method leads to the exact value by adopting a order of power   . In fact, for = 1, the exact variance is obtained for   43 , while for = 4, the exact variance is obtained for   25 . Nevertheless, in both cases, the exact variance can be obtained by selecting different weighting functions corresponding to more pairs of and . In this way, the GGPL method is equivalent to different “true” linearizations of the original non-linear system. In Figs. 5 and 6 the approximate solutions in term of variance normalized to the exact value obtained with the GGPL method for 0  1.4 and 0  1 are illustrated. In particular, in Fig. 5 the previous case in which the GPL underestimates the exact variance is reported; in this case, for fixed 0  1, the exact variance is found for . Moreover, from Fig. 6, relative to the case in which the GPL overestimates the exact

S. Lacquaniti, G. Ricciardi / International Journal of Non-Linear Mechanics 41 (2006) 1191 – 1205

The proposed GGPL method is equivalent to replacing the FPK equation of the system (82) with an equivalent one relative to the following linear oscillator:

1.8  = 1.0  = 0.8  = 0.6  = 0.4  = 0.2 GPL

1.6

2 X2 /X,ex

1.4

1203

¯ = W¯ 1 (t) + W¯ 2 (t), X¨ + ¯ X˙ + kX

(84)

0.8

where the unknown parameters are the equivalent coefficients ¯ 11 and Q ¯ 22 . Since the ¯ and k¯ and the equivalent strengths Q velocity response of the non-linear system is Gaussian, the GGPL method, will be applied in order to find an improved solution for the variance of the displacement. The following weighting functions have been used:

0.6

w1 (X1 , X2 ) = |X1 | sign(X1 )X2 ,

1.2 1

0

0.2

0.4

0.6

0.8

1

1.2

1.4



Fig. 6. Normalized variance of the response process for a van der Pol oscillator with parameters: = 3.65,  = 1.4, Q = 12 and Q33 = 20. Comparison of the solutions obtained by GPL and GGPL method versus for different values of .

variance, two different behaviors can be evidenced: for fixed 0.2  0.4 the proposed method leads to the exact solution for a value , while for 0.6  1.0 the exact solution is obtained for a value > . As evidenced in the cases examined, the exact solution can be obtained by means of different combinations of the weighting functions. According to Elishakoff [19,20] who, obtained the exact variance of the displacement for a Duffing oscillator for a specific set of parameters by means of the energy based stochastic linearization, it is possible to affirm that there is not a unique “true” linearization valid for all cases, but different “true” linearizations exist which lead to the exact variance, characterized by different equivalent coefficients. 8.2. A non-linear Duffing-type oscillator under additive and multiplicative excitations As a second example, let us consider the non-linear oscillator whose equation of motion is given by: ˙ 2 X˙ + εX 3 X¨ + X˙ + (X + X) ˙ 1 (t) + W2 (t), = (aX + bX)W

(82)

where the excitations W1 (t) and W2 (t) are uncorrelated Gaussian white noise with zero mean and strengths Q11 and Q22 , respectively. This example is performed by assuming that: b = a,

=

  Q11 2 Q11 2 a + b . Q22 2

(83)

In these conditions, the system belongs to the class of the generalized stationary potential [2], and the exact probability density function of the response is known. The following parameters have been selected: a =  = 0.1,  = 1.0, and ε = 10.

w2 (X1 , X2 ) = X1 X2 , w11 (X1 , X2 ) = w22 (X1 , X2 ) = X22 ,

(85)

where 0  1. By using the expressions (72a–c), the following equivalent parameters are obtained:  k¯ = 2 EG [X22 ] + 2ε 1 + EG [X12 ], 2 ¯ =  + EG [X12 ] + 3 2 EG [X22 ] − 21 b2 Q11 ,

(86a,b)

¯ 11 = (a 2 EG [X12 ] + b2 EG [X22 ])Q11 , Q

(86c,d)

¯ 22 = Q22 . Q

Since for the equivalent linear system (84) the variance of the displacement and velocity are given by EG [X12 ] =

¯ 22 ¯ 11 + Q Q , 2¯ k¯

EG [X22 ] =

¯ 22 ¯ 11 + Q Q 2¯

(87a,b)

by substituting the equivalent parameters (86b–d) in Eq. (87b), 2 we obtain the exact variance of the velocity 2,ex X2 = EG [X2 ] = 2.8571, that coincide with the solution of the proposed GPL (or GSL) method with = 1. From the remaining Eq. (87a) and by using the expressions (86a,c) will be possible to determine the approximate solution in terms of the variance of the displacement. In Figs. 7(a)–(d), the percentage errors e for the variance of the displacement are depicted for different levels of the excitations. From these figures, it appears that, for each level of the additive and multiplicative excitation, the GGPL method leads to small errors (< 7%) for values of < 0.6. On the contrary, in Figs. 7(a) and 7(b) it is possible to note that for > 0.8 we obtain higher errors (> 10%) for low level of the multiplicative excitation only, while for high level these errors are contained within 8%. Moreover, in Figs. 7(c) and 7(d) is shown that the errors are greater than 10% for each level of the multiplicative excitations when > 0.8. It appears that it is possible to obtain the exact variance by choosing an appropriate value of the coefficient . In particular, in Fig. 7(a) is shown that for low level of the multiplicative excitation the exact solution is obtained for low values of ( 0.2), while for high level of the multiplicative excitation it is need to choice a greater coefficient ( 0.7). In Fig. 7(b) we note that for Q11 = 10.0, the exact solution is obtained for 0.45. Moreover, for high levels of the additive excitations, as

1204

S. Lacquaniti, G. Ricciardi / International Journal of Non-Linear Mechanics 41 (2006) 1191 – 1205

14

14 12

10

8

e (%)

e (%)

10

Q22 = 0.1 6

8

Q22 = 0.1

6

4

4

2

2

0

0 0

0.2

0.4

(c)

0.6

0.8

1

0

0.2

0.4

(d)



0.6

0.8

1

0.6

0.8

1

14

14 Q11 = 0.1 Q11 = 1.0 Q11 = 5.0 Q11 = 10.0

10 8

Q11 = 0.1 Q11 = 1.0 Q11 = 5.0 Q11 = 10.0

12 10 e (%)

12

e (%)

Q11 = 0.1 Q11 = 1.0 Q11 = 5.0 Q11 = 10.0

12

Q11 = 0.1 Q11 = 1.0 Q11 = 5.0 Q11 = 10.0

Q22 = 5.0

6

8

Q22 = 10.0

6

4

4

2

2 0

0 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4



Fig. 7. Percentage error for the variance of the displacement of the Duffing-type oscillator versus for different values of the intensity Q11 of the multiplicative excitation and for assigned intensity Q22 of the additive excitation: (a) Q22 = 0.1; (b) Q22 = 1.0; (c) Q22 = 5.0; (d) Q22 = 10.0.

shown in Figs. 7(c) and 7(d), the exact variance is obtained for low values of ( 0.2–0.3) for each level of the multiplicative excitations. Also for this example, the application of the GGPL method affrims that there does not exist a unique linearization method valid in general for all non-linear systems, but it is possible to make different linearizations for the specific cases, by choosing in a suitable way the weighting functions used in order to obtain the equivalent coefficients. 9. Conclusions A new version of the stochastic linearization (SL) method has been proposed. This method, called “Probabilistic Linearization” (PL), (following the terminology of Polidori and Beck [17]), based on the linearization of the FPK equation of the original system, represents a general approach to produce approximate solution in the study of non-linear systems subjected to additive and multiplicative random excitations. The application of the PL method allows one to obtain the unknown parameters of an equivalent linear system, whose FPK equation is close to the original one. The explicit expressions

of the equivalent coefficients have been obtained by means of the weighted residual method. The solutions achieved with the proposed method depend on the weighting functions used and on the probability density function adopted in order to calculate the expressions of the equivalent coefficients. By assuming the Gaussian probability density function and by choosing in a suitable way the weighting functions, two different versions of the PL method have been illustrated: the GPL and the GGPL method. The numerical analyses performed on a van der Pol-type and a Duffing-type oscillators under additive and multiplicative excitations allowed to obtain different solutions in terms of variance, by varying the weighting functions. In particular, for each specific case, several combinations of the weighting functions, chosen in opportune way, furnish the same exact value of the variance, even though they lead to different equivalent coefficients. References [1] T.K. Caughey, F. Ma, The exact steady-state solution of a class of nonlinear stochastic systems, Int. J. Non-Linear Mech. 17 (3) (1982) 137–142.

S. Lacquaniti, G. Ricciardi / International Journal of Non-Linear Mechanics 41 (2006) 1191 – 1205 [2] Y. Lin, G.C. Cai, Probabilistic Structural Dynamics: Advanced Theory and Applications, McGraw-Hill, New York, 1995. [3] Y.K. Lin, G.C. Cai, Exact stationary-response solution for second-order nonlinear systems under parametric and external white-noise excitations: Part II, ASME J. Appl. Mech. 55 (3) (1988) 702–705. [4] R.C. Booton, The analysis of nonlinear control systems with random inputs, in: Proceedings of the Symposium on Nonlinear Circuit Analysis, vol. 2, Polytechnic Institute of Brooklyn, 1953. [5] I.E. Kazakov, Approximate probabilistic analysis of the accuracy of operation of essentially nonlinear systems, Autom. Remote Control 17 (1956) 423–450. [6] T.K. Caughey, Equivalent linearization techniques, J. Acoust. Soc. Am. 35 (1963) 1706–1711. [7] J.B. Roberts, P.D. Spanos, Random Vibration and Statistical Linearization, Wiley, New York, 1990. [8] L. Socha, Linearization in analysis of nonlinear stochastic systems: recent results—part I: theory, Trans. ASME 58 (2005) 178–205. [9] L. Socha, T.T. Soong, Linearization in analysis of nonlinear stochastic systems, Appl. Mech. Rev. 44 (1991) 399–422. [10] H. Risken, The Fokker–Planck equation, Methods of Solution and Application, Springer, Berlin, 1989. [11] C.W. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences, Springer, Berlin, 1983. [12] R.I. Chang, G.E. Young, Methods and Gaussian criterion for statistical linearization of stochastic parametrically and externally excited nonlinear systems, J. Appl. Meth. 56 (1989) 179–186. [13] A. Bruckner, Y.K. Lin, Generalization of the equivalent linearization method for non-linear random vibration problems, Int. J. Non-Linear Mech. 22 (3) (1987) 227–235. [14] J. Kottalam, K. Lindberg, B. West, Statistical replacement for systems with delta-correlated fluctuations, J. Statist. Phys. 42 (1986) 979–1008. [15] W.F. Wu, Comparison of Gaussian closure technique and equivalent linearization method, Probab. Eng. Mech. 2 (1) (1987) 2–8. [16] G. Falsone, Stochastic linearization of MDOF systems under parametric excitations, Int. J. Non-Linear Mech. 27 (6) (1992) 1025–1037. [17] D.C. Polidori, J.L. Beck, Approximate solutions for non-linear random vibration problems, Probab. Eng. Mech. 11 (1996) 179–185. [18] G.Q. Cai, Y.K. Lin, I. Elishakoff, A new approximate solution technique for randomly excited non-linear oscillators, part II, Int. J. Non Linear Mech. 27 (6) (1992) 969–979.

1205

[19] I. Elishakoff, Method of stochastic linearization revised and improved, in: P.D. Spanos, C.A. Brebbia (Eds.), Fifth International Computational Stochastic Methods, Southampton, 1991, pp. 101–111. [20] I. Elishakoff, Stochastic linearization technique: a new interpretation and a selective review, Shock Vib. Dig. 32 (2000) 179–188. [21] F. Kozin, The method of statistical linearization for non-linear stochastic vibrations, in: F. Ziegler, G.I. Schueller (Eds.), IUTAM Symposium, Nonlinear Stochastic Dynamic Systems, Springer, Berlin, 1988, pp. 45–66. [22] E. Wong, M. Zakai, On the relation between ordinary and stochastic differential equations, Int. J. Eng. Sci. 3 (1965) 213–229. [23] T.K. Caughey, J.K. Dienes, Analysis of a nonlinear first-order system with a white noise input, J. Appl. Phys. 23 (11) (1961) 2476–2479. [24] T.K. Caughey, Nonlinear theory of random vibrations, Adv. Appl. Mech. 11 (1971) 209–253. [25] T.K. Caughey, On the response of a class of nonlinear oscillator to stochastic excitation, in: Proceedings of Colloquium International du Centre Nat. de la Rech. Scientifique, vol. 148, Marseilles, 1964, pp. 393–402. [26] R. Graham, H. Haken, Generalized thermo-dynamics potential for Markov systems in detailed balance and far from thermal equilibrium, Z. Phys. 243 (1971) 289–302. [27] M.F. Dimentberg, An exact solution to a certain nonlinear random vibration problem, Int. J. Non-Linear Mech. 17 (4) (1982) 231–236. [28] Y. Yong, Y.K. Lin, Exact stationary-response solution for second-order nonlinear systems under parametric and external white-noise excitations, ASME J. Appl. Mech. 54 (2) (1987) 414–418. [29] J.J. Beaman, J.K. Hendrick, Improved statistical linearization for analysis and control of nonlinear stochastic systems: part I: an extended statistical linearization technique, Trans. ASME 103 (1981) 14–21. [30] H.J. Pradlwarter, Non-Gaussian linearization, An efficient tool to analyze nonlinear MDOF-systems, J. Nucl. Eng. Des. 128 (1991) 175–192. [31] R.J. Chang, Non-Gaussian linearization method for stochastic parametrically and externally excited nonlinear systems, Trans. ASME 114 (1992) 20–26. [32] G. Ricciardi, A non-Gaussian stochastic linearization method, Prob. Eng. Mech. 22 (1) (2007) 1–11. [33] S.H. Crandall, Heuristic and equivalent linearization techniques for random vibration of nonlinear oscillators, in: Proceedings of Eighth International Conference on Nonlinear Oscillators, vol. 1, Academica, Prague, 1978, pp. 211–226.