Multiresolution Analysis for Stochastic Finite Element Problems with Wavelet-Based Karhunen-Loève Expansion

Multiresolution analysis for problems involving random parameter fields is considered. The random field is discretized by a Karhunen-Loève expansion. The eigenfunctions involved in this representation are computed by a wavelet expansion. The wavelet expansion allows to control the spatial resolution of the problem. Fine and coarse scales are defined, and the fine scales are taken into account by projection operators. The influence of the truncation level for the wavelet expansion on the computed reliability is documented.


Introduction
Reliability analysis for problems involving random parameter fields is concerned with the solution of stochastic elliptic boundary value problems.For the solution of stochastic boundary value problems, the random parameter field has to be discretized.This can be accomplished either by a Karhunen-Loève expansion or a projection on a polynomial basis 1, 2 .As a result of this discretization, the random parameter field can be approximated by an expression containing only a finite number of random variables.The reliability problem reduces then to the computation of failure probabilities with respect to a finite number of random variables.
After introduction of the discretization, the stochastic boundary value problem becomes equivalent to a deterministic one, for which the approximation on the physical domain, and on the stochastic domain can be treated differently.The deterministic boundary value problem can be solved by standard approximation methods.On the physical domain, finite element FE approximations are prevalent.Galerkin projection methods or collocation schemes on the stochastic domain give rise to stochastic Galerkin methods 1 and stochastic collocation schemes 3-5 .Their theoretical foundation has been laid in 6-9 , where local and global polynomial chaos expansions were investigated and where a priori error estimates have been proved for a fixed number of terms of the Karhunen-Loève expansion.Approximate solutions of the stochastic boundary value problem can be viewed as local stochastic response surfaces 10 that depend on three parameters: a discretization parameter for the physical domain, a discretization parameter for the stochastic domain, and the discretization level of the input random field.
Although the focus of many investigations concerning the approximate solution of stochastic boundary value problems lies on the computation of the first-and second-order moments of the solution, already, in 1 , reliability computation techniques were described, that are based on series representations of the response distribution, on the reliability index, or on Monte Carlo simulation methods.In most of the papers that followed, for example, 11, 12 and references therein, either the reliability index, or simulation methods has been employed.Alternatively, once the random field is discretized, sampling-based response surfaces 13, 14 can be applied as well.Except for the stochastic Galerkin method, all methods provide nonintrusive algorithms, which allow combining the solution procedure with repetitive runs of an FE solver for deterministic problems.
In this paper, the focus is on the discretization of the random field itself.The discretization of the random field may enforce a rather fine FE mesh and thus increase the computational effort.The Karhunen-Loève expansion of the random field requires the solution of an eigenvalue problem for the determination of the expansion functions.Recently, a wavelet expansion of the eigenfunctions has been introduced 15 .The advantage of wavelet bases resides in the fact that they are of localized compact support, which lead to sparse representations of functions and integral operators.Moreover, the discrete wavelet transform provides an efficient means to solve the integral equation related to the determination of the eigenfunctions.Wavelet expansions have also been studied in the context of simulation of random fields 16, 17 and for the solution of stochastic dynamic systems by polynomial chaos expansion 18 .
Here, the wavelet representation of the eigenfunctions is introduced into the stochastic FE approximation procedure in order to control the spatial resolution of the eigenfunctions.This leads to a multiresolution approximation scheme, where finer scales are taken into account by projection operators.The truncation level of the wavelet expansion has been identified as an additional parameter of the metamodel.As the multiresolution analysis is based on the eigenfunctions involved in the Karhunen-Loève expansion of the random field, it can be combined with stochastic Galerkin, stochastic collocation, or sampling-based methods.In this paper, only the latter are considered in the examples.
The paper is organized as follows: the next section gives an introduction to multiresolution analysis.Section 3 discusses random field discretization by Karhunen-Loève expansion with emphasis on wavelet approximations of the eigenfunctions.Section 4 provides basic information on the approximate solution of the stochastic boundary value problem.Section 5 explains the proposed multiresolution approximation, and Section 6 gives some details on reliability assessment.In Section 7, the proposed procedure is validated by examples and relative errors for the failure probability are given.Finally, conclusions are drawn in Section 8.

