On the generalization of probabilistic transformation method

Applied Mathematics and Computation 190 (2007) 1284–1289 www.elsevier.com/locate/amc

On the generalization of probabilistic transformation method Seifedine Kadry LaMI-IFMA, France

Abstract The probabilistic transformation method with the finite element analysis is a new technique to solve random differential equation. The advantage of this technique is finding the ‘‘exact’’ expression of the probability density function of the solution when the probability density function of the input is known. However the disadvantage is due to the characteristics (geometrics and materials) of the analyzed structure included in the random differential equation. In this paper, a developed formula is used to generalize this technique by obtaining the ‘‘exact’’ joint probability density function of the solution in any situations, as well as the proposed technique for the non-linear case. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Probabilistic transformation method; FEA; Random matrix; Stochastic differential equation

1. Introduction For several decades, the theory of probability has been used in the theory of random differential equations to model the random properties. The solution of a Stochastic Differential Equation is obtained when evaluating the statistical characteristic of the solution process like the mean, standard deviation, high order moments and the most important characteristic ‘‘the probability density function’’ (p.d.f.). The probabilistic transformation method together with the finite element analysis is a new technique for evaluating the p.d.f. of the random differential equation. As opposed to the stochastic finite element method or perturbation based methods, no series expansion is involved in this expression. The power of this technique is demonstrated in the context of Stochastic Differential Equation [1,2], in the reliability analysis [3] and in the random eigenvalue problem [4]. The main problem of this technique is limited and applicable only for some special cases such as non-linearity, complexity and the heterogeneous in the materials of the structure. In this article, the idea is to evaluate in exact form the joint density function of the solution of random differential equation for any complex structure, based on a developed formula that gives the p.d.f of the inverse of stiffness matrix. A technique also was proposed in the non-linear case for the classic technique.

E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2007.02.009

S. Kadry / Applied Mathematics and Computation 190 (2007) 1284–1289

1285

Fig. 1. Algorithm of PTM-FEM.

2. PTM-FEM technique The probabilistic transformation methods (PTM) evaluate the probability density function (PDF) of the system output by multiplying the input PDF by the Jacobian of the inverse function. The idea of PTM is based on the following formulae [5]: fU ðuÞ

¼ fP ðpÞ  jJ p;U j  1    ¼ fP ðpÞ  ou ouðuÞ;

where p is the input parameter, u is the response (solution) and u1 ðuÞ is the inverse transformation, which can be determined either analytically or numerically. The PTM-FEM technique is a combination between the finite element method and the probabilistic transformation method. The PTM-FEM algorithm is illustrated in Fig. 1. 3. Limitation and extension of PTM-FEM The principal mathematical condition in the probabilistic transformation must be one-to-one. Problems frequently arise when we wish to find the probability distribution of the random function Y ¼ U ðX Þ when X is a continuous random variable and the transformation is not one-to-one. Therefore the two major problems or limitations of PTM are: 1. the transformation function (U) should be bijective, 2. the determinant of Jacobian should be ‘‘not null’’. To avoid these limitations of PTM, we first treat a system in n input and only one output y 1 ¼ uðx1 ; x2 ; . . . ; xn Þ. We can proceed as follows:   U : Rn ! Rn y 1 ¼ uðx1 ; x2 ; . . . ; xn Þ; let with U ðX Þ ¼ y i ¼ xi ; i ¼ 2; . . . ; n: Y ! Y ¼ U ðX Þ

1286

S. Kadry / Applied Mathematics and Computation 190 (2007) 1284–1289

Theorem 1. The function U (above) is reversible, i.e. U 1 9. Proof. The function u  ou  ox    1 0 1 0   jJ j ¼   0 1       0   At least, i exist as

ou oxi

(above) is reversible only of the determinant of the Jacobian is not null.   oxoun   0   ou 6¼ 0:   ¼  ox1  1 0 0 1

6¼ 0 i.e., U 1 9. h

So the algorithm becomes: 1. Use the FEM to express the solution or the response of stochastic differential equation in terms of the input, i.e. the load, young’s modulus and the length of the analyzed structure. 2. If the number of random variables equal to A then add A  1 fictitious variables to transform the input to be equal to the output. 3. Calculate the determinant of the Jacobian using Theorem 1. R 4. Calculate fy 1 ðy 1 Þ ¼ Dðx1 Þ jJ jfx ðxÞ dxi , Dðxi Þ: domain of xi. 5. End.

