# Generalized Neumann Expansion and Its Application in Stochastic Finite Element Methods

Authors: Xiangyu Wang; Song Cen; Chenfeng Li
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

An acceleration technique, termed generalized Neumann expansion (GNE), is presented for evaluating the responses of uncertain systems. The GNE method, which solves stochastic linear algebraic equations arising in stochastic finite element analysis, is easy to implement and is of high efficiency. The convergence condition of the new method is studied, and a rigorous error estimator is proposed to evaluate the upper bound of the relative error of a given GNE solution. It is found that the third-order GNE solution is sufficient to achieve a good accuracy even when the variation of the source stochastic field is relatively high. The relationship between the GNE method, the perturbation method, and the standard Neumann expansion method is also discussed. Based on the links between these three methods, quantitative error estimations for the perturbation method and the standard Neumann method are obtained for the first time in the probability context.


---

## Body

## 1. Introduction

Over the past few decades, the stochastic finite element method (SFEM) has attracted increasing attention from both academia and industries. The objective of SFEM is to analyze uncertainty propagation in practical engineering systems that are affected by various random factors, for example, heterogeneous materials and seismic loading. The current research of SFEM has two main focuses: probabilistic modeling of random media and efficient numerical solution of stochastic equations. Properties of heterogeneous materials (e.g., rocks and composites) vary randomly through the medium and are generally modeled with stochastic fields. Various probabilistic methods [1–7] have been developed to quantitatively describe the random material properties, which include the spectral representation method [1, 2], the Karhunen-Loève expansion method [3, 4], and the Fourier-Karhunen-Loève expansion method [5, 6], to name a few. In these methods, the random material properties are represented by a linear expansion such that they can be directly inserted into the governing equations of the physical system, which facilitates further numerical analysis. The resulting equations are stochastic equations, whose solution is the second main focus of SFEM research and also the focus of this work.For static and steady-state problems, the resulting SFEM equations have the following form [8–13]:
(1)(A0+A1α1+A2α2+⋯+Amαm)u(ω)=f(ω),
where A1,…,Am are n×n symmetric and deterministic matrices, u(ω) the unknown stochastic nodal vector, and f(ω) the nodal force vector. If the material properties are modeled as Gaussian fields, then α1,α2,…,αm are independent Gaussian random variables. The uncertainties from the right-hand side of (1) are relatively easy to handle because they are not directly coupled with the random matrix from the left-hand side.A number of methods have been developed for the solution of (1), including Monte Carlo methods, the perturbation method, the Neumann expansion method, and various projection methods. The Monte Carlo methods [14–19] first generate a number of samples of the random material properties, then solve each sample with the deterministic FEM solver, and finally evaluate the statistics from these deterministic solutions. The Monte Carlo methods are simple, flexible, friendly to parallel computing, and not restricted by the number of random variables [20]. The computational cost of Monte Carlo methods is generally high especially for large-scale problems, although there exist a number of efficient sampling techniques such as important sampling [14, 19], subset simulation [15, 17], and line sampling [18, 19]. Owning to its generosity and accuracy when sufficient sampling can be ensured, the Monte Carlo method is often taken as the reference method to verify other SFEM techniques. The perturbation method [11, 21–25] expresses the stiffness matrix, the load vector, and the unknown nodal displacement vector in terms of Taylor expansion with respect to the input random variables, and the Taylor expansion coefficients corresponding to the unknown nodal displacement are then solved by equating terms of the same order from both sides of (1). Its accuracy decreases when the random fluctuation increases, and there has not been a quantitative conclusion for the relationship between the solution error and the level of perturbation, making it difficult to confidently deploy the perturbation method in practical problems. The Neumann expansion method [9, 26, 27] is an iterative process for the solution of linear algebraic equations. Because the spectral radius of the iterative matrix must be smaller than one, it is often considered not suitable for problems with large random variations. However, the fluctuation range allowed by the Neumann expansion method is generally higher than the perturbation methods, because it is easier to adopt higher order expansions in the former one. It should be noted that the stationary iterative procedure employed in Neumann expansion is not as efficient as the preconditioned conjugate gradient procedure employed in Monte Carlo methods. Therefore, the computational efficiency of the Neumann expansion method is inferior to direct Monte Carlo simulations. In the projection methods [3, 4, 12, 28–35], the unknown nodal displacements are projected onto a set of orthogonal polynomials or other orthogonal bases, after which the Galerkin approach is employed to solve the projection coefficients. However, the projection methods typically result in an equation system much larger than the deterministic counterpart, and thus these methods are often limited to cases with a small set of base functions and a few independent random variables. In order to handle problems with more random variables, the stochastic collocation methods [36–42] have been developed, in which a reduced random basis is used in the projection operation and the projection coefficients are determined by approximating the precomputed Monte Carlo solutions. However, the stochastic collocation methods are believed to be less rigorous and their accuracy is poorer than the full projection methods using the Galerkin approach. Other variants of the projection procedures can be found in [13, 43, 44]. The fatal drawback of the projection types of solvers is their low efficiency [45], because the size of the base function set grows exponentially with the number of random variables and the order of the polynomials used. Moreover, the accuracy of the projection methods is criticized by [46, 47]. For linear stochastic problems, a joint diagonalization solution strategy was proposed in [8–10], in which the stochastic solution is obtained in a closed form by jointly diagonalizing the matrix family in (1). This method is not sensitive to the number of random variables in the system, while there has not been a mathematically rigorous proof for its accuracy.We propose a generalized Neumann expansion method for the solution of (1). The new GNE method is many times more efficient than the standard Neumann method, which is compromised by increased memory demand. In the context of probability, a mathematically rigorous error estimation is established for the new GNE scheme. It is also proved that for (1), the GNE scheme, the standard Neumann expansion method, and the perturbation method are equivalent. Hence, the stochastic error estimation obtained here also shed a light on the standard Neumann method and the perturbation method, which has remained as an outstanding challenge in the literature. The rest of the paper is organized as follows. The GNE scheme and its error estimation are derived in Section 2, after which it is validated and examined in Section 3 via a comprehensive numerical test. Concluding remarks are given in Section 4.


## 2. Generalized Neumann Expansion Method

