# Hybrid Stochastic Finite Element Method for Mechanical Vibration Problems

Authors: Harri Hakula; Mikael Laaksonen
Journal: Shock and Vibration
 (2015)

---

## Abstract

We present and analyze a new hybrid stochastic finite element method for solving
eigenmodes of structures with random geometry and random elastic modulus.
The fundamental assumption is that the smallest eigenpair is well defined
over the whole stochastic parameter space.
The geometric uncertainty is resolved using collocation and random material
models using Galerkin method at each collocation point. The response statistics,
expectation and variance of the smallest eigenmode, are computed in numerical
experiments. The hybrid approach is superior to alternatives in practical
cases where the number of random parameters used to describe geometric uncertainty
is much smaller than that of the material models.


---

## Body

## 1. Introduction

In standard engineering models many physical quantities such as material parameters are taken to be constant, even though their statistical nature is well understood. Similarly, assumptions of geometric constants, such as thickness of a structure, are not realistic due to manufacturing imperfections. In a detailed report on a state-of-the-art verification and validation process comparing modern simulations with the set of experiments performed in the Oak Ridge National Laboratory in the early 70s, Szabo and Muntges report discrepancies of over 20% in some quantities of interest [1]. These discrepancies are attributed to machining imperfections not accounted for in the computations. Also, in important nonlinear problems such as buckling of a shell, it is known that variation between manufactured specimens has a profound effect in the actual performance [2]. This suggests that a stochastic dimension should be added to the models.The modern era of uncertainty quantification starts with the works of Babuska et al. [3, 4] and the ETH-group led by Schwab and Gittelson (e.g., [5]) with provably faster convergence rates than the standard Monte Carlo methods. Although the so-called stochastic finite elements had been studied for a relative long time before (the classic reference is Ghanem and Spanos [6]), their application was thus limited to highly specific cases.The solution methods can broadly speaking be divided intointrusive andnonintrusive ones. The same division applies to stochastic eigenvalue problems (SEVPs) as well. SEVPs have attracted a lot of attention recently and various algorithms have been suggested for computing approximate eigenpairs [7–10], specially the power iteration [11] and inverse iteration [12].In this paper our focus is on effects of material models and manufacturing imperfections of geometric nature. It should be noted that in the context of this paper it is assumed that the problems are positive definite and the eigenpair of interest is the ground state, that is, the one with the smallest eigenvalue which, in theory, can be a double eigenvalue. Our experimental setup is nonsymmetric and, thus, a spectral gap exists and the inverse iteration converges to the desired eigenmode.In stochastic eigenvalue problems one must address two central issues that do not arise in stochastic source problems: first, the eigenmodes are defined only up to a sign and, second, the eigenmodes must be normalized over the whole parameter space; that is, every realization must be normalized in the same way. Here our quantity of interest is the eigenvalue and, therefore, we do not necessarily have to fix the signs. The normalization is handled by solution of a nonlinear system of equations as in [12].The main result of this paper is the new hybrid algorithm which combines a nonintrusive outer loop (collocation) with an intrusive inner one (Galerkin). The randomness in geometry is resolved using collocation and that of materials with Galerkinat every collocation point. The model problem is an idealized engine bed vibration problem, where a 2D hollow square is clamped from one side (see Figure 1). The inner cavity is assumed to be of fixed shape but random orientation and Young’s modulus is taken to be random in the normal direction into the body. The choice of the geometric uncertainty models the case of sand casting with moulds of fixed shape. The hybrid approach combines the geometric flexibility of collocation with much faster computation times of the Galerkin approach used in the inner loop and, thus, combines the best of both worlds.Figure 1
Two realisations of the domain with meshes. The nominal domain is the reference. The rotation and scaling of the cavity as well as the direction of the uncertainty in the Young’s modulus is indicated.
(a)
Nominal domain
(b)
Perturbed domainThe rest of the paper is organized as follows: first the abstract problem is introduced in Section2. The representation of the random input is given in Section 3. Higher order finite elements are described in Section 4. Section 5 is central to our discussion; both solution methods, collocation and spectral inverse iteration, are presented. Numerical experiments are discussed in Section 6 and finally conclusions are drawn in Section 7.


## 2. Problem Statement

We start by presenting the deterministic system of Navier’s equations of elasticity and the corresponding weak form. We then extend this to the stochastic setting in order to cover the case of uncertain domain and modulus of elasticity. We aim to compute statistics of the smallest eigenvalue of the resulting stochastic system. The geometry of the computational domain is illustrated in Figure1.


### 2.1. Deterministic Formulation

We use the Navier equations of elasticity to model the system of interest. Find the eigenvalueκ, the displacement field u=(u1,u2), and the symmetric stress tensor σ=σiji,j=12, such that(1)σ=λdiv⁡uI+2μϵu,inD,-div⁡σ=κu,inD,u=0,on∂DD,σ·n=0,on∂DN,where ∂D=∂DD∪∂DN is a partitioned boundary of D. The Lamé constants are(2)λ=Eν1+ν1-2ν,μ=E21+ν,with E and ν being Young’s modulus and Poisson’s ratio, respectively. Further, I is the identity tensor, n denotes the outward unit normal to ∂DN, and the strain tensor is(3)ϵu=12∇u+∇uT.The vector-valued tensor divergence is(4)div⁡σ=∑j=12∂σij∂xji=12.This formulation assumes a constitutive relation corresponding to linear isotropic elasticity with stresses and strains related by Hooke’s generalized law:(5)σv=σ11σ22σ12=Dλ,μϵ11ϵ22ϵ12=Dλ,μϵv.


### 2.2. Weak Formulation

Let us introduce a function spaceVD={v:v∈H1(D),v∂DD=0} and define the eigenproblem: find the eigenpair (κ,u)∈R×(VD×VD) such that(6)au,v=κu,vD,∀v∈VD×VD,where the bilinear form(7)au,v=∫Dσu:ϵvdxis the integrated tensor contraction(8)σ:ϵ=∑i,j=12σijϵij,and ·,·D denotes the mass matrix(9)u,vD=∫Du1v1+u2v2dx.As usual, we define the kinematic relation(10)
ϵ
v
u
=
∂
∂
x
1
0
0
∂
∂
x
2
∂
∂
x
2
∂
∂
x
1
u
1
u
2and specify the constitutive matrix(11)Dλ,μ=λ+2μλ0λλ+2μ000μfor the purpose of rewriting the bilinear form as(12)au,v=∫DϵvuTDλ,μϵvvdx.


### 2.3. Stochastic Extension of the Eigenproblem

