# Multiscale Discretization Scheme Based on the Rayleigh Quotient Iterative Method for the Steklov Eigenvalue Problem

Authors: Hai Bi; Yidu Yang
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

This paper discusses efficient numerical methods for the Steklov eigenvalue problem and establishes a new multiscale discretization scheme and an adaptive algorithm based on the Rayleigh quotient iterative method. The efficiency of these schemes is analyzed theoretically, and the constants appeared in the error estimates are also analyzed elaborately. Finally, numerical experiments are provided to support the theory.


---

## Body

## 1. Introduction

Steklov eigenvalue problems have several deep applications both in physical and mechanical fields. For instance, they are found in the study of surface waves (see [1]), in the analysis of stability of mechanical oscillators immersed in a viscous fluid (see [2]), and in the study of the vibration modes of a structure in contact with an incompressible fluid (see, e.g. [3]). Thus, numerical methods for Steklov eigenvalue problems have attracted more and more scholars’ attention in recent years, for example, see [4–12].However, in practical applications, it is a challenging task to adopt efficient numerical methods to reduce the computational costs, such as the CPU time and the storage requirement, without any loss of the accuracy of numerical approximation. A two-grid finite element discretization scheme is considered one of these efficient methods. This discretization technique was first introduced by Xu [13, 14] for nonsymmetric and nonlinear elliptic problems, and due to its outstanding performance in computation, it has been successfully applied and further investigated for many other problems, for example, Poisson eigenvalue equations and integral equations in [15], nonlinear eigenvalue problems in [16], Schro¨dinger equation in [17, 18], Stokes equations in [19], and so forth. As for the Steklov eigenvalue problem, Li and Yang [20], Bi and Yang [21] established conforming and nonconforming finite element two-grid discretization scheme, respectively, and obtained satisfactory results. So far, to the best of our knowledge, there have been developed two kinds of two-grid discretization scheme for eigenvalue problems. The first type is contributed by Xu and Zhou [15], as just mentioned, and the other is proposed by Yang and Bi in [22], which has been applied to elliptic eigenvalue problems [22] and eigenvalue problem of electric field [23]. Compared with the first type, the second one is based on the shifted-inverse power method which makes itself more efficient, and starting from it, one can construct multiscale discretization scheme and parallel algorithm. So, this paper aims to establish a multi-scale discretization scheme based on the shifted-inverse power method, and an adaptive algorithm for the Steklov eigenvalue problem. The special work of this paper is as follows.(1) It is well known that the Rayleigh quotient iteration method is a basic approach for solving matrix eigenvalue problems (see Algorithm  27.3 in [24]). The multi-scale discretization scheme established in this paper is a combination of the finite element method and the Rayleigh quotient iteration method. We prove that this scheme is efficient, and analyze elaborately that the constants appeared in the error estimates are independent of mesh size and the iterative time. It can be seen easily that a two-grid discretization scheme is a special case of the multi-scale algorithm by taking the iterative time equal to 2 (see Scheme 2). From the theoretical analysis (see Theorem 3.2) we know that when taking H=𝒪(h3), the resulted solution can maintain an asymptotically optimal accuracy, while with the schemes in [20, 21], the resulted solution maintains an asymptotically optimal accuracy by taking H=𝒪(h).(2) The multi-scale discretization scheme is actually an iterative process. Since the approximate eigenvalue sequence obtained by the scheme converges to the exact eigenvalue (see Theorem3.1, Corollary 3.4) and the constants appeared in its error estimates are independent of the mesh size and the iterative time, we can stop the iteration when the error between two neighboring approximate eigenvalues is less or equal to a given error tolerance. Thus, we establish an adaptive algorithm (see Scheme 3). With this adaptive scheme, first we solve an eigenvalue problem on a coarse grid, in each step after that we only need to solve a linear algebraic system on a fine grid. Compared with the existing adaptive method (e.g., see Algorithm  1 in [11]) which computes an eigenvalue problem in each step, our approach reduces the computational complexity.The rest of this paper is organized as follows. In the subsequent section, some preliminaries needed in this paper are presented. In Section3, a multi-scale and a two-grid discretization scheme are proposed, and the error estimates are also given. In Section 4, an adaptive algorithm is established. Finally, numerical experiments are provided to support our theoretical analysis.


## 2. Preliminaries

