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

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.

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 Section 2, 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: −Δ𝑢+𝑢=0,inΩ,𝜕𝑢𝜕𝑛=𝑓,on𝜕Ω, (1) where Ω⊂𝑅² is a rectangular domain and 𝜕/𝜕𝑛 is the outward normal derivative on 𝜕Ω. We denote the Sobolev space with norm ‖⋅‖_𝑚 and seminorm |⋅|_𝑚 by 𝐻^𝑚(Ω). Let ‖⋅‖₀ and ‖⋅‖_0,𝜕Ω be the norms in the space 𝐿²(Ω) and 𝐿²(𝜕Ω), respectively.

Throughout this paper, 𝐶 denotes a generic positive constant independent of the polynomial degrees 𝑁, which may not be the same at each occurrence.

The variational problem associated with (1) is given by the following.

Find 𝑢∈𝐻¹(Ω), such that 𝑎(𝑢,𝑣)=𝑏𝑓,𝑣,∀𝑣∈𝐻¹(Ω), (2) where 𝑎(𝑢,𝑣)=_Ω∇𝑢⋅∇𝑣+_Ω𝑢𝑣,∀𝑢,𝑣∈𝐻¹(Ω),𝑏𝑓,𝑣=_𝜕Ω𝑓𝑣,∀𝑓,𝑣∈𝐿²(𝜕Ω). (3) Obviously, the bilinear form 𝑎(⋅,⋅) is continuous and 𝐻¹(Ω)-elliptic, that is, there exists constants 𝐴 and 𝐵 independent of 𝑢, 𝑣, such that |𝑎(𝑢,𝑣)|≤𝐴‖𝑢‖₁‖𝑣‖₁,∀𝑢,𝑣∈𝐻¹(Ω),𝑎(𝑣,𝑣)≥𝐵‖𝑣‖²₁,∀𝑢,𝑣∈𝐻¹(Ω). (4)

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 𝑓∈𝐿²(𝜕Ω), the solution of problem (1) 𝑢∈𝐻^1+𝑟/2(Ω), and ‖𝑢‖_1+𝑟/2≤𝑀₁𝑓_0,𝜕Ω; (5) furthermore, if 𝑓∈𝐻^1/2(𝜕Ω), then 𝑢∈𝐻^1+𝑟(Ω) and ‖𝑢‖_1+𝑟≤𝑀₂𝑓_1/2,𝜕Ω, (6) where 𝑀₁, 𝑀₂ are constants independent of 𝑓, and 𝑟=1 if Ω is convex and 𝑟<𝜋/𝜔 if Ω is concave (with 𝜔 being the largest inner angle of Ω).

Define 𝑆_𝑁(Ω)=span𝑥^𝑖𝑦^𝑗∣0≤𝑖,𝑗≤𝑁. (7)

The spectral-Galerkin approximation to (2) is the following.

Find 𝑢_𝑁∈𝑆_𝑁(Ω), such that 𝑎𝑢_𝑁,𝑣_𝑁=𝑏𝑓,𝑣_𝑁,∀𝑣_𝑁∈𝑆_𝑁(Ω). (8)

For 𝑚≥3/2 and each integer 𝑁≥1, we define the interpolation operator 𝐼_𝑁∶𝐻^𝑚(Ω)⟶𝑆_𝑁(Ω) (9) with the LGL interpolation nodes in Ω. We quote from [18] the following interpolation estimates for spectral method at the tensor-product LGL points 𝑢−𝐼_𝑁𝑢₁≤𝐶(𝑚)𝑁^1−𝑚‖𝑢‖_𝑚. (10)

Let 𝑢∈𝐻^𝑚(Ω) and 𝑢_𝑁 be the solution of (2) and (8), respectively. Then we derive from Céa Lemma that 𝑢_𝑁−𝑢₁≤𝐶(𝑚)𝑁^1−𝑚‖𝑢‖_𝑚. (11) 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): 𝑢_𝑁−𝑢_0,𝜕Ω≤𝐶𝑁^{1−𝑚−𝑟/2}‖𝑢‖_𝑚. (12)

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: −Δ𝑢+𝑢=0,inΩ,𝜕𝑢𝜕𝑛=𝜆𝑢,on𝜕Ω, (13) where Ω⊂𝑅² is a rectangular domain.

The variational problem associated with (13) is given by the following.

Find 𝜆∈𝑅, 0≠𝑢∈𝐻¹(Ω), such that 𝑎(𝑢,𝑣)=𝜆𝑏(𝑢,𝑣),∀𝑣∈𝐻¹(Ω),‖𝑢‖_0,𝜕Ω=1. (14)

