# Spectral Method with the Tensor-Product Nodal Basis for the Steklov Eigenvalue Problem

Authors: Xuqing Zhang; Yidu Yang; Hai Bi
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This paper discusses spectral method with the tensor-product nodal basis at the Legendre-Gauss-Lobatto points for solving the Steklov eigenvalue problem. A priori error estimates of spectral method are discussed, and based on the work of Melenk and Wohlmuth (2001), a posterior error estimator of the residual type is given and analyzed. In addition, this paper combines the shifted-inverse iterative method and spectral method to establish an efficient scheme. Finally, numerical experiments with MATLAB program are reported.


---

## Body

## 1. Introduction

Since Steklov eigenvalue problems have important physical background and many applications, for instance, they appear in the analysis of stability of mechanical oscillators immersed in a viscous fluid (see [1] and the references therein), in the analysis of the antiplane shearing on a system of collinear faults under slip-dependent friction law (see [2]), in the study of surface waves (see [3]), in the study of the vibration modes of a structure in contact with an incompressible fluid (see [4]) and vibration of a pendulum (see [5]), and in the eigen-oscillations of mechanical systems with boundary conditions containing frequency [6], and thus, numerical methods for solving Steklov eigenvalue problems have received increasing attention in recent years. Bramble and Osborn [7] studied the Galerkin approximation of a Steklov eigenvalue problem of nonself-adjoint second order elliptic operators in smooth domain, Andreev and Todorov [8] discussed the isoparametric finite element method for the approximation of the Steklov eigenvalue problem of second-order self-adjoint elliptic differential operators, Armentano and Padra [9] introduced and analyzed the conforming linear finite element approximation of the Steklov eigenvalue problem in a bounded polygonal domain, Alonso and Russo [10], Yang et al. [11], and Li et al. [12] studied nonconforming finite elements approximation of the Steklov eigenvalue problem, Li and Yang [13] and Bi and Yang [14] discussed a two-grid method of the conforming and non-conforming finite element method, respectively, Li et al. [15] studied the extrapolation and superconvergence of the Steklov eigenvalue problem, Tang et al. [16] studied the boundary element approximation, and Cao et al. [17] discussed multiscale asymptotic method for Steklov eigenvalue equations in composite medias. However, to the best of our knowledge, there have been no reports on spectral method for Steklov eigenvalue problems.Spectral method is an important numerical method for solving differential equations developed after finite-difference and finite-element methods. Because spectral method has the characteristics of superior accuracy, it is widely used in the field of meteorology, physics, and mechanics (see [18, 19]). This paper discusses spectral method with the tensor-product nodal basis at the Legendre-Gauss-Lobatto (LGL for short) points for solving the Steklov eigenvalue problem. We use spectral approximation theory (see [20]) to give a priori error estimates. Based on the work of [21], we discuss a posterior error estimator of the residual type. Moreover, inspired by the work of [22] this paper combines the shifted-inverse iterative method and spectral method to propose an efficient computation scheme. Finally, by using program with Matlab, we implement the numerical experiments and get satisfactory results.The rest of the paper is organized as follows. In Section2, some preliminaries needed in this paper are presented. In Section 3, we discuss a priori and a posterior error estimates for spectral method. In Section 4, we establish an efficient scheme combining the shifted-inverse iteration and spectral method. In Section 5, the numerical experiments on the square domain are reported.


## 2. Preliminaries

We consider the following boundary value problem:(1)-Δu+u=0,inΩ,∂u∂n=f,on∂Ω,
where Ω⊂R2 is a rectangular domain and ∂/∂n is the outward normal derivative on ∂Ω. We denote the Sobolev space with norm ∥·∥m and seminorm |·|m by Hm(Ω). Let ∥·∥0 and ∥·∥0,∂Ω be the norms in the space L2(Ω) and L2(∂Ω), respectively.Throughout this paper,C denotes a generic positive constant independent of the polynomial degrees N, which may not be the same at each occurrence.The variational problem associated with (1) is given by the following.Findu∈H1(Ω), such that
(2)a(u,v)=b(f,v),∀v∈H1(Ω),
where
(3)a(u,v)=∫Ω‍∇u·∇v+∫Ω‍uv,∀u,v∈H1(Ω),b(f,v)=∫∂Ω‍fv,∀f,v∈L2(∂Ω).
Obviously, the bilinear form a(·,·) is continuous and H1(Ω)-elliptic, that is, there exists constants A and B independent of u, v, such that
(4)|a(u,v)|≤A∥u∥1∥v∥1,∀u,v∈H1(Ω),a(v,v)≥B∥v∥12,∀u,v∈H1(Ω).From the Lax-Milgram theorem we know that there exists a unique solution to (2).The following regularity results (see [4]) will be needed in the sequel: for f∈L2(∂Ω), the solution of problem (1) u∈H1+r/2(Ω), and
(5)∥u∥1+r/2≤M1∥f∥0,∂Ω;
furthermore, if f∈H1/2(∂Ω), then u∈H1+r(Ω) and
(6)∥u∥1+r≤M2∥f∥1/2,∂Ω,
where M1, M2 are constants independent of f, and r=1 if Ω is convex and r<π/ω if Ω is concave (with ω being the largest inner angle of Ω).Define(7)SN(Ω)=span{xiyj∣0≤i,j≤N}.The spectral-Galerkin approximation to (2) is the following.FinduN∈SN(Ω), such that
(8)a(uN,vN)=b(f,vN),∀vN∈SN(Ω).Form≥3/2 and each integer N≥1, we define the interpolation operator
(9)IN:Hm(Ω)⟶SN(Ω)
with the LGL interpolation nodes in Ω. We quote from [18] the following interpolation estimates for spectral method at the tensor-product LGL points
(10)∥u-INu∥1≤C(m)N1-m∥u∥m.Letu∈Hm(Ω) and uN be the solution of (2) and (8), respectively. Then we derive from Céa Lemma that
(11)∥uN-u∥1≤C(m)N1-m∥u∥m.
Notice that (11) is also valid for special method with modal basis (see [18, 19]).Using Aubin-Nitsche technique, we deduce the a priori error estimate from (5) and (10)-(11):
(12)∥uN-u∥0,∂Ω≤CN1-m-r/2∥u∥m.


