- Email: [email protected]
Monte Carlo simulation in the stochastic analysis of non-linear systems under external stationary Poisson white noise input
Available online at www.sciencedirect.com
International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
Monte Carlo simulation in the stochastic analysis of non-linear systems under external stationary Poisson white noise input G. Muscolino ∗ , G. Ricciardi, P. Cacciola Dipartimento di Costruzioni e Tecnologie Avanzate, University of Messina, Salita Sperone 31, S. Agata, I-98166 Messina, Italy
Abstract A method for the evaluation of the probability density function (p.d.f.) of the response process of non-linear systems under external stationary Poisson white noise excitation is presented. The method takes advantage of the great accuracy of the Monte Carlo simulation (MCS) in evaluating the 1rst two moments of the response process by considering just few samples. The quasi-moment neglect closure is used to close the in1nite hierarchy of the moment di4erential equations of the response process. Moreover, in order to determine the higher order statistical moments of the response, the second-order probabilistic information given by MCS in conjunction with the quasi-moment neglect closure leads to a set of linear di4erential equations. The quasi-moments up to a given order are used as partial probabilistic information on the response process in order to 1nd the p.d.f. by means of the C-type Gram–Charlier series expansion. ? 2002 Elsevier Science Ltd. All rights reserved. Keywords: Non-linear stochastic dynamics; Stationary Poisson process; Monte Carlo simulation; Non-Gaussian probability density function; Quasi-moment neglect closure
1. Introduction A broad class of random excitations on structures can be modeled as Gaussian white noise or 1ltered Gaussian white noise processes. Then both analytical and numerical methods have been recently developed to determine the response of structures driven by these types of excitations. However, when excitations are characterized by impulsive loading, these representations may not be satisfactory. Moving loads traveling on a bridge [1,2], earthquake shake on structures [3] and wave action on ships [4] are examples of this type of excitations. ∗
Corresponding author. Tel.: +390-90-676-5620; fax: +39090-395022. E-mail address: [email protected] (G. Muscolino).
In these cases the Poisson white noise and the 1ltered Poisson white noise processes are more realistic models [5]. Methods for the evaluation of the response of linear and non-linear systems under non-Gaussian delta-correlated processes have been recently proposed by several authors. Some of these are based on the theory of Markov processes and determine the probability density function (p.d.f.) of the response process as solution of the generalized Kolmogorov equation [6]. Numerical methods have been developed to solve this di4erential equation such as the Petrov–Galerkin method [7] and the path integral method [8,9]. Recent works on this subject utilize the characteristic function to describe the response process. These methods require the solution of the partial di4erential equation for the characteristic function obtained by
0020-7462/03/$ - see front matter ? 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 4 6 2 ( 0 2 ) 0 0 0 7 2 - 0
1270
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
taking the Fourier transform of the generalized Kolmogorov equation [10] or starting by the generalized Itˆo’s di4erential rule for Poisson white noise input [11]. Several approaches based on the moment equation method have been proposed [12,13]. Unfortunately, these methods su4er from the problem associated with the approximation techniques used to close the in1nite hierarchy of equations. The Monte Carlo simulation (MCS) constitutes the universal alternative method of solution when the system under study is large and strongly non-linear. This method requires the generation of samples of the input process, the deterministic analysis for each of them and the evaluation of the resultant response statistics, which are used to develop estimates of various characteristics of the response process. This approach can provide accurate solutions for a large class of stochastic problems [14]. Unfortunately, the disadvantages of MCS method are that it is generally time-consuming, when the statistics of higher order or the p.d.f. of the response process are of interest. Moreover, the MCS is rather ineIcient for estimating small probabilities such as those in the tail of the response distribution, which are particularly important for the reliability of structural systems. However, the great accuracy of MCS in estimating lower order statistics of the response has been recently exploited in stochastic dynamics of non-linear systems [15,16]. Following this strategy, in this paper an alternative approach to evaluate the stationary p.d.f. of non-linear systems subjected to external stationary Poisson white noise input is presented. The proposed approach takes advantage of: (a) the great accuracy of the MCS in evaluating the 1rst two moments of the response process with a small number of input samples; (b) the non-Gaussian closure (NGC) technique in conjunction with the moment equation approach, which leads to a set of di4erential equations in which the non-linearity is due to the presence of mean and variance of the response process [17]. As a consequence of the evaluation of the mean and variance by MCS, the di4erential equations of moments of higher order than two become a closed set of linear di4erential equations. In the stationary case, the set of linear di4erential equations becomes an algebraic one. Once a good approximation of the statistical moments up to a given order
is achieved, the p.d.f. of the response process can be easily determined. In this paper the so-called C-type Gram–Charlier (CGC) series expansion [18] is used, whose coeIcients are determined as linear functions of quasi-moments of the response process [19,20]. 2. Preliminary considerations 2.1. One-dimensional systems For clarity’s sake and in order to evidence the advantages of the proposed formulation, a non-linear one-dimensional system is 1rst studied. Indeed, as it will be shown later, the multi-dimensional case is formally similar to the one-dimensional case if the Kronecker algebra is used. Let us consider the di4erential equation describing a one-dimensional system subjected to a purely external Poisson white noise process, given by ˙ = a[Z(t); t] + g(t)(t); Z(t)
(1)
where Z(t) is the response process, the dot means time di4erentiation, a[Z(t); t] is a deterministic non-linear function of the response process Z(t) and g(t) is a deterministic time-dependent function. In Eq. (1), (t) is a scalar Poisson white noise process [21], that is a sequence of impulses with random amplitude and arriving at random times, given by (t) =
N (t)
Yk (t − tk );
(2)
k=1
the amplitudes {Yk } being a family of random variables, mutually independent and independent of the time instants tk , with prescribed distribution pY (y),
(t) the Dirac’s delta function, tk the random time arrivals. In Eq. (2), N (t) is a counting Poisson process with parameter . The probabilistic characterization of the Poisson white noise process (t) can be given in terms of its correlation functions [22] Cr [(t1 ); (t2 ); : : : ; (tr )] = E[Y r ] (t1 − t2 ) (t1 − t3 ) · · · (t1 − tr );
(3)
where E[ • ] means stochastic average. It follows that (t) is a non-normal delta-correlated process, with intensity coeIcients related to the probabilistic
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
1271
characterization of the random variables Yk and of the counting Poisson process N (t). The Poisson white noise process (t) is not mean square Riemann integrable and, consequently, Eq. (1) has no traditional mathematical meaning. The latter can be considered to be formally equivalent to the following generalized Itˆo-type stochastic di4erential equation:
Since (dZ)r is an in1nitesimal of order dt, it follows that the series expansion of the increment N(Z r ) has to be written in the following form [1,23]: r 1 k r d (Z ) = rZ r−1 dZ N(Z r ) = k!
dZ(t) = a[Z(t); t] dt + g(t) d(t);
where Lr; k = r!=[k!(r − k)!]. Eq. (11) is obtained observing that neglecting terms of order higher than dt the following relationships hold:
(4)
where d(t) is an increment of the compound Poisson process (t), de1ned by (t) =
N (t)
Yk U (t − tk );
(5)
k=1
U (t) being the unit step function. The compound Poisson process is characterized by the following correlation functions [1]: Cr [(t1 ); (t2 ); : : : ; (tr )] = E[Y r ] min(t1 ; t2 ; : : : ; tr ):
(6)
If we proceed formally, that is in terms of correlation functions, it can be shown that the Poisson white noise process can be considered the derivative of the compound Poisson process. In fact, the correlation functions of an increment of the compound Poisson process are given by the following relationships: Cr [d(t1 ); d(t2 ); : : : ; d(tr )] = E[Y r ] (t1 − t2 ) (t1 − t3 ) · · ·
(t1 − tr ) dt1 dt2 : : : dtr :
(7)
Accounting for the relationships between correlation functions and moments, putting t1 = t2 = · · · = tr = t and neglecting in1nitesimals of order higher than dt, we obtain mr [d(t)] ≡ kr [d(t)] = E[Y r ] dt;
(8)
where mr [ • ] means rth statistical moment and kr [ • ] means the rth cumulant. Moreover, by taking into account the non-anticipating property, the following relationships hold: E[(Z(t))k (d(t))r ] = E[Y r ]mk [Z(t)] dt:
(9)
Neglecting in1nitesimals of order higher than dt, we have (dZ)r = gr (t)(d)r ;
r ¿ 2:
(10)
k=1
+
r
Lr; k Z r−k gk (t)(d)k ;
(11)
k=2
d(Z r ) = rZ r−1 dZ; d k (Z r ) = k!Lr; k Z r−k (dZ)k = k!Lr; k Z r−k gk (t)(d)k ; k = 2; 3; : : : ; r: (12) Using the stochastic average of both sides of Eq. (11), taking into account Eqs. (8) – (10) and dividing the result by dt, we obtain the di4erential equation governing the evolution of the rth moment of the response m˙ r [Z(t)] = rE[Z r−1 (t)a[Z(t); t]] r + Lr; k mr−k [Z(t)]gk (t) E[Y k ]:
(13)
k=1
Notice that the stationary rth order moment of the response can be evaluated by solving the following algebraic equation: rE[Z r−1 a[Z]] r + Lr; k mr−k [Z]gk E[Y k ] = 0;
(14)
k=1
obtained from Eq. (13) observing that this equation is not time dependent. Moreover, since a zero mean Gaussian white process can be considered as delta-correlated up to the second order, for this input the sum in Eqs. (13) and (14) need to be performed only for k = 2. In this case q2 = E[Y 2 ] is usually de1ned as the strength of the white process. 2.2. Multi-dimensional systems Let us consider a N -degree-of-freedom structural system, whose nodal displacements are the elements of vector u(t). By introducing the vector of state variables Z (t) = [uT (t) u˙T (t)]T , the equation of motion of the non-linear N -DOF structural systems subjected
1272
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
to external normal or non-normal excitations can be written as follows:
multi-dimensional system for the Poisson white noise process are obtained as follows:
Z˙ (t) = a[Z (t); t] + g(t)(t)
m˙ r [Z (t)] = Qr; 2N E[Z [r−1] (t) ⊗ a[Z (t); t]] r [k] +
E[Y k ]: gr; k (t) mr−k [Z (t)] ⊗ I2N
(15)
or alternatively dZ (t) = a[Z (t); t] dt + g(t) d(t);
(16)
where Z ; a, and g are three time dependent vectors of order 2N . By using the Kronecker algebra [24], the series expansion of the increment N(Z [r] ) can be written as follows: N(Z [r] ) = (∇TZ ⊗ Z [r] ) dZ r T 1 [r] [k] + (∇[k] Z ⊗ Z ) dZ k!
(17)
k=1
In this equation the symbol ⊗ means Kronecker product; I2N is the identity matrix of order 2N; ∇TZ is the following row di4erential vector: @ @ @ T (18) ∇Z = ··· @Z1 @Z2 @Z2N the exponent in square bracket means Kronecker power, that is, Z [2] = Z ⊗ Z ;
Z [3] = Z ⊗ Z ⊗ Z ; · · ·
(19)
and gr; k (t) =
1 Qr; 2N (Qr−1; 2N ⊗ I2N ) · · · k! [k−1] [r−k] (Qr−k+1; 2N ⊗ I2N )(I2N ⊗ g [k] (t)): (20)
In the previous equations the matrix Qr; 2N is a matrix of order (2N )r × (2N )r de1ned as follows [25]: Qr; 2N =
r−1
E(2N )r−s ; 2N s ;
(22)
Notice that the stationary rth moment vector mr [Z ] of the response can be evaluated by solving the following algebraic equation: Qr; 2N E[Z [r−1] ⊗ a[Z ]] r [k] + gr; k (mr−k [Z ] ⊗ I2N ) E[Y k ] = 0;
(23)
k=1
k=2
= Qr; 2N (Z [r−1] ⊗ a[Z ; t]) dt r [k] + gr; k (t)(Z [r−k] ⊗ I2N )(d)k :
k=1
(21)
s=0
where Eq; p denotes a permutation matrix of order (qp)×(qp) consisting of q×p arrays of p×q dimensional elementary submatrices E is which have one on the (i; s)th position and zero in all other positions. By use of the stochastic average of both sides of Eq. (17) and dividing the result by dt, the di4erential equations governing the evolution of the moments of a
obtained from Eq. (22) assuming that this equation is not time dependent. Moreover, since a zero mean Gaussian white noise process can be considered as a delta-correlated process up to the second order, for this input the sum in Eqs. (22) and (23) needs to be performed only for k = 2. 3. Non-Gaussian closures for polynomial non-linearity 3.1. One-dimensional systems For simplicity’s sake, let us assume that the deterministic non-linear function a[Z(t); t] is a polynomial of pth order: p a[Z(t); t] = aj Z j (t): (24) j=1
By means of this equation, we can rewrite Eq. (13) as follows: p m˙ r [Z(t)] = r aj (t)mr+j−1 [Z(t)] j=1 r
+
Lr; k mr−k [Z(t)]gk (t) E[Y r ]:
(25)
k=1
As a consequence of the system non-linearity in the 1rst summation of the right member of Eq. (25) moments until the (r + p − 1)th order appear. It follows that the di4erential equation governing the evolution of statistics of order r will now include statistics of higher order than r. As an example, if
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
cubic non-linearity is considered, the di4erential Eq. (25) contains moments of (r + 2)th order. The numerical solution of Eq. (25) can be obtained only if a Gaussian or NGC is performed [17]. The most used closure techniques, such as the cumulant and the quasi-moment neglect closure, lead to a set of non-linear equations. This fact represents the major drawback in the numerical applications related to a non-uniqueness of solutions, as emphasized by Wojtkiewicz et al. [26]. In this paper, a new closure technique is proposed, based on the modi1cation of the well-known quasi-moment neglect closure and taking advantage on the great accuracy of the MCS in the evaluation of the second-order statistics of the response. The proposed closure technique leads to a set of linear equations for the evaluation of higher order statistical moments. In fact, the relationships between quasi-moments and statistical moments can be written as follows [17]: s−2
bs [Z] = ms [Z] + (−1)s
Bs; i Zs−i mi [Z]
i=0;2 (s=even) i=3;5 (s=odd)
+ (−1)s ×
s
Bs; i Zs−i
i=0;2 (s=even) i=3;5 (s=odd)
i r=1
i! r mi−r [Z] Z [Z] ; (−1) r!(i − r)! r
s ¿ 2; (26) in which bs [Z] are the quasi-moments of the response process and Bs; s = (−1)s ;
Bs; i = (1 + i)Bs−1; i+1 − Bs−1; i−1 : (27)
If the mean Z and the standard deviation Z are approximated by MCS, that is Z = m1 ∼ (28) = ˜Z ; Z = m2 − m21 ∼ = ˜Z ; and if a quasi-moments closure of order R is adopted, the statistical moments of order higher than R can be evaluated from Eq. (26) assuming the quasi-moments of order higher than R equal to zero (bs = 0, s ¿ R)
1273
and obtaining the following linear relationships: 1 + (−1)s [(Bs; 0 + Bs; 2 ) ˜sZ ] ms [Z] = (−1)s+1 2 s−2 s+1 Bs; i ˜Zs−i mi [Z] + (−1) i=4;6 (s=even) i=3;5 (s=odd) s
+ (−1)s+1
×
i−3
i=4;6 (s=even) i=3;5 (s=odd)
(−1)r
r=1 i
Bs; i ˜Zs−i
i! mi−r [Z] ˜rZ r!(i − r)!
i! ˜i r!(i − r)! Z r=i−2 i(i − 1) 2 i−2 ; + (−1)i−2 ˜Z ˜Z 2
+
(−1)r
(s ¿ R): (29)
By substituting Eq. (29) into Eq. (25), we obtain a set of R − 2 linear di4erential equations, whose unknowns are the statistical moments of the process Z(t) of order greater than two, up to order R. The unique solution of the set of linear di4erential equations can then be obtained by traditional numerical procedures. Moreover, since the number of realization required to approximate accurately the 1rst two statistical moments of the response is relatively small, the proposed method is very competitive to 1nd additional probabilistic information in terms of moments of higher order than two. This information is best exploited in order to evaluate the p.d.f. of the response process, as explained in the next section. This procedure allows avoiding the large number of realizations required by the Monte Carlo method to accurately produce the p.d.f. 3.2. Multi-dimensional systems Let us consider the non-linearity to be in polynomial form, this happens, for example, in the case of geometrically non-linear structural systems modeled by using the 1nite element method (see, e.g., [27,28]). This formulation could be used also in structural problems in which the vector a[Z (t); t] can be conveniently expanded by Taylor’s series expansion or when an
1274
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
equivalent non-linearization is employed [29]. The series representation of the vector a[Z (t); t] by adopting Kronecker algebra can be written as follows: p a[Z (t); t] = Aj Z [ j] (t); (30)
×(#i; i−2 Vec( ˜ Z ) ⊗ ˜Z[i−2] ) i ; (s ¿ R); (−1)j #i; j ˜[i] + Z j=i−2
j=0
where Aj are matrices of order 2N × (2N )j given as 1 T[ j] (31) Aj = [∇Z ⊗ a[Z (t); t]Z =0 : j! For polynomial non-linearity, it is possible to evaluate the statistical moments of rth order of the response by solving the following set of di4erential equations (3 6 r 6 R): p m˙ r [Z (t)] = Aj; r mr+j−1 [Z (t)] j=1
+
r
[k] gr; k (t)(mr−k [Z (t)] ⊗ I2N )
where ˜Z and ˜ Z are the mean vector and the cross-covariance matrix evaluated by MCS and Vec() means vectorialized form of the matrix in parentheses, i.e. a column formed by all columns of the matrix in parentheses in such way that these columns are written below each other. The singular matrices Pj; 2N are given as P1; 2N = I2N ;
Pj; 2N = Qj; 2N (Qj−1; 2N ⊗ I2N ) · · ·
where (r times);
(33)
the symbol ⊕ being the Kronecker summation. It follows that for polynomial non-linearity, in order to evaluate the moments of the response process until Rth order, it is necessary to know moments of higher order than R. However, by using a quasi-moments NGC of Rth order it is possible to relate the moments higher order than R with moments until Rth order by means of the following linear relationship [17,27]: 1 + (−1)s 1=2 mz = (−1)s+1 ( ˜ Z )[s] 2
−1=2 [2] ˜ ˜ × [Bs; 0 + Bs; 2 ( Z ) Vec( Z )] + (−1)s+1 s−2 1=2 [s] −1=2 ˜ ×( Z ) Bs; i ( ˜ Z )[i] mi i=4;6 (s=even) i=3;5 (s=odd)
+
s i=4;6 (s=even) i=3;5 (s=odd)
−1=2 [i]
Bs; i ( ˜ Z
)
Pi; 2N i!
i−3 × (−1)j (#i; j mi−j ⊗ ˜[ j] ) + (−1)i−2 j=1
[ j−2] (Q2; 2N ⊗ I2N )
(32)
Aj; r = Aj ⊕ Aj ⊕ · · · ⊕ Aj
Z
P2; 2N = Q2; 2N ;
P3; 2N = Q3; 2N (Q2; 2N ⊗ I2N );
k=1
× E[Y k ];
(34)
(35)
and Bs; i can be evaluated in the recursive form as follows [17,27]: [s] Bs; s = (−1)s I2N ;
Bs; i = (I2N ⊗ Bs−1; i+1 )Ri+1 − (I2N ⊗ Bs−1; i−1 ) (36) with [i−1] Ri = Rˆ i [Vec(I2N ) ⊗ I2N ]; [i] [i−2] Rˆ i = I2N ⊗ (I2N + E2N; 2N ⊗ I2N [i−3] + E2N; (2N )2 ⊗ I2N + · · · + E2N; (2N )i−1 ):
(37)
In Eq. (34), the coeIcient #i; j can be evaluated recursively as [27] #i; j = #i−1; j + #i−1; j−1 ;
#i; 0 = #i; i = 1:
(38)
By substituting Eq. (34) into Eq. (32), a set of linear di4erential equations whose unknowns are the statistical moments of the response until the Rth order is obtained. Notice that the proposed approach only requires the MCS estimation until the second-order statistics, because the quasi-moment neglect closure is adopted. In order to avoid the use of a closure technique, the hierarchy of moment di4erential equations could be
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
closed by estimating via MCS the moments of order greater than two also, as proposed in [16]. For higher order polynomial non-linearities, this fact can demand a larger number of samples in MCS, in order to better estimate higher order moments. For example, if the order of non-linearity is p and the statistical moments up to Rth order are of interest, the Rth order moment equation contains unknown moments until the order R + p − 1. To close the hierarchy of moment equations without the use of any closure technique, lower order moments up to (p − 1)th order must be estimated by MCS. On the contrary, if the quasi-moment neglect closure is adopted, only moments of the 1rstand second-order are required by MCS, regardless of the order of the non-linearity. However, the use of a closure technique introduces additional approximations with respect to that given by MCS in estimating the 1rst- and second-order moments. A further remark is necessary concerning the sensitivity analysis required to give the rate of change of the higher order moments with the lower 1rst- and second-order moments, as proposed in [16]. This investigation is of great interest in order to know how the error in estimating the 1rst- and second-order moments propagates over the higher order moments. In fact, the error of the higher order moments obtained by the proposed method has two sources: (1) the statistical uncertainty in the estimated 1rst- and second-order moments by MCS; and (2) the mathematical structure of the relationships given by the closure technique adopted, which relate the higher order moments to the lower ones, that is the non-linear functional dependence of the resultant moment equations with respect to the 1rst- and second-order moments. An exhaustive treatment of this problem should be performed. However, the error resulting from the estimation of the 1rst- and second-order statistical moments by MCS can be assumed reasonably negligible with respect to that given by the approximations due to the adopted quasi-moment neglect closure. 4. Approximate probability density function of the response 4.1. One-dimensional systems The problem to 1nd the p.d.f. of a random variable starting from partial probabilistic information
1275
in terms of statistical moments is generally solved by using the expansions of p.d.f. of the response based on the so-called A-type Gram–Charlier series [17,18]. The coeIcient of this series expansion are the quasi-moments of the random variable, that can be easily evaluated in terms of statistical moments by means of Eq. (26). Unfortunately, the accuracy and convergence of this series expansion is not guaranteed, indeed this expansion sometimes leads to negative values at low probability regions which do not have physical meaning, making the reliability analysis greatly inaccurate. To avoid this de1ciency, some non-negative probability functions of exponential form have been recently proposed [19,20], based on the C-type Gram– Charlier (CGC) series expansion. The latter is more rapidly convergent than the A-type one and ensures a good accuracy of p.d.f., with positive values over the whole region. Recently, the relationships between the quasi-moments of the response and the coeIcients of the C-type series expansion have been determined, overcoming the major drawback in using this expansion in practical structural problems. The C-type series expansion can be written as follows: ∞ 1 pZ (z) = N exp j [Z]Hj (z ∗ ) ; (39) j j=1
where N is a normalizing constant, Hj (z ∗ ) is the jth Hermite polynomial and z ∗ = (z − ˜Z )= ˜Z is the standardized variable. Notice that the exponential form of this series ensures that the p.d.f. is always positive; moreover, only few terms of this series permit to approximate accurately those probability densities which strongly deviate from the Gaussian condition. Expression (39) requires the knowledge of the time-dependent coeIcients j [Z], that are related to the quasi-moments of the response process, as will be shown later. In order to evaluate the coeIcients j [Z] let us consider the expansion (39) truncated at M th term. By performing an integration by part, it can easily be seen that the following equality holds: ∞ dpZ (z) Hi−1 (z ∗ ) dz d z∗ −∞ ∞ dHi−1 (z ∗ ) =− pZ (z) d z: (40) d z∗ −∞
1276
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
By recalling the mathematical de1nition of stochastic average and by substituting in the 1rst integral of Eq. (40) the probability density function (39) truncated at M th term and by taking into account the following properties of the Hermite polynomials: bi [Z] Ci [Z] = E[Hi (Z ∗ )] = i ; ˜Z d Hi (z ∗ ) = iHi−1 (z ∗ ); d z∗
Once an approximation of the unknowns coeIcients i is obtained, the normalization of the approximate p.d.f. of the response process is required. It has to be emphasized that the approximation of the p.d.f. with M terms requires the evaluation of moments and quasi-moments of the response process up to R = (2M − 2)th order.
Hi+1 (z ∗ ) = z ∗ Hi (z ∗ ) − iHi−1 (z ∗ );
In this section the CGC series expansion is adopted to approximate the multi-dimensional p.d.f. of the response. However, to use this series it is needed to use the explicit relationships between the quasi-moments and the coeIcients of this series recently obtained [20] for multi-dimensional uncorrelated random process. To this aim, let us write the CGC series expansion in the space of the vector process Z truncated at the M th term as follows: M 1 pZ (z) = N exp H T (z ∗ )j [Z ] ; (47) j j
∗
∗
Hi−1 (z )Hj−1 (z ) i+j−2
=
1 &i−1; j−1; s Hs (z ∗ ); s!
(41)
s=0
where Ci [Z] are the so-called Hermite moments and &p; q; r is de1ned as &p; q; r = '(# + 1) '(( − # + 1) #=p; q; r
if p + q + r = even
( = (p + q + r)=2 ¿ p; q; r; 0 elsewhere;
(42)
' being the Gamma function, it is possible to rewrite Eq. (40) as follows: i+j−2 M 1 &i−1; j−1; s Cs [Z] j [Z] s! j=1
s=0
4.2. Multi-dimensional systems
j=1
where j [Z ] is the jth coeIcient vector of the C-type series expansion of order (2N )j , Hj (z ∗ ) is a vector of order (2N )j of the multi-dimensional Hermite polynomials (see Ref. [20]) and z ∗ is the standardized variable vector de1ned as follows: −1=2
z ∗ = ˜ Z
(z − ˜Z ):
(48)
Eq. (43) represents a set of M algebraic equations, that can be written in the following matrix form:
The coeIcients j [Z ] of the CGC series can be evalM uated as the solution of a system of m = i=1 (2N )i algebraic equations, that can be written in the following matrix form [20]:
X = x;
X = ;
= − (i − 1)Ci−2 [Z];
(i = 1; 2; : : : ; M ):
(43)
(44)
where is the M -vector of unknown coeIcients of the truncated CGC series expansion, while the M elements of vector x and the M (M + 1)=2 elements of the symmetric matrix X can be written as follows: xi = −(i − 1)Ci−2 [Z]; i+j−2 1 Xij = &i−1; j−1; s Cs [Z]: s!
(45)
s=0
It is useful to remember that x1 = 0;
b0 [Z] = 1;
b1 [Z] = 0;
b2 [Z] = 0:
(46)
where
X11 X12 X21 X22 X = . .. .. . XN 1 XN 2 1 2 = . ; .. N
(49) · · · X1N · · · X2N . .. . .. · · · XNN
;
=
1 2 .. .
;
N
(50)
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
X being a symmetric matrix, of order m × m whose elements Xij (of order (2N )i × (2N )j ) are given as T Xij = E[Hi−1 (z ∗ )Hj−1 (z ∗ )] ⊗ I2N
= i−1; j−1 ⊗ I2N :
(51)
It can easily be shown that Xji = XijT ; moreover, i−1; j−1 (of order (2N )i−1 × (2N )j−1 ) can be related to the Hermite moment vectors −1=2 Cj [Z ] = ( ˜ Z )[j] bj [Z ]
(52)
as follows: Vec{i−1; j−1 } = E[Hj−1 (z ∗ ) ⊗ Hi−1 (z ∗ )] i+j−2
=
1 i−1; j−1; k Ck [Z ]; k!
(53)
k=0
where the matrix p; q; r of order (2N )p+q × (2N )r is given as 1 p; q; r = [Hq (x) ⊗ Hp (x)]HrT (x) (2+)N R2N
1 ×exp − xT x d x1 d x2 : : : d x2N ; (54) 2 which can be evaluated by using Eq. (42) and taking into account that an element of the Hermite polynomial vectors Hi (z ∗ ) is the product of one-dimensional Hermite polynomials of the type Hj (zr∗ )Hk (zs∗ ) · · · Hl (zt∗ ), with j + k + · · · + l = i. Finally, in Eq. (49) ith subvector i of vector can be obtained as follows: T i = −Vec[(Ci−2 [Z ] ⊗ I2N )Qi−1; 2N ]:
(55)
It has to be emphasized that the coeIcient matrix X of the linear algebraic system given by Eq. (49) is singular. This is due to the fact that, by applying Kronecker algebra, some equations in (49) and consequently some elements of the unknown vectors j are equal between each other. To overcome this drawback, one can reduce the system equation (49) by eliminating both the repeated equations and repeated unknown coeIcients. Moreover, as stated in [20] once an approximation of the unknown vector coeIcients j is obtained, the normalization of the approximate p.d.f. given by Eq. (47) requires a certain caution; in fact, by inspection of this equation, it is obvious that the exponential function can violate the limit condition limz→∞pZ (z) = 0.
1277
In this case, an accurate normalization of the approximate p.d.f. can be performed by an ‘ad hoc’ truncation of the exponential function in the range of interest, which can be chosen on the basis of probabilistic information given by the estimated 1rst- and second-order moments. 5. Application 5.1. Example 1 In order to validate the proposed method, let us consider as an application the following non-linear half-oscillator: ˙ = −-Z(t) − .Z 3 (t) + (t); Z(t)
(56)
where - and . are positive constant and (t) is a Poisson white process, with the impulse occurrence de1ned by the parameter . The impulse amplitude is assumed as a standardized Gaussian random variable Y . Then, the moments of odd order are zero and the moments of even order are given by E[Y 2k ] = 1 × 3 × · · · × (2k − 1) = (2k − 1)!!: (57) Taking into account that a[Z(t); t] = −-Z − .Z 3 is an odd function, the moments of odd order of the response process are zero and only the moments of even order must be evaluated. The stationary problem was solved for the following parameter set: - = 1, . = 1, = 30. The approximate variance ˜2Z = 2:289 of the response process has been evaluated by MCS with 100 samples only and exploiting the ergodicity of the response process. The obtained value is a very good approximation of the assumed exact value determined by MCS with 109 samples. In order to obtain the stationary statistical moments of the response, the proposed procedure applied to this example leads to the following set of algebraic equations: m˙ r [Z(t)] = rmr [Z(t)] + rmr+2 [Z(t)] r + Lr; k mr−k [Z(t)] E[Y r ] = 0; k=2;4
r = 4; 6; : : : :
(58)
The proposed procedure has been performed, by considering di4erent values of R, that de1ne the order of the NGC. The moment equations until Rth
1278
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
15
Table 1 Stationary moments of the response evaluated by the proposed method for di4erent order of NGC compared with assumed exact values by MCS (109 samples)
m4[Z]
10
5 MCS NGC (R=6) NGC (R=12)
0 0
2
4
(a)
6
8
10
t 150
m8 [Z]
m10 [Z]
R=4 e (%)
13.74 (8.96)
176.1 (58.50)
—
—
R=6 e (%)
13.08 (3.73)
112.4 (1.17)
715 (48.34)
—
R=8 e (%)
12.85 (1.90)
112.7 (1.44)
1441 (4.12)
30,591 (31.52)
R = 10 e (%)
12.76 (1.19)
112.7 (1.44)
1435 (3.68)
24,460 (5.16)
R = 12 e (%)
12.73 (0.95)
112.8 (1.53)
1432 (3.47)
24,440 (5.07)
Exact (MCS)
12.61
111.1
1384
23,260
50 0.25
MCS NGC (R=6) NGC (R=12)
0 (b)
m6 [Z]
0
2
4
6
8
0.2
10
t
p Z (z )
m6[Z]
100
m4 [Z]
Fig. 1. Comparison of the proposed method with the non-stationary MCS results for di4erent order of NGC: (a) fourth-order moment of the response process and (b) sixth-order moment of the response process.
0.15 0.1 MCS CGC (M=4) CGC (M=6) CGC (M=12)
0.05 0 -5
0
2.5
5
z 1
0.1
p Z (z)
order have been considered. By setting to zero the quasi-moment of order (R + 2), the moment of (R + 2)th order has been approximated. The resulting linear algebraic moment equations have been solved and the stationary values of the approximate moments have been compared with non-stationary MCS results with 10,000 samples. In Figs. 1a and b the fourth-order and the sixth-order approximate stationary moments are plotted for different values of R, revealing the good approximation of results compared with non-stationary solutions by simulation (with zero initial condition). In Table 1, the stationary moments of the response process up to 10th order evaluated by the proposed method for di4erent closure orders are compared with MCS results with 109 samples, assumed as exact
-2.5
(a)
0.01 MCS CGC (M=4) CGC (M=6) CGC (M=12)
0.001
0.0001 -5
(b)
-2.5
0
2.5
5
z
Fig. 2. Stationary p.d.f. of the response process for di4erent values of M (CGC) and comparison with MCS: (a) decimal scale and (b) logarithmic scale.
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
solutions. This comparison reveals the good performances of the method and the accuracy increases as the closure order increases. The percentage error used is de1ned as |mk − mek | e(%) = × 100; (59) mek
0.025 0.02
D(M)
0.015
where the apex e means the exact value. Once the moments of the response process are evaluated, the p.d.f. has been approximated by the proposed method that utilizes the CGC series expansion. The coeIcients of this series have been determined by solving a linear algebraic system in terms of the quasi-moment of the response process. The analysis has been performed for di4erent values of M , that is the quasi-moments (or moments) up to R=(2M −2)th order are required.
0.01 0.005 0 2
4
6
8
10
12
M Fig. 3. Kullback–Leibler cross-entropy measure D between the approximated stationary p.d.f. (evaluated by the proposed method) and the stationary p.d.f. by MCS (with 109 samples) for di4erent values of M .
0.2
0.008
0.16
E[u2u2]
0.006
E[u4]
1279
0.004
0.12 0.08
0.002
MCS NGC (R=6) NGC (R=10)
0 0
2
4
6
(a)
8
MCS NGC (R=6) NGC (R=10)
0.04 0
10
0
2
4
(b)
t
6
8
10
t
30
E[u 4]
20
10 MCS NGC (R=6) NGC (R=10)
0 0 (c)
2
4
6
8
10
t
Fig. 4. Comparison of the proposed method with the non-stationary MCS results for di4erent order of NGC: fourth-order moments of the response process (a–c).
1280
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283 2.5
0.0016
2
E[u2u 4]
E[u 6]
0.0012
0.0008
0.0004
MCS NGC (R=6) NGC (R=10)
0 2
4
6
8
10
t
0
2
4
(c)
6
8
10
t
0.03
600
0.02
400
E[u 6]
E[u 4u 2]
MCS NGC (R=6) NGC (R=10)
0.5
200
0.01 MCS NGC (R=6) NGC (R=10)
0 0
(b)
1
0 0
(a)
1.5
2
4
6
8
t
MCS NGC (R=6) NGC (R=10)
0 0
10
(d)
2
4
6
8
10
t
Fig. 5. Comparison of the proposed method with the non-stationary MCS results for di4erent order of NGC: sixth-order moments of the response process (a–d).
In Fig. 2a the stationary p.d.f. of the response process evaluated by the proposed method for di4erent values of M are plotted and set in comparison with MCS with 109 samples. Note that even if few probabilistic information in terms of quasi-moments are considered, that is for M = 4; 6, satisfactory results are obtained. Better results are obtained by setting M =12, shown by comparison with simulation. In this case the method is capable to capture the bimodal shape of the probability distribution. The good accuracy in estimating small probabilities in the tail of the response distribution is shown in Fig. 2b, where the probability densities obtained by the proposed method are compared with MCS in logarithmic scale. In order to investigate the convergence of the CGC series expansion as M increases, a measure of diver-
gence between the approximate stationary probability distributions and the exact one (by MCS with 109 samples) is introduced [30]: D(M ) = D[pZM (z); pZe (z)] M +∞ pZ (z) M d z; = pZ (z) log pZe (z) −∞
(60)
which is the Kullback–Leibler cross-entropy measure. In Eq. (60), pZM (z) means the Gram–Charlier series approximation with M terms, while pZe (z) means the assumed exact solution by MCS with 109 samples. In Fig. 3 the cross-entropy (or divergence) D(M ) as a function of M is reported, revealing the fast convergence of the Gram–Charlier series approximation to the exact distribution.
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283 2
0.25
MCS CGC (M=2) CGC (M=4) CGC (M=6)
1.6
MCS CGC (M=2) CGC (M=4)
0.2 0.15 . p(u)
1.2
p (u)
1281
0.8
0.1
0.4
0.05
0
0
-1
-0.5
0
0.5
1
u
(a)
-9
-6
-3
(a)
0 .
3
6
9
6
9
u 1
1 0.1
0.1 .
0.001
0.001 MCS CGC (M=2) CGC (M=4) CGC (M=6)
0.0001 1E-005
(b)
0.01
p(u)
p (u)
0.01
MCS CGC (M=2) CGC (M=4)
0.0001 1E-005 -9
-1
-0.5
0
0.5
1
(b)
-6
-3
0 .
3
u
u
Fig. 6. Marginal p.d.f. p(u) of the displacement u and the velocity u˙ for di4erent values of M (CGC) and comparison with MCS: (a) decimal scale and (b) logarithmic scale.
5.2. Example 2 Let us consider the DuIng oscillator, which is a single-degree-of-freedom system with linear damping and cubic non-linear spring, governed by the following equation of motion: uW + 22!u˙ + (1 + .u)!2 u = ((t);
(61)
where ! is the natural frequency of the undamped linear structure (with . = 0), 2 is the damping ratio, . is a parameter representing the level of non-linearity and (t) is a Poisson white process, de1ned in the previous example. According to the proposed procedure the approximate variance of the displacement and of the velocity has been evaluated by MCS with 1000 samples only and exploiting the ergodicity of the response process, giving ˜2u = 0:056 and ˜2u˙ = 2:962, respectively, assuming ! = 2, 2 = 0:2, . = 100, ( = 0:4 and = 30. Introducing the state variable vector Z T = [u; u], ˙ the
Fig. 7. Marginal p.d.f. p(u) ˙ of the displacement u and the velocity u˙ for di4erent values of M (CGC) and comparison with MCS: (a) decimal scale and (b) logarithmic scale.
statistical moments of the response have been evaluated by applying the proposed procedure for several closure orders. Note that, since the input intensities of odd order are zeros and because of the antisymmetric characteristics of the non-linearity, the stationary odd moments of the response are zero. Figs. 4a–c and 5a–d show fourth- and sixth-order approximate stationary moments of the response for di4erent values of R, which are compared with non-stationary MCS results with 10,000 samples. The comparison shows the good accuracy of the method as the closure order increases. Once the moments of the response are evaluated the stationary joint p.d.f. of the displacement and the velocity have been approximated by using the CGC series expansion given by Eq. (47) truncated at M th term. The coeIcients of this series have been evaluated by solving Eq. (49). To the end of comparing the results obtained through this those given by approach with MCS the marginal p.d.f. p(u) and p(u) ˙
1282
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
are plotted in Figs. 6a, b and Figs. 7a, b for decimal and logarithmic scales. The logarithmic plots show the good accuracy in estimating small probabilities in the tail of the response distribution, which are particularly important for reliability of structural systems. 6. Conclusions A method for evaluating in approximate form the probability density function of non-linear multi-degrees of freedom systems under non-normal delta correlated input has been presented. The proposed method, based on the C-type Gram–Charlier series expansion of probability density function, requires the following steps: (a) the evaluation of mean vector and covariance matrix of non-linear systems by means of Monte Carlo Simulation; (b) the use of the quasi-moment non-Gaussian closure in order to truncate the in1nite hierarchy of statistical moments which appear in equations for polynomial non-linearities; (c) the evaluation of moments higher order than two by solving a set of linear equations; and (d) the evaluation of the probability density function of the response process by the knowledge of the statistical moments. The proposed approach takes advantage of: (a) the great accuracy of the Monte Carlo simulation in evaluating the 1rst two moments of the response process with a small number of input samples; (b) the possibility to obtain statistical moments of higher order by solving a set of linear equations; and (c) the positivity of the probability density function in all the range of interest and the capacity to capture with great accuracy the tails of the density. References [1] G. Ricciardi, Random vibration of beam under moving loads, J. Eng. Mech. ASCE 120 (11) (1994) 2361–2380. [2] C.C. Tung, Random response of highway bridges to vehicle loads, J. Eng. Mech. Div. ASCE 93 (5) (1967) 79–94. [3] Y.K. Lin, Application of non-stationary shot noise in the study of system response to non-stationary excitations, J. Appl. Mech. ASME 30 (1963) 555–558. [4] J.B. Roberts, System response to random impulse, J. Sound Vib. 24 (1) (1972) 23–34. [5] M. Grigoriu, Applied Non-Gaussian Processes, Prentice-Hall, Englewood Cli4s, NJ, 1995.
[6] I.I. Gikhman, A.W. Skorohod, Stochastic Di4erential Equations, Springer, Berlin, 1972. [7] H.U. KWoylWuoglu, S.R.K. Nielsen, R. Iwankiewicz, Reliability of non-linear oscillators subjected to Poisson driven impulses, J. Sound Vib. 176 (1994) 19–33. [8] H.U. KWoylWuoglu, S.R.K. Nielsen, A.S. Cakmak, Fast cell-to-cell mapping (path integration) for non-linear white noise and Poisson driven systems, Struct. Saf. 17 (1995) 151–165. [9] H.U. KWoylWuoglu, S.R.K. Nielsen, R. Iwankiewicz, Response and reliability of Poisson driven systems by path integration, J. Eng. Mech. ASCE 121 (1995) 117–129. [10] A. Tylikowski, W. Marowski, Vibration of a non-linear single degree of freedom system due to Poissonian impulse excitation, Int. J. Non-Linear Mech. 21 (1986) 229–238. [11] M. Grigoriu, A partial di4erential equation for the characteristic function of the response of non-linear systems to additive Poisson white noise, J. Sound Vib. 198 (1996) 193–202. [12] M. Di Paola, G. Falsone, Itˆo and Stratonovich integrals for delta-correlated processes, Prob. Eng. Mech. 8 (1993) 197– 208. [13] M. Di Paola, M. Vasta, Stochastic integro-di4erential and di4erential equations of non-linear systems excited by parametric Poisson pulses, Int. J. Non-Linear Mech. 8 (1997) 855–862. [14] G.I. Schueller, (Guest Ed.), A state-of-the-art report on computational stochastic mechanics—simulation techniques (Special issue IASSAR Report), Prob. Eng. Mech. 12(4) (1997) 203–229. [15] M. Vasta, G.I. Schueller, Phase space reduction in stochastic dynamics, J. Eng. Mech. (ASCE) 126 (2000) 626–632. [16] M. Grigoriu, Moment closure by Monte Carlo simulation and moment sensitivity factors, Int. J. Non-Linear Mech. 34 (1999) 739–748. [17] G. Muscolino, Response of non-linear structural systems under Gaussian or non-Gaussian 1ltered input, in: F. Casciati (Ed.), Dynamic Motion, Chaotic and Stochastic Behaviour, Springer, Wien, 1993, pp. 203–299. [18] N.C. Hampl, G.I. Schueller, Probability density of the response of non-linear structures under stochastic dynamic excitation, Prob. Eng. Mech. 4 (1996) 2–9. [19] G. Muscolino, G. Ricciardi, M. Vasta, Stationary and non-stationary probability density function for non-linear oscillators, Int. J. Non-linear Mech. 32 (1997) 1051–1064. [20] G. Muscolino, G. Ricciardi, Probability density function of MDOF structural systems under non-normal delta-correlated inputs, Comput. Meth. Appl. Mech. Eng. 168 (1999) 121– 133. [21] M. Grigoriu, White noise processes, J. Eng. Mech. ASCE 113 (5) (1987) 757–765. [22] R.L. Stratonovich, Topics in the Theory of Random Noise, Vol. 1, Gordon and Breach Science Publishers, New York, 1963. [23] M. Di Paola, G. Falsone, Stochastic dynamics of non-linear systems driven by non-normal delta correlated process, J. Appl. Mech. ASME 60 (1993) 141–148.
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283 [24] M. Fiedler, Special Matrices and their Applications in Numerical Mathematics, Martinus Nijho4, Dordrecht, The Netherlands, 1986. [25] F. Ma, Extension of second moment analysis to vector valued and matrix valued functions, Int. J. Non-Linear Mech. 22 (1987) 251–260. [26] S.F. Wojtkiewicz, B.F. Spencer Jr., L.A. Bergman, On the cumulant-neglect closure methods in stochastic dynamics, Int. J. Non-Linear Mech. 31 (5) (1996) 657–684.
1283
[27] M. Di Paola, G. Falsone, G. Muscolino, Random analysis of geometrically non-linear FE modelled structures under seismic actions, Struct. Saf. 8 (1990) 209–222. [28] M. Di Paola, G. Muscolino, Di4erential moment equations of FE modelled structures with geometrical non-linearities, Int. J. Non-Linear Mech. 22 (1990) 363–373. [29] J.B. Roberts, P.D. Spanos, Random Vibration and Statistical Linearization, Wiley, New York, 1989. [30] S. Kullback, Information Theory and Statistics, Wiley, New York, 1959.
International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
Monte Carlo simulation in the stochastic analysis of non-linear systems under external stationary Poisson white noise input G. Muscolino ∗ , G. Ricciardi, P. Cacciola Dipartimento di Costruzioni e Tecnologie Avanzate, University of Messina, Salita Sperone 31, S. Agata, I-98166 Messina, Italy
Abstract A method for the evaluation of the probability density function (p.d.f.) of the response process of non-linear systems under external stationary Poisson white noise excitation is presented. The method takes advantage of the great accuracy of the Monte Carlo simulation (MCS) in evaluating the 1rst two moments of the response process by considering just few samples. The quasi-moment neglect closure is used to close the in1nite hierarchy of the moment di4erential equations of the response process. Moreover, in order to determine the higher order statistical moments of the response, the second-order probabilistic information given by MCS in conjunction with the quasi-moment neglect closure leads to a set of linear di4erential equations. The quasi-moments up to a given order are used as partial probabilistic information on the response process in order to 1nd the p.d.f. by means of the C-type Gram–Charlier series expansion. ? 2002 Elsevier Science Ltd. All rights reserved. Keywords: Non-linear stochastic dynamics; Stationary Poisson process; Monte Carlo simulation; Non-Gaussian probability density function; Quasi-moment neglect closure
1. Introduction A broad class of random excitations on structures can be modeled as Gaussian white noise or 1ltered Gaussian white noise processes. Then both analytical and numerical methods have been recently developed to determine the response of structures driven by these types of excitations. However, when excitations are characterized by impulsive loading, these representations may not be satisfactory. Moving loads traveling on a bridge [1,2], earthquake shake on structures [3] and wave action on ships [4] are examples of this type of excitations. ∗
Corresponding author. Tel.: +390-90-676-5620; fax: +39090-395022. E-mail address: [email protected] (G. Muscolino).
In these cases the Poisson white noise and the 1ltered Poisson white noise processes are more realistic models [5]. Methods for the evaluation of the response of linear and non-linear systems under non-Gaussian delta-correlated processes have been recently proposed by several authors. Some of these are based on the theory of Markov processes and determine the probability density function (p.d.f.) of the response process as solution of the generalized Kolmogorov equation [6]. Numerical methods have been developed to solve this di4erential equation such as the Petrov–Galerkin method [7] and the path integral method [8,9]. Recent works on this subject utilize the characteristic function to describe the response process. These methods require the solution of the partial di4erential equation for the characteristic function obtained by
0020-7462/03/$ - see front matter ? 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 4 6 2 ( 0 2 ) 0 0 0 7 2 - 0
1270
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
taking the Fourier transform of the generalized Kolmogorov equation [10] or starting by the generalized Itˆo’s di4erential rule for Poisson white noise input [11]. Several approaches based on the moment equation method have been proposed [12,13]. Unfortunately, these methods su4er from the problem associated with the approximation techniques used to close the in1nite hierarchy of equations. The Monte Carlo simulation (MCS) constitutes the universal alternative method of solution when the system under study is large and strongly non-linear. This method requires the generation of samples of the input process, the deterministic analysis for each of them and the evaluation of the resultant response statistics, which are used to develop estimates of various characteristics of the response process. This approach can provide accurate solutions for a large class of stochastic problems [14]. Unfortunately, the disadvantages of MCS method are that it is generally time-consuming, when the statistics of higher order or the p.d.f. of the response process are of interest. Moreover, the MCS is rather ineIcient for estimating small probabilities such as those in the tail of the response distribution, which are particularly important for the reliability of structural systems. However, the great accuracy of MCS in estimating lower order statistics of the response has been recently exploited in stochastic dynamics of non-linear systems [15,16]. Following this strategy, in this paper an alternative approach to evaluate the stationary p.d.f. of non-linear systems subjected to external stationary Poisson white noise input is presented. The proposed approach takes advantage of: (a) the great accuracy of the MCS in evaluating the 1rst two moments of the response process with a small number of input samples; (b) the non-Gaussian closure (NGC) technique in conjunction with the moment equation approach, which leads to a set of di4erential equations in which the non-linearity is due to the presence of mean and variance of the response process [17]. As a consequence of the evaluation of the mean and variance by MCS, the di4erential equations of moments of higher order than two become a closed set of linear di4erential equations. In the stationary case, the set of linear di4erential equations becomes an algebraic one. Once a good approximation of the statistical moments up to a given order
is achieved, the p.d.f. of the response process can be easily determined. In this paper the so-called C-type Gram–Charlier (CGC) series expansion [18] is used, whose coeIcients are determined as linear functions of quasi-moments of the response process [19,20]. 2. Preliminary considerations 2.1. One-dimensional systems For clarity’s sake and in order to evidence the advantages of the proposed formulation, a non-linear one-dimensional system is 1rst studied. Indeed, as it will be shown later, the multi-dimensional case is formally similar to the one-dimensional case if the Kronecker algebra is used. Let us consider the di4erential equation describing a one-dimensional system subjected to a purely external Poisson white noise process, given by ˙ = a[Z(t); t] + g(t)(t); Z(t)
(1)
where Z(t) is the response process, the dot means time di4erentiation, a[Z(t); t] is a deterministic non-linear function of the response process Z(t) and g(t) is a deterministic time-dependent function. In Eq. (1), (t) is a scalar Poisson white noise process [21], that is a sequence of impulses with random amplitude and arriving at random times, given by (t) =
N (t)
Yk (t − tk );
(2)
k=1
the amplitudes {Yk } being a family of random variables, mutually independent and independent of the time instants tk , with prescribed distribution pY (y),
(t) the Dirac’s delta function, tk the random time arrivals. In Eq. (2), N (t) is a counting Poisson process with parameter . The probabilistic characterization of the Poisson white noise process (t) can be given in terms of its correlation functions [22] Cr [(t1 ); (t2 ); : : : ; (tr )] = E[Y r ] (t1 − t2 ) (t1 − t3 ) · · · (t1 − tr );
(3)
where E[ • ] means stochastic average. It follows that (t) is a non-normal delta-correlated process, with intensity coeIcients related to the probabilistic
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
1271
characterization of the random variables Yk and of the counting Poisson process N (t). The Poisson white noise process (t) is not mean square Riemann integrable and, consequently, Eq. (1) has no traditional mathematical meaning. The latter can be considered to be formally equivalent to the following generalized Itˆo-type stochastic di4erential equation:
Since (dZ)r is an in1nitesimal of order dt, it follows that the series expansion of the increment N(Z r ) has to be written in the following form [1,23]: r 1 k r d (Z ) = rZ r−1 dZ N(Z r ) = k!
dZ(t) = a[Z(t); t] dt + g(t) d(t);
where Lr; k = r!=[k!(r − k)!]. Eq. (11) is obtained observing that neglecting terms of order higher than dt the following relationships hold:
(4)
where d(t) is an increment of the compound Poisson process (t), de1ned by (t) =
N (t)
Yk U (t − tk );
(5)
k=1
U (t) being the unit step function. The compound Poisson process is characterized by the following correlation functions [1]: Cr [(t1 ); (t2 ); : : : ; (tr )] = E[Y r ] min(t1 ; t2 ; : : : ; tr ):
(6)
If we proceed formally, that is in terms of correlation functions, it can be shown that the Poisson white noise process can be considered the derivative of the compound Poisson process. In fact, the correlation functions of an increment of the compound Poisson process are given by the following relationships: Cr [d(t1 ); d(t2 ); : : : ; d(tr )] = E[Y r ] (t1 − t2 ) (t1 − t3 ) · · ·
(t1 − tr ) dt1 dt2 : : : dtr :
(7)
Accounting for the relationships between correlation functions and moments, putting t1 = t2 = · · · = tr = t and neglecting in1nitesimals of order higher than dt, we obtain mr [d(t)] ≡ kr [d(t)] = E[Y r ] dt;
(8)
where mr [ • ] means rth statistical moment and kr [ • ] means the rth cumulant. Moreover, by taking into account the non-anticipating property, the following relationships hold: E[(Z(t))k (d(t))r ] = E[Y r ]mk [Z(t)] dt:
(9)
Neglecting in1nitesimals of order higher than dt, we have (dZ)r = gr (t)(d)r ;
r ¿ 2:
(10)
k=1
+
r
Lr; k Z r−k gk (t)(d)k ;
(11)
k=2
d(Z r ) = rZ r−1 dZ; d k (Z r ) = k!Lr; k Z r−k (dZ)k = k!Lr; k Z r−k gk (t)(d)k ; k = 2; 3; : : : ; r: (12) Using the stochastic average of both sides of Eq. (11), taking into account Eqs. (8) – (10) and dividing the result by dt, we obtain the di4erential equation governing the evolution of the rth moment of the response m˙ r [Z(t)] = rE[Z r−1 (t)a[Z(t); t]] r + Lr; k mr−k [Z(t)]gk (t) E[Y k ]:
(13)
k=1
Notice that the stationary rth order moment of the response can be evaluated by solving the following algebraic equation: rE[Z r−1 a[Z]] r + Lr; k mr−k [Z]gk E[Y k ] = 0;
(14)
k=1
obtained from Eq. (13) observing that this equation is not time dependent. Moreover, since a zero mean Gaussian white process can be considered as delta-correlated up to the second order, for this input the sum in Eqs. (13) and (14) need to be performed only for k = 2. In this case q2 = E[Y 2 ] is usually de1ned as the strength of the white process. 2.2. Multi-dimensional systems Let us consider a N -degree-of-freedom structural system, whose nodal displacements are the elements of vector u(t). By introducing the vector of state variables Z (t) = [uT (t) u˙T (t)]T , the equation of motion of the non-linear N -DOF structural systems subjected
1272
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
to external normal or non-normal excitations can be written as follows:
multi-dimensional system for the Poisson white noise process are obtained as follows:
Z˙ (t) = a[Z (t); t] + g(t)(t)
m˙ r [Z (t)] = Qr; 2N E[Z [r−1] (t) ⊗ a[Z (t); t]] r [k] +
E[Y k ]: gr; k (t) mr−k [Z (t)] ⊗ I2N
(15)
or alternatively dZ (t) = a[Z (t); t] dt + g(t) d(t);
(16)
where Z ; a, and g are three time dependent vectors of order 2N . By using the Kronecker algebra [24], the series expansion of the increment N(Z [r] ) can be written as follows: N(Z [r] ) = (∇TZ ⊗ Z [r] ) dZ r T 1 [r] [k] + (∇[k] Z ⊗ Z ) dZ k!
(17)
k=1
In this equation the symbol ⊗ means Kronecker product; I2N is the identity matrix of order 2N; ∇TZ is the following row di4erential vector: @ @ @ T (18) ∇Z = ··· @Z1 @Z2 @Z2N the exponent in square bracket means Kronecker power, that is, Z [2] = Z ⊗ Z ;
Z [3] = Z ⊗ Z ⊗ Z ; · · ·
(19)
and gr; k (t) =
1 Qr; 2N (Qr−1; 2N ⊗ I2N ) · · · k! [k−1] [r−k] (Qr−k+1; 2N ⊗ I2N )(I2N ⊗ g [k] (t)): (20)
In the previous equations the matrix Qr; 2N is a matrix of order (2N )r × (2N )r de1ned as follows [25]: Qr; 2N =
r−1
E(2N )r−s ; 2N s ;
(22)
Notice that the stationary rth moment vector mr [Z ] of the response can be evaluated by solving the following algebraic equation: Qr; 2N E[Z [r−1] ⊗ a[Z ]] r [k] + gr; k (mr−k [Z ] ⊗ I2N ) E[Y k ] = 0;
(23)
k=1
k=2
= Qr; 2N (Z [r−1] ⊗ a[Z ; t]) dt r [k] + gr; k (t)(Z [r−k] ⊗ I2N )(d)k :
k=1
(21)
s=0
where Eq; p denotes a permutation matrix of order (qp)×(qp) consisting of q×p arrays of p×q dimensional elementary submatrices E is which have one on the (i; s)th position and zero in all other positions. By use of the stochastic average of both sides of Eq. (17) and dividing the result by dt, the di4erential equations governing the evolution of the moments of a
obtained from Eq. (22) assuming that this equation is not time dependent. Moreover, since a zero mean Gaussian white noise process can be considered as a delta-correlated process up to the second order, for this input the sum in Eqs. (22) and (23) needs to be performed only for k = 2. 3. Non-Gaussian closures for polynomial non-linearity 3.1. One-dimensional systems For simplicity’s sake, let us assume that the deterministic non-linear function a[Z(t); t] is a polynomial of pth order: p a[Z(t); t] = aj Z j (t): (24) j=1
By means of this equation, we can rewrite Eq. (13) as follows: p m˙ r [Z(t)] = r aj (t)mr+j−1 [Z(t)] j=1 r
+
Lr; k mr−k [Z(t)]gk (t) E[Y r ]:
(25)
k=1
As a consequence of the system non-linearity in the 1rst summation of the right member of Eq. (25) moments until the (r + p − 1)th order appear. It follows that the di4erential equation governing the evolution of statistics of order r will now include statistics of higher order than r. As an example, if
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
cubic non-linearity is considered, the di4erential Eq. (25) contains moments of (r + 2)th order. The numerical solution of Eq. (25) can be obtained only if a Gaussian or NGC is performed [17]. The most used closure techniques, such as the cumulant and the quasi-moment neglect closure, lead to a set of non-linear equations. This fact represents the major drawback in the numerical applications related to a non-uniqueness of solutions, as emphasized by Wojtkiewicz et al. [26]. In this paper, a new closure technique is proposed, based on the modi1cation of the well-known quasi-moment neglect closure and taking advantage on the great accuracy of the MCS in the evaluation of the second-order statistics of the response. The proposed closure technique leads to a set of linear equations for the evaluation of higher order statistical moments. In fact, the relationships between quasi-moments and statistical moments can be written as follows [17]: s−2
bs [Z] = ms [Z] + (−1)s
Bs; i Zs−i mi [Z]
i=0;2 (s=even) i=3;5 (s=odd)
+ (−1)s ×
s
Bs; i Zs−i
i=0;2 (s=even) i=3;5 (s=odd)
i r=1
i! r mi−r [Z] Z [Z] ; (−1) r!(i − r)! r
s ¿ 2; (26) in which bs [Z] are the quasi-moments of the response process and Bs; s = (−1)s ;
Bs; i = (1 + i)Bs−1; i+1 − Bs−1; i−1 : (27)
If the mean Z and the standard deviation Z are approximated by MCS, that is Z = m1 ∼ (28) = ˜Z ; Z = m2 − m21 ∼ = ˜Z ; and if a quasi-moments closure of order R is adopted, the statistical moments of order higher than R can be evaluated from Eq. (26) assuming the quasi-moments of order higher than R equal to zero (bs = 0, s ¿ R)
1273
and obtaining the following linear relationships: 1 + (−1)s [(Bs; 0 + Bs; 2 ) ˜sZ ] ms [Z] = (−1)s+1 2 s−2 s+1 Bs; i ˜Zs−i mi [Z] + (−1) i=4;6 (s=even) i=3;5 (s=odd) s
+ (−1)s+1
×
i−3
i=4;6 (s=even) i=3;5 (s=odd)
(−1)r
r=1 i
Bs; i ˜Zs−i
i! mi−r [Z] ˜rZ r!(i − r)!
i! ˜i r!(i − r)! Z r=i−2 i(i − 1) 2 i−2 ; + (−1)i−2 ˜Z ˜Z 2
+
(−1)r
(s ¿ R): (29)
By substituting Eq. (29) into Eq. (25), we obtain a set of R − 2 linear di4erential equations, whose unknowns are the statistical moments of the process Z(t) of order greater than two, up to order R. The unique solution of the set of linear di4erential equations can then be obtained by traditional numerical procedures. Moreover, since the number of realization required to approximate accurately the 1rst two statistical moments of the response is relatively small, the proposed method is very competitive to 1nd additional probabilistic information in terms of moments of higher order than two. This information is best exploited in order to evaluate the p.d.f. of the response process, as explained in the next section. This procedure allows avoiding the large number of realizations required by the Monte Carlo method to accurately produce the p.d.f. 3.2. Multi-dimensional systems Let us consider the non-linearity to be in polynomial form, this happens, for example, in the case of geometrically non-linear structural systems modeled by using the 1nite element method (see, e.g., [27,28]). This formulation could be used also in structural problems in which the vector a[Z (t); t] can be conveniently expanded by Taylor’s series expansion or when an
1274
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
equivalent non-linearization is employed [29]. The series representation of the vector a[Z (t); t] by adopting Kronecker algebra can be written as follows: p a[Z (t); t] = Aj Z [ j] (t); (30)
×(#i; i−2 Vec( ˜ Z ) ⊗ ˜Z[i−2] ) i ; (s ¿ R); (−1)j #i; j ˜[i] + Z j=i−2
j=0
where Aj are matrices of order 2N × (2N )j given as 1 T[ j] (31) Aj = [∇Z ⊗ a[Z (t); t]Z =0 : j! For polynomial non-linearity, it is possible to evaluate the statistical moments of rth order of the response by solving the following set of di4erential equations (3 6 r 6 R): p m˙ r [Z (t)] = Aj; r mr+j−1 [Z (t)] j=1
+
r
[k] gr; k (t)(mr−k [Z (t)] ⊗ I2N )
where ˜Z and ˜ Z are the mean vector and the cross-covariance matrix evaluated by MCS and Vec() means vectorialized form of the matrix in parentheses, i.e. a column formed by all columns of the matrix in parentheses in such way that these columns are written below each other. The singular matrices Pj; 2N are given as P1; 2N = I2N ;
Pj; 2N = Qj; 2N (Qj−1; 2N ⊗ I2N ) · · ·
where (r times);
(33)
the symbol ⊕ being the Kronecker summation. It follows that for polynomial non-linearity, in order to evaluate the moments of the response process until Rth order, it is necessary to know moments of higher order than R. However, by using a quasi-moments NGC of Rth order it is possible to relate the moments higher order than R with moments until Rth order by means of the following linear relationship [17,27]: 1 + (−1)s 1=2 mz = (−1)s+1 ( ˜ Z )[s] 2
−1=2 [2] ˜ ˜ × [Bs; 0 + Bs; 2 ( Z ) Vec( Z )] + (−1)s+1 s−2 1=2 [s] −1=2 ˜ ×( Z ) Bs; i ( ˜ Z )[i] mi i=4;6 (s=even) i=3;5 (s=odd)
+
s i=4;6 (s=even) i=3;5 (s=odd)
−1=2 [i]
Bs; i ( ˜ Z
)
Pi; 2N i!
i−3 × (−1)j (#i; j mi−j ⊗ ˜[ j] ) + (−1)i−2 j=1
[ j−2] (Q2; 2N ⊗ I2N )
(32)
Aj; r = Aj ⊕ Aj ⊕ · · · ⊕ Aj
Z
P2; 2N = Q2; 2N ;
P3; 2N = Q3; 2N (Q2; 2N ⊗ I2N );
k=1
× E[Y k ];
(34)
(35)
and Bs; i can be evaluated in the recursive form as follows [17,27]: [s] Bs; s = (−1)s I2N ;
Bs; i = (I2N ⊗ Bs−1; i+1 )Ri+1 − (I2N ⊗ Bs−1; i−1 ) (36) with [i−1] Ri = Rˆ i [Vec(I2N ) ⊗ I2N ]; [i] [i−2] Rˆ i = I2N ⊗ (I2N + E2N; 2N ⊗ I2N [i−3] + E2N; (2N )2 ⊗ I2N + · · · + E2N; (2N )i−1 ):
(37)
In Eq. (34), the coeIcient #i; j can be evaluated recursively as [27] #i; j = #i−1; j + #i−1; j−1 ;
#i; 0 = #i; i = 1:
(38)
By substituting Eq. (34) into Eq. (32), a set of linear di4erential equations whose unknowns are the statistical moments of the response until the Rth order is obtained. Notice that the proposed approach only requires the MCS estimation until the second-order statistics, because the quasi-moment neglect closure is adopted. In order to avoid the use of a closure technique, the hierarchy of moment di4erential equations could be
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
closed by estimating via MCS the moments of order greater than two also, as proposed in [16]. For higher order polynomial non-linearities, this fact can demand a larger number of samples in MCS, in order to better estimate higher order moments. For example, if the order of non-linearity is p and the statistical moments up to Rth order are of interest, the Rth order moment equation contains unknown moments until the order R + p − 1. To close the hierarchy of moment equations without the use of any closure technique, lower order moments up to (p − 1)th order must be estimated by MCS. On the contrary, if the quasi-moment neglect closure is adopted, only moments of the 1rstand second-order are required by MCS, regardless of the order of the non-linearity. However, the use of a closure technique introduces additional approximations with respect to that given by MCS in estimating the 1rst- and second-order moments. A further remark is necessary concerning the sensitivity analysis required to give the rate of change of the higher order moments with the lower 1rst- and second-order moments, as proposed in [16]. This investigation is of great interest in order to know how the error in estimating the 1rst- and second-order moments propagates over the higher order moments. In fact, the error of the higher order moments obtained by the proposed method has two sources: (1) the statistical uncertainty in the estimated 1rst- and second-order moments by MCS; and (2) the mathematical structure of the relationships given by the closure technique adopted, which relate the higher order moments to the lower ones, that is the non-linear functional dependence of the resultant moment equations with respect to the 1rst- and second-order moments. An exhaustive treatment of this problem should be performed. However, the error resulting from the estimation of the 1rst- and second-order statistical moments by MCS can be assumed reasonably negligible with respect to that given by the approximations due to the adopted quasi-moment neglect closure. 4. Approximate probability density function of the response 4.1. One-dimensional systems The problem to 1nd the p.d.f. of a random variable starting from partial probabilistic information
1275
in terms of statistical moments is generally solved by using the expansions of p.d.f. of the response based on the so-called A-type Gram–Charlier series [17,18]. The coeIcient of this series expansion are the quasi-moments of the random variable, that can be easily evaluated in terms of statistical moments by means of Eq. (26). Unfortunately, the accuracy and convergence of this series expansion is not guaranteed, indeed this expansion sometimes leads to negative values at low probability regions which do not have physical meaning, making the reliability analysis greatly inaccurate. To avoid this de1ciency, some non-negative probability functions of exponential form have been recently proposed [19,20], based on the C-type Gram– Charlier (CGC) series expansion. The latter is more rapidly convergent than the A-type one and ensures a good accuracy of p.d.f., with positive values over the whole region. Recently, the relationships between the quasi-moments of the response and the coeIcients of the C-type series expansion have been determined, overcoming the major drawback in using this expansion in practical structural problems. The C-type series expansion can be written as follows: ∞ 1 pZ (z) = N exp j [Z]Hj (z ∗ ) ; (39) j j=1
where N is a normalizing constant, Hj (z ∗ ) is the jth Hermite polynomial and z ∗ = (z − ˜Z )= ˜Z is the standardized variable. Notice that the exponential form of this series ensures that the p.d.f. is always positive; moreover, only few terms of this series permit to approximate accurately those probability densities which strongly deviate from the Gaussian condition. Expression (39) requires the knowledge of the time-dependent coeIcients j [Z], that are related to the quasi-moments of the response process, as will be shown later. In order to evaluate the coeIcients j [Z] let us consider the expansion (39) truncated at M th term. By performing an integration by part, it can easily be seen that the following equality holds: ∞ dpZ (z) Hi−1 (z ∗ ) dz d z∗ −∞ ∞ dHi−1 (z ∗ ) =− pZ (z) d z: (40) d z∗ −∞
1276
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
By recalling the mathematical de1nition of stochastic average and by substituting in the 1rst integral of Eq. (40) the probability density function (39) truncated at M th term and by taking into account the following properties of the Hermite polynomials: bi [Z] Ci [Z] = E[Hi (Z ∗ )] = i ; ˜Z d Hi (z ∗ ) = iHi−1 (z ∗ ); d z∗
Once an approximation of the unknowns coeIcients i is obtained, the normalization of the approximate p.d.f. of the response process is required. It has to be emphasized that the approximation of the p.d.f. with M terms requires the evaluation of moments and quasi-moments of the response process up to R = (2M − 2)th order.
Hi+1 (z ∗ ) = z ∗ Hi (z ∗ ) − iHi−1 (z ∗ );
In this section the CGC series expansion is adopted to approximate the multi-dimensional p.d.f. of the response. However, to use this series it is needed to use the explicit relationships between the quasi-moments and the coeIcients of this series recently obtained [20] for multi-dimensional uncorrelated random process. To this aim, let us write the CGC series expansion in the space of the vector process Z truncated at the M th term as follows: M 1 pZ (z) = N exp H T (z ∗ )j [Z ] ; (47) j j
∗
∗
Hi−1 (z )Hj−1 (z ) i+j−2
=
1 &i−1; j−1; s Hs (z ∗ ); s!
(41)
s=0
where Ci [Z] are the so-called Hermite moments and &p; q; r is de1ned as &p; q; r = '(# + 1) '(( − # + 1) #=p; q; r
if p + q + r = even
( = (p + q + r)=2 ¿ p; q; r; 0 elsewhere;
(42)
' being the Gamma function, it is possible to rewrite Eq. (40) as follows: i+j−2 M 1 &i−1; j−1; s Cs [Z] j [Z] s! j=1
s=0
4.2. Multi-dimensional systems
j=1
where j [Z ] is the jth coeIcient vector of the C-type series expansion of order (2N )j , Hj (z ∗ ) is a vector of order (2N )j of the multi-dimensional Hermite polynomials (see Ref. [20]) and z ∗ is the standardized variable vector de1ned as follows: −1=2
z ∗ = ˜ Z
(z − ˜Z ):
(48)
Eq. (43) represents a set of M algebraic equations, that can be written in the following matrix form:
The coeIcients j [Z ] of the CGC series can be evalM uated as the solution of a system of m = i=1 (2N )i algebraic equations, that can be written in the following matrix form [20]:
X = x;
X = ;
= − (i − 1)Ci−2 [Z];
(i = 1; 2; : : : ; M ):
(43)
(44)
where is the M -vector of unknown coeIcients of the truncated CGC series expansion, while the M elements of vector x and the M (M + 1)=2 elements of the symmetric matrix X can be written as follows: xi = −(i − 1)Ci−2 [Z]; i+j−2 1 Xij = &i−1; j−1; s Cs [Z]: s!
(45)
s=0
It is useful to remember that x1 = 0;
b0 [Z] = 1;
b1 [Z] = 0;
b2 [Z] = 0:
(46)
where
X11 X12 X21 X22 X = . .. .. . XN 1 XN 2 1 2 = . ; .. N
(49) · · · X1N · · · X2N . .. . .. · · · XNN
;
=
1 2 .. .
;
N
(50)
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
X being a symmetric matrix, of order m × m whose elements Xij (of order (2N )i × (2N )j ) are given as T Xij = E[Hi−1 (z ∗ )Hj−1 (z ∗ )] ⊗ I2N
= i−1; j−1 ⊗ I2N :
(51)
It can easily be shown that Xji = XijT ; moreover, i−1; j−1 (of order (2N )i−1 × (2N )j−1 ) can be related to the Hermite moment vectors −1=2 Cj [Z ] = ( ˜ Z )[j] bj [Z ]
(52)
as follows: Vec{i−1; j−1 } = E[Hj−1 (z ∗ ) ⊗ Hi−1 (z ∗ )] i+j−2
=
1 i−1; j−1; k Ck [Z ]; k!
(53)
k=0
where the matrix p; q; r of order (2N )p+q × (2N )r is given as 1 p; q; r = [Hq (x) ⊗ Hp (x)]HrT (x) (2+)N R2N
1 ×exp − xT x d x1 d x2 : : : d x2N ; (54) 2 which can be evaluated by using Eq. (42) and taking into account that an element of the Hermite polynomial vectors Hi (z ∗ ) is the product of one-dimensional Hermite polynomials of the type Hj (zr∗ )Hk (zs∗ ) · · · Hl (zt∗ ), with j + k + · · · + l = i. Finally, in Eq. (49) ith subvector i of vector can be obtained as follows: T i = −Vec[(Ci−2 [Z ] ⊗ I2N )Qi−1; 2N ]:
(55)
It has to be emphasized that the coeIcient matrix X of the linear algebraic system given by Eq. (49) is singular. This is due to the fact that, by applying Kronecker algebra, some equations in (49) and consequently some elements of the unknown vectors j are equal between each other. To overcome this drawback, one can reduce the system equation (49) by eliminating both the repeated equations and repeated unknown coeIcients. Moreover, as stated in [20] once an approximation of the unknown vector coeIcients j is obtained, the normalization of the approximate p.d.f. given by Eq. (47) requires a certain caution; in fact, by inspection of this equation, it is obvious that the exponential function can violate the limit condition limz→∞pZ (z) = 0.
1277
In this case, an accurate normalization of the approximate p.d.f. can be performed by an ‘ad hoc’ truncation of the exponential function in the range of interest, which can be chosen on the basis of probabilistic information given by the estimated 1rst- and second-order moments. 5. Application 5.1. Example 1 In order to validate the proposed method, let us consider as an application the following non-linear half-oscillator: ˙ = −-Z(t) − .Z 3 (t) + (t); Z(t)
(56)
where - and . are positive constant and (t) is a Poisson white process, with the impulse occurrence de1ned by the parameter . The impulse amplitude is assumed as a standardized Gaussian random variable Y . Then, the moments of odd order are zero and the moments of even order are given by E[Y 2k ] = 1 × 3 × · · · × (2k − 1) = (2k − 1)!!: (57) Taking into account that a[Z(t); t] = −-Z − .Z 3 is an odd function, the moments of odd order of the response process are zero and only the moments of even order must be evaluated. The stationary problem was solved for the following parameter set: - = 1, . = 1, = 30. The approximate variance ˜2Z = 2:289 of the response process has been evaluated by MCS with 100 samples only and exploiting the ergodicity of the response process. The obtained value is a very good approximation of the assumed exact value determined by MCS with 109 samples. In order to obtain the stationary statistical moments of the response, the proposed procedure applied to this example leads to the following set of algebraic equations: m˙ r [Z(t)] = rmr [Z(t)] + rmr+2 [Z(t)] r + Lr; k mr−k [Z(t)] E[Y r ] = 0; k=2;4
r = 4; 6; : : : :
(58)
The proposed procedure has been performed, by considering di4erent values of R, that de1ne the order of the NGC. The moment equations until Rth
1278
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
15
Table 1 Stationary moments of the response evaluated by the proposed method for di4erent order of NGC compared with assumed exact values by MCS (109 samples)
m4[Z]
10
5 MCS NGC (R=6) NGC (R=12)
0 0
2
4
(a)
6
8
10
t 150
m8 [Z]
m10 [Z]
R=4 e (%)
13.74 (8.96)
176.1 (58.50)
—
—
R=6 e (%)
13.08 (3.73)
112.4 (1.17)
715 (48.34)
—
R=8 e (%)
12.85 (1.90)
112.7 (1.44)
1441 (4.12)
30,591 (31.52)
R = 10 e (%)
12.76 (1.19)
112.7 (1.44)
1435 (3.68)
24,460 (5.16)
R = 12 e (%)
12.73 (0.95)
112.8 (1.53)
1432 (3.47)
24,440 (5.07)
Exact (MCS)
12.61
111.1
1384
23,260
50 0.25
MCS NGC (R=6) NGC (R=12)
0 (b)
m6 [Z]
0
2
4
6
8
0.2
10
t
p Z (z )
m6[Z]
100
m4 [Z]
Fig. 1. Comparison of the proposed method with the non-stationary MCS results for di4erent order of NGC: (a) fourth-order moment of the response process and (b) sixth-order moment of the response process.
0.15 0.1 MCS CGC (M=4) CGC (M=6) CGC (M=12)
0.05 0 -5
0
2.5
5
z 1
0.1
p Z (z)
order have been considered. By setting to zero the quasi-moment of order (R + 2), the moment of (R + 2)th order has been approximated. The resulting linear algebraic moment equations have been solved and the stationary values of the approximate moments have been compared with non-stationary MCS results with 10,000 samples. In Figs. 1a and b the fourth-order and the sixth-order approximate stationary moments are plotted for different values of R, revealing the good approximation of results compared with non-stationary solutions by simulation (with zero initial condition). In Table 1, the stationary moments of the response process up to 10th order evaluated by the proposed method for di4erent closure orders are compared with MCS results with 109 samples, assumed as exact
-2.5
(a)
0.01 MCS CGC (M=4) CGC (M=6) CGC (M=12)
0.001
0.0001 -5
(b)
-2.5
0
2.5
5
z
Fig. 2. Stationary p.d.f. of the response process for di4erent values of M (CGC) and comparison with MCS: (a) decimal scale and (b) logarithmic scale.
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
solutions. This comparison reveals the good performances of the method and the accuracy increases as the closure order increases. The percentage error used is de1ned as |mk − mek | e(%) = × 100; (59) mek
0.025 0.02
D(M)
0.015
where the apex e means the exact value. Once the moments of the response process are evaluated, the p.d.f. has been approximated by the proposed method that utilizes the CGC series expansion. The coeIcients of this series have been determined by solving a linear algebraic system in terms of the quasi-moment of the response process. The analysis has been performed for di4erent values of M , that is the quasi-moments (or moments) up to R=(2M −2)th order are required.
0.01 0.005 0 2
4
6
8
10
12
M Fig. 3. Kullback–Leibler cross-entropy measure D between the approximated stationary p.d.f. (evaluated by the proposed method) and the stationary p.d.f. by MCS (with 109 samples) for di4erent values of M .
0.2
0.008
0.16
E[u2u2]
0.006
E[u4]
1279
0.004
0.12 0.08
0.002
MCS NGC (R=6) NGC (R=10)
0 0
2
4
6
(a)
8
MCS NGC (R=6) NGC (R=10)
0.04 0
10
0
2
4
(b)
t
6
8
10
t
30
E[u 4]
20
10 MCS NGC (R=6) NGC (R=10)
0 0 (c)
2
4
6
8
10
t
Fig. 4. Comparison of the proposed method with the non-stationary MCS results for di4erent order of NGC: fourth-order moments of the response process (a–c).
1280
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283 2.5
0.0016
2
E[u2u 4]
E[u 6]
0.0012
0.0008
0.0004
MCS NGC (R=6) NGC (R=10)
0 2
4
6
8
10
t
0
2
4
(c)
6
8
10
t
0.03
600
0.02
400
E[u 6]
E[u 4u 2]
MCS NGC (R=6) NGC (R=10)
0.5
200
0.01 MCS NGC (R=6) NGC (R=10)
0 0
(b)
1
0 0
(a)
1.5
2
4
6
8
t
MCS NGC (R=6) NGC (R=10)
0 0
10
(d)
2
4
6
8
10
t
Fig. 5. Comparison of the proposed method with the non-stationary MCS results for di4erent order of NGC: sixth-order moments of the response process (a–d).
In Fig. 2a the stationary p.d.f. of the response process evaluated by the proposed method for di4erent values of M are plotted and set in comparison with MCS with 109 samples. Note that even if few probabilistic information in terms of quasi-moments are considered, that is for M = 4; 6, satisfactory results are obtained. Better results are obtained by setting M =12, shown by comparison with simulation. In this case the method is capable to capture the bimodal shape of the probability distribution. The good accuracy in estimating small probabilities in the tail of the response distribution is shown in Fig. 2b, where the probability densities obtained by the proposed method are compared with MCS in logarithmic scale. In order to investigate the convergence of the CGC series expansion as M increases, a measure of diver-
gence between the approximate stationary probability distributions and the exact one (by MCS with 109 samples) is introduced [30]: D(M ) = D[pZM (z); pZe (z)] M +∞ pZ (z) M d z; = pZ (z) log pZe (z) −∞
(60)
which is the Kullback–Leibler cross-entropy measure. In Eq. (60), pZM (z) means the Gram–Charlier series approximation with M terms, while pZe (z) means the assumed exact solution by MCS with 109 samples. In Fig. 3 the cross-entropy (or divergence) D(M ) as a function of M is reported, revealing the fast convergence of the Gram–Charlier series approximation to the exact distribution.
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283 2
0.25
MCS CGC (M=2) CGC (M=4) CGC (M=6)
1.6
MCS CGC (M=2) CGC (M=4)
0.2 0.15 . p(u)
1.2
p (u)
1281
0.8
0.1
0.4
0.05
0
0
-1
-0.5
0
0.5
1
u
(a)
-9
-6
-3
(a)
0 .
3
6
9
6
9
u 1
1 0.1
0.1 .
0.001
0.001 MCS CGC (M=2) CGC (M=4) CGC (M=6)
0.0001 1E-005
(b)
0.01
p(u)
p (u)
0.01
MCS CGC (M=2) CGC (M=4)
0.0001 1E-005 -9
-1
-0.5
0
0.5
1
(b)
-6
-3
0 .
3
u
u
Fig. 6. Marginal p.d.f. p(u) of the displacement u and the velocity u˙ for di4erent values of M (CGC) and comparison with MCS: (a) decimal scale and (b) logarithmic scale.
5.2. Example 2 Let us consider the DuIng oscillator, which is a single-degree-of-freedom system with linear damping and cubic non-linear spring, governed by the following equation of motion: uW + 22!u˙ + (1 + .u)!2 u = ((t);
(61)
where ! is the natural frequency of the undamped linear structure (with . = 0), 2 is the damping ratio, . is a parameter representing the level of non-linearity and (t) is a Poisson white process, de1ned in the previous example. According to the proposed procedure the approximate variance of the displacement and of the velocity has been evaluated by MCS with 1000 samples only and exploiting the ergodicity of the response process, giving ˜2u = 0:056 and ˜2u˙ = 2:962, respectively, assuming ! = 2, 2 = 0:2, . = 100, ( = 0:4 and = 30. Introducing the state variable vector Z T = [u; u], ˙ the
Fig. 7. Marginal p.d.f. p(u) ˙ of the displacement u and the velocity u˙ for di4erent values of M (CGC) and comparison with MCS: (a) decimal scale and (b) logarithmic scale.
statistical moments of the response have been evaluated by applying the proposed procedure for several closure orders. Note that, since the input intensities of odd order are zeros and because of the antisymmetric characteristics of the non-linearity, the stationary odd moments of the response are zero. Figs. 4a–c and 5a–d show fourth- and sixth-order approximate stationary moments of the response for di4erent values of R, which are compared with non-stationary MCS results with 10,000 samples. The comparison shows the good accuracy of the method as the closure order increases. Once the moments of the response are evaluated the stationary joint p.d.f. of the displacement and the velocity have been approximated by using the CGC series expansion given by Eq. (47) truncated at M th term. The coeIcients of this series have been evaluated by solving Eq. (49). To the end of comparing the results obtained through this those given by approach with MCS the marginal p.d.f. p(u) and p(u) ˙
1282
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283
are plotted in Figs. 6a, b and Figs. 7a, b for decimal and logarithmic scales. The logarithmic plots show the good accuracy in estimating small probabilities in the tail of the response distribution, which are particularly important for reliability of structural systems. 6. Conclusions A method for evaluating in approximate form the probability density function of non-linear multi-degrees of freedom systems under non-normal delta correlated input has been presented. The proposed method, based on the C-type Gram–Charlier series expansion of probability density function, requires the following steps: (a) the evaluation of mean vector and covariance matrix of non-linear systems by means of Monte Carlo Simulation; (b) the use of the quasi-moment non-Gaussian closure in order to truncate the in1nite hierarchy of statistical moments which appear in equations for polynomial non-linearities; (c) the evaluation of moments higher order than two by solving a set of linear equations; and (d) the evaluation of the probability density function of the response process by the knowledge of the statistical moments. The proposed approach takes advantage of: (a) the great accuracy of the Monte Carlo simulation in evaluating the 1rst two moments of the response process with a small number of input samples; (b) the possibility to obtain statistical moments of higher order by solving a set of linear equations; and (c) the positivity of the probability density function in all the range of interest and the capacity to capture with great accuracy the tails of the density. References [1] G. Ricciardi, Random vibration of beam under moving loads, J. Eng. Mech. ASCE 120 (11) (1994) 2361–2380. [2] C.C. Tung, Random response of highway bridges to vehicle loads, J. Eng. Mech. Div. ASCE 93 (5) (1967) 79–94. [3] Y.K. Lin, Application of non-stationary shot noise in the study of system response to non-stationary excitations, J. Appl. Mech. ASME 30 (1963) 555–558. [4] J.B. Roberts, System response to random impulse, J. Sound Vib. 24 (1) (1972) 23–34. [5] M. Grigoriu, Applied Non-Gaussian Processes, Prentice-Hall, Englewood Cli4s, NJ, 1995.
[6] I.I. Gikhman, A.W. Skorohod, Stochastic Di4erential Equations, Springer, Berlin, 1972. [7] H.U. KWoylWuoglu, S.R.K. Nielsen, R. Iwankiewicz, Reliability of non-linear oscillators subjected to Poisson driven impulses, J. Sound Vib. 176 (1994) 19–33. [8] H.U. KWoylWuoglu, S.R.K. Nielsen, A.S. Cakmak, Fast cell-to-cell mapping (path integration) for non-linear white noise and Poisson driven systems, Struct. Saf. 17 (1995) 151–165. [9] H.U. KWoylWuoglu, S.R.K. Nielsen, R. Iwankiewicz, Response and reliability of Poisson driven systems by path integration, J. Eng. Mech. ASCE 121 (1995) 117–129. [10] A. Tylikowski, W. Marowski, Vibration of a non-linear single degree of freedom system due to Poissonian impulse excitation, Int. J. Non-Linear Mech. 21 (1986) 229–238. [11] M. Grigoriu, A partial di4erential equation for the characteristic function of the response of non-linear systems to additive Poisson white noise, J. Sound Vib. 198 (1996) 193–202. [12] M. Di Paola, G. Falsone, Itˆo and Stratonovich integrals for delta-correlated processes, Prob. Eng. Mech. 8 (1993) 197– 208. [13] M. Di Paola, M. Vasta, Stochastic integro-di4erential and di4erential equations of non-linear systems excited by parametric Poisson pulses, Int. J. Non-Linear Mech. 8 (1997) 855–862. [14] G.I. Schueller, (Guest Ed.), A state-of-the-art report on computational stochastic mechanics—simulation techniques (Special issue IASSAR Report), Prob. Eng. Mech. 12(4) (1997) 203–229. [15] M. Vasta, G.I. Schueller, Phase space reduction in stochastic dynamics, J. Eng. Mech. (ASCE) 126 (2000) 626–632. [16] M. Grigoriu, Moment closure by Monte Carlo simulation and moment sensitivity factors, Int. J. Non-Linear Mech. 34 (1999) 739–748. [17] G. Muscolino, Response of non-linear structural systems under Gaussian or non-Gaussian 1ltered input, in: F. Casciati (Ed.), Dynamic Motion, Chaotic and Stochastic Behaviour, Springer, Wien, 1993, pp. 203–299. [18] N.C. Hampl, G.I. Schueller, Probability density of the response of non-linear structures under stochastic dynamic excitation, Prob. Eng. Mech. 4 (1996) 2–9. [19] G. Muscolino, G. Ricciardi, M. Vasta, Stationary and non-stationary probability density function for non-linear oscillators, Int. J. Non-linear Mech. 32 (1997) 1051–1064. [20] G. Muscolino, G. Ricciardi, Probability density function of MDOF structural systems under non-normal delta-correlated inputs, Comput. Meth. Appl. Mech. Eng. 168 (1999) 121– 133. [21] M. Grigoriu, White noise processes, J. Eng. Mech. ASCE 113 (5) (1987) 757–765. [22] R.L. Stratonovich, Topics in the Theory of Random Noise, Vol. 1, Gordon and Breach Science Publishers, New York, 1963. [23] M. Di Paola, G. Falsone, Stochastic dynamics of non-linear systems driven by non-normal delta correlated process, J. Appl. Mech. ASME 60 (1993) 141–148.
G. Muscolino et al. / International Journal of Non-Linear Mechanics 38 (2003) 1269 – 1283 [24] M. Fiedler, Special Matrices and their Applications in Numerical Mathematics, Martinus Nijho4, Dordrecht, The Netherlands, 1986. [25] F. Ma, Extension of second moment analysis to vector valued and matrix valued functions, Int. J. Non-Linear Mech. 22 (1987) 251–260. [26] S.F. Wojtkiewicz, B.F. Spencer Jr., L.A. Bergman, On the cumulant-neglect closure methods in stochastic dynamics, Int. J. Non-Linear Mech. 31 (5) (1996) 657–684.
1283
[27] M. Di Paola, G. Falsone, G. Muscolino, Random analysis of geometrically non-linear FE modelled structures under seismic actions, Struct. Saf. 8 (1990) 209–222. [28] M. Di Paola, G. Muscolino, Di4erential moment equations of FE modelled structures with geometrical non-linearities, Int. J. Non-Linear Mech. 22 (1990) 363–373. [29] J.B. Roberts, P.D. Spanos, Random Vibration and Statistical Linearization, Wiley, New York, 1989. [30] S. Kullback, Information Theory and Statistics, Wiley, New York, 1959.