Multiresolution Analysis
A multiresolution analysis is a sequence of subspaces For every j ∈ Z, denote with W j the orthogonal complement of V j in V j−1 : where It can be shown 19 that there exists a wavelet ψ x , such that ψ j,k x 2 −j/2 ψ 2 −j x − k is an orthonormal basis for W j .Thus, any function f ∈ L 2 R can be decomposed into where P m is the projection operator on V m .As we are interested in functions with compact support, the second sum in 2.5 is finite.For m 0, we have after truncation at level −ν where the error e x is an element of W − ν 1 .Note that all sums are finite, if f x has compact support.
The same representation can be obtained for functions L 2 R d with d ∈ N.For example, for d 2, the functions Mathematical Problems in Engineering constitute an orthonormal basis of V j and the functions

2.8
where superscript h stands for horizontal, υ for vertical, and d for diagonal translation of the unidimensional wavelet ψ x , constitute an orthonormal basis of W j .

Karhunen-Lo ève Discretization of Random Fields
Let D be a convex bounded open set in R n and Ω, F, P a complete probability space, where Ω is the set of outcomes, F the σ-field of events, and P : F → 0 : 1 a probability measure.We consider a random field α : D × Ω → R that has a continuous and squareintegrable covariance function where the expectation operator E α x Ω α x, ω dP ω denotes the mean value of the random field.It is assumed that α x, ω is bounded and coercive, that is, there exist positive constants a min , a max , such that Due to the properties of the covariance function, the operator T : is compact and self-adjoint and thus admits a spectrum of decreasing nonnegative eigenvalues x, y w i x dx λ i w i y .

3.4
The corresponding eigenfunctions {w i x } ∞ i 1 are orthonormal in L 2 D .The random variables given by Mathematical Problems in Engineering 5 are uncorrelated but in general not independent , have zero mean and unit variance, and allow to represent the random field by the Karhunen-Loève expansion Conditions for stronger convergence properties are given in 8 .The Karhunen-Loève expansion is usually truncated by retaining only the first M terms.In order to keep the computational effort small, a fast decay of the spectrum of 3.3 is important.It is shown in 21 that fast eigenvalue decay corresponds to smoothness of the covariance function.
The solution of the Fredholm integral equation 3.4 can be computed by means of the wavelet basis introduced before cf., e.g., 15 .For the kth eigenfunction w k x , an approximation with a finite number of basis functions Ψ i x , i 1, 2, . . ., n, of V 0 and W j , 0 ≤ j ≤ −ν, is given by 3.7 A Galerkin technique can be applied in order to compute the eigenvalue λ k and the coefficients of the normalized D w 2 k x dx 1 approximate eigenfunctions.To this end, the representation 3.7 is inserted into the Fredholm integral equation 3.4 and the equation is multiplied with a basis function Ψ j y .After integration over the domain D, one obtains the algebraic eigenvalue problem where the coefficients on the left hand side can be obtained from discrete wavelet transforms of the covariance function C x, y .More details and convergence studies for this approximation of the eigenfunctions may be found in 15 .Here, the decomposition is introduced in order to obtain a coarse approximation by taking scales until a certain level −μ > −ν into account and a fine approximation scales from −μ − 1 to −ν .

Stochastic Linear Elliptic Boundary Value Problems
Consider the following model problem with stochastic operator and deterministic input function on D × Ω: find u : D × Ω → R, such that P almost surely: It is assumed that the deterministic input function g x is square-integrable.