## 3. Spectral Method and Its Error Estimates for the Steklov Eigenvalue Problem

### 3.1. Spectral Method for the Steklov Eigenvalue Problem

We consider the following Steklov eigenvalue problem:(13)-Δu+u=0,inΩ,∂u∂n=λu,on∂Ω,
where Ω⊂R2 is a rectangular domain.The variational problem associated with (13) is given by the following.Findλ∈R, 0≠u∈H1(Ω), such that
(14)a(u,v)=λb(u,v),∀v∈H1(Ω),∥u∥0,∂Ω=1.Let(λ,u) be an eigenpair of (14); then from (6) we know that u∈H1+r(Ω).The spectral-Galerkin approximation to (14) is as the follows.FindλN∈R, 0≠uN∈SN(Ω), such that
(15)a(uN,vN)=λNb(uN,vN),∀vN∈SN(Ω),∥uN∥0,∂Ω=1.LetT:L2(∂Ω)→H1(Ω) and TN:L2(∂Ω)→SN(Ω) be the solution operators defined by
(16)a(Tf,v)=b(f,v),∀v∈H1(Ω),a(TNf,v)=b(f,v),∀vN∈SN(Ω).
Obviously, T:H1(Ω)→H1(Ω) is completely continuous and TN:H1(Ω)→SN(Ω) is a finite ranked operator. Equations (14) and (15) have the following equivalent operator forms, respectively:
(17)Tu=λ-1u,TNuN=λN-1uN.


### 3.2. A Prior Error Estimates of Spectral Method for the Steklov Eigenvalue Problem

Let(λ,u) be an eigenpair of (14), let V(λ) be the space spanned by all eigenfunctions corresponding to the eigenvalue λ, and let VN(λ) be the direct sum of the eigenspace corresponding to the eigenvalues λN of (15) that converge to λ.Denote(18)εN=εN(λ)=supu∈V(λ),∥u∥1=1infv∈SN(Ω)∥u-v∥1,ηa(N)=supf∈H1(Ω),∥f∥1=1infv∈SN(Ω)∥Tf-v∥1,ρN=supf∈L2(∂Ω),∥f∥0,∂Ω=1infv∈SN(Ω)∥Tf-v∥1,
where εN characterizes the degree of the space SN(Ω) approximating V(λ), ηa(N) and ρN characterize the degrees of SN(Ω) approximating the generalized solution Tf of (2).Lemma 1.
LetV(λ)⊂Hm(Ω); then
(19)εN≤CN1-msupu∈V(λ),∥u∥1=1∥u∥m,ηa(N)≤CN-r,ρN≤CN-r/2.Proof.
We deduce from the interpolation estimate (10) that
(20)εN≤supu∈V(λ),∥u∥1=1∥u-INu∥1≤CN1-msupu∈V(λ),∥u∥1=1∥u∥m.
Using (10) and (6) we have
(21)ηa(N)≤supf∈H1(Ω),∥f∥1=1∥Tf-INTf∥1≤CN-rsupf∈H1(Ω),∥f∥1=1∥Tf∥1+r≤CN-r.
From (10) and regularity estimate (5) we get
(22)ρN≤supf∈L2(∂Ω),∥f∥0,∂Ω=1∥Tf-INTf∥1≤CN-r/2supf∈L2(∂Ω),∥f∥0,∂Ω=1∥Tf∥1+r/2,Ω≤CN-r/2.
The proof is completed.Theorem 2.
LetV(λ)⊂Hm(Ω); then
(23)λN-λ≤C{N1-msupu∈V(λ),∥u∥1=1∥u∥m}2;
let uN∈VN(λ) and ∥uN∥1=1; then there exists u∈V(λ), such that
(24)∥uN-u∥1≤CN1-msupu∈V(λ),∥u∥1=1∥u∥m,(25)∥uN-u∥0,∂Ω≤CN1-m-r/2supu∈V(λ),∥u∥1=1∥u∥m;
let u∈V(λ) and ∥u∥1=1; then there exists uN∈VN(λ), such that
(26)∥uN-u∥1≤CN1-msupu∈V(λ),∥u∥1=1∥u∥m.Proof.
From [20, 23] we have that
(27)λN-λ≤CεN2;
let uN∈VN(λ) and ∥uN∥1=1; then there exists u∈V(λ), such that
(28)∥uN-u∥1≤CεN,∥uN-u∥0,∂Ω≤CρNεN;
let u∈V(λ) and ∥u∥1=1; then there exists uN∈VN(λ), such that
(29)∥uN-u∥1≤CεN.
Combining Lemma 1 we obtain the desired results.


### 3.3. A Posterior Error Estimates for the Steklov Eigenvalue Problem