We introduce two probability spaces(Ω1,F1,P1) and (Ω2,F2,P2) to model randomness of the domain and the elastic modulus, respectively. Here Ωi denotes the set of outcomes, Fi is a σ-algebra of events, and Pi is a probability measure on Ωi for each i=1,2. Furthermore we set(13)Ω,F,P=⨂i=12⁡Ωi,Fi,Pito be the underlying product probability space.The geometry of the computational domainD is illustrated in Figure 1. We assume that the radius of the cavity R:Ω1→R and the angle θ:Ω1→R are random variables satisfying Rmin⁡≤R(ω1)≤Rmax⁡ and θmin⁡≤θ(ω1)≤θmax⁡ for some suitable constants Rmin⁡,Rmax⁡,θmin⁡,θmax⁡>0. Thus we may write D=D(ω1)=D(R(ω1),θ(ω1)). In the numerical examples we set R(ω1) and θ(ω1) to be uniformly distributed random variables.Suppose that for eachω1∈Ω1 the Young modulus is a random field on the domain D(ω1). By this we mean that E=E(ω1,·,·):Ω2×D(ω1)→R for every ω1∈Ω1. We assume that E(ω1,ω2,x) is strictly positive and bounded; that is, for some constants Emin⁡,Emax⁡>0 it holds that(14)Emin⁡≤Eω1,ω2,·≤Emax⁡inDω1.The stochastic extension of eigenproblem (6) is now given by the following: find the eigenvalue κ:Ω→R and the eigenfunction u:Ω→VD(ω1)×VD(ω1) such that(15)aω;uω1,ω2,·,v=κω1,ω2uω1,ω2,·,vDω1,∀v∈VDω1×VDω1holds for P-almost every ω=(ω1,ω2)∈Ω. Here a(ω;u,v) is the stochastic equivalent of the bilinear form a(u,v) and is given by(16)aω;u,v=∫Dω1σω,u:ϵvdx=∫Dω1ϵvuT·Dλω1,ω2,·,μω1,ω2,·ϵvvdx.The Lamé constants λ(ω1,ω2,x) and μ(ω1,ω2,x) are simply random fields defined by (2) with E=E(ω1,ω2,x).In most cases we are mainly interested in computing statistical moments of the solution. Suppose that we have a random variablev(ω1,ω2). The expected value or mean of v(ω1,ω2) is defined to be(17)Ev=∫Ωvω1,ω2dPωand its variance is(18)Var⁡v=Ev-Ev2.In this paper our aim is to compute the expected value and variance of the smallest eigenvalue of system (15).


## 2.1. Deterministic Formulation

We use the Navier equations of elasticity to model the system of interest. Find the eigenvalueκ, the displacement field u=(u1,u2), and the symmetric stress tensor σ=σiji,j=12, such that(1)σ=λdiv⁡uI+2μϵu,inD,-div⁡σ=κu,inD,u=0,on∂DD,σ·n=0,on∂DN,where ∂D=∂DD∪∂DN is a partitioned boundary of D. The Lamé constants are(2)λ=Eν1+ν1-2ν,μ=E21+ν,with E and ν being Young’s modulus and Poisson’s ratio, respectively. Further, I is the identity tensor, n denotes the outward unit normal to ∂DN, and the strain tensor is(3)ϵu=12∇u+∇uT.The vector-valued tensor divergence is(4)div⁡σ=∑j=12∂σij∂xji=12.This formulation assumes a constitutive relation corresponding to linear isotropic elasticity with stresses and strains related by Hooke’s generalized law:(5)σv=σ11σ22σ12=Dλ,μϵ11ϵ22ϵ12=Dλ,μϵv.


## 2.2. Weak Formulation

Let us introduce a function spaceVD={v:v∈H1(D),v∂DD=0} and define the eigenproblem: find the eigenpair (κ,u)∈R×(VD×VD) such that(6)au,v=κu,vD,∀v∈VD×VD,where the bilinear form(7)au,v=∫Dσu:ϵvdxis the integrated tensor contraction(8)σ:ϵ=∑i,j=12σijϵij,and ·,·D denotes the mass matrix(9)u,vD=∫Du1v1+u2v2dx.As usual, we define the kinematic relation(10)
ϵ
v
u
=
∂
∂
x
1
0
0
∂
∂
x
2
∂
∂
x
2
∂
∂
x
1
u
1
u
2and specify the constitutive matrix(11)Dλ,μ=λ+2μλ0λλ+2μ000μfor the purpose of rewriting the bilinear form as(12)au,v=∫DϵvuTDλ,μϵvvdx.


## 2.3. Stochastic Extension of the Eigenproblem

We introduce two probability spaces(Ω1,F1,P1) and (Ω2,F2,P2) to model randomness of the domain and the elastic modulus, respectively. Here Ωi denotes the set of outcomes, Fi is a σ-algebra of events, and Pi is a probability measure on Ωi for each i=1,2. Furthermore we set(13)Ω,F,P=⨂i=12⁡Ωi,Fi,Pito be the underlying product probability space.The geometry of the computational domainD is illustrated in Figure 1. We assume that the radius of the cavity R:Ω1→R and the angle θ:Ω1→R are random variables satisfying Rmin⁡≤R(ω1)≤Rmax⁡ and θmin⁡≤θ(ω1)≤θmax⁡ for some suitable constants Rmin⁡,Rmax⁡,θmin⁡,θmax⁡>0. Thus we may write D=D(ω1)=D(R(ω1),θ(ω1)). In the numerical examples we set R(ω1) and θ(ω1) to be uniformly distributed random variables.Suppose that for eachω1∈Ω1 the Young modulus is a random field on the domain D(ω1). By this we mean that E=E(ω1,·,·):Ω2×D(ω1)→R for every ω1∈Ω1. We assume that E(ω1,ω2,x) is strictly positive and bounded; that is, for some constants Emin⁡,Emax⁡>0 it holds that(14)Emin⁡≤Eω1,ω2,·≤Emax⁡inDω1.The stochastic extension of eigenproblem (6) is now given by the following: find the eigenvalue κ:Ω→R and the eigenfunction u:Ω→VD(ω1)×VD(ω1) such that(15)aω;uω1,ω2,·,v=κω1,ω2uω1,ω2,·,vDω1,∀v∈VDω1×VDω1holds for P-almost every ω=(ω1,ω2)∈Ω. Here a(ω;u,v) is the stochastic equivalent of the bilinear form a(u,v) and is given by(16)aω;u,v=∫Dω1σω,u:ϵvdx=∫Dω1ϵvuT·Dλω1,ω2,·,μω1,ω2,·ϵvvdx.The Lamé constants λ(ω1,ω2,x) and μ(ω1,ω2,x) are simply random fields defined by (2) with E=E(ω1,ω2,x).In most cases we are mainly interested in computing statistical moments of the solution. Suppose that we have a random variablev(ω1,ω2). The expected value or mean of v(ω1,ω2) is defined to be(17)Ev=∫Ωvω1,ω2dPωand its variance is(18)Var⁡v=Ev-Ev2.In this paper our aim is to compute the expected value and variance of the smallest eigenvalue of system (15).


## 3. Spectral Representations

For purposes of numerical treatment we need to approximate the random Young modulus using only a finite number of random variables. A natural way of achieving this is to use the truncated Karhunen-Loève expansion. Once a finite-dimensional approximation has been obtained we may express our solution in a generalized polynomial chaos basis and apply solution schemes based on the stochastic Galerkin finite element method; see [6] for a thorough introduction to this approach. Instead of using the Karhunen-Loève expansion one could also fix a spatial basis and then estimate the polynomial chaos coefficients straight from measurement data; see, for instance, [13, 14].


### 3.1. Karhunen-Loève Expansion

