# Numerical Procedures for Random Differential Equations

Authors: Mohamed Ben Said; Lahcen Azrar; Driss Sarsri
Journal: Journal of Applied Mathematics
 (2018)

---

## Abstract

Some methodological approaches based on generalized polynomial chaos for linear differential equations with random parameters following various types of distribution laws are proposed. Mainly, an internal random coefficients method ‘IRCM’ is elaborated for a large number of random parameters. A procedure to build a new polynomial chaos basis and a connection between the one-dimensional and multidimensional polynomials are developed. This allows handling easily random parameters with various laws. A compact matrix formulation is given and the required matrices and scalar products are explicitly presented. For random excitations with an arbitrary number of uncertain variables, the IRCM is couplet to the superposition method leading to successive random differential equations with the same main random operator and right-hand sides depending only on one random parameter. This methodological approach leads to equations with a reduced number of random variables and thus to a large reduction of CPU time and memory required for the numerical solution. The conditional expectation method is also elaborated for reference solutions as well as the Monte-Carlo procedure. The applicability and effectiveness of the developed methods are demonstrated by some numerical examples.


---

## Body

## 1. Introduction

Stochastic and random differential equations constitute a growing field of great scientific interest. There are mainly three categories of random differential equations. The first and the simplest class is one where only the initial conditions are random. The second class is characterized by the presence of random nonhomogeneous or input terms and the third one is the differential equations with random coefficients. To deal with errors and uncertainties, random coefficients have been increasingly used in the last few decades.This paper focuses on the combined second and third classes because this type of equations offers a natural and rational approach for mathematical modeling of many physical phenomena. The last decades have witnessed an enormous effort in the fields of parameters uncertainty and random or stochastic differential processes. This is due to the fact that any physical system contains uncertainties and its real phenomena may be modeled by stochastic differential equations with random or stochastic process coefficients. These equations take into account the approximate knowledge of the numerical values of the physical parameters on which the system depends and have been a matter of intensive investigation.A number of techniques are available for uncertainty sensitivity and propagation such as Monte-Carlo procedure [1, 2], sensitivity analysis methods [3], and polynomial chaos [4, 5] among others. Monte-Carlo (MC) method has been the mainstream uncertainty quantification technique for decades. It is the most used method and is valid for a wide range of problems. However, it is very computationally expensive since it requires a large number of simulations using the full model.An alternative approach is based on the expansion of the response in terms of a series of polynomials that are orthogonal with respect to mean value operations. Polynomial chaos was first introduced by Wiener [6] where Hermite polynomials were used to model stochastic processes with Gaussian random variables. A number of other expansions have been proposed in the literature for representing non-Gaussian process [7, 8]. Recent review papers by Stefanou [9] and by Schuëller and Pradlwarter [10] summarized the assessment of the past and current status of the procedure for stochastic structural analysis.This polynomial representation provides a framework suitable for computational simulation and then widespread in mathematical and numerical analysis of many engineering problems. Various problems have been solved based on this approximation such as solution of stochastic differential equations [11], linear structural dynamics [4, 5], and nonlinear random vibration [12, 13]; soil-structure interaction [14], structural reliability [15], and identification [16, 17]. More recently, Trcala used polynomial chaos for nonlinear diffusion problems of moisture transfer in wood [18]. The accuracy of the PC approximation has been evaluated by Field and Grigoriu [19]. A convergence of the decomposition of the solution into the polynomials chaos is studied by Dvurecenskij et al. [20] and the conditions associated with the distribution function of the random vector appearing in the solution for a convergence toward the solution are given by Ernst et al. [21].The polynomial chaos has been used in many finite elements problems [4]. Accurate discrete modeling of complex industrial structures leads to a large finite element model. To reduce the CPU time, reducing the order of the models is very useful. Component mode synthesis (CMS) is a well-established method for efficiently constructing models to analyze the dynamics of large and complex structures that are often described by separate substructure (or components) models. Sarsri et al. [22] have coupled the CMS methods with the projection chaos polynomials methods in the first and second orders to compute the frequency transfer functions of stochastic structures. This coupling methodological approach has been used by Sarsri and Azrar in time domain [23] as well as a coupling with the perturbation method [24].The polynomial chaos methods are well suited for the random differential equations, RDE, with a very few number of random variables defining their main coefficients. It is well known that if the number of the considered random variables increases, the needed number of unknowns to be determined for solving the random systems increases very rapidly with the degree of the polynomials. Thus, for accurate solution, the CPU time and memory required may be prohibitive. This greatly limits these methods to random differential equations with very few numbers of random parameters.An alternative approach called internal random coefficients method (IRCM) is developed in this paper. A careful presentation is given in the frame of higher-order random differential equations. This method is based on generalized polynomial chaos and the superposition principle. It can be used to solve random differential equations with a large number of random variables and an input right-hand side decomposed in an arbitrary number of random coefficients. The considered random parameters may follow various distribution laws.A procedure to build a new polynomial chaos basis and a connection between the one-dimensional and multidimensional polynomials is established. Different distribution laws can be easily considered. Based on the superposition principle, the random differential equation with an input depending on several random variables is decomposed on a sequence of RDE with the same main random operator and reduced right-hand sides. A series of RDEs with reduced number of random variables have thus to be solved based on the generalized polynomial chaos decomposition. The global system is then solved by a drastic reduction of the CPU time and memory space. For the sake of comparison, the conditional expectation method is developed for the considered random differential equations as well as the Monte-Carlo method. The applicability and effectiveness of the presented methodological approach have been demonstrated by numerically solving various examples.


## 2. Mathematical Formulation

In this work, various methodological approaches are elaborated to solve higher-order initial value problems with linear and nonlinear random variables subjected to a random input right-hand side. For this aim, the following stochastic differential equation is considered:(1)
L
·
U
t
,
ω
=
F
twith deterministic initial conditions, where U is the stochastic process response and L is a linear random operator of order m defined by(2)
L
=
∑
i
=
0
m
C
i
ω
1
d
i
d
t
i
.It is assumed that the random coefficientsC
i depend on the random vector ω
1 which is defined in a probability space (
Ω
1
,
F
1
,
P
1
)
. The input right-hand side, F, is assumed to be dependent on the random vector ω
2 that is defined in the probability space (Ω
2
,
F
2
,
P
2). Explicit expression of F
(
t
) is given later.For numerical solution of (1), a new procedure based on the general polynomial chaos, GPC, expansion procedure is elaborated. Herein, the classical GPC procedure is reviewed in a clear manner.In the present work, the random variables, component of the vector(
ω
1
,
ω
2
), are assumed to be independent but may have general distinct distributions f
j. If f
j are classical distributions, such as normal, uniform, gamma, and beta, the associated known polynomial chaos can be used. Otherwise, the procedure, developed in this paper, will be used to build the needed polynomial basis. Explicit expressions of this basis for general cases are given. In addition, the number of random variables (
n
+
1
) and the differential order m are arbitrary.The concept of internal random coefficients is introduced and combined with the superposition principle and generalized polynomial chaos expansion. For the sake of comparison, conditional expectation and Monte-Carlo methods are also elaborated.


### 2.1. General Polynomial Chaos Formulation

#### 2.1.1. General Formulation