Let (𝜆,𝑢) be an eigenpair of (14); then from (6) we know that 𝑢∈𝐻^1+𝑟(Ω).

The spectral-Galerkin approximation to (14) is as the follows.

Find 𝜆_𝑁∈𝑅, 0≠𝑢_𝑁∈𝑆_𝑁(Ω), such that 𝑎𝑢_𝑁,𝑣_𝑁=𝜆_𝑁𝑏𝑢_𝑁,𝑣_𝑁,∀𝑣_𝑁∈𝑆_𝑁(Ω),𝑢_𝑁0,𝜕Ω=1. (15)

Let 𝑇∶𝐿²(𝜕Ω)→𝐻¹(Ω) and 𝑇_𝑁∶𝐿²(𝜕Ω)→𝑆_𝑁(Ω) be the solution operators defined by 𝑎𝑇𝑓,𝑣=𝑏𝑓,𝑣,∀𝑣∈𝐻¹(Ω),𝑎𝑇_𝑁𝑓,𝑣=𝑏𝑓,𝑣,∀𝑣_𝑁∈𝑆_𝑁(Ω). (16) Obviously, 𝑇∶𝐻¹(Ω)→𝐻¹(Ω) is completely continuous and 𝑇_𝑁∶𝐻¹(Ω)→𝑆_𝑁(Ω) is a finite ranked operator. Equations (14) and (15) have the following equivalent operator forms, respectively: 𝑇𝑢=𝜆⁻¹𝑇𝑢,_𝑁𝑢_𝑁=𝜆_𝑁⁻¹𝑢_𝑁. (17)

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

Let (𝜆,𝑢) be an eigenpair of (14), let 𝑉(𝜆) be the space spanned by all eigenfunctions corresponding to the eigenvalue 𝜆, and let 𝑉_𝑁(𝜆) be the direct sum of the eigenspace corresponding to the eigenvalues 𝜆_𝑁 of (15) that converge to 𝜆.

Denote 𝜀_𝑁=𝜀_𝑁(𝜆)=sup_{‖𝑢‖𝑢∈𝑉(𝜆),1=1}inf_{𝑣∈𝑆𝑁(Ω)}‖𝑢−𝑣‖₁,𝜂_𝑎(𝑁)=sup_{𝑓∈𝐻1𝑓(Ω),1=1}inf_{𝑣∈𝑆𝑁(Ω)}𝑇𝑓−𝑣₁,𝜌_𝑁=sup_{𝑓∈𝐿2𝑓(𝜕Ω),0,𝜕Ω=1}inf_{𝑣∈𝑆𝑁(Ω)}𝑇𝑓−𝑣₁, (18) where 𝜀_𝑁 characterizes the degree of the space 𝑆_𝑁(Ω) approximating 𝑉(𝜆), 𝜂_𝑎(𝑁) and 𝜌_𝑁 characterize the degrees of 𝑆_𝑁(Ω) approximating the generalized solution 𝑇𝑓 of (2).

Lemma 1. Let 𝑉(𝜆)⊂𝐻^𝑚(Ω); then 𝜀_𝑁≤𝐶𝑁^1−𝑚sup_{(𝜆)𝑢∈𝑉,‖𝑢‖1=1}‖𝑢‖_𝑚,𝜂_𝑎(𝑁)≤𝐶𝑁^−𝑟,𝜌_𝑁≤𝐶𝑁^−𝑟/2. (19)

Proof. We deduce from the interpolation estimate (10) that 𝜀_𝑁≤sup_{(𝜆),‖𝑢‖𝑢∈𝑉1=1}𝑢−𝐼_𝑁𝑢₁≤𝐶𝑁^1−𝑚sup_{(𝜆),‖𝑢‖𝑢∈𝑉1=1}‖𝑢‖_𝑚. (20) Using (10) and (6) we have 𝜂_𝑎(𝑁)≤sup_{𝑓∈𝐻1(Ω),𝑓1=1}𝑇𝑓−𝐼_𝑁𝑇𝑓₁≤𝐶𝑁^−𝑟sup_{𝑓∈𝐻1(Ω),𝑓1=1}𝑇𝑓_1+𝑟≤𝐶𝑁^−𝑟. (21) From (10) and regularity estimate (5) we get 𝜌_𝑁≤sup_{𝑓∈𝐿2(),𝑓𝜕Ω0,𝜕Ω=1}𝑇𝑓−𝐼_𝑁𝑇𝑓₁≤𝐶𝑁^−𝑟/2sup_{𝑓∈𝐿2(),𝑓𝜕Ω0,𝜕Ω=1}𝑇𝑓_1+𝑟/2,Ω≤𝐶𝑁^−𝑟/2. (22) The proof is completed.