Consider a linear elastic statics system with Young’s modulus modeled as a random field. By using the Karhunen-Loève expansion [3, 4] or the Fourier-Karhunen-Loève expansion [5, 6], random Young’s modulus can be expressed as
(2)E(x,ω)=E0(x)+∑iλiαi(ω)ψi(x),
where E(x,ω) denotes random Young’s modulus at point x, E0(x) is the expectation of E(x,ω), αi(ω) are uncorrelated random variables, and λi and ψi(x) are deterministic eigen values and eigen functions corresponding to the covariance function of the source field E(x,ω). If E(x,ω) is modeled as a Gaussian field, αi(ω) become independent standard Gaussian random variables. For the sake of simplicity, we drop the random symbol ω and for the rest of the paper use αi to denote the random variables in (2).Following the standard finite element (FE) discretization procedure, a stochastic system of linear algebraic equations can be built from (2)
(3)(A0+∑iαiAi)u=f.
In the above equation, A0 is the conventional FE stiffness matrix, which is symmetric, positive definite, and sparse while Ai are symmetric and sparse matrices. The deterministic matrices Ai share the same entry pattern as the stiffness matrix A0, but they are usually not positive definite.


### 2.1. Related Methods

The Monte Carlo methods first generate a set of samples ofαi in (3), and each sample of αi corresponds to a deterministic equation, which can be solved by the standard FEM method.In the perturbation procedure, the unknown random vectoru is expanded as [21, 25, 26]
(4)u=u0+∑iαiui+∑i∑jαiαjuij+⋯.
Substituting (4) into (3) yields
(5)(A0+∑iαiAi)(u0+∑iαiui+∑i∑jαiαjuij+⋯)=f,
from which the components u0, ui, uij of u are solved as follows:
(6)u0=A0-1f,ui=-A0-1Aiu0,uij=A0-1AiA0-1Aju0.
It should be noted that, in the stochastic context, there has not been a rigorous convergence criterion for the perturbation scheme described above.In the Neumann expansion method, the solution of (3) is written as
(7)u=(A0+ΔA)-1f=(I+B)-1A0-1f,
where ΔA=∑iαiAi, B=A0-1ΔA and I is the identity matrix. According to Neumann expansion
(8)(I+B)-1=I-B+B2-B3+⋯,
solution (7) can be simplified as
(9)u=u0-Bu0+B2u0-B3u0+⋯,
where u0=A0-1f. The spectral radius of matrix B must be smaller than 1 to ensure the convergence of (9).


### 2.2. The Generalized Neumann Expansion Method

The Neumann expansion (8) can be generalized to the inverse of multiple matrices:
(10)(I+∑iαiBi)-1=I-∑iαiBi+∑i∑jαiαjBiBj-⋯.
The generalized Neumann expansion (GNE) in (10) can be proved as follows:
(11)(I+∑iαiBi)(I+∑iαiBi)-1=I-∑iαiBi+∑i∑jαiαjBiBj-⋯+α1B1-∑iα1αiB1Bi+⋯+α2B2-∑iα2αiB2Bi+⋯⋮+αnBn-∑iαnαiBnBi+⋯=I,
where n is the total number of Bi.The Neumann expansion (8) and the generalized Neumann expansion (10) are equivalent in mathematics and only differ in organizing the matrix operations. This difference can lead to a significant acceleration in computing the solution of (3). Specifically, let Bi=A0-1Ai, and applying the GNE in (10) to (3) gives
(12)u=(A0+∑iαiAi)-1f=(I+∑iαiA0-1Ai)-1A0-1fu=u0-∑iαiA0-1Aiu0+∑i∑jαjαiA0-1AjA0-1Aiu0-⋯u=u0-∑iαiui+∑i∑jαiαjuij-⋯ in which
(13)ui=A0-1Aiu0,uij=A0-1AiA0-1Aju0.
Equation (12) is the GNE solution.


### 2.3. Convergence Analysis and Error Estimation

The convergence condition of the GNE (10) is essentially the same as the Neumann expansion (8), that is, ρ(∑αiBi)<1, where ρ(·) denotes the spectral radius of a matrix. Hence, for the GNE solution in (12), the convergence condition is
(14)ρ(∑αiΑ0-1Ai)<1.
It should be noted that
(15)ρ(∑αiΑ0-1Ai)=ρ(∑αiC-1AiC-1)=∥∑αiC-1AiC-1∥m2,
where ∥·∥m2 denotes the spectral norm of a matrix, and the symmetric and positive definite matrix C is defined as
(16)A0=C2.
As a type of matrix norm, the spectral norm has the following properties:
(17)∥P+Q∥m2≤∥P∥m2+∥Q∥m2,∥PQ∥m2≤∥P∥m2∥Q∥m2,∥hP∥m2=h∥P∥m2,
where P and Q are general matrices and h is a positive constant. Hence,
(18)∥∑αiC-1AiC-1∥m2≤∥α1C-1A1C-1∥m2+∥α2C-1A2C-1∥m2+⋯∥αnC-1AnC-1∥m2,ρ(∑αiC-1AiC-1)≤ρ(α1C-1A1C-1)+ρ(α2C-1A2C-1)+⋯ρ(αnC-1AnC-1),ρ(∑αiΑ0-1Ai)≤ρ(α1Α0-1A1)+ρ(α2Α0-1A2)+⋯ρ(αnΑ0-1An)≤nmax(αi)max(ρ(Α0-1Ai)).
Based on the above relation, a sufficient but not necessary condition for convergence is
(19)nmax(αi)max(ρ(Α0-1Ai))<1.
The maximum spectral radius ρmax=max(ρ(Α0-1Ai)) can be obtained from A0 and Ai. Equation (19) can be used to verify the convergence of GNE in advance. When αi are independent standard Gaussian random variables, their absolute values are smaller than 3 with the probability of 99.7%. In this case, (19) can be further simplified:
(20)nmax(αi)max(ρ(Α0-1Ai))≈3nρmax<1,ρmax<13n.
Considering the kth-order GNE (10) or Neumann expansion (8), their error can be estimated as follows:
(21)∥(I+B)-1-(I-B+B2+⋯+(-1)kBk)∥m2=∥Bk+1(I+B)-1∥m2≤∥Bk+1∥m2∥(I+B)-1∥m2=σBk+11-σB,
where σB is the spectral norm of B, which is assumed to satisfy
(22)σB<1.
The exact solution of (3) is
(23)u=(I+∑iαiA0-1Ai)-1A0-1f=(I+B)-1A0-1f,Cu=(I+∑iαiC-1AiC-1)-1C-1fCu=(I+D)-1C-1f,
where B=∑iαiA0-1Ai, D=∑iαiC-1AiC-1. The kth-order GNE or Neumann expansion solution is given by
(24)u~=(I-B+B2+⋯+(-1)kBk)A0-1f,Cu~=(I-D+D2+⋯+(-1)kDk)C-1f.
The norm of the solution error is defined as
(25)∥Cu-Cu~∥2=∥[(I+D)-1-(I-D+D2+⋯+(-1)kDk)]C-1f∥2=∥Dk+1(I+D)-1C-1f∥2=∥Dk+1Cu∥2≤∥Dk+1∥m2∥Cu∥2=σDk+1∥Cu∥2=ρBk+1∥Cu∥2,
where ρB is the spectral radius of B.The relative error of the solution is defined as(26)RE=(u-u~)TA0(u-u~)uTA0uRE=(u-u~)TC2(u-u~)uTC2u=∥Cu-Cu~∥2∥Cu∥2≤ρBk+1.
Comparing (21) and (26), it can be seen that the accuracy for the solution vector u is even better than the accuracy for the matrix inverse (I+B)-1.