The Karhunen-Loève expansion is a representation of a random field as a linear combination of the eigenfunctions of the associated covariance operator. Truncating the resulting series we obtain a finite-dimensional approximation of the original random field. The Karhunen-Loève expansion is the optimal choice among linear expansions in the sense that it minimizes the mean square truncation error [6].For a fixedω1∈Ω1 let E(ω2,x)=E(ω1,ω2,x) be a random field on D=D(ω1) with mean field(19)E-x=∫Ω2Eω2,xdPω2and a covariance function(20)CEx1,x2=∫Ω2Eω2,x1-E-x1·Eω2,x2-E-x2dPω2that is symmetric and positive definite. The associated covariance operator(21)CEvx1≔∫DCEx1,x2vx2dx2is compact as an operator CE:L2(D)→L2(D) and therefore admits a countable number of positive eigenvalues λmm=1∞ and associated eigenfunctions ϕmm=1∞. The eigenvalues λmm=1∞ tend to zero as m tends to infinity and the eigenfunctions ϕmm=1∞ form an orthonormal basis of L2(D). The eigenpairs may be computed numerically using, for example, fast multipole methods [15].The Karhunen-Loève expansion of the random fieldE(ω2,x) is now given by(22)Eω2,x=E-x+∑m=1∞λmϕmxymω2,where ymm=1∞ are centered and orthogonal but not necessarily independent random variables.In numerical computations we replace the random fieldE(ω2,x) with an approximation obtained by truncating series (22) after M terms. The truncation error depends upon the decay of the Karhunen-Loève eigenvalues and eigenfunctions. See [16] for bounds on this decay.Remark 1.
It is essential for numerical algorithms that the random variablesymm=1∞ are mutually independent, even though this is not in general implied by their orthogonality. A solution to this issue has been proposed in [3], where suitable auxiliary density functions have been introduced. In the numerical examples we assume a representation of the random field with respect to a set of independent uniformly distributed random variables. Other types of random variables (e.g., Gaussian) could be considered as well.


### 3.2. Legendre Chaos

We employ the generalized polynomial chaos framework which essentially means representing our solution on a basis of orthogonal polynomials. In our case we assume the input random variables to be uniformly distributed and the polynomial chaos basis is therefore given by tensorized Legendre polynomials. The use of orthogonal polynomials as a basis allows us to apply stochastic Galerkin based methods and ensures optimal convergence of these methods; see [17]. In this paper we aim to solve the eigenmodes using a method of spectral inverse iteration introduced in [12].Assume thaty(ω2)=(y1(ω2),…,yM(ω2)) is a vector of mutually independent random variables that are uniformly distributed in the interval [-1,1]. The expected value or mean of a random variable v(y) is now given by(23)E2v=∫Ω2vyω2dPω2=∫-1,1Mvy2-Mdy.Here we have used a subscript to indicate that the expected value is taken over the sample space Ω2 only. We associate each multi-index α=(α1,…,αM) with the multivariate Legendre polynomial(24)Λαy≔∏m=1MLαmym,where Lp(x) denotes the univariate Legendre polynomial of degree p. We assume that the polynomials are normalized so that E2[Λα2]=1 for all α∈N0M.The system{Λα(y)∣α∈N0M} forms an orthonormal basis of L2(-1,1M). This allows us to express any square integrable random variable v(y) in a polynomial chaos expansion(25)vy=∑α∈N0MvαΛαy,where the coefficients are given by vα=E2[vΛα].For numerical computations we truncate expansion (25). To this end we choose a finite set of multi-indices A⊂N0M and approximate v(y) with the series(26)vAy=∑α∈AvαΛαy.Here the choice of the set A is key to obtaining accurate approximations and is currently a topic of active research. In our numerical examples we consider sets that for a fixed level L∈R are given by AM,L={α∈N0M∣∑m=1Mαm/λm≤L}, where λmm=1∞ are the Karhunen-Loève eigenvalues. More sophisticated adaptive schemes for selecting the multi-index sets have been presented in [18–20].The most important reason for using the polynomial chaos basis is that it allows fast computation of any expectations involved. This is because we may evaluate integrals over the stochastic domain analytically and are thus able to avoid numerical integration in high dimensions. Orthogonality of the Legendre polynomials yields(27)E2Λα=δα0,E2ΛαΛβ=δαβfor all multi-indices α,β∈N0M. For notational convenience we set(28)cmβγ≔E2ymΛαΛβ,m∈N,α,β∈N0M,cαβγ≔E2ΛαΛβΛγ,α,β,γ∈N0M.These coefficients are easy to evaluate analytically; see the appendix in [12]. It is also worth noting that most of the coefficients are evaluated to zero. This fact can be exploited in order to speed up the solution process.


## 3.1. Karhunen-Loève Expansion

The Karhunen-Loève expansion is a representation of a random field as a linear combination of the eigenfunctions of the associated covariance operator. Truncating the resulting series we obtain a finite-dimensional approximation of the original random field. The Karhunen-Loève expansion is the optimal choice among linear expansions in the sense that it minimizes the mean square truncation error [6].For a fixedω1∈Ω1 let E(ω2,x)=E(ω1,ω2,x) be a random field on D=D(ω1) with mean field(19)E-x=∫Ω2Eω2,xdPω2and a covariance function(20)CEx1,x2=∫Ω2Eω2,x1-E-x1·Eω2,x2-E-x2dPω2that is symmetric and positive definite. The associated covariance operator(21)CEvx1≔∫DCEx1,x2vx2dx2is compact as an operator CE:L2(D)→L2(D) and therefore admits a countable number of positive eigenvalues λmm=1∞ and associated eigenfunctions ϕmm=1∞. The eigenvalues λmm=1∞ tend to zero as m tends to infinity and the eigenfunctions ϕmm=1∞ form an orthonormal basis of L2(D). The eigenpairs may be computed numerically using, for example, fast multipole methods [15].The Karhunen-Loève expansion of the random fieldE(ω2,x) is now given by(22)Eω2,x=E-x+∑m=1∞λmϕmxymω2,where ymm=1∞ are centered and orthogonal but not necessarily independent random variables.In numerical computations we replace the random fieldE(ω2,x) with an approximation obtained by truncating series (22) after M terms. The truncation error depends upon the decay of the Karhunen-Loève eigenvalues and eigenfunctions. See [16] for bounds on this decay.Remark 1.
It is essential for numerical algorithms that the random variablesymm=1∞ are mutually independent, even though this is not in general implied by their orthogonality. A solution to this issue has been proposed in [3], where suitable auxiliary density functions have been introduced. In the numerical examples we assume a representation of the random field with respect to a set of independent uniformly distributed random variables. Other types of random variables (e.g., Gaussian) could be considered as well.


## 3.2. Legendre Chaos

We employ the generalized polynomial chaos framework which essentially means representing our solution on a basis of orthogonal polynomials. In our case we assume the input random variables to be uniformly distributed and the polynomial chaos basis is therefore given by tensorized Legendre polynomials. The use of orthogonal polynomials as a basis allows us to apply stochastic Galerkin based methods and ensures optimal convergence of these methods; see [17]. In this paper we aim to solve the eigenmodes using a method of spectral inverse iteration introduced in [12].Assume thaty(ω2)=(y1(ω2),…,yM(ω2)) is a vector of mutually independent random variables that are uniformly distributed in the interval [-1,1]. The expected value or mean of a random variable v(y) is now given by(23)E2v=∫Ω2vyω2dPω2=∫-1,1Mvy2-Mdy.Here we have used a subscript to indicate that the expected value is taken over the sample space Ω2 only. We associate each multi-index α=(α1,…,αM) with the multivariate Legendre polynomial(24)Λαy≔∏m=1MLαmym,where Lp(x) denotes the univariate Legendre polynomial of degree p. We assume that the polynomials are normalized so that E2[Λα2]=1 for all α∈N0M.The system{Λα(y)∣α∈N0M} forms an orthonormal basis of L2(-1,1M). This allows us to express any square integrable random variable v(y) in a polynomial chaos expansion(25)vy=∑α∈N0MvαΛαy,where the coefficients are given by vα=E2[vΛα].For numerical computations we truncate expansion (25). To this end we choose a finite set of multi-indices A⊂N0M and approximate v(y) with the series(26)vAy=∑α∈AvαΛαy.Here the choice of the set A is key to obtaining accurate approximations and is currently a topic of active research. In our numerical examples we consider sets that for a fixed level L∈R are given by AM,L={α∈N0M∣∑m=1Mαm/λm≤L}, where λmm=1∞ are the Karhunen-Loève eigenvalues. More sophisticated adaptive schemes for selecting the multi-index sets have been presented in [18–20].The most important reason for using the polynomial chaos basis is that it allows fast computation of any expectations involved. This is because we may evaluate integrals over the stochastic domain analytically and are thus able to avoid numerical integration in high dimensions. Orthogonality of the Legendre polynomials yields(27)E2Λα=δα0,E2ΛαΛβ=δαβfor all multi-indices α,β∈N0M. For notational convenience we set(28)cmβγ≔E2ymΛαΛβ,m∈N,α,β∈N0M,cαβγ≔E2ΛαΛβΛγ,α,β,γ∈N0M.These coefficients are easy to evaluate analytically; see the appendix in [12]. It is also worth noting that most of the coefficients are evaluated to zero. This fact can be exploited in order to speed up the solution process.