LetHs(Ω) and Hs(∂Ω) denote Sobolev spaces on Ω and ∂Ω with real order s, respectively. The norms in Hs(Ω) and Hs(∂Ω) are denoted by ∥·∥s and ∥·∥s,∂Ω, respectively.Consider the model problem(2.1)-Δu+u=0inΩ,∂u∂n=λuon∂Ω,
where Ω⊂R2 is a polygonal domain and ∂u/∂n is the outward normal derivative on ∂Ω.The weak form of (2.1) is given by the following: find λ∈R, 0≠u∈H1(Ω), such that
(2.2)a(u,v)=λb(u,v),∀v∈H1(Ω),
where
(2.3)a(u,v)=∫Ω∇u∇v+uvdx,b(u,v)=∫∂Ωuvds.
It is easy to know that a(·,·) is a symmetric, continuous, and H1(Ω)-elliptic bilinear form on H1(Ω)×H1(Ω). So, we use a(·,·) and ∥·∥a=a(·,·)=∥·∥1 as the inner product and norm on H1(Ω), respectively. We use b(·,·) and ∥·∥b=b(·,·) as the inner product and norm on L2(∂Ω), respectively.From [25], we know that (2.2) has a countable sequence of positive eigenvalues
(2.4)0<λ1≤λ2≤λ3≤⋯→+∞
(here each eigenvalue occurs as many times as given by its multiplicity), and the corresponding eigenfunctions uk∈H1+r(Ω)(k=1,2,…), where r=1 if Ω is convex, and r<(π/ω) which can be arbitrarily close to π/ω when Ω is concave (with ω being the largest inner angle of Ω, ω<2π).Letπh be a regular triangulation of Ω with the mesh diameter h, and let Sh be a piecewise polynomial space of degree m (m≥1) defined on πh.The conforming finite element approximation of (2.2) is the following: find λh∈R, 0≠uh∈Sh⊂H1(Ω), such that
(2.5)a(uh,v)=λhb(uh,v),∀v∈Sh.It is well known that (2.5) has a finite sequence of eigenvalues
(2.6)0<λ1,h≤λ2,h≤⋯λNh,h(Nh=dimSh),
and their corresponding eigenfunctions are uk,h(k=1,2,⋯,Nh).Consider the following source problem (2.7) associated with (2.2) and the approximate source problem (2.8) associated with (2.5), respectively.Findu∈H1(Ω), such that
(2.7)a(u,v)=b(f,v),∀v∈H1(Ω).
Find uh∈Sh, such that
(2.8)a(uh,v)=b(f,v),∀v∈Sh.Several regular estimates of the Steklov eigenvalue problem are presented in the following lemma, which will be used in the sequel.Lemma 2.1.
Iff∈L2(∂Ω), then there exists a unique solution u∈H1+(r/2)(Ω) to (2.7) and
(2.9)∥u∥1+(r/2)≤Cp∥f∥0,∂Ω.
If f∈H1/2(∂Ω), then there exists a unique solution u∈H1+r(Ω) to (2.7) and
(2.10)∥u∥1+r≤Cp∥f∥1/2,∂Ω.
If f∈L2(∂Ω), then there exists a unique solution u∈H1(Ω) to (2.7) and
(2.11)∥u∥a≤Cp∥f∥0,∂Ω,
where Cp appeared in (2.9), (2.10), and (2.11) stands for a positive constant but may take different values.Proof.
As regards the proof of (2.9) and (2.10), see [26]. Since a(·,·) is H1(Ω)-elliptic, from (2.7), it is easy to know that (2.11) is valid.Thus, from (2.7), we can define the operator A:L2(∂Ω)→H1(Ω) by
(2.12)a(Af,v)=b(f,v),∀v∈H1(Ω).Similarly, from (2.8), we define the operator Ah:L2(∂Ω)→Sh⊂H1(Ω) by
(2.13)a(Ahf,v)=b(f,v),∀v∈Sh.
It is obvious that A:H1(Ω)→H1(Ω) is a self-adjoint operator. In fact, for any u,v∈H1(Ω), a(Au,v)=b(u,v)=b(v,u)=a(Av,u)=a(u,Av). Analogously, Ah is also a self-adjoint operator. Observe that Af and Ahf are the exact solution and the finite element solution of (2.7), respectively, and a(Af-Ahf,v)=0,∀v∈Sh⊂H1(Ω). Define the Ritz-Galerkin projection operator Ph:H1(Ω)→Sh by
(2.14)a(u-Phu,v)=0,∀u∈H1(Ω),∀v∈Sh.
Thus, ∀f∈H1(Ω),
(2.15)a(Ahf-Ph(Af),v)=a(Ahf-Af+Af-Ph(Af),v)=0,∀v∈Sh.
Therefore, Ahf=PhAf,∀f∈H1(Ω), then Ah=PhA. From Lemma 2.1 and the interpolation error estimate, we have
(2.16)∥Ah-A∥a=supg∈H1(Ω)∥(Ah-A)g∥a∥g∥a=supg∈H1(Ω)∥PhAg-Ag∥a∥g∥a≤supg∈H1(Ω)CIhr∥Ag∥1+r∥g∥a≤supg∈H1(Ω)CIhr∥g∥a∥g∥a=CIhr→0(h→0),
where CI is the interpolation constant. It is clear that Ah is a finite rank operator, then A is a completely continuous operator. From [25, 27], we know that (2.2) and (2.5) have the following equivalent operator forms, respectively:
(2.17)Au=μu,Ahuh=μhuh,
where μ=1/λ,μh=1/λh. In this paper, μk, and μk,h, λk and λk,h are all called eigenvalues.Suppose that the algebraic multiplicity ofμk is equal to q, μk=μk+1=⋯=μk+q-1. Let M(μk) be the space spanned by all eigenfunctions corresponding to μk of A, and let Mh(μk) be the direct sum of eigenspaces corresponding to all eigenvalues of Ah that converge to μk. Let M^(μk)={v:v∈M(μk),∥v∥a=1}, M^h(μk)={v:v∈Mh(μk),∥v∥a=1}. We also write M(λk)=M(μk), Mh(λk)=Mh(μk), M^(λk)=M^(μk), and M^h(λk)=M^h(μk).Let(2.18)δh(λk)=supw∈M^(λk)infv∈Sh∥w-v∥a.From Lemma2.1, we know that M^(λk)⊂H1+r(Ω), then by the interpolation error estimate, we have
(2.19)δh(λk)≤CIhrsupw∈M^(λk)∥w∥1+r.Lemma 2.2.
Letλk,h and λk be the kth eigenvalue of (2.5) and (2.2), respectively. Then
(2.20)λk≤λk,h≤λk+C1δh2(λk),
for any eigenfunction uk,h corresponding to λk,h, satisfying ∥uk,h∥a=1, there exists uk∈M^(λk) such that
(2.21)∥uk,h-uk∥a≤C2δh(λk),
and for any uk∈M^(λk), there exists uh∈Mh(λk) such that
(2.22)∥uh-uk∥a≤C3δh(λk),
where Ci,i=1,2,3 are constants independent of h.Proof.
By the argument in [25, 27], we can obtain the desired results.The following lemma states a crucial property (but straightforward) of eigenvalue and eigenfunction approximation.Lemma 2.3.
Let(λ,u) be an eigenpair of (2.2), then for any w∈H1(Ω),∥w∥b≠0, the Rayleigh quotient a(w,w)/∥w∥b2 satisfies
(2.23)a(w,w)∥w∥b2-λ=∥w-u∥a2∥w∥b2-λ∥w-u∥b2∥w∥b2.Proof.
See, for instance, Lemma  9.1 in [27] for details.We also need the following basic estimate of shifted-inverse power method (see Theorem  3.2 in [22]).Lemma 2.4.
Let(μ0,u0) be an approximation for (μk,uk), where μ0 is not an eigenvalue of Ah, and u0∈Sh with ∥u0∥a=1. Suppose that maxk≤j≤k+q-1|(μj,h-μk,h)/(μ0-μj,h)|≤1/2, dist(u0,Mh(μk))≤1/2, |μ0-μj,h|≥ρ/2 for j≠k,k+1,…,k+q-1, and u∈Sh,ukh∈Sh satisfy
(2.24)(μ0-Ah)u=u0,ukh=u∥u∥a.Then(2.25)dist(ukh,M^h(μk))≤16ρ|μ0-μk,h|dist(u0,Mh(μk)),
where ρ=minμj≠μk|μj-μk| be the separation constant of the eigenvalue μk.For the multi-scale discretization scheme established in this paper, the conditions of Lemma2.4 are satisfied which will be verified in the proof of Theorem 3.1.