Theorem 2.1 in [24] gave the following results.Lemma 3.
There exists a positive constantC and a linear bounded operator I:L1(Ω)→SN(Ω)⊂L1(Ω), such that for all v∈H1(Ω) and an arbitrary edge e of  ∂Ω, and there holds
(30)∥v-Iv∥0≤C1N∥∇v∥0,∥v-Iv∥0,e≤C1N∥∇v∥0.Analysis of the error indicators requires polynomial inverse estimates in weighted Sobolev spaces. Define on the intervalS=(0,1) the weight function ΦS∶=x(1-x); then there holds the following polynomial inverse estimates.Lemma 4.
Let-1<α<β, σ∈[0,1]; then there exists C1, C2=C(α,β), and C3=Cσ such that for every univariate polynomial πN of degree N, there holds
(31)∫S‍ΦS(x)(πN′(x))2dx≤C1N2∫S‍πN2(x)dx,(32)∫S‍ΦSα(x)πN2(x)dx≤C2N2(β-α)∫S‍ΦSβπN2(x)dx,(33)∫S‍ΦS2σ(x)(πN′(x))2dx≤C3N2(2-σ)∫S‍ΦSσπN2(x)dx.Proof.
See Lemma 2.4 in [21].Consider the domainκ=(0,1)2 and the weight function Φκ∶=dist(x,∂κ). From [21] we have the two-dimensional polynomial inverse estimates.Lemma 5.
Let-1<α<β, σ∈[0,1]; then there exists C1, C2=C(α,β), and C3=Cσ such that for all polynomials πN∈SN(Ω) there hold
(34)∫κ‍Φκ|∇πN|2dxdy≤C1N2∫κ‍|πN|2dxdy,(35)∫κ‍ΦκαπN2dxdy≤C2N2(β-α)∫κ‍ΦκβπN2dxdy,(36)∫κ‍Φκ2σ|∇πN|2dxdy≤C3N2(2-σ)∫κ‍ΦκσπN2dxdy.Proof.
See Lemma 2.5 in [21].Lemma 6.
Letα∈(1/2,1], set e^∶=(0,1)×{0}, and let Φe^ be defined as the same as ΦS, that is, Φe^=x(1-x). Then there exists Cα>0, such that for every ε∈(0,1) and every univariate polynomial πN of degree N, there exists an extension ve^∈H1(κ), such that the following holds:
(37)ve^|e^=πN·Φe^α,ve^|∂κ∖e^=0,(38)∥ve^∥0,κ2≤Cαε∥πNΦe^α/2∥0,e^2,(39)∥∇ve^∥0,κ2≤Cα(εN2(2-α)+ε-1)∥πNΦe^α/2∥0,e^2.Proof.
See Lemma 2.6 in [21].Remark 7.
Because there has an affine transformation from rectangular domain to the unit square, Lemma5 is suited to rectangular domain, and Lemmas 4 and 6 are valid for all edges of rectangular.Define the global error indicatorηα, α∈[0,1], which is expressed as the sum of two terms:
(40)ηα2∶=ηα;B2+ηα;E2,
where the first term ηα;B2 is the weighted internal residual given by
(41)ηα;B2≔1N2∥(ΔuN-uN)ΦΩα/2∥02,
and the second term ηα;E2 is the weighted edge residual given by
(42)ηα;E2∶=1N∥(λNuN-∂uN∂n)Φeα/2∥0,∂Ω2,
where the weight functions ΦΩ and Φe are scaled transformations of the weight functions Φκ and Φe^, respectively.By using the proof of [21], we can get the following results.Theorem 8.
Letα∈[0,1]. Then there exists C>0 independent of N, such that
(43)∥u-uN∥12≤C∥λu-λNuN∥0,∂Ω2+CN2α(ηα;B2+ηα;E2).Proof.
We denotew∶=u-uN-I(u-uN) where I is an interpolation operator satisfying Lemma 3. We derive from H1(Ω)-ellipticity of a(·,·) that
(44)C∥u-uN∥12≤a(u-uN,u-uN)=a(u-uN,w)+a(u-uN,I(u-uN))=∫∂Ω‍λuw-a(uN,w)+∫∂Ω‍(λu-λNuN)I(u-uN).
Since
(45)a(uN,w)=∫Ω‍∇uN·∇w+∫Ω‍uNw=-∫Ω‍ΔuNw+∫∂Ω‍∂uN∂nw+∫Ω‍uNw,
thus
(46)C∥u-uN∥12≤∫∂Ω‍(λu-λNuN)w+∫∂Ω‍(λNuN-∂uN∂n)w+∫Ω‍(ΔuN-uN)w+∫∂Ω‍(λu-λNuN)I(u-uN)=∫∂Ω‍(λu-λNuN)(u-uN)+∫∂Ω‍(λNuN-∂uN∂n)w+∫Ω‍(ΔuN-uN)w.
By using Cauchy-Schwartz inequality, the trace theorem, and Lemma 3, we get
(47)∫∂Ω‍(λu-λNuN)(u-uN)≤∥λu-λNuN∥0,∂Ω∥u-uN∥1,∫∂Ω‍(λNuN-∂uN∂n)w≤C∥λNuN-∂uN∂n∥0,∂Ω1N∥∇(u-uN)∥0,∫Ω‍(ΔuN-uN)w≤C∥ΔuN-uN∥01N∥∇(u-uN)∥0.
Plugging (47) into (46) gives
(48)C∥u-uN∥1≤∥λu-λNuN∥0,∂Ω+C1N∥λNuN-∂uN∂n∥0,∂Ω+C1N∥ΔuN-uN∥0.
Setting α=0, β=α in (32) and (35), respectively, we obtain η0;B≤CNαηα;B and η0;E≤CNαηα;E, and thus
(49)C∥u-uN∥1≤∥λu-λNuN∥0,∂Ω+CNαηα;E+CNαηα;B.
The proof is completed.Lemma 9.
Letα∈[0,1]. Then there exists C>0 independent of N, such that
(50)ηα;B2≤CN2(1-α)∥uN-u∥12.Proof.
DenotevΩ∶=(uN-ΔuN)ΦΩα∈H01(Ω), α∈(0,1]; then
(51)∥vΩΦΩ-α/2∥02=∫Ω‍(uN-ΔuN)vΩ=∫Ω‍uNvΩ+∫Ω‍∇uN·∇vΩ-∫∂Ω‍∂uN∂nvΩ=a(uN-u,vΩ)≤C∥uN-u∥1|vΩ|1.
Now, we consider|vΩ|1. Using (35)-(36), we find for α>1/2 that
(52)|vΩ|12=∥∇((uN-ΔuN)ΦΩα)∥02≤2∫Ω‍ΦΩ2α|∇(uN-ΔuN)|2+2∫Ω‍|∇ΦΩα|2(uN-ΔuN)2≤CN2(2-α)∫Ω‍ΦΩα(uN-ΔuN)2+C∫Ω‍ΦΩ2(α-1)(uN-ΔuN)2≤CN2(2-α)∫Ω‍ΦΩα(uN-ΔuN)2=CN2(2-α)∥vΩΦΩ-α/2∥02,
noticing that (35) is valid since the assumption α>1/2 implies 2α-2>-1, hence we take β=α, α=2(α-1) in (35) and the above third inequality holds. By the definition ∥vΩΦΩ-α/2∥0=Nηα;B, we obtain
(53)ηα;B≤CN1-α∥uN-u∥1.
To get an upper bound of the error indicator ηα;B in the case of 0≤α≤1/2, using (35) we find for β>1/2 that
(54)Nηα;B=∥(uN-ΔuN)ΦΩα/2∥0≤CNβ-α∥(uN-ΔuN)ΦΩβ/2∥0≤CNβ-α·N·N1-β∥uN-u∥1.
Therefore,
(55)ηα;B2≤CN2(1-α)∥uN-u∥12.
We obtain the desired result immediately.Lemma 10.
Letα∈[0,1], ε>0. Then there exists C>0 independent of N, such that
(56)ηα;E2≤CNmax{1-2α+2ε,0}ηα;E2≤×{N∥uN-u∥12+1N∥λNuN-λu∥0,∂Ω2}.Proof.
In order to obtain an upper bound for the error indicatorηα;E, we will use weight function and a proper extension operator. Let e be an edge of Ω; we construct a function we∈H1(Ω) and we|e=(λNuN-∂uN/∂n)Φeα.
Letve=(λNuN-∂uN/∂n)Φeα. We extend ve to ∂Ω by using Lemma 6 where the polynomial πN corresponds to λNuN-(∂uN/∂n), and from (37) we know that ve is zero on ∂Ω∖e. We extend we to ∂Ω such that we=0 on ∂Ω∖e; then
(57)∥(λNuN-∂uN∂n)Φeα/2∥0,e2=∫∂Ω‍(λNuN-∂uN∂n)we=∫∂Ω‍(λNuN-λu)we+a(u-uN,we)+∫Ω‍(uN-ΔuN)we≤∥λNuN-λu∥0,∂Ω∥we∥0,∂Ω+C∥u-uN∥1∥we∥1+∥uN-ΔuN∥0∥we∥0.
Now we estimate ∥we∥0,∂Ω, ∥we∥0, ∥we∥1 in turn. First,
(58)∥we∥0,∂Ω2≤C∥(λNuN-∂uN∂n)Φeα/2∥0,e2.
In the case of α∈(1/2,1] we have
(59)∥we∥02≤Cε∥(λNuN-∂uN∂n)Φeα/2∥0,e2,∥we∥12≤C(εN2(2-α)+ε-1)∥(λNuN-∂uN∂n)Φeα/2∥0,e2.
Combining (58) and (59), we get
(60)∥(λNuN-∂uN∂n)Φeα/2∥0,e≤C∥λNuN-λu∥0,∂Ω+C[εN2(2-α)+ε-1]1/2×∥u-uN∥1+Cε1/2∥uN-ΔuN∥0.
By the definition of ηα;E2 and setting α=0 in Lemma 9, we have that
(61)ηα;E2≤C1N∥λNuN-λu∥0,∂Ω2ηα;E2≤+C{1N[εN2(2-α)+ε-1]+N3ε}∥u-uN∥12.
Setting ε=1/N2 yields the assertion for α>1/2 in (61).
For the caseα∈[0,1/2], setting β=1/2+ε and using (35), we find
(62)ηα;E=1N∥(λNuN-∂uN∂n)Φeα/2∥0,∂Ω≤C1NNβ-α∥(λNuN-∂uN∂n)Φeβ/2∥0,∂Ω=CNβ-αηβ;E;
therefore,
(63)ηα;E2≤CN1+2ε-2αηβ;E2.
From the above discussion we complete the proof.Combining Lemmas9 and 10, we obtain the following theorem.Theorem 11.
Letα∈[0,1]. Then there exists C>0 independent of N, such that
(64)ηα2≤CNmax{1-2α+2ε,0}ηα2≤×{N∥uN-u∥12+1N∥λNuN-λu∥0,∂Ω2}.Remark 12.
The term∥λNuN-λu∥0,∂Ω2 is a term of higher order. In fact, for any v∈H1(Ω), we have
(65)∫∂Ω‍(λNuN-λu)v=∫∂Ω‍λ(uN-u)v+∫∂Ω‍(λN-λ)uNv≤(λ∥uN-u∥0,∂Ω+(λN-λ)∥uN∥0,∂Ω)∥v∥0,∂Ω.
Thus, taking v=λNuN-λu and using the a priori estimates (23) and (25), we conclude that
(66)∥λNuN-λu∥0,∂Ω≤CN1-m-r/2supu∈V(λ),∥u∥1=1∥u∥m∥λNuN-λu∥0,∂Ω≤+C{N1-msupu∈V(λ),∥u∥1=1∥u∥m}2.