Mathematical Problems in Engineering
We are interested in the probability that a functional F u of the solution u x, ω exceeds a threshold F 0 , that is, we want to evaluate the integral 1/2 , the space L 2 P Ω of square integrable random variables and the tensor product space H 1 0 D ⊗ L 2 P Ω of H 1 0 Dvalued random fields with finite second-order moments, equipped with the inner product The variational formulation of the stochastic linear elliptic boundary value problem 4.1 then The assumptions on the random field α x, ω guarantee the continuity and coercivity of the bilinear form in 4.4 .Existence and uniqueness of a solution to 4.4 follow from the Lax-Milgram lemma.For a prescribed, uniformly bounded random field α x, ω , the random variables ξ i ω in 3.6 would be dependent non-Gaussian random variables whose joint distribution function is very difficult to identify.If, on the other hand, independent but bounded distributions are prescribed for ξ i ω , i 1, 2, . . ., M, the random field α x, ω is not necessarily bounded for M → ∞.Thus, one is left with Gaussian distributions for ξ i ω and α x, ω , with transformations of Gaussian random fields 22 or with some situations, where nonnegative distributions for ξ i ω lead to meaningful e.g., Erlang-distributions for α x, ω .Some authors 11, 23 consider models that contain a finite number of random variables as parameters of a boundary value problem.For this kind of problems, Babuška et al. 8 investigate convergence properties of several approximation schemes.The examples studied in this paper belong to this class of problems.It is henceforth assumed that the random field α x, ω is given by a truncated sum where ξ i ω are continuous and independent random variables with zero mean and unit variance, where Γ i ξ i Ω are bounded intervals in R and the parameters λ i and functions w i x are the first M eigenvalues and corresponding eigenfunctions of the operator 3.3 .In the examples, ξ i ω are Gaussian random variables truncated at ±3σ, σ 1.The representation of α M x, ω by a finite number of random variables allows to consider a deterministic auxiliary problem: where From the expression for u x, y , one obtains u x, ω by replacing the vector y with the vector of random variables ξ i ω , i 1, 2, . . ., M.
This equation is discretized on finite dimensional approximation spaces for the physical domain.For H 1 0 D , a family of standard FE approximation spaces X h ⊂ H 1 0 D of continuous piecewise linear functions in a regular triangulation T h of D with mesh parameter h is considered.
Denote with N i x , i 1, 2, . . ., N, a basis of X h ⊂ H 1 0 D .The solution u x, y is approximated by For a fixed value y, the unknown coefficients u i y , i 1, 2, . . ., N, can be computed from the solution of the FE problem where the matrices K 0 , K s , s 1, 2, . . ., M and the vector g are given by and u y is the vector containing the nodal displacements u i y , i 1, 2, . . ., N, for the fixed value y.The matrices K s , s 1, 2, . . ., M, can be interpreted, for example, as FE stiffness matrices for a spatial variation of elastic properties.

Multiresolution Approximation
By means of multiresolution analysis, each of the matrices K s in 4.10 can be decomposed as K s Given a fixed value for y, the decomposition of K s leads to a decomposition of the global stiffness matrix: where . Also the solution u is decomposed into u c and u f , where u c is represented on a coarse mesh and u f on a fine mesh.With the help of an interpolation matrix P, the solution u is written as u Pu c u f .Inserting these decompositions into equation 5.1 leads to This equation is split into 3 Multiplication of the first equation with P T projects this equation onto the coarse mesh:

5.5
For u f , a coarse mesh approximation is computed in the following manner: ii u f is represented on the coarse mesh by u f , thus u f P u f .
After multiplication of 5.4 with P T , one then obtains for the coarse mesh approximation u f of u f .Thus, the linear system of equations in 5.5 and 5.6 has to be resolved only on the reduced set of degrees of freedom given by the coarse mesh, which saves a considerable amount of CPU time.The interpolation matrix P can be generated from FE interpolation functions.