Theorem 2. Let 𝑉(𝜆)⊂𝐻^𝑚(Ω); then 𝜆_𝑁−𝜆≤𝐶𝑁^1−𝑚sup_{(𝜆),‖𝑢‖𝑢∈𝑉1=1}‖𝑢‖_𝑚²; (23) let 𝑢_𝑁∈𝑉_𝑁(𝜆) and ‖𝑢_𝑁‖₁=1; then there exists 𝑢∈𝑉(𝜆), such that 𝑢_𝑁−𝑢₁≤𝐶𝑁^1−𝑚sup_{(𝜆),‖𝑢‖𝑢∈𝑉1=1}‖𝑢‖_𝑚,𝑢_𝑁−𝑢_0,𝜕Ω≤𝐶𝑁^{1−𝑚−𝑟/2}sup_{(𝜆),‖𝑢‖𝑢∈𝑉1=1}‖𝑢‖_𝑚; (24)(25) let 𝑢∈𝑉(𝜆) and ‖𝑢‖₁=1; then there exists 𝑢_𝑁∈𝑉_𝑁(𝜆), such that 𝑢_𝑁−𝑢₁≤𝐶𝑁^1−𝑚sup_{‖𝑢‖𝑢∈𝑉(𝜆),1=1}‖𝑢‖_𝑚. (26)

Proof. From [20, 23] we have that 𝜆_𝑁−𝜆≤𝐶𝜀²_𝑁; (27) let 𝑢_𝑁∈𝑉_𝑁(𝜆) and ‖𝑢_𝑁‖₁=1; then there exists 𝑢∈𝑉(𝜆), such that 𝑢_𝑁−𝑢₁≤𝐶𝜀_𝑁,𝑢_𝑁−𝑢_0,𝜕Ω≤𝐶𝜌_𝑁𝜀_𝑁; (28) let 𝑢∈𝑉(𝜆) and ‖𝑢‖₁=1; then there exists 𝑢_𝑁∈𝑉_𝑁(𝜆), such that 𝑢_𝑁−𝑢₁≤𝐶𝜀_𝑁. (29) 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 constant 𝐶 and a linear bounded operator 𝐼∶𝐿¹(Ω)→𝑆_𝑁(Ω)⊂𝐿¹(Ω), such that for all 𝑣∈𝐻¹(Ω) and an arbitrary edge 𝑒 of  𝜕Ω, and there holds ‖𝑣−𝐼𝑣‖₀1≤𝐶𝑁‖∇𝑣‖₀,‖𝑣−𝐼𝑣‖_0,𝑒1≤𝐶√𝑁‖∇𝑣‖₀. (30)

Analysis of the error indicators requires polynomial inverse estimates in weighted Sobolev spaces. Define on the interval 𝑆=(0,1) the weight function Φ_𝑆∶=𝑥(1−𝑥); then there holds the following polynomial inverse estimates.

Lemma 4. Let −1<𝛼<𝛽, 𝜎∈[0,1]; then there exists 𝐶₁, 𝐶₂=𝐶(𝛼,𝛽), and 𝐶₃=𝐶_𝜎 such that for every univariate polynomial 𝜋_𝑁 of degree 𝑁, there holds _𝑆Φ_𝑆(𝑥)𝜋^′_𝑁(𝑥)²𝑑𝑥≤𝐶₁𝑁²_𝑆𝜋²_𝑁(𝑥)𝑑x,_𝑆Φ^𝛼_𝑆(𝑥)𝜋²_𝑁(𝑥)𝑑𝑥≤𝐶₂𝑁^{2(𝛽−𝛼)}_𝑆Φ^𝛽_𝑆𝜋²_𝑁(𝑥)𝑑𝑥,_𝑆Φ_𝑆^2𝜎(𝑥)𝜋^′_𝑁(𝑥)²𝑑𝑥≤𝐶₃𝑁^2(2−𝜎)_𝑆Φ^𝜎_𝑆𝜋²_𝑁(𝑥)𝑑𝑥. (31)(32)(33)

Proof. See Lemma 2.4 in [21].

Consider the domain 𝜅=(0,1)² and the weight function Φ_𝜅∶=dist(𝑥,𝜕𝜅). From [21] we have the two-dimensional polynomial inverse estimates.

