Generalized Neumann Expansion and Its Application in Stochastic Finite Element Methods

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.

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) where 𝐀₁,…,𝐀_𝑚 are 𝑛×𝑛 symmetric and deterministic matrices, 𝐮(𝜔) the unknown stochastic nodal vector, and 𝐟(𝜔) the nodal force vector. If the material properties are modeled as Gaussian fields, then 𝛼₁,𝛼₂,…,𝛼_𝑚 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) where 𝐸(𝐱,𝜔) denotes random Young’s modulus at point 𝐱, 𝐸₀(𝐱) is the expectation of 𝐸(𝐱,𝜔), 𝛼_𝑖(𝜔) are uncorrelated random variables, and 𝜆_𝑖 and 𝜓_𝑖(𝐱) are deterministic eigen values and eigen functions corresponding to the covariance function of the source field 𝐸(𝐱,𝜔). If 𝐸(𝐱,𝜔) is modeled as a Gaussian field, 𝛼_𝑖(𝜔) become independent standard Gaussian random variables. For the sake of simplicity, we drop the random symbol 𝜔 and for the rest of the paper use 𝛼_𝑖 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) In the above equation, 𝐀₀ is the conventional FE stiffness matrix, which is symmetric, positive definite, and sparse while 𝐀_𝑖 are symmetric and sparse matrices. The deterministic matrices 𝐀_𝑖 share the same entry pattern as the stiffness matrix 𝐀₀, but they are usually not positive definite.

2.1. Related Methods

The Monte Carlo methods first generate a set of samples of 𝛼_𝑖 in (3), and each sample of 𝛼_𝑖 corresponds to a deterministic equation, which can be solved by the standard FEM method.