## 4. Higher Order Finite Elements

Here we give a short overview of thehp-FEM following closely the one in [21]. In the p-version of the FEM the polynomial degree p is used to control the accuracy of the solution while keeping the mesh fixed in contrast to the h-version or standard finite element method, where the polynomial degree is constant but the mesh size varies. The hp-method simply combines the h- and p-refinements.In our setting we can use topologically fixed meshes with high-order elements withp=4 and p=8 to ensure sufficient resolution of the deterministic part of the solution.In the following one way to construct ap-type quadrilateral element is given. The construction of triangles follows similar lines. First of all, the choice of shape functions is not unique. We use the so-called hierarchic integrated Legendre shape functions.Legendre polynomials of degreen can be defined by a recursion formula(29)n+1Pn+1x-2n+1xPnx+nPn-1x=0,where P0(x)=1 and P1(x)=x.The derivatives can similarly be computed by using the recursion(30)1-x2Pn′x=-nxPnx+nPn-1x.The integrated Legendre polynomials are defined forξ∈[-1,1] as(31)ψnξ=2n-12∫-1ξPn-1tdt,n=2,3,…,and can be rewritten as linear combinations of Legendre polynomials(32)ψnξ=122n-1Pnξ-Pn-2ξ,n=2,3,….The normalizing coefficients are chosen so that(33)∫-11dψiξdξdψjξdξdξ=δij,i,j≥2.Using these polynomials we can now define the shape functions for a quadrilateral reference element over the domain[-1,1]×[-1,1]. The shape functions are divided into three categories: nodal shape functions, side modes, and internal modes.There are four nodal shape functions:(34)N1ξ,η=141-ξ1-η,N2ξ,η=141+ξ1-η,N3ξ,η=141+ξ1+η,N4ξ,η=141-ξ1+η,which taken alone define the standard four-node quadrilateral finite element. There are 4(p-1) side modes associated with the sides of a quadrilateral (p≥2), with i=2,…,p,(35)Ni1ξ,η=121-ηψiξ,Ni2ξ,η=121+ξψiη,Ni3ξ,η=121+ηψiη,Ni4ξ,η=121-ξψiξ.For the internal modes we choose the (p-1)(p-1) shape functions:(36)Ni,j0ξ,η=ψiξψjη,i=2,…,p,j=2,…,p.The internal shape functions are often referred to as bubble-functions.Note that some additional book-keeping is necessary. The Legendre polynomials have the propertyPn(-x)=-1nPn(x). This means that every edge must be globally parameterized in the same way in both elements where it belongs.


### 4.1. Curved Boundaries

Since we want to use fixed mesh topologies even with perturbed domains, it is important to represent curved boundary segments accurately. The linear blending function method of Gordon and Hall [22] is our choice for this purpose.In the general case all sides of an element can be curved, but in our case only one side is—as in Figure1. We assume that every side is parameterized:(37)x=xit,y=yit,-1≤t≤1,where i=1,…,number  of  sides. Using capital letters as coordinates of the corner points (Xi,Yi), we can write the mapping for the global x-coordinates of a quadrilateral as(38)x=121-ηx1ξ+121+ξx2η+121+ηx3ξ+121-ξx4η-141-ξ1-ηX1-141+ξ1+ηX2-141+ξ1+ηX3-141-ξ1+ηX4and symmetrically for the y-coordinate. Note that if the side parameterizations represent straight edges, the mapping simplifies to the standard bilinear mapping of quadrilaterals.


## 4.1. Curved Boundaries

Since we want to use fixed mesh topologies even with perturbed domains, it is important to represent curved boundary segments accurately. The linear blending function method of Gordon and Hall [22] is our choice for this purpose.In the general case all sides of an element can be curved, but in our case only one side is—as in Figure1. We assume that every side is parameterized:(37)x=xit,y=yit,-1≤t≤1,where i=1,…,number  of  sides. Using capital letters as coordinates of the corner points (Xi,Yi), we can write the mapping for the global x-coordinates of a quadrilateral as(38)x=121-ηx1ξ+121+ξx2η+121+ηx3ξ+121-ξx4η-141-ξ1-ηX1-141+ξ1+ηX2-141+ξ1+ηX3-141-ξ1+ηX4and symmetrically for the y-coordinate. Note that if the side parameterizations represent straight edges, the mapping simplifies to the standard bilinear mapping of quadrilaterals.


## 5. The Solution Method

We will present a hybrid method of stochastic finite elements for solving the stochastic eigenvalue problem (15). More precisely, we apply stochastic collocation to resolve the dependency on the geometry and the stochastic Galerkin scheme for computing the effects of the random field. The method of stochastic collocation allows for an easy implementation, since the solution statistics are sampled from an ensemble of deterministic solutions. This makes the method attractive from the point of view of random geometry. On the other hand, if a sufficiently accurate parametric representation for the input random field exists, then the use of stochastic Galerkin based methods is well motivated by their high rate of convergence—especially if the number of parameters is large. In this paper we apply a method of spectral inverse iteration introduced in [12] to compute the eigenpair of interest at each collocation point.


### 5.1. Discretization

Let us consider the stochastic eigenproblem (15) for a fixed ω1∈Ω1. We replace the Young modulus E=E(ω1,·,·) with the expression(39)EMω2,x=E-x+∑m=1Mλmϕmxymω2obtained by truncating the Karhunen-Loève expansion (22) after M terms. Plugging (39) into (2) we may write the constitutive matrix (11) in the form(40)Dλω1,ω2,x,μω1,ω2,x=D0ω1,x+∑m=1MDmω1,xymω2.We assume that ymm=1M is a set of mutually independent and uniformly distributed random variables with supports scaled to the interval [-1,1]. The sample space Ω2 is thus parameterized by the vector y=(y1,…,yM)⊂-1,1M. With ω1∈Ω1 being fixed we may now interpret the eigenvalue κ as a function of y and the eigenfunction u as a function of y and x.Next we address the spatial discretization of (15) on the fixed domain D=D(ω1). Let φkk=1N be a set of global basis functions for the approximation space of VD×VD. Here we employ the p-version of the finite element method so that the components of each φk are some shape functions (34), (35), or (36) from Section 4. From (15) we obtain the deterministic Galerkin projection(41)∫DϵvuTD0ϵvφk+∑m=1MymϵvuTDmϵvφkdx=κ∫Du·φkdx,k=1,…,N.We define the stiffness matrix A(m) for each m=0,…,M by setting(42)Aijm=∫DϵvφiTDmϵvφjdxand the mass matrix M via(43)Mij=∫Dφi·φjdx.For a fixed ω1∈Ω1 we have now reduced (15) to a stochastic matrix eigenvalue problem: find a random variable κ(ω1,·):-1,1M→R and a random vector u(ω1,·):-1,1M→RN such that(44)A0+∑m=1MAmymuy=κyMuy.