Reliability Assessment
Once the algebraic problems are solved and the correction to the coarse scale solution is computed, an approximation for u x, ω has been obtained.The approximation quality depends on the following parameters: i partition of the physical domain and the stochastic domain, ii truncation of the Karhunen-Loève expansion M , iii wavelet scales μ, ν .
The wavelet scales μ, ν appear as an additional parameter that influences the approximation quality.
For solving reliability problems, 4.8 yields a functional relationship between the input random variables and u x, y .It is then possible to compute the most probable point of failure MPP , that is, the point ξ ∈ Γ with F u x, ξ F 0 with lowest Euclidean norm.The norm of the MPP may serve as a control variable for the approximation quality and for the adaptation of the parameters mentioned above.
The MPP may also be useful for the evaluation of the integral in 4.2 by means of variance reduced Monte Carlo simulation importance sampling .To this end, a sampling density p y is introduced by shifting the original probability density function p y of the random variables ξ i ω , i 1, . . ., M, to the previously obtained MPP, and 4.2 is approximated by p y j p y j p y j , where the sampling points y j , j 1, 2, . . ., N s , are generated according to p y .

Example 1: Clamped Square Plate
The first example deals with a standard problem for stochastic FE techniques, a clamped thin square plate under uniform in-plane tension q cf. 1 .The problem is depicted in Figure 1.
The plate has unit length l.The product of Young's modulus and the thickness of the plate are assumed to be an isotropic normal random field with covariance function with standard deviation σ 0.2 and unit mean value.Poisson's ratio is set to 0.3.A first investigation is concerned with the efficiency of the multiresolution analysis.To this end, relative errors for the maximum longitudinal displacement have been computed by comparing the coarse mesh-coarse level μ −1 solution and the coarse solution corrected by the fine level solution ν −4 to the fine mesh-fine level reference solution i.e., μ 0, ν −4 .Dividing the relative error of the coarse scale solution by the relative Mathematical Problems in Engineering  error of the corrected solution yields an improvement factor that indicates the efficiency of the multiresolution correction.For a coarse mesh of 2 × 2 quadrilateral elements only and a fine mesh of 4 × 4 quadrilateral elements, Figures 2 and 3 display the improvement factors for a varying number of Karhunen-Loève expansion terms and for a variation of the correlation length l c l x l y , respectively.Here and in the following, the fine mesh is generated from the coarse mesh by halving the edge length of the elements.For Figure 2, the correlation length was l c 0.5, and, for Figure 3, 12 Karhunen-Loève expansion terms have been retained.The figures reveal that the multiresolution scheme is more efficient if fluctuations become more important, that is for a large number of Karhunen-Loève expansion terms and low correlation lengths.
Next, the influence of the multiresolution correction on the prediction of failure has been investigated.A threshold value of 1.5 is assumed for the maximum longitudinal displacement, and failure occurs if the maximum longitudinal displacement exceeds this threshold value.
The random field has been discretized by a Karhunen-Loève expansion with M 4 terms for a correlation length l c 1.The eigenfunctions were computed either on a coarse 8 × 8 elements or a fine 16 × 16 elements mesh.A reference solution MPP has been computed on the fine mesh including up to seven scales for the wavelet expansion of the eigenfunctions, and a reference result for the failure probability has been obtained on the fine mesh from importance sampling with 30000 samples.These results allow to evaluate the  error due to the multiresolution scheme.To this end, relative errors for the Euclidean norm of the MPP and the failure probability with respect to the reference results have been computed from the multiresolution scheme with one coarse level μ −1 and up to seven fine levels.Figure 4 presents the development of the relative errors for the norm of the design point.Figure 5 displays the relative error of the failure probability.From the figure, an exponential decrease of the error with respect to the number of fine scales involved can be deduced.The coarse mesh solution needed less than half of the CPU time of the fine scale solution.