In the perturbation procedure, the unknown random vector 𝐮 is expanded as [21, 25, 26] 𝐮=𝐮₀+_𝑖𝛼_𝑖𝐮_𝑖+_𝑖𝑗𝛼_𝑖𝛼_𝑗𝐮_𝑖𝑗+⋯. (4) Substituting (4) into (3) yields 𝐀₀+_𝑖𝛼_𝑖𝐀_𝑖𝐮₀+_𝑖𝛼_𝑖𝐮_𝑖+_𝑖𝑗𝛼_𝑖𝛼_𝑗𝐮_𝑖𝑗+⋯=𝐟, (5) from which the components 𝐮₀, 𝐮_𝑖, 𝐮_𝑖𝑗 of 𝐮 are solved as follows: 𝐮₀=𝐀₀⁻¹𝐮𝐟,_𝑖=−𝐀₀⁻¹𝐀_𝑖𝐮₀,𝐮_𝑖𝑗=𝐀₀⁻¹𝐀_𝑖𝐀₀⁻¹𝐀_𝑗𝐮₀. (6) 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) where Δ𝐀=∑_𝑖𝛼_𝑖𝐀_𝑖, 𝐁=𝐀₀⁻¹Δ𝐀 and 𝐈 is the identity matrix. According to Neumann expansion (𝐈+𝐁)⁻¹=𝐈−𝐁+𝐁²−𝐁³+⋯, (8) solution (7) can be simplified as 𝐮=𝐮₀−𝐁𝐮₀+𝐁²𝐮₀−𝐁³𝐮₀+⋯, (9) where 𝐮₀=𝐀₀⁻¹𝐟. The spectral radius of matrix 𝐁 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) The generalized Neumann expansion (GNE) in (10) can be proved as follows: 𝐈+_𝑖𝛼_𝑖𝐁_𝑖𝐈+_𝑖𝛼_𝑖𝐁_𝑖⁻¹=𝐈−_𝑖𝛼_𝑖𝐁_𝑖+_𝑖𝑗𝛼_𝑖𝛼_𝑗𝐁_𝑖𝐁_𝑗−⋯+𝛼₁𝐁₁−_𝑖𝛼₁𝛼_𝑖𝐁₁𝐁_𝑖+⋯+𝛼₂𝐁₂−_𝑖𝛼₂𝛼_𝑖𝐁₂𝐁_𝑖⋮+⋯+𝛼_𝑛𝐁_𝑛−_𝑖𝛼_𝑛𝛼_𝑖𝐁_𝑛𝐁_𝑖+⋯=𝐈, (11) where 𝑛 is the total number of 𝐁_𝑖.

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 𝐁_𝑖=𝐀₀⁻¹𝐀_𝑖, and applying the GNE in (10) to (3) gives 𝐮=𝐀₀+_𝑖𝛼_𝑖𝐀_𝑖⁻¹𝐟=𝐈+_𝑖𝛼_𝑖𝐀₀⁻¹𝐀_𝑖⁻¹𝐀₀⁻¹𝐟=𝐮₀−_𝑖𝛼_𝑖𝐀₀⁻¹𝐀_𝑖𝐮₀+_𝑖𝑗𝛼_𝑗𝛼_𝑖𝐀₀⁻¹𝐀_𝑗𝐀₀⁻¹𝐀_𝑖𝐮₀−⋯=𝐮₀−_𝑖𝛼_𝑖𝐮_𝑖+_𝑖𝑗𝛼_𝑖𝛼_𝑗𝐮_𝑖𝑗−⋯ (12) in which 𝐮_𝑖=𝐀₀⁻¹𝐀_𝑖𝐮₀,𝐮_𝑖𝑗=𝐀₀⁻¹𝐀_𝑖𝐀₀⁻¹𝐀_𝑗𝐮₀. (13) 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, 𝜌(∑𝛼_𝑖𝐁_𝑖)<1, where 𝜌(⋅) denotes the spectral radius of a matrix. Hence, for the GNE solution in (12), the convergence condition is 𝜌𝛼_𝑖𝚨₀⁻¹𝐀_𝑖<1. (14) It should be noted that 𝜌𝛼_𝑖𝚨₀⁻¹𝐀_𝑖=𝜌𝛼_𝑖𝐂⁻¹𝐀_𝑖𝐂⁻¹=𝛼_𝑖𝐂⁻¹𝐀_𝑖𝐂⁻¹_𝑚2, (15) where ‖⋅‖_𝑚2 denotes the spectral norm of a matrix, and the symmetric and positive definite matrix 𝐂 is defined as 𝐀₀=𝐂². (16) As a type of matrix norm, the spectral norm has the following properties: ‖𝐏+𝐐‖_𝑚2≤‖𝐏‖_𝑚2+‖𝐐‖_𝑚2,‖𝐏𝐐‖_𝑚2≤‖𝐏‖_𝑚2‖𝐐‖_𝑚2,‖ℎ𝐏‖_𝑚2=ℎ‖𝐏‖_𝑚2, (17) where 𝐏 and 𝐐 are general matrices and ℎ is a positive constant. Hence, 𝛼_𝑖𝐂⁻¹𝐀_𝑖𝐂⁻¹_𝑚2≤𝛼₁𝐂⁻¹𝐀₁𝐂⁻¹_𝑚2+𝛼₂𝐂⁻¹𝐀₂𝐂⁻¹_𝑚2+⋯𝛼_𝑛𝐂⁻¹𝐀_𝑛𝐂⁻¹_𝑚2,𝜌𝛼_𝑖𝐂⁻¹𝐀_𝑖𝐂⁻¹≤𝜌𝛼₁𝐂⁻¹𝐀₁𝐂⁻¹+𝜌𝛼₂𝐂⁻¹𝐀₂𝐂⁻¹+⋯𝜌𝛼_𝑛𝐂⁻¹𝐀_𝑛𝐂⁻¹,𝜌𝛼_𝑖𝚨₀⁻¹𝐀_𝑖≤𝜌𝛼₁𝚨₀⁻¹𝐀₁+𝜌𝛼₂𝚨₀⁻¹𝐀₂+⋯𝜌𝛼_𝑛𝚨₀⁻¹𝐀_𝑛≤𝑛max𝛼_𝑖max𝜌𝚨₀⁻¹𝐀_𝑖. (18) Based on the above relation, a sufficient but not necessary condition for convergence is 𝑛max𝛼_𝑖max𝜌𝚨₀⁻¹𝐀_𝑖<1. (19) The maximum spectral radius 𝜌_max=max(𝜌(𝚨₀⁻¹𝐀_𝑖)) can be obtained from 𝐀₀ and 𝐀_𝑖. Equation (19) can be used to verify the convergence of GNE in advance. When 𝛼_𝑖 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: 𝑛max𝛼_𝑖max𝜌𝚨₀⁻¹𝐀_𝑖≈3𝑛𝜌_max𝜌<1,_max<1.3𝑛 (20) Considering the 𝑘th-order GNE (10) or Neumann expansion (8), their error can be estimated as follows: (𝐈+𝐁)⁻¹−𝐈−𝐁+𝐁²+⋯+(−1)^𝑘𝐁^𝑘_𝑚2=𝐁^𝑘+1(𝐈+𝐁)⁻¹_𝑚2≤𝐁^𝑘+1_𝑚2(𝐈+𝐁)⁻¹_𝑚2=𝜎_𝐵^𝑘+11−𝜎_𝐵, (21) where 𝜎_𝐵 is the spectral norm of 𝐁, which is assumed to satisfy 𝜎_𝐵<1. (22) The exact solution of (3) is 𝐮=𝐈+_𝑖𝛼_𝑖𝐀₀⁻¹𝐀_𝑖⁻¹𝐀₀⁻¹𝐟=(𝐈+𝐁)⁻¹𝐀₀⁻¹𝐟,𝐂𝐮=𝐈+_𝑖𝛼_𝑖𝐂⁻¹𝐀_𝑖𝐂⁻¹⁻¹𝐂⁻¹𝐟=(𝐈+𝐃)⁻¹𝐂⁻¹𝐟, (23) where 𝐁=∑_𝑖𝛼_𝑖𝐀₀⁻¹𝐀_𝑖, 𝐃=∑_𝑖𝛼_𝑖𝐂⁻¹𝐀_𝑖𝐂⁻¹. The 𝑘th-order GNE or Neumann expansion solution is given by 𝐮=𝐈−𝐁+𝐁²+⋯+(−1)^𝑘𝐁^𝑘𝐀₀⁻¹𝐂𝐟,𝐮=𝐈−𝐃+𝐃²+⋯+(−1)^𝑘𝐃^𝑘𝐂⁻¹𝐟. (24) The norm of the solution error is defined as ‖𝐂𝐮−𝐂𝐮‖₂=(𝐈+𝐃)⁻¹−𝐈−𝐃+𝐃²+⋯+(−1)^𝑘𝐃^𝑘𝐂⁻¹𝐟₂=𝐃^𝑘+1(𝐈+𝐃)⁻¹𝐂⁻¹𝐟₂=𝐃^𝑘+1𝐂𝐮₂≤𝐃^𝑘+1_𝑚2‖𝐂𝐮‖₂=𝜎_𝐷^𝑘+1‖𝐂𝐮‖₂=𝜌_𝐵^𝑘+1‖𝐂𝐮‖₂, (25) where 𝜌_𝐵 is the spectral radius of 𝐁.