## 3. Multi-Scale Discretization Scheme

In this section, we combine the finite element method with the Rayleigh quotient iteration method and establish a multi-scale discretization scheme. Let{πhi}1l be a family of regular meshes, hi-1≫hi, and let {Shi}1l be the conforming finite element spaces defined on {πhi}1l, and let πH=πh1, SH=Sh1, πh=πhl, and Sh=Shl.Scheme 1 (multi-scale discretization scheme).Step 1.
Solve (2.2) on the πH: find λk,H∈R,uk,H∈SH such that ∥uk,H∥a=1 and
(3.1)a(uk,H,v)=λk,Hb(uk,H,v),∀v∈SH.Step 2.
Execute the assignments:ukh1⇐uk,H, λkh1⇐λk,H, i⇐2.Step 3.
Solve a linear system on theπhi: find u~∈Shi such that
(3.2)a(u~,v)-λkhi-1b(u~,v)=b(ukhi-1,v),∀v∈Shi.
And set ukhi=u~/(∥u~∥a).Step 4.
Compute the Rayleigh quotient(3.3)λkhi=a(ukhi,ukhi)b(ukhi,ukhi).Step 5.
Ifi=l, then output (λkhl,ukhl), that is, (λkh,ukh), stop. Else, i⇐i+1, and return to Step 3.Theorem 3.1.
Let(λkhl,ukhl) be obtained by Scheme 1. Assume that (λkhl-1,ukhl-1) approximate an eigenpair, (λk,u-), of (2.2), u-∈M^(λk), and λkhl-1-λk is a small quantity of lower order than λk,hl-λk. Then there exists uk∈M(λk) such that
(3.4)∥ukhl-uk∥a≤48ρC4C5λk|λkhl-1-λk|dist(ukhl-1,M^(λk))+3C2qδhl(λk),(3.5)|λkhl-λk|≤2λk(1+C6λk)∥ukhl-uk∥a2,l≥2,
where the constants C4,C5, and C6 are stated in (3.6), (3.11), and (3.29), respectively.Proof.
We use Lemma2.4 to complete the proof. First, we will verify that the conditions of Lemma 2.4 are satisfied.
Selectμ0=1/λkhl-1, and u0=λkhl-1Ahlukhl-1/∥λkhl-1Ahlukhl-1∥a. From (2.16), we know that ∥Ahl-A∥a→0(hl→0), then there exists a constant C4 independent of hl and l such that
(3.6)∥Ahlf∥a≤C4∥f∥a,∀f∈H1(Ω).
Thus, by the assumption, we deduce that
(3.7)∥λkhl-1Ahlukhl-1-u-∥a≤∥λkhl-1Ahlukhl-1-λkhl-1Ahlu-+λkhl-1Ahlu--λkAhlu-+λkAhlu--λkAu-∥a≤C4λkhl-1∥ukhl-1-u-∥a+C4|λkhl-1-λk|·∥u-∥a+λk∥(Ahl-A)u-∥a.
Note that in any normed space, for any nonzero u,v∈Sh, there holds
(3.8)∥u∥u∥-v∥v∥∥≤2∥u-v∥∥u∥,∥u∥u∥-v∥v∥∥≤2∥u-v∥∥v∥.
Hence, we have
(3.9)dist(u0,M^(λk))≤∥u0-u-∥u-∥a∥a≤2∥λkhl-1Ahlukhl-1-u-∥a≤2{C4λkhl-1∥ukhl-1-u-∥a+C4|λkhl-1-λk|·∥u-∥a+λk∥(Ahl-A)u-∥a}.
Using the triangle inequality and (2.22), we get
(3.10)dist(u0,Mhl(λk))≤dist(u0,M^(λk))+C3δhl(λk).
It follows from (2.20) that λk,hl→λk(hl→0), then by the assumption, we have
(3.11)|μ0-μk,hl|=|λkhl-1-λk+λk-λk,hlλk,hlλkhl-1|≤C5|λkhl-1-λk|,
where C5 is a constant independent of hl. By (2.13), we see that Step 3 in Scheme 1 is equivalent to the following:
(3.12)a(u~,v)-λkhl-1a(Ahlu~,v)=a(Ahlukhl-1,v),∀v∈Sh,
u
k
h
l
=
u
~
/
∥
u
~
∥
a, that is,
(3.13)(1λkhl-1-Ahl)u~=1λkhl-1Ahlukhl-1,ukhl=u~∥u~∥a.
Notice that(1/λkhl-1)Ahlukhl-1=∥(1/λkhl-1)Ahlukhl-1∥au0 differs from u0 by only a constant, then Step 3 is equivalent to
(3.14)(1λkhl-1-Ahl)u~=u0,ukhl=u~∥u~∥a.
In fact, it is obvious that ukhl obtained by the above two formulae are the same. When hl-1 is small enough, noting that hl≪hl-1, from (3.10) and (3.9), we get
(3.15)dist(u0,Mhl(μk))≤12.
From (2.20), having in mind that λk=λk+1=⋯=λk+q-1, we know that
(3.16)λj,hl-λk≤C1δhl2(λk)(j=k,k+1,…,k+q-1),
then
(3.17)|μj,hl-μk,hl|=|λk,hl-λj,hlλk,hlλj,hl|=|λk,hl-λk+λj-λj,hlλk,hlλj,hl|.
Combining (3.17) with (3.11) and noting that the quantity on the right-hand side of (3.17) is an infinitesimal of higher order comparing with λkhl-1-λk, we derive
(3.18)maxk≤j≤k+q-1|μj,hl-μk,hlμ0-μj,hl|≤12.
Since ρ is the separation constant, hl-1 is small enough, and hl≪hl-1, there holds
(3.19)|μ0-μj,hl|≥ρ2,j≠k,k+1,…,k+q-1.
From the above arguments, we see that the conditions of Lemma 2.4 hold.
Next, we will prove that (3.4) and (3.5) are valid.
Substituting (3.10) and (3.11) into (2.25), we obtain
(3.20)dist(ukhl,M^hl(μk))≤16ρC5|λkhl-1-λk|(dist(u0,M^(λk))+C3δhl(λk)).
Let the eigenvectors {uj,hl}kk+q-1 be an orthonormal basis of Mhl(λk). Note that
(3.21)dist(ukhl,Mhl(λk))=∥ukhl-∑j=kk+q-1‍a(ukhl,uj,hl)uj,hl∥a.
Let
(3.22)u*=∑j=kk+q-1‍a(ukhl,uj,hl)uj,hl,
it follows directly from (3.20) that
(3.23)∥ukhl-u*∥a≤16ρC5|λkhl-1-λk|(dist(u0,M^(λk))+C3δhl(λk)).
By Lemma 2.2, there exists {uj0}kk+q-1⊂M^(λk) making uj,hl-uj0 satisfy (2.21). Let
(3.24)uk=∑j=kk+q-1‍a(ukhl,uj,hl)uj0,
then uk∈M(λk). Using (2.21), we deduce that
(3.25)∥u*-uk∥a=∥∑j=kk+q-1a(ukhl,uj,hl)(uj,hl-uj0)∥a≤(∑j=kk+q-1∥uj,hl-uj0∥a2)1/2≤C2∑j=kk+q-1δhl(λj)≤C2qδhl(λk).
Combining (3.23) with the above inequality, we have
(3.26)∥ukhl-uk∥a≤16ρC5|λkhl-1-λk|(dist(u0,M^(λk))+C3δhl(λk))+C2qδhl(λk)≤16ρC5|λkhl-1-λk|dist(u0,M^(λk))+2C2qδhl(λk).
It is obvious that there exists uk'∈M^(λk) such that
(3.27)∥ukhl-1-uk'∥a=dist(ukhl-1,M^(λk)),λkAuk'=uk'∈M^(λk), and ∥(Ahl-A)uk'∥a≤(1/λk)δhl(λk).
Sinceu0=λkhl-1Ahlukhl-1/∥λkhl-1Ahlukhl-1∥a, from (2.23) and (3.8), we derive
(3.28)dist(u0,M^(λk))≤2dist(λkhl-1Ahlukhl-1,M^(λk))≤2∥λkhl-1Ahlukhl-1-λkAuk'∥a≤2∥(λkhl-1-λk)Ahlukhl-1∥a+2∥λkAhl(ukhl-1-uk′)∥a+2∥λk(Ahl-A)uk'∥a≤2C4|λkhl-1-λk|+2C4λkdist(ukhl-1,M^(λk))+2δhl(λk)≤3C4λkdist(ukhl-1,M^(λk))+2δhl(λk).
Substituting (3.28) into (3.26), we obtain (3.4). Equation (3.4) indicates that ukhl converges to uk in the sense of norm ∥·∥a, then from the trace theorem, we know that ukhl converges to uk in the sense of norm ∥·∥b; thus, 1/∥ukhl∥b2→1/∥uk∥b2=λk. Therefore, when hl is small enough, we have 1/∥ukhl∥b2≤2λk. Picking w=ukhl in (2.23), we get
(3.29)|λkhl-λk|≤∥ukhl-uk∥a2∥ukhl∥b2+λk∥ukhl-uk∥b2∥ukhl∥b2≤2λk∥ukhl-uk∥a2+2λk2∥ukhl-uk∥b2≤2λk(1+C6λk)∥ukhl-uk∥a2,
where the last inequality in the above holds due to the trace theorem, and the constant C6 is independent of hl.Setl=2 and denote H=h1,h=h2, then we immediately get the following two-grid discretization scheme based on the shifted-inverse power method.Scheme 2 (two-grid discretization scheme).Step 1.
Solve (2.2) on a coarse grid πH: find λk,H∈R,uk,H∈SH, such that ∥uk,H∥a=1 and
(3.30)a(uk,H,v)=λk,Hb(uk,H,v),∀v∈SH.Step 2.
Solve a linear system on a fine gridπh: find u~∈Sh, such that
(3.31)a(u~,v)-λk,Hb(u~,v)=b(uk,H,v),∀v∈Sh.
Set ukh=u~/∥u~∥a.Step 3.
Compute the Rayleigh quotient(3.32)λkh=a(ukh,ukh)b(ukh,ukh).Then from Theorem3.1, we have the following error estimates for Scheme 2.Theorem 3.2.
Let(λkh,ukh) be obtained by Scheme 2. Then when H is properly small and h≪H, there exists uk∈M(λk) such that
(3.33)∥ukh-uk∥a≤48ρC4C5λk(λk,H-λk)dist(uk,H,M^(λk))+3C2qδh(λk),(3.34)|λkh-λk|≤2λk(1+C6λk)∥ukh-uk∥a2.Proof.
Since(λk,H,uk,H) is obtained by Step 1 in Scheme 2, that is, (λk,H,uk,H) is a solution of (2.5) on the coarse grid πH, it certainly approximates an eigenpair, (λk,u-), of (2.2), u-∈M^(λk). Note that λk,h is an eigenvalue of (2.5) on the fine grid πh, so when H is properly small and h≪H, it is valid that λk,H-λk is a small quantity of lower order than λk,h-λk. The above arguments imply that the conditions of Theorem 3.1 hold for l=2. Therefore, the desired results follow from Theorem 3.1.Remark 3.3.
It follows directly from (2.20) that
(3.35)λk,H-λk≤C1δH2(λk).And from (2.21), we know that dist(uk,H,M^(λk))≤C2δH(λk), then from Theorem 3.2, we get
(3.36)∥ukh-uk∥a≤48ρC1C2C4C5λkδH3(λk)+3C2qδh(λk).
Substituting (2.19) into (3.36), we obtain
(3.37)∥ukh-uk∥a≤48ρC2C4C5λkC1CI3H3r(supw∈M^(λk)∥w∥1+r)3+3C2qCIhrsupw∈M^(λk)∥w∥1+r.
From (3.37) and (3.34), we can see that when we take H=𝒪(h3), the resulted solution, (λkh,ukh), can maintain an asymptotically optimal accuracy.Corollary 3.4.
Suppose thathi=hi-1ti,ti∈(1,3),i=2,3,…,supiti<3. By using the triangle linear conforming element with regular meshes, when h1, that is, H is properly small, there exists uk∈M(λk) such that the following error estimates for Scheme 1 hold:
(3.38)∥ukhl-uk∥a≤6qC2CIhlrsupw∈M^(λk)∥w∥1+r,(3.39)|λkhl-λk|≤2λk(1+C6λk)∥ukhl-uk∥a2,l≥2.Proof.
We prove the results by using induction.
Forl=2, we estimate λk,H-λk and dist(uk,H,M^(λk)) in (3.33), respectively.
Combining (3.35) and (2.19) yields
(3.40)λk,H-λk≤C1CI2H2r(supw∈M^(λk)∥w∥1+r)2.
From Lemma 2.2, and (3.8) we have
(3.41)dist(uk,H,M^(λk))≤C2δH(λk)≤C2CIHrsupw∈M^(λk)∥w∥1+r.
Substituting the above estimates into (3.33), we get
(3.42)∥ukh-uk∥a≤48ρC1C2C4C5λkCI3H3r(supw∈M^(λk)∥w∥1+r)3+3qC2CIhrsupw∈M^(λk)∥w∥1+r.
By the assumption, h=Ht2,t2∈(1,3), we can see that the first term on the right-hand side of (3.42) is a small quantity of higher order than the second term, hence,
(3.43)∥ukh-uk∥a≤6qC2CIhrsupw∈M^(λk)∥w∥1+r.
Suppose that (3.38), (3.39) hold for l-1, that is, there holds
(3.44)∥ukhl-1-uk∥a≤6qC2CIhl-1rsupw∈M^(λk)∥w∥1+r,(3.45)|λkhl-1-λk|≤2λk(1+C6λk)∥ukhl-1-uk∥a2,l≥2.
By (3.8), we get
(3.46)dist(ukhl-1,M^(λk))≤2∥ukhl-1-uk∥a.
Substituting (3.45) and (3.46) into (3.4), and noting (3.44), we deduce,
(3.47)∥ukhl-uk∥a≤192ρC4C5λk2(1+C6λk)(6qC2CIsupw∈M^(λk)∥w∥1+r)3hl-13r+3C2qδhl(λk)≤192ρC4C5λk2(1+C6λk)(6qC2CIsupw∈M^(λk)∥w∥1+r)3hl-13r+3qC2CIhlrsupw∈M^(λk)∥w∥1+r≤6qC2CIhlrsupw∈M^(λk)∥w∥1+r,
that is, (3.38) holds. The inequality (3.39) follows directly from (3.29).From Theorem3.1 and Corollary 3.4, we can see that the error is independent of the mesh diameters and the iteration time, and Corollary 3.4 also shows that with Scheme 1, by using the linear element, the approximate eigenvalue sequence {λkhl} converges to {λk} as l→∞, and the convergence is quick. Hence, we have the reason to believe that |λkhl-λkhl+1|<|λkhl-λkhl-1|. A simple calculation shows that
(3.48)|λkhl-λk|≤|λkhl-λkhl+1|+|λkhl+1-λk|≤|λkhl-λkhl-1|+|λkhl+1-λk|.Note that in (3.48) |λkhl+1-λk| is a small quantity of higher order. So we can use |λkhl-λkhl-1|<ε (ε is a given tolerate error) as a criteria to stop the iteration. Scheme 1 is actually an iterative process. It is natural to establish an adaptive algorithm based on Scheme 1.Scheme 3 (adaptive algorithm).
Give an error toleranceε, an initial coarse grid πH, and a grid πh2 derived from πH.Step 1.
Solve (2.2) on the πH: find λk,H∈R,uk,H∈SH, such that ∥uk,H∥a=1 and
(3.49)a(uk,H,v)=λk,Hb(uk,H,v),∀v∈SH.Step 2.
Execute the assignments:ukh1⇐uk,H, λkh1⇐λk,H, i⇐2.Step 3.
Solve a linear system on theπhi: find u~∈Shi, such that
(3.50)a(u~,v)-λkhi-1b(u~,v)=b(ukhi-1,v),∀v∈Shi.
And set ukhi=u~/∥u~∥a.Step 4.
Compute the Rayleigh quotient(3.51)λkhi=a(ukhi,ukhi)b(ukhi,ukhi).Step 5.
If|λkhi-λkhi-1|<ε, then output (λkhi,ukhi), stop. Else, return to Step 6.Step 6.
Determineti+1 and hi+1=hiti+1, then refine πhi to get πhi+1. Set i⇐i+1, and return to Step 3.