For general purpose, let us consider random vectorsω
1 and ω
2, presented in the following forms:(3)
ω
1
=
ξ
0
,
ξ
1
,
…
,
ξ
p
ω
2
=
ξ
p
+
1
,
ξ
p
+
2
,
…
,
ξ
n
,where ξ
i for i
=
0 to n are random variables defined from the probabilistic field Ω
i
,
F
i
,
P
i to R. The random vectors ω
1 and ω
2 are assumed to be independent and gathered in the vector ω:(4)
ω
=
ω
1
,
ω
2
.We assume that the random vector ω has a distribution function with respect to the Lebesgue measure denoted by f. L
f
2
R
n
+
1 denotes the set of square-integrable functions with respect to the weight measure f:(5)
L
f
2
R
n
+
1
=
H
:
∫
R
n
+
1
H
2
x
f
x
d
x
<
+
∞with the following associated inner product:(6)
H
,
G
=
∫
R
n
+
1
H
x
G
x
f
x
d
x
.Let us note that the distribution f may be Gaussian or non-Gaussian. In the present analysis, various types of distribution functions may be considered.The general polynomial chaos associated with the random vectorω is denoted by ϕ
k
;
k
∈
N. These polynomials coincide with the orthogonal polynomials associated with the inner product defined in (6) and verify(7)
ϕ
k
,
ϕ
j
=
α
k
δ
k
j
,where α
k are given by(8)
α
k
=
ϕ
k
,
ϕ
k
.The solutionU of the main equation (1) is a time stochastic process depending on the random vectorω and decomposed in the polynomial chaos basis ϕ
k
;
k
∈
N:(9)
U
t
,
ω
=
∑
k
≥
0
ϕ
k
ω
U
k
t
.For general purpose, the random coefficientsC
i are assumed to depend linearly and nonlinearly on the random variables ξ
l; 0
≤
l
≤
p and written in the following general form:(10)
C
i
ω
1
=
∑
k
=
0
q
i
∑
α
=
k
c
α
i
∏
l
=
0
p
ξ
l
α
l
,in which α
=
(
α
0
,
…
,
α
p
), α
j
∈
N and α
=
∑
j
=
0
p
α
j.The right-hand side of (1), F
(
t
), is assumed to be a time-dependent random function that depends linearly on ω
2 and expressed by(11)
F
t
,
ω
2
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
g
j
t
,where g
j
(
t
) and a
0
j
,
a
1
j are considered deterministic function and constants.Note that a more general right-hand side excitation can be decomposed in the form (11) using the Karhunen-Loéve expansion [25].Based on a reduced decomposition using the(
N
+
1
) first terms, the stochastic process U
t
,
ω can be approximated by(12)
U
N
t
,
ω
≅
∑
l
=
0
N
ϕ
l
ω
U
l
t
,where ϕ
l
(
ω
) is a multidimensional general polynomial chaos depending on the random vector ω
=
(
ξ
0
,
ξ
1
,
…
,
ξ
n
).The insertion of these expressions in (1) leads to the following m
t
h order random differential equation:(13)
∑
i
=
0
m
∑
l
=
0
N
∑
s
=
0
q
i
∑
α
=
s
c
α
i
∏
d
=
0
p
ξ
d
α
d
ϕ
l
ω
d
i
U
l
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
g
j
t
.Projecting this equation with respect toϕ
k for k
=
0 to N, the following deterministic differential system is then obtained:(14)
∑
i
=
0
m
∑
l
=
0
N
∑
s
=
0
q
i
∑
α
=
s
c
α
i
∏
d
=
0
p
ξ
d
α
d
ϕ
l
,
ϕ
k
d
i
U
l
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
,
ϕ
k
g
j
t
.Note that the first- and mainly the second-order differential equations of the above kind,m
=
1 or 2, have been investigated by many authors when the number n of random variables is too small. When the random variables ξ
j are Gaussian, Hermite-chaos polynomials in ]-
∞
,
+
∞[ are used in [4]. This standard approach is very often used in structural dynamics. Various other works are elaborated when ξ
j are uniform, gamma, or beta and thus Legendre-chaos in [
a
,
b
], Laguerre-chaos in [
0
,
+
∞
[, and Jacobi-chaos in [
a
,
b
] are, respectively, used [8].The expansion on polynomial basis of the vectorω is related to the polynomial basis associated wih each random variable ξ
j in the case of independent variables. This relation is clarified and the correspondence is clearly established herein. The procedure allowing clarifying the inner product used in the differential system (14) is established and an explicit simple procedure is given, in the next subsection.Firstly, this procedure is established in the next paragraphs for independent random variables, based on the relationships between the random variablesξ
j and random vector ω. Secondly, new variables are introduced that are not necessarily independent and the procedure is established for general cases.


#### 2.1.2. Condensed Formulation

In order to formulate the considered problem in a condensed form, the following mathematical developments will be used. As the variables(
ξ
i
,
0
≤
i
≤
n
) are assumed to be pairwise independent, the joint distribution function f is then given by(15)
f
ω
=
∏
k
=
0
n
f
k
ξ
k
,where f
k, k
=
0 to n are functions of marginal distributions associated with each variable ξ
k.The general polynomial chaos associated with each variableξ
k is denoted by ϕ
k
j
,
j
≥
0. These polynomials coincide with the orthogonal polynomials associated with the inner product defined in L
f
k
2
(
R
) with respect to the weight function f
k:(16)
L
f
k
2
R
=
h
:
∫
R
h
2
x
f
k
x
d
x
<
+
∞with the associated inner product given by(17)
h
,
g
k
=
∫
R
h
x
g
x
f
k
x
d
x
.The set of orthogonal polynomials satisfies the orthogonal conditions:(18)
ϕ
k
i
,
ϕ
k
j
k
=
α
i
k
δ
i
j
,where(19)
α
i
k
=
ϕ
k
i
,
ϕ
k
i
k
.In order to make a correspondence between the set of polynomial chaos associated with each variableξ
i and that associated with the random vector ω, the following total order is introduced on the set N
n
+
1 by(20)
∀
i
0
,
…
,
i
n
,
j
0
,
…
,
j
n
∈
N
n
+
1
×
N
n
+
1
,
i
0
,
…
,
i
n
=
j
0
,
…
,
j
n
⟺
i
k
=
j
k
,
∀
k
∈
0
,
n
i
0
,
…
,
i
n
<
j
0
,
…
,
j
n
i
f
∑
k
=
0
n
i
k
<
∑
k
=
0
n
j
k
o
r
i
f
∑
k
=
0
n
i
k
=
∑
k
=
0
n
j
k
∃
l
∈
0
,
n
,
i
0
,
…
,
i
l
-
1
=
j
0
,
…
,
j
l
-
1
,
i
l
>
j
l
.This order allows considering a bijectionψ
n from N
n
+
1 to N, defined by(21)
ψ
n
:
N
n
+
1
⟶
N
i
0
,
…
,
i
n
⟼
ψ
n
i
0
,
…
,
i
n
=
i
.This bijection relates each elementi
0
,
…
,
i
n in N
n
+
1 by its order i defined in (20). Let L be a nonzero integer; the integer N used in the decomposition (12) is taken as(22)
N
=
ψ
n
0,0
,
…
,
0
,
L
.The choice of N allows decomposing the solution in a set of polynomial chaos associated with the vector ω of degree less than or equal to L. So, for all integers k between 0 and N, there is a single k
0
,
…
,
k
n
∈
N
n
+
1 such that(23)
k
=
ψ
n
k
0
,
…
,
k
n
,
(24)
∑
j
=
0
n
k
j
≤
L
.This one-to-one correspondence allows writing the multidimensional polynomial chaos associated with the random vectorω as a function of the one-dimensional polynomial chaos corresponding to each variable ξ
i by(25)
ϕ
k
ω
=
∏
d
=
0
n
ϕ
d
k
d
,
∀
k
∈
0
,
N
,where (k
0
,
…
,
k
n) is the multi-index associated with the integer k, introduced by (23).For all integersi and h between 0 and n the square matrix β
i
h of order L
+
1 is defined by(26)
β
i
h
l
+
1
,
k
+
1
=
ξ
i
h
ϕ
i
l
,
ϕ
i
k
i
∀
l
,
k
∈
0
,
L
2
.Let℘
L
i
=
R
L
ξ
i, the set of polynomials of degree less than or equal to L. Then, the set of the classical polynomial chaos B
C
P
L
i
=
ϕ
i
k
,
k
∈
0
,
L is an orthogonal basis of ℘
L
i for the inner product defined by (17). Let B
c
L
i
=
1
,
ξ
i
,
…
,
ξ
i
L be the canonical basis of ℘
L
i and P
L
i be the passage matrix from the canonical basis B
C
L
i to the chaos basis B
C
P
L
i. This matrix can be obtained in a standard way, using the Gram-Schmidt procedure or a recursive method.Let the vectors(27)
V
C
i
L
=
1
,
…
,
ξ
i
L
T
,
V
C
P
i
L
=
ϕ
i
0
,
…
,
ϕ
i
L
T
.Then, one has(28)
V
C
P
i
L
=
P
L
i
T
V
C
i
L
.Forh
∈
N and i between 0 and n, the moment m
i
h of order h of the random variable ξ
i is defined by(29)
m
i
h
=
E
ξ
i
h
=
∫
-
∞
+
∞
x
h
f
i
x
d
x
.The square matrices of order(
L
+
1
), B
i
L
h, are defined by(30)
B
i
L
h
l
+
1
,
k
+
1
=
m
i
h
+
k
+
l
∀
l
,
k
∈
0
,
L
2
.This gives explicitly matrixβ
i
h by(31)
β
i
h
=
P
L
i
T
B
i
L
h
P
L
i
.Letk and l be two integers between 0 and N. From (23) we have, respectively, unique elements k
0
,
…
,
k
n and l
0
,
…
,
l
n in N
n
+
1 such that k
=
ψ
n
k
0
,
…
,
k
n and l
=
ψ
n
l
0
,
…
,
l
n. Rewriting ϕ
k using expression (26), the following expression is explicitly obtained:(32)
∏
d
=
0
p
ξ
d
α
d
ϕ
l
,
ϕ
k
=
∏
d
=
0
p
β
d
α
d
l
d
+
1
,
k
d
+
1
∏
d
=
p
+
1
n
β
d
0
l
d
+
1
,
k
d
+
1
.The expression in the right-hand side of (14) is thus given by(33)
a
0
j
+
a
1
j
ξ
j
,
ϕ
k
=
a
0
j
β
j
0
1
,
k
j
+
1
+
a
1
j
β
j
1
1
,
k
j
+
1
∏
d
=
0
d
≠
j
n
β
d
0
1
,
k
d
+
1
.For independent random variables, a relationship between the multidimensional and associated one-dimensional generalized polynomials is established. This leads to explicit and closed forms of the used scalar product and needed terms to be numerically computed. Using these relationships resulted in the following deterministic differential system:(34)
∑
i
=
0
m
∑
l
=
0
N
∑
s
=
0
q
i
∑
α
=
s
c
α
i
∏
d
=
0
p
β
d
α
d
l
d
+
1
,
k
d
+
1
∏
d
=
p
+
1
n
β
d
0
l
d
+
1
,
k
d
+
1
d
i
U
l
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
β
j
0
1
,
k
j
+
1
+
a
1
j
β
j
1
1
,
k
j
+
1
∏
d
=
0
d
≠
j
n
β
d
0
1
,
k
d
+
1
g
j
t
.For a compact formulation and using notation (26), the square matrices T
i of order (
N
+
1
) and the time-dependent vector D of dimension (
N
+
1
) are introduced for all integers k, l
∈
[
0
,
N
] by(35)
T
i
l
+
1
,
k
+
1
=
∑
s
=
0
q
i
∑
α
=
s
c
α
i
∏
d
=
0
p
β
d
α
d
l
d
+
1
,
k
d
+
1
∏
d
=
p
+
1
n
β
d
0
l
d
+
1
,
k
d
+
1
,
(36)
D
k
+
1
t
=
∑
j
=
p
+
1
n
a
0
j
β
j
0
1
,
k
j
+
1
+
a
1
j
β
j
1
1
,
k
j
+
1
∏
d
=
0
d
≠
j
n
β
d
0
1
,
k
d
+
1
g
j
t
.Using these notations, differential system (34) is rewritten in the following closed form system:(37)
∑
i
=
0
m
∑
l
=
0
N
T
i
l
+
1
,
k
+
1
d
i
U
l
t
d
t
i
=
D
k
+
1
t
;
k
=
0
,
…
,
N
.Using the presented methodological approach, the truncated solution (12) can be numerically obtained. Its mean and variance are given by(38)
m
N
t
=
U
0
t
,
V
N
t
=
∑
k
=
1
N
α
k
U
k
2
t
.Equation (14) is usually given for Hermite polynomials when random variables are Gaussian. This kind of projection is classically done and many authors follow this procedure.In this paper, the random variablesξ
i may follow various types of laws. The presented generalized formalism allows one to handle multilaws by using the canonical basis and the orthogonalization principle with respect to a scalar product associated with the distribution function of the random vector. Explicit and closed forms of the used scalar products and needed terms to be numerically computed are given.It should to be noted that the major problem of this classical decomposition into polynomials chaos expansion is that the number of unknowns to estimate increases very rapidly when the degree of the polynomial chaos and the number of random parameters increase. More clearly, for n random variables, the number of unknown coefficients in the polynomial chaos of orders less than or equal toL is n
+
L
L. Table 1 presents the numbers of the needed unknown terms for various L and n. This very fast growth of dimensionality is the main limitation of this classical approach.Table 1
Number of needed unknown terms in the polynomial chaos.
n
L
5
10
20
30
5
252
3,003
53,130
32,4632
10
3,003
184,756
30,045,015
847,660,528
20
53,130
30,045,015
137,846,528,820
47,129,212,243,960
30
324,632
847,660,528
47,129,212,243,960
11,826,458,156,4861,424To overcome this drawback, a concept of internal random coefficients is introduced herein. This allows reducing drastically the number of random variables, especially if the initially considered number of random variables is large.


### 2.2. Internal Random Coefficients Method (IRCM)

Let us recall thatm is the derivative order of the random operator L and p is the number of the random variables defining L.Ifp
≤
m, then better use the methodological approach presented in the previous sections.Ifm
<
p, then the coefficients C
i will be considered as the new random variables and called here internal random coefficients. In fact, the following procedure will be of big interest when the number p of random variables ξ
j is too large with respect to the differential order m. A large reduction will result based on the procedures developed herein.Let us consider the worst case:m
<
p. The coefficients C
i
0
≤
i
≤
m and the random vector ω
2 can be gathered in a new random vector:(39)
Y
1
=
C
0
,
C
1
,
…
,
C
m
,
ω
2
,which depends on the initial random vector ω
=
(
ω
1
,
ω
2
). It has to be noted that the coefficients C
i are not independent of each other but independent of ω
2.In this section, it is assumed that the distribution functionf of the random vector ω is a continuous function in R
n
+
1 and has a compact support denoted by K
1.This hypothesis ensures the quadratic convergence of the series formed by polynomials chaos to the solution [20, 21]. Further, we assume that there exists a diffeomorphism h from K
1 to K
2
=
h
K
1, defined by(40)
h
:
K
1
⟶
K
2
ω
⟼
h
ω
.Such that, for alli between 0 and m, one has(41)
h
i
ω
=
C
iand, for all i between p
+
1 and n, one has(42)
h
i
ω
=
ξ
i
.The other components ofh, from m
+
1 to p, are chosen from the components of ω
1 to complete the construction of h.Let us putY
=
h
(
ω
). The random vector Y
, so defined, has a distribution function with respect to the Lebesgue measure, denoted by g. This function is given over the distribution function associated with the random vector ω by (43)
g
y
=
f
h
-
1
y
J
h
-
1
y
,where J
h
-
1
y is the determinant of the Jacobean matrix of h
-
1
y.Letg
1 be the distribution function related to the considered Lebesgue measure of the random vector Y
1 given by(44)
g
1
y
1
=
∫
R
p
-
m
g
y
d
y
m
+
1
⋯
d
y
p
,where y
1
=
(
y
0
,
…
,
y
m
,
y
p
+
1
,
…
,
y
n
) is an observation of the random vector Y
1.For normalized parameters, the following reduced variables are introduced:(45)
X
i
1
=
C
i
-
E
C
i
σ
i
F
o
r
0
≤
i
≤
m
,where σ
i
=
σ
(
C
i
) and(46)
X
1
=
X
11
,
ω
2
;
X
11
=
X
0
1
,
…
,
X
m
1
.Then, for, all integers i between 0 and m, one has(47)
C
i
=
E
C
i
+
σ
i
X
i
1
.Using these new expressions, the main random equation (1) is reduced to the following simplified random differential equation:(48)
∑
i
=
0
m
E
C
i
+
σ
i
X
i
1
d
i
U
t
,
X
1
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
g
j
⁡
t
.It has to be noted here that the number of random variables is reduced and the resulting random differential equation is thus easier to handle.Recall that the coefficients and the random variables, used in the right-hand side of (48), are independent of each other as well as of X
i
1. This is due to the fact that the random vector X
11 depends only on the random vector ω
1. X
11 and ω
2 are then independent.Letg
11 and g
12 be, respectively, the distribution functions of the random vectors X
11 and ω
2 with respect to the associated Lebesgue measures. Taking into account the independence of random variables ξ
i
p
+
1
≤
i
≤
n one has (49)
g
12
x
p
+
1
,
…
,
x
n
=
∏
k
=
p
+
1
n
f
k
x
k
.The decomposition of the solutionU
t
,
X
1 of the random equation (48), according to the multidimensional polynomial chaos basis, is given by (50)
U
t
,
X
1
=
∑
k
≥
0
Ψ
k
1
X
1
u
k
t
.UsingM
t
h first terms of the series (50), the process U
t
,
X
1 can be approximated by(51)
U
M
t
,
X
1
≅
∑
k
=
0
M
Ψ
k
1
X
1
u
k
t
.Inserting these expressions in (48), the following random differential equation resulted:(52)
∑
k
=
0
M
∑
i
=
0
m
E
C
i
+
σ
i
1
X
i
1
Ψ
k
1
X
1
d
i
u
k
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
g
j
t
.Projecting this equation with respect toΨ
l
1 for l
=
0, M leads to(53)
∑
k
=
0
M
∑
i
=
0
m
E
C
i
+
σ
i
1
X
i
1
Ψ
k
1
,
Ψ
l
1
d
i
u
k
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
,
Ψ
l
1
g
j
t
.LetΨ
j
11
,
j
∈
N be polynomial chaos associated with the random vector X
11
.Using bijection (21), there exist two elements in N
n
-
p
+
m
+
1
k
0
,
…
,
k
n
+
m
-
p and l
0
,
…
,
l
n
+
m
-
p, such that k
0
,
…
,
k
n
+
m
-
p
=
ψ
n
+
m
-
p
-
1
k
,
l
0
,
…
,
l
n
+
m
-
p
=
ψ
n
+
m
-
p
-
1
l for all integers k, l and(54)
k
1
=
ψ
m
k
0
,
…
,
k
m
,
l
1
=
ψ
m
l
0
,
…
,
l
m
.The multidimensional polynomial chaosΨ
k
1 associated with the random vector X
1 can be expressed by the polynomial chaos associated with the random vector X
11 and the polynomial chaos associated with the random variables (
ξ
j
,
p
+
1
≤
j
≤
n
) using the notations (23) and (54) by(55)
Ψ
k
1
X
1
=
Ψ
k
1
11
X
11
∏
q
=
p
+
1
n
ϕ
q
k
q
+
m
-
p
ξ
q
.Inserting this relationship into (53), the following random differential equation resulted: (56)
∑
k
=
0
M
∑
i
=
0
m
E
C
i
+
σ
i
1
X
i
1
Ψ
k
1
11
X
11
∏
q
=
p
+
1
n
ϕ
q
k
q
+
m
-
p
ξ
q
,
Ψ
l
1
11
X
11
∏
q
=
p
+
1
n
ϕ
q
l
q
+
m
-
p
ξ
q
d
i
u
k
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
,
Ψ
l
1
11
X
11
∏
q
=
p
+
1
n
ϕ
q
l
q
+
m
-
p
ξ
b
g
j
t
.Based on the independence of the random vectorX
11 and the random variables (
ξ
j
,
p
+
1
≤
j
≤
n
), the inner products in (56) are rearranged as(57)
E
C
i
+
σ
i
1
X
i
1
Ψ
k
1
11
∏
q
=
p
+
1
n
ϕ
q
k
q
+
m
-
p
ξ
q
,
Ψ
l
1
11
∏
q
=
p
+
1
n
ϕ
q
l
q
+
m
-
p
ξ
q
=
m
i
1
Ψ
k
1
11
,
Ψ
l
1
11
11
+
σ
i
1
X
i
1
Ψ
k
1
11
,
Ψ
l
1
11
11
∏
q
=
p
+
1
n
α
k
q
+
m
-
p
q
δ
k
q
+
m
-
p
l
q
+
m
-
p
a
0
j
+
a
1
j
ξ
j
,
Ψ
l
1
11
∏
q
=
p
+
1
n
ϕ
q
l
q
+
m
-
p
ξ
q
=
δ
0
l
1
∏
q
=
p
+
1
q
≠
j
n
δ
0
l
q
+
m
-
p
a
0
j
δ
0
l
j
+
m
-
p
+
a
1
j
δ
j
+
1
+
m
-
p
l
j
+
m
-
p
,where δ
p
q is the Kronecker symbol.LetB
c be the canonical basis of R
X
0
,
…
,
X
m. This basis can be ordered in the following form:(58)
B
c
=
∏
k
=
0
m
X
k
α
k
-
α
k
+
1
:
w
h
e
r
e
α
0
=
q
,
0
≤
α
k
+
1
≤
α
k
,
α
m
+
1
=
0
;
q
∈
N
.The orthogonalization of the basisB
c gives the polynomial chaos Ψ
k
11
,
k
∈
N associated with the random vector X
11. This orthogonalization is performed using a classical method such as the Gram-Schmidt procedure in R
X
0
,
…
,
X
m and the following inner product: (59)
H
,
K
11
=
∫
R
m
+
1
H
x
K
x
g
11
x
d
x
∀
H
,
K
∈
R
X
11
2
.LetB
1 be the tensor defined for all α
∈
N
m
+
1, α
=
(
α
0
,
…
,
α
m
) by(60)
B
1
α
=
E
∏
k
=
0
m
X
k
α
k
=
∫
R
m
+
1
∏
k
=
0
m
x
k
α
k
g
11
x
d
x
.This integral can be computed using the distribution function of the vectorω
1 and the diffeomorphism h by(61)
B
1
α
=
∫
R
p
+
1
∏
k
=
0
m
h
k
ω
1
-
m
k
σ
k
α
k
f
1
ω
1
d
ω
1
.The integerM
, considered in decomposition (51), is chosen such that the considered polynomials have the degree less than or equal to L and expressed by(62)
M
=
ψ
n
+
m
-
p
0,0
,
…
,
0
,
L
∀
L
∈
N
.LetB
c
L be the canonical basis of R
L
X
0
,
…
,
X
m, then(63)
B
c
L
=
∏
k
=
0
m
X
k
α
k
-
α
k
+
1
:
where
α
0
=
q
,
0
≤
α
k
+
1
≤
α
k
,
α
m
+
1
=
0
;
q
∈
N
,
q
≤
L
.Let us denote byP
1 the passage matrix from the canonical basis B
c
L to the polynomial chaos basis associated with the vector X
1 of degree less than or equal to L. Using the following notations,(64)
N
L
=
ψ
m
0,0
,
…
,
0
,
L
,
Q
l
=
∏
i
=
0
m
X
i
l
i
w
h
e
r
e
l
0
,
…
,
l
m
=
ψ
m
-
1
l
for
0
≤
l
≤
N
L
,
V
c
11
=
Q
0
,
…
,
Q
N
L
T
,
V
C
P
11
=
Ψ
0
11
,
…
,
Ψ
N
L
11
T
,one gets(65)
V
C
P
11
=
P
1
T
V
c
11
.The matricesΥ, Τ
i and D, defined for all integers k and l between 0 and N
L, are introduced for 0
≤
i
≤
m, by(66)
Υ
k
+
1
,
l
+
1
=
Ψ
k
11
,
Ψ
l
11
11
=
α
k
11
δ
k
l
(67)
Τ
i
k
+
1
,
l
+
1
=
X
i
1
Ψ
k
11
,
Ψ
l
11
11
(68)
D
k
+
1
,
l
+
1
=
B
1
ψ
m
-
1
k
+
ψ
m
-
1
l
for
0
≤
k
≤
N
L
,
0
≤
l
≤
N
L
.Let(
e
i
)
0
≤
i
≤
m, (
m
+
1
) elements in N
m
+
1 such that(69)
e
0
=
1,0
,
…
,
0
;
e
1
=
0,1
,
…
,
0
;
⋮
e
m
=
0,0
,
…
,
1
.The matricesD
i of order (N
L
+
1) for 0
≤
i
≤
m are defined by (70)
D
i
k
+
1
,
l
+
1
=
B
1
ψ
m
-
1
k
+
ψ
m
-
1
l
+
e
i
for
0
≤
k
≤
N
L
,
0
≤
l
≤
N
L
.This allows defining the matricesΥ and Τ
i by (71)
Υ
=
P
1
D
P
1
T
,
Τ
i
=
P
1
D
i
P
1
T
.The random differential system (52) is thus reduced to the following deterministic one:(72)
∑
k
=
0
M
∑
i
=
0
m
E
C
i
Υ
k
1
+
1
,
l
1
+
1
+
σ
i
1
Τ
i
k
1
+
1
,
l
1
+
1
∏
q
=
p
+
1
n
α
k
q
+
m
-
p
q
δ
k
q
+
m
-
p
l
q
+
m
-
p
d
i
u
k
t
d
t
i
=
∑
j
=
p
+
1
n
δ
0
l
1
∏
q
=
p
+
1
q
≠
j
n
δ
0
l
q
+
m
-
p
a
0
j
δ
0
l
j
+
m
-
p
+
a
1
j
δ
j
+
1
+
m
-
p
l
j
+
m
-
p
g
j
t
,where k
1 and l
1 are integers related to k and l and introduced in (54).Finally, the last differential system is rewritten in the following condensed matrix form:(73)
∑
k
=
0
M
∑
i
=
0
m
H
i
k
1
+
1
,
l
1
+
1
d
i
u
k
t
d
t
i
=
E
l
t
,where (74)
H
i
k
1
+
1
,
l
1
+
1
=
E
C
i
Υ
k
1
+
1
,
l
1
+
1
+
σ
i
1
Τ
i
k
1
+
1
,
l
1
+
1
∏
q
=
p
+
1
n
α
k
q
+
m
-
p
q
δ
k
q
+
m
-
p
l
q
+
m
-
p
.
(75)
E
l
t
=
∑
j
=
p
+
1
n
δ
0
l
1
∏
q
=
p
+
1
q
≠
j
n
δ
0
l
q
+
m
-
p
a
0
j
δ
0
l
j
+
m
-
p
+
a
1
j
δ
j
+
1
+
m
-
p
l
j
+
m
-
p
g
j
tBased on the presented mathematical procedures, the truncated solution (51) can be numerically obtained.An equivalent form of (73) has been introduced in (37) based on original random variables. These two equations give the approximate solution ofU. The coefficients of these two equations are given, respectively, by the inner products in (57) for (73) and (32)-(33) for (37).It should be noted, on one hand, that these coefficients are hard to be determined for (37) but easier for (73). On the other hand, the number of the random variables used in (73) is largely reduced. This reduces the number of unknowns used to determine the solutionU and greatly simplifies the required numerical computation.Note that when uncertainties come from random parameters of the system parameters, they can be efficiently handled by the previous mathematical procedures. When the right-hand side excitation depends on many random variables, the previous internal coefficients procedure can be largely improved by the so-called superposition method.


### 2.3. Superposition Method

Thanks to the linear behavior of the considered system, the superposition principle can be applied to the standard random differential equation (13) and to the resulting equation from the internal random coefficients procedure (52). In both cases, the considered random right-hand side function F
(
t
) depends linearly on (
n
-
p
) random variables ξ
j. The superposition method consists of solving (
n
-
p
) differential equations with elementary right-hand sides given by(76)
L
-
·
U
t
,
ω
=
a
0
j
+
a
1
j
ξ
j
g
j
⁡
t
for
j
=
p
+
1
,
…
,
n
,where L
- may be the left-hand side differential operator of (13) or (52). In fact, this superposition is adding only one random variable ξ
j to the system and the resulting random vectors are (
ω
1
,
ξ
j
) for (13) and (
X
11
,
ξ
j
) for (52) for each j.As the internal random coefficient procedure leads to a reduced number of random variables, the superposition principle is applied to (52). Our focus is then on the resulting random differential equation: (77)
∑
i
=
0
m
E
C
i
+
σ
i
1
X
i
1
d
i
U
t
,
X
1
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
g
j
⁡
t
,with the deterministic initials conditions given by(78)
d
i
U
t
d
t
i
t
=
0
=
u
0
i
for
0
≤
i
≤
m
-
1
.The polynomial chaos associated with the random vectorω
2 is considered and denoted by ϕ
k
2
;
k
∈
N. The solution U is in this case a process that depends on t, on the random vector X
11, and on the random vector ω
2. The solution is then written in the following form: (79)
U
t
,
X
1
=
U
t
,
X
11
,
ω
2
.Since random vectorsX
11 and ω
2 are independent, the factorization by the polynomial chaos associated with the random vector ω
2 is used in the decomposition (50). The solution can be rewritten as follows:(80)
U
t
,
X
11
,
ω
2
=
∑
k
≥
0
U
k
t
,
X
11
ϕ
k
2
ω
2
.The insertion of this decomposition in the differential equation (77) and the projection of the equation and of the initial conditions over each polynomial chaos ϕ
l
2 for l
≥
0 lead to the following simplified system of random differential equations:Fork
=
0,(81a)
∑
i
=
0
m
E
C
i
+
σ
i
1
X
i
1
d
i
U
0
t
,
X
11
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
g
j
⁡
t
d
i
U
0
t
,
X
11
d
t
i
t
=
0
=
u
0
i
for
0
≤
i
≤
m
-
1
.
For1
≤
k
≤
n
-
p,(81b)
∑
i
=
0
m
E
C
i
+
σ
i
1
X
i
1
d
i
U
0
t
,
X
11
d
t
i
=
a
1
p
+
k
g
p
+
k
⁡
t
d
i
U
k
t
,
X
11
d
t
i
t
=
0
=
0
for
0
≤
i
≤
m
-
1
.
Fork
≥
n
-
p
+
1,(81c)
∑
i
=
0
m
E
C
i
+
σ
i
1
X
i
1
d
i
U
k
t
,
X
11
d
t
i
=
0
d
i
U
k
t
,
X
11
d
t
i
t
=
0
=
0
for
0
≤
i
≤
m
-
1
.It has to be noted that the casek
≥
n
-
p
+
1 leads to null initial conditions and null right-hand side. In this case, the solution is then null.The solution given in (80) is then reduced to (82)
U
t
,
X
11
,
ω
2
=
∑
k
=
0
n
-
p
ϕ
k
2
ω
2
U
k
t
,
X
11
.For0
≤
k
≤
n
-
p and considering the decomposition of the solution U
k
(
t
,
X
11
) according to the polynomials chaos associated with the random vector X
11, one gets(83)
U
k
t
,
X
11
=
∑
d
=
0
N
L
u
k
d
t
Ψ
d
11
X
11
.Based on this decomposition, the differential equations (81a)–(81c) are then reduced to the following:Fork
=
0,(84a)
∑
i
=
0
m
∑
d
=
0
N
L
Ψ
d
11
E
C
i
+
σ
i
1
X
i
1
d
i
u
0
d
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
g
j
⁡
t
d
i
u
0
d
t
,
X
11
d
t
i
t
=
0
=
u
0
d
i
for
0
≤
i
≤
m
-
1
.
For1
≤
k
≤
n
-
p,(84b)
∑
i
=
0
m
∑
d
=
0
N
L
Ψ
d
11
E
C
i
+
σ
i
1
X
i
1
d
i
u
k
d
t
d
t
i
=
a
1
p
+
k
g
p
+
k
t
d
i
u
k
d
t
,
X
11
d
t
i
t
=
0
=
0
for
0
≤
i
≤
m
-
1
.Projecting this expression over each polynomial chaosΨ
q
11 for q
=
0 to N
L, one gets (85)
∑
i
=
0
m
∑
d
=
0
N
L
Ψ
d
11
E
C
i
+
σ
i
1
X
i
1
,
Ψ
q
11
d
i
u
0
d
t
d
t
i
=
δ
0
q
∑
j
=
p
+
1
n
a
0
j
g
j
tand, for 1
≤
k
≤
n
-
p,(86)
∑
i
=
0
m
∑
d
=
0
N
L
Ψ
d
11
E
C
i
+
σ
i
1
X
i
1
,
Ψ
q
11
d
i
u
k
d
t
d
t
i
=
δ
k
q
a
1
p
+
k
g
p
+
k
⁡
tUsing the matricesΥ and Τ
i introduced in (71), the last equations become the following.Forq
=
0 to N
L,(87)
∑
i
=
0
m
∑
d
=
0
N
L
E
C
i
Υ
d
1
+
1
,
q
1
+
1
+
σ
i
1
Τ
i
d
1
+
1
,
q
1
+
1
d
i
u
0
d
t
d
t
i
=
δ
0
q
∑
j
=
p
+
1
n
a
0
j
g
j
⁡
tand, for 1
≤
k
≤
n
-
p,(88)
∑
i
=
0
m
∑
d
=
0
N
L
E
C
i
Υ
d
1
+
1
,
q
1
+
1
+
σ
i
1
Τ
i
d
1
+
1
,
q
1
+
1
d
i
u
k
d
t
d
t
i
=
δ
kq
a
1
p
+
k
g
p
+
k
⁡
t
,where the integers d
1 and q
1 are related to d and q by notations (54).These deterministic systems are then rewritten in the following compact matrix form.Forq
=
0 to N
L,(89)
∑
i
=
0
m
∑
d
=
0
N
L
F
i
d
1
+
1
,
q
1
+
1
d
i
u
0
d
t
d
t
i
=
δ
0
q
∑
j
=
p
+
1
n
a
0
j
g
j
⁡
t
d
i
u
0
d
t
,
X
11
d
t
i
t
=
0
=
u
0
d
i
for
0
≤
i
≤
m
-
1and, for 1
≤
k
≤
n
-
p,(90)
∑
i
=
0
m
∑
d
=
0
N
L
F
i
d
1
+
1
,
q
1
+
1
d
i
u
k
d
t
d
t
i
=
a
1
p
+
k
g
p
+
k
⁡
t
δ
k
q
d
i
u
k
d
t
,
X
11
d
t
i
t
=
0
=
0
for
0
≤
i
≤
m
-
1
,where(91)
F
i
=
E
C
i
Υ
+
σ
i
1
Τ
i
.The final truncated solution is then given by(92a)
U
t
=
∑
k
=
0
n
-
p
ϕ
k
2
ω
2
∑
d
=
0
N
L
u
k
d
t
Ψ
d
11
X
11
.Its mean and variance are given:(92b)
m
N
t
=
u
00
t
,
V
N
t
=
∑
k
=
1
n
-
p
∑
d
=
0
N
α
d
11
u
k
d
2
t
+
∑
d
=
1
N
α
d
11
u
0
d
2
t
.This solution is the main approximate solution of the random initial value problem (1) elaborated herein. Following (89) and (90), one has to solve NL(
n
-
p
+
1
) simple initial values differential systems. These systems can be solved by standard methods such as the Runge-Kutta method.The main advantages of this methodological approach are as follows:
(i)
The deterministic differential systems have the same left-hand side that is built only once and in a compact matrix form.
(ii)
Simple right-hand sides are obtained for an arbitrary number of random variables defining the excitation.
(iii)
A large reduction of the number of unknowns based on polynomial chaos can result.
(iv)
The CPU time and memory needed can be largely reduced.
(v)
The numerical solution can be done in parallel manner and particularly for a large number of uncertain parameters.


### 2.4. Conditional Expectation Method

For the sake of comparison, a methodological approach based on the conditional expectation method is also elaborated here to solve (1). Remember that the solution U of (1) is a stochastic process that depends on the random vectors ω
1 and ω
2 and time t. Let ω
0
=
(
ω
10
,
ω
20
) be an observation of the random vector ω
=
(
ω
1
,
ω
2
). The conditional expectation of U such that ω
=
ω
0, noted by E
U
∣
ω
=
ω
0, is the solution of the deterministic equation:(93)
L
0
u
t
,
ω
0
=
F
t
,
ω
20
,where L
0 is a deterministic operator defined by(94)
L
0
=
∑
i
=
0
m
C
i
ω
10
d
i
d
t
i
.Like in the first section, the random vectorω is assumed to have a density function relative to the Lebesgue measure denoted by f. The moments of order p (p
∈
N
∗) of the solution of (1) are then given by the Bayes formulas:(95)
m
p
t
=
∫
R
n
+
1
E
U
∣
ω
=
ω
0
p
f
ω
0
d
ω
0
.The numerical computation of the integral (95) needs the solution of (93) for any observation vector ω
0. The analytical solution can be obtained in general cases. Assume that ω
0
∈
Ω, where Ω is a bounded domain in R
n
+
1.The integral can be approximated by the Gauss-Legendre quadrature method. Let the setsω
-
0
i
,
1
≤
i
≤
q and p
-
i
,
1
≤
i
≤
q be Gauss points and associated Gauss weights. The moments m
p
t can be approximated by (96)
m
p
t
≅
∑
i
=
1
q
p
-
i
E
U
∣
ω
=
ω
-
0
i
p
f
ω
0
i
.The main advantage of this method is that it allows obtaining a general solution that is considered here as a reference solution as well as the obtained one based on the Monte-Carlo method. But, the inconvenience of these two methods is that a very large CPU time is needed for accurate numerical solutions.


## 2.1. General Polynomial Chaos Formulation

### 2.1.1. General Formulation

For general purpose, let us consider random vectorsω
1 and ω
2, presented in the following forms:(3)
ω
1
=
ξ
0
,
ξ
1
,
…
,
ξ
p
ω
2
=
ξ
p
+
1
,
ξ
p
+
2
,
…
,
ξ
n
,where ξ
i for i
=
0 to n are random variables defined from the probabilistic field Ω
i
,
F
i
,
P
i to R. The random vectors ω
1 and ω
2 are assumed to be independent and gathered in the vector ω:(4)
ω
=
ω
1
,
ω
2
.We assume that the random vector ω has a distribution function with respect to the Lebesgue measure denoted by f. L
f
2
R
n
+
1 denotes the set of square-integrable functions with respect to the weight measure f:(5)
L
f
2
R
n
+
1
=
H
:
∫
R
n
+
1
H
2
x
f
x
d
x
<
+
∞with the following associated inner product:(6)
H
,
G
=
∫
R
n
+
1
H
x
G
x
f
x
d
x
.Let us note that the distribution f may be Gaussian or non-Gaussian. In the present analysis, various types of distribution functions may be considered.The general polynomial chaos associated with the random vectorω is denoted by ϕ
k
;
k
∈
N. These polynomials coincide with the orthogonal polynomials associated with the inner product defined in (6) and verify(7)
ϕ
k
,
ϕ
j
=
α
k
δ
k
j
,where α
k are given by(8)
α
k
=
ϕ
k
,
ϕ
k
.The solutionU of the main equation (1) is a time stochastic process depending on the random vectorω and decomposed in the polynomial chaos basis ϕ
k
;
k
∈
N:(9)
U
t
,
ω
=
∑
k
≥
0
ϕ
k
ω
U
k
t
.For general purpose, the random coefficientsC
i are assumed to depend linearly and nonlinearly on the random variables ξ
l; 0
≤
l
≤
p and written in the following general form:(10)
C
i
ω
1
=
∑
k
=
0
q
i
∑
α
=
k
c
α
i
∏
l
=
0
p
ξ
l
α
l
,in which α
=
(
α
0
,
…
,
α
p
), α
j
∈
N and α
=
∑
j
=
0
p
α
j.The right-hand side of (1), F
(
t
), is assumed to be a time-dependent random function that depends linearly on ω
2 and expressed by(11)
F
t
,
ω
2
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
g
j
t
,where g
j
(
t
) and a
0
j
,
a
1
j are considered deterministic function and constants.Note that a more general right-hand side excitation can be decomposed in the form (11) using the Karhunen-Loéve expansion [25].Based on a reduced decomposition using the(
N
+
1
) first terms, the stochastic process U
t
,
ω can be approximated by(12)
U
N
t
,
ω
≅
∑
l
=
0
N
ϕ
l
ω
U
l
t
,where ϕ
l
(
ω
) is a multidimensional general polynomial chaos depending on the random vector ω
=
(
ξ
0
,
ξ
1
,
…
,
ξ
n
).The insertion of these expressions in (1) leads to the following m
t
h order random differential equation:(13)
∑
i
=
0
m
∑
l
=
0
N
∑
s
=
0
q
i
∑
α
=
s
c
α
i
∏
d
=
0
p
ξ
d
α
d
ϕ
l
ω
d
i
U
l
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
g
j
t
.Projecting this equation with respect toϕ
k for k
=
0 to N, the following deterministic differential system is then obtained:(14)
∑
i
=
0
m
∑
l
=
0
N
∑
s
=
0
q
i
∑
α
=
s
c
α
i
∏
d
=
0
p
ξ
d
α
d
ϕ
l
,
ϕ
k
d
i
U
l
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
,
ϕ
k
g
j
t
.Note that the first- and mainly the second-order differential equations of the above kind,m
=
1 or 2, have been investigated by many authors when the number n of random variables is too small. When the random variables ξ
j are Gaussian, Hermite-chaos polynomials in ]-
∞
,
+
∞[ are used in [4]. This standard approach is very often used in structural dynamics. Various other works are elaborated when ξ
j are uniform, gamma, or beta and thus Legendre-chaos in [
a
,
b
], Laguerre-chaos in [
0
,
+
∞
[, and Jacobi-chaos in [
a
,
b
] are, respectively, used [8].The expansion on polynomial basis of the vectorω is related to the polynomial basis associated wih each random variable ξ
j in the case of independent variables. This relation is clarified and the correspondence is clearly established herein. The procedure allowing clarifying the inner product used in the differential system (14) is established and an explicit simple procedure is given, in the next subsection.Firstly, this procedure is established in the next paragraphs for independent random variables, based on the relationships between the random variablesξ
j and random vector ω. Secondly, new variables are introduced that are not necessarily independent and the procedure is established for general cases.


### 2.1.2. Condensed Formulation

In order to formulate the considered problem in a condensed form, the following mathematical developments will be used. As the variables(
ξ
i
,
0
≤
i
≤
n
) are assumed to be pairwise independent, the joint distribution function f is then given by(15)
f
ω
=
∏
k
=
0
n
f
k
ξ
k
,where f
k, k
=
0 to n are functions of marginal distributions associated with each variable ξ
k.The general polynomial chaos associated with each variableξ
k is denoted by ϕ
k
j
,
j
≥
0. These polynomials coincide with the orthogonal polynomials associated with the inner product defined in L
f
k
2
(
R
) with respect to the weight function f
k:(16)
L
f
k
2
R
=
h
:
∫
R
h
2
x
f
k
x
d
x
<
+
∞with the associated inner product given by(17)
h
,
g
k
=
∫
R
h
x
g
x
f
k
x
d
x
.The set of orthogonal polynomials satisfies the orthogonal conditions:(18)
ϕ
k
i
,
ϕ
k
j
k
=
α
i
k
δ
i
j
,where(19)
α
i
k
=
ϕ
k
i
,
ϕ
k
i
k
.In order to make a correspondence between the set of polynomial chaos associated with each variableξ
i and that associated with the random vector ω, the following total order is introduced on the set N
n
+
1 by(20)
∀
i
0
,
…
,
i
n
,
j
0
,
…
,
j
n
∈
N
n
+
1
×
N
n
+
1
,
i
0
,
…
,
i
n
=
j
0
,
…
,
j
n
⟺
i
k
=
j
k
,
∀
k
∈
0
,
n
i
0
,
…
,
i
n
<
j
0
,
…
,
j
n
i
f
∑
k
=
0
n
i
k
<
∑
k
=
0
n
j
k
o
r
i
f
∑
k
=
0
n
i
k
=
∑
k
=
0
n
j
k
∃
l
∈
0
,
n
,
i
0
,
…
,
i
l
-
1
=
j
0
,
…
,
j
l
-
1
,
i
l
>
j
l
.This order allows considering a bijectionψ
n from N
n
+
1 to N, defined by(21)
ψ
n
:
N
n
+
1
⟶
N
i
0
,
…
,
i
n
⟼
ψ
n
i
0
,
…
,
i
n
=
i
.This bijection relates each elementi
0
,
…
,
i
n in N
n
+
1 by its order i defined in (20). Let L be a nonzero integer; the integer N used in the decomposition (12) is taken as(22)
N
=
ψ
n
0,0
,
…
,
0
,
L
.The choice of N allows decomposing the solution in a set of polynomial chaos associated with the vector ω of degree less than or equal to L. So, for all integers k between 0 and N, there is a single k
0
,
…
,
k
n
∈
N
n
+
1 such that(23)
k
=
ψ
n
k
0
,
…
,
k
n
,
(24)
∑
j
=
0
n
k
j
≤
L
.This one-to-one correspondence allows writing the multidimensional polynomial chaos associated with the random vectorω as a function of the one-dimensional polynomial chaos corresponding to each variable ξ
i by(25)
ϕ
k
ω
=
∏
d
=
0
n
ϕ
d
k
d
,
∀
k
∈
0
,
N
,where (k
0
,
…
,
k
n) is the multi-index associated with the integer k, introduced by (23).For all integersi and h between 0 and n the square matrix β
i
h of order L
+
1 is defined by(26)
β
i
h
l
+
1
,
k
+
1
=
ξ
i
h
ϕ
i
l
,
ϕ
i
k
i
∀
l
,
k
∈
0
,
L
2
.Let℘
L
i
=
R
L
ξ
i, the set of polynomials of degree less than or equal to L. Then, the set of the classical polynomial chaos B
C
P
L
i
=
ϕ
i
k
,
k
∈
0
,
L is an orthogonal basis of ℘
L
i for the inner product defined by (17). Let B
c
L
i
=
1
,
ξ
i
,
…
,
ξ
i
L be the canonical basis of ℘
L
i and P
L
i be the passage matrix from the canonical basis B
C
L
i to the chaos basis B
C
P
L
i. This matrix can be obtained in a standard way, using the Gram-Schmidt procedure or a recursive method.Let the vectors(27)
V
C
i
L
=
1
,
…
,
ξ
i
L
T
,
V
C
P
i
L
=
ϕ
i
0
,
…
,
ϕ
i
L
T
.Then, one has(28)
V
C
P
i
L
=
P
L
i
T
V
C
i
L
.Forh
∈
N and i between 0 and n, the moment m
i
h of order h of the random variable ξ
i is defined by(29)
m
i
h
=
E
ξ
i
h
=
∫
-
∞
+
∞
x
h
f
i
x
d
x
.The square matrices of order(
L
+
1
), B
i
L
h, are defined by(30)
B
i
L
h
l
+
1
,
k
+
1
=
m
i
h
+
k
+
l
∀
l
,
k
∈
0
,
L
2
.This gives explicitly matrixβ
i
h by(31)
β
i
h
=
P
L
i
T
B
i
L
h
P
L
i
.Letk and l be two integers between 0 and N. From (23) we have, respectively, unique elements k
0
,
…
,
k
n and l
0
,
…
,
l
n in N
n
+
1 such that k
=
ψ
n
k
0
,
…
,
k
n and l
=
ψ
n
l
0
,
…
,
l
n. Rewriting ϕ
k using expression (26), the following expression is explicitly obtained:(32)
∏
d
=
0
p
ξ
d
α
d
ϕ
l
,
ϕ
k
=
∏
d
=
0
p
β
d
α
d
l
d
+
1
,
k
d
+
1
∏
d
=
p
+
1
n
β
d
0
l
d
+
1
,
k
d
+
1
.The expression in the right-hand side of (14) is thus given by(33)
a
0
j
+
a
1
j
ξ
j
,
ϕ
k
=
a
0
j
β
j
0
1
,
k
j
+
1
+
a
1
j
β
j
1
1
,
k
j
+
1
∏
d
=
0
d
≠
j
n
β
d
0
1
,
k
d
+
1
.For independent random variables, a relationship between the multidimensional and associated one-dimensional generalized polynomials is established. This leads to explicit and closed forms of the used scalar product and needed terms to be numerically computed. Using these relationships resulted in the following deterministic differential system:(34)
∑
i
=
0
m
∑
l
=
0
N
∑
s
=
0
q
i
∑
α
=
s
c
α
i
∏
d
=
0
p
β
d
α
d
l
d
+
1
,
k
d
+
1
∏
d
=
p
+
1
n
β
d
0
l
d
+
1
,
k
d
+
1
d
i
U
l
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
β
j
0
1
,
k
j
+
1
+
a
1
j
β
j
1
1
,
k
j
+
1
∏
d
=
0
d
≠
j
n
β
d
0
1
,
k
d
+
1
g
j
t
.For a compact formulation and using notation (26), the square matrices T
i of order (
N
+
1
) and the time-dependent vector D of dimension (
N
+
1
) are introduced for all integers k, l
∈
[
0
,
N
] by(35)
T
i
l
+
1
,
k
+
1
=
∑
s
=
0
q
i
∑
α
=
s
c
α
i
∏
d
=
0
p
β
d
α
d
l
d
+
1
,
k
d
+
1
∏
d
=
p
+
1
n
β
d
0
l
d
+
1
,
k
d
+
1
,
(36)
D
k
+
1
t
=
∑
j
=
p
+
1
n
a
0
j
β
j
0
1
,
k
j
+
1
+
a
1
j
β
j
1
1
,
k
j
+
1
∏
d
=
0
d
≠
j
n
β
d
0
1
,
k
d
+
1
g
j
t
.Using these notations, differential system (34) is rewritten in the following closed form system:(37)
∑
i
=
0
m
∑
l
=
0
N
T
i
l
+
1
,
k
+
1
d
i
U
l
t
d
t
i
=
D
k
+
1
t
;
k
=
0
,
…
,
N
.Using the presented methodological approach, the truncated solution (12) can be numerically obtained. Its mean and variance are given by(38)
m
N
t
=
U
0
t
,
V
N
t
=
∑
k
=
1
N
α
k
U
k
2
t
.Equation (14) is usually given for Hermite polynomials when random variables are Gaussian. This kind of projection is classically done and many authors follow this procedure.In this paper, the random variablesξ
i may follow various types of laws. The presented generalized formalism allows one to handle multilaws by using the canonical basis and the orthogonalization principle with respect to a scalar product associated with the distribution function of the random vector. Explicit and closed forms of the used scalar products and needed terms to be numerically computed are given.It should to be noted that the major problem of this classical decomposition into polynomials chaos expansion is that the number of unknowns to estimate increases very rapidly when the degree of the polynomial chaos and the number of random parameters increase. More clearly, for n random variables, the number of unknown coefficients in the polynomial chaos of orders less than or equal toL is n
+
L
L. Table 1 presents the numbers of the needed unknown terms for various L and n. This very fast growth of dimensionality is the main limitation of this classical approach.Table 1
Number of needed unknown terms in the polynomial chaos.
n
L
5
10
20
30
5
252
3,003
53,130
32,4632
10
3,003
184,756
30,045,015
847,660,528
20
53,130
30,045,015
137,846,528,820
47,129,212,243,960
30
324,632
847,660,528
47,129,212,243,960
11,826,458,156,4861,424To overcome this drawback, a concept of internal random coefficients is introduced herein. This allows reducing drastically the number of random variables, especially if the initially considered number of random variables is large.


## 2.1.1. General Formulation

For general purpose, let us consider random vectorsω
1 and ω
2, presented in the following forms:(3)
ω
1
=
ξ
0
,
ξ
1
,
…
,
ξ
p
ω
2
=
ξ
p
+
1
,
ξ
p
+
2
,
…
,
ξ
n
,where ξ
i for i
=
0 to n are random variables defined from the probabilistic field Ω
i
,
F
i
,
P
i to R. The random vectors ω
1 and ω
2 are assumed to be independent and gathered in the vector ω:(4)
ω
=
ω
1
,
ω
2
.We assume that the random vector ω has a distribution function with respect to the Lebesgue measure denoted by f. L
f
2
R
n
+
1 denotes the set of square-integrable functions with respect to the weight measure f:(5)
L
f
2
R
n
+
1
=
H
:
∫
R
n
+
1
H
2
x
f
x
d
x
<
+
∞with the following associated inner product:(6)
H
,
G
=
∫
R
n
+
1
H
x
G
x
f
x
d
x
.Let us note that the distribution f may be Gaussian or non-Gaussian. In the present analysis, various types of distribution functions may be considered.The general polynomial chaos associated with the random vectorω is denoted by ϕ
k
;
k
∈
N. These polynomials coincide with the orthogonal polynomials associated with the inner product defined in (6) and verify(7)
ϕ
k
,
ϕ
j
=
α
k
δ
k
j
,where α
k are given by(8)
α
k
=
ϕ
k
,
ϕ
k
.The solutionU of the main equation (1) is a time stochastic process depending on the random vectorω and decomposed in the polynomial chaos basis ϕ
k
;
k
∈
N:(9)
U
t
,
ω
=
∑
k
≥
0
ϕ
k
ω
U
k
t
.For general purpose, the random coefficientsC
i are assumed to depend linearly and nonlinearly on the random variables ξ
l; 0
≤
l
≤
p and written in the following general form:(10)
C
i
ω
1
=
∑
k
=
0
q
i
∑
α
=
k
c
α
i
∏
l
=
0
p
ξ
l
α
l
,in which α
=
(
α
0
,
…
,
α
p
), α
j
∈
N and α
=
∑
j
=
0
p
α
j.The right-hand side of (1), F
(
t
), is assumed to be a time-dependent random function that depends linearly on ω
2 and expressed by(11)
F
t
,
ω
2
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
g
j
t
,where g
j
(
t
) and a
0
j
,
a
1
j are considered deterministic function and constants.Note that a more general right-hand side excitation can be decomposed in the form (11) using the Karhunen-Loéve expansion [25].Based on a reduced decomposition using the(
N
+
1
) first terms, the stochastic process U
t
,
ω can be approximated by(12)
U
N
t
,
ω
≅
∑
l
=
0
N
ϕ
l
ω
U
l
t
,where ϕ
l
(
ω
) is a multidimensional general polynomial chaos depending on the random vector ω
=
(
ξ
0
,
ξ
1
,
…
,
ξ
n
).The insertion of these expressions in (1) leads to the following m
t
h order random differential equation:(13)
∑
i
=
0
m
∑
l
=
0
N
∑
s
=
0
q
i
∑
α
=
s
c
α
i
∏
d
=
0
p
ξ
d
α
d
ϕ
l
ω
d
i
U
l
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
g
j
t
.Projecting this equation with respect toϕ
k for k
=
0 to N, the following deterministic differential system is then obtained:(14)
∑
i
=
0
m
∑
l
=
0
N
∑
s
=
0
q
i
∑
α
=
s
c
α
i
∏
d
=
0
p
ξ
d
α
d
ϕ
l
,
ϕ
k
d
i
U
l
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
,
ϕ
k
g
j
t
.Note that the first- and mainly the second-order differential equations of the above kind,m
=
1 or 2, have been investigated by many authors when the number n of random variables is too small. When the random variables ξ
j are Gaussian, Hermite-chaos polynomials in ]-
∞
,
+
∞[ are used in [4]. This standard approach is very often used in structural dynamics. Various other works are elaborated when ξ
j are uniform, gamma, or beta and thus Legendre-chaos in [
a
,
b
], Laguerre-chaos in [
0
,
+
∞
[, and Jacobi-chaos in [
a
,
b
] are, respectively, used [8].The expansion on polynomial basis of the vectorω is related to the polynomial basis associated wih each random variable ξ
j in the case of independent variables. This relation is clarified and the correspondence is clearly established herein. The procedure allowing clarifying the inner product used in the differential system (14) is established and an explicit simple procedure is given, in the next subsection.Firstly, this procedure is established in the next paragraphs for independent random variables, based on the relationships between the random variablesξ
j and random vector ω. Secondly, new variables are introduced that are not necessarily independent and the procedure is established for general cases.


## 2.1.2. Condensed Formulation

In order to formulate the considered problem in a condensed form, the following mathematical developments will be used. As the variables(
ξ
i
,
0
≤
i
≤
n
) are assumed to be pairwise independent, the joint distribution function f is then given by(15)
f
ω
=
∏
k
=
0
n
f
k
ξ
k
,where f
k, k
=
0 to n are functions of marginal distributions associated with each variable ξ
k.The general polynomial chaos associated with each variableξ
k is denoted by ϕ
k
j
,
j
≥
0. These polynomials coincide with the orthogonal polynomials associated with the inner product defined in L
f
k
2
(
R
) with respect to the weight function f
k:(16)
L
f
k
2
R
=
h
:
∫
R
h
2
x
f
k
x
d
x
<
+
∞with the associated inner product given by(17)
h
,
g
k
=
∫
R
h
x
g
x
f
k
x
d
x
.The set of orthogonal polynomials satisfies the orthogonal conditions:(18)
ϕ
k
i
,
ϕ
k
j
k
=
α
i
k
δ
i
j
,where(19)
α
i
k
=
ϕ
k
i
,
ϕ
k
i
k
.In order to make a correspondence between the set of polynomial chaos associated with each variableξ
i and that associated with the random vector ω, the following total order is introduced on the set N
n
+
1 by(20)
∀
i
0
,
…
,
i
n
,
j
0
,
…
,
j
n
∈
N
n
+
1
×
N
n
+
1
,
i
0
,
…
,
i
n
=
j
0
,
…
,
j
n
⟺
i
k
=
j
k
,
∀
k
∈
0
,
n
i
0
,
…
,
i
n
<
j
0
,
…
,
j
n
i
f
∑
k
=
0
n
i
k
<
∑
k
=
0
n
j
k
o
r
i
f
∑
k
=
0
n
i
k
=
∑
k
=
0
n
j
k
∃
l
∈
0
,
n
,
i
0
,
…
,
i
l
-
1
=
j
0
,
…
,
j
l
-
1
,
i
l
>
j
l
.This order allows considering a bijectionψ
n from N
n
+
1 to N, defined by(21)
ψ
n
:
N
n
+
1
⟶
N
i
0
,
…
,
i
n
⟼
ψ
n
i
0
,
…
,
i
n
=
i
.This bijection relates each elementi
0
,
…
,
i
n in N
n
+
1 by its order i defined in (20). Let L be a nonzero integer; the integer N used in the decomposition (12) is taken as(22)
N
=
ψ
n
0,0
,
…
,
0
,
L
.The choice of N allows decomposing the solution in a set of polynomial chaos associated with the vector ω of degree less than or equal to L. So, for all integers k between 0 and N, there is a single k
0
,
…
,
k
n
∈
N
n
+
1 such that(23)
k
=
ψ
n
k
0
,
…
,
k
n
,
(24)
∑
j
=
0
n
k
j
≤
L
.This one-to-one correspondence allows writing the multidimensional polynomial chaos associated with the random vectorω as a function of the one-dimensional polynomial chaos corresponding to each variable ξ
i by(25)
ϕ
k
ω
=
∏
d
=
0
n
ϕ
d
k
d
,
∀
k
∈
0
,
N
,where (k
0
,
…
,
k
n) is the multi-index associated with the integer k, introduced by (23).For all integersi and h between 0 and n the square matrix β
i
h of order L
+
1 is defined by(26)
β
i
h
l
+
1
,
k
+
1
=
ξ
i
h
ϕ
i
l
,
ϕ
i
k
i
∀
l
,
k
∈
0
,
L
2
.Let℘
L
i
=
R
L
ξ
i, the set of polynomials of degree less than or equal to L. Then, the set of the classical polynomial chaos B
C
P
L
i
=
ϕ
i
k
,
k
∈
0
,
L is an orthogonal basis of ℘
L
i for the inner product defined by (17). Let B
c
L
i
=
1
,
ξ
i
,
…
,
ξ
i
L be the canonical basis of ℘
L
i and P
L
i be the passage matrix from the canonical basis B
C
L
i to the chaos basis B
C
P
L
i. This matrix can be obtained in a standard way, using the Gram-Schmidt procedure or a recursive method.Let the vectors(27)
V
C
i
L
=
1
,
…
,
ξ
i
L
T
,
V
C
P
i
L
=
ϕ
i
0
,
…
,
ϕ
i
L
T
.Then, one has(28)
V
C
P
i
L
=
P
L
i
T
V
C
i
L
.Forh
∈
N and i between 0 and n, the moment m
i
h of order h of the random variable ξ
i is defined by(29)
m
i
h
=
E
ξ
i
h
=
∫
-
∞
+
∞
x
h
f
i
x
d
x
.The square matrices of order(
L
+
1
), B
i
L
h, are defined by(30)
B
i
L
h
l
+
1
,
k
+
1
=
m
i
h
+
k
+
l
∀
l
,
k
∈
0
,
L
2
.This gives explicitly matrixβ
i
h by(31)
β
i
h
=
P
L
i
T
B
i
L
h
P
L
i
.Letk and l be two integers between 0 and N. From (23) we have, respectively, unique elements k
0
,
…
,
k
n and l
0
,
…
,
l
n in N
n
+
1 such that k
=
ψ
n
k
0
,
…
,
k
n and l
=
ψ
n
l
0
,
…
,
l
n. Rewriting ϕ
k using expression (26), the following expression is explicitly obtained:(32)
∏
d
=
0
p
ξ
d
α
d
ϕ
l
,
ϕ
k
=
∏
d
=
0
p
β
d
α
d
l
d
+
1
,
k
d
+
1
∏
d
=
p
+
1
n
β
d
0
l
d
+
1
,
k
d
+
1
.The expression in the right-hand side of (14) is thus given by(33)
a
0
j
+
a
1
j
ξ
j
,
ϕ
k
=
a
0
j
β
j
0
1
,
k
j
+
1
+
a
1
j
β
j
1
1
,
k
j
+
1
∏
d
=
0
d
≠
j
n
β
d
0
1
,
k
d
+
1
.For independent random variables, a relationship between the multidimensional and associated one-dimensional generalized polynomials is established. This leads to explicit and closed forms of the used scalar product and needed terms to be numerically computed. Using these relationships resulted in the following deterministic differential system:(34)
∑
i
=
0
m
∑
l
=
0
N
∑
s
=
0
q
i
∑
α
=
s
c
α
i
∏
d
=
0
p
β
d
α
d
l
d
+
1
,
k
d
+
1
∏
d
=
p
+
1
n
β
d
0
l
d
+
1
,
k
d
+
1
d
i
U
l
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
β
j
0
1
,
k
j
+
1
+
a
1
j
β
j
1
1
,
k
j
+
1
∏
d
=
0
d
≠
j
n
β
d
0
1
,
k
d
+
1
g
j
t
.For a compact formulation and using notation (26), the square matrices T
i of order (
N
+
1
) and the time-dependent vector D of dimension (
N
+
1
) are introduced for all integers k, l
∈
[
0
,
N
] by(35)
T
i
l
+
1
,
k
+
1
=
∑
s
=
0
q
i
∑
α
=
s
c
α
i
∏
d
=
0
p
β
d
α
d
l
d
+
1
,
k
d
+
1
∏
d
=
p
+
1
n
β
d
0
l
d
+
1
,
k
d
+
1
,
(36)
D
k
+
1
t
=
∑
j
=
p
+
1
n
a
0
j
β
j
0
1
,
k
j
+
1
+
a
1
j
β
j
1
1
,
k
j
+
1
∏
d
=
0
d
≠
j
n
β
d
0
1
,
k
d
+
1
g
j
t
.Using these notations, differential system (34) is rewritten in the following closed form system:(37)
∑
i
=
0
m
∑
l
=
0
N
T
i
l
+
1
,
k
+
1
d
i
U
l
t
d
t
i
=
D
k
+
1
t
;
k
=
0
,
…
,
N
.Using the presented methodological approach, the truncated solution (12) can be numerically obtained. Its mean and variance are given by(38)
m
N
t
=
U
0
t
,
V
N
t
=
∑
k
=
1
N
α
k
U
k
2
t
.Equation (14) is usually given for Hermite polynomials when random variables are Gaussian. This kind of projection is classically done and many authors follow this procedure.In this paper, the random variablesξ
i may follow various types of laws. The presented generalized formalism allows one to handle multilaws by using the canonical basis and the orthogonalization principle with respect to a scalar product associated with the distribution function of the random vector. Explicit and closed forms of the used scalar products and needed terms to be numerically computed are given.It should to be noted that the major problem of this classical decomposition into polynomials chaos expansion is that the number of unknowns to estimate increases very rapidly when the degree of the polynomial chaos and the number of random parameters increase. More clearly, for n random variables, the number of unknown coefficients in the polynomial chaos of orders less than or equal toL is n
+
L
L. Table 1 presents the numbers of the needed unknown terms for various L and n. This very fast growth of dimensionality is the main limitation of this classical approach.Table 1
Number of needed unknown terms in the polynomial chaos.
n
L
5
10
20
30
5
252
3,003
53,130
32,4632
10
3,003
184,756
30,045,015
847,660,528
20
53,130
30,045,015
137,846,528,820
47,129,212,243,960
30
324,632
847,660,528
47,129,212,243,960
11,826,458,156,4861,424To overcome this drawback, a concept of internal random coefficients is introduced herein. This allows reducing drastically the number of random variables, especially if the initially considered number of random variables is large.


## 2.2. Internal Random Coefficients Method (IRCM)

Let us recall thatm is the derivative order of the random operator L and p is the number of the random variables defining L.Ifp
≤
m, then better use the methodological approach presented in the previous sections.Ifm
<
p, then the coefficients C
i will be considered as the new random variables and called here internal random coefficients. In fact, the following procedure will be of big interest when the number p of random variables ξ
j is too large with respect to the differential order m. A large reduction will result based on the procedures developed herein.Let us consider the worst case:m
<
p. The coefficients C
i
0
≤
i
≤
m and the random vector ω
2 can be gathered in a new random vector:(39)
Y
1
=
C
0
,
C
1
,
…
,
C
m
,
ω
2
,which depends on the initial random vector ω
=
(
ω
1
,
ω
2
). It has to be noted that the coefficients C
i are not independent of each other but independent of ω
2.In this section, it is assumed that the distribution functionf of the random vector ω is a continuous function in R
n
+
1 and has a compact support denoted by K
1.This hypothesis ensures the quadratic convergence of the series formed by polynomials chaos to the solution [20, 21]. Further, we assume that there exists a diffeomorphism h from K
1 to K
2
=
h
K
1, defined by(40)
h
:
K
1
⟶
K
2
ω
⟼
h
ω
.Such that, for alli between 0 and m, one has(41)
h
i
ω
=
C
iand, for all i between p
+
1 and n, one has(42)
h
i
ω
=
ξ
i
.The other components ofh, from m
+
1 to p, are chosen from the components of ω
1 to complete the construction of h.Let us putY
=
h
(
ω
). The random vector Y
, so defined, has a distribution function with respect to the Lebesgue measure, denoted by g. This function is given over the distribution function associated with the random vector ω by (43)
g
y
=
f
h
-
1
y
J
h
-
1
y
,where J
h
-
1
y is the determinant of the Jacobean matrix of h
-
1
y.Letg
1 be the distribution function related to the considered Lebesgue measure of the random vector Y
1 given by(44)
g
1
y
1
=
∫
R
p
-
m
g
y
d
y
m
+
1
⋯
d
y
p
,where y
1
=
(
y
0
,
…
,
y
m
,
y
p
+
1
,
…
,
y
n
) is an observation of the random vector Y
1.For normalized parameters, the following reduced variables are introduced:(45)
X
i
1
=
C
i
-
E
C
i
σ
i
F
o
r
0
≤
i
≤
m
,where σ
i
=
σ
(
C
i
) and(46)
X
1
=
X
11
,
ω
2
;
X
11
=
X
0
1
,
…
,
X
m
1
.Then, for, all integers i between 0 and m, one has(47)
C
i
=
E
C
i
+
σ
i
X
i
1
.Using these new expressions, the main random equation (1) is reduced to the following simplified random differential equation:(48)
∑
i
=
0
m
E
C
i
+
σ
i
X
i
1
d
i
U
t
,
X
1
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
g
j
⁡
t
.It has to be noted here that the number of random variables is reduced and the resulting random differential equation is thus easier to handle.Recall that the coefficients and the random variables, used in the right-hand side of (48), are independent of each other as well as of X
i
1. This is due to the fact that the random vector X
11 depends only on the random vector ω
1. X
11 and ω
2 are then independent.Letg
11 and g
12 be, respectively, the distribution functions of the random vectors X
11 and ω
2 with respect to the associated Lebesgue measures. Taking into account the independence of random variables ξ
i
p
+
1
≤
i
≤
n one has (49)
g
12
x
p
+
1
,
…
,
x
n
=
∏
k
=
p
+
1
n
f
k
x
k
.The decomposition of the solutionU
t
,
X
1 of the random equation (48), according to the multidimensional polynomial chaos basis, is given by (50)
U
t
,
X
1
=
∑
k
≥
0
Ψ
k
1
X
1
u
k
t
.UsingM
t
h first terms of the series (50), the process U
t
,
X
1 can be approximated by(51)
U
M
t
,
X
1
≅
∑
k
=
0
M
Ψ
k
1
X
1
u
k
t
.Inserting these expressions in (48), the following random differential equation resulted:(52)
∑
k
=
0
M
∑
i
=
0
m
E
C
i
+
σ
i
1
X
i
1
Ψ
k
1
X
1
d
i
u
k
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
g
j
t
.Projecting this equation with respect toΨ
l
1 for l
=
0, M leads to(53)
∑
k
=
0
M
∑
i
=
0
m
E
C
i
+
σ
i
1
X
i
1
Ψ
k
1
,
Ψ
l
1
d
i
u
k
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
,
Ψ
l
1
g
j
t
.LetΨ
j
11
,
j
∈
N be polynomial chaos associated with the random vector X
11
.Using bijection (21), there exist two elements in N
n
-
p
+
m
+
1
k
0
,
…
,
k
n
+
m
-
p and l
0
,
…
,
l
n
+
m
-
p, such that k
0
,
…
,
k
n
+
m
-
p
=
ψ
n
+
m
-
p
-
1
k
,
l
0
,
…
,
l
n
+
m
-
p
=
ψ
n
+
m
-
p
-
1
l for all integers k, l and(54)
k
1
=
ψ
m
k
0
,
…
,
k
m
,
l
1
=
ψ
m
l
0
,
…
,
l
m
.The multidimensional polynomial chaosΨ
k
1 associated with the random vector X
1 can be expressed by the polynomial chaos associated with the random vector X
11 and the polynomial chaos associated with the random variables (
ξ
j
,
p
+
1
≤
j
≤
n
) using the notations (23) and (54) by(55)
Ψ
k
1
X
1
=
Ψ
k
1
11
X
11
∏
q
=
p
+
1
n
ϕ
q
k
q
+
m
-
p
ξ
q
.Inserting this relationship into (53), the following random differential equation resulted: (56)
∑
k
=
0
M
∑
i
=
0
m
E
C
i
+
σ
i
1
X
i
1
Ψ
k
1
11
X
11
∏
q
=
p
+
1
n
ϕ
q
k
q
+
m
-
p
ξ
q
,
Ψ
l
1
11
X
11
∏
q
=
p
+
1
n
ϕ
q
l
q
+
m
-
p
ξ
q
d
i
u
k
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
,
Ψ
l
1
11
X
11
∏
q
=
p
+
1
n
ϕ
q
l
q
+
m
-
p
ξ
b
g
j
t
.Based on the independence of the random vectorX
11 and the random variables (
ξ
j
,
p
+
1
≤
j
≤
n
), the inner products in (56) are rearranged as(57)
E
C
i
+
σ
i
1
X
i
1
Ψ
k
1
11
∏
q
=
p
+
1
n
ϕ
q
k
q
+
m
-
p
ξ
q
,
Ψ
l
1
11
∏
q
=
p
+
1
n
ϕ
q
l
q
+
m
-
p
ξ
q
=
m
i
1
Ψ
k
1
11
,
Ψ
l
1
11
11
+
σ
i
1
X
i
1
Ψ
k
1
11
,
Ψ
l
1
11
11
∏
q
=
p
+
1
n
α
k
q
+
m
-
p
q
δ
k
q
+
m
-
p
l
q
+
m
-
p
a
0
j
+
a
1
j
ξ
j
,
Ψ
l
1
11
∏
q
=
p
+
1
n
ϕ
q
l
q
+
m
-
p
ξ
q
=
δ
0
l
1
∏
q
=
p
+
1
q
≠
j
n
δ
0
l
q
+
m
-
p
a
0
j
δ
0
l
j
+
m
-
p
+
a
1
j
δ
j
+
1
+
m
-
p
l
j
+
m
-
p
,where δ
p
q is the Kronecker symbol.LetB
c be the canonical basis of R
X
0
,
…
,
X
m. This basis can be ordered in the following form:(58)
B
c
=
∏
k
=
0
m
X
k
α
k
-
α
k
+
1
:
w
h
e
r
e
α
0
=
q
,
0
≤
α
k
+
1
≤
α
k
,
α
m
+
1
=
0
;
q
∈
N
.The orthogonalization of the basisB
c gives the polynomial chaos Ψ
k
11
,
k
∈
N associated with the random vector X
11. This orthogonalization is performed using a classical method such as the Gram-Schmidt procedure in R
X
0
,
…
,
X
m and the following inner product: (59)
H
,
K
11
=
∫
R
m
+
1
H
x
K
x
g
11
x
d
x
∀
H
,
K
∈
R
X
11
2
.LetB
1 be the tensor defined for all α
∈
N
m
+
1, α
=
(
α
0
,
…
,
α
m
) by(60)
B
1
α
=
E
∏
k
=
0
m
X
k
α
k
=
∫
R
m
+
1
∏
k
=
0
m
x
k
α
k
g
11
x
d
x
.This integral can be computed using the distribution function of the vectorω
1 and the diffeomorphism h by(61)
B
1
α
=
∫
R
p
+
1
∏
k
=
0
m
h
k
ω
1
-
m
k
σ
k
α
k
f
1
ω
1
d
ω
1
.The integerM
, considered in decomposition (51), is chosen such that the considered polynomials have the degree less than or equal to L and expressed by(62)
M
=
ψ
n
+
m
-
p
0,0
,
…
,
0
,
L
∀
L
∈
N
.LetB
c
L be the canonical basis of R
L
X
0
,
…
,
X
m, then(63)
B
c
L
=
∏
k
=
0
m
X
k
α
k
-
α
k
+
1
:
where
α
0
=
q
,
0
≤
α
k
+
1
≤
α
k
,
α
m
+
1
=
0
;
q
∈
N
,
q
≤
L
.Let us denote byP
1 the passage matrix from the canonical basis B
c
L to the polynomial chaos basis associated with the vector X
1 of degree less than or equal to L. Using the following notations,(64)
N
L
=
ψ
m
0,0
,
…
,
0
,
L
,
Q
l
=
∏
i
=
0
m
X
i
l
i
w
h
e
r
e
l
0
,
…
,
l
m
=
ψ
m
-
1
l
for
0
≤
l
≤
N
L
,
V
c
11
=
Q
0
,
…
,
Q
N
L
T
,
V
C
P
11
=
Ψ
0
11
,
…
,
Ψ
N
L
11
T
,one gets(65)
V
C
P
11
=
P
1
T
V
c
11
.The matricesΥ, Τ
i and D, defined for all integers k and l between 0 and N
L, are introduced for 0
≤
i
≤
m, by(66)
Υ
k
+
1
,
l
+
1
=
Ψ
k
11
,
Ψ
l
11
11
=
α
k
11
δ
k
l
(67)
Τ
i
k
+
1
,
l
+
1
=
X
i
1
Ψ
k
11
,
Ψ
l
11
11
(68)
D
k
+
1
,
l
+
1
=
B
1
ψ
m
-
1
k
+
ψ
m
-
1
l
for
0
≤
k
≤
N
L
,
0
≤
l
≤
N
L
.Let(
e
i
)
0
≤
i
≤
m, (
m
+
1
) elements in N
m
+
1 such that(69)
e
0
=
1,0
,
…
,
0
;
e
1
=
0,1
,
…
,
0
;
⋮
e
m
=
0,0
,
…
,
1
.The matricesD
i of order (N
L
+
1) for 0
≤
i
≤
m are defined by (70)
D
i
k
+
1
,
l
+
1
=
B
1
ψ
m
-
1
k
+
ψ
m
-
1
l
+
e
i
for
0
≤
k
≤
N
L
,
0
≤
l
≤
N
L
.This allows defining the matricesΥ and Τ
i by (71)
Υ
=
P
1
D
P
1
T
,
Τ
i
=
P
1
D
i
P
1
T
.The random differential system (52) is thus reduced to the following deterministic one:(72)
∑
k
=
0
M
∑
i
=
0
m
E
C
i
Υ
k
1
+
1
,
l
1
+
1
+
σ
i
1
Τ
i
k
1
+
1
,
l
1
+
1
∏
q
=
p
+
1
n
α
k
q
+
m
-
p
q
δ
k
q
+
m
-
p
l
q
+
m
-
p
d
i
u
k
t
d
t
i
=
∑
j
=
p
+
1
n
δ
0
l
1
∏
q
=
p
+
1
q
≠
j
n
δ
0
l
q
+
m
-
p
a
0
j
δ
0
l
j
+
m
-
p
+
a
1
j
δ
j
+
1
+
m
-
p
l
j
+
m
-
p
g
j
t
,where k
1 and l
1 are integers related to k and l and introduced in (54).Finally, the last differential system is rewritten in the following condensed matrix form:(73)
∑
k
=
0
M
∑
i
=
0
m
H
i
k
1
+
1
,
l
1
+
1
d
i
u
k
t
d
t
i
=
E
l
t
,where (74)
H
i
k
1
+
1
,
l
1
+
1
=
E
C
i
Υ
k
1
+
1
,
l
1
+
1
+
σ
i
1
Τ
i
k
1
+
1
,
l
1
+
1
∏
q
=
p
+
1
n
α
k
q
+
m
-
p
q
δ
k
q
+
m
-
p
l
q
+
m
-
p
.
(75)
E
l
t
=
∑
j
=
p
+
1
n
δ
0
l
1
∏
q
=
p
+
1
q
≠
j
n
δ
0
l
q
+
m
-
p
a
0
j
δ
0
l
j
+
m
-
p
+
a
1
j
δ
j
+
1
+
m
-
p
l
j
+
m
-
p
g
j
tBased on the presented mathematical procedures, the truncated solution (51) can be numerically obtained.An equivalent form of (73) has been introduced in (37) based on original random variables. These two equations give the approximate solution ofU. The coefficients of these two equations are given, respectively, by the inner products in (57) for (73) and (32)-(33) for (37).It should be noted, on one hand, that these coefficients are hard to be determined for (37) but easier for (73). On the other hand, the number of the random variables used in (73) is largely reduced. This reduces the number of unknowns used to determine the solutionU and greatly simplifies the required numerical computation.Note that when uncertainties come from random parameters of the system parameters, they can be efficiently handled by the previous mathematical procedures. When the right-hand side excitation depends on many random variables, the previous internal coefficients procedure can be largely improved by the so-called superposition method.


## 2.3. Superposition Method

Thanks to the linear behavior of the considered system, the superposition principle can be applied to the standard random differential equation (13) and to the resulting equation from the internal random coefficients procedure (52). In both cases, the considered random right-hand side function F
(
t
) depends linearly on (
n
-
p
) random variables ξ
j. The superposition method consists of solving (
n
-
p
) differential equations with elementary right-hand sides given by(76)
L
-
·
U
t
,
ω
=
a
0
j
+
a
1
j
ξ
j
g
j
⁡
t
for
j
=
p
+
1
,
…
,
n
,where L
- may be the left-hand side differential operator of (13) or (52). In fact, this superposition is adding only one random variable ξ
j to the system and the resulting random vectors are (
ω
1
,
ξ
j
) for (13) and (
X
11
,
ξ
j
) for (52) for each j.As the internal random coefficient procedure leads to a reduced number of random variables, the superposition principle is applied to (52). Our focus is then on the resulting random differential equation: (77)
∑
i
=
0
m
E
C
i
+
σ
i
1
X
i
1
d
i
U
t
,
X
1
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
+
a
1
j
ξ
j
g
j
⁡
t
,with the deterministic initials conditions given by(78)
d
i
U
t
d
t
i
t
=
0
=
u
0
i
for
0
≤
i
≤
m
-
1
.The polynomial chaos associated with the random vectorω
2 is considered and denoted by ϕ
k
2
;
k
∈
N. The solution U is in this case a process that depends on t, on the random vector X
11, and on the random vector ω
2. The solution is then written in the following form: (79)
U
t
,
X
1
=
U
t
,
X
11
,
ω
2
.Since random vectorsX
11 and ω
2 are independent, the factorization by the polynomial chaos associated with the random vector ω
2 is used in the decomposition (50). The solution can be rewritten as follows:(80)
U
t
,
X
11
,
ω
2
=
∑
k
≥
0
U
k
t
,
X
11
ϕ
k
2
ω
2
.The insertion of this decomposition in the differential equation (77) and the projection of the equation and of the initial conditions over each polynomial chaos ϕ
l
2 for l
≥
0 lead to the following simplified system of random differential equations:Fork
=
0,(81a)
∑
i
=
0
m
E
C
i
+
σ
i
1
X
i
1
d
i
U
0
t
,
X
11
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
g
j
⁡
t
d
i
U
0
t
,
X
11
d
t
i
t
=
0
=
u
0
i
for
0
≤
i
≤
m
-
1
.
For1
≤
k
≤
n
-
p,(81b)
∑
i
=
0
m
E
C
i
+
σ
i
1
X
i
1
d
i
U
0
t
,
X
11
d
t
i
=
a
1
p
+
k
g
p
+
k
⁡
t
d
i
U
k
t
,
X
11
d
t
i
t
=
0
=
0
for
0
≤
i
≤
m
-
1
.
Fork
≥
n
-
p
+
1,(81c)
∑
i
=
0
m
E
C
i
+
σ
i
1
X
i
1
d
i
U
k
t
,
X
11
d
t
i
=
0
d
i
U
k
t
,
X
11
d
t
i
t
=
0
=
0
for
0
≤
i
≤
m
-
1
.It has to be noted that the casek
≥
n
-
p
+
1 leads to null initial conditions and null right-hand side. In this case, the solution is then null.The solution given in (80) is then reduced to (82)
U
t
,
X
11
,
ω
2
=
∑
k
=
0
n
-
p
ϕ
k
2
ω
2
U
k
t
,
X
11
.For0
≤
k
≤
n
-
p and considering the decomposition of the solution U
k
(
t
,
X
11
) according to the polynomials chaos associated with the random vector X
11, one gets(83)
U
k
t
,
X
11
=
∑
d
=
0
N
L
u
k
d
t
Ψ
d
11
X
11
.Based on this decomposition, the differential equations (81a)–(81c) are then reduced to the following:Fork
=
0,(84a)
∑
i
=
0
m
∑
d
=
0
N
L
Ψ
d
11
E
C
i
+
σ
i
1
X
i
1
d
i
u
0
d
t
d
t
i
=
∑
j
=
p
+
1
n
a
0
j
g
j
⁡
t
d
i
u
0
d
t
,
X
11
d
t
i
t
=
0
=
u
0
d
i
for
0
≤
i
≤
m
-
1
.
For1
≤
k
≤
n
-
p,(84b)
∑
i
=
0
m
∑
d
=
0
N
L
Ψ
d
11
E
C
i
+
σ
i
1
X
i
1
d
i
u
k
d
t
d
t
i
=
a
1
p
+
k
g
p
+
k
t
d
i
u
k
d
t
,
X
11
d
t
i
t
=
0
=
0
for
0
≤
i
≤
m
-
1
.Projecting this expression over each polynomial chaosΨ
q
11 for q
=
0 to N
L, one gets (85)
∑
i
=
0
m
∑
d
=
0
N
L
Ψ
d
11
E
C
i
+
σ
i
1
X
i
1
,
Ψ
q
11
d
i
u
0
d
t
d
t
i
=
δ
0
q
∑
j
=
p
+
1
n
a
0
j
g
j
tand, for 1
≤
k
≤
n
-
p,(86)
∑
i
=
0
m
∑
d
=
0
N
L
Ψ
d
11
E
C
i
+
σ
i
1
X
i
1
,
Ψ
q
11
d
i
u
k
d
t
d
t
i
=
δ
k
q
a
1
p
+
k
g
p
+
k
⁡
tUsing the matricesΥ and Τ
i introduced in (71), the last equations become the following.Forq
=
0 to N
L,(87)
∑
i
=
0
m
∑
d
=
0
N
L
E
C
i
Υ
d
1
+
1
,
q
1
+
1
+
σ
i
1
Τ
i
d
1
+
1
,
q
1
+
1
d
i
u
0
d
t
d
t
i
=
δ
0
q
∑
j
=
p
+
1
n
a
0
j
g
j
⁡
tand, for 1
≤
k
≤
n
-
p,(88)
∑
i
=
0
m
∑
d
=
0
N
L
E
C
i
Υ
d
1
+
1
,
q
1
+
1
+
σ
i
1
Τ
i
d
1
+
1
,
q
1
+
1
d
i
u
k
d
t
d
t
i
=
δ
kq
a
1
p
+
k
g
p
+
k
⁡
t
,where the integers d
1 and q
1 are related to d and q by notations (54).These deterministic systems are then rewritten in the following compact matrix form.Forq
=
0 to N
L,(89)
∑
i
=
0
m
∑
d
=
0
N
L
F
i
d
1
+
1
,
q
1
+
1
d
i
u
0
d
t
d
t
i
=
δ
0
q
∑
j
=
p
+
1
n
a
0
j
g
j
⁡
t
d
i
u
0
d
t
,
X
11
d
t
i
t
=
0
=
u
0
d
i
for
0
≤
i
≤
m
-
1and, for 1
≤
k
≤
n
-
p,(90)
∑
i
=
0
m
∑
d
=
0
N
L
F
i
d
1
+
1
,
q
1
+
1
d
i
u
k
d
t
d
t
i
=
a
1
p
+
k
g
p
+
k
⁡
t
δ
k
q
d
i
u
k
d
t
,
X
11
d
t
i
t
=
0
=
0
for
0
≤
i
≤
m
-
1
,where(91)
F
i
=
E
C
i
Υ
+
σ
i
1
Τ
i
.The final truncated solution is then given by(92a)
U
t
=
∑
k
=
0
n
-
p
ϕ
k
2
ω
2
∑
d
=
0
N
L
u
k
d
t
Ψ
d
11
X
11
.Its mean and variance are given:(92b)
m
N
t
=
u
00
t
,
V
N
t
=
∑
k
=
1
n
-
p
∑
d
=
0
N
α
d
11
u
k
d
2
t
+
∑
d
=
1
N
α
d
11
u
0
d
2
t
.This solution is the main approximate solution of the random initial value problem (1) elaborated herein. Following (89) and (90), one has to solve NL(
n
-
p
+
1
) simple initial values differential systems. These systems can be solved by standard methods such as the Runge-Kutta method.The main advantages of this methodological approach are as follows:
(i)
The deterministic differential systems have the same left-hand side that is built only once and in a compact matrix form.
(ii)
Simple right-hand sides are obtained for an arbitrary number of random variables defining the excitation.
(iii)
A large reduction of the number of unknowns based on polynomial chaos can result.
(iv)
The CPU time and memory needed can be largely reduced.
(v)
The numerical solution can be done in parallel manner and particularly for a large number of uncertain parameters.


## 2.4. Conditional Expectation Method

For the sake of comparison, a methodological approach based on the conditional expectation method is also elaborated here to solve (1). Remember that the solution U of (1) is a stochastic process that depends on the random vectors ω
1 and ω
2 and time t. Let ω
0
=
(
ω
10
,
ω
20
) be an observation of the random vector ω
=
(
ω
1
,
ω
2
). The conditional expectation of U such that ω
=
ω
0, noted by E
U
∣
ω
=
ω
0, is the solution of the deterministic equation:(93)
L
0
u
t
,
ω
0
=
F
t
,
ω
20
,where L
0 is a deterministic operator defined by(94)
L
0
=
∑
i
=
0
m
C
i
ω
10
d
i
d
t
i
.Like in the first section, the random vectorω is assumed to have a density function relative to the Lebesgue measure denoted by f. The moments of order p (p
∈
N
∗) of the solution of (1) are then given by the Bayes formulas:(95)
m
p
t
=
∫
R
n
+
1
E
U
∣
ω
=
ω
0
p
f
ω
0
d
ω
0
.The numerical computation of the integral (95) needs the solution of (93) for any observation vector ω
0. The analytical solution can be obtained in general cases. Assume that ω
0
∈
Ω, where Ω is a bounded domain in R
n
+
1.The integral can be approximated by the Gauss-Legendre quadrature method. Let the setsω
-
0
i
,
1
≤
i
≤
q and p
-
i
,
1
≤
i
≤
q be Gauss points and associated Gauss weights. The moments m
p
t can be approximated by (96)
m
p
t
≅
∑
i
=
1
q
p
-
i
E
U
∣
ω
=
ω
-
0
i
p
f
ω
0
i
.The main advantage of this method is that it allows obtaining a general solution that is considered here as a reference solution as well as the obtained one based on the Monte-Carlo method. But, the inconvenience of these two methods is that a very large CPU time is needed for accurate numerical solutions.


## 3. Applications

For the sake of clarity, the developed methodological approaches are explicitly presented for the most standard initial value problem. Let us consider the following second-order differential equation modeling the dynamical forced behavior of Euler-Bernoulli beams with the following uncertain parameters: thicknessh, mass density ρ, and Young modulus E:(97)
C
2
U
¨
t
+
C
1
U
˙
t
+
C
0
U
t
=
F
t
,
U
0
=
U
O
,
U
˙
0
=
U
˙
0where(98a)
C
2
=
ρ
h
C
1
=
α
C
2
+
β
C
0
C
0
=
c
0
E
h
+
c
1
E
h
3
(98b)
F
t
=
∑
i
=
3
n
a
0
i
+
a
1
i
ξ
i
sin
⁡
θ
i
t
,in which c
0, c
1, α, β, a
0
i, a
1
i, θ
i, U
0, and U
˙
0 are deterministic constants. The considered physical random variables ρ, h, and E are written in the following forms:(98c)
ρ
=
ρ
0
+
ξ
0
ρ
1
h
=
h
0
+
ξ
1
h
1
E
=
E
0
+
ξ
2
E
1
,where (
ξ
i
,
0
≤
i
≤
n
) are random variables with zero means and one variance, two by two independent, and possessing distribution functions relative to the Lebesgue measure, respectively, denoted by f
i, 0
≤
i
≤
n.It has to be noted that the internal random coefficientsC
0, C
1, and C
2 depend nonlinearly on the random variables ξ
0, ξ
1, and ξ
2 and n is an arbitrary integer. For a clear presentation, the three methodological approaches, elaborated here, are developed for the solution of (97).


### 3.1. Generalized Polynomial Chaos Method

The random vectorω
1
,
ω
1
=
(
ξ
0
,
ξ
1
,
ξ
2
), the order of linear random operator L is m
=
2, and the internal coefficients in this case are C
0
ω
1, C
1
ω
1, C
2
ω
1. Note that these coefficients depend nonlinearly on ξ
0, ξ
1, and ξ
2 and the order of nonlinearity is 2 for C
2 and 4 for C
0 and C
1. Based on the notations of the first section, these coefficients can be expanded as follows:(99)
C
0
ω
1
=
c
000
0
+
c
010
0
ξ
1
+
c
001
0
ξ
2
+
c
020
0
ξ
1
2
+
c
011
0
ξ
1
ξ
2
+
c
030
0
ξ
1
3
+
c
021
0
ξ
1
2
ξ
2
+
c
031
0
ξ
1
3
ξ
2
C
1
ω
1
=
c
000
1
+
c
100
1
ξ
0
+
c
010
1
ξ
1
+
c
001
1
ξ
2
+
c
110
1
ξ
0
ξ
1
+
c
020
1
ξ
1
2
+
c
011
1
ξ
1
ξ
2
+
c
030
1
ξ
1
3
+
c
021
1
ξ
1
2
ξ
2
+
c
031
1
ξ
1
3
ξ
2
,
C
2
ω
1
=
c
000
2
+
c
100
2
ξ
0
+
c
010
2
ξ
1
+
c
110
2
ξ
0
ξ
1
,where these coefficients are given by(100)
c
000
0
=
c
0
E
0
h
0
+
c
1
E
0
h
0
3
,
c
010
0
=
c
0
E
0
h
1
+
3
c
1
E
0
h
0
2
h
1
,
c
001
0
=
c
0
E
1
h
0
+
c
1
E
1
h
0
3
,
c
020
0
=
3
c
1
E
0
h
0
h
1
2
,
c
011
0
=
c
0
E
1
h
1
+
3
c
1
E
1
h
0
2
h
1
,
c
003
0
=
c
1
E
0
h
1
3
,
c
021
0
=
3
c
1
E
1
h
0
h
1
2
,
c
031
0
=
c
1
E
1
h
1
3
;
c
000
2
=
ρ
0
h
0
,
c
100
2
=
ρ
1
h
0
,
c
010
2
=
ρ
0
h
1
,
c
110
2
=
ρ
1
h
1
;
c
000
1
=
α
c
000
2
+
β
c
000
0
,
c
100
1
=
α
c
100
2
,
c
010
1
=
α
c
010
2
+
β
c
010
0
,
c
001
1
=
β
c
001
0
,
c
110
1
=
α
c
110
2
,
c
020
1
=
β
c
020
0
,
c
011
1
=
β
c
011
0
,
c
003
1
=
β
c
003
0
,
c
021
1
=
β
c
021
0
,
c
031
1
=
β
c
031
0
.Letk and l
=
0 to N, and, using the notation (23), one gets k
0
,
…
,
k
n
=
ψ
n
-
1
(
k
) and l
0
,
…
,
l
n
=
ψ
n
-
1
(
l
), where k
0
,
…
,
k
n and l
0
,
…
,
l
n are in N
n
+
1.The matricesT
i for i
=
0 to 2, introduced in (35), are then given by the following.Fork, l
=
0
:
N,(101)
T
0
l
+
1
,
k
+
1
=
c
000
0
∏
d
=
0
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
010
0
β
11
l
1
+
1
,
k
1
+
1
∏
d
=
0
d
≠
1
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
001
0
β
21
l
2
+
1
,
k
2
+
1
∏
d
=
0
d
≠
2
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
020
0
β
12
l
1
+
1
,
k
1
+
1
∏
d
=
0
d
≠
1
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
011
0
β
11
l
1
+
1
,
k
1
+
1
β
21
l
2
+
1
,
k
2
+
1
∏
d
=
0
d
≠
1,2
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
030
0
β
13
l
1
+
1
,
k
1
+
1
∏
d
=
0
d
≠
1
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
021
0
β
12
l
1
+
1
,
k
1
+
1
β
21
l
2
+
1
,
k
2
+
1
∏
d
=
0
d
≠
1,2
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
031
0
β
13
l
1
+
1
,
k
1
+
1
β
21
l
2
+
1
,
k
2
+
1
∏
d
=
0
d
≠
1,2
n
β
d
0
l
d
+
1
,
k
d
+
1
T
2
l
+
1
,
k
+
1
=
c
000
2
∏
d
=
0
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
100
2
β
01
l
1
+
1
,
k
1
+
1
∏
d
=
1
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
010
2
β
11
l
1
+
1
,
k
1
+
1
∏
d
=
0
d
≠
1
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
110
2
β
01
l
1
+
1
,
k
1
+
1
β
11
l
1
+
1
,
k
1
+
1
∏
d
=
2
n
β
d
0
l
d
+
1
,
k
d
+
1
,
T
1
=
α
T
0
+
β
T
2
.The variable vectorD of dimension (
N
+
1
), introduced in (36), is given by the following.Fork
=
0 to N,(102)
D
k
+
1
t
=
∑
j
=
3
n
a
0
j
β
j
0
1
,
k
j
+
1
+
a
1
j
β
j
1
1
,
k
j
+
1
∏
d
=
0
d
≠
j
n
β
d
0
1
,
k
d
+
1
s
i
n
⁡
θ
j
t
.The passage from the index for0
≤
l
≤
N to the multi-indices (l
0
,
l
1
,
…
,
l
n), used for the determination the matrices T
0, T
1, and T
2 needs the bijection ψ
n. This bijection is given with Algorithm 1.Algorithm 1
Read L and n
L non null integer
t
=
0
.
For
l
0
=
0
:
L
;
l
1
=
0
:
l
0
;
l
2
=
0
:
l
1
;
…
;
l
n
=
0
:
l
n
-
1
.
ψ
n
(
l
0
-
l
1
,
l
1
-
l
2
,
…
,
l
n
-
1
-
l
n
,
l
n
)
=
t
+
1
;
t
=
t
+
1
;
end.The determination ofT
0, T
1, and T
2, in this case, needs the computation of the inner products ξ
i
p
ϕ
k
,
ϕ
l for p
≤
3, by (32).If the random coefficientsξ
i are not identically distributed then ξ
i
p
ϕ
k
,
ϕ
l have to be computed for each i
=
0 to n. Otherwise, these scalar products have to be computed only once. This difficulty as well as the number of unknowns can be largely reduced by using the internal random coefficients method.


### 3.2. Internal Random Coefficients Method

Based on this method, the internal coefficients of (97) are rewritten in the following form:(103)
C
0
ω
1
=
E
C
0
+
σ
0
1
X
0
1
C
2
ω
1
=
E
C
0
+
σ
2
1
X
2
1
C
1
ω
1
=
α
C
0
ω
1
+
β
C
2
ω
1
.The random vectorX
11, introduced by this method, is then given by(104)
X
11
=
X
0
1
,
X
2
1
.The number of random variables used here is(
n
-
1
) variables, two variables for the coefficients C
j, and (
n
-
3
) variables for the right-hand side.Fork
,
l
=
0
,
…
,
N
L; and let k
1, l
1 be the integers associated with integers k, l, introduced by the notations (23)–(54); one has(105)
k
0
,
…
,
k
n
-
1
=
ψ
n
-
1
-
1
k
;
l
0
,
…
,
l
n
-
1
=
ψ
n
-
1
-
1
l
,
k
1
=
ψ
2
k
0
,
k
1
,
l
1
=
ψ
2
l
0
,
l
1
.The matricesH
i, introduced in (74), are given in this case by(106)
H
0
k
1
,
l
1
=
E
C
0
Υ
k
1
+
1
,
l
1
+
1
+
σ
0
1
Τ
0
k
1
+
1
,
l
1
+
1
∏
d
=
2
n
-
1
β
d
0
l
d
+
1
,
k
d
+
1
H
2
k
1
,
l
1
=
E
C
1
Υ
k
1
+
1
,
l
1
+
1
+
σ
1
1
Τ
2
k
1
+
1
,
l
1
+
1
∏
d
=
2
n
-
1
β
d
0
l
d
+
1
,
k
d
+
1
,
H
1
=
α
H
0
+
β
H
2and the variable vector E of dimension (
N
L
+
1
), introduced in (75), is given by(107)
E
l
t
=
∑
j
=
3
n
δ
0
l
1
∏
q
=
3
q
≠
j
n
δ
0
l
q
-
1
a
0
j
δ
0
l
j
-
1
+
a
1
j
δ
j
l
j
-
1
sin
⁡
θ
j
t
l
=
1
,
N
L
+
1
.The number of random variables,n, is arbitrary considered. For more clarity about the needed operations, let us consider n
=
5. Based on the classical polynomial chaos and if the decomposition of the solution into the polynomials chaos of degree less than or equal to seven is considered then the needed unknowns to be determined are of number 13
7
=
1716, while the internal random coefficient method requires 12
7
=
792 unknowns to be determined. The number unknowns of is thus diminished by half. In addition the determination of the tensors T
0, T
1, T
2 needs the computation of the scalar products ∏
d
=
0
p
ξ
d
α
d
ψ
k
,
ψ
l where ∑
d
=
0
p
α
d
≤
q
i and q
i is the degree of a random polynomial coefficient C
i
(
ω
1
) given by relationship (10), while the tensors H
0, H
1, H
2 need only the calculation of the inner products X
i
1
p
Ψ
k
1
,
Ψ
l
1 for p
=
0 or 1. This demonstrates an extra reduction in the computation effort. For more reduction in the CPU time, the internal random coefficient method is coupled with the superposition method.


### 3.3. Superposition Method

Here, (97) is replaced with an equivalent system with (
n
-
2
) equations; each equation is a random equation depending on the random vector X
11 introduced in (46). These equations are given by(108a)
∑
i
=
0
2
C
i
ω
1
d
i
U
0
t
,
X
11
d
t
i
=
∑
j
=
3
n
a
0
j
sin
⁡
θ
j
twith the initial conditions given by d
U
0
t
,
X
11
/
d
t
t
=
0
=
u
˙
0
.
U
0
0
,
X
11
=
u
0,And, for1
≤
k
≤
n
-
2,(108b)
∑
i
=
0
2
C
i
ω
1
d
i
U
k
t
,
X
11
d
t
i
=
a
1
2
+
k
sin
⁡
θ
2
+
k
twith the initial conditions given by(109)
d
U
k
t
,
X
11
d
t
t
=
0
=
0
,
U
k
0
,
X
11
=
0
.The associated matrices, introduced in (91), are given by(110)
F
0
=
E
C
0
Υ
+
σ
0
1
Τ
0
F
2
=
E
C
2
Υ
+
σ
2
1
Τ
2
,
F
1
=
α
F
0
+
β
F
2
.These matrices are replaced with their expressions in systems (108a) and (108b). In this case, the determination of the tensorsF
0,F
1,F
2 needs the determination of inner products X
i
11
p
Ψ
k
1
11
,
Ψ
l
1
11 for p
=
0 or 1 where k
1 and l1 are the integers associated with integers k and l given by notation (54). In this case, the polynomial chaos Ψ
k
1
11 associated with the random vector X
11 is simple. The solution of the deterministic system defined by (86)-(87) determines an approximate solution U of the initial random equation (97).This system largely reduces the number of unknowns to be determined for the approximate solutionU and makes the use of polynomial chaos possible even when the classical projection method is unable to determine the solution. The CPU time and memory required to solve this problem by classical projection methods are significantly reduced by this method.It is well known that, for a stochastic excitation and with Karhunen-Loéve expansion, the classical projection methods are limited to a small number of random variables used in the decomposition. This will affect the accuracy of the obtained solution. This kind of problems is completely solved the by the presented method.


## 3.1. Generalized Polynomial Chaos Method

The random vectorω
1
,
ω
1
=
(
ξ
0
,
ξ
1
,
ξ
2
), the order of linear random operator L is m
=
2, and the internal coefficients in this case are C
0
ω
1, C
1
ω
1, C
2
ω
1. Note that these coefficients depend nonlinearly on ξ
0, ξ
1, and ξ
2 and the order of nonlinearity is 2 for C
2 and 4 for C
0 and C
1. Based on the notations of the first section, these coefficients can be expanded as follows:(99)
C
0
ω
1
=
c
000
0
+
c
010
0
ξ
1
+
c
001
0
ξ
2
+
c
020
0
ξ
1
2
+
c
011
0
ξ
1
ξ
2
+
c
030
0
ξ
1
3
+
c
021
0
ξ
1
2
ξ
2
+
c
031
0
ξ
1
3
ξ
2
C
1
ω
1
=
c
000
1
+
c
100
1
ξ
0
+
c
010
1
ξ
1
+
c
001
1
ξ
2
+
c
110
1
ξ
0
ξ
1
+
c
020
1
ξ
1
2
+
c
011
1
ξ
1
ξ
2
+
c
030
1
ξ
1
3
+
c
021
1
ξ
1
2
ξ
2
+
c
031
1
ξ
1
3
ξ
2
,
C
2
ω
1
=
c
000
2
+
c
100
2
ξ
0
+
c
010
2
ξ
1
+
c
110
2
ξ
0
ξ
1
,where these coefficients are given by(100)
c
000
0
=
c
0
E
0
h
0
+
c
1
E
0
h
0
3
,
c
010
0
=
c
0
E
0
h
1
+
3
c
1
E
0
h
0
2
h
1
,
c
001
0
=
c
0
E
1
h
0
+
c
1
E
1
h
0
3
,
c
020
0
=
3
c
1
E
0
h
0
h
1
2
,
c
011
0
=
c
0
E
1
h
1
+
3
c
1
E
1
h
0
2
h
1
,
c
003
0
=
c
1
E
0
h
1
3
,
c
021
0
=
3
c
1
E
1
h
0
h
1
2
,
c
031
0
=
c
1
E
1
h
1
3
;
c
000
2
=
ρ
0
h
0
,
c
100
2
=
ρ
1
h
0
,
c
010
2
=
ρ
0
h
1
,
c
110
2
=
ρ
1
h
1
;
c
000
1
=
α
c
000
2
+
β
c
000
0
,
c
100
1
=
α
c
100
2
,
c
010
1
=
α
c
010
2
+
β
c
010
0
,
c
001
1
=
β
c
001
0
,
c
110
1
=
α
c
110
2
,
c
020
1
=
β
c
020
0
,
c
011
1
=
β
c
011
0
,
c
003
1
=
β
c
003
0
,
c
021
1
=
β
c
021
0
,
c
031
1
=
β
c
031
0
.Letk and l
=
0 to N, and, using the notation (23), one gets k
0
,
…
,
k
n
=
ψ
n
-
1
(
k
) and l
0
,
…
,
l
n
=
ψ
n
-
1
(
l
), where k
0
,
…
,
k
n and l
0
,
…
,
l
n are in N
n
+
1.The matricesT
i for i
=
0 to 2, introduced in (35), are then given by the following.Fork, l
=
0
:
N,(101)
T
0
l
+
1
,
k
+
1
=
c
000
0
∏
d
=
0
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
010
0
β
11
l
1
+
1
,
k
1
+
1
∏
d
=
0
d
≠
1
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
001
0
β
21
l
2
+
1
,
k
2
+
1
∏
d
=
0
d
≠
2
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
020
0
β
12
l
1
+
1
,
k
1
+
1
∏
d
=
0
d
≠
1
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
011
0
β
11
l
1
+
1
,
k
1
+
1
β
21
l
2
+
1
,
k
2
+
1
∏
d
=
0
d
≠
1,2
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
030
0
β
13
l
1
+
1
,
k
1
+
1
∏
d
=
0
d
≠
1
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
021
0
β
12
l
1
+
1
,
k
1
+
1
β
21
l
2
+
1
,
k
2
+
1
∏
d
=
0
d
≠
1,2
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
031
0
β
13
l
1
+
1
,
k
1
+
1
β
21
l
2
+
1
,
k
2
+
1
∏
d
=
0
d
≠
1,2
n
β
d
0
l
d
+
1
,
k
d
+
1
T
2
l
+
1
,
k
+
1
=
c
000
2
∏
d
=
0
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
100
2
β
01
l
1
+
1
,
k
1
+
1
∏
d
=
1
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
010
2
β
11
l
1
+
1
,
k
1
+
1
∏
d
=
0
d
≠
1
n
β
d
0
l
d
+
1
,
k
d
+
1
+
c
110
2
β
01
l
1
+
1
,
k
1
+
1
β
11
l
1
+
1
,
k
1
+
1
∏
d
=
2
n
β
d
0
l
d
+
1
,
k
d
+
1
,
T
1
=
α
T
0
+
β
T
2
.The variable vectorD of dimension (
N
+
1
), introduced in (36), is given by the following.Fork
=
0 to N,(102)
D
k
+
1
t
=
∑
j
=
3
n
a
0
j
β
j
0
1
,
k
j
+
1
+
a
1
j
β
j
1
1
,
k
j
+
1
∏
d
=
0
d
≠
j
n
β
d
0
1
,
k
d
+
1
s
i
n
⁡
θ
j
t
.The passage from the index for0
≤
l
≤
N to the multi-indices (l
0
,
l
1
,
…
,
l
n), used for the determination the matrices T
0, T
1, and T
2 needs the bijection ψ
n. This bijection is given with Algorithm 1.Algorithm 1
Read L and n
L non null integer
t
=
0
.
For
l
0
=
0
:
L
;
l
1
=
0
:
l
0
;
l
2
=
0
:
l
1
;
…
;
l
n
=
0
:
l
n
-
1
.
ψ
n
(
l
0
-
l
1
,
l
1
-
l
2
,
…
,
l
n
-
1
-
l
n
,
l
n
)
=
t
+
1
;
t
=
t
+
1
;
end.The determination ofT
0, T
1, and T
2, in this case, needs the computation of the inner products ξ
i
p
ϕ
k
,
ϕ
l for p
≤
3, by (32).If the random coefficientsξ
i are not identically distributed then ξ
i
p
ϕ
k
,
ϕ
l have to be computed for each i
=
0 to n. Otherwise, these scalar products have to be computed only once. This difficulty as well as the number of unknowns can be largely reduced by using the internal random coefficients method.


## 3.2. Internal Random Coefficients Method

Based on this method, the internal coefficients of (97) are rewritten in the following form:(103)
C
0
ω
1
=
E
C
0
+
σ
0
1
X
0
1
C
2
ω
1
=
E
C
0
+
σ
2
1
X
2
1
C
1
ω
1
=
α
C
0
ω
1
+
β
C
2
ω
1
.The random vectorX
11, introduced by this method, is then given by(104)
X
11
=
X
0
1
,
X
2
1
.The number of random variables used here is(
n
-
1
) variables, two variables for the coefficients C
j, and (
n
-
3
) variables for the right-hand side.Fork
,
l
=
0
,
…
,
N
L; and let k
1, l
1 be the integers associated with integers k, l, introduced by the notations (23)–(54); one has(105)
k
0
,
…
,
k
n
-
1
=
ψ
n
-
1
-
1
k
;
l
0
,
…
,
l
n
-
1
=
ψ
n
-
1
-
1
l
,
k
1
=
ψ
2
k
0
,
k
1
,
l
1
=
ψ
2
l
0
,
l
1
.The matricesH
i, introduced in (74), are given in this case by(106)
H
0
k
1
,
l
1
=
E
C
0
Υ
k
1
+
1
,
l
1
+
1
+
σ
0
1
Τ
0
k
1
+
1
,
l
1
+
1
∏
d
=
2
n
-
1
β
d
0
l
d
+
1
,
k
d
+
1
H
2
k
1
,
l
1
=
E
C
1
Υ
k
1
+
1
,
l
1
+
1
+
σ
1
1
Τ
2
k
1
+
1
,
l
1
+
1
∏
d
=
2
n
-
1
β
d
0
l
d
+
1
,
k
d
+
1
,
H
1
=
α
H
0
+
β
H
2and the variable vector E of dimension (
N
L
+
1
), introduced in (75), is given by(107)
E
l
t
=
∑
j
=
3
n
δ
0
l
1
∏
q
=
3
q
≠
j
n
δ
0
l
q
-
1
a
0
j
δ
0
l
j
-
1
+
a
1
j
δ
j
l
j
-
1
sin
⁡
θ
j
t
l
=
1
,
N
L
+
1
.The number of random variables,n, is arbitrary considered. For more clarity about the needed operations, let us consider n
=
5. Based on the classical polynomial chaos and if the decomposition of the solution into the polynomials chaos of degree less than or equal to seven is considered then the needed unknowns to be determined are of number 13
7
=
1716, while the internal random coefficient method requires 12
7
=
792 unknowns to be determined. The number unknowns of is thus diminished by half. In addition the determination of the tensors T
0, T
1, T
2 needs the computation of the scalar products ∏
d
=
0
p
ξ
d
α
d
ψ
k
,
ψ
l where ∑
d
=
0
p
α
d
≤
q
i and q
i is the degree of a random polynomial coefficient C
i
(
ω
1
) given by relationship (10), while the tensors H
0, H
1, H
2 need only the calculation of the inner products X
i
1
p
Ψ
k
1
,
Ψ
l
1 for p
=
0 or 1. This demonstrates an extra reduction in the computation effort. For more reduction in the CPU time, the internal random coefficient method is coupled with the superposition method.


## 3.3. Superposition Method

Here, (97) is replaced with an equivalent system with (
n
-
2
) equations; each equation is a random equation depending on the random vector X
11 introduced in (46). These equations are given by(108a)
∑
i
=
0
2
C
i
ω
1
d
i
U
0
t
,
X
11
d
t
i
=
∑
j
=
3
n
a
0
j
sin
⁡
θ
j
twith the initial conditions given by d
U
0
t
,
X
11
/
d
t
t
=
0
=
u
˙
0
.
U
0
0
,
X
11
=
u
0,And, for1
≤
k
≤
n
-
2,(108b)
∑
i
=
0
2
C
i
ω
1
d
i
U
k
t
,
X
11
d
t
i
=
a
1
2
+
k
sin
⁡
θ
2
+
k
twith the initial conditions given by(109)
d
U
k
t
,
X
11
d
t
t
=
0
=
0
,
U
k
0
,
X
11
=
0
.The associated matrices, introduced in (91), are given by(110)
F
0
=
E
C
0
Υ
+
σ
0
1
Τ
0
F
2
=
E
C
2
Υ
+
σ
2
1
Τ
2
,
F
1
=
α
F
0
+
β
F
2
.These matrices are replaced with their expressions in systems (108a) and (108b). In this case, the determination of the tensorsF
0,F
1,F
2 needs the determination of inner products X
i
11
p
Ψ
k
1
11
,
Ψ
l
1
11 for p
=
0 or 1 where k
1 and l1 are the integers associated with integers k and l given by notation (54). In this case, the polynomial chaos Ψ
k
1
11 associated with the random vector X
11 is simple. The solution of the deterministic system defined by (86)-(87) determines an approximate solution U of the initial random equation (97).This system largely reduces the number of unknowns to be determined for the approximate solutionU and makes the use of polynomial chaos possible even when the classical projection method is unable to determine the solution. The CPU time and memory required to solve this problem by classical projection methods are significantly reduced by this method.It is well known that, for a stochastic excitation and with Karhunen-Loéve expansion, the classical projection methods are limited to a small number of random variables used in the decomposition. This will affect the accuracy of the obtained solution. This kind of problems is completely solved the by the presented method.


## 4. Numerical Results

The dynamical system defined in (97) is considered. To analyze the effect of the uncertain parameters, the free system (
F
(
t
)
=
0
) is first considered. For free systems, two variants of the internal random coefficients method are used. The first one consists of using two random coefficients C
0 and C
2 called here IRCM2; the second one consists of dividing by C
2 and then using only one random coefficient (
C
2
/
C
0
). This variant is referred to here as IRCM1.


### 4.1. Free Dynamics Response

For numerical analysis, the following values are used:(111)
ρ
0
=
7800
,
ρ
1
=
α
0
ρ
0
,
h
0
=
5
∗
10
-
5
,
h
1
=
α
1
h
0
,
E
0
=
21
∗
10
10
,
E
1
=
α
2
E
0
,
α
=
0.01
,
β
=
0.001
.The random variablesξ
0, ξ
1, ξ
2 are assumed to follow various laws and α
0, α
1, and α
2 are the ratios of the standard deviation and the mean of the random variables ρ, h, and E.Figures1(a) and 1(b) present the standard deviation and mean of the random solution U, where ξ
i are random variables with normalized truncated exponential laws in [
0,10
], with 5% of standard deviation for parameters (
α
0
=
α
1
=
α
2
=
5
%). The solution U
(
t
) is computed using the chaos methods associated with the random variables ξ
i, the IRCM1, IRCM2, the conditional expectation method, and the Monte-Carlo method with 1000 simples.Figure 1
(a) Standard deviation ofU when ξ
i follow the truncated exponential law in [
0,10
] and σ
i = 5%. (b) Means of U when ξ
i follow the truncated exponential law in [
0,10
] and σ
i = 5%.
(a)
(b)In the case of the classical chaos method, solutionU is expanded in the polynomial chaos of degree less than or equal to 7 and these polynomials are dependent on the three random variables ξ
i for i
=
0, 1, and 2. For the first internal random coefficient method, IRCM1, the coefficients of the random equation (97) are divided by the random variables C
0
(
ω
1
) and the random equation depends only on one new random variable η
=
C
2
ω
1
/
C
0
ω
1. The solution U is then expanded into the polynomials chaos associated with the random variable η of degree less than or equal to 12. For the second internal random coefficient method (IRCM2), the random equation depends on the random variables C
0
ω
1
,
C
2
ω
2. Solution Uis here expanded into the polynomials chaos associated with these two random variables of degree less than or equal to 2 and also less than 7. The number of unknowns used by the classical chaos method is 120 variables for degrees ≤ 7. The ICRM2 needs 6 unknowns for degrees ≤ 2 and 36 unknowns for degrees ≤ 7 while the IRCM1 needs only 13 unknowns for degrees ≤ 12 as presented in Table 2. This simple case demonstrates clearly how the IRCM reduces the number of unknowns.Table 2
Number of coefficients to compute for IRCM1-IRCM2 and PC at the considered degrees.
Method
IRCM1
IRCM2
IRCM2
Chaos
Degree
12
2
7
7
Unknowns
13
6
36
120The mean and standard deviation of the solutionU obtained by IRCM1, IRCM2, chaos, conditional expectation, and Monte-Carlo methods, presented in Figures 1(a) and 1(b), demonstrate a good accuracy between the predictions obtained by these methods.Figures2(a) and 2(b) present the standard deviation and mean of the solution U, when the random variables ξ
i are uniform in [
-
4,4
] with 5% standard deviation for each parameter (
α
0
=
α
1
=
α
2
=
5
%
). Solution U is computed using the previous methods. For the case of chaos method, the solution U is expanded into polynomials chaos of degree less than or equal to seven and these polynomials are dependent on the three random variables ξ
i for i
=
0, 1, and 2. For the ICRM1, solution U is expanded into the polynomial chaos associated with the random variable η of degree less than or equal to 12. For ICRM2, solution U is expanded into the polynomial chaos associated with two random variables of degree less than or equal to 7.Figure 2
(a) Standard deviation ofU when ξ
i follow the uniform law in [
-
4,4
] and σ
i = 5%. (b) Means of U when ξ
i follow the uniform law in [
-
4,4
] and σ
i = 5%.
(a)
(b)For large standard deviations, these methods are also tested. Figures3(a) and 3(b) present the standard deviation and mean of the solution U, when the random variables ξ
i are uniform in [-
4,4
] with 20% standard deviation (α
0
=
α
1
=
α
2
=
20
%). For the case of chaos method, solution U is expanded in polynomial chaos of degree less than or equal to 7 and these polynomials are dependent on the three random variables ξ
i for i
=
0, 1, and 2. ICRM1 and ICRM2 are also used with polynomial of degree ≤ 7. The number of coefficients to compute is presented in Table 3. Results obtained by the conditional expectation and Monte-Carlo methods are also presented. It is observed that, even for large standard deviation, 20%, a good agreement is observed between the predictions obtained by these methods.Table 3
Number of coefficients to compute for IRCM1-IRCM2 and PC at degree 7 for the considered number of variables.
Method
IRCM1
IRCM2
Chaos
Degree
7
7
7
Unknowns
8
36
120Figure 3
(a) Standard deviation ofU when ξ
i follow the uniform law in [
-
4,4
] and σ
i = 20%. (b) Means of U when ξ
i follow the uniform law in [
-
4,4
] and σ
i = 20%.
(a)
(b)Figures4(a) and 4(b) present the standard deviation and mean of solution U, where the random variables ξ
i follow a truncated exponential law in [
0
,
10
] with 20% standard deviation (α
0
=
α
1
=
α
2
=
20
%). The degree of the used polynomials is ≤7 for chaos, ≤12 for IRCM1, and ≤2 and ≤7 for IRCM2 and the number of coefficients to compute is presented in Table 4. It is demonstrated that, for the mean value, all the methods agree with each other but there are some discrepancies for the standard deviation. The IRCM2 is better than IRCM1 and particularly for large standard deviations. The division procedure, used in IRCM1, reduces the number of random parameters by one but introduces additional nonlinear random effect. This effect affects the accuracy of the IRCM1, as clearly presented in Figure 3.Table 4
Number of coefficients to compute for IRCM1-IRCM2 and PC at different degrees.
Method
IRCM1
IRCM2
IRCM2
Chaos
Degree
12
2
7
7
Unknowns
13
6
36
120Figure 4
(a) Standard deviation ofU when ξ
i follow the truncated exponential law in [
0,10
] and σ
i = 20%. (b) Means of U when ξ
i follow the truncated exponential law in [
0,10
] and σ
i = 20%.
(a)
(b)


### 4.2. Dynamical System with a Random Excitation

To test the flexibility of the presented methodological approach to handle dynamical systems with random variables following different laws, (97) is considered with a random excitation given by F
(
t
)
=
a
0
+
a
1
ξ
3
sin
⁡
(
ω
0
t
). a
0, a
1, and ω
0 are considered deterministic. The random coefficients ξ
0, ξ
1, ξ
2, and ξ
3 are assumed here to follow different laws.Figures5(a) and 5(b) present the mean and standard deviation of U
(
t
) obtained by the presented methods. ξ
0 and ξ
2 are truncated normal variables in [
-
4,4
], ξ
1 is a uniform variable in [
-
1,1
], and ξ
3 is a truncated exponential variable in [
0
,
10
] with 5% of standard deviations.Figure 5
(a) Standard deviation of solution whereξ
0 and ξ
2 are truncated normal variables in [
-
4,4
], ξ
1 is a uniform variable in [
-
1,1
], and ξ
3 is a truncated exponential variable in [
0,10
] and σ
i = 5%. (b) Means of solution where ξ
0 and ξ
2 are truncated normal variables in [
-
4,4
], ξ
1 is a uniform variable in [
-
1,1
], and ξ
3 is a truncated exponential variable in [
0,10
] and σ
i = 5%.
(a)
(b)The presented results are obtained by the chaos method where solutionU is expanded in the polynomials chaos of degree less than or equal to 2 and 5 and these polynomials depend on the four random variables ξ
i for i
=
0, 1, 2, and 3. For internal random coefficients method, solutionU is expanded in polynomial chaos associated with the two random variables C
0
ω
1, C
2
ω
1 of degree less than or equal to 2 and 5. The results predicted by the IRCM are better than those obtained by the chaos method. In addition, the number of unknowns used by the chaos method is 15 for polynomial chaos of degrees ≤ 2 and 126 for polynomial chaos of degree ≤ 5 while by the IRCM combined to the superposition method we have to solve a system with two differential equations and each one with 6 unknowns for expansion into polynomials of degrees ≤ 2 and 21 unknowns for expansion into the polynomial of degree ≤ 7 as presented in Table 5.Table 5
Number of coefficients to compute for IRCM and PC at degrees 2 and 5.
Method
IRCM
IRCM
Chaos
Chaos
Degree
2
5
2
5
Unknowns
12
42
15
126To analyze the standard deviation effect on the obtained results, the parameters of the system and of the excitation are random variables with 20% standard deviations. The obtained results are given in Figures6(a) and 6(b) for mean and standard deviation. It is shown that chaos of degree ≤ 2 does not predict good results in comparison with reference results, obtained by conditional expectation and Monte-Carlo 2000 methods. The IRCM with degree ≤ 2 predicts better results than chaos with degree ≤ 5 and IRCM with degree ≤ 5 coincides perfectly with the reference results. The random coefficients ξ
0, ξ
1, ξ
2, and ξ
3 are assumed here to follow different laws. The required CPU time for the used methods is presented in Table 6. These results demonstrate clearly the efficiency and accuracy of the proposed IRCM.Table 6
CPU times to compute for IRCM, PC at degrees 12, MC (2000), and conditional method.
Method
IRCM d° 12
PC d° 12
Conditional
MC (2000)
CPUT
6.52
1.47
e
+
002
2.98
e
+
003
36.43Figure 6
(a) Standard deviation of solution whereξ
0 and ξ
2 are truncated normal variables in [
-
4
,
4
], ξ
1 is a uniform variable in [
-
1,1
], and ξ
3 is a truncated exponential variable in [
0,10
] and σ
i = 20%. (b) Means of solution where ξ
0 and ξ
2 are truncated normal variables in [
-
4,4
], ξ
1 is a uniform variable in [
-
1,1
], and ξ
3 is a truncated exponential variable in [
0,10
] and σ
i = 20%.
(a)
(b)The error between the used methods is clearly demonstrated in Figures7, 8, and 9. It is clearly seen that the error is lower for IRCM than for PC related to the conditional solution assumed to be the reference solution.Figure 7
Absolute error between PC of degree 5 and IRCM of degree 5 and reference methods whereξ
0 and ξ
2 are truncated normal variables in [
-
4,4
], ξ
1 is a uniform variable in [
-
1,1
], and ξ
3 is a truncated exponential variable in [
0,10
] and σ
i = 20%.Figure 8
Absolute error between PC of degree 12 and IRCM of degree 12 and reference methods whereξ
0 and ξ
2 are truncated normal variables in [
-
4,4
], ξ
1 is a uniform variable in [
-
1,1
] and ξ
3 is a truncated exponential variable in [
0,10
] and σ
i = 20%.Figure 9
Absolute error between PC of degree 2 and IRCM of degree 2 MC (2000) and reference methods whereξ
0 and ξ
2 are truncated normal variables in [
-
4,4
], ξ
1 is a uniform variable in [
-
1,1
], and ξ
3 is a truncated exponential variable in [
0,10
] and σ
i = 5%.


## 4.1. Free Dynamics Response

For numerical analysis, the following values are used:(111)
ρ
0
=
7800
,
ρ
1
=
α
0
ρ
0
,
h
0
=
5
∗
10
-
5
,
h
1
=
α
1
h
0
,
E
0
=
21
∗
10
10
,
E
1
=
α
2
E
0
,
α
=
0.01
,
β
=
0.001
.The random variablesξ
0, ξ
1, ξ
2 are assumed to follow various laws and α
0, α
1, and α
2 are the ratios of the standard deviation and the mean of the random variables ρ, h, and E.Figures1(a) and 1(b) present the standard deviation and mean of the random solution U, where ξ
i are random variables with normalized truncated exponential laws in [
0,10
], with 5% of standard deviation for parameters (
α
0
=
α
1
=
α
2
=
5
%). The solution U
(
t
) is computed using the chaos methods associated with the random variables ξ
i, the IRCM1, IRCM2, the conditional expectation method, and the Monte-Carlo method with 1000 simples.Figure 1
(a) Standard deviation ofU when ξ
i follow the truncated exponential law in [
0,10
] and σ
i = 5%. (b) Means of U when ξ
i follow the truncated exponential law in [
0,10
] and σ
i = 5%.
(a)
(b)In the case of the classical chaos method, solutionU is expanded in the polynomial chaos of degree less than or equal to 7 and these polynomials are dependent on the three random variables ξ
i for i
=
0, 1, and 2. For the first internal random coefficient method, IRCM1, the coefficients of the random equation (97) are divided by the random variables C
0
(
ω
1
) and the random equation depends only on one new random variable η
=
C
2
ω
1
/
C
0
ω
1. The solution U is then expanded into the polynomials chaos associated with the random variable η of degree less than or equal to 12. For the second internal random coefficient method (IRCM2), the random equation depends on the random variables C
0
ω
1
,
C
2
ω
2. Solution Uis here expanded into the polynomials chaos associated with these two random variables of degree less than or equal to 2 and also less than 7. The number of unknowns used by the classical chaos method is 120 variables for degrees ≤ 7. The ICRM2 needs 6 unknowns for degrees ≤ 2 and 36 unknowns for degrees ≤ 7 while the IRCM1 needs only 13 unknowns for degrees ≤ 12 as presented in Table 2. This simple case demonstrates clearly how the IRCM reduces the number of unknowns.Table 2
Number of coefficients to compute for IRCM1-IRCM2 and PC at the considered degrees.
Method
IRCM1
IRCM2
IRCM2
Chaos
Degree
12
2
7
7
Unknowns
13
6
36
120The mean and standard deviation of the solutionU obtained by IRCM1, IRCM2, chaos, conditional expectation, and Monte-Carlo methods, presented in Figures 1(a) and 1(b), demonstrate a good accuracy between the predictions obtained by these methods.Figures2(a) and 2(b) present the standard deviation and mean of the solution U, when the random variables ξ
i are uniform in [
-
4,4
] with 5% standard deviation for each parameter (
α
0
=
α
1
=
α
2
=
5
%
). Solution U is computed using the previous methods. For the case of chaos method, the solution U is expanded into polynomials chaos of degree less than or equal to seven and these polynomials are dependent on the three random variables ξ
i for i
=
0, 1, and 2. For the ICRM1, solution U is expanded into the polynomial chaos associated with the random variable η of degree less than or equal to 12. For ICRM2, solution U is expanded into the polynomial chaos associated with two random variables of degree less than or equal to 7.Figure 2
(a) Standard deviation ofU when ξ
i follow the uniform law in [
-
4,4
] and σ
i = 5%. (b) Means of U when ξ
i follow the uniform law in [
-
4,4
] and σ
i = 5%.
(a)
(b)For large standard deviations, these methods are also tested. Figures3(a) and 3(b) present the standard deviation and mean of the solution U, when the random variables ξ
i are uniform in [-
4,4
] with 20% standard deviation (α
0
=
α
1
=
α
2
=
20
%). For the case of chaos method, solution U is expanded in polynomial chaos of degree less than or equal to 7 and these polynomials are dependent on the three random variables ξ
i for i
=
0, 1, and 2. ICRM1 and ICRM2 are also used with polynomial of degree ≤ 7. The number of coefficients to compute is presented in Table 3. Results obtained by the conditional expectation and Monte-Carlo methods are also presented. It is observed that, even for large standard deviation, 20%, a good agreement is observed between the predictions obtained by these methods.Table 3
Number of coefficients to compute for IRCM1-IRCM2 and PC at degree 7 for the considered number of variables.
Method
IRCM1
IRCM2
Chaos
Degree
7
7
7
Unknowns
8
36
120Figure 3
(a) Standard deviation ofU when ξ
i follow the uniform law in [
-
4,4
] and σ
i = 20%. (b) Means of U when ξ
i follow the uniform law in [
-
4,4
] and σ
i = 20%.
(a)
(b)Figures4(a) and 4(b) present the standard deviation and mean of solution U, where the random variables ξ
i follow a truncated exponential law in [
0
,
10
] with 20% standard deviation (α
0
=
α
1
=
α
2
=
20
%). The degree of the used polynomials is ≤7 for chaos, ≤12 for IRCM1, and ≤2 and ≤7 for IRCM2 and the number of coefficients to compute is presented in Table 4. It is demonstrated that, for the mean value, all the methods agree with each other but there are some discrepancies for the standard deviation. The IRCM2 is better than IRCM1 and particularly for large standard deviations. The division procedure, used in IRCM1, reduces the number of random parameters by one but introduces additional nonlinear random effect. This effect affects the accuracy of the IRCM1, as clearly presented in Figure 3.Table 4
Number of coefficients to compute for IRCM1-IRCM2 and PC at different degrees.
Method
IRCM1
IRCM2
IRCM2
Chaos
Degree
12
2
7
7
Unknowns
13
6
36
120Figure 4
(a) Standard deviation ofU when ξ
i follow the truncated exponential law in [
0,10
] and σ
i = 20%. (b) Means of U when ξ
i follow the truncated exponential law in [
0,10
] and σ
i = 20%.
(a)
(b)


## 4.2. Dynamical System with a Random Excitation

To test the flexibility of the presented methodological approach to handle dynamical systems with random variables following different laws, (97) is considered with a random excitation given by F
(
t
)
=
a
0
+
a
1
ξ
3
sin
⁡
(
ω
0
t
). a
0, a
1, and ω
0 are considered deterministic. The random coefficients ξ
0, ξ
1, ξ
2, and ξ
3 are assumed here to follow different laws.Figures5(a) and 5(b) present the mean and standard deviation of U
(
t
) obtained by the presented methods. ξ
0 and ξ
2 are truncated normal variables in [
-
4,4
], ξ
1 is a uniform variable in [
-
1,1
], and ξ
3 is a truncated exponential variable in [
0
,
10
] with 5% of standard deviations.Figure 5
(a) Standard deviation of solution whereξ
0 and ξ
2 are truncated normal variables in [
-
4,4
], ξ
1 is a uniform variable in [
-
1,1
], and ξ
3 is a truncated exponential variable in [
0,10
] and σ
i = 5%. (b) Means of solution where ξ
0 and ξ
2 are truncated normal variables in [
-
4,4
], ξ
1 is a uniform variable in [
-
1,1
], and ξ
3 is a truncated exponential variable in [
0,10
] and σ
i = 5%.
(a)
(b)The presented results are obtained by the chaos method where solutionU is expanded in the polynomials chaos of degree less than or equal to 2 and 5 and these polynomials depend on the four random variables ξ
i for i
=
0, 1, 2, and 3. For internal random coefficients method, solutionU is expanded in polynomial chaos associated with the two random variables C
0
ω
1, C
2
ω
1 of degree less than or equal to 2 and 5. The results predicted by the IRCM are better than those obtained by the chaos method. In addition, the number of unknowns used by the chaos method is 15 for polynomial chaos of degrees ≤ 2 and 126 for polynomial chaos of degree ≤ 5 while by the IRCM combined to the superposition method we have to solve a system with two differential equations and each one with 6 unknowns for expansion into polynomials of degrees ≤ 2 and 21 unknowns for expansion into the polynomial of degree ≤ 7 as presented in Table 5.Table 5
Number of coefficients to compute for IRCM and PC at degrees 2 and 5.
Method
IRCM
IRCM
Chaos
Chaos
Degree
2
5
2
5
Unknowns
12
42
15
126To analyze the standard deviation effect on the obtained results, the parameters of the system and of the excitation are random variables with 20% standard deviations. The obtained results are given in Figures6(a) and 6(b) for mean and standard deviation. It is shown that chaos of degree ≤ 2 does not predict good results in comparison with reference results, obtained by conditional expectation and Monte-Carlo 2000 methods. The IRCM with degree ≤ 2 predicts better results than chaos with degree ≤ 5 and IRCM with degree ≤ 5 coincides perfectly with the reference results. The random coefficients ξ
0, ξ
1, ξ
2, and ξ
3 are assumed here to follow different laws. The required CPU time for the used methods is presented in Table 6. These results demonstrate clearly the efficiency and accuracy of the proposed IRCM.Table 6
CPU times to compute for IRCM, PC at degrees 12, MC (2000), and conditional method.
Method
IRCM d° 12
PC d° 12
Conditional
MC (2000)
CPUT
6.52
1.47
e
+
002
2.98
e
+
003
36.43Figure 6
(a) Standard deviation of solution whereξ
0 and ξ
2 are truncated normal variables in [
-
4
,
4
], ξ
1 is a uniform variable in [
-
1,1
], and ξ
3 is a truncated exponential variable in [
0,10
] and σ
i = 20%. (b) Means of solution where ξ
0 and ξ
2 are truncated normal variables in [
-
4,4
], ξ
1 is a uniform variable in [
-
1,1
], and ξ
3 is a truncated exponential variable in [
0,10
] and σ
i = 20%.
(a)
(b)The error between the used methods is clearly demonstrated in Figures7, 8, and 9. It is clearly seen that the error is lower for IRCM than for PC related to the conditional solution assumed to be the reference solution.Figure 7
Absolute error between PC of degree 5 and IRCM of degree 5 and reference methods whereξ
0 and ξ
2 are truncated normal variables in [
-
4,4
], ξ
1 is a uniform variable in [
-
1,1
], and ξ
3 is a truncated exponential variable in [
0,10
] and σ
i = 20%.Figure 8
Absolute error between PC of degree 12 and IRCM of degree 12 and reference methods whereξ
0 and ξ
2 are truncated normal variables in [
-
4,4
], ξ
1 is a uniform variable in [
-
1,1
] and ξ
3 is a truncated exponential variable in [
0,10
] and σ
i = 20%.Figure 9
Absolute error between PC of degree 2 and IRCM of degree 2 MC (2000) and reference methods whereξ
0 and ξ
2 are truncated normal variables in [
-
4,4
], ξ
1 is a uniform variable in [
-
1,1
], and ξ
3 is a truncated exponential variable in [
0,10
] and σ
i = 5%.


## 5. Conclusion

Methodological approaches based on generalized polynomial chaos for random differential equation with an arbitrary number of random parameters following different types of distribution laws are developed. These random parameters can be linear or nonlinear. The conditional expectation method for the considered random differential equations is well clarified. The obtained results based on this method are considered to be exact results and used as reference results as well as those obtained by the Monte-Carlo method.A procedure to build a new polynomial chaos basis and a connection between the one-dimensional and multidimensional polynomials is well-established. Firstly, this procedure is given for independent variables and then generalized for variables that may be dependent. This procedure permits considering a basis given by the general polynomial chaos associated with the considered parameters and the solution is expanded into this new basis. A new methodological approach called internal random coefficient method is developed. Compact matrix representations are elaborated and the required matrices, vectors, and scalar products are explicitly given. These projection methods, combined with the superposition method, lead to large CPU time and memory space reductions required for the solution. It is demonstrated that the developed methods are efficient and accurate and these methods may lead to spectacular memory space and CPU time reductions and particularly when a large number of random parameters are considered.



---

*Source: 7403745-2018-05-21.xml*