## 3.1. Spectral Method for the Steklov Eigenvalue Problem

We consider the following Steklov eigenvalue problem:(13)-Δu+u=0,inΩ,∂u∂n=λu,on∂Ω,
where Ω⊂R2 is a rectangular domain.The variational problem associated with (13) is given by the following.Findλ∈R, 0≠u∈H1(Ω), such that
(14)a(u,v)=λb(u,v),∀v∈H1(Ω),∥u∥0,∂Ω=1.Let(λ,u) be an eigenpair of (14); then from (6) we know that u∈H1+r(Ω).The spectral-Galerkin approximation to (14) is as the follows.FindλN∈R, 0≠uN∈SN(Ω), such that
(15)a(uN,vN)=λNb(uN,vN),∀vN∈SN(Ω),∥uN∥0,∂Ω=1.LetT:L2(∂Ω)→H1(Ω) and TN:L2(∂Ω)→SN(Ω) be the solution operators defined by
(16)a(Tf,v)=b(f,v),∀v∈H1(Ω),a(TNf,v)=b(f,v),∀vN∈SN(Ω).
Obviously, T:H1(Ω)→H1(Ω) is completely continuous and TN:H1(Ω)→SN(Ω) is a finite ranked operator. Equations (14) and (15) have the following equivalent operator forms, respectively:
(17)Tu=λ-1u,TNuN=λN-1uN.