Lemma 5. Let −1<𝛼<𝛽, 𝜎∈[0,1]; then there exists 𝐶₁, 𝐶₂=𝐶(𝛼,𝛽), and 𝐶₃=𝐶_𝜎 such that for all polynomials 𝜋_𝑁∈𝑆_𝑁(Ω) there hold _𝜅Φ_𝜅∇𝜋_𝑁²𝑑𝑥𝑑𝑦≤𝐶₁𝑁²_𝜅𝜋_𝑁²𝑑𝑥𝑑𝑦,_𝜅Φ^𝛼_𝜅𝜋²_𝑁𝑑𝑥𝑑𝑦≤𝐶₂𝑁^{2(𝛽−𝛼)}_𝜅Φ^𝛽_𝜅𝜋²_𝑁𝑑𝑥𝑑𝑦,_𝜅Φ_𝜅^2𝜎∇𝜋_𝑁²𝑑𝑥𝑑𝑦≤𝐶₃𝑁^2(2−𝜎)_𝜅Φ^𝜎_𝜅𝜋²_𝑁𝑑𝑥𝑑𝑦. (34)(35)(36)

Proof. See Lemma 2.5 in [21].

Lemma 6. Let 𝛼∈(1/2,1], set 𝑒̂∶=(0,1)×{0}, and let Φ𝑒̂ be defined as the same as Φ_𝑆, that is, Φ𝑒̂=𝑥(1−𝑥). Then there exists 𝐶_𝛼>0, such that for every 𝜀∈(0,1) and every univariate polynomial 𝜋_𝑁 of degree 𝑁, there exists an extension 𝑣𝑒̂∈𝐻¹(𝜅), such that the following holds: 𝑣𝑒̂𝑒̂=𝜋_𝑁⋅Φ𝑒^𝛼,̂𝑣𝑒̂_𝜕𝜅⧵𝑒̂𝑣𝑒̂=0,_20,𝜅≤𝐶_𝛼𝜀𝜋_𝑁Φ̂^𝛼/2𝑒²_0,𝑒̂,∇𝑣𝑒̂²_0,𝜅≤𝐶_𝛼𝜀𝑁^2(2−𝛼)+𝜀⁻¹𝜋_𝑁Φ𝑒̂^𝛼/22_0,𝑒.̂ (37)(38)(39)

Proof. See Lemma 2.6 in [21].

Remark 7. Because there has an affine transformation from rectangular domain to the unit square, Lemma 5 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) where the first term 𝜂²_𝛼;𝐵 is the weighted internal residual given by 𝜂²_𝛼;𝐵1∶=𝑁²Δ𝑢_𝑁−𝑢_𝑁Φ_Ω^𝛼/22₀, (41) and the second term 𝜂²_𝛼;𝐸 is the weighted edge residual given by 𝜂²_𝛼;𝐸1∶=𝑁𝜆_𝑁𝑢_𝑁−𝜕𝑢_𝑁Φ𝜕𝑛^𝑒𝛼/22_0,𝜕Ω, (42) where the weight functions Φ_Ω and Φ_𝑒 are scaled transformations of the weight functions Φ_𝜅 and Φ𝑒̂, respectively.

By using the proof of [21], we can get the following results.

Theorem 8. Let 𝛼∈[0,1]. Then there exists 𝐶>0 independent of 𝑁, such that 𝑢−𝑢_𝑁²₁≤𝐶𝜆𝑢−𝜆_𝑁𝑢_𝑁²_0,𝜕Ω+𝐶𝑁^2𝛼𝜂²_𝛼;𝐵+𝜂²_𝛼;𝐸. (43)