### 5.2. Stochastic Collocation

In our problem of interest the geometry of the computational domain depends on the random variablesR(ω1) and θ(ω1). Thus, these variables give us a parametrization of the sample space Ω1. We present an anisotropic collocation operator and apply it for computing the dependency of our solution on the parameters. Implementation of a sparse grid collocation operator (see [4, 23]) would also be possible but since we only consider collocation in two dimensions we are restricted to full tensor collocation in the following. In our case we rely on the Gauss-Legendre quadrature, since the random variables R(ω1) and θ(ω1) are uniformly distributed. Furthermore, for notational simplicity we assume these to be scaled so that R,θ∈[-1,1].Denote byzknk=0n the abscissae of the Gauss-Legendre quadrature of order n and by wknk=0n the associated quadrature weights. The quadrature points zknk=0n are the zeros of the univariate Legendre polynomial of degree n+1. The one-dimensional Lagrange interpolation operators In with respect to the interpolation points zknk=0n are defined via(45)Invz=∑k=0nvzknlknzforn≥0.Here lknk=0n are the Lagrange basis polynomials of degree n. For information on orthogonal polynomials and computation of the quadrature weights see [24].The full tensor Lagrange interpolation operator is defined as the tensor product of the univariate operators. In our two-dimensional case this is(46)IγF=In1R⊗In2θ,where In1(R) and In2(θ) denote the one-dimensional interpolation operators in the variables R(ω1) and θ(ω1), respectively, and γ={n1,n2} are the orders of quadrature. Applied on a random variable v(R,θ) the full tensor interpolation gives(47)IγFvR,θ=∑k1=0n1∑k2=0n2vzk1n1,zk2n2lk1n1Rlk2n2θ.The expected value of a random variablev(R,θ) is given by(48)E1v=∫Ω1vRω1,θω1dPω1=∫-1,12vR,θ2-2dRdθ.Here we have again used a subscript to designate that the expectation is taken over the sample space Ω1. From (47) we obtain(49)E1IγFv=∑k1=0n1∑k2=0n2vzk1n1,zk2n2wk1n1wk2n2.Hence we may compute the expected value of the solution once we have solved the underlying problem at each collocation point.


### 5.3. Spectral Inverse Iteration

We choose to solve (44) at each collocation point using a method of spectral inverse iteration. This method has been introduced in [12] and there its convergence towards the eigenpair with the smallest eigenvalue has been verified in the case of an elliptic operator with random coefficients. However, assuming that the eigenpair of interest is of multiplicity one, the method is in fact applicable for a much wider range of problems.We may consider the spectral inverse iteration as a stochastic extension of the deterministic inverse iteration which has been given in Algorithm2 for reference. For the generalized eigenvalue problem of a stiffness matrix A and a mass matrix M, the deterministic inverse iteration converges to the eigenpair for which the eigenvalue is closest to the given parameter κ^.Algorithm 2 (deterministic inverse iteration).
Fixκ^∈R and tol>0. Let u(0) be a normalized initial guess for the eigenvector and set κ(0)=u0TAu(0). For k=1,2,… do the following:(1)
Solve(A-κ^M)v=Mu(k-1) for v.(2)
Setu(k)=v/v.(3)
Setκ(k)=ukTAu(k).(4)
Stop ifu(k)-u(k-1)<tol and return (κ,u)=(κ(k),u(k)).The idea in the spectral inverse iteration is to interpret each equation involved in Algorithm2 in a stochastic sense using Galerkin projections with respect to a given spectral basis. To this end we start by fixing a finite polynomial chaos basis {Λα(y)∣α∈A⊂N0M} and proceed by formulating the respective Galerkin projections. Here we only recap the essentials of the process as it has already been thoroughly explained in [12].For someκ^∈R let us consider the equation(50)Ay-κ^Mvy=Muythat arises from the first step of Algorithm 2. Here we assume the matrix A to take the form(51)A=A0+∑m=1MAmymas in (44). We replace the random vectors u(y) and v(y) with the truncated polynomial chaos expansions(52)uAy=∑β∈AuβΛβy,vAy=∑β∈AvβΛβy.The Galerkin projection of (50) on the basis Λαα∈A is now defined as(53)E2A-κ^MvAΛα=E2MuAΛα,α∈A.Using the orthogonality relations (27) this reduces to the linear system(54)Ψη=ζ,where ζ, η, and Ψ are block vectors and matrices characterized by ζα=Muα,[ηβ]=vβ and(55)Ψαβ=A0-κ^Mδαβ+∑m=1MAmcmαβ.Let us next consider the normalization step in Algorithm2. For a vector v(y) we aim to compute w(y) such that(56)wy=vyvy.If we set n(y)=v(y) then the respective Galerkin projection on the basis Λαα∈A is defined as(57)E2nA2Λα=E2vA2Λα,α∈Aand reduces to(58)∑β∈A∑γ∈Anβnγcαβγ=∑β∈A∑γ∈AvβTMvγcαβγ,α∈A.This is a nonlinear system of equations for the coefficients nαα∈A and can be solved using, for example, Newton’s method with the initial guess nα=vα(y). The Jacobian is given explicitly by(59)Jδϵ=2nϵ∑γ∈Acδϵγ.The Galerkin projection of (56) is now defined as(60)E2wAnAΛα=E2vAΛα,α∈Aand reduces to the block linear system(61)Φθ=η,where ηα=vα, [θβ]=wβ, and(62)Φαβij=δij∑γ∈Anγcαβγ.For computing the eigenvalue we need to be able to evaluate the Rayleigh quotient(63)ry=uyTAyuy,for a vector u(y). As explained in [12] we may approximate r(y) in the Galerkin sense to obtain the formula(64)rα=∑β∈A∑γ∈AuβTA0uγ+∑δ∈A∑m=1MuβTAmuδcmγδcαβγfor the Galerkin coefficients rαα∈A.Algorithm3 now formulates the method of spectral inverse iteration. Similar to the deterministic case, the choice of the parameter κ^ affects the speed of convergence. However, if κ^ is chosen too close to the actual eigenvalue, the iteration might not converge as (50) becomes singular for some y∈-1,1M. We refer to [12] for more information. In our numerical examples we set κ^=0 which ensures convergence to the eigenpair with the smallest eigenvalue as long as the matrix A(y) is positive definite for all y∈-1,1M.Algorithm 3 (spectral inverse iteration for stochastic eigenvalue problems).
Fixκ^∈R and tol>0. Let u(0)=uα0α∈A be a normalized initial guess for the eigenvector and set κ(0)=rαu0α∈A. For k=1,2,… do the following:(1)
Withu=uαk-1α∈A solve system (54) for v=vαα∈A.(2)
Solven=nαα∈A from the nonlinear system (58). Solve system (61) and set u(k)=wαα∈A.(3)
Setκ(k)=rαukα∈A.(4)
Stop if∑α∈Auα(k)-uα(k-1)<tol and return (κ,u)=(κ(k),u(k)).


### 5.4. Response Statistics

We may easily calculate statistical moments of a random variablev(R,θ,y) once we know its polynomial chaos expansion(65)vAR,θ,y=∑α∈AvαR,θΛαyat each collocation point. By using the orthogonality relations (27) we obtain(66)E2vAR,θ,·=v0R,θ,E2vAR,θ,·2=∑α∈AvαR,θ2.The first and second moments of v(R,θ,y) are now given by(67)Ev=E1E2v,Ev2=E1E2v2.These can be computed by applying formula (49) on (66). From these we obtain approximations for the expected value and variance of v(R,θ,y).