## 4. Numerical Experiments

Example 4.1.
We compute the first four approximate eigenvalues of (2.1) with the triangle linear finite element by using Scheme 2 on Ω-=[0,1]×[0,1] and [0,1]×[0,1/2]∪[0,1/2]×[1/2,1], respectively, [9] and [28] have proved that the non-conforming EQ1rot element can provide the lower bound for the exact eigenvalues, and the minimum-maximum principle ensures that conforming finite elements can give the upper bound for the exact eigenvalues. So we compute the following ranges for the first four exact eigenvalues of (2.1) by the linear triangle element and the EQ1rot element.WhenΩ-=[0,1]×[0,1] the first four exact eigenvalues
(4.1)λ1≈0.24007909,λ2≈1.492303,λ3≈1.492303,λ4≈2.08265;
when Ω-=[0,1]×[0,1/2]∪[0,1/2]×[1/2,1] the first four exact eigenvalues
(4.2)λ1≈0.18296424,λ2≈0.8936,λ3≈1.68860,λ4≈3.2179.
When compute with Scheme 2, we adopt a uniform isosceles right triangulation along three directions to obtain the coarse grid πH, and refine the coarse grid in a uniform way (each triangle is divided into four congruent subtriangles) repeatedly to obtain the fine grid πh. Then we compute the approximate eigenvalues with MATLAB7.1. The numerical results are shown in Tables 1 and 2.Table 1
Numerical eigenvalues on the square domain[0,1]×[0,1] by Scheme 2: set h=2/512,h=𝒪(Ht),t∈(1,3].
H
h
λ
1
h
λ
2
h
λ
3
h
λ
4
h
2
/
8
2
/
512
0.2400791223
1.4923055029
1.4923059934
2.0826626453
2
/
16
2
/
512
0.2400791223
1.4923054995
1.4923059872
2.0826625905
2
/
32
2
/
512
0.2400791223
1.4923054994
1.4923059871
2.0826625901Table 2
Numerical eigenvalues on the L-shaped domain[0,1]×[0,1/2]∪[0,1/2]×[1/2,1] by Scheme 2:  set h=2/512,h=𝒪(Ht),t∈(1,3].
H
h
λ
1
h
λ
2
h
λ
3
h
λ
4
h
2
/
8
2
/
512
0.1829642799
0.8937364938
1.6886068115
3.2179250611
2
/
16
2
/
512
0.1829642799
0.8937364004
1.6886067424
3.2179003466
2
/
32
2
/
512
0.1829642799
0.8937363984
1.6886067421
3.2179002027Remark 4.2.
In Example4.1, we set the diameter of fine grid h=2/512 and H=𝒪(ht),t∈(1,3]. From Table 1 we can see that the approximate eigenvalues obtain an asymptotically optimal accuracy. While with the two-grid discretization scheme in [20], to achieve the same accuracy the relationship between the diameter of coarse grid and that of fine grid should satisfy H=𝒪(h2). For example, we carry out the computation on the same grids with h=2/512 and H=2/8, by our scheme the resulted solution λ1h=0.2400791223 which has 7 significant digits, and by the scheme in [20] the resulted solution λ~1h=0.2400791543 with 6 significant digits. Note that the computational complexities of these two methods are almost the same which indicates that our method is also efficient.Remark 4.3.
We program by using the finite element package of Chen [29] to solve Example 4.1. It costs 19.745s by using Scheme 2, with H=2/8 and h=2/512, to get the first approximate eigenvalue λ1h=0.2400791223 of (2.1) on the square [0,1]×[0,1]. To illustrate the efficiency of our scheme we compute the first approximate eigenvalue on the fine grid with mesh diameter h=2/512 directly, and it costs 75.267s to get the same approximation. We also compare the computation times of our approach and direct calculation on the L-shaped domain [0,1]×[0,1/2]∪[0,1/2]×[1/2,1]. By using Scheme 2 with H=2/8 and h=2/512, we spend 10.779s to obtain λ1h=0.1829642799; while computing on the fine grid with mesh size h=2/512 directly it costs 31.986s. From these comparisons we can see that our scheme is very efficient.Example 4.4.
We compute the first four approximate eigenvalues of (2.1) with the triangle linear finite element by using Scheme 3 on [0,1]×[0,1] and [0,1]×[0,1/2]∪[0,1/2]×[1/2,1], respectively. The results are listed in Tables 3 and 4, respectively. And in Tables 3 and 4, the symbol “-" means that we have obtained the approximate eigenvalues which meet the accuracy requirements and the iteration stops here.Table 3
Numerical eigenvalues on the square domainΩ-=[0,1]×[0,1] by Scheme 3: set h1=H=1/8,h2=h12,hi=hi-1/2,i=3,4,….
ε
=
5
×
1
0
-
7
ε
=
5
×
1
0
-
5
ε
=
5
×
1
0
-
5
ε
=
5
×
1
0
-
5
l
λ
1
h
l
λ
2
h
l
λ
3
h
l
λ
4
h
l
1
0.2402262809
1.5014059516
1.5032096484
2.1452661239
2
0.2400814379
1.4924542699
1.4924853962
2.0836406391
3
0.2400796738
1.4923409581
1.4923487554
2.0828955752
4
0.2400792326
1.4923125932
1.4923145438
2.0827091950
5
—
—
—
2.0826625901
6
—
—
—
—Table 4
Numerical eigenvalues on the L-shaped domainΩ-=[0,1]×[0,1/2]∪[0,1/2]×[1/2,1] by Scheme 3:  set h1=H=1/8,h2=h12,hi=hi-1/2,i=3,4,….
ε
=
5
×
1
0
-
7
ε
=
5
×
1
0
-
4
ε
=
5
×
1
0
-
5
ε
=
5
×
1
0
-
4
l
λ
1
h
l
λ
2
h
l
λ
3
h
l
λ
4
h
l
1
0.1831328879
0.9211580591
1.7114443522
3.3645126465
2
0.1829669801
0.8951626591
1.6889966003
3.2204599675
3
0.1829649244
0.8942351070
1.6887001313
3.2185057817
4
0.1829644089
0.8938763111
1.6886254813
3.2180214091
5
0.1829642799
—
1.6886067421
—
6
—
—
—
—Remark 4.5.
In the research of numerical methods for Steklov eigenvalue problems, the boundary element method and the finite element method have been studied. However, we have not seen any literatures on the finite difference method, a major numerical method, for Steklov eigenvalue problems. To apply the finite difference method to Steklov eigenvalue problems would be a new and interesting issue which is our next work.


## 5. Concluding Remarks

This paper discusses the Steklov eigenvalue problem, and establishes a new multi-scale discretization scheme and an adaptive algorithm. We prove that our approach is efficient. With the adaptive scheme, first we solve an eigenvalue problem on a coarse grid, in each step after that we only need to solve a linear algebraic system on a fine grid. Compared with the existing adaptive method which computes an eigenvalue problem in each step, our approach reduces the computational complexity. However, the a posteriori error indicator we used is a global estimate, thus, to establish a local a posteriori error estimate and improve our adaptive algorithm is our next goal.



---

*Source: 487207-2012-09-23.xml*