## 3.2. A Prior Error Estimates of Spectral Method for the Steklov Eigenvalue Problem

Let(λ,u) be an eigenpair of (14), let V(λ) be the space spanned by all eigenfunctions corresponding to the eigenvalue λ, and let VN(λ) be the direct sum of the eigenspace corresponding to the eigenvalues λN of (15) that converge to λ.Denote(18)εN=εN(λ)=supu∈V(λ),∥u∥1=1infv∈SN(Ω)∥u-v∥1,ηa(N)=supf∈H1(Ω),∥f∥1=1infv∈SN(Ω)∥Tf-v∥1,ρN=supf∈L2(∂Ω),∥f∥0,∂Ω=1infv∈SN(Ω)∥Tf-v∥1,
where εN characterizes the degree of the space SN(Ω) approximating V(λ), ηa(N) and ρN characterize the degrees of SN(Ω) approximating the generalized solution Tf of (2).Lemma 1.
LetV(λ)⊂Hm(Ω); then
(19)εN≤CN1-msupu∈V(λ),∥u∥1=1∥u∥m,ηa(N)≤CN-r,ρN≤CN-r/2.Proof.
We deduce from the interpolation estimate (10) that
(20)εN≤supu∈V(λ),∥u∥1=1∥u-INu∥1≤CN1-msupu∈V(λ),∥u∥1=1∥u∥m.
Using (10) and (6) we have
(21)ηa(N)≤supf∈H1(Ω),∥f∥1=1∥Tf-INTf∥1≤CN-rsupf∈H1(Ω),∥f∥1=1∥Tf∥1+r≤CN-r.
From (10) and regularity estimate (5) we get
(22)ρN≤supf∈L2(∂Ω),∥f∥0,∂Ω=1∥Tf-INTf∥1≤CN-r/2supf∈L2(∂Ω),∥f∥0,∂Ω=1∥Tf∥1+r/2,Ω≤CN-r/2.
The proof is completed.Theorem 2.
LetV(λ)⊂Hm(Ω); then
(23)λN-λ≤C{N1-msupu∈V(λ),∥u∥1=1∥u∥m}2;
let uN∈VN(λ) and ∥uN∥1=1; then there exists u∈V(λ), such that
(24)∥uN-u∥1≤CN1-msupu∈V(λ),∥u∥1=1∥u∥m,(25)∥uN-u∥0,∂Ω≤CN1-m-r/2supu∈V(λ),∥u∥1=1∥u∥m;
let u∈V(λ) and ∥u∥1=1; then there exists uN∈VN(λ), such that
(26)∥uN-u∥1≤CN1-msupu∈V(λ),∥u∥1=1∥u∥m.Proof.
From [20, 23] we have that
(27)λN-λ≤CεN2;
let uN∈VN(λ) and ∥uN∥1=1; then there exists u∈V(λ), such that
(28)∥uN-u∥1≤CεN,∥uN-u∥0,∂Ω≤CρNεN;
let u∈V(λ) and ∥u∥1=1; then there exists uN∈VN(λ), such that
(29)∥uN-u∥1≤CεN.
Combining Lemma 1 we obtain the desired results.


## 3.3. A Posterior Error Estimates for the Steklov Eigenvalue Problem