4. Non-linearity technique In the complex mechanical structure, we always face the problem of non-linearity in the stiffness matrix. The difficulty then arises when the non-linear element is random. The method proposed in this article is to linearize the element of stiffness matrix by using the Variable Changes Technique which in turn allows us to apply the PTM-FEM. Some times the non-linearity present in the numerator of the element of stiffness 2 matrix; for example k i ¼ EIL ; here the non-linearity is in inertia moment I2 which is supposed as a random variable and its probability density function (p.d.f) is known. To linearize I2 we suppose M ¼ I 2 , then we apply the PTM Technique to get the p.d.f of M. We then apply easily the PTM to the stiffness element k i ¼ EM which L is linear in terms of random variable M, some times the non-linearity appears in the denominator of the term of stiffness matrix; for example k i ¼  EI . Here the random variable is the lengths (l1 ; l2 Þ of the analyzed 1 1 l1 þl2

structure. We apply the same technique and suppose for example M ¼ l11 and N ¼ l12 . Then we apply PTM to find the p.d.f of M and N again we suppose L ¼ M þ N and we find the p.d.f of L and finally we suppose T ¼ L1 and we evaluate its p.d.f. In this case, the term of stiffness become linear in terms of the random variable k i ¼ EIT then easily we can apply the PTM to find the p.d.f of ki. 5. Generalization of PTM-FEM The PTM-FEM is a new powerful technique that has been developed for calculating the most important factor in probability which is the probability density function as all the other probabilistic characteristics are calculated from p.d.f.; if we look at the literature we would not find any method to evaluate the p.d.f in a deterministic form. For example, the SFEM, which is the famous method in this domain, can give us only the first and second moment and in some cases the moment order 3 and 4 only. Usually the inverse of stiffness matrix is non-linear specially for a complex matrix. As we see, though the PTM-FEM technique is powerful and good yet it is very difficult to generalize due to the non-homogenize of the structure of stiffness matrix which in turn should be well studied due to the non-homogenize of mechanical structure. The application of PTM-FEM should be studied problem by problem, i.e. each structure has its own PTM-FEM. For this reason, to generalize this technique, we should work on the whole stiffness matrix with the p.d.f of the stiffness matrix and not for the element of stiffness matrix as before. In this section, we

S. Kadry / Applied Mathematics and Computation 190 (2007) 1284–1289

1287

developed a new formula which is very important in the calculation in close form of the p.d.f of stiffness matrix. In linear elastic static, the displacement formulation, by using FEM, gives the structural equilibrium in the form ½KfU g ¼ fF g;

ð1Þ

where [K] is the structure stiffness matrix, {U} is the vector of Nodal Displacements and {F} is the vector of Applied Nodal Forces. Suppose the matrix probability density function of the non-singular stiffness matrix K is given by pK ðKÞ : Rnn ! R: We are interested in the p.d.f of H ¼ K 1 2 Rnn

ð2Þ

because in general from (1), U ¼ K 1 F and the main problem is to find the p.d.f of K1. This implies that we want to obtain the joint probability density functions of all the elements of H. The elements of H are complicated non-linear functions of the elements of K. Therefore, even if the elements of K have Gaussian Distribution, the joint distribution of the elements of H will be non-Gaussian. Also note that H may not have any banded structure even if K is of banded nature. The key step to obtain the p.d.f of the inverse of a random matrix is the calculation of the Jacobian of the non-liner matrix transformation of the above equation. Now, our proposed formula is ðnþ1Þ

J ¼ jKj

:

ð3Þ

This expression of the Jacobian of the inverse matrix transformation plays the central role in the determination of the response of a linear stochastic system. Thereafter, we proof this formula for n ¼ 1 and 2 then we proof it for all n. For a single-degree-of-freedom system ðn ¼ 1Þ, Eq. (1) reduces to ku ¼ f where k; u; f 2 R. Suppose the probability density function (may be non-Gaussian) of the random variable k is given by pk ðkÞ and we are interested in deriving the p.d.f of h ¼ k 1 : The Jacobian of the above transformation   oh 2 J ¼   ¼ jkj : ok Using the Jacobian, the probability density function of h can be obtained using PTM technique, as follows:   1 2 ph ðhÞ ¼ jJ jpk ðkÞ ! ph ðhÞ ¼ jhj pk : h The final step is to multiply the p.d.f of h ¼ k 1 by the load vector which is supposed to be deterministic in this case to obtain the p.d.f of u. i.e. pu ðuÞ ¼ ph ðhÞ  f . For n ¼ 2, suppose the 2  2 stiffness symmetric matrix is given by  a b K¼ : b c Here a, b and c are real scalar variables which can be random. From the direct matrix inversion we get  c b 1 1 H ¼K ¼ : ac  b2 b a

1288

S. Kadry / Applied Mathematics and Computation 190 (2007) 1284–1289

We observe that the structure of the inverse of stiffness matrix represents complex coupled non-linear transformations of the original (random) variables in the stiffness matrix. This is why it is difficult to say that the PTM of univariate case is always applicable so the best and the one deterministic solution to work on is the probability of the stiffness matrix itself. Differentiating the upper triangular part of the A matrix, with respect to the independent elements of the K matrix one, has  2 2 3  c bc b2      1 6 7 ¼ J ¼ oH 4 2bc ac  b2 2ab 5 oK  ðacb2 Þ  b2 ab a2         ac  b2 2ab   2bc ac  b2  2  2bc 2ab  c2 bc b       ¼     2 2 2 2 6 2 6 2 6 2  ab a ab ðac  b Þ ðac  b Þ b ðac  b Þ b a ¼ ¼

1 ðac  b2 Þ6 1 ðac  b2 Þ

6

½a3 c3  a2 b2 c2 þ 2a2 b2 c2 þ 2a2 b2 c2  2ab4 c  2ab4 c þ ab4 c þ b6  3

3

ðac  b2 Þ ¼ ðac  b2 Þ

ð2þ1Þ

¼ ðac  b2 Þ

¼ jKj

3

General proof The matrix differential of (2) oH ¼ K 1 oKK 1 :

ð4Þ

It is sometimes convenient when we work with matrix differentiation to rearrange the elements of an m  n matrix A ¼ faij g, in the form on an mn-dimensional column vector. The conventional way of doing so is to successively stack the first, second, . . . , nth columns a1 ; a2 ; . . . ; an of A under each other, giving the mn-dimensional column vector: ða1 a2    an ÞT , this mn-dimensional column vector is referred to as the vec of A. A special case when the matrix A is symmetric we need to stack only the elements that are on or below the diagonal. In this case, vec becomes vec h (h for half). Example: if A  a11 a12 A¼ a21 a22 2 3 2 3 a11 a11 6 a12 7 7 4 5 ) vecðAÞ ¼ 6 4 a21 5, and if A symmetric, i.e. a12 ¼ a21 , ) vec hðAÞ ¼ a12 : a22 a22 Applying [6, Theorem 16.2.1] to Eq. (4) vecðoH Þ ¼ vecðK 1  K 1 Þ vecðoKÞ;