The relative error of the solution is defined as RE=(𝐮−𝐮)^𝑇𝐀₀(𝐮−𝐮)𝐮^𝑇𝐀₀𝐮=(𝐮−𝐮)^𝑇𝐂²(𝐮−𝐮)𝐮^𝑇𝐂²𝐮=‖𝐂𝐮−𝐂𝐮‖₂‖𝐂𝐮‖₂≤𝜌_𝐵^𝑘+1. (26) Comparing (21) and (26), it can be seen that the accuracy for the solution vector 𝐮 is even better than the accuracy for the matrix inverse (𝐈+𝐁)⁻¹.

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) the three methods are not equivalent. Specifically, applying the perturbation method to the above equation yields 𝐀₀+_𝑖𝛼_𝑖𝐀_𝑖+_𝑖𝑗𝛼_𝑖𝛼_𝑗𝐀_𝑖𝑗×𝐮+⋯₀+_𝑖𝛽_𝑖𝐮_𝑖+_𝑖𝑗𝛽_𝑖𝛽_𝑗𝐮_𝑖𝑗+⋯=𝐟₀+_𝑖𝛾_𝑖𝐟_𝑖+_𝑖𝑗𝛾_𝑖𝛾_𝑗𝐟_𝑖𝑗+⋯. (28) Equation (28) can be reorganized as𝐀₀𝐮₀+∑_𝑖𝛼_𝑖𝐀_𝑖𝐮₀+∑_𝑖∑_𝑗𝛼_𝑖𝛼_𝑗𝐀_𝑖𝑗𝐮₀⋯+∑_𝑖𝐀₀𝛽_𝑖𝐮_𝑖+∑_𝑖∑_𝑗𝛼_𝑖𝐀_𝑖𝛽_𝑗𝐮_𝑗+∑_𝑖∑_𝑗∑_𝑘𝛼_𝑖𝛼_𝑗𝐀_𝑖𝑗𝛽_𝑘𝐮_𝑘⋯+∑_𝑖∑_𝑗𝐀₀𝛽_𝑖𝛽_𝑗𝐮_𝑖𝑗+∑_𝑖∑_𝑗∑_𝑘𝛼_𝑖𝐀_𝑖𝛽_𝑗𝛽_𝑘𝐮_𝑗𝑘+∑_𝑖∑_𝑗∑_𝑘∑_𝑙𝛼_𝑖𝛼_𝑗𝐀_𝑖𝑗𝛽_𝑘𝛽_𝑙𝐮_𝑘𝑙⋯⋯⋯⋯⋱=𝐟₀+_𝑖𝛾_𝑖𝐟_𝑖+_𝑖𝑗𝛾_𝑖𝛾_𝑗𝐟_𝑖𝑗+⋯. (29)If 𝛼_𝑖=𝛽_𝑖=𝛾_𝑖, then the perturbation components 𝐮₀,𝐮_𝑖,𝐮_𝑖𝑗,…,𝐮_{𝑖1𝑖2⋯𝑖𝑛}of 𝐮 can be solved as 𝐀₀𝐮₀=𝐟₀,𝐀₀𝐮_𝑖+𝐀_𝑖𝐮₀=𝐟_𝑖,𝐀₀𝐮_𝑖𝑗+𝐀_𝑖𝐮_𝑗+𝐀_𝑖𝑗𝐮₀=𝐟_𝑖𝑗,⋮𝐀₀𝐮_{𝑖1𝑖2⋯𝑖𝑛}+𝐀_𝑖1𝐮_{𝑖2𝑖3⋯𝑖𝑛}+𝐀_𝑖1𝑖2𝐮_{𝑖3𝑖4⋯𝑖𝑛}+𝐀_{𝑖1𝑖2⋯𝑖𝑛}𝐮₀=𝐟_{𝑖1𝑖2⋯𝑖𝑛}. (30) Solving (27) with the Neumann expansion gives 𝐀⁻¹=𝐀₀+Δ𝐀⁻¹=(𝐈+𝐁)⁻¹𝐀₀⁻¹=𝐈−𝐁+𝐁²−𝐁³+⋯𝐀₀⁻¹, (31) where 𝐁=𝐀₀⁻¹Δ𝐀,Δ𝐀=∑_𝑖𝛼_𝑖𝐀_𝑖+∑_𝑖∑_𝑗𝛼_𝑖𝛼_𝑗𝐀_𝑖𝑗+⋯. Then 𝐮=𝐈−𝐁+𝐁²−𝐁³+⋯𝐀₀⁻¹𝐟₀+Δ𝐟, (32) in which Δ𝐟=_𝑖𝛾_𝑖𝐟_𝑖+_𝑖𝑗𝛾_𝑖𝛾_𝑗𝐟_𝑖𝑗+⋯. (33) Because higher order entries are included in 𝐁, the term 𝐁𝐀₀⁻¹(𝐟₀+Δ𝐟) (the first-order term of the Neumann solution in (32)) does not correspond to the first-order entries of 𝐮. 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 𝐮=𝐀₀+_𝑖𝑖𝛼_𝑖1𝐀_𝑖1+_𝑖2𝑖1𝛼_𝑖1𝛼_𝑖2𝐀_𝑖1𝑖2+⋯+_𝑖𝑛⋯_𝑖2𝑖1𝛼_{𝑖1𝑖2⋯𝑖𝑛}𝐀_{𝑖1𝑖2⋯𝑖𝑛}⁻¹𝐟=𝐈+_𝑖𝑖𝛼_𝑖1𝐁_𝑖1+_𝑖2𝑖1𝛼_𝑖1𝛼_𝑖2𝐁_𝑖1𝑖2+⋯+_𝑖𝑛⋯_𝑖2𝑖1𝛼_{𝑖1𝑖2⋯𝑖𝑛}𝐁_{𝑖1𝑖2⋯𝑖𝑛}⁻¹𝐀₀⁻¹𝐟=𝐈+_𝑖𝛽_𝑖𝐂_𝑖⁻¹𝐀₀⁻¹𝐟=𝐈−_𝑖𝛽_𝑖𝐂_𝑖+_𝑖𝑗𝛽_𝑖𝛽_𝑗𝐂_𝑖𝐂_𝑗−_𝑖𝑗𝑘𝛽_𝑖𝛽_𝑗𝛽_𝑘𝐂_𝑖𝐂_𝑗𝐂_𝑘+⋯𝐀₀⁻¹𝐟, (34) where 𝐟=𝐟₀+∑_𝑖𝛾_𝑖𝐟_𝑖+∑_𝑖∑_𝑗𝛾_𝑖𝛾_𝑗𝐟_𝑖𝑗+⋯, 𝐁_𝑖=𝐀₀⁻¹𝐀_𝑖, 𝐁_𝑖𝑗=𝐀₀⁻¹𝐀_𝑖𝑗, and other 𝐁 terms can be similarly formed. Matrices 𝐂_𝑖 are determined by realignment of 𝐁_𝑖, 𝐁_𝑖𝑗,… and random variables 𝛽_𝑖 are formed by rearrangement of 𝛼_𝑖1,𝛼_𝑖1𝑖2,…,𝛼_{𝑖1𝑖2⋯𝑖𝑛}.

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) The first-order perturbation solution is 𝐮₀=𝐀₀⁻¹𝐟₀,𝐮_𝑖=𝐀₀⁻¹𝐟_𝑖−𝐀₀⁻¹𝐀_𝑖𝐮₀,𝐮=𝐮₀+_𝑖𝛼_𝑖𝐮_𝑖=𝐀₀⁻¹𝐟₀+_𝑖𝛼_𝑖𝐀₀⁻¹𝐟_𝑖−_𝑖𝛼_𝑖𝐀₀⁻¹𝐀_𝑖𝐮₀, (36) and the first-order GNE solution is 𝐮=𝐀₀⁻¹−_𝑖𝛼_𝑖𝐀₀⁻¹𝐀_𝑖𝐀₀⁻¹𝐟₀+_𝑖𝛼_𝑖𝐟_𝑖=𝐀₀⁻¹𝐟₀−_𝑖𝛼_𝑖𝐀₀⁻¹𝐀_𝑖𝐮₀+_𝑖𝛼_𝑖𝐀₀⁻¹𝐟_𝑖−_𝑖𝑗𝛼_𝑖𝛼_𝑗𝐀₀⁻¹𝐀_𝑖𝐟_𝑗. (37) Comparing (36) and (37), it can be seen that the GNE solution contains an extra second-order term −∑_𝑖∑_𝑗𝛼_𝑖𝛼_𝑗𝐀₀⁻¹𝐀_𝑖𝐟_𝑗.

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) Denote 𝐮₀=𝐀₀⁻¹𝐟₀, 𝐁_𝑖=𝐀₀⁻¹𝐀_𝑖, and 𝐁_𝑖𝑗=𝐀₀⁻¹𝐀_𝑖𝑗; then the first-order perturbation solution can be expressed as 𝐮₀=𝐀₀⁻¹𝐟₀,𝐮_𝑖=−_𝑖𝐀₀⁻¹𝐀_𝑖𝐮₀=−_𝑖𝐁_𝑖𝐮₀,𝐮=𝐮₀+_𝑖𝛼_𝑖𝐮_𝑖=𝐮₀−_𝑖𝛼_𝑖𝐁_𝑖𝐮₀, (39) while the first-order GNE solution is 𝐮=𝐈+_𝑖𝛼_𝑖𝐀₀⁻¹𝐀_𝑖+_𝑖𝑗𝛼_𝑖𝛼_𝑗𝐀₀⁻¹𝐀_𝑖𝑗⁻¹𝐮₀=𝐈+_𝑖𝛼_𝑖𝐁_𝑖+_𝑖𝑗𝛼_𝑖𝛼_𝑗𝐁_𝑖𝑗⁻¹𝐮₀=𝐈−_𝑖𝛼_𝑖𝐁_𝑖−_𝑖𝑗𝛼_𝑖𝛼_𝑗𝐁_𝑖𝑗𝐮₀. (40) Comparing (39) and (40), the second-order terms −∑_𝑖∑_𝑗𝛼_𝑖𝛼_𝑗𝐀₀⁻¹𝐀_𝑖𝑗𝐮₀ in the GNE solution do not appear in the perturbation solution (39).