## 5.1. Discretization

Let us consider the stochastic eigenproblem (15) for a fixed ω1∈Ω1. We replace the Young modulus E=E(ω1,·,·) with the expression(39)EMω2,x=E-x+∑m=1Mλmϕmxymω2obtained by truncating the Karhunen-Loève expansion (22) after M terms. Plugging (39) into (2) we may write the constitutive matrix (11) in the form(40)Dλω1,ω2,x,μω1,ω2,x=D0ω1,x+∑m=1MDmω1,xymω2.We assume that ymm=1M is a set of mutually independent and uniformly distributed random variables with supports scaled to the interval [-1,1]. The sample space Ω2 is thus parameterized by the vector y=(y1,…,yM)⊂-1,1M. With ω1∈Ω1 being fixed we may now interpret the eigenvalue κ as a function of y and the eigenfunction u as a function of y and x.Next we address the spatial discretization of (15) on the fixed domain D=D(ω1). Let φkk=1N be a set of global basis functions for the approximation space of VD×VD. Here we employ the p-version of the finite element method so that the components of each φk are some shape functions (34), (35), or (36) from Section 4. From (15) we obtain the deterministic Galerkin projection(41)∫DϵvuTD0ϵvφk+∑m=1MymϵvuTDmϵvφkdx=κ∫Du·φkdx,k=1,…,N.We define the stiffness matrix A(m) for each m=0,…,M by setting(42)Aijm=∫DϵvφiTDmϵvφjdxand the mass matrix M via(43)Mij=∫Dφi·φjdx.For a fixed ω1∈Ω1 we have now reduced (15) to a stochastic matrix eigenvalue problem: find a random variable κ(ω1,·):-1,1M→R and a random vector u(ω1,·):-1,1M→RN such that(44)A0+∑m=1MAmymuy=κyMuy.


## 5.2. Stochastic Collocation

In our problem of interest the geometry of the computational domain depends on the random variablesR(ω1) and θ(ω1). Thus, these variables give us a parametrization of the sample space Ω1. We present an anisotropic collocation operator and apply it for computing the dependency of our solution on the parameters. Implementation of a sparse grid collocation operator (see [4, 23]) would also be possible but since we only consider collocation in two dimensions we are restricted to full tensor collocation in the following. In our case we rely on the Gauss-Legendre quadrature, since the random variables R(ω1) and θ(ω1) are uniformly distributed. Furthermore, for notational simplicity we assume these to be scaled so that R,θ∈[-1,1].Denote byzknk=0n the abscissae of the Gauss-Legendre quadrature of order n and by wknk=0n the associated quadrature weights. The quadrature points zknk=0n are the zeros of the univariate Legendre polynomial of degree n+1. The one-dimensional Lagrange interpolation operators In with respect to the interpolation points zknk=0n are defined via(45)Invz=∑k=0nvzknlknzforn≥0.Here lknk=0n are the Lagrange basis polynomials of degree n. For information on orthogonal polynomials and computation of the quadrature weights see [24].The full tensor Lagrange interpolation operator is defined as the tensor product of the univariate operators. In our two-dimensional case this is(46)IγF=In1R⊗In2θ,where In1(R) and In2(θ) denote the one-dimensional interpolation operators in the variables R(ω1) and θ(ω1), respectively, and γ={n1,n2} are the orders of quadrature. Applied on a random variable v(R,θ) the full tensor interpolation gives(47)IγFvR,θ=∑k1=0n1∑k2=0n2vzk1n1,zk2n2lk1n1Rlk2n2θ.The expected value of a random variablev(R,θ) is given by(48)E1v=∫Ω1vRω1,θω1dPω1=∫-1,12vR,θ2-2dRdθ.Here we have again used a subscript to designate that the expectation is taken over the sample space Ω1. From (47) we obtain(49)E1IγFv=∑k1=0n1∑k2=0n2vzk1n1,zk2n2wk1n1wk2n2.Hence we may compute the expected value of the solution once we have solved the underlying problem at each collocation point.


## 5.3. Spectral Inverse Iteration

We choose to solve (44) at each collocation point using a method of spectral inverse iteration. This method has been introduced in [12] and there its convergence towards the eigenpair with the smallest eigenvalue has been verified in the case of an elliptic operator with random coefficients. However, assuming that the eigenpair of interest is of multiplicity one, the method is in fact applicable for a much wider range of problems.We may consider the spectral inverse iteration as a stochastic extension of the deterministic inverse iteration which has been given in Algorithm2 for reference. For the generalized eigenvalue problem of a stiffness matrix A and a mass matrix M, the deterministic inverse iteration converges to the eigenpair for which the eigenvalue is closest to the given parameter κ^.Algorithm 2 (deterministic inverse iteration).
Fixκ^∈R and tol>0. Let u(0) be a normalized initial guess for the eigenvector and set κ(0)=u0TAu(0). For k=1,2,… do the following:(1)
Solve(A-κ^M)v=Mu(k-1) for v.(2)
Setu(k)=v/v.(3)
Setκ(k)=ukTAu(k).(4)
Stop ifu(k)-u(k-1)<tol and return (κ,u)=(κ(k),u(k)).The idea in the spectral inverse iteration is to interpret each equation involved in Algorithm2 in a stochastic sense using Galerkin projections with respect to a given spectral basis. To this end we start by fixing a finite polynomial chaos basis {Λα(y)∣α∈A⊂N0M} and proceed by formulating the respective Galerkin projections. Here we only recap the essentials of the process as it has already been thoroughly explained in [12].For someκ^∈R let us consider the equation(50)Ay-κ^Mvy=Muythat arises from the first step of Algorithm 2. Here we assume the matrix A to take the form(51)A=A0+∑m=1MAmymas in (44). We replace the random vectors u(y) and v(y) with the truncated polynomial chaos expansions(52)uAy=∑β∈AuβΛβy,vAy=∑β∈AvβΛβy.The Galerkin projection of (50) on the basis Λαα∈A is now defined as(53)E2A-κ^MvAΛα=E2MuAΛα,α∈A.Using the orthogonality relations (27) this reduces to the linear system(54)Ψη=ζ,where ζ, η, and Ψ are block vectors and matrices characterized by ζα=Muα,[ηβ]=vβ and(55)Ψαβ=A0-κ^Mδαβ+∑m=1MAmcmαβ.Let us next consider the normalization step in Algorithm2. For a vector v(y) we aim to compute w(y) such that(56)wy=vyvy.If we set n(y)=v(y) then the respective Galerkin projection on the basis Λαα∈A is defined as(57)E2nA2Λα=E2vA2Λα,α∈Aand reduces to(58)∑β∈A∑γ∈Anβnγcαβγ=∑β∈A∑γ∈AvβTMvγcαβγ,α∈A.This is a nonlinear system of equations for the coefficients nαα∈A and can be solved using, for example, Newton’s method with the initial guess nα=vα(y). The Jacobian is given explicitly by(59)Jδϵ=2nϵ∑γ∈Acδϵγ.The Galerkin projection of (56) is now defined as(60)E2wAnAΛα=E2vAΛα,α∈Aand reduces to the block linear system(61)Φθ=η,where ηα=vα, [θβ]=wβ, and(62)Φαβij=δij∑γ∈Anγcαβγ.For computing the eigenvalue we need to be able to evaluate the Rayleigh quotient(63)ry=uyTAyuy,for a vector u(y). As explained in [12] we may approximate r(y) in the Galerkin sense to obtain the formula(64)rα=∑β∈A∑γ∈AuβTA0uγ+∑δ∈A∑m=1MuβTAmuδcmγδcαβγfor the Galerkin coefficients rαα∈A.Algorithm3 now formulates the method of spectral inverse iteration. Similar to the deterministic case, the choice of the parameter κ^ affects the speed of convergence. However, if κ^ is chosen too close to the actual eigenvalue, the iteration might not converge as (50) becomes singular for some y∈-1,1M. We refer to [12] for more information. In our numerical examples we set κ^=0 which ensures convergence to the eigenpair with the smallest eigenvalue as long as the matrix A(y) is positive definite for all y∈-1,1M.Algorithm 3 (spectral inverse iteration for stochastic eigenvalue problems).
Fixκ^∈R and tol>0. Let u(0)=uα0α∈A be a normalized initial guess for the eigenvector and set κ(0)=rαu0α∈A. For k=1,2,… do the following:(1)
Withu=uαk-1α∈A solve system (54) for v=vαα∈A.(2)
Solven=nαα∈A from the nonlinear system (58). Solve system (61) and set u(k)=wαα∈A.(3)
Setκ(k)=rαukα∈A.(4)
Stop if∑α∈Auα(k)-uα(k-1)<tol and return (κ,u)=(κ(k),u(k)).