Proof. We denote 𝑤∶=𝑢−𝑢_𝑁−𝐼(𝑢−𝑢_𝑁) where 𝐼 is an interpolation operator satisfying Lemma 3. We derive from 𝐻¹(Ω)-ellipticity of 𝑎(⋅,⋅) that 𝐶𝑢−𝑢_𝑁²₁≤𝑎𝑢−𝑢_𝑁,𝑢−𝑢_𝑁=𝑎𝑢−𝑢_𝑁,𝑤+𝑎𝑢−𝑢_𝑁,𝐼𝑢−𝑢_𝑁=_𝜕Ω𝜆𝑢𝑤−𝑎𝑢_𝑁+,𝑤_𝜕Ω𝜆𝑢−𝜆_𝑁𝑢_𝑁𝐼𝑢−𝑢_𝑁. (44) Since 𝑎𝑢_𝑁,𝑤=_Ω∇𝑢_𝑁⋅∇𝑤+_Ω𝑢_𝑁𝑤=−_ΩΔ𝑢_𝑁𝑤+_𝜕Ω𝜕𝑢_𝑁𝜕𝑛𝑤+_Ω𝑢_𝑁𝑤, (45) thus 𝐶𝑢−𝑢_𝑁²₁≤_𝜕Ω𝜆𝑢−𝜆_𝑁𝑢_𝑁𝑤+_𝜕Ω𝜆_𝑁𝑢_𝑁−𝜕𝑢_𝑁𝑤+𝜕𝑛_ΩΔ𝑢_𝑁−𝑢_𝑁𝑤+_𝜕Ω𝜆𝑢−𝜆_𝑁𝑢_𝑁𝐼𝑢−𝑢_𝑁=_𝜕Ω𝜆𝑢−𝜆_𝑁𝑢_𝑁𝑢−𝑢_𝑁+_𝜕Ω𝜆_𝑁𝑢_𝑁−𝜕𝑢_𝑁𝜕𝑛𝑤+_ΩΔ𝑢_𝑁−𝑢_𝑁𝑤. (46) By using Cauchy-Schwartz inequality, the trace theorem, and Lemma 3, we get _𝜕Ω𝜆𝑢−𝜆_𝑁𝑢_𝑁𝑢−𝑢_𝑁≤𝜆𝑢−𝜆_𝑁𝑢_𝑁0,𝜕Ω𝑢−𝑢_𝑁1,_𝜕Ω𝜆_𝑁𝑢_𝑁−𝜕𝑢_𝑁𝑤𝜆𝜕𝑛≤𝐶_𝑁𝑢_𝑁−𝜕𝑢_𝑁𝜕𝑛_0,𝜕Ω1√𝑁∇𝑢−𝑢_𝑁0,_ΩΔ𝑢_𝑁−𝑢_𝑁𝑤≤𝐶Δ𝑢_𝑁−𝑢_𝑁01𝑁∇𝑢−𝑢_𝑁0. (47) Plugging (47) into (46) gives 𝐶𝑢−𝑢_𝑁1≤𝜆𝑢−𝜆_𝑁𝑢_𝑁0,𝜕Ω1+𝐶√𝑁𝜆_𝑁𝑢_𝑁−𝜕𝑢_𝑁𝜕𝑛_0,𝜕Ω1+𝐶𝑁Δ𝑢_𝑁−𝑢_𝑁0. (48) Setting 𝛼=0, 𝛽=𝛼 in (32) and (35), respectively, we obtain 𝜂_0;𝐵≤𝐶𝑁^𝛼𝜂_𝛼;𝐵 and 𝜂_0;𝐸≤𝐶𝑁^𝛼𝜂_𝛼;𝐸, and thus 𝐶𝑢−𝑢_𝑁1≤𝜆𝑢−𝜆_𝑁𝑢_𝑁0,𝜕Ω+𝐶𝑁^𝛼𝜂_𝛼;𝐸+𝐶𝑁^𝛼𝜂_𝛼;𝐵. (49) The proof is completed.

Lemma 9. Let 𝛼∈[0,1]. Then there exists 𝐶>0 independent of 𝑁, such that 𝜂²_𝛼;𝐵≤𝐶𝑁^2(1−𝛼)𝑢_𝑁−𝑢²₁. (50)

Proof. Denote 𝑣_Ω∶=(𝑢_𝑁−Δ𝑢_𝑁)Φ^𝛼_Ω∈𝐻¹₀(Ω), 𝛼∈(0,1]; then 𝑣_ΩΦ_Ω^−𝛼/22₀=_Ω𝑢_𝑁−Δ𝑢_𝑁𝑣_Ω=_Ω𝑢_𝑁𝑣_Ω+_Ω∇𝑢_𝑁⋅∇𝑣_Ω−_𝜕Ω𝜕𝑢_𝑁𝑣𝜕𝑛^Ω=𝑎𝑢_𝑁−𝑢,𝑣_Ω≤𝐶𝑢_𝑁−𝑢₁𝑣_Ω1. (51)
Now, we consider |𝑣_Ω|₁. Using (35)-(36), we find for 𝛼>1/2 that 𝑣_Ω²₁=∇𝑢_𝑁−Δ𝑢_𝑁Φ^𝛼_Ω²₀≤2_ΩΦ_Ω^2𝛼∇𝑢_𝑁−Δ𝑢_𝑁²+2_Ω∇Φ^𝛼_Ω²𝑢_𝑁−Δ𝑢_𝑁²≤𝐶𝑁^2(2−𝛼)_ΩΦ^𝛼_Ω𝑢_𝑁−Δ𝑢_𝑁²+𝐶_ΩΦ_Ω^2(𝛼−1)𝑢_𝑁−Δ𝑢_𝑁²≤𝐶𝑁^2(2−𝛼)_ΩΦ^𝛼_Ω𝑢_𝑁−Δ𝑢_𝑁²=𝐶𝑁^2(2−𝛼)𝑣_ΩΦ_Ω^−𝛼/22₀, (52) 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 ‖𝑣_ΩΦ_Ω^−𝛼/2‖₀=𝑁𝜂_𝛼;𝐵, we obtain 𝜂_𝛼;𝐵≤𝐶𝑁^1−𝛼𝑢_𝑁−𝑢₁. (53) To get an upper bound of the error indicator 𝜂_𝛼;𝐵 in the case of 0≤𝛼≤1/2, using (35) we find for 𝛽>1/2 that 𝑁𝜂_𝛼;𝐵=𝑢_𝑁−Δ𝑢_𝑁Φ_Ω^𝛼/2₀≤𝐶𝑁^𝛽−𝛼𝑢_𝑁−Δ𝑢_𝑁Φ_Ω^𝛽/2₀≤𝐶𝑁^𝛽−𝛼⋅𝑁⋅𝑁^1−𝛽𝑢_𝑁−𝑢₁. (54) Therefore, 𝜂²_𝛼;𝐵≤𝐶𝑁^2(1−𝛼)𝑢_𝑁−𝑢²₁. (55) We obtain the desired result immediately.