### 2.4. Solutions for More General Stochastic Algebraic Equations

For the solution of (3), the GNE solution in (12) shares the same form as the perturbation solution in (5). Hence, for solving (3), the GNE method, the perturbation method, and the standard Neumann expansion method are identical. The error estimation obtained in (26) also holds for the perturbation approach.However, for the following more general stochastic algebraic equations [26]:
(27)(A0+∑iαiAi+∑i∑jαiαjAij+⋯)u=f0+∑iγifi+∑i∑jγiγjfij+⋯,
the three methods are not equivalent. Specifically, applying the perturbation method to the above equation yields
(28)(A0+∑iαiAi+∑i∑jαiαjAij+⋯)×(u0+∑iβiui+∑i∑jβiβjuij+⋯)=f0+∑iγifi+∑i∑jγiγjfij+⋯.
Equation (28) can be reorganized as(29)A0u0+∑iαiAiu0+∑i∑jαiαjAiju0⋯+∑iA0βiui+∑i∑jαiAiβjuj+∑i∑j∑kαiαjAijβkuk⋯+∑i∑jA0βiβjuij+∑i∑j∑kαiAiβjβkujk+∑i∑j∑k∑lαiαjAijβkβlukl⋯⋯⋯⋯⋱=f0+∑iγifi+∑i∑jγiγjfij+⋯.If αi=βi=γi, then the perturbation components u0,ui,uij,…,ui1i2⋯inof u can be solved as
(30)A0u0=f0,A0ui+Aiu0=fi,A0uij+Aiuj+Aiju0=fij,⋮A0ui1i2⋯in+Ai1ui2i3⋯in+Ai1i2ui3i4⋯in+Ai1i2⋯inu0=fi1i2⋯in.
Solving (27) with the Neumann expansion gives
(31)A-1=(A0+ΔA)-1A-1=(I+B)-1A0-1=(I-B+B2-B3+⋯)A0-1,
where B=A0-1ΔA,ΔA=∑iαiAi+∑i∑jαiαjAij+⋯. Then
(32)u=(I-B+B2-B3+⋯)A0-1(f0+Δf),
in which
(33)Δf=∑iγifi+∑i∑jγiγjfij+⋯.
Because higher order entries are included in B, the term BA0-1(f0+Δf) (the first-order term of the Neumann solution in (32)) does not correspond to the first-order entries of u. Therefore, the Neumann solution in (32) and the perturbation solution in (30) are different for the general stochastic equation (27).Solving (27) with the GNE approach gives
(34)u=(A0+∑iiαi1Ai1+∑i2∑i1αi1αi2Ai1i2+⋯+∑in⋯∑i2∑i1αi1i2⋯inAi1i2⋯in)-1fu=(I+∑iiαi1Bi1+∑i2∑i1αi1αi2Bi1i2+⋯+∑in⋯∑i2∑i1αi1i2⋯inBi1i2⋯in)-1A0-1fu=(I+∑iβiCi)-1A0-1fu=(I-∑iβiCi+∑i∑jβiβjCiCj-∑i∑j∑kβiβjβkCiCjCk+⋯)A0-1f,
where f=f0+∑iγifi+∑i∑jγiγjfij+⋯, Bi=A0-1Ai, Bij=A0-1Aij, and other B terms can be similarly formed. Matrices Ci are determined by realignment of Bi, Bij,… and random variables βi are formed by rearrangement of αi1,αi1i2,…,αi1i2⋯in.In summary, when solving (27), the GNE method and the Neumann method are identical, while they are not equivalent to the perturbation procedure. However, if the higher order terms in the GNE or Neumann expansion solutions are omitted, they degenerate into the perturbation solution. Two simplified examples are given here to demonstrate the difference between the GNE and the perturbation method.


#### 2.4.1. Example I

Only first-order entries are included in both sides of (27) such that the equation becomes
(35)(A0+∑iαiAi)u=f0+∑iαifi.
The first-order perturbation solution is
(36)u0=A0-1f0,ui=A0-1fi-A0-1Aiu0,u=u0+∑iαiuiu=A0-1f0+∑iαiA0-1fi-∑iαiA0-1Aiu0,
and the first-order GNE solution is
(37)u=(A0-1-∑iαiA0-1AiA0-1)(f0+∑iαifi)u=A0-1f0-∑iαiA0-1Aiu0+∑iαiA0-1fiu=-∑i∑jαiαjA0-1Aifj.
Comparing (36) and (37), it can be seen that the GNE solution contains an extra second-order term -∑i∑jαiαjA0-1Aifj.


#### 2.4.2. Example II

Up to second-order terms are included in the left-hand side of (27), while only the constant term is included in the right-hand side. In this case, the equation becomes
(38)(A0+∑iαiAi+∑i∑jαiαjAij)u=f0.
Denote u0=A0-1f0, Bi=A0-1Ai, and Bij=A0-1Aij; then the first-order perturbation solution can be expressed as
(39)u0=A0-1f0,ui=-∑iA0-1Aiu0=-∑iBiu0,u=u0+∑iαiui=u0-∑iαiBiu0,
while the first-order GNE solution is
(40)u=(I+∑iαiA0-1Ai+∑i∑jαiαjA0-1Aij)-1u0u=(I+∑iαiBi+∑i∑jαiαjBij)-1u0u=(I-∑iαiBi-∑i∑jαiαjBij)u0.
Comparing (39) and (40), the second-order terms -∑i∑jαiαjA0-1Aiju0 in the GNE solution do not appear in the perturbation solution (39).


