- Email: [email protected]
A multilevel finite element method for Fredholm integral eigenvalue problems
Accepted Manuscript A multilevel finite element method for Fredholm integral eigenvalue problems
Hehu Xie, Tao Zhou
PII: DOI: Reference:
S0021-9991(15)00642-7 http://dx.doi.org/10.1016/j.jcp.2015.09.043 YJCPH 6146
To appear in:
Journal of Computational Physics
Received date: Revised date: Accepted date:
9 November 2014 16 July 2015 24 September 2015
Please cite this article in press as: H. Xie, T. Zhou, A multilevel finite element method for Fredholm integral eigenvalue problems, J. Comput. Phys. (2015), http://dx.doi.org/10.1016/j.jcp.2015.09.043
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A multilevel finite element method for Fredholm integral eigenvalue problems Hehu Xie∗ and Tao Zhou†
Abstract In this work, we proposed a multigrid finite element (MFE) method for solving the Fredholm integral eigenvalue problems. The main motivation for such studies is to compute the Karhunen-Lo`eve expansions of random fields, which plays an important role in the applications of uncertainty quantification. In our MFE framework, solving the eigenvalue problem is converted to doing a series of integral iterations and eigenvalue solving in the coarsest mesh. Then, any exist efficient integration scheme can be used for the associated integration process. The error estimates are provided, and the computational complexity is analyzed. It is noticed that the total computational work of our method is comparable with a single integration step in the finest mesh. Several numerical experiments are presented to validate the efficiency of the proposed numerical method. Keywords. Uncertainty quantification, Karhunen-L`oeve expansion, Fredholm eigenvalue problem, multigrid finite element. AMS subject classifications. 65N30, 65N25, 65L15, 65B99.
1
Introduction
Eigenvalue problems involving the Fredholm integral equation arises in many areas of sciences and engineering (see e.g. [6, 9, 11, 13] and references therein). Here, we focus on its particular application for Karhunen-L`oeve approximation of random fields, which plays an important role in uncertainty quantification (UQ) [2, 10]. In UQ, a typical problem is to numerically solve a PDE with random field coefficients. Then, to transform the original problem into a PDE with finite random parameters, one usually adopts the truncated Karhunen-L`oeve expansion of the random field ∗
LSEC, NCMIS, Institute of Computational Mathematics, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China. Email: [email protected]. † LSEC, Institute of Computational Mathematics, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China. Email: [email protected].
1
[16,32]. In such framework, the random field is represented in terms of uncorrelated random variables and eigenfunctions of a homogeneous Fredholm integral equation of the second kind, whose kernel is the covariance function of the random field. Thus, to efficiently solve the original problem, a key issue is how to efficiently build the truncated Karhunen-L`oeve expansion, by solving a Fredholm integral eigenvalue problem. Despite the large amount of the literature on numerical methods for Fredholm integral equations (see e.g. [1, 7] and references therein), there are only a few numerical results for Fredholm eigenvalue problems (see e.g. [12, 22, 25]). From the UQ application point of view (especially for the computations of Karhunen-L`oeve approximations), one can refer to [21] for a finite element method combined with the fast multipole methods, to [19] for a wavelet Galerkin arpproach using Haar basis functions, to [18] for the finite cell method, to [17] for a spectral element approximation, to name a few. In this paper, we propose a multigrid finite element method for solving the Fredholm integral eigenvalue problem with the idea coming from [15, 27, 28]. In our framework, solving the eigenvalue problem is converted to a series of integral processes and eigenvalue solving on the coarsest mesh. Then, any exist efficient numerical integral scheme can be used for the integral processes. The corresponding error estimates are given, and the computational complexity is analyzed. It is noticed that the total computational work of our method is comparable with the integral scheme. Several numerical experiments are presented to validate the efficiency of the method. The rest of the paper is organized as follows. In section 2, we give a brief review of the Karhunen-L`oeve expansion of random fields, and introduce its applications in UQ. In section 3, a multigrid finite element method is proposed to solve the corresponding Fredholm integral eigenvalue problem, moreover, error estimates and complexity analysis will be provided. Several numerical examples are given in Section 4 to illustrate the efficiency of the proposed method. We finally give some conclusions in section 5.
2
Karhunen-L` oeve expansion of random fields
In this section, we shall give a brief review of the Karhunen-L`oeve expansion of random fields, and introduce its applications in UQ. To this end, let (Ω, Σ, P) be a probability space, where Ω is the space of events, Σ the corresponding σ-algebra and P the probability measure. A spatial dependent random field is a P -measurable map a(x, ω) : Ω → L∞ (D) (where D ⊂ Rd is the physical domain). The mean of the random field a(x, ω) is defined by a ¯(x) = E[a(x, ·)] = a(x, ω)dP(ω). (2.1) Ω
2
We assume that all the random fields considered in this paper are square integrable (i.e., the so-called second order random fields), that is, E[a2 (x, ·)]dx < +∞. (2.2) D
We define the covariance function of a space-dependent random field a(x, ω) as ¯(x)) (a(y, ·) − a ¯(y)) , x, y ∈ D. (2.3) Va (x, y) = E (a(x, ·) − a Given such a random field a(x, ω), prediction of for example a concentration φ(x, ω) may requires solution of a stochastic partial differential equation such as (here we take the elliptic problem as an example) −∇ · (a(x, ω)∇x φ(x, ω)) = f (x),
∀x ∈ D.
(2.4)
It is well known that any second order random field a(x, ω), with continuous and positive definite covariance function Va : D × D → R, can be represented as an infinite sum of random variables, by means of the Karhunen-L`oeve expansion [10]. To this end, we introduce the compact and self-adjoint operator Ka : L2 (D) → L2 (D), which is defined by Va (x, ·)u(x)dx, ∀v ∈ L2 (D). (2.5) Ka u(·) = D
Then, consider the sequence of non-negative decreasing eigenvalues of Ka , {λi }∞ i=1 , ∞ and the corresponding sequence of orthonormal eigenfunctions, {vi }i=1 , satisfying ui uj dx = δij for i, j ∈ N + , (2.6) Ka ui = λi ui , D
where δij is the Kronecker symbol. In addition, define the mutually uncorrelated real random variables 1 Yi (ω) := √ (2.7) a(x, ω) − a ¯(x) ui (x)dx, i = 1, 2, ... λi D with zero mean and unit variance, i.e. E[Yi ] = 0 and E[Yi Yj ] = δij for i, j ∈ N + . Then, the N -term truncated Karhunen-L`oeve expansion of the random field a(x, ω), which we denote by aN , is defined by ¯(x) + aN (x, ω) = a
N λi Yi (ω)ui (x), ∀N ∈ N + .
(2.8)
i=1
By Mercer’s theorem [20], it follows that lim sup E
N →∞ x∈D
a(x, ·) − aN (x, ·)
2
= lim sup N →∞ x∈D
3
∞ i=N +1
λi u2i (x) = 0.
(2.9)
Note that such N -term Karhunen-L`oeve truncation is optimal in the mean square sense [10]. By using the truncated Karhunen-L`oeve expansion, one converts the original stochastic elliptic problem (2.4) into the following elliptic problem with N random parameters: −∇ · (a(x, Y)∇x φ(x, Y)) = f (x),
∀x ∈ D,
(2.10)
with Y being the random vector Y = (Y1 , ...YN ). For the above finite dimensional stochastic problem, many approaches have been developed to get the quantity of interests of the solution φ(x, Y)). For example, the stochastic Galerkin method which proposed originally in [10, 26] and further developed in [30], the stochastic collocation method [16, 24, 29], the polynomial regression [23, 33], etc. By the discussions above, a key issue is to efficiently get the truncated KarhunenL`oeve expansion, i.e., efficiently solve the corresponding Fredholm integral eigenvalue problem: Ka u(x) = Va (x, y)u(y)dy = λu(x), ∀x ∈ D. (2.11) D
The purpose of this work is to propose a multigrid finite element method for solving (2.11). To this end, we remark that regularity results of the eigenfunctions and decay properties of the eigenvalues have been investigated in [21]. Particularly, it is proved in [21] that the eigenfunctions inherited the regularity of the corresponding covariance function, and furthermore, the decay properties of the eigenvalues is also depend on the regularity of the covariance function. More precisely, we have the following propositions [21]: Proposition 2.1. Let Va ∈ L2 (D × D) be a symmetric kernel defining the compact and nonnegative integral operator 2 2 Ka : L (D) → L (D), Ka u(x) = Va (x, y)u(y)dy. (2.12) D
If Va is piecewise analytic (see a detailed definition in [21]) on D × D and {λm }m≥1 is the eigenvalue sequence of its associated operator Ka , then there exist constants c1 , c2 > 0 depending only on Va such that 0 ≤ λm ≤ c1 e−c2 m
1/d
,
∀m ≥ 1.
If Va is piecewise H k,0 := H k ⊗ L2 on D × D with k > 0, then there exists a constant c3 > 0 depending only on Va such that λm ≤ c3 m−k/d ,
∀m ≥ 1.
Proposition 2.2. Assume that the covariance kernel Va (x, y) of a random field a(x, ω) is piecewise analytic (piecewise H p,q ) on D × D (for rigorous definitions, see 4
[21]). Then the Karhunen-L`oeve eigenfunctions {vm (x)}m≥1 of a(x, ω) are piecewise analytic (piecewise H p ). Moreover, for D ⊂ Rd a bounded domain and Va piecewise smooth on D ×D, such that the each piece Dj have the uniform cone property, we denote by {λm , vm }m≥1 the eigenpair sequence of the associated integral operator Ka via (2.12), such that
vm L2 = 1 for m ≥ 1. Then for any s > 0 and any multi-index α ∈ N d there exists cV,s,α > 0 such that
∂ α vm L∞ (Dj ) ≤ cV,s,α λ−s m ,
∀m ≥ 1.
For more properties of the Karhunen-L`oeve expansion of random fields, please refer to [21]. We also remark that the Karhunen-L`oeve expansion is not the only way to approximate random fields, other approaches can also be considered, see e.g. [5] for a Fourier-based decomposition, where the authors uses trigonometric polynomials as basis functions in the physical space.
3
Fredholm integral eigenvalue problem and its finite element discretization
In this section, we shall introduce a multigrid finite element (MFE) approach for solving the following Fredholm integral eigenvalue problems (FIEPs): Ka u = λu,
(3.1)
where the integral operator Ka is defined as Ka u := Va (x, y)u(y)dy,
(3.2)
D
with Va (x, y) being some covariance functions of certain random fields. We denote by V the space L2 (D) in this paper. The weak form of (3.1) is: Find (λ, u) ∈ R × V such that a(u, u) = 1, where
and λa(u, v) = b(u, v),
(3.3)
uvdx,
a(u, v) =
v ∈ V,
Ka uvdx.
b(u, v) =
D
D
We define the corresponding norm · a and semi-norm · b as:
v 2a := a(v, v),
v 2b := b(v, v). 5
(3.4)
For the eigenvalue λ, there exists the following Rayleigh quotient expression (see, e.g., [3, 4, 31]) λ=
b(u, u) . a(u, u)
(3.5)
We shall use the finite element discretizition for problem (3.3). To do this, let us define the finite element approximations of the problem (3.3). First we generate a shape-regular decomposition of the computing domain D ⊂ Rd (d = 2, 3) into triangles or rectangles for d = 2 (tetrahedrons or hexahedrons for d = 3). The diameter of a cell K ∈ Th is denoted by hK . The mesh diameter h describes the maximum diameter of all cells K ∈ Th . Based on the mesh Th , we construct a finite element space Vh ⊂ V as a piecewise polynomial space. Then we define the approximation for the eigenpair (λ, u) of (3.3) by the finite element method as: ¯ h , u¯h ) ∈ R × Vh such that Find (λ a(¯ uh , u¯h ) = 1,
¯ h a(¯ and λ uh , vh ) = b(¯ uh , vh ),
∀vh ∈ Vh .
(3.6)
For (3.6), the following Rayleigh quotient expression for λh holds [3, 4, 31]: uh , u¯h ) ¯ h = b(¯ λ . a(¯ uh , u¯h )
(3.7)
We assume that Vh ⊂ V is a family of finite-dimensional spaces that satisfy the following assumption: lim inf w − v a = 0,
h→0 v∈Vh
∀w ∈ V.
(3.8)
Let Ph be the finite element projection operator of V onto Vh that defined by a(w − Ph w, vh ) = 0,
∀w ∈ V and ∀vh ∈ Vh .
(3.9)
Obviously, we have
Ph w a ≤ w a ,
∀w ∈ V.
(3.10)
as h → 0.
(3.11)
For any w ∈ V , by (3.8) we have
w − Ph w a → 0, Define ηa (h) as ηa (h) =
sup
inf T f − v a ,
f ∈V,f a =1 v∈Vh
6
(3.12)
where the operator T : V → V is defined as ∀f ∈ V and ∀v ∈ V.
a(T f, v) = b(f, v),
(3.13)
As we know, the operator T : V → V is compact. Then the results about the eigenvalues of compact operators can be used here. In order to derive the error estimate of eigenpair approximation in the weak semi-norm · b , we need the following error estimates of the finite element projection operator Ph . Lemma 3.1. ( [3, Lemma 3.3 and Lemma 3.4]) ηa (h) → 0,
as h → 0,
(3.14)
and
w − Ph w b ηa (h) w − Ph w a ,
∀w ∈ V.
(3.15)
Let M (λi ) denote the unit ball in the eigenfunction set corresponding to the eigenvalue λi which is defined by M (λi ) = w ∈ V : w is an eigenfunction of (3.3) corresponding (3.16) to λi and w a = 1 . Then we define δh (λi ) =
sup
inf w − vh a .
w∈M (λi ) vh ∈Vh
(3.17)
In this paper, we are concerned with the finite positive eigenvalues of (3.3). It is well known that [3,4,8] give a thorough theoretical analysis for the eigenvalue problems associated with the compact operators. Based on these results, the following error estimates holds for the eigenpair approximations by the finite element method. Proposition 3.1. ( [3, Lemma 3.7, (3.29b)], [4, P. 699] and [8]) (i) For any eigenfunction approximation u¯i,h of (3.6) (i = 1, 2, · · · , Nh ), there is an eigenfunction ui of (3.3) corresponding to λi such that ui a = 1 and
ui − u¯i,h a ≤ Ci δh (λi ).
(3.18)
ui − u¯i,h b ≤ Ci ηa (h)δh (λi ).
(3.19)
Furthermore,
(ii) For each eigenvalue, we have ¯ j,h − λi | ≤ Ci δ 2 (λi ), |λ h
(3.20)
for j = i, · · · , i + q − 1. Here and hereafter Ci is some constant depending on i but independent of the mesh size h. 7
For generality, we set the multiplicity of our desired eigenvalue is q. It means λi = · · · = λi+q−1 . We use (λi,h , ui,h ), · · · , (λi+q−1,h , ui+q−1,h ) to denote the eigenpair approximations for the eigenvalues λi = · · · = λi+q−1 and their corresponding eigenfunction space M (λi ). Let (3.21) Mh (λi ) = span ui,h , · · · , ui+q−1,h . For two linear spaces A and B, we define
Θ(A, B) =
sup
inf w − v a , Φ(A, B) =
w∈A,wa =1 v∈B
sup
inf w − v b .
w∈A,wa =1 v∈B
We define the gaps between M (λi ) and Mh (λi ) in · a as
h (λi ), M (λi )) , (λi ), Mh (λi )), Θ(M Θ(M (λi ), Mh (λi )) = max Θ(M
(3.22)
and in · b as
h (λi ), M (λi )) . (λi ), Mh (λi )), Φ(M Φ(M (λi ), Mh (λi )) = max Φ(M
4
(3.23)
Multigrid finite element method
In this section, we shall introduce a multigrid finite element (MFE) approach for solving Fredholm integral eigenvalue problem (3.3). In order to do multigrid scheme, we first generate a coarse mesh TH with the mesh size H and the coarse linear finite element space VH is defined on the mesh TH . Then we define a sequence of triangulations Thk of D ⊂ Rd determined as follows. Suppose Th1 (produced from TH by regular refinements) is given and let Thk be obtained from Thk−1 via regular refinement (produce β d subelements) such that hk ≈
1 hk−1 , β
where the integer β denotes the refinement index and larger than 1 (always equals 2 in the text book). Based on this sequence of meshes, we construct the corresponding linear finite element spaces such that V H ⊆ V h1 ⊂ V h2 ⊂ · · · ⊂ V hn and the following relation of approximation accuracy holds 1 γ δhk−1 (λi ), k = 2, · · · , n, ηa (H) δh1 (λi ), δhk (λi ) ≈ β
(4.1)
(4.2)
where the index γ > 1 is determined by the degree of the finite element space Vh and the regularity of the concerned eigenfunctions. 8
4.1
One correction step with iterative method
In this section, we present a type of one correction step to improve the accuracy of the current eigenvalue and eigenfunction approximations. This correction method contains doing the integral iterative step (Sloan iteration step [14,22]) and an eigenvalue problem on a coarse finite element space. Assume we have obtained the eigenpair approximations (λj,hk , uj,hk ) ∈ R × Vhk for i+q−1 j = i, · · · , i + q − 1, where eigenvalues {λj,h }j=i converge to the desired eigenvalue λi of (3.3). Now we introduce a correction step to improve the accuracy of the current eigenpair approximation {λj,hk , uj,hk }i+q−1 j=i . Algorithm 4.1. One Correction Step 1. For j = i, · · · , i + q − 1 Do:. Obtain u j,hk+1 ∈ Vhk+1 by the Sloan iteration u j,hk+1 = Ka,hk+1 uj,hk ,
(4.3)
where Ka,hk+1 is the discretization of Ka on the space Vhk+1 . End Do 2. Define a new finite element space VH,hk+1 = VH +span{ ui,hk+1 , · · · , u i+q−1,hk+1 } and solve the following eigenvalue problem: Find (λj,hk+1 , uj,hk+1 ) ∈ R × VH,hk+1 such that a(uj,hk+1 , uj,hk+1 ) = 1 and a(uj,hk+1 , vH,hk+1 ) = λj,hk+1 b(uj,hk+1 , vH,hk+1 ),
∀vH,hk+1 ∈ VH,hk+1 .
(4.4)
We choose the corresponding eigenpair {λj,hk+1 , uj,hk+1 }i+q−1 as the output of the j=i One Correction Step. In order to simplify the notation and summarize the above two steps, we define {λj,hk+1 , uj,hk+1 }i+q−1 = Correction(VH , {λj,hk , uj,hk }i+q−1 j=i j=i , Vhk+1 ), where VH denotes the background coarse finite element space, {λj,hk , uj,hk }i+q−1 are j=i the given eigenpair approximations, Vhk+1 denotes the computing space. Remark 4.1. The construction of VH,hk+1 is significant for our designation of the multigrid method. The good properties of VH,hk+1 are that it not only keeps the accuracy of { uj,hk+1 }i+q−1 j=i , but also has the duality argument (see the proof of Theorem 4.1). Theorem 4.1. Assume there exists a real number εhk (λi ) such that the given eigenpairs {λj,hk , uj,hk }i+q−1 in Algorithm 4.1 have the following error estimates j=i Θ(M (λi ), Mhk (λi )) εhk (λi ), 9
(4.5)
Φ(M (λi ), Mhk (λi )) ηa (H)εhk (λi ), |λi − λj,hk | ε2hk (λi ),
(4.6) (4.7)
for j = i, · · · , i + q − 1. Then after one correction step, the resultant eigenpair approximation {λhk+1 , uhk+1 }i+q−1 have the following error estimates j=i Θ(M (λi ), Mhk+1 (λi )) εhk+1 (λi ), Φ(M (λi ), Mhk+1 (λi )) ηa (H)εhk+1 (λi ), |λi − λj,hk+1 | ε2hk+1 (λi ),
(4.8) (4.9) (4.10)
where εhk+1 (λi ) := ηa (H)εhk (λi ) + δhk+1 (λi ) and j = i, · · · , i + q − 1. i+q−1 Proof. From (4.5), (4.6) and (4.7), we know there exist an orthogonal basis {uj }j=i of M (λi ) such that
uj − uj,hk a εhk (λi ),
uj − uj,hk b ηa (H)εhk (λi ), |λi − λj,hk | ε2hk (λi ).
(4.11) (4.12) (4.13)
From problems (3.9), (3.3), and (4.3), and inequalities (4.11), (4.12), and (4.13), the following estimates hold for j = i, · · · , i + q − 1
λ−1 j,hk+1 − Phk+1 uj 2a a(λ−1 j,hk+1 − Phk+1 uj , λ−1 j,hk+1 − Phk+1 uj ) i u i u i u −1 −1 j,hk+1 − Phk+1 uj ) = λi b(uj,hk − uj , λi u −1 −1 j,hk+1 − Phk+1 uj a λi uj,hk − uj b λi u −1 ηa (H)εhk (λi ) λi u j,hk+1 − Phk+1 uj a . Then we have
λ−1 j,hk+1 − Phk+1 uj a ηa (H)εhk (λi ), i u
(4.14)
for j = i, · · · , i + q − 1. Combining (4.14) and the error estimate of finite element projection
uj − Phk+1 uj a δhk+1 (λi ), we have
λ−1 j,hk+1 − uj a ηa (H)εhk (λi ) + δhk+1 (λi ), i u
(4.15)
for j = i, · · · , i + q − 1. i+q−1 Now we come to estimate the error for the eigenpair solution {λj,hk+1 , uj,hk+1 }j=i of problem (4.4). Based on the error estimate theory of eigenvalue problem by finite
10
element method (see, e.g., [3, 4] and Proposition 3.1), (4.15), and the definition of the space VH,hk+1 , the following estimates hold Θ(M (λi ), Mhk+1 (λi ))
sup
inf
w∈M (λi ) vH,hk+1 ∈VH,hk+1
w − vH,hk+1 a
max{ λ−1 i,hk+1 − ui a , · · · , λ−1 i+q−1,hk+1 − ui+q−1 a } i u i u (4.16) ηa (H)εhk (λi ) + δhk+1 (λi ), and Φ(M (λi ), Mhk+1 (λi ))
sup
inf
w∈M (λi ) vH,hk+1 ∈VH,hk+1
ηa (H) sup
w − vH,hk+1 −a
inf
w∈M (λi ) vH,hk+1 ∈VH,hk+1
w − vH,hk+1 a
ηa (H)εhk+1 (λi ),
(4.17)
where ηa (H) =
sup
inf
f ∈W,f a =1 v∈VH,hk+1
T f − v a ≤ ηa (H).
(4.18)
From (4.16), (4.17), and (4.18), we can obtain (4.8) and (4.9). The estimate (4.10) can be derived by (3.20) and (4.8).
4.2
Multigrid scheme for the eigenvalue problem
In this section, we introduce a multigrid scheme based on the One Correction Step defined in Algorithm 4.1. This type of multigrid method can obtain the optimal accuracy as same as solving the eigenvalue problem directly in the finest finite element space. Algorithm 4.2. Eigenvalue Multigrid Scheme 1. Construct a series of nested finite element spaces Vh1 , Vh2 , · · · , Vhn such that (4.1) and (4.2) hold. 2. Solve the following eigenvalue problem: Find (λh1 , uh1 ) ∈ R × Vh1 such that b(uh1 , uh1 ) = 1 and λh1 a(uh1 , vh1 ) = b(uh1 , vh1 ),
∀vh1 ∈ Vh1 .
(4.19)
Choose q eigenpairs {λj,h1 , uj,h1 }i+q−1 which approximate to the desired eigenj=i value λi and its eigenspace to do the following correction steps.
11
3. Do k = 1, · · · , n − 1 Obtain new eigenpair approximations {λj,hk+1 , uj,hk+1 }i+q−1 ∈ R × Vhk+1 by j=i Algorithm 4.1 i+q−1 i+q−1 {λj,hk+1 , uj,hk+1 }j=i = Correction(VH , {λj,hk , uj,hk }j=i , Vhk+1 ).
End Do ∈ R × Vhn . Finally, we obtain q eigenpair approximations {λj,hn , uj,hn }i+q−1 j=i Theorem 4.2. After implementing Algorithm 4.2, the resultant eigenpair approxii+q−1 mations {λj,hn , uj,hn }j=i has the following error estimates Θ(M (λi ), Mhn (λi )) δhn (λi ), Φ(M (λi ), Mhn (λi )) ηa (H)δhn (λi ), |λj,hn − λi | δh2n (λi ), j = i, · · · , i + q − 1,
(4.20) (4.21) (4.22)
with the condition Cβ γ ηa (H) < 1, for some constant C. Proof. First we have Θ(M (λi ), Mh1 (λi )) δh1 (λ), Φ(M (λi ), Mh1 (λi )) ηa (h1 )δh1 (λi ) ηa (H)δh1 (λi ), |λj,h1 − λi | δh21 (λi ), j = i, · · · , i + q − 1.
(4.23) (4.24) (4.25)
Let εh1 (λi ) := δh1 (λi ). From (4.23)-(4.25) and Theorem 4.1, we have εhk+1 (λj ) ηa (H)εhk (λj ) + δhk+1 (λj ),
for 1 ≤ k ≤ n − 1.
(4.26)
Then by recursive relation and based on the proof in Theorem 4.1, (4.2), (4.26) and Cβηa (H) < 1, we have εhn (λi ) ηa (H)εhn−1 (λi ) + δhn (λi ) ηa (H)2 εhn−2 (λi ) + ηa (H)δhn−1 (λi ) + δhn (λi ) n n n−k ηa (H) δhk (λi ) = (β γ ηa (H))n−k δhn (λi ) k=1
k=1
β γ ηa (H) δh (λi ) 1 − β γ ηa (H) n δhn (λi ).
(4.27)
This is the estimate (4.20). From the proof of Theorem 4.1 and (4.20), we can obtain the desired results (4.21) and (4.22).
12
4.3
Work estimate of eigenvalue multigrid scheme
In this section, we turn our attention to the estimate of computational work for Algorithm 4.2. We will show that Algorithm 4.2 makes solving the eigenvalue problem need almost the same work as doing one integration process on the finest mesh. First, we define the dimension of each level linear finite element space as Nk := dimVhk . Then we have Nk ≈
1 d(n−k) β
Nn ,
k = 1, 2, · · · , n.
(4.28)
Theorem 4.3. Assume the eigenvalue problem solved in the coarse spaces VH and Vh1 need work O(MH ) and O(Mh1 ), respectively, and the work of Sloan iteration in each level space Vhk is O(Nkα ) (α ≥ 1) for j = i, · · · , i + q − 1 and k = 2, 3, · · · , n. We do Sloan iteration in Step 1 of Algorithm 4.1 in parallel and solve the eigenvalue problem in Step 2 of Algorithm 4.1 in one processor. Then the most work involved in each processor of Algorithm 4.2 is O(Nnα + MH log(Nk ) + Mh1 ). Furthermore, the complexity will be O(Nnα ) provided MH Nnα and Mh1 ≤ Nnα . Proof. Let Wk denote the work in each processor in the k-th finite element space Vhk for the k-th One Correction Step. Then with the definition of Algorithm 4.1, we have Wk = O(Nkα + MH ).
(4.29)
Iterating (4.29) and using the fact (4.28), we obtain Total Work =
n
n α Wk = O Mh1 + Nk + MH
k=1
k=2
= O Mh1 + (n − 1)MH + = O Mh1 + (n − 1)MH + =
Nkα
k=2
O(Nnα
n
n 1 αd(n−k) k=2
+ MH log Nn + Mh1 ).
β
Nnα
(4.30)
This is the desired result O(Nnα + MH log Nn + Mh1 ) and the one O(Nnα ) can be obtained by the conditions MH Nnα and Mh1 ≤ Nnα . Remark 4.2. From Theorem 4.3, we can find Algorithm 4.2 can make the complexity of the Fredhold integral eigenvalue problems become the complexity of the integral process on the finest mesh. If we have an efficient integral scheme with α = 1, we can make the computational complexity of Algorithm 4.2 become optimal.
13
5
Numerical examples
In this section, we shall provide with some numerical examples to illustrate the efficiency of our MFE method, for computing the Karhunen-L`oeve expansion of random fields.
5.1
The Gaussian covariance function
Our first example is the following FIEP with Gaussian covariance function, (x−y)2 Va (x, y)v(y)dy = λv(x), with Va (x, y) = e− L , x, y ∈ D.
(5.1)
D
It is well know that the Gaussian covariance function is differentiable everywhere, and the eigenvalue sequence admits a fast decay rate. In this example, the linear finite element is adopted. We fist set the domain D to be (−1, 1) × (−1, 1). The corresponding numerical results are shown in Figure 1 (left). A salient advantage of using finite element methods is that it is easy to handle general domain. We next consider a L-shape domain D = (−1, 1) × (−1, 1)\[0, 1] × [0, 1] for (5.1). The numerical results are also shown in Figure 1 (right). From Figure 1, we can find that Algorithm 4.2 can reach the fourth order convergence rate which validate our theoretical results. Eigenvalues errors on unit square domain
−3
Eigenvalues errors on L−shape domain
−2
10
10 Σ10 |λ −λ | j=1 j j,h
Σ15 |λ −λ | j=1 j
slope=−2
−4
10
10
−5
−4
10 Errors
Errors
10
−6
10
−7
−5
10
−6
10
10
−8
−7
10
10
−9
10
j,h
slope=−2
−3
−8
1
10
2
10
3
10
4
10
Number of elements
10
1
10
2
10
3
10
4
10
Number of elements
Figure 1: Errors for the first 10 eigenvalues by the multi-grid way: the left subfigure is on the unit square domain and the right subfigure is on the L-shape domain, number of elements ≈ 1/h2
14
5.2
The exponential covariance function
We next consider the FIEP with an exponential covariance function, i.e., Va (x, y)v(y)dy = λv(x), with Va (x, y) = e−c|x−y| , x, y ∈ D.
(5.2)
D
The above exponential kernel is often used on benchmark tests, since closed formulas are available for its eigenvalues and eigenfunctions. For example, let D ≡ (−1, 1), a random field a(x, ω) with exponential covariance function given by (5.2) admits the following Karhunen-L`oeve expansion [10]: a(x, ω) =
∞ ξn λn fn (x) + ξn∗ λ∗n fn∗ (x) ,
(5.3)
n=1
where 2c 2c , λ∗n = ∗ 2 , 2 +c (μn ) + c2 cos(μn x) cos(μ∗n x) , fn∗ (x) = . fn (x) = sin(2μ∗n ) n) 1 + 1 + sin(2μ ∗ 2μn 2μn λn =
μ2n
(5.4) (5.5)
The parameters μn , μ∗n are solutions of the following problems μ∗ + ctan(μ∗ ) = 0.
c − μtan(μ) = 0,
(5.6)
In this example, we choose c = 1. The piecewise constant finite element method is adopted to compute the eigenvalue problem (5.2). The corresponding numerical results are shown in Figure 2 which also confirm the theoretical second order convergence.
5.3
A non-stationary case: the Wiener-Levy process
Finally, we consider the covariance function of the Wiener-Levy process which reads Va (x, y) = min{x, y},
x, y ∈ (0, 1) × (0, 1).
(5.7)
This example is provided to demonstrate that our method can be used as a simulation tool for non-stationary processes. The resulting normalized eigenfunctions and eigenvalues of the corresponding FIEP are [10]: x √ 4 √ , f (x) = 2sin . (5.8) λn = 2 n π (2n + 1)2 λn In this example, the eigenvalue problem (5.7) is discretized by the piecewise constant finite element method. We investigate the numerical results for the coarse 15
Eigenpair errors for H=1/20
−3
Eigenpair errors for H=1/40
−4
10
10 Σ11 |λ −λ | j=1 j j,h Σ11 |λ −λdir| j=1 j j,h
Σ11 |λ −λ | j=1 j
slope=−2
−4
10
j,h
Σ11 |λ −λdir| j=1 j j,h slope=−2 −5
Errors
Errors
10 −5
10
−6
10 −6
10
−7
10
−7
1
2
10
10
10
3
10
1
2
10
10
Number of elements
3
10
Number of elements
Figure 2: Errors for the first 11 eigenvalues by the multi-grid way with H = 1/20 (left) or H = 1/40 (right): λj,h is obtained by the multilevel method and λdir j,h denotes the approximation by the direct eigenvalue solving, number of elements ≈ 1/h mesh H = 1/20 and H = 1/40 for the biggest 11 eigenvalues and the associated first 5 eigenfunctions. The results are shown in Figure 3 (left for H = 1/20 and right for H = 1/40). It can be seen that the convergence orders of eigenvalue and eigenfunction approximations correspond to the theoretical prediction, respectively. Eigenpair errors for H=1/20
0
10
−1
−1
10
10
−2
−2
10 Errors
10 Errors
Eigenpair errors for H=1/40
0
10
−3
10
−4
−3
10
−4
10
10
Σ11 |λ −λ |
Σ11 |λ −λ | j=1 j j,h −5
10
1
10
j,h
Σ5j=1||uj−uj,h||0
−5
10
slope=−2 slope=−1
slope=−2 slope=−1
−6
10
j=1 j
Σ5 ||u −u || j=1 j j,h 0
−6
2
10
3
10
Number of elements
10
1
10
2
10
3
10
Number of elements
Figure 3: Errors for the first 11 eigenvalues and the first 5 eigenfunctions by the multi-grid way with H = 1/20 (left) or H = 1/40 (right), number of elements ≈ 1/h
16
6
Concluding remarks
In this paper, we propose a type of multilevel finite element method for the Fredholm integral eigenvalue problems based on the multilevel correction scheme. In this multilevel finite element method, we only need to solve integral eigenvalue problems in the coarsest finite element space and the Sloan iteration for the integral equation. The multilevel correction idea can also be extended to the collocation method for the integral eigenvalue problems which may be our future work.
Acknowledgment This first author is supported in part by the National Science Foundation of China (NSFC 91330202, 11371026, 11001259, 11201501, 2011CB309703), the National Center for Mathematics and Interdisciplinary Science (CAS), and the President Foundation of AMSS-CAS. The second author is supported by the National Science Foundation of China (No.91130003 and No.11201461).
References [1] K. Atkinson, The Numerical Solution of Integral Equations of the Second Kind, Cambridge University Press, Cambridge, 1997. [2] I. Babuˇska, F. Nobile and R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM J. Numer. Anal., 45 (2007), 1005-1034. [3] I. Babuˇska and J. E. Osborn, Finite element-Galerkin approximation of the eigenvalues and eigenvectors of selfadjoint problems, Math. Comp. 52 (1989), 275-297. [4] I. Babuˇska and J. E. Osborn, Eigenvalue Problems, In Handbook of Numerical Analysis, Vol. II, (Eds. P. G. Lions and Ciarlet P.G.), Finite Element Methods (Part 1), North-Holland, Amsterdam, 641-787, 1991. [5] J. B¨ack, F. Nobile, L. Tamellini, and R. Tempone, On the optimal polynomial approximation of stochastic PDEs by Galerkin and Collocation methods, Math. Mod. Methods Appl. Sci., 22(9) (2012), 1-33. [6] H. Brunner, A. Iserles, S. P. Nørsett, The spectral problem for a class of highly oscillatory Fredholm integral operators, IMA J. Numer. Anal., 30(1) (2010), 108-130.
17
[7] H. Brunner, Collocation Methods for Volterra Integral and Related Functional Equations, Cambridge University Press, Cambridge, 2004. [8] F. Chatelin, Spectral Approximation of Linear Operators, Academic Press Inc, New York, 1983. [9] F. R. DiNapoli, F. H. Middleton, Surface-wave propagation in a continuously stratified medium, J. Acoust. Soc. Am., 39(5A) (1966), 899-903. [10] R. Ghanem, P. Spanos, Stochastic Finite Element: A Spectral Aproach, Springer, New York, 1991. [11] M. Girolami, Orthogonal series density estimation and the kernel eigenvalue problem, Neural Comput., 14(3) (2002), 669-688. [12] H. B. Keller, On the accuracy of finite difference approximations to the eigenvalues of differential and integral operators, Numer. Math., 7(5) (1965), 412-419. [13] K. Khare, Sampling theorem, bandlimited integral kernels and inverse problems, Inverse Problems, 23(4) (2007), 1395-1416. [14] Q. Lin, Some problems concerning approximate solutions of operator equations, Acta Math. Sinica, 22 (1979), 219-230 (Chinese). [15] Q. Lin, H. Xie, A multi-level correction scheme for eigenvalue problems, Math. Comput., http://dx.doi.org/10.1090/S0025-5718-2014-02825-1. [16] F. Nobile, R. Tempone and C. Webster, A sparse grid stochastic collocation method for partial differential equations with random input data, SIAM J. Numer. Anal., 46(5) (2008), 2309-2345 [17] S. P. Oliveira and J. S. Azevedob, Spectral element approximation of Fredholm integral eigenvalue problems, J. Comput. Appl. Math., 257 (2014), 46-56. [18] J. Parvizian, A. D¨ uster and E. Rank, Finite cell method, Computational Mechanics, 41 (2007), 121-133. [19] K. K. Phoon, S.P. Huang, S. T. Quek, Implementation of Karhunen–Lo`eve expansion for simulation using a wavelet-Galerkin scheme, Probab. Eng. Mech., 17(3) (2002), 293-303. [20] F. Riesz and B. Sz-Nagy, Functional Analysis, Dover, 1990. [21] C. Schwab and R. A. Todor, Karhunen-Lo`eve approximation of random fields by generalized fast multipole methods, J. Comput. Phys., 217 (2006), 100-122. [22] I. Sloan, Iterated Galerkin method for eigenvalue problems, SIAM J. Numer. Anal., 13(5) (1976), 753-760. 18
[23] T. Tang and T. Zhou, Discrete least square projection in unbounded domain with random evaluations and its application to parametric uncertainty quantification, SIAM J. Sci. Comput., 36(5) (2014), A2272-A2295. [24] T. Tang and T. Zhou, Convergence analysis for stochastic collocation methods to scalar hyperbolic equations, Commun. Comput. Phys., 8 (2010), 226-248. [25] R. Todor, Robust eigenvalue computation for smoothing operators, SIAM J. Numer. Anal., 44(2) (2006), 865-878. [26] N. Wiener, The homogeneous chaos, Am. J. Math., 60 (1938), 897-936. [27] H. Xie, A type of multilevel method for the Steklov eigenvalue problem, IMA J. Numer. Anal., 34(2) (2014), 592-608. [28] H. Xie, A multigrid method for eigenvalue problem, J. Comput. Phys., 274 (2014), 550-561. [29] D. Xiu and J. S. Hesthaven, High-order collocation methods for differential equations with random inputs, SIAM J. Sci. Comput, 27 (2005), 1118-1139. [30] D. Xiu and G. E. Karniadakis, The Wiener-Askey polynomial chaos for stochastic differential equations, SIAM J. Sci. Comput., 24(2) (2002), 619-644. [31] J. Xu and A. Zhou, A two-grid discretization scheme for eigenvalue problems, Math. Comput., 70(233) (2001), 17-25. [32] T. Zhou and T. Tang, Galerkin methods for stochastic hyperbolic problems using Bi-orthogonal polynomials, J. Sci. Comput., 51 (2012), 274-292. [33] T. Zhou, A. Narayan and Z. Xu, Multivariate discrete least-squares approximations with a new type of collocation grid, SIAM J. Sci. Comput., 36(5) (2014), A2401-A2422.
19
Hehu Xie, Tao Zhou
PII: DOI: Reference:
S0021-9991(15)00642-7 http://dx.doi.org/10.1016/j.jcp.2015.09.043 YJCPH 6146
To appear in:
Journal of Computational Physics
Received date: Revised date: Accepted date:
9 November 2014 16 July 2015 24 September 2015
Please cite this article in press as: H. Xie, T. Zhou, A multilevel finite element method for Fredholm integral eigenvalue problems, J. Comput. Phys. (2015), http://dx.doi.org/10.1016/j.jcp.2015.09.043
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A multilevel finite element method for Fredholm integral eigenvalue problems Hehu Xie∗ and Tao Zhou†
Abstract In this work, we proposed a multigrid finite element (MFE) method for solving the Fredholm integral eigenvalue problems. The main motivation for such studies is to compute the Karhunen-Lo`eve expansions of random fields, which plays an important role in the applications of uncertainty quantification. In our MFE framework, solving the eigenvalue problem is converted to doing a series of integral iterations and eigenvalue solving in the coarsest mesh. Then, any exist efficient integration scheme can be used for the associated integration process. The error estimates are provided, and the computational complexity is analyzed. It is noticed that the total computational work of our method is comparable with a single integration step in the finest mesh. Several numerical experiments are presented to validate the efficiency of the proposed numerical method. Keywords. Uncertainty quantification, Karhunen-L`oeve expansion, Fredholm eigenvalue problem, multigrid finite element. AMS subject classifications. 65N30, 65N25, 65L15, 65B99.
1
Introduction
Eigenvalue problems involving the Fredholm integral equation arises in many areas of sciences and engineering (see e.g. [6, 9, 11, 13] and references therein). Here, we focus on its particular application for Karhunen-L`oeve approximation of random fields, which plays an important role in uncertainty quantification (UQ) [2, 10]. In UQ, a typical problem is to numerically solve a PDE with random field coefficients. Then, to transform the original problem into a PDE with finite random parameters, one usually adopts the truncated Karhunen-L`oeve expansion of the random field ∗
LSEC, NCMIS, Institute of Computational Mathematics, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China. Email: [email protected]. † LSEC, Institute of Computational Mathematics, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China. Email: [email protected].
1
[16,32]. In such framework, the random field is represented in terms of uncorrelated random variables and eigenfunctions of a homogeneous Fredholm integral equation of the second kind, whose kernel is the covariance function of the random field. Thus, to efficiently solve the original problem, a key issue is how to efficiently build the truncated Karhunen-L`oeve expansion, by solving a Fredholm integral eigenvalue problem. Despite the large amount of the literature on numerical methods for Fredholm integral equations (see e.g. [1, 7] and references therein), there are only a few numerical results for Fredholm eigenvalue problems (see e.g. [12, 22, 25]). From the UQ application point of view (especially for the computations of Karhunen-L`oeve approximations), one can refer to [21] for a finite element method combined with the fast multipole methods, to [19] for a wavelet Galerkin arpproach using Haar basis functions, to [18] for the finite cell method, to [17] for a spectral element approximation, to name a few. In this paper, we propose a multigrid finite element method for solving the Fredholm integral eigenvalue problem with the idea coming from [15, 27, 28]. In our framework, solving the eigenvalue problem is converted to a series of integral processes and eigenvalue solving on the coarsest mesh. Then, any exist efficient numerical integral scheme can be used for the integral processes. The corresponding error estimates are given, and the computational complexity is analyzed. It is noticed that the total computational work of our method is comparable with the integral scheme. Several numerical experiments are presented to validate the efficiency of the method. The rest of the paper is organized as follows. In section 2, we give a brief review of the Karhunen-L`oeve expansion of random fields, and introduce its applications in UQ. In section 3, a multigrid finite element method is proposed to solve the corresponding Fredholm integral eigenvalue problem, moreover, error estimates and complexity analysis will be provided. Several numerical examples are given in Section 4 to illustrate the efficiency of the proposed method. We finally give some conclusions in section 5.
2
Karhunen-L` oeve expansion of random fields
In this section, we shall give a brief review of the Karhunen-L`oeve expansion of random fields, and introduce its applications in UQ. To this end, let (Ω, Σ, P) be a probability space, where Ω is the space of events, Σ the corresponding σ-algebra and P the probability measure. A spatial dependent random field is a P -measurable map a(x, ω) : Ω → L∞ (D) (where D ⊂ Rd is the physical domain). The mean of the random field a(x, ω) is defined by a ¯(x) = E[a(x, ·)] = a(x, ω)dP(ω). (2.1) Ω
2
We assume that all the random fields considered in this paper are square integrable (i.e., the so-called second order random fields), that is, E[a2 (x, ·)]dx < +∞. (2.2) D
We define the covariance function of a space-dependent random field a(x, ω) as ¯(x)) (a(y, ·) − a ¯(y)) , x, y ∈ D. (2.3) Va (x, y) = E (a(x, ·) − a Given such a random field a(x, ω), prediction of for example a concentration φ(x, ω) may requires solution of a stochastic partial differential equation such as (here we take the elliptic problem as an example) −∇ · (a(x, ω)∇x φ(x, ω)) = f (x),
∀x ∈ D.
(2.4)
It is well known that any second order random field a(x, ω), with continuous and positive definite covariance function Va : D × D → R, can be represented as an infinite sum of random variables, by means of the Karhunen-L`oeve expansion [10]. To this end, we introduce the compact and self-adjoint operator Ka : L2 (D) → L2 (D), which is defined by Va (x, ·)u(x)dx, ∀v ∈ L2 (D). (2.5) Ka u(·) = D
Then, consider the sequence of non-negative decreasing eigenvalues of Ka , {λi }∞ i=1 , ∞ and the corresponding sequence of orthonormal eigenfunctions, {vi }i=1 , satisfying ui uj dx = δij for i, j ∈ N + , (2.6) Ka ui = λi ui , D
where δij is the Kronecker symbol. In addition, define the mutually uncorrelated real random variables 1 Yi (ω) := √ (2.7) a(x, ω) − a ¯(x) ui (x)dx, i = 1, 2, ... λi D with zero mean and unit variance, i.e. E[Yi ] = 0 and E[Yi Yj ] = δij for i, j ∈ N + . Then, the N -term truncated Karhunen-L`oeve expansion of the random field a(x, ω), which we denote by aN , is defined by ¯(x) + aN (x, ω) = a
N λi Yi (ω)ui (x), ∀N ∈ N + .
(2.8)
i=1
By Mercer’s theorem [20], it follows that lim sup E
N →∞ x∈D
a(x, ·) − aN (x, ·)
2
= lim sup N →∞ x∈D
3
∞ i=N +1
λi u2i (x) = 0.
(2.9)
Note that such N -term Karhunen-L`oeve truncation is optimal in the mean square sense [10]. By using the truncated Karhunen-L`oeve expansion, one converts the original stochastic elliptic problem (2.4) into the following elliptic problem with N random parameters: −∇ · (a(x, Y)∇x φ(x, Y)) = f (x),
∀x ∈ D,
(2.10)
with Y being the random vector Y = (Y1 , ...YN ). For the above finite dimensional stochastic problem, many approaches have been developed to get the quantity of interests of the solution φ(x, Y)). For example, the stochastic Galerkin method which proposed originally in [10, 26] and further developed in [30], the stochastic collocation method [16, 24, 29], the polynomial regression [23, 33], etc. By the discussions above, a key issue is to efficiently get the truncated KarhunenL`oeve expansion, i.e., efficiently solve the corresponding Fredholm integral eigenvalue problem: Ka u(x) = Va (x, y)u(y)dy = λu(x), ∀x ∈ D. (2.11) D
The purpose of this work is to propose a multigrid finite element method for solving (2.11). To this end, we remark that regularity results of the eigenfunctions and decay properties of the eigenvalues have been investigated in [21]. Particularly, it is proved in [21] that the eigenfunctions inherited the regularity of the corresponding covariance function, and furthermore, the decay properties of the eigenvalues is also depend on the regularity of the covariance function. More precisely, we have the following propositions [21]: Proposition 2.1. Let Va ∈ L2 (D × D) be a symmetric kernel defining the compact and nonnegative integral operator 2 2 Ka : L (D) → L (D), Ka u(x) = Va (x, y)u(y)dy. (2.12) D
If Va is piecewise analytic (see a detailed definition in [21]) on D × D and {λm }m≥1 is the eigenvalue sequence of its associated operator Ka , then there exist constants c1 , c2 > 0 depending only on Va such that 0 ≤ λm ≤ c1 e−c2 m
1/d
,
∀m ≥ 1.
If Va is piecewise H k,0 := H k ⊗ L2 on D × D with k > 0, then there exists a constant c3 > 0 depending only on Va such that λm ≤ c3 m−k/d ,
∀m ≥ 1.
Proposition 2.2. Assume that the covariance kernel Va (x, y) of a random field a(x, ω) is piecewise analytic (piecewise H p,q ) on D × D (for rigorous definitions, see 4
[21]). Then the Karhunen-L`oeve eigenfunctions {vm (x)}m≥1 of a(x, ω) are piecewise analytic (piecewise H p ). Moreover, for D ⊂ Rd a bounded domain and Va piecewise smooth on D ×D, such that the each piece Dj have the uniform cone property, we denote by {λm , vm }m≥1 the eigenpair sequence of the associated integral operator Ka via (2.12), such that
vm L2 = 1 for m ≥ 1. Then for any s > 0 and any multi-index α ∈ N d there exists cV,s,α > 0 such that
∂ α vm L∞ (Dj ) ≤ cV,s,α λ−s m ,
∀m ≥ 1.
For more properties of the Karhunen-L`oeve expansion of random fields, please refer to [21]. We also remark that the Karhunen-L`oeve expansion is not the only way to approximate random fields, other approaches can also be considered, see e.g. [5] for a Fourier-based decomposition, where the authors uses trigonometric polynomials as basis functions in the physical space.
3
Fredholm integral eigenvalue problem and its finite element discretization
In this section, we shall introduce a multigrid finite element (MFE) approach for solving the following Fredholm integral eigenvalue problems (FIEPs): Ka u = λu,
(3.1)
where the integral operator Ka is defined as Ka u := Va (x, y)u(y)dy,
(3.2)
D
with Va (x, y) being some covariance functions of certain random fields. We denote by V the space L2 (D) in this paper. The weak form of (3.1) is: Find (λ, u) ∈ R × V such that a(u, u) = 1, where
and λa(u, v) = b(u, v),
(3.3)
uvdx,
a(u, v) =
v ∈ V,
Ka uvdx.
b(u, v) =
D
D
We define the corresponding norm · a and semi-norm · b as:
v 2a := a(v, v),
v 2b := b(v, v). 5
(3.4)
For the eigenvalue λ, there exists the following Rayleigh quotient expression (see, e.g., [3, 4, 31]) λ=
b(u, u) . a(u, u)
(3.5)
We shall use the finite element discretizition for problem (3.3). To do this, let us define the finite element approximations of the problem (3.3). First we generate a shape-regular decomposition of the computing domain D ⊂ Rd (d = 2, 3) into triangles or rectangles for d = 2 (tetrahedrons or hexahedrons for d = 3). The diameter of a cell K ∈ Th is denoted by hK . The mesh diameter h describes the maximum diameter of all cells K ∈ Th . Based on the mesh Th , we construct a finite element space Vh ⊂ V as a piecewise polynomial space. Then we define the approximation for the eigenpair (λ, u) of (3.3) by the finite element method as: ¯ h , u¯h ) ∈ R × Vh such that Find (λ a(¯ uh , u¯h ) = 1,
¯ h a(¯ and λ uh , vh ) = b(¯ uh , vh ),
∀vh ∈ Vh .
(3.6)
For (3.6), the following Rayleigh quotient expression for λh holds [3, 4, 31]: uh , u¯h ) ¯ h = b(¯ λ . a(¯ uh , u¯h )
(3.7)
We assume that Vh ⊂ V is a family of finite-dimensional spaces that satisfy the following assumption: lim inf w − v a = 0,
h→0 v∈Vh
∀w ∈ V.
(3.8)
Let Ph be the finite element projection operator of V onto Vh that defined by a(w − Ph w, vh ) = 0,
∀w ∈ V and ∀vh ∈ Vh .
(3.9)
Obviously, we have
Ph w a ≤ w a ,
∀w ∈ V.
(3.10)
as h → 0.
(3.11)
For any w ∈ V , by (3.8) we have
w − Ph w a → 0, Define ηa (h) as ηa (h) =
sup
inf T f − v a ,
f ∈V,f a =1 v∈Vh
6
(3.12)
where the operator T : V → V is defined as ∀f ∈ V and ∀v ∈ V.
a(T f, v) = b(f, v),
(3.13)
As we know, the operator T : V → V is compact. Then the results about the eigenvalues of compact operators can be used here. In order to derive the error estimate of eigenpair approximation in the weak semi-norm · b , we need the following error estimates of the finite element projection operator Ph . Lemma 3.1. ( [3, Lemma 3.3 and Lemma 3.4]) ηa (h) → 0,
as h → 0,
(3.14)
and
w − Ph w b ηa (h) w − Ph w a ,
∀w ∈ V.
(3.15)
Let M (λi ) denote the unit ball in the eigenfunction set corresponding to the eigenvalue λi which is defined by M (λi ) = w ∈ V : w is an eigenfunction of (3.3) corresponding (3.16) to λi and w a = 1 . Then we define δh (λi ) =
sup
inf w − vh a .
w∈M (λi ) vh ∈Vh
(3.17)
In this paper, we are concerned with the finite positive eigenvalues of (3.3). It is well known that [3,4,8] give a thorough theoretical analysis for the eigenvalue problems associated with the compact operators. Based on these results, the following error estimates holds for the eigenpair approximations by the finite element method. Proposition 3.1. ( [3, Lemma 3.7, (3.29b)], [4, P. 699] and [8]) (i) For any eigenfunction approximation u¯i,h of (3.6) (i = 1, 2, · · · , Nh ), there is an eigenfunction ui of (3.3) corresponding to λi such that ui a = 1 and
ui − u¯i,h a ≤ Ci δh (λi ).
(3.18)
ui − u¯i,h b ≤ Ci ηa (h)δh (λi ).
(3.19)
Furthermore,
(ii) For each eigenvalue, we have ¯ j,h − λi | ≤ Ci δ 2 (λi ), |λ h
(3.20)
for j = i, · · · , i + q − 1. Here and hereafter Ci is some constant depending on i but independent of the mesh size h. 7
For generality, we set the multiplicity of our desired eigenvalue is q. It means λi = · · · = λi+q−1 . We use (λi,h , ui,h ), · · · , (λi+q−1,h , ui+q−1,h ) to denote the eigenpair approximations for the eigenvalues λi = · · · = λi+q−1 and their corresponding eigenfunction space M (λi ). Let (3.21) Mh (λi ) = span ui,h , · · · , ui+q−1,h . For two linear spaces A and B, we define
Θ(A, B) =
sup
inf w − v a , Φ(A, B) =
w∈A,wa =1 v∈B
sup
inf w − v b .
w∈A,wa =1 v∈B
We define the gaps between M (λi ) and Mh (λi ) in · a as
h (λi ), M (λi )) , (λi ), Mh (λi )), Θ(M Θ(M (λi ), Mh (λi )) = max Θ(M
(3.22)
and in · b as
h (λi ), M (λi )) . (λi ), Mh (λi )), Φ(M Φ(M (λi ), Mh (λi )) = max Φ(M
4
(3.23)
Multigrid finite element method
In this section, we shall introduce a multigrid finite element (MFE) approach for solving Fredholm integral eigenvalue problem (3.3). In order to do multigrid scheme, we first generate a coarse mesh TH with the mesh size H and the coarse linear finite element space VH is defined on the mesh TH . Then we define a sequence of triangulations Thk of D ⊂ Rd determined as follows. Suppose Th1 (produced from TH by regular refinements) is given and let Thk be obtained from Thk−1 via regular refinement (produce β d subelements) such that hk ≈
1 hk−1 , β
where the integer β denotes the refinement index and larger than 1 (always equals 2 in the text book). Based on this sequence of meshes, we construct the corresponding linear finite element spaces such that V H ⊆ V h1 ⊂ V h2 ⊂ · · · ⊂ V hn and the following relation of approximation accuracy holds 1 γ δhk−1 (λi ), k = 2, · · · , n, ηa (H) δh1 (λi ), δhk (λi ) ≈ β
(4.1)
(4.2)
where the index γ > 1 is determined by the degree of the finite element space Vh and the regularity of the concerned eigenfunctions. 8
4.1
One correction step with iterative method
In this section, we present a type of one correction step to improve the accuracy of the current eigenvalue and eigenfunction approximations. This correction method contains doing the integral iterative step (Sloan iteration step [14,22]) and an eigenvalue problem on a coarse finite element space. Assume we have obtained the eigenpair approximations (λj,hk , uj,hk ) ∈ R × Vhk for i+q−1 j = i, · · · , i + q − 1, where eigenvalues {λj,h }j=i converge to the desired eigenvalue λi of (3.3). Now we introduce a correction step to improve the accuracy of the current eigenpair approximation {λj,hk , uj,hk }i+q−1 j=i . Algorithm 4.1. One Correction Step 1. For j = i, · · · , i + q − 1 Do:. Obtain u j,hk+1 ∈ Vhk+1 by the Sloan iteration u j,hk+1 = Ka,hk+1 uj,hk ,
(4.3)
where Ka,hk+1 is the discretization of Ka on the space Vhk+1 . End Do 2. Define a new finite element space VH,hk+1 = VH +span{ ui,hk+1 , · · · , u i+q−1,hk+1 } and solve the following eigenvalue problem: Find (λj,hk+1 , uj,hk+1 ) ∈ R × VH,hk+1 such that a(uj,hk+1 , uj,hk+1 ) = 1 and a(uj,hk+1 , vH,hk+1 ) = λj,hk+1 b(uj,hk+1 , vH,hk+1 ),
∀vH,hk+1 ∈ VH,hk+1 .
(4.4)
We choose the corresponding eigenpair {λj,hk+1 , uj,hk+1 }i+q−1 as the output of the j=i One Correction Step. In order to simplify the notation and summarize the above two steps, we define {λj,hk+1 , uj,hk+1 }i+q−1 = Correction(VH , {λj,hk , uj,hk }i+q−1 j=i j=i , Vhk+1 ), where VH denotes the background coarse finite element space, {λj,hk , uj,hk }i+q−1 are j=i the given eigenpair approximations, Vhk+1 denotes the computing space. Remark 4.1. The construction of VH,hk+1 is significant for our designation of the multigrid method. The good properties of VH,hk+1 are that it not only keeps the accuracy of { uj,hk+1 }i+q−1 j=i , but also has the duality argument (see the proof of Theorem 4.1). Theorem 4.1. Assume there exists a real number εhk (λi ) such that the given eigenpairs {λj,hk , uj,hk }i+q−1 in Algorithm 4.1 have the following error estimates j=i Θ(M (λi ), Mhk (λi )) εhk (λi ), 9
(4.5)
Φ(M (λi ), Mhk (λi )) ηa (H)εhk (λi ), |λi − λj,hk | ε2hk (λi ),
(4.6) (4.7)
for j = i, · · · , i + q − 1. Then after one correction step, the resultant eigenpair approximation {λhk+1 , uhk+1 }i+q−1 have the following error estimates j=i Θ(M (λi ), Mhk+1 (λi )) εhk+1 (λi ), Φ(M (λi ), Mhk+1 (λi )) ηa (H)εhk+1 (λi ), |λi − λj,hk+1 | ε2hk+1 (λi ),
(4.8) (4.9) (4.10)
where εhk+1 (λi ) := ηa (H)εhk (λi ) + δhk+1 (λi ) and j = i, · · · , i + q − 1. i+q−1 Proof. From (4.5), (4.6) and (4.7), we know there exist an orthogonal basis {uj }j=i of M (λi ) such that
uj − uj,hk a εhk (λi ),
uj − uj,hk b ηa (H)εhk (λi ), |λi − λj,hk | ε2hk (λi ).
(4.11) (4.12) (4.13)
From problems (3.9), (3.3), and (4.3), and inequalities (4.11), (4.12), and (4.13), the following estimates hold for j = i, · · · , i + q − 1
λ−1 j,hk+1 − Phk+1 uj 2a a(λ−1 j,hk+1 − Phk+1 uj , λ−1 j,hk+1 − Phk+1 uj ) i u i u i u −1 −1 j,hk+1 − Phk+1 uj ) = λi b(uj,hk − uj , λi u −1 −1 j,hk+1 − Phk+1 uj a λi uj,hk − uj b λi u −1 ηa (H)εhk (λi ) λi u j,hk+1 − Phk+1 uj a . Then we have
λ−1 j,hk+1 − Phk+1 uj a ηa (H)εhk (λi ), i u
(4.14)
for j = i, · · · , i + q − 1. Combining (4.14) and the error estimate of finite element projection
uj − Phk+1 uj a δhk+1 (λi ), we have
λ−1 j,hk+1 − uj a ηa (H)εhk (λi ) + δhk+1 (λi ), i u
(4.15)
for j = i, · · · , i + q − 1. i+q−1 Now we come to estimate the error for the eigenpair solution {λj,hk+1 , uj,hk+1 }j=i of problem (4.4). Based on the error estimate theory of eigenvalue problem by finite
10
element method (see, e.g., [3, 4] and Proposition 3.1), (4.15), and the definition of the space VH,hk+1 , the following estimates hold Θ(M (λi ), Mhk+1 (λi ))
sup
inf
w∈M (λi ) vH,hk+1 ∈VH,hk+1
w − vH,hk+1 a
max{ λ−1 i,hk+1 − ui a , · · · , λ−1 i+q−1,hk+1 − ui+q−1 a } i u i u (4.16) ηa (H)εhk (λi ) + δhk+1 (λi ), and Φ(M (λi ), Mhk+1 (λi ))
sup
inf
w∈M (λi ) vH,hk+1 ∈VH,hk+1
ηa (H) sup
w − vH,hk+1 −a
inf
w∈M (λi ) vH,hk+1 ∈VH,hk+1
w − vH,hk+1 a
ηa (H)εhk+1 (λi ),
(4.17)
where ηa (H) =
sup
inf
f ∈W,f a =1 v∈VH,hk+1
T f − v a ≤ ηa (H).
(4.18)
From (4.16), (4.17), and (4.18), we can obtain (4.8) and (4.9). The estimate (4.10) can be derived by (3.20) and (4.8).
4.2
Multigrid scheme for the eigenvalue problem
In this section, we introduce a multigrid scheme based on the One Correction Step defined in Algorithm 4.1. This type of multigrid method can obtain the optimal accuracy as same as solving the eigenvalue problem directly in the finest finite element space. Algorithm 4.2. Eigenvalue Multigrid Scheme 1. Construct a series of nested finite element spaces Vh1 , Vh2 , · · · , Vhn such that (4.1) and (4.2) hold. 2. Solve the following eigenvalue problem: Find (λh1 , uh1 ) ∈ R × Vh1 such that b(uh1 , uh1 ) = 1 and λh1 a(uh1 , vh1 ) = b(uh1 , vh1 ),
∀vh1 ∈ Vh1 .
(4.19)
Choose q eigenpairs {λj,h1 , uj,h1 }i+q−1 which approximate to the desired eigenj=i value λi and its eigenspace to do the following correction steps.
11
3. Do k = 1, · · · , n − 1 Obtain new eigenpair approximations {λj,hk+1 , uj,hk+1 }i+q−1 ∈ R × Vhk+1 by j=i Algorithm 4.1 i+q−1 i+q−1 {λj,hk+1 , uj,hk+1 }j=i = Correction(VH , {λj,hk , uj,hk }j=i , Vhk+1 ).
End Do ∈ R × Vhn . Finally, we obtain q eigenpair approximations {λj,hn , uj,hn }i+q−1 j=i Theorem 4.2. After implementing Algorithm 4.2, the resultant eigenpair approxii+q−1 mations {λj,hn , uj,hn }j=i has the following error estimates Θ(M (λi ), Mhn (λi )) δhn (λi ), Φ(M (λi ), Mhn (λi )) ηa (H)δhn (λi ), |λj,hn − λi | δh2n (λi ), j = i, · · · , i + q − 1,
(4.20) (4.21) (4.22)
with the condition Cβ γ ηa (H) < 1, for some constant C. Proof. First we have Θ(M (λi ), Mh1 (λi )) δh1 (λ), Φ(M (λi ), Mh1 (λi )) ηa (h1 )δh1 (λi ) ηa (H)δh1 (λi ), |λj,h1 − λi | δh21 (λi ), j = i, · · · , i + q − 1.
(4.23) (4.24) (4.25)
Let εh1 (λi ) := δh1 (λi ). From (4.23)-(4.25) and Theorem 4.1, we have εhk+1 (λj ) ηa (H)εhk (λj ) + δhk+1 (λj ),
for 1 ≤ k ≤ n − 1.
(4.26)
Then by recursive relation and based on the proof in Theorem 4.1, (4.2), (4.26) and Cβηa (H) < 1, we have εhn (λi ) ηa (H)εhn−1 (λi ) + δhn (λi ) ηa (H)2 εhn−2 (λi ) + ηa (H)δhn−1 (λi ) + δhn (λi ) n n n−k ηa (H) δhk (λi ) = (β γ ηa (H))n−k δhn (λi ) k=1
k=1
β γ ηa (H) δh (λi ) 1 − β γ ηa (H) n δhn (λi ).
(4.27)
This is the estimate (4.20). From the proof of Theorem 4.1 and (4.20), we can obtain the desired results (4.21) and (4.22).
12
4.3
Work estimate of eigenvalue multigrid scheme
In this section, we turn our attention to the estimate of computational work for Algorithm 4.2. We will show that Algorithm 4.2 makes solving the eigenvalue problem need almost the same work as doing one integration process on the finest mesh. First, we define the dimension of each level linear finite element space as Nk := dimVhk . Then we have Nk ≈
1 d(n−k) β
Nn ,
k = 1, 2, · · · , n.
(4.28)
Theorem 4.3. Assume the eigenvalue problem solved in the coarse spaces VH and Vh1 need work O(MH ) and O(Mh1 ), respectively, and the work of Sloan iteration in each level space Vhk is O(Nkα ) (α ≥ 1) for j = i, · · · , i + q − 1 and k = 2, 3, · · · , n. We do Sloan iteration in Step 1 of Algorithm 4.1 in parallel and solve the eigenvalue problem in Step 2 of Algorithm 4.1 in one processor. Then the most work involved in each processor of Algorithm 4.2 is O(Nnα + MH log(Nk ) + Mh1 ). Furthermore, the complexity will be O(Nnα ) provided MH Nnα and Mh1 ≤ Nnα . Proof. Let Wk denote the work in each processor in the k-th finite element space Vhk for the k-th One Correction Step. Then with the definition of Algorithm 4.1, we have Wk = O(Nkα + MH ).
(4.29)
Iterating (4.29) and using the fact (4.28), we obtain Total Work =
n
n α Wk = O Mh1 + Nk + MH
k=1
k=2
= O Mh1 + (n − 1)MH + = O Mh1 + (n − 1)MH + =
Nkα
k=2
O(Nnα
n
n 1 αd(n−k) k=2
+ MH log Nn + Mh1 ).
β
Nnα
(4.30)
This is the desired result O(Nnα + MH log Nn + Mh1 ) and the one O(Nnα ) can be obtained by the conditions MH Nnα and Mh1 ≤ Nnα . Remark 4.2. From Theorem 4.3, we can find Algorithm 4.2 can make the complexity of the Fredhold integral eigenvalue problems become the complexity of the integral process on the finest mesh. If we have an efficient integral scheme with α = 1, we can make the computational complexity of Algorithm 4.2 become optimal.
13
5
Numerical examples
In this section, we shall provide with some numerical examples to illustrate the efficiency of our MFE method, for computing the Karhunen-L`oeve expansion of random fields.
5.1
The Gaussian covariance function
Our first example is the following FIEP with Gaussian covariance function, (x−y)2 Va (x, y)v(y)dy = λv(x), with Va (x, y) = e− L , x, y ∈ D.
(5.1)
D
It is well know that the Gaussian covariance function is differentiable everywhere, and the eigenvalue sequence admits a fast decay rate. In this example, the linear finite element is adopted. We fist set the domain D to be (−1, 1) × (−1, 1). The corresponding numerical results are shown in Figure 1 (left). A salient advantage of using finite element methods is that it is easy to handle general domain. We next consider a L-shape domain D = (−1, 1) × (−1, 1)\[0, 1] × [0, 1] for (5.1). The numerical results are also shown in Figure 1 (right). From Figure 1, we can find that Algorithm 4.2 can reach the fourth order convergence rate which validate our theoretical results. Eigenvalues errors on unit square domain
−3
Eigenvalues errors on L−shape domain
−2
10
10 Σ10 |λ −λ | j=1 j j,h
Σ15 |λ −λ | j=1 j
slope=−2
−4
10
10
−5
−4
10 Errors
Errors
10
−6
10
−7
−5
10
−6
10
10
−8
−7
10
10
−9
10
j,h
slope=−2
−3
−8
1
10
2
10
3
10
4
10
Number of elements
10
1
10
2
10
3
10
4
10
Number of elements
Figure 1: Errors for the first 10 eigenvalues by the multi-grid way: the left subfigure is on the unit square domain and the right subfigure is on the L-shape domain, number of elements ≈ 1/h2
14
5.2
The exponential covariance function
We next consider the FIEP with an exponential covariance function, i.e., Va (x, y)v(y)dy = λv(x), with Va (x, y) = e−c|x−y| , x, y ∈ D.
(5.2)
D
The above exponential kernel is often used on benchmark tests, since closed formulas are available for its eigenvalues and eigenfunctions. For example, let D ≡ (−1, 1), a random field a(x, ω) with exponential covariance function given by (5.2) admits the following Karhunen-L`oeve expansion [10]: a(x, ω) =
∞ ξn λn fn (x) + ξn∗ λ∗n fn∗ (x) ,
(5.3)
n=1
where 2c 2c , λ∗n = ∗ 2 , 2 +c (μn ) + c2 cos(μn x) cos(μ∗n x) , fn∗ (x) = . fn (x) = sin(2μ∗n ) n) 1 + 1 + sin(2μ ∗ 2μn 2μn λn =
μ2n
(5.4) (5.5)
The parameters μn , μ∗n are solutions of the following problems μ∗ + ctan(μ∗ ) = 0.
c − μtan(μ) = 0,
(5.6)
In this example, we choose c = 1. The piecewise constant finite element method is adopted to compute the eigenvalue problem (5.2). The corresponding numerical results are shown in Figure 2 which also confirm the theoretical second order convergence.
5.3
A non-stationary case: the Wiener-Levy process
Finally, we consider the covariance function of the Wiener-Levy process which reads Va (x, y) = min{x, y},
x, y ∈ (0, 1) × (0, 1).
(5.7)
This example is provided to demonstrate that our method can be used as a simulation tool for non-stationary processes. The resulting normalized eigenfunctions and eigenvalues of the corresponding FIEP are [10]: x √ 4 √ , f (x) = 2sin . (5.8) λn = 2 n π (2n + 1)2 λn In this example, the eigenvalue problem (5.7) is discretized by the piecewise constant finite element method. We investigate the numerical results for the coarse 15
Eigenpair errors for H=1/20
−3
Eigenpair errors for H=1/40
−4
10
10 Σ11 |λ −λ | j=1 j j,h Σ11 |λ −λdir| j=1 j j,h
Σ11 |λ −λ | j=1 j
slope=−2
−4
10
j,h
Σ11 |λ −λdir| j=1 j j,h slope=−2 −5
Errors
Errors
10 −5
10
−6
10 −6
10
−7
10
−7
1
2
10
10
10
3
10
1
2
10
10
Number of elements
3
10
Number of elements
Figure 2: Errors for the first 11 eigenvalues by the multi-grid way with H = 1/20 (left) or H = 1/40 (right): λj,h is obtained by the multilevel method and λdir j,h denotes the approximation by the direct eigenvalue solving, number of elements ≈ 1/h mesh H = 1/20 and H = 1/40 for the biggest 11 eigenvalues and the associated first 5 eigenfunctions. The results are shown in Figure 3 (left for H = 1/20 and right for H = 1/40). It can be seen that the convergence orders of eigenvalue and eigenfunction approximations correspond to the theoretical prediction, respectively. Eigenpair errors for H=1/20
0
10
−1
−1
10
10
−2
−2
10 Errors
10 Errors
Eigenpair errors for H=1/40
0
10
−3
10
−4
−3
10
−4
10
10
Σ11 |λ −λ |
Σ11 |λ −λ | j=1 j j,h −5
10
1
10
j,h
Σ5j=1||uj−uj,h||0
−5
10
slope=−2 slope=−1
slope=−2 slope=−1
−6
10
j=1 j
Σ5 ||u −u || j=1 j j,h 0
−6
2
10
3
10
Number of elements
10
1
10
2
10
3
10
Number of elements
Figure 3: Errors for the first 11 eigenvalues and the first 5 eigenfunctions by the multi-grid way with H = 1/20 (left) or H = 1/40 (right), number of elements ≈ 1/h
16
6
Concluding remarks
In this paper, we propose a type of multilevel finite element method for the Fredholm integral eigenvalue problems based on the multilevel correction scheme. In this multilevel finite element method, we only need to solve integral eigenvalue problems in the coarsest finite element space and the Sloan iteration for the integral equation. The multilevel correction idea can also be extended to the collocation method for the integral eigenvalue problems which may be our future work.
Acknowledgment This first author is supported in part by the National Science Foundation of China (NSFC 91330202, 11371026, 11001259, 11201501, 2011CB309703), the National Center for Mathematics and Interdisciplinary Science (CAS), and the President Foundation of AMSS-CAS. The second author is supported by the National Science Foundation of China (No.91130003 and No.11201461).
References [1] K. Atkinson, The Numerical Solution of Integral Equations of the Second Kind, Cambridge University Press, Cambridge, 1997. [2] I. Babuˇska, F. Nobile and R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM J. Numer. Anal., 45 (2007), 1005-1034. [3] I. Babuˇska and J. E. Osborn, Finite element-Galerkin approximation of the eigenvalues and eigenvectors of selfadjoint problems, Math. Comp. 52 (1989), 275-297. [4] I. Babuˇska and J. E. Osborn, Eigenvalue Problems, In Handbook of Numerical Analysis, Vol. II, (Eds. P. G. Lions and Ciarlet P.G.), Finite Element Methods (Part 1), North-Holland, Amsterdam, 641-787, 1991. [5] J. B¨ack, F. Nobile, L. Tamellini, and R. Tempone, On the optimal polynomial approximation of stochastic PDEs by Galerkin and Collocation methods, Math. Mod. Methods Appl. Sci., 22(9) (2012), 1-33. [6] H. Brunner, A. Iserles, S. P. Nørsett, The spectral problem for a class of highly oscillatory Fredholm integral operators, IMA J. Numer. Anal., 30(1) (2010), 108-130.
17
[7] H. Brunner, Collocation Methods for Volterra Integral and Related Functional Equations, Cambridge University Press, Cambridge, 2004. [8] F. Chatelin, Spectral Approximation of Linear Operators, Academic Press Inc, New York, 1983. [9] F. R. DiNapoli, F. H. Middleton, Surface-wave propagation in a continuously stratified medium, J. Acoust. Soc. Am., 39(5A) (1966), 899-903. [10] R. Ghanem, P. Spanos, Stochastic Finite Element: A Spectral Aproach, Springer, New York, 1991. [11] M. Girolami, Orthogonal series density estimation and the kernel eigenvalue problem, Neural Comput., 14(3) (2002), 669-688. [12] H. B. Keller, On the accuracy of finite difference approximations to the eigenvalues of differential and integral operators, Numer. Math., 7(5) (1965), 412-419. [13] K. Khare, Sampling theorem, bandlimited integral kernels and inverse problems, Inverse Problems, 23(4) (2007), 1395-1416. [14] Q. Lin, Some problems concerning approximate solutions of operator equations, Acta Math. Sinica, 22 (1979), 219-230 (Chinese). [15] Q. Lin, H. Xie, A multi-level correction scheme for eigenvalue problems, Math. Comput., http://dx.doi.org/10.1090/S0025-5718-2014-02825-1. [16] F. Nobile, R. Tempone and C. Webster, A sparse grid stochastic collocation method for partial differential equations with random input data, SIAM J. Numer. Anal., 46(5) (2008), 2309-2345 [17] S. P. Oliveira and J. S. Azevedob, Spectral element approximation of Fredholm integral eigenvalue problems, J. Comput. Appl. Math., 257 (2014), 46-56. [18] J. Parvizian, A. D¨ uster and E. Rank, Finite cell method, Computational Mechanics, 41 (2007), 121-133. [19] K. K. Phoon, S.P. Huang, S. T. Quek, Implementation of Karhunen–Lo`eve expansion for simulation using a wavelet-Galerkin scheme, Probab. Eng. Mech., 17(3) (2002), 293-303. [20] F. Riesz and B. Sz-Nagy, Functional Analysis, Dover, 1990. [21] C. Schwab and R. A. Todor, Karhunen-Lo`eve approximation of random fields by generalized fast multipole methods, J. Comput. Phys., 217 (2006), 100-122. [22] I. Sloan, Iterated Galerkin method for eigenvalue problems, SIAM J. Numer. Anal., 13(5) (1976), 753-760. 18
[23] T. Tang and T. Zhou, Discrete least square projection in unbounded domain with random evaluations and its application to parametric uncertainty quantification, SIAM J. Sci. Comput., 36(5) (2014), A2272-A2295. [24] T. Tang and T. Zhou, Convergence analysis for stochastic collocation methods to scalar hyperbolic equations, Commun. Comput. Phys., 8 (2010), 226-248. [25] R. Todor, Robust eigenvalue computation for smoothing operators, SIAM J. Numer. Anal., 44(2) (2006), 865-878. [26] N. Wiener, The homogeneous chaos, Am. J. Math., 60 (1938), 897-936. [27] H. Xie, A type of multilevel method for the Steklov eigenvalue problem, IMA J. Numer. Anal., 34(2) (2014), 592-608. [28] H. Xie, A multigrid method for eigenvalue problem, J. Comput. Phys., 274 (2014), 550-561. [29] D. Xiu and J. S. Hesthaven, High-order collocation methods for differential equations with random inputs, SIAM J. Sci. Comput, 27 (2005), 1118-1139. [30] D. Xiu and G. E. Karniadakis, The Wiener-Askey polynomial chaos for stochastic differential equations, SIAM J. Sci. Comput., 24(2) (2002), 619-644. [31] J. Xu and A. Zhou, A two-grid discretization scheme for eigenvalue problems, Math. Comput., 70(233) (2001), 17-25. [32] T. Zhou and T. Tang, Galerkin methods for stochastic hyperbolic problems using Bi-orthogonal polynomials, J. Sci. Comput., 51 (2012), 274-292. [33] T. Zhou, A. Narayan and Z. Xu, Multivariate discrete least-squares approximations with a new type of collocation grid, SIAM J. Sci. Comput., 36(5) (2014), A2401-A2422.
19