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

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.


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 , Schr ödinger equation in 17, 18 , Stokes equations in 19 , and so forth.As for the Steklov eigenvalue

Preliminaries
Let H s Ω and H s ∂Ω denote Sobolev spaces on Ω and ∂Ω with real order s, respectively.The norms in H s Ω and H s ∂Ω are denoted by • s and • s,∂Ω , respectively.
Consider the model problem −Δu u 0 inΩ, ∂u ∂n λu on ∂Ω, 2.1 where Ω ⊂ R 2 is a polygonal domain and ∂u/∂n is the outward normal derivative on ∂Ω.

Mathematical Problems in Engineering
The weak form of 2.1 is given by the following: find λ ∈ R, 0 / u ∈ H 1  From 25 , we know that 2.2 has a countable sequence of positive eigenvalues here each eigenvalue occurs as many times as given by its multiplicity , and the corresponding eigenfunctions u k ∈ H 1 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 S h 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 It is well known that 2.5 has a finite sequence of eigenvalues and their corresponding eigenfunctions are Consider the following source problem 2.7 associated with 2.2 and the approximate source problem 2.8 associated with 2.5 , respectively. Find Several regular estimates of the Steklov eigenvalue problem are presented in the following lemma, which will be used in the sequel.Lemma 2.1.If f ∈ L 2 ∂Ω , then there exists a unique solution u ∈ H 1 r/2 Ω to 2.7 and 2.9 If f ∈ H 1/2 ∂Ω , then there exists a unique solution u ∈ H 1 r Ω to 2.7 and

2.10
If f ∈ L 2 ∂Ω , then there exists a unique solution u ∈ H 1 Ω to 2.7 and where C p 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 H 1 Ω -elliptic, from 2.7 , it is easy to know that 2.11 is valid.
Thus, from 2.7 , we can define the operator A : Similarly, from 2.8 , we define the operator

2.13
It is obvious that A : Analogously, A h is also a self-adjoint operator.Observe that Af and A h f are the exact solution and the finite element solution of 2.7 , respectively, and a Af − A h f, v 0, ∀v ∈ S h ⊂ H 1 Ω .Define the Ritz-Galerkin projection operator P h :

2.15
Therefore, A h f P h Af, ∀f ∈ H 1 Ω , then A h P h A. From Lemma 2.1 and the interpolation error estimate, we have where C I is the interpolation constant.It is clear that A h 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: 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 M h μ k be the direct sum of eigenspaces corresponding to all eigenvalues of A h that converge to

2.18
From Lemma 2.1, we know that M λ k ⊂ H 1 r Ω , then by the interpolation error estimate, we have

2.19
Lemma 2.2.Let λ k,h and λ k be the kth eigenvalue of 2.5 and 2.2 , respectively.Then for any eigenfunction u k,h corresponding to λ k,h , satisfying u k,h a 1, there exists and for any u k ∈ M λ k , there exists u h ∈ M h λ k such that where C i , 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.

2.23
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 , u 0 be an approximation for μ k , u k , where μ 0 is not an eigenvalue of A h , and u 0 ∈ S h with u 0 a 1. Suppose that 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 Lemma 2.4 are satisfied which will be verified in the proof of Theorem 3.1.

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 {π h i } l 1 be a family of regular meshes, h i−1 h i , and let {S h i } l 1 be the conforming finite element spaces defined on {π h i } l 1 , and let π H π h 1 , S H S h 1 , π h π h l , and S h S h l .Scheme 1 multi-scale discretization scheme .
Step 1. Solve 2.2 on the π H : find Step 2. Execute the assignments: Step 3. Solve a linear system on the Step 4. Compute the Rayleigh quotient

3.3
Step 5.If i l, then output λ h l k , u h l k , that is, λ h k , u h k , stop.Else, i ⇐ i 1, and return to Step 3.
where the constants C 4 , C 5 , and C 6 are stated in 3.6 , 3.11 , and 3.29 , respectively.
Proof.We use Lemma 2.4 to complete the proof.First, we will verify that the conditions of Lemma 2.4 are satisfied.Select μ 0 1/λ h l−1 k , and 16 , we know that A h l − A a → 0 h l → 0 , then there exists a constant C 4 independent of h l and l such that Thus, by the assumption, we deduce that