Lemma 10. Let 𝛼∈[0,1], 𝜀>0. Then there exists 𝐶>0 independent of 𝑁, such that 𝜂²_𝛼;𝐸≤𝐶𝑁^{max{1−2𝛼+2𝜀,0}}×𝑁𝑢_𝑁−𝑢²₁+1𝑁𝜆_𝑁𝑢_𝑁−𝜆𝑢²_0,𝜕Ω. (56)

Proof. In order to obtain an upper bound for the error indicator 𝜂_𝛼;𝐸, we will use weight function and a proper extension operator. Let 𝑒 be an edge of Ω; we construct a function 𝑤_𝑒∈𝐻¹(Ω) and 𝑤_𝑒|_𝑒=(𝜆_𝑁𝑢_𝑁−𝜕𝑢_𝑁/𝜕𝑛)Φ^𝛼_𝑒.
Let 𝑣_𝑒=(𝜆_𝑁𝑢_𝑁−𝜕𝑢_𝑁/𝜕𝑛)Φ^𝛼_𝑒. We extend 𝑣_𝑒 to 𝜕Ω by using Lemma 6 where the polynomial 𝜋_𝑁 corresponds to 𝜆_𝑁𝑢_𝑁−(𝜕𝑢_𝑁/𝜕𝑛), and from (37) we know that 𝑣_𝑒 is zero on 𝜕Ω⧵𝑒. We extend 𝑤_𝑒 to 𝜕Ω such that 𝑤_𝑒=0 on 𝜕Ω⧵𝑒; then 𝜆_𝑁𝑢_𝑁−𝜕𝑢_𝑁Φ𝜕𝑛^𝑒𝛼/22_0,𝑒=_𝜕Ω𝜆_𝑁𝑢_𝑁−𝜕𝑢_𝑁𝑤𝜕𝑛^𝑒=_𝜕Ω𝜆_𝑁𝑢_𝑁−𝜆𝑢𝑤_𝑒+𝑎𝑢−𝑢_𝑁,𝑤_𝑒+_Ω𝑢_𝑁−Δ𝑢_𝑁𝑤_𝑒≤𝜆_𝑁𝑢_𝑁−𝜆𝑢_0,𝜕Ω𝑤_𝑒0,𝜕Ω+𝐶𝑢−𝑢_𝑁1𝑤_𝑒1+𝑢_𝑁−Δ𝑢_𝑁0𝑤_𝑒0. (57) Now we estimate ‖𝑤_𝑒‖_0,𝜕Ω, ‖𝑤_𝑒‖₀, ‖𝑤_𝑒‖₁ in turn. First, 𝑤_𝑒²_0,𝜕Ω≤𝐶𝜆_𝑁𝑢_𝑁−𝜕𝑢_𝑁Φ𝜕𝑛^𝑒𝛼/22_0,𝑒. (58) In the case of 𝛼∈(1/2,1] we have 𝑤_𝑒²₀≤𝐶𝜀𝜆_𝑁𝑢_𝑁−𝜕𝑢_𝑁Φ𝜕𝑛^𝑒𝛼/22_0,𝑒,𝑤_𝑒²₁≤𝐶𝜀𝑁^2(2−𝛼)+𝜀⁻¹𝜆_𝑁𝑢_𝑁−𝜕𝑢_𝑁Φ𝜕𝑛^𝑒𝛼/22_0,𝑒. (59) Combining (58) and (59), we get 𝜆_𝑁𝑢_𝑁−𝜕𝑢_𝑁Φ𝜕𝑛^𝑒𝛼/2_0,𝑒≤𝐶𝜆_𝑁𝑢_𝑁−𝜆𝑢_0,𝜕Ω+𝐶𝜀𝑁^2(2−𝛼)+𝜀^−11/2×𝑢−𝑢_𝑁1+𝐶𝜀^1/2𝑢_𝑁−Δ𝑢_𝑁0. (60) By the definition of 𝜂²_𝛼;𝐸 and setting 𝛼=0 in Lemma 9, we have that 𝜂²_𝛼;𝐸1≤𝐶𝑁𝜆_𝑁𝑢_𝑁−𝜆𝑢²_0,𝜕Ω1+𝐶𝑁𝜀𝑁^2(2−𝛼)+𝜀⁻¹+𝑁³𝜀𝑢−𝑢_𝑁²₁. (61) Setting 𝜀=1/𝑁² yields the assertion for 𝛼>1/2 in (61).
For the case 𝛼∈[0,1/2], setting 𝛽=1/2+𝜀 and using (35), we find 𝜂_𝛼;𝐸=1√𝑁𝜆_𝑁𝑢_𝑁−𝜕𝑢_𝑁Φ𝜕𝑛^𝑒𝛼/2_0,𝜕Ω1≤𝐶√𝑁𝑁^𝛽−𝛼𝜆_𝑁𝑢_𝑁−𝜕𝑢_𝑁Φ𝜕𝑛^𝑒𝛽/2_0,𝜕Ω=𝐶𝑁^𝛽−𝛼𝜂_𝛽;𝐸; (62) therefore, 𝜂²_𝛼;𝐸≤𝐶𝑁^{1+2𝜀−2𝛼}𝜂²_𝛽;𝐸. (63) From the above discussion we complete the proof.