ð5Þ

where  is the Kronecker symbol. We can easily demonstrate that Eq. (4) become vec hðoH Þ ¼ ½H n ðK 1  K 1 ÞGn  vec hðoKÞ;

ð6Þ

where Gn (is called the duplication matrix) can be described in terms of its rows or columns. For i P j, the ½ðj  1Þn þ ith and ½ði  1Þn þ jth rows of Gn equal the ½ðj  1Þðn  j=2Þ þ ith row of I nðnþ1Þ=2 , that is, they are equal to the ½nðn þ 1Þ=2-dimensional row vector whose ½ðj  1Þðn  j=2Þ þ ith element is 1 and whose remaining elements are 0. and (for i P j), the ½ðj  1Þðn  j=2Þ þ ith column of Gn is an n2 dimensional column vector whose ½ðj  1Þn þ ith and ½ði  1Þn þ jth elements are 1 and whose remaining n2  1 (if i ¼ j) or n2  2 (if i > j) elements are 0. Hn is the inverse matrix of Gn. The Jacobian associated with the above transformation is simply the determinant of the matrix H n ðK 1  K 1 ÞGn , again by using [6, Theorem 16.4.2] one obtains

S. Kadry / Applied Mathematics and Computation 190 (2007) 1284–1289

jH n ðK 1  K 1 ÞGn j ¼ jKj

nþ1

1289

:

This is the generalization of the simple case of the PTM-FEM. This expression of the Jacobian of the inverse matrix transformation plays the central role in the determination of the response of a linear stochastic system. In the series of papers of Soize [7,8], the p.d.f of stiffness matrix has a Gamma distribution and the mean of the stiffness matrix can be evaluated from the stochastic finite element method or from Monte-Carlo simulation. Once, the p.d.f of the inverse of stiffness matrix is known, then using (3) the calculation of the p.d.f of the joint probability density function of the response is straightforward fU ðuÞ ¼ jJ j

ðnþ1Þ

 fK ðkÞ:

6. Conclusion In this paper, the generalization of the probabilistic transformation method with finite element analysis is developed. The main disadvantage of the PTM-FEM technique is applicable on a special case, especially when the analyzed structure is homogenize, i.e. when the form of the stiffness matrix is simple. So the idea of generalization is to work with the whole matrix and not with its element. The main difficulty to work with the whole stiffness matrix is in finding the p.d.f of the inverse of stiffness matrix which becomes solved by using our developed formulae. References [1] M. El-Tawil, W. El-Tahhan, A. Hussein, A proposed technique of SFEM on solving ordinary random differential equation, J. Appl. Math. Comput. (2005) 35–47. [2] S. Kadry, R. Younes. Etude probabiliste d’un syste`me me´canique a` parame`tres incertains par une technique base´e sur la me´thode de transformation, in: Proceeding of CANCAM 2004. [3] S. Kadry, A. Chateauneuf, K. El-Tawil, Random eigenvalue problem of Stochastic Systems, in: Proceedings of the 8th International Conference on Computational Structure Technology, Civil-Comp Press, Spain, 2006. [4] S. Kadry, A. Chateauneuf, K. El-Tawil, One-dimensional transformation method in reliability analysis, in: Proceedings of the 8th International Conference on Computational Structure Technology, Civil-Comp Press, Spain, 2006. [5] A. Papoulis, Probability, Random Variables and Stochastic Processes, fourth ed., McGraw-Hill, Boston, USA, 2002. [6] A. Gupta, D. Nagar, Matrix Variate Distributions, Monographs & Surveys in Pure & Applied Mathematics, Chapman & Hall/CRC, London, 2000. [7] C. Soize, Random matrix theory and non-parametric model of random uncertainties in vibration analysis, J. Sound Vib. 263 (4) (2003) 893–916. [8] C. Soize, Non-gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators, Comput. Methods Appl. Mech. Eng. 195 (1-3) (2006) 26–64.