## 2.1. Related Methods

The Monte Carlo methods first generate a set of samples ofαi in (3), and each sample of αi corresponds to a deterministic equation, which can be solved by the standard FEM method.In the perturbation procedure, the unknown random vectoru is expanded as [21, 25, 26]
(4)u=u0+∑iαiui+∑i∑jαiαjuij+⋯.
Substituting (4) into (3) yields
(5)(A0+∑iαiAi)(u0+∑iαiui+∑i∑jαiαjuij+⋯)=f,
from which the components u0, ui, uij of u are solved as follows:
(6)u0=A0-1f,ui=-A0-1Aiu0,uij=A0-1AiA0-1Aju0.
It should be noted that, in the stochastic context, there has not been a rigorous convergence criterion for the perturbation scheme described above.In the Neumann expansion method, the solution of (3) is written as
(7)u=(A0+ΔA)-1f=(I+B)-1A0-1f,
where ΔA=∑iαiAi, B=A0-1ΔA and I is the identity matrix. According to Neumann expansion
(8)(I+B)-1=I-B+B2-B3+⋯,
solution (7) can be simplified as
(9)u=u0-Bu0+B2u0-B3u0+⋯,
where u0=A0-1f. The spectral radius of matrix B must be smaller than 1 to ensure the convergence of (9).


## 2.2. The Generalized Neumann Expansion Method

The Neumann expansion (8) can be generalized to the inverse of multiple matrices:
(10)(I+∑iαiBi)-1=I-∑iαiBi+∑i∑jαiαjBiBj-⋯.
The generalized Neumann expansion (GNE) in (10) can be proved as follows:
(11)(I+∑iαiBi)(I+∑iαiBi)-1=I-∑iαiBi+∑i∑jαiαjBiBj-⋯+α1B1-∑iα1αiB1Bi+⋯+α2B2-∑iα2αiB2Bi+⋯⋮+αnBn-∑iαnαiBnBi+⋯=I,
where n is the total number of Bi.The Neumann expansion (8) and the generalized Neumann expansion (10) are equivalent in mathematics and only differ in organizing the matrix operations. This difference can lead to a significant acceleration in computing the solution of (3). Specifically, let Bi=A0-1Ai, and applying the GNE in (10) to (3) gives
(12)u=(A0+∑iαiAi)-1f=(I+∑iαiA0-1Ai)-1A0-1fu=u0-∑iαiA0-1Aiu0+∑i∑jαjαiA0-1AjA0-1Aiu0-⋯u=u0-∑iαiui+∑i∑jαiαjuij-⋯ in which
(13)ui=A0-1Aiu0,uij=A0-1AiA0-1Aju0.
Equation (12) is the GNE solution.


## 2.3. Convergence Analysis and Error Estimation

The convergence condition of the GNE (10) is essentially the same as the Neumann expansion (8), that is, ρ(∑αiBi)<1, where ρ(·) denotes the spectral radius of a matrix. Hence, for the GNE solution in (12), the convergence condition is
(14)ρ(∑αiΑ0-1Ai)<1.
It should be noted that
(15)ρ(∑αiΑ0-1Ai)=ρ(∑αiC-1AiC-1)=∥∑αiC-1AiC-1∥m2,
where ∥·∥m2 denotes the spectral norm of a matrix, and the symmetric and positive definite matrix C is defined as
(16)A0=C2.
As a type of matrix norm, the spectral norm has the following properties:
(17)∥P+Q∥m2≤∥P∥m2+∥Q∥m2,∥PQ∥m2≤∥P∥m2∥Q∥m2,∥hP∥m2=h∥P∥m2,
where P and Q are general matrices and h is a positive constant. Hence,
(18)∥∑αiC-1AiC-1∥m2≤∥α1C-1A1C-1∥m2+∥α2C-1A2C-1∥m2+⋯∥αnC-1AnC-1∥m2,ρ(∑αiC-1AiC-1)≤ρ(α1C-1A1C-1)+ρ(α2C-1A2C-1)+⋯ρ(αnC-1AnC-1),ρ(∑αiΑ0-1Ai)≤ρ(α1Α0-1A1)+ρ(α2Α0-1A2)+⋯ρ(αnΑ0-1An)≤nmax(αi)max(ρ(Α0-1Ai)).
Based on the above relation, a sufficient but not necessary condition for convergence is
(19)nmax(αi)max(ρ(Α0-1Ai))<1.
The maximum spectral radius ρmax=max(ρ(Α0-1Ai)) can be obtained from A0 and Ai. Equation (19) can be used to verify the convergence of GNE in advance. When αi are independent standard Gaussian random variables, their absolute values are smaller than 3 with the probability of 99.7%. In this case, (19) can be further simplified:
(20)nmax(αi)max(ρ(Α0-1Ai))≈3nρmax<1,ρmax<13n.
Considering the kth-order GNE (10) or Neumann expansion (8), their error can be estimated as follows:
(21)∥(I+B)-1-(I-B+B2+⋯+(-1)kBk)∥m2=∥Bk+1(I+B)-1∥m2≤∥Bk+1∥m2∥(I+B)-1∥m2=σBk+11-σB,
where σB is the spectral norm of B, which is assumed to satisfy
(22)σB<1.
The exact solution of (3) is
(23)u=(I+∑iαiA0-1Ai)-1A0-1f=(I+B)-1A0-1f,Cu=(I+∑iαiC-1AiC-1)-1C-1fCu=(I+D)-1C-1f,
where B=∑iαiA0-1Ai, D=∑iαiC-1AiC-1. The kth-order GNE or Neumann expansion solution is given by
(24)u~=(I-B+B2+⋯+(-1)kBk)A0-1f,Cu~=(I-D+D2+⋯+(-1)kDk)C-1f.
The norm of the solution error is defined as
(25)∥Cu-Cu~∥2=∥[(I+D)-1-(I-D+D2+⋯+(-1)kDk)]C-1f∥2=∥Dk+1(I+D)-1C-1f∥2=∥Dk+1Cu∥2≤∥Dk+1∥m2∥Cu∥2=σDk+1∥Cu∥2=ρBk+1∥Cu∥2,
where ρB is the spectral radius of B.The relative error of the solution is defined as(26)RE=(u-u~)TA0(u-u~)uTA0uRE=(u-u~)TC2(u-u~)uTC2u=∥Cu-Cu~∥2∥Cu∥2≤ρBk+1.
Comparing (21) and (26), it can be seen that the accuracy for the solution vector u is even better than the accuracy for the matrix inverse (I+B)-1.