## 5.4. Response Statistics

We may easily calculate statistical moments of a random variablev(R,θ,y) once we know its polynomial chaos expansion(65)vAR,θ,y=∑α∈AvαR,θΛαyat each collocation point. By using the orthogonality relations (27) we obtain(66)E2vAR,θ,·=v0R,θ,E2vAR,θ,·2=∑α∈AvαR,θ2.The first and second moments of v(R,θ,y) are now given by(67)Ev=E1E2v,Ev2=E1E2v2.These can be computed by applying formula (49) on (66). From these we obtain approximations for the expected value and variance of v(R,θ,y).


## 6. Numerical Experiments

Let us consider a problem with uncertain domain and uncertain field. The computational domain isD=[-2,2]×[-2,2]∖C, where the cavity C=[-R,R]×[-R,R]∪annulus  with  diameter=1/4 which is rotated θ radians by the origin. Here R=U([1/2,3/5]) and θ=U([0,π/12]) are random with uniform distribution. The bottom edge y=-2 is clamped; that is, all displacements are inhibited.In Figure2 the effect of domain perturbation on the modes is illustrated. Notice that if one wanted to compute the statistics of the modes, it would be necessary to map every realisation onto the nominal domain via conformal mappings, for instance.Figure 2
Case A: effect of perturbation on eigenmodes. Contour lines of displacement.
(a)
Contour plot:u1 (R=4/5)
(b)
Contour plot:u2 (R=4/5)
(c)
Contour plot:u1 (R=4/5, θ=π/12)
(d)
Contour plot:u2 (R=4/5, θ=π/12)The (constant) material parameters are Poisson ratioν=1/3 and the mean field of the Young modulus E-=2.069·1011 MPa. We have used the convention, however, where the reported results are normalized by the value of E-.We will proceed in two steps by first letting only the geometry vary before adding uncertainty to the field as well. The cases are(A)
uncertain domain with deterministic field, solved with collocation,(B)
uncertain domain with uncertain field, solved with the new hybrid method.In both cases the nominal domain is the same (see Figure 1). In the absence of exact solutions overkill solutions have been used in convergence graphs unless otherwise specified. We have also verified the results by comparing to a full Monte Carlo simulation of 620000 draws.


### 6.1. Case A

The values of the smallest eigenvalue over the parameter space are shown in Figure3 using both contour and surface plots.Figure 3
Case A: eigenvalues over the parameter space: the standard deviation isσ=0.0008.
(a)
Contour plot:R (x-axis) versus θ (y-axis)
(b)
Surface plot:R (x-axis) versus θ (y-axis)We consider a collocation sequence, where the collocation points are taken to span all 16 cases from1×1 to 4×4 points over the R×θ-parameter space. The relative convergence of the smallest eigenvalue over the whole collocation sequence is shown in Figure 4. We have used constant polynomial order p=8 (1728 degrees of freedom) throughout and p=16 (6528 degrees of freedom) for the reference computation.Figure 4
Case A: convergence in eigenvalue: relative error versus collocation sequences in log-scale. Dashed line: fixed deterministic discretization in reference. Solid line: higher resolution reference. Plotmarkers: uniform collocation grids.
(a)
Expected value
(b)
VarianceThe dashed line shows the convergence in the case when the deterministic discretization of the reference solution is used also in collocation. The exponential convergence in expectation is due to convergence in parameter space only. The poorly performing grids in Figure4 are the anisotropic ones with more points for θ. This is due to relative effects of R and θ not being exactly balanced (see Figure 3). The solid line indicates that already at rule (2,2) the spatial discretization used is not sufficient when the reference has been computed with much higher resolution. The convergence stalls completely since the deterministic error dominates. Convergence in variance is slower, as suggested by theory.


### 6.2. Case B

We proceed to consider a fixed collocation grid with4×4 collocation points for the parameters R and θ. Assuming normalization of the mean field (E-=1) we set the coefficients in the Karhunen-Loève expansion of the Young modulus (39) to be ϕm(x)=cos⁡(mπd(x)), where d(x) is the distance from x to the inner boundary and λm=m+1-3. Thus the series can be considered to converge at an algebraic rate of s=-1.5. We set the truncation parameter to be M=3. Using p=4 and 25 multi-indices, every linear system has 12000 degrees of freedom.At each collocation point we employ the spectral inverse iteration given in Algorithm3 with κ^=0. For the polynomial chaos basis we consider multi-index sets of the form AM,L={α∈N0M∣∑m=1Mαm/λm≤L}. Varying the level L∈R results in multi-index sets of different size. The overkill solution used as a reference has been computed with 64 multi-indices. Convergence of the relative error of the smallest eigenvalue towards the overkill solution as a function of the size of the multi-index set is illustrated in Figure 5. We observe algebraic rates of convergence for both the expected value and the variance.Figure 5
Case B: convergence in eigenvalue in the case of the4×4 collocation grid: size of the multi-index sets versus relative error in the hybrid solution in loglog-scale. Plotmarkers: multi-index sets of different size. Dashed line: error between the overkill solution and the Monte Carlo solution.
(a)
Expected value
(b)
VarianceThe dashed horizontal line in Figure5 illustrates the relative error between the overkill solution and the Monte Carlo solution. The statistics for the eigenvalue are(68)Eκ≈1.64286·10-2,Var⁡κ≈9.19957·10-7for the overkill solution and(69)Eκ≈1.64236·10-2,Var⁡κ≈9.11820·10-7for the Monte Carlo solution. The central limit theorem suggests that the standard deviation of the Monte Carlo mean is 1.21271·10-6, which is of the same size range as the error between the overkill solution and the Monte Carlo solution. Comparison with standard collocation is more complicated, since the hybrid method includes a collocation part. For an isotropic grid aligned for the geometric uncertainty dimensions, (R,θ,M1,M2,M3)=(4,4,8,8,8), that is, 8 points per every term in the Karhunen-Loève expansion, we get statistics:(70)Eκ≈1.64237·10-2,Var⁡κ≈9.12691·10-7within expected accuracy as well.Algorithm3 also gives us the eigenvector associated with the smallest eigenvalue. More precisely we obtain the spectral components uαα∈AM,L at each collocation point. In Figure 6 we have illustrated the means of the sorted L2 and H1 energies of these spectral components. At each of the 16 collocation points we would observe results similar to those presented in Figure 6. Theory suggests that the asymptotic convergence rates for these energies should be r1=s-1/2 and r2=s-1 for the L2 norms and H1 seminorms, respectively [25, 26]. In our example we observe much faster convergence than these asymptotic rates would predict. However, the observed convergence rates differ by one half, which is in accordance with the theory.Figure 6
Case B: energies of the spectral componentsuαα∈AM,L of the eigenfunction in the case of the 4×4 collocation grid. Dashed line: algebraic convergence rates of −3.8 and −3.3 for the L2 norm and H1 seminorm, respectively. Plotmarkers: sorted means of the energies in loglog-scale.
(a)
L
2 energies
(b)
H
1 energiesThe Monte Carlo reference results have been computed from 620000 draws using 432 single core CPU hours on Intel Xeon (2009, 2.8 GHz). The anisotropic collocation took two and the hybrid algorithm for the highest reported number of multi-indices (64) exactly one single core CPU hour on the same machine. Naturally, in iterative algorithms there are many tolerances with which to tune the performance, but these reported times can be taken as representative.Remark 4.
Notice that in more elaborate models for geometric uncertainty with much higher dimensions of the parameter space the asymptotic convergence rates will apply to the collocation method. For lower dimensional cases it is possible to converge at faster rates, of course, but in real-life applications one should consider the modeling error as well. Similarly, for the Galerkin scheme the size of multi-index sets is too small for the energies of the spectral components to exhibit asymptotic behaviour.