Combining Lemmas 9 and 10, we obtain the following theorem.

Theorem 11. Let 𝛼∈[0,1]. Then there exists 𝐶>0 independent of 𝑁, such that 𝜂²_𝛼≤𝐶𝑁^{max{1−2𝛼+2𝜀,0}}×𝑁𝑢_𝑁−𝑢²₁+1𝑁𝜆_𝑁𝑢_𝑁−𝜆𝑢²_0,𝜕Ω. (64)

Remark 12. The term ‖𝜆_𝑁𝑢_𝑁−𝜆𝑢‖²_0,𝜕Ω is a term of higher order. In fact, for any 𝑣∈𝐻¹(Ω), we have _𝜕Ω𝜆_𝑁𝑢_𝑁−𝜆𝑢𝑣=_𝜕Ω𝜆𝑢_𝑁−𝑢𝑣+_𝜕Ω𝜆_𝑁−𝜆𝑢_𝑁𝑣≤𝜆𝑢_𝑁−𝑢_0,𝜕Ω+𝜆_𝑁−𝜆𝑢_𝑁0,𝜕Ω‖𝑣‖_0,𝜕Ω. (65) Thus, taking 𝑣=𝜆_𝑁𝑢_𝑁−𝜆𝑢 and using the a priori estimates (23) and (25), we conclude that 𝜆_𝑁𝑢_𝑁−𝜆𝑢_0,𝜕Ω≤𝐶𝑁^{1−𝑚−𝑟/2}sup_{(𝜆),‖𝑢‖𝑢∈𝑉1=1}‖𝑢‖_𝑚+𝐶𝑁^1−𝑚sup_{(𝜆),‖𝑢‖𝑢∈𝑉1=1}‖𝑢‖_𝑚². (66)

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).

Step 1. Solve (13) in 𝑆_𝑁1(Ω): Find 𝜆_𝑁1∈𝑅, 𝑢_𝑁1∈𝑆_𝑁1(Ω), such that 𝑎𝑢_𝑁1,𝑣=𝜆_𝑁1𝑏𝑢_𝑁1,𝑣,∀𝑣∈𝑆_𝑁1(Ω). (67)

Step 2. Solve a linear system in 𝑆_𝑁2(Ω)(𝑁₂>𝑁₁). Find 𝑢^𝑁2∈𝑆_𝑁2(Ω), such that 𝑎𝑢^𝑁2,𝑣−𝜆_𝑁1𝑏𝑢^𝑁2,𝑣=𝑏𝑢_𝑁1,𝑣,∀𝑣∈𝑆_𝑁2(Ω). (68)

Step 3. Compute the Rayleigh quotient 𝜆^𝑁2=𝑎𝑢^𝑁2,𝑢^𝑁2𝑏𝑢^𝑁2,𝑢^𝑁2. (69)