## 2.4. Solutions for More General Stochastic Algebraic Equations

For the solution of (3), the GNE solution in (12) shares the same form as the perturbation solution in (5). Hence, for solving (3), the GNE method, the perturbation method, and the standard Neumann expansion method are identical. The error estimation obtained in (26) also holds for the perturbation approach.However, for the following more general stochastic algebraic equations [26]:
(27)(A0+∑iαiAi+∑i∑jαiαjAij+⋯)u=f0+∑iγifi+∑i∑jγiγjfij+⋯,
the three methods are not equivalent. Specifically, applying the perturbation method to the above equation yields
(28)(A0+∑iαiAi+∑i∑jαiαjAij+⋯)×(u0+∑iβiui+∑i∑jβiβjuij+⋯)=f0+∑iγifi+∑i∑jγiγjfij+⋯.
Equation (28) can be reorganized as(29)A0u0+∑iαiAiu0+∑i∑jαiαjAiju0⋯+∑iA0βiui+∑i∑jαiAiβjuj+∑i∑j∑kαiαjAijβkuk⋯+∑i∑jA0βiβjuij+∑i∑j∑kαiAiβjβkujk+∑i∑j∑k∑lαiαjAijβkβlukl⋯⋯⋯⋯⋱=f0+∑iγifi+∑i∑jγiγjfij+⋯.If αi=βi=γi, then the perturbation components u0,ui,uij,…,ui1i2⋯inof u can be solved as
(30)A0u0=f0,A0ui+Aiu0=fi,A0uij+Aiuj+Aiju0=fij,⋮A0ui1i2⋯in+Ai1ui2i3⋯in+Ai1i2ui3i4⋯in+Ai1i2⋯inu0=fi1i2⋯in.
Solving (27) with the Neumann expansion gives
(31)A-1=(A0+ΔA)-1A-1=(I+B)-1A0-1=(I-B+B2-B3+⋯)A0-1,
where B=A0-1ΔA,ΔA=∑iαiAi+∑i∑jαiαjAij+⋯. Then
(32)u=(I-B+B2-B3+⋯)A0-1(f0+Δf),
in which
(33)Δf=∑iγifi+∑i∑jγiγjfij+⋯.
Because higher order entries are included in B, the term BA0-1(f0+Δf) (the first-order term of the Neumann solution in (32)) does not correspond to the first-order entries of u. Therefore, the Neumann solution in (32) and the perturbation solution in (30) are different for the general stochastic equation (27).Solving (27) with the GNE approach gives
(34)u=(A0+∑iiαi1Ai1+∑i2∑i1αi1αi2Ai1i2+⋯+∑in⋯∑i2∑i1αi1i2⋯inAi1i2⋯in)-1fu=(I+∑iiαi1Bi1+∑i2∑i1αi1αi2Bi1i2+⋯+∑in⋯∑i2∑i1αi1i2⋯inBi1i2⋯in)-1A0-1fu=(I+∑iβiCi)-1A0-1fu=(I-∑iβiCi+∑i∑jβiβjCiCj-∑i∑j∑kβiβjβkCiCjCk+⋯)A0-1f,
where f=f0+∑iγifi+∑i∑jγiγjfij+⋯, Bi=A0-1Ai, Bij=A0-1Aij, and other B terms can be similarly formed. Matrices Ci are determined by realignment of Bi, Bij,… and random variables βi are formed by rearrangement of αi1,αi1i2,…,αi1i2⋯in.In summary, when solving (27), the GNE method and the Neumann method are identical, while they are not equivalent to the perturbation procedure. However, if the higher order terms in the GNE or Neumann expansion solutions are omitted, they degenerate into the perturbation solution. Two simplified examples are given here to demonstrate the difference between the GNE and the perturbation method.


### 2.4.1. Example I

Only first-order entries are included in both sides of (27) such that the equation becomes
(35)(A0+∑iαiAi)u=f0+∑iαifi.
The first-order perturbation solution is
(36)u0=A0-1f0,ui=A0-1fi-A0-1Aiu0,u=u0+∑iαiuiu=A0-1f0+∑iαiA0-1fi-∑iαiA0-1Aiu0,
and the first-order GNE solution is
(37)u=(A0-1-∑iαiA0-1AiA0-1)(f0+∑iαifi)u=A0-1f0-∑iαiA0-1Aiu0+∑iαiA0-1fiu=-∑i∑jαiαjA0-1Aifj.
Comparing (36) and (37), it can be seen that the GNE solution contains an extra second-order term -∑i∑jαiαjA0-1Aifj.


### 2.4.2. Example II

Up to second-order terms are included in the left-hand side of (27), while only the constant term is included in the right-hand side. In this case, the equation becomes
(38)(A0+∑iαiAi+∑i∑jαiαjAij)u=f0.
Denote u0=A0-1f0, Bi=A0-1Ai, and Bij=A0-1Aij; then the first-order perturbation solution can be expressed as
(39)u0=A0-1f0,ui=-∑iA0-1Aiu0=-∑iBiu0,u=u0+∑iαiui=u0-∑iαiBiu0,
while the first-order GNE solution is
(40)u=(I+∑iαiA0-1Ai+∑i∑jαiαjA0-1Aij)-1u0u=(I+∑iαiBi+∑i∑jαiαjBij)-1u0u=(I-∑iαiBi-∑i∑jαiαjBij)u0.
Comparing (39) and (40), the second-order terms -∑i∑jαiαjA0-1Aiju0 in the GNE solution do not appear in the perturbation solution (39).


## 2.4.1. Example I