Theorem 2.1 in [24] gave the following results.Lemma 3.
There exists a positive constantC and a linear bounded operator I:L1(Ω)→SN(Ω)⊂L1(Ω), such that for all v∈H1(Ω) and an arbitrary edge e of  ∂Ω, and there holds
(30)∥v-Iv∥0≤C1N∥∇v∥0,∥v-Iv∥0,e≤C1N∥∇v∥0.Analysis of the error indicators requires polynomial inverse estimates in weighted Sobolev spaces. Define on the intervalS=(0,1) the weight function ΦS∶=x(1-x); then there holds the following polynomial inverse estimates.Lemma 4.
Let-1<α<β, σ∈[0,1]; then there exists C1, C2=C(α,β), and C3=Cσ such that for every univariate polynomial πN of degree N, there holds
(31)∫S‍ΦS(x)(πN′(x))2dx≤C1N2∫S‍πN2(x)dx,(32)∫S‍ΦSα(x)πN2(x)dx≤C2N2(β-α)∫S‍ΦSβπN2(x)dx,(33)∫S‍ΦS2σ(x)(πN′(x))2dx≤C3N2(2-σ)∫S‍ΦSσπN2(x)dx.Proof.
See Lemma 2.4 in [21].Consider the domainκ=(0,1)2 and the weight function Φκ∶=dist(x,∂κ). From [21] we have the two-dimensional polynomial inverse estimates.Lemma 5.
Let-1<α<β, σ∈[0,1]; then there exists C1, C2=C(α,β), and C3=Cσ such that for all polynomials πN∈SN(Ω) there hold
(34)∫κ‍Φκ|∇πN|2dxdy≤C1N2∫κ‍|πN|2dxdy,(35)∫κ‍ΦκαπN2dxdy≤C2N2(β-α)∫κ‍ΦκβπN2dxdy,(36)∫κ‍Φκ2σ|∇πN|2dxdy≤C3N2(2-σ)∫κ‍ΦκσπN2dxdy.Proof.
See Lemma 2.5 in [21].Lemma 6.
Letα∈(1/2,1], set e^∶=(0,1)×{0}, and let Φe^ be defined as the same as ΦS, that is, Φe^=x(1-x). Then there exists Cα>0, such that for every ε∈(0,1) and every univariate polynomial πN of degree N, there exists an extension ve^∈H1(κ), such that the following holds:
(37)ve^|e^=πN·Φe^α,ve^|∂κ∖e^=0,(38)∥ve^∥0,κ2≤Cαε∥πNΦe^α/2∥0,e^2,(39)∥∇ve^∥0,κ2≤Cα(εN2(2-α)+ε-1)∥πNΦe^α/2∥0,e^2.Proof.
See Lemma 2.6 in [21].Remark 7.
Because there has an affine transformation from rectangular domain to the unit square, Lemma5 is suited to rectangular domain, and Lemmas 4 and 6 are valid for all edges of rectangular.Define the global error indicatorηα, α∈[0,1], which is expressed as the sum of two terms:
(40)ηα2∶=ηα;B2+ηα;E2,
where the first term ηα;B2 is the weighted internal residual given by
(41)ηα;B2≔1N2∥(ΔuN-uN)ΦΩα/2∥02,
and the second term ηα;E2 is the weighted edge residual given by
(42)ηα;E2∶=1N∥(λNuN-∂uN∂n)Φeα/2∥0,∂Ω2,
where the weight functions ΦΩ and Φe are scaled transformations of the weight functions Φκ and Φe^, respectively.By using the proof of [21], we can get the following results.Theorem 8.
Letα∈[0,1]. Then there exists C>0 independent of N, such that
(43)∥u-uN∥12≤C∥λu-λNuN∥0,∂Ω2+CN2α(ηα;B2+ηα;E2).Proof.
We denotew∶=u-uN-I(u-uN) where I is an interpolation operator satisfying Lemma 3. We derive from H1(Ω)-ellipticity of a(·,·) that
(44)C∥u-uN∥12≤a(u-uN,u-uN)=a(u-uN,w)+a(u-uN,I(u-uN))=∫∂Ω‍λuw-a(uN,w)+∫∂Ω‍(λu-λNuN)I(u-uN).
Since
(45)a(uN,w)=∫Ω‍∇uN·∇w+∫Ω‍uNw=-∫Ω‍ΔuNw+∫∂Ω‍∂uN∂nw+∫Ω‍uNw,
thus
(46)C∥u-uN∥12≤∫∂Ω‍(λu-λNuN)w+∫∂Ω‍(λNuN-∂uN∂n)w+∫Ω‍(ΔuN-uN)w+∫∂Ω‍(λu-λNuN)I(u-uN)=∫∂Ω‍(λu-λNuN)(u-uN)+∫∂Ω‍(λNuN-∂uN∂n)w+∫Ω‍(ΔuN-uN)w.
By using Cauchy-Schwartz inequality, the trace theorem, and Lemma 3, we get
(47)∫∂Ω‍(λu-λNuN)(u-uN)≤∥λu-λNuN∥0,∂Ω∥u-uN∥1,∫∂Ω‍(λNuN-∂uN∂n)w≤C∥λNuN-∂uN∂n∥0,∂Ω1N∥∇(u-uN)∥0,∫Ω‍(ΔuN-uN)w≤C∥ΔuN-uN∥01N∥∇(u-uN)∥0.
Plugging (47) into (46) gives
(48)C∥u-uN∥1≤∥λu-λNuN∥0,∂Ω+C1N∥λNuN-∂uN∂n∥0,∂Ω+C1N∥ΔuN-uN∥0.
Setting α=0, β=α in (32) and (35), respectively, we obtain η0;B≤CNαηα;B and η0;E≤CNαηα;E, and thus
(49)C∥u-uN∥1≤∥λu-λNuN∥0,∂Ω+CNαηα;E+CNαηα;B.
The proof is completed.Lemma 9.
Letα∈[0,1]. Then there exists C>0 independent of N, such that
(50)ηα;B2≤CN2(1-α)∥uN-u∥12.Proof.
DenotevΩ∶=(uN-ΔuN)ΦΩα∈H01(Ω), α∈(0,1]; then
(51)∥vΩΦΩ-α/2∥02=∫Ω‍(uN-ΔuN)vΩ=∫Ω‍uNvΩ+∫Ω‍∇uN·∇vΩ-∫∂Ω‍∂uN∂nvΩ=a(uN-u,vΩ)≤C∥uN-u∥1|vΩ|1.
Now, we consider|vΩ|1. Using (35)-(36), we find for α>1/2 that
(52)|vΩ|12=∥∇((uN-ΔuN)ΦΩα)∥02≤2∫Ω‍ΦΩ2α|∇(uN-ΔuN)|2+2∫Ω‍|∇ΦΩα|2(uN-ΔuN)2≤CN2(2-α)∫Ω‍ΦΩα(uN-ΔuN)2+C∫Ω‍ΦΩ2(α-1)(uN-ΔuN)2≤CN2(2-α)∫Ω‍ΦΩα(uN-ΔuN)2=CN2(2-α)∥vΩΦΩ-α/2∥02,
noticing that (35) is valid since the assumption α>1/2 implies 2α-2>-1, hence we take β=α, α=2(α-1) in (35) and the above third inequality holds. By the definition ∥vΩΦΩ-α/2∥0=Nηα;B, we obtain
(53)ηα;B≤CN1-α∥uN-u∥1.
To get an upper bound of the error indicator ηα;B in the case of 0≤α≤1/2, using (35) we find for β>1/2 that
(54)Nηα;B=∥(uN-ΔuN)ΦΩα/2∥0≤CNβ-α∥(uN-ΔuN)ΦΩβ/2∥0≤CNβ-α·N·N1-β∥uN-u∥1.
Therefore,
(55)ηα;B2≤CN2(1-α)∥uN-u∥12.
We obtain the desired result immediately.Lemma 10.
Letα∈[0,1], ε>0. Then there exists C>0 independent of N, such that
(56)ηα;E2≤CNmax{1-2α+2ε,0}ηα;E2≤×{N∥uN-u∥12+1N∥λNuN-λu∥0,∂Ω2}.Proof.
In order to obtain an upper bound for the error indicatorηα;E, we will use weight function and a proper extension operator. Let e be an edge of Ω; we construct a function we∈H1(Ω) and we|e=(λNuN-∂uN/∂n)Φeα.
Letve=(λNuN-∂uN/∂n)Φeα. We extend ve to ∂Ω by using Lemma 6 where the polynomial πN corresponds to λNuN-(∂uN/∂n), and from (37) we know that ve is zero on ∂Ω∖e. We extend we to ∂Ω such that we=0 on ∂Ω∖e; then
(57)∥(λNuN-∂uN∂n)Φeα/2∥0,e2=∫∂Ω‍(λNuN-∂uN∂n)we=∫∂Ω‍(λNuN-λu)we+a(u-uN,we)+∫Ω‍(uN-ΔuN)we≤∥λNuN-λu∥0,∂Ω∥we∥0,∂Ω+C∥u-uN∥1∥we∥1+∥uN-ΔuN∥0∥we∥0.
Now we estimate ∥we∥0,∂Ω, ∥we∥0, ∥we∥1 in turn. First,
(58)∥we∥0,∂Ω2≤C∥(λNuN-∂uN∂n)Φeα/2∥0,e2.
In the case of α∈(1/2,1] we have
(59)∥we∥02≤Cε∥(λNuN-∂uN∂n)Φeα/2∥0,e2,∥we∥12≤C(εN2(2-α)+ε-1)∥(λNuN-∂uN∂n)Φeα/2∥0,e2.
Combining (58) and (59), we get
(60)∥(λNuN-∂uN∂n)Φeα/2∥0,e≤C∥λNuN-λu∥0,∂Ω+C[εN2(2-α)+ε-1]1/2×∥u-uN∥1+Cε1/2∥uN-ΔuN∥0.
By the definition of ηα;E2 and setting α=0 in Lemma 9, we have that
(61)ηα;E2≤C1N∥λNuN-λu∥0,∂Ω2ηα;E2≤+C{1N[εN2(2-α)+ε-1]+N3ε}∥u-uN∥12.
Setting ε=1/N2 yields the assertion for α>1/2 in (61).
For the caseα∈[0,1/2], setting β=1/2+ε and using (35), we find
(62)ηα;E=1N∥(λNuN-∂uN∂n)Φeα/2∥0,∂Ω≤C1NNβ-α∥(λNuN-∂uN∂n)Φeβ/2∥0,∂Ω=CNβ-αηβ;E;
therefore,
(63)ηα;E2≤CN1+2ε-2αηβ;E2.
From the above discussion we complete the proof.Combining Lemmas9 and 10, we obtain the following theorem.Theorem 11.
Letα∈[0,1]. Then there exists C>0 independent of N, such that
(64)ηα2≤CNmax{1-2α+2ε,0}ηα2≤×{N∥uN-u∥12+1N∥λNuN-λu∥0,∂Ω2}.Remark 12.
The term∥λNuN-λu∥0,∂Ω2 is a term of higher order. In fact, for any v∈H1(Ω), we have
(65)∫∂Ω‍(λNuN-λu)v=∫∂Ω‍λ(uN-u)v+∫∂Ω‍(λN-λ)uNv≤(λ∥uN-u∥0,∂Ω+(λN-λ)∥uN∥0,∂Ω)∥v∥0,∂Ω.
Thus, taking v=λNuN-λu and using the a priori estimates (23) and (25), we conclude that
(66)∥λNuN-λu∥0,∂Ω≤CN1-m-r/2supu∈V(λ),∥u∥1=1∥u∥m∥λNuN-λu∥0,∂Ω≤+C{N1-msupu∈V(λ),∥u∥1=1∥u∥m}2.