Example 2: Soil-Structure Interaction
The second example, previously introduced in 24 , considers a soil structure interaction problem.In contrast to the first example, it deals with an anisotropic autocorrelation function of the random field and an anisotropic FE mesh.The settlement of a foundation width 2B 10 m , represented by a uniform pressure of 0.2 MPa, on an elastic soil layer of thickness Relative error (%)   30 m lying on a rigid substratum, is investigated.A plane strain deformation of the soil is assumed with linear elastic material properties cf. Figure 6 .The Young's modulus is assumed to be a random field with mean value 50 MPa, 20% coefficient of variation, and an exponential correlation function as in 7.1 but with correlation length l x 250 m, l y 100 m.Poisson's ratio is set to 0.3.Six truncated Gaussian random variables M 6 for the representation of the random field have been considered.

Mathematical Problems in Engineering
Starting from the center of the foundation, the right half of the soil layer is meshed by quadrilateral finite elements until a length of L 60 m, as indicated in Figure 6.The coarse mesh consists of 24 elements, while the fine mesh is divided into 96 elements.The reliability problem deals with the displacement of point A in Figure 6, situated at the center of the foundation and on the surface of the soil, where the vertical displacements of the corresponding deterministic problem i.e., considering only the mean value of Young's modulus attain their maximum.Failure is defined as exceedance of a limit value of 0.1 m for the vertical displacement at point A. Figures 7 and 8 summarize the behavior of the relative errors for the Euclidean norm of the MPP and for the failure probability, respectively.As in the previous example, μ has been set to −1 and ν has been varied.It can be seen that, while the error in the norm is relatively small, rather large errors for the probability of failure may occur.As in the previous example, the relative error for the failure probability decreases  exponentially with the wavelet level −ν; however, this decrease starts only after taking a certain number of fine scales 4 and more into account.

Conclusions
In this paper, wavelet-based Karhunen-Loève expansion has been extended to multiresolution analysis of stochastic FE problems.The proposed procedure allows to include fine scale corrections of the FE solution by means of coarse mesh computations.Improvement due to these corrections has been highlighted by examples, and it could be seen that the relative error of the failure probability decreases exponentially with the wavelet level.
Although the procedure has been combined in this paper to importance samplingbased computation of failure probabilities, it applies to stochastic Galerkin and stochastic collocation methods as well, because it is solely based on the decomposition of the Karhunen-Loève eigenfunctions.However, it relies on linear problems.
Typical applications for multiresolution analysis in the context of stochastic FE methods arise in multiscale analysis of material properties, especially if there is scale coupling and length scales that correspond to the length scales of the representative volume element have to be taken into account.For this kind of application, the limitation of wavelets to rather simple geometries, as those in the examples, will not be a drawback, since simple geometries are frequently used on the micro-and mesoscale and more complex geometries could be embedded into simple ones for the solution of the Fredholm integral equation.

4 . 2 where
χ I • , the indicator function, assumes the value 1 in the interval I and vanishes elsewhere.The variational formulation of the stochastic boundary value problem necessitates the introduction of the Sobolev space H 1 0 D of functions having generalized derivatives in L 2 D and vanishing on the boundary ∂D with norm u H 1 0 D D |∇u| 2 dx contains the low-frequency content, expressed by the Mathematical Problems in Engineering multiresolution scheme until the level −μ, while K s f takes the high frequency content of the levels −μ − 1 until −ν into account.

Figure 1 :
Figure 1: Thin square plate under uniform in-plane tension.

Figure 2 :
Figure 2: Efficiency of the multiresolution scheme with respect to the number of Karhunen-Loève expansion terms.

Figure 3 :
Figure 3: Efficiency of the multiresolution scheme with respect to the correlation length.

Figure 4 :
Figure 4: Relative error of the norm for the MPP.

Figure 5 :
Figure 5: Relative error of the failure probability.

Figure 7 :
Figure 7: Relative error of the norm for the MPP.

Figure 8 :
Figure 8: Relative error of the failure probability.