Only first-order entries are included in both sides of (27) such that the equation becomes
(35)(A0+∑iαiAi)u=f0+∑iαifi.
The first-order perturbation solution is
(36)u0=A0-1f0,ui=A0-1fi-A0-1Aiu0,u=u0+∑iαiuiu=A0-1f0+∑iαiA0-1fi-∑iαiA0-1Aiu0,
and the first-order GNE solution is
(37)u=(A0-1-∑iαiA0-1AiA0-1)(f0+∑iαifi)u=A0-1f0-∑iαiA0-1Aiu0+∑iαiA0-1fiu=-∑i∑jαiαjA0-1Aifj.
Comparing (36) and (37), it can be seen that the GNE solution contains an extra second-order term -∑i∑jαiαjA0-1Aifj.


## 2.4.2. Example II

Up to second-order terms are included in the left-hand side of (27), while only the constant term is included in the right-hand side. In this case, the equation becomes
(38)(A0+∑iαiAi+∑i∑jαiαjAij)u=f0.
Denote u0=A0-1f0, Bi=A0-1Ai, and Bij=A0-1Aij; then the first-order perturbation solution can be expressed as
(39)u0=A0-1f0,ui=-∑iA0-1Aiu0=-∑iBiu0,u=u0+∑iαiui=u0-∑iαiBiu0,
while the first-order GNE solution is
(40)u=(I+∑iαiA0-1Ai+∑i∑jαiαjA0-1Aij)-1u0u=(I+∑iαiBi+∑i∑jαiαjBij)-1u0u=(I-∑iαiBi-∑i∑jαiαjBij)u0.
Comparing (39) and (40), the second-order terms -∑i∑jαiαjA0-1Aiju0 in the GNE solution do not appear in the perturbation solution (39).


## 3. Numerical Examples

### 3.1. Example Description

As shown in Figure1, a plain strain problem with an orthohexagonal cross section is considered here. The length of each edge of the hexagon is 0.4 m. The Poisson ratio of the material is set as 0.2, and Young’s modulus is assumed to be a Gaussian field with a mean value of 30 GPa and a covariance function defined as
(41)R(τ1,τ2)=c1e-τ12/c2-τ22/c2,
in which τ1=x1-x2 and τ2=y1-y2 are the coordinate differences and c1 and c2 are adjustable constants. As shown in Figure 2, the cross section is discretized with triangular linear elements, where each side is divided into 8 equal parts. There are in total 384 elements and 217 nodes.Figure 1
Cross section of the random material sample.Figure 2
Finite element mesh.The relative error is defined as(42)RE=∥CuGNE-CuMC∥2∥CuMC∥2×100%,
in which uMC is the nodal displacement vector solved with the Monte Carlo method and uGNE is the GNE solution. Unless stated explicitly, the numerical testing is performed for 4340 samples generated from the Gaussian field defined in (41).The convergence condition of the GNE method has been rigorously proved. The main purpose of this example is to examine the error distribution and its relationship with various random parameters, for which a number of numerical tests are performed in the following subsections.


### 3.2. Error Estimation

The third-order GNE solution is adopted. The constantsc1 and c2 in (41) are varied to examine the error estimator in (26). The constant c1 is the variance of the Gaussian field, and the constant c2 controls the shape of the covariance function. The effective correlation length Leff of the covariance function R(τ) is defined by
(43)R(Leff)R(0)=0.01.
For any two points in the medium, if the distance between them is larger than Leff, their material properties are only weakly correlated. The effective correlation length of the Gaussian field defined in (41) is
(44)R(τ)R(0)=e-τ2/c2=0.01⟹Leff=2.146c2.
In the first set of tests, the constant c1 is varied to change the variance of the Gaussian field. The parameter settings are
(45)c1=0.01,c2=3.2E-2,Leff=0.38,c1=0.04,c2=3.2E-2,Leff=0.38.
In the second set of tests, the constant c2 is varied to change the shape of the covariance function. The parameter settings are
(46)c1=0.01,c2=0.8E-2,Leff=0.19,c1=0.01,c2=2.4E-2,Leff=0.33,c1=0.01,c2=4.0E-2,Leff=0.43.


#### 3.2.1. The Influence ofc1

Corresponding toc1=0.01 and c1=0.04, the standard deviation σ of the stochastic field is 0.1 and 0.2, respectively. A total number of 4340 samples are generated for each case. For each realization, both the Monte Carlo solution and the third-order GNE solution are computed. The relationship between the spectral radius and the relative error is drawn in Figures 3 and 4, in which the red crosses show the actual relative errors corresponding to each sample and the blue line shows the theoretical error estimation. It can been seen that for both cases, the theoretical error estimator in (26) holds accurately and its effectiveness is not affected by the change of the variance controlled by c1.Figure 3
The relation between the relative error and the spectral radius.Figure 4
The relation between the relative error and the spectral radius.


#### 3.2.2. The Influence ofc2

The constantc2 is adjusted to take three different values 0.8×10-2, 2.4×10-2, and 4.0×10-2, while the constant c1 remains unchanged. Following (41), the effective correlation length grows with the increase of c2. Again, 4340 samples are generated for each stochastic field, and for each realization both the Monte Carlo solution and the third-order GNE solution are computed. Figure 5 plots the relationship between spectral radius and the relative error, which verifies that, for all cases, the theoretical error estimation in (26) holds. The changing of c2, which changes the variation frequency of the stochastic field, does not affect the error estimation.Figure 5
The relation between the relative error and the spectral radius.As indicated in (44), the effective correlation length Leff is directly associated with the parameter c2. To investigate the influence of Leff on the relative error, the empirical cumulative distribution functions (ECDF) of relative errors are plotted in Figure 6 for different Leff cases. Comparing the three ECDF curves at different effective correlation lengths, it is found that with the decrease of the effective correlation length, the accuracy of the GNE solution improves.Figure 6
ECDF of the relative errors.


### 3.3. The Relative Error for Other Types of Stochastic Fields