Let (𝜆_𝑁1,𝑢_𝑁1) be the 𝑘th eigenpair of (67) then (𝜆^𝑁2,𝑢^𝑁2) is regarded as the 𝑘th approximate eigenpair of (13).

Theorem 14. Suppose that 𝑉(𝜆)⊂𝐻^𝑚(Ω) and 𝑁₁ is properly large. Let (𝜆^𝑁2,𝑢^𝑁2) be the 𝑘th approximate eigenpair obtained by Scheme 13; then there exists the 𝑘th eigenpair (𝜆,𝑢) of (13), such that 𝑢^𝑁2−𝑢₁≤𝐶𝑁₁^{−𝑟+3(1−𝑚)}+𝑁₂^1−𝑚,𝜆^𝑁2−𝜆≤𝐶𝑁₁^{−2𝑟+6(1−𝑚)}+𝑁₂^2(1−𝑚). (70)

Proof. From the proof of Theorem 4.1 in [22] we can deduce that 𝑢^𝑁2−𝑢₁≤𝐶𝜂_𝑎𝑁₁𝜀³_𝑁1+𝜀_𝑁2,𝜆^𝑁2−𝜆≤𝐶𝜂_𝑎𝑁₁²𝜀⁶_𝑁1+𝜀²_𝑁2, (71) which together with Lemma 1 yields (70).

Theorem 14 tells us, when 𝑁₁<𝑁₂≤𝑁₁^{(𝑟+3(𝑚−1))/(𝑚−1)}, the degree of precision (DOP) of 𝜆^𝑁2 and 𝑢^𝑁2 obtained by Scheme 13 is 𝑂(𝑁₂^2(1−𝑚)) and 𝑂(𝑁₂^1−𝑚), respectively; that is to say, (𝜆^𝑁2,𝑢^𝑁2) have the same DOP with (𝜆_𝑁2,𝑢_𝑁2).

5. Numerical Experiments

Consider the Steklov eigenvalue problem (13) where Ω=[0,1]².

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 {ℎ_𝑘}^𝑁+1_𝑘=1 be the Lagrange basis polynomials associated with the LGL points; then 𝑆_𝑁(Ω)=span{𝜙_𝑖∶𝜙_𝑖=ℎ_𝑖1(𝑥)ℎ_𝑗1(𝑦),𝑖1,𝑗1=1,2,…,𝑁+1,𝑖=1,2,…,(𝑁+1)²}, and let 𝑢_𝑁=∑_𝑖𝑢_𝑖𝜙_𝑖; then (15) is equivalent to the following linear system: 𝐴𝑈=𝜆_𝑁𝐵𝑈, (72) where 𝐴 is called stiffness matrix, 𝐵 is called mass matrix, and 𝑈=𝑢₁,𝑢₂,…,𝑢_(𝑁+1)^2𝑇,𝐴=𝑎_𝑖,𝑗,𝑎_𝑖,𝑗=_Ω∇𝜙_𝑖⋅∇𝜙_𝑗+_Ω𝜙_𝑖𝜙_𝑗,1≤𝑖,𝑗≤(𝑁+1)²,𝐵=𝑏_𝑖,𝑗,𝑏_𝑖,𝑗=_𝜕Ω𝜙_𝑖𝜙_𝑗,1≤𝑖,𝑗≤(𝑁+1)². (73) We program with Matlab to solve (72), and the numerical results are presented in Table 1.

Table 1 

The eigenvalues of LGL spectral method.

In the tables in this section, DOF denotes the degree of freedom (number of nodes).

From Table 1 we can see that 𝜆₁𝜆≈0.240079085427,₂=𝜆₃≈1.492303134533,𝜆₄≈2.082647054033. (74)

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.

From Table 2 we see that 𝜆₁≈0.2400791,𝜆₂=𝜆₃𝜆≈1.49230,₄≈2.0826470. (75)

We use Scheme 13 to solve (13), and let 𝜆^𝑁2_𝑖 be the 𝑖th approximate eigenvalue (𝑖=1,2,4), and the numerical results are presented in Table 3.

Table 3 

The eigenvalues of shifted-inverse iterative method.

From Table 3 we find the eigenvalues obtained by Scheme 13 have the same accuracy with 𝜆_𝑖,𝑁2 which is computed directly in 𝑆_𝑁2(Ω) 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 𝜆₁𝜆≈0.240079085427,₂=𝜆₃≈1.492303134533,𝜆₄≈2.082647054033. (76) We think these approximate eigenvalues have high accuracy.

Acknowledgment

This work is supported by the National Natural Science Foundation of China (Grant no. 11201093).