## 4. Spectral Method Based on the Shifted-Inverse Iteration

High efficient finite element schemes have been researched further in recent years; for example, see [22, 25–31]. Based on these works, this paper combines the shifted-inverse iteration and spectral method to propose an efficient scheme.Scheme 13 (spectral method based on the shifted-inverse iteration).
z
z
zStep 1. Solve (13) in SN1(Ω): Find λN1∈R, uN1∈SN1(Ω), such that
(67)a(uN1,v)=λN1b(uN1,v),∀v∈SN1(Ω).Step 2. Solve a linear system in SN2(Ω)(N2>N1). Find uN2∈SN2(Ω), such that
(68)a(uN2,v)-λN1b(uN2,v)=b(uN1,v),∀v∈SN2(Ω).Step 3. Compute the Rayleigh quotient
(69)λN2=a(uN2,uN2)b(uN2,uN2).Let(λN1,uN1) be the kth eigenpair of (67) then (λN2,uN2) is regarded as the kth approximate eigenpair of (13).Theorem 14.
Suppose thatV(λ)⊂Hm(Ω) and N1 is properly large. Let (λN2,uN2) be the kth approximate eigenpair obtained by Scheme 13; then there exists the kth eigenpair (λ,u) of (13), such that
(70)∥uN2-u∥1≤C(N1-r+3(1-m)+N21-m),|λN2-λ|≤C(N1-2r+6(1-m)+N22(1-m)).Proof.
From the proof of Theorem 4.1 in [22] we can deduce that
(71)∥uN2-u∥1≤C(ηa(N1)εN13+εN2),|λN2-λ|≤C(ηa(N1)2εN16+εN22),
which together with Lemma 1 yields (70).Theorem14 tells us, when N1<N2≤N1(r+3(m-1))/(m-1), the degree of precision (DOP) of λN2 and uN2 obtained by Scheme 13 is O(N22(1-m)) and O(N21-m), respectively; that is to say, (λN2,uN2) have the same DOP with (λN2,uN2).