In the last two subsections, the error estimation analysis has taken the Gaussian assumption; that is, all random variablesαi in (3) are assumed to be independent Gaussian random variables. In this subsection, we investigate the relative error of the GNE method for other types of stochastic fields. Specifically, the random variables αi in (3) are assumed to have the two-point distribution, the binomial distribution, and the β distribution, respectively. For the case of two-point distribution, values {-0.3E0,0.3E0} are assigned with equal probability, that is, 0.5. For the case of binomial distribution, ten times Bernoulli trials with the probability 0.5 are adopted, and the integer values between 0 and 10 are linearly projected onto [-0.3E0,0.3E0]. For the case of β distribution, the probability density function is defined as
(47)f(x∣a,b)=1B(a,b)xa-1(1-x)b-1I(0,1)(x),
where B(a,b) is a β function, and the filter function I(0,1)(x) is defined as
(48)I(0,1)(x)={1x∈(0,1)0x∉(0,1).
Three kinds of β distributions are considered, with a=b=0.75, a=b=1, and a=b=4, respectively. The interval (0,1) is linearly mapped to [-0.3E0,0.3E0].Corresponding to these five distributions, the empirical cumulative distribution functions of relative errors are plotted in Figure7. Comparing the ECDF curves at different effective correlation lengths, the same trend is observed as the Gaussian case. That is, the relative error decreases with the drop of effective correlation length.ECDF of the relative error.
(a)
Two-point distribution
(b)
Binomial distribution
(c)
β distribution with a=b=0.75
(d)
β distribution with a=b=1
(e)
β distribution with a=b=4


### 3.4. Efficiency Analysis

Using the same computing platform, the Monte Carlo method, the Neumann expansion method, and the GNE method are employed to solve the same set of solutions (corresponding to 4340 samples), and the computation time is shown in Table1. The number of random variables in the stochastic equation system differs depending on the effective correlation length. A fixed mesh is adopted in all of the computation such that the dimension of a single equation system is 434×434.Table 1
Computation time for 4340 samples (unit: s).
Number of random variables
10
11
16
21
39
Monte Carlo
40.78
43.26
54.02
65.55
107.13
GNE
First order
0.44
0.44
0.48
0.53
0.69
Second order
0.79
0.83
1.05
1.28
2.21
Third order
1.30
1.43
2.29
3.58
13.33
Neumann
First order
26.74
28.96
40.43
52.23
93.90
Second order
28.63
30.88
42.28
53.70
95.22
Third order
30.28
32.57
44.05
55.49
97.08If the mesh density is doubled, then the dimension of each matrix is1634×1634. The number of random variables remains the same in each case. The corresponding computation time is shown in Table 2.Table 2
Computation time for 16340 samples (unit: s).
Number of random variables
11
16
21
Monte Carlo
6560.01
7871.31
9118.39
GNEM
First order
13.31
14.85
16.39
Second order
16.06
19.23
22.50
Third order
24.55
41.86
70.97
Neumann
First order
3325.30
4851.50
6017.91
Second order
3639.38
5010.79
6314.89
Third order
3881.43
5473.71
6609.92The two data tables indicate that the computational efficiency of the GNE method (up to the third order) is much higher than Neumann expansion and Monte Carlo procedures. However, the memory requirement for higher order GNE schemes increases rapidly causing the computational efficiency to decrease.


## 3.1. Example Description

As shown in Figure1, a plain strain problem with an orthohexagonal cross section is considered here. The length of each edge of the hexagon is 0.4 m. The Poisson ratio of the material is set as 0.2, and Young’s modulus is assumed to be a Gaussian field with a mean value of 30 GPa and a covariance function defined as
(41)R(τ1,τ2)=c1e-τ12/c2-τ22/c2,
in which τ1=x1-x2 and τ2=y1-y2 are the coordinate differences and c1 and c2 are adjustable constants. As shown in Figure 2, the cross section is discretized with triangular linear elements, where each side is divided into 8 equal parts. There are in total 384 elements and 217 nodes.Figure 1
Cross section of the random material sample.Figure 2
Finite element mesh.The relative error is defined as(42)RE=∥CuGNE-CuMC∥2∥CuMC∥2×100%,
in which uMC is the nodal displacement vector solved with the Monte Carlo method and uGNE is the GNE solution. Unless stated explicitly, the numerical testing is performed for 4340 samples generated from the Gaussian field defined in (41).The convergence condition of the GNE method has been rigorously proved. The main purpose of this example is to examine the error distribution and its relationship with various random parameters, for which a number of numerical tests are performed in the following subsections.


## 3.2. Error Estimation

The third-order GNE solution is adopted. The constantsc1 and c2 in (41) are varied to examine the error estimator in (26). The constant c1 is the variance of the Gaussian field, and the constant c2 controls the shape of the covariance function. The effective correlation length Leff of the covariance function R(τ) is defined by
(43)R(Leff)R(0)=0.01.
For any two points in the medium, if the distance between them is larger than Leff, their material properties are only weakly correlated. The effective correlation length of the Gaussian field defined in (41) is
(44)R(τ)R(0)=e-τ2/c2=0.01⟹Leff=2.146c2.
In the first set of tests, the constant c1 is varied to change the variance of the Gaussian field. The parameter settings are
(45)c1=0.01,c2=3.2E-2,Leff=0.38,c1=0.04,c2=3.2E-2,Leff=0.38.
In the second set of tests, the constant c2 is varied to change the shape of the covariance function. The parameter settings are
(46)c1=0.01,c2=0.8E-2,Leff=0.19,c1=0.01,c2=2.4E-2,Leff=0.33,c1=0.01,c2=4.0E-2,Leff=0.43.


### 3.2.1. The Influence ofc1

Corresponding toc1=0.01 and c1=0.04, the standard deviation σ of the stochastic field is 0.1 and 0.2, respectively. A total number of 4340 samples are generated for each case. For each realization, both the Monte Carlo solution and the third-order GNE solution are computed. The relationship between the spectral radius and the relative error is drawn in Figures 3 and 4, in which the red crosses show the actual relative errors corresponding to each sample and the blue line shows the theoretical error estimation. It can been seen that for both cases, the theoretical error estimator in (26) holds accurately and its effectiveness is not affected by the change of the variance controlled by c1.Figure 3
The relation between the relative error and the spectral radius.Figure 4
The relation between the relative error and the spectral radius.


### 3.2.2. The Influence ofc2

The constantc2 is adjusted to take three different values 0.8×10-2, 2.4×10-2, and 4.0×10-2, while the constant c1 remains unchanged. Following (41), the effective correlation length grows with the increase of c2. Again, 4340 samples are generated for each stochastic field, and for each realization both the Monte Carlo solution and the third-order GNE solution are computed. Figure 5 plots the relationship between spectral radius and the relative error, which verifies that, for all cases, the theoretical error estimation in (26) holds. The changing of c2, which changes the variation frequency of the stochastic field, does not affect the error estimation.Figure 5
The relation between the relative error and the spectral radius.As indicated in (44), the effective correlation length Leff is directly associated with the parameter c2. To investigate the influence of Leff on the relative error, the empirical cumulative distribution functions (ECDF) of relative errors are plotted in Figure 6 for different Leff cases. Comparing the three ECDF curves at different effective correlation lengths, it is found that with the decrease of the effective correlation length, the accuracy of the GNE solution improves.Figure 6
ECDF of the relative errors.


## 3.2.1. The Influence ofc1

Corresponding toc1=0.01 and c1=0.04, the standard deviation σ of the stochastic field is 0.1 and 0.2, respectively. A total number of 4340 samples are generated for each case. For each realization, both the Monte Carlo solution and the third-order GNE solution are computed. The relationship between the spectral radius and the relative error is drawn in Figures 3 and 4, in which the red crosses show the actual relative errors corresponding to each sample and the blue line shows the theoretical error estimation. It can been seen that for both cases, the theoretical error estimator in (26) holds accurately and its effectiveness is not affected by the change of the variance controlled by c1.Figure 3
The relation between the relative error and the spectral radius.Figure 4
The relation between the relative error and the spectral radius.


## 3.2.2. The Influence ofc2

The constantc2 is adjusted to take three different values 0.8×10-2, 2.4×10-2, and 4.0×10-2, while the constant c1 remains unchanged. Following (41), the effective correlation length grows with the increase of c2. Again, 4340 samples are generated for each stochastic field, and for each realization both the Monte Carlo solution and the third-order GNE solution are computed. Figure 5 plots the relationship between spectral radius and the relative error, which verifies that, for all cases, the theoretical error estimation in (26) holds. The changing of c2, which changes the variation frequency of the stochastic field, does not affect the error estimation.Figure 5
The relation between the relative error and the spectral radius.As indicated in (44), the effective correlation length Leff is directly associated with the parameter c2. To investigate the influence of Leff on the relative error, the empirical cumulative distribution functions (ECDF) of relative errors are plotted in Figure 6 for different Leff cases. Comparing the three ECDF curves at different effective correlation lengths, it is found that with the decrease of the effective correlation length, the accuracy of the GNE solution improves.Figure 6
ECDF of the relative errors.


## 3.3. The Relative Error for Other Types of Stochastic Fields

In the last two subsections, the error estimation analysis has taken the Gaussian assumption; that is, all random variablesαi in (3) are assumed to be independent Gaussian random variables. In this subsection, we investigate the relative error of the GNE method for other types of stochastic fields. Specifically, the random variables αi in (3) are assumed to have the two-point distribution, the binomial distribution, and the β distribution, respectively. For the case of two-point distribution, values {-0.3E0,0.3E0} are assigned with equal probability, that is, 0.5. For the case of binomial distribution, ten times Bernoulli trials with the probability 0.5 are adopted, and the integer values between 0 and 10 are linearly projected onto [-0.3E0,0.3E0]. For the case of β distribution, the probability density function is defined as
(47)f(x∣a,b)=1B(a,b)xa-1(1-x)b-1I(0,1)(x),
where B(a,b) is a β function, and the filter function I(0,1)(x) is defined as
(48)I(0,1)(x)={1x∈(0,1)0x∉(0,1).
Three kinds of β distributions are considered, with a=b=0.75, a=b=1, and a=b=4, respectively. The interval (0,1) is linearly mapped to [-0.3E0,0.3E0].Corresponding to these five distributions, the empirical cumulative distribution functions of relative errors are plotted in Figure7. Comparing the ECDF curves at different effective correlation lengths, the same trend is observed as the Gaussian case. That is, the relative error decreases with the drop of effective correlation length.ECDF of the relative error.
(a)
Two-point distribution
(b)
Binomial distribution
(c)
β distribution with a=b=0.75
(d)
β distribution with a=b=1
(e)
β distribution with a=b=4


## 3.4. Efficiency Analysis

Using the same computing platform, the Monte Carlo method, the Neumann expansion method, and the GNE method are employed to solve the same set of solutions (corresponding to 4340 samples), and the computation time is shown in Table1. The number of random variables in the stochastic equation system differs depending on the effective correlation length. A fixed mesh is adopted in all of the computation such that the dimension of a single equation system is 434×434.Table 1
Computation time for 4340 samples (unit: s).
Number of random variables
10
11
16
21
39
Monte Carlo
40.78
43.26
54.02
65.55
107.13
GNE
First order
0.44
0.44
0.48
0.53
0.69
Second order
0.79
0.83
1.05
1.28
2.21
Third order
1.30
1.43
2.29
3.58
13.33
Neumann
First order
26.74
28.96
40.43
52.23
93.90
Second order
28.63
30.88
42.28
53.70
95.22
Third order
30.28
32.57
44.05
55.49
97.08If the mesh density is doubled, then the dimension of each matrix is1634×1634. The number of random variables remains the same in each case. The corresponding computation time is shown in Table 2.Table 2
Computation time for 16340 samples (unit: s).
Number of random variables
11
16
21
Monte Carlo
6560.01
7871.31
9118.39
GNEM
First order
13.31
14.85
16.39
Second order
16.06
19.23
22.50
Third order
24.55
41.86
70.97
Neumann
First order
3325.30
4851.50
6017.91
Second order
3639.38
5010.79
6314.89
Third order
3881.43
5473.71
6609.92The two data tables indicate that the computational efficiency of the GNE method (up to the third order) is much higher than Neumann expansion and Monte Carlo procedures. However, the memory requirement for higher order GNE schemes increases rapidly causing the computational efficiency to decrease.


## 4. Conclusions

An acceleration approach of Neumann expansion is proposed for solving stochastic linear equations. The new method, termed generalized Neumann expansion, generalizes Neumann expansion from one matrix to multiple matrices. The GNE method is equivalent to Neumann expansion for general stochastic equations, and for stochastic linear equations, both methods are also equivalent to the perturbation method. According to this equivalence relation, a rigorous error estimation is obtained, which is applicable to all the three methods. Numerical tests show that the third-order GNE expansion can provide good accuracy even for relatively large random fluctuations.



---

*Source: 325025-2013-09-23.xml*