3.7
Note that in any normed space, for any nonzero u, v ∈ S h , there holds Hence, we have 3.9 Using the triangle inequality and 2.22 , we get It follows from 2.20 that λ k,h l → λ k h l → 0 , then by the assumption, we have where C 5 is a constant independent of h l .By 2.13 , we see that Step 3 in Scheme 1 is equivalent to the following: In fact, it is obvious that u h l k obtained by the above two formulae are the same.When h l−1 is small enough, noting that h l h l−1 , from 3.10 and 3.9 , we get From 2.20 , having in mind that 3.17 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 Since ρ is the separation constant, h l−1 is small enough, and h l h l−1 , there holds

3.19
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.21
Let it follows directly from 3.20 that By Lemma 2.2, there exists {u 0 j }

3.25
Combining 3.23 with the above inequality, we have

3.26
It is obvious that there exists u k ∈ M λ k such that a , from 2.23 and 3.8 , we derive

3.28
Substituting 3.28 into 3.26 , we obtain 3.4 .Equation 3.4 indicates that u h l k converges to u k in the sense of norm • a , then from the trace theorem, we know that u h l k converges to u k in the sense of norm

3.29
where the last inequality in the above holds due to the trace theorem, and the constant C 6 is independent of h l .
Set l 2 and denote H h 1 , h h 2 , then we immediately get the following two-grid discretization scheme based on the shifted-inverse power method.
Step 1. Solve 2.2 on a coarse grid π H : find λ k,H ∈ R, u k,H ∈ S H , such that u k,H a 1 and

3.30
Step 2. Solve a linear system on a fine grid π h : find u ∈ S h , such that

3.31
Set u h k u/ u a .
Step 3. Compute the Rayleigh quotient

3.32
Then from Theorem 3.1, we have the following error estimates for Scheme 2.

3.34
Proof.Since λ k,H , u k,H is obtained by Step 1 in Scheme 2, that is, λ k,H , u k,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.

3.35
And from 2.21 , we know that dist u k,H , M λ k ≤ C 2 δ H λ k , then from Theorem 3.2, we get

3.37
From 3.37 and 3.34 , we can see that when we take H O 3 √ h , the resulted solution, λ h k , u h k , can maintain an asymptotically optimal accuracy.Corollary 3.4.Suppose that h i h t i i−1 , t i ∈ 1, 3 , i 2, 3, . . ., sup i t i < 3.By using the triangle linear conforming element with regular meshes, when h 1 , that is, H is properly small, there exists u k ∈ M λ k such that the following error estimates for Scheme 1 hold:

3.39
Proof.We prove the results by using induction.
For l 2, we estimate λ k,H − λ k and dist u k,H , M λ k in 3.33 , respectively.Combining 3.35 and 2.19 yields

3.40
From Lemma 2.2, and 3.8 we have

3.41
Substituting the above estimates into 3.33 , we get

3.42
By the assumption, h H t 2 , t 2 ∈ 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
Suppose that 3.38 , 3.39 hold for l − 1, that is, there holds

3.46
Mathematical Problems in Engineering Substituting 3.45 and 3.46 into 3.4 , and noting 3.44 , we deduce, 3.47 that is, 3.38 holds.The inequality 3.39 follows directly from 3.29 .
From Theorem 3.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 {λ h l k } converges to {λ k } as l → ∞, and the convergence is quick.Hence, we have the reason to believe that

3.48
Note that in 3.48 |λ h l 1 k − λ k | is a small quantity of higher order.So we can use |λ h l k − λ h l−1 k | < ε ε 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 π h 2 derived from π H .
Step 1. Solve 2.2 on the π H : find λ k,H ∈ R, u k,H ∈ S H , such that u k,H a 1 and

3.49
Step 2. Execute the assignments: Step 3. Solve a linear system on the π h i : find u ∈ S h i , such that

3.50
And set u h i k u/ u a .Step 4. Compute the Rayleigh quotient
Step 6. Determine t i 1 and h i 1 h t i 1 i , then refine π h i to get π h i 1 .Set i ⇐ i 1, and return to Step 3.

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 EQ rot 1 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 EQ rot 1 element.
When Ω 0, 1 × 0, 1 the first four exact eigenvalues 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  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 λ h 1 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.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

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.
It is easy to know that a •, • is a symmetric, continuous, and H 1 Ω -elliptic bilinear form onH 1 Ω × H 1 Ω .So,we use a •, • and • a a •, • • 1 as the inner product and norm on H 1 Ω , respectively.We use b •, • and • b b •, • as the inner product and norm on L 2 ∂Ω , respectively.

Table 3 :
Numerical eigenvalues on the square domain Ω 0, 1 × 0, 1 by Scheme 3: set h 1 H 1/8, h 2 Note that the computational complexities of these two methods are almost the same which indicates that our method is also efficient.
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 O 2 √ h .For example, we carry out the computation on the same grids with h