## 5. Numerical Experiments

Consider the Steklov eigenvalue problem (13) where Ω=[0,1]2.We use spectral method with tensor-product nodal basis at the LGL points (LGL spectral method) to solve (13). We construct two-dimensional basis functions by using the tensor product of one-dimensional basis functions. More precisely, let {hk}k=1N+1 be the Lagrange basis polynomials associated with the LGL points; then SN(Ω)=span{ϕi:ϕi=hi1(x)hj1(y),i1,j1=1,2,…,N+1,i=1,2,…,(N+1)2}, and let uN=∑iuiϕi; then (15) is equivalent to the following linear system:
(72)AU=λNBU,
where A is called stiffness matrix, B is called mass matrix, and
(73)U=(u1,u2,…,u(N+1)2)T,A=(ai,j),ai,j=∫Ω‍∇ϕi·∇ϕj+∫Ω‍ϕiϕj,1≤i,j≤(N+1)2,B=(bi,j),bi,j=∫∂Ω‍ϕiϕj,1≤i,j≤(N+1)2.
We program with Matlab to solve (72), and the numerical results are presented in Table 1.Table 1
The eigenvalues of LGL spectral method.
N
DOF
λ
1
λ
2
λ
4
3
16
0.24007936507936
1.49252393896074
2.08264705882353
4
25
0.24007908542791
1.49230719085030
2.08264705882353
5
36
0.24007908542792
1.49230316584831
2.08264705403320
6
49
0.24007908542722
1.49230313468790
2.08264705403320
7
64
0.24007908542723
1.49230313453462
2.08264705403320
8
81
0.24007908542723
1.49230313453360
2.08264705403320
9
100
0.24007908542722
1.49230313453358
2.08264705403318
10
121
0.24007908542721
1.49230313453357
2.08264705403318
11
144
0.24007908542722
1.49230313453359
2.08264705403321
12
169
0.24007908542721
1.49230313453357
2.08264705403318
13
196
0.24007908542723
1.49230313453360
2.08264705403323
14
225
0.24007908542722
1.49230313453357
2.08264705403318
15
256
0.24007908542724
1.49230313453359
2.08264705403320
16
289
0.24007908542723
1.49230313453359
2.08264705403319
17
324
0.24007908542728
1.49230313453361
2.08264705403321
18
361
0.24007908542725
1.49230313453359
2.08264705403319
19
400
0.24007908542721
1.49230313453356
2.08264705403317
20
441
0.24007908542718
1.49230313453349
2.08264705403311
21
484
0.24007908542705
1.49230313453341
2.08264705403305
22
529
0.24007908542721
1.49230313453356
2.08264705403317
23
576
0.24007908542727
1.49230313453364
2.08264705403326
24
625
0.24007908542700
1.49230313453339
2.08264705403301
25
676
0.24007908542723
1.49230313453356
2.08264705403317
27
784
0.24007908542747
1.49230313453375
2.08264705403338
29
900
0.24007908542735
1.49230313453373
2.08264705403335
30
961
0.24007908542716
1.49230313453351
2.08264705403311
31
1024
0.24007908542735
1.49230313453364
2.08264705403327
34
1225
0.24007908542735
1.49230313453361
2.08264705403323In the tables in this section, DOF denotes the degree of freedom (number of nodes).From Table1 we can see that
(74)λ1≈0.240079085427,λ2=λ3≈1.492303134533,λ4≈2.082647054033.We also use the bilinear finite element to solve (13), and the numerical results are listed in Table 2.Table 2
The eigenvalues of the bilinear finite element.
h
DOF
λ
1
λ
2
λ
4
1/8
81
0.24022930645534
1.49577435506528
2.08269976804992
1/16
289
0.24011662732054
1.49317640860702
2.08266035463158
1/32
1089
0.24008847006616
1.49252179566534
2.08265038681122
1/64
4225
0.24008143153493
1.49235782124229
2.08264788770449
1/128
16641
0.24007967195002
1.49231680754965
2.08264726248054
1/150
22801
0.24007951252139
1.49231309100913
2.08264720582238
1/200
40401
0.24007932566556
1.49230873510836
2.08264713941521
1/256
66049
0.24007923205310
1.49230655286905
2.08264710614529
1/300
90601
0.24007919219220
1.49230562369281
2.08264709197866
1/400
160801
0.24007914547718
1.49230453468660
2.08264707537565
1/512
263169
0.24007912210763
1.49230398913452
2.08264706706984From Table2 we see that
(75)λ1≈0.2400791,λ2=λ3≈1.49230,λ4≈2.0826470.We use Scheme13 to solve (13), and let λiN2 be the ith approximate eigenvalue (i=1,2,4), and the numerical results are presented in Table 3.Table 3
The eigenvalues of shifted-inverse iterative method.
N
1
N
2
λ
1
N
2
λ
2
N
2
λ
4
N
2
4
5
0.24007908542792
1.49230316584831
2.08264705403320
4
6
0.24007908542722
1.49230313468790
2.08264705403320
4
7
0.24007908542722
1.49230313453463
2.08264705403320
4
8
0.24007908542723
1.49230313453360
2.08264705403320From Table3 we find the eigenvalues obtained by Scheme 13 have the same accuracy with λi,N2 which is computed directly in SN2(Ω) by LGL spectral method.Since we do not know the exact eigenvalues of the Steklov eigenvalue problem, many research studies reported the approximate eigenvalues; for example, see [15, 22]. This paper uses the LGL spectral method to solve (13), and in Table 1 we give
(76)λ1≈0.240079085427,λ2=λ3≈1.492303134533,λ4≈2.082647054033.
We think these approximate eigenvalues have high accuracy.



---

*Source: 650530-2013-10-12.xml*