## 6.1. Case A

The values of the smallest eigenvalue over the parameter space are shown in Figure3 using both contour and surface plots.Figure 3
Case A: eigenvalues over the parameter space: the standard deviation isσ=0.0008.
(a)
Contour plot:R (x-axis) versus θ (y-axis)
(b)
Surface plot:R (x-axis) versus θ (y-axis)We consider a collocation sequence, where the collocation points are taken to span all 16 cases from1×1 to 4×4 points over the R×θ-parameter space. The relative convergence of the smallest eigenvalue over the whole collocation sequence is shown in Figure 4. We have used constant polynomial order p=8 (1728 degrees of freedom) throughout and p=16 (6528 degrees of freedom) for the reference computation.Figure 4
Case A: convergence in eigenvalue: relative error versus collocation sequences in log-scale. Dashed line: fixed deterministic discretization in reference. Solid line: higher resolution reference. Plotmarkers: uniform collocation grids.
(a)
Expected value
(b)
VarianceThe dashed line shows the convergence in the case when the deterministic discretization of the reference solution is used also in collocation. The exponential convergence in expectation is due to convergence in parameter space only. The poorly performing grids in Figure4 are the anisotropic ones with more points for θ. This is due to relative effects of R and θ not being exactly balanced (see Figure 3). The solid line indicates that already at rule (2,2) the spatial discretization used is not sufficient when the reference has been computed with much higher resolution. The convergence stalls completely since the deterministic error dominates. Convergence in variance is slower, as suggested by theory.


## 6.2. Case B

We proceed to consider a fixed collocation grid with4×4 collocation points for the parameters R and θ. Assuming normalization of the mean field (E-=1) we set the coefficients in the Karhunen-Loève expansion of the Young modulus (39) to be ϕm(x)=cos⁡(mπd(x)), where d(x) is the distance from x to the inner boundary and λm=m+1-3. Thus the series can be considered to converge at an algebraic rate of s=-1.5. We set the truncation parameter to be M=3. Using p=4 and 25 multi-indices, every linear system has 12000 degrees of freedom.At each collocation point we employ the spectral inverse iteration given in Algorithm3 with κ^=0. For the polynomial chaos basis we consider multi-index sets of the form AM,L={α∈N0M∣∑m=1Mαm/λm≤L}. Varying the level L∈R results in multi-index sets of different size. The overkill solution used as a reference has been computed with 64 multi-indices. Convergence of the relative error of the smallest eigenvalue towards the overkill solution as a function of the size of the multi-index set is illustrated in Figure 5. We observe algebraic rates of convergence for both the expected value and the variance.Figure 5
Case B: convergence in eigenvalue in the case of the4×4 collocation grid: size of the multi-index sets versus relative error in the hybrid solution in loglog-scale. Plotmarkers: multi-index sets of different size. Dashed line: error between the overkill solution and the Monte Carlo solution.
(a)
Expected value
(b)
VarianceThe dashed horizontal line in Figure5 illustrates the relative error between the overkill solution and the Monte Carlo solution. The statistics for the eigenvalue are(68)Eκ≈1.64286·10-2,Var⁡κ≈9.19957·10-7for the overkill solution and(69)Eκ≈1.64236·10-2,Var⁡κ≈9.11820·10-7for the Monte Carlo solution. The central limit theorem suggests that the standard deviation of the Monte Carlo mean is 1.21271·10-6, which is of the same size range as the error between the overkill solution and the Monte Carlo solution. Comparison with standard collocation is more complicated, since the hybrid method includes a collocation part. For an isotropic grid aligned for the geometric uncertainty dimensions, (R,θ,M1,M2,M3)=(4,4,8,8,8), that is, 8 points per every term in the Karhunen-Loève expansion, we get statistics:(70)Eκ≈1.64237·10-2,Var⁡κ≈9.12691·10-7within expected accuracy as well.Algorithm3 also gives us the eigenvector associated with the smallest eigenvalue. More precisely we obtain the spectral components uαα∈AM,L at each collocation point. In Figure 6 we have illustrated the means of the sorted L2 and H1 energies of these spectral components. At each of the 16 collocation points we would observe results similar to those presented in Figure 6. Theory suggests that the asymptotic convergence rates for these energies should be r1=s-1/2 and r2=s-1 for the L2 norms and H1 seminorms, respectively [25, 26]. In our example we observe much faster convergence than these asymptotic rates would predict. However, the observed convergence rates differ by one half, which is in accordance with the theory.Figure 6
Case B: energies of the spectral componentsuαα∈AM,L of the eigenfunction in the case of the 4×4 collocation grid. Dashed line: algebraic convergence rates of −3.8 and −3.3 for the L2 norm and H1 seminorm, respectively. Plotmarkers: sorted means of the energies in loglog-scale.
(a)
L
2 energies
(b)
H
1 energiesThe Monte Carlo reference results have been computed from 620000 draws using 432 single core CPU hours on Intel Xeon (2009, 2.8 GHz). The anisotropic collocation took two and the hybrid algorithm for the highest reported number of multi-indices (64) exactly one single core CPU hour on the same machine. Naturally, in iterative algorithms there are many tolerances with which to tune the performance, but these reported times can be taken as representative.Remark 4.
Notice that in more elaborate models for geometric uncertainty with much higher dimensions of the parameter space the asymptotic convergence rates will apply to the collocation method. For lower dimensional cases it is possible to converge at faster rates, of course, but in real-life applications one should consider the modeling error as well. Similarly, for the Galerkin scheme the size of multi-index sets is too small for the energies of the spectral components to exhibit asymptotic behaviour.


## 7. Conclusions

We have demonstrated the practicality of the new hybrid algorithm in a practical, yet idealized, application. The combination of collocation for the geometric uncertainty and Galerkin for that of materials allows us to take advantage of the higher computational efficiency of the Galerkin approach whilst keeping mesh generation simple. Indeed, in this paper we have kept the mesh topologically fixed. The observed convergence rates for the smallest eigenvalue are in line with theoretical predictions.The main two remaining issues for future research are the mapping of the eigenmodes to the nominal domain and handling of double eigenvalues in the general case. In 2D we are confident that this extension can be readily done, but in 3D there are still open mathematical questions on how to connect the nominal domain and the perturbed ones.



---

*Source: 812069-2015-07-27.xml*