3. Numerical Examples

3.1. Example Description

As shown in Figure 1, 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 𝑅𝜏₁,𝜏₂=𝑐₁𝑒^{−𝜏12/𝑐2−𝜏22/𝑐2}, (41) in which 𝜏₁=𝑥₁−𝑥₂ and 𝜏₂=𝑦₁−𝑦₂ are the coordinate differences and 𝑐₁ and 𝑐₂ 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 RE=𝐂𝐮_GNE−𝐂𝐮_MC2𝐂𝐮_MC2×100%, (42) in which 𝐮_MC is the nodal displacement vector solved with the Monte Carlo method and 𝐮_GNE 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 constants 𝑐₁ and 𝑐₂ in (41) are varied to examine the error estimator in (26). The constant 𝑐₁ is the variance of the Gaussian field, and the constant 𝑐₂ controls the shape of the covariance function. The effective correlation length 𝐿_eff of the covariance function 𝑅(𝜏) is defined by 𝑅𝐿_eff𝑅(0)=0.01. (43) For any two points in the medium, if the distance between them is larger than 𝐿_eff, their material properties are only weakly correlated. The effective correlation length of the Gaussian field defined in (41) is 𝑅(𝜏)𝑅(0)=𝑒^{−𝜏2/𝑐2}=0.01⟹𝐿_eff=2.146√𝑐₂. (44) In the first set of tests, the constant 𝑐₁ is varied to change the variance of the Gaussian field. The parameter settings are 𝑐₁=0.01,𝑐₂=3.2𝐸−2,𝐿_eff𝑐=0.38,₁=0.04,𝑐₂=3.2𝐸−2,𝐿_eff=0.38. (45) In the second set of tests, the constant 𝑐₂ is varied to change the shape of the covariance function. The parameter settings are 𝑐₁=0.01,𝑐₂=0.8𝐸−2,𝐿_eff𝑐=0.19,₁=0.01,𝑐₂=2.4𝐸−2,𝐿_eff𝑐=0.33,₁=0.01,𝑐₂=4.0𝐸−2,𝐿_eff=0.43. (46)

3.2.1. The Influence of 𝑐₁

Corresponding to 𝑐₁=0.01 and 𝑐₁=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 𝑐₁.

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 of 𝑐₂

The constant 𝑐₂ is adjusted to take three different values 0.8×10⁻², 2.4×10⁻², and 4.0×10⁻², while the constant 𝑐₁ remains unchanged. Following (41), the effective correlation length grows with the increase of 𝑐₂. 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 𝑐₂, 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 𝐿_eff is directly associated with the parameter 𝑐₂. To investigate the influence of 𝐿_eff on the relative error, the empirical cumulative distribution functions (ECDF) of relative errors are plotted in Figure 6 for different 𝐿_eff 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 𝛼_𝑖 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 𝛼_𝑖 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.3𝐸₀,0.3𝐸₀} 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.3𝐸₀,0.3𝐸₀]. For the case of 𝛽 distribution, the probability density function is defined as 𝑓()=1𝑥∣𝑎,𝑏𝐵()𝑥𝑎,𝑏^𝑎−1(1−𝑥)^𝑏−1𝐼_(0,1)(𝑥), (47) where 𝐵(𝑎,𝑏) is a 𝛽 function, and the filter function 𝐼_(0,1)(𝑥) is defined as 𝐼_(0,1)(𝑥)=()().1𝑥∈0,10𝑥∉0,1 (48) Three kinds of 𝛽 distributions are considered, with 𝑎=𝑏=0.75, 𝑎=𝑏=1, and 𝑎=𝑏=4, respectively. The interval (0,1) is linearly mapped to [−0.3𝐸₀,0.3𝐸₀].

Corresponding to these five distributions, the empirical cumulative distribution functions of relative errors are plotted in Figure 7. 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.

(a) Two-point distribution

(b) Binomial distribution

(c) 𝛽 distribution with 𝑎=𝑏=0.75

(d) 𝛽 distribution with 𝑎=𝑏=1

(e) 𝛽 distribution with 𝑎=𝑏=4

Figure 7 

ECDF of the relative error.

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 Table 1. 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).

If the mesh density is doubled, then the dimension of each matrix is 1634×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).

The 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.

Conflict of Interests

The authors do not have any conflict of interests with the content of the paper.

Acknowledgments

The authors would like to acknowledge the financial supports of the National Natural Science Foundation of China (11272181), the Specialized Research Fund for the Doctoral Program of Higher Education of China (20120002110080), and the National Basic Research Program of China (Project no. 2010CB832701). The authors would also like to thank the support from British Council and Chinese Scholarship Council through the support of an award for Sino-UK Higher Education Research Partnership for PhD Studies.