A multigrid discretization scheme based on the shifted inverse iteration for the Steklov eigenvalue problem in inverse scattering

Jiali Xie; Hai Bi

doi:10.1515/math-2022-0607

Article Open Access

A multigrid discretization scheme based on the shifted inverse iteration for the Steklov eigenvalue problem in inverse scattering

Jiali Xie and Hai Bi

Published/Copyright: July 28, 2023

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

$Open Mathematics$

From the journal Open Mathematics Volume 21 Issue 1

Abstract

Numerical methods for computing Steklov eigenvalues have attracted the attention of academia for their important physical background and wide applications. In this article we discuss the multigrid discretization scheme based on the shifted inverse iteration for the Steklov eigenvalue problem in inverse scattering, and give the error estimation of the proposed scheme. In addition, on the basis of the a posteriori error indicator, we design an adaptive multigrid algorithm. Finally, we present numerical examples to show the efficiency of the proposed scheme.

Keywords: Steklov eigenvalue problem in inverse scattering; multigrid discretization; the shifted-inverse iteration; error estimation; adaptive algorithm

MSC 2010: 65N25; 65N30

1 Introduction

Steklov eigenvalues have many applications. They are extensively encountered in surface waves, mechanical oscillators immersed in a viscous fluid, the vibration modes of a structure in contact with an incompressible fluid (see [1–3] and the references therein), etc. In recent years, Steklov eigenvalues have been introduced in scattering problems for the inhomogeneous medium [4]. By combining the eigenvalues measured from far field data with the Steklov eigenvalues obtained by numerical computations, the index of refraction associated with the inhomogeneous medium can be reconstructed using optimization methods and other methods. Due to the wide applications [4–6], numerical methods for computing Steklov eigenvalues have attracted the attention of academia. For instance, for the Steklov eigenvalue problem in inverse scattering, Yang et al. [7] studied the nonconforming Crouzeix-Raviart element and gave the convergence analysis for discrete eigenvalues and eigenfunctions; Liu et al. [8] explored the spectral indicator method; Zhang et al. [9] proposed a multigrid correction scheme and provided the error estimates of eigenvalues and eigenfunctions; Zhang et al. [10] studied the a posteriori error estimates and adaptive algorithm; Meng and Mei [11] discussed the discontinuous Galerkin method; Ren et al. [12] studied spectral method, etc.

The two-grid discretization of finite element method was first introduced by Xu [13,14]. Because the two-grid discretization technique can significantly improve the accuracy of approximations and save computational costs, the two-grid discretization and multigrid discretization developed based on it are widely used to solve eigenvalue problems, such as biharmonic eigenvalue problem [15], semi-linear eigenvalue problem [16], quantum eigenvalue problem [17], Stokes eigenvalue problem [18,19], Maxwell eigenvalue problem [20], 2 m order elliptic eigenvalue problem [21], and so on. Based on the two-grid discretization established in the literature [13,14], another two-grid discretizations based on the shifted inverse iteration has been developed [22]. The two-grid and multigrid discretizations have been based on the shifted inverse iteration can be regarded as a combination of the finite element method and the shifted inverse iteration method for solving matrix eigenvalue problems (see [22–24]). Later, the multigrid method based on shift inverse iteration has been successfully applied to the Laplace eigenvalue problem [25], Maxwell eigenvalue problem [26,27], Stokes eigenvalue problem [28], general self-adjoint eigenvalue problems and integral operator eigenvalue problem [29] and so on.

In this article, we study the multigrid discretization scheme based on the shifted inverse iteration for the Steklov eigenvalue problem in inverse scattering. In the next section, we present some results of conforming element approximations for the problem, and in Section 3, we establish the multigrid scheme. Since the bilinear form in the weak formulation is not coercive, we extend an fundamental lemma in the existing literature to prove a basic tool (Lemma 3.3) to analyze the multigrid scheme and give the error estimation. In addition, on the basis of the a posteriori error indicator, we design an adaptive multigrid algorithm in Section 4. Finally, we present numerical examples to show the efficiency of the multigrid discretization scheme.

2 Preliminaries

Let Ω ⊂ R 2 be a bounded polygonal domain with Lipschitz boundary ∂ Ω . Let H t ( Ω ) denote the usual Sobolev space on Ω with order t , and ‖ ⋅ ‖ t , Ω is the norm on H t ( Ω ) . H 0 ( Ω ) = L 2 ( Ω ) . Let H t ( ∂ Ω ) be the Sobolev space on ∂ Ω with order t equipped with the norm ‖ ⋅ ‖ t , ∂ Ω .

The Steklov eigenvalue problem in inverse scattering states as follows: Find λ ∈ C and a nontrivial function u ∈ H 1 ( Ω ) satisfying

(2.1) Δ u + k 2 n ( x ) u = 0 in Ω , ∂ u ∂ ν + λ u = 0 on ∂ Ω ,

where ν is the unit outward normal to ∂ Ω , k is the wave number, n ( x ) = n 1 ( x ) + i n 2 ( x ) k is the index of refraction, i = − 1 , n 1 ( x ) > 0 , and n 2 ( x ) ≥ 0 are bounded and smooth functions.

Denote ( u , v ) = ∫ Ω u v ¯ d x , ⟨ f , g ⟩ = ∫ ∂ Ω f g ¯ d s , and define a continuous bilinear form

a ( u , v ) ≔ ( ∇ u , ∇ v ) − k 2 ( n u , v ) , ∀ u , v ∈ H 1 ( Ω ) .

For any g ∈ H 1 ( Ω ) , ⟨ f , g ⟩ has a continuous extension to f ∈ H − 1 2 ( ∂ Ω ) such that ⟨ f , g ⟩ is continuous on H − 1 2 ( ∂ Ω ) × H 1 2 ( ∂ Ω ) .

Then the weak formulation of (2.1) is to find ( λ , u ) ∈ C × H 1 ( Ω ) , u ≠ 0 , such that

(2.2) a ( u , v ) = − λ ⟨ u , v ⟩ , ∀ v ∈ H 1 ( Ω ) .

When n ( x ) is a real valued function, it is obvious that a ( ⋅ , ⋅ ) is symmetric, and in this case, the eigenvalues of (2.2) λ ∈ R . In this article, we focus ourselves to the case of n ( x ) being real valued function.

From [8], we know that a ( ⋅ , ⋅ ) satisfies Gärding’s inequality, i.e., there exists constant K < ∞ , α 0 > 0 such that

Re { a ( v , v ) } + K ‖ v ‖ 0 , Ω 2 ≥ α 0 ‖ v ‖ 1 , Ω 2 , ∀ v ∈ H 1 ( Ω ) .

Suppose that K is a sufficiently large positive constant, and define the bilinear form

a ˜ ( u , v ) ≔ a ( u , v ) + K ( u , v ) = ( ∇ u , ∇ v ) − k 2 ( n u , v ) + K ( u , v ) , u , v ∈ H 1 ( Ω ) .

It is ease to verify that a ˜ ( u , v ) is H 1 ( Ω ) -coercive [8].

Let π h = { τ } be a regular mesh of Ω , h ( x ) be the diameter of element τ that containing x , and h = h Ω = max x ∈ Ω h ( x ) be the diameter of π h . Let V h ⊂ C ( Ω ¯ ) be the space of piecewise polynomials of degree m ( m ≥ 1 ) defined on π h and V h B ≔ V h ∣ ∂ Ω be the restriction of V h on ∂ Ω .

Suppose that the finite element spaces in this article satisfy the following condition [30]: There exists m ≥ 1 , such that for ψ ∈ H 1 + s ( Ω ) , it is valid that

inf v ∈ V h ( ∥ h − 1 ( ψ − v ) ∥ 0 , Ω + ‖ ψ − v ‖ 1 , Ω ) ≤ C h s ‖ ψ ‖ 1 + s , Ω , 0 ≤ s ≤ m .

The conforming finite element discretization of (2.2) is to find ( λ h , u h ) ∈ R × V h , u h ≠ 0 such that

(2.3) a ( u h , v h ) = − λ h ⟨ u h , v h ⟩ , ∀ v h ∈ V h .

Consider the following source problem (2.4) associated with the eigenvalue problem (2.2), and the approximate source problem (2.5) associated with (2.3):

Given g ∈ H − 1 2 ( ∂ Ω ) , find ϕ ∈ H 1 ( Ω ) satisfying

(2.4) a ( ϕ , v ) = ⟨ g , v ⟩ , ∀ v ∈ H 1 ( Ω ) ;

Find ϕ h ∈ V h satisfying

(2.5) a ( ϕ h , v ) = ⟨ g , v ⟩ , ∀ v ∈ V h .

Introduce the following Neumann eigenvalue problem:

(2.6) Δ ϕ + k 2 n ( x ) ϕ = 0 in Ω , ∂ ϕ ∂ v = 0 on ∂ Ω .

When k 2 is not an eigenvalue of (2.6), from Fredholm alternative (see Section 5.3 in [31], or Lemma 1 in [7]), we know that for g ∈ H − 1 2 ( ∂ Ω ) , (2.4) admits the unique solution ϕ ∈ H 1 ( Ω ) satisfying

(2.7) ‖ ϕ ‖ 1 , Ω ≤ C ‖ g ‖ − 1 2 , ∂ Ω .

So we can define the solution operator A : H − 1 2 ( ∂ Ω ) → H 1 ( Ω ) by

a ( A g , v ) = ⟨ g , v ⟩ , ∀ v ∈ H 1 ( Ω ) ,

and a Neumann-to-Dirichlet mapping T : H − 1 2 ( ∂ Ω ) → H 1 2 ( ∂ Ω ) satisfying

T g = ( A g ) ∣ ∂ Ω ,

i.e., T is the restriction of A on ∂ Ω .

Then, (2.2) has the following equivalent operator form:

T u = − 1 λ u .

Similarly, from (2.5), we can define the discrete solution operator A h : H − 1 2 ( ∂ Ω ) → V h by

(2.8) a ( A h g , v ) = ⟨ g , v ⟩ , ∀ v ∈ V h ,

and define T h : H − 1 2 ( ∂ Ω ) → V h B satisfying

T h g = ( A h g ) ∣ ∂ Ω .

Thus, (2.3) has the equivalent operator form:

T h u h = − 1 λ h u h .

In this article, we always assume that k 2 is not an eigenvalue of (2.6).

According to [32] or Proposition 4.1 in [33], we have the following regularity estimates as for the source problem (2.4).

Lemma 2.1

Let ϕ be the solution of (2.4). If g ∈ L 2 ( ∂ Ω ) , then ϕ ∈ H 1 + r 2 ( Ω ) and

‖ ϕ ‖ 1 + r 2 , Ω ≤ C R ‖ g ‖ 0 , ∂ Ω .

If g ∈ H 1 2 ( ∂ Ω ) , then ϕ ∈ H 1 + r ( Ω ) and

(2.9) ‖ ϕ ‖ 1 + r , Ω ≤ C R ‖ g ‖ 1 2 , ∂ Ω ,

where C R is the regularity constant, r = 1 when Ω is a convex domain, and when Ω is nonconvex, r < π θ and r can be arbitrarily close to π θ with θ being the largest inner angle of Ω .

Lemma 2.1 ensures that the eigenfunction of (2.2) u ∈ H 1 + r ( Ω ) .

Define the projection operator P h : H 1 ( Ω ) → V h by

(2.10) a ( ϕ − P h ϕ , v ) = 0 , ∀ v ∈ V h .

Then, for any g ∈ H − 1 2 ( ∂ Ω ) , we have

a ( A h g − P h ( A g ) , v ) = a ( A h g − A g + A g − P h ( A g ) , v ) = a ( A h g − A g , v ) + a ( A g − P h ( A g ) , v ) = 0 , ∀ v ∈ V h ,

thus, A h g = P h A g , ∀ g ∈ H − 1 2 ( ∂ Ω ) , and hence, A h = P h A .

Denote ρ h = sup f ∈ L 2 ( ∂ Ω ) , ‖ f ‖ 0 , ∂ Ω = 1 inf v ∈ V h ‖ A f − v ‖ 1 , Ω .

From [8] and the Aubin-Nitsche technique, we have the following results.

Lemma 2.2

Let ϕ be the solution of (2.4). If ϕ ∈ H 1 + t ( Ω ) ( t ≥ r ) , then

(2.11) ‖ ϕ − P h ϕ ‖ 1 , Ω ≤ C inf v ∈ V h ‖ ϕ − v ‖ 1 , Ω ,

(2.12) ‖ ϕ − P h ϕ ‖ 0 , ∂ Ω ≤ C ρ h ‖ ϕ − P h ϕ ‖ 1 , Ω .

Proof

From Theorem 3.1 in [8], we know that (2.11) holds. Next we will prove (2.12).

By (2.10) and the Aubin-Nitsche technique, for any v ∈ V h , we deduce that

∥ ϕ − P h ϕ ∥ 0 , ∂ Ω = sup g ∈ L 2 ( ∂ Ω ) , g ≠ 0 ∣ ⟨ ϕ − P h ϕ , g ⟩ ∣ ‖ g ‖ 0 , ∂ Ω = sup g ∈ L 2 ( ∂ Ω ) , g ≠ 0 ∣ a ( A g , ϕ − P h ϕ ) ∣ ‖ g ‖ 0 , ∂ Ω = sup g ∈ L 2 ( ∂ Ω ) , v ∈ V h ( Ω ) ∣ a ( ϕ − P h ϕ , A g − v ) ∣ ‖ g ‖ 0 , ∂ Ω ≤ C sup g ∈ L 2 ( ∂ Ω ) , v ∈ V h ( Ω ) ∥ ϕ − P h ϕ ∥ 1 , Ω ‖ A g − v ‖ 1 , Ω ‖ g ‖ 0 , ∂ Ω ≤ C ρ h ∥ ϕ − P h ϕ ∥ 1 , Ω .

The proof is completed.□

Lemma 2.3

For any ϕ ∈ H 1 ( Ω ) , there holds

(2.13) ‖ P h ϕ ‖ 1 , Ω ≤ C ‖ ϕ ‖ 1 , Ω .

Proof

See Lemma 2.5 in [34].□

Note that (2.7) can be stated as ‖ A g ‖ 1 , Ω ≤ C ‖ g ‖ − 1 2 , ∂ Ω , then

‖ A g ‖ 1 , Ω ≤ C ‖ g ‖ − 1 2 , ∂ Ω ≤ C ‖ g ‖ 0 , ∂ Ω .

From A h = P h A and (2.13), we have

(2.14) ‖ A h g ‖ 1 , Ω ≤ C ‖ A g ‖ 1 , Ω ≤ C ‖ g ‖ 0 , ∂ Ω ,

and from the trace theorem, it is valid that

‖ A h g ‖ 1 , Ω ≤ C C t r ‖ g ‖ 1 , Ω .

Lemma 2.4

When h → 0 , ∥ T h − T ∥ L 2 ( ∂ Ω ) → L 2 ( ∂ Ω ) → 0 and T is compact.

Proof

Referring to Theorem 3.2 in [8], we can prove the conclusions.□

Suppose that { λ p } and { λ p , h } are enumerations of the eigenvalues of (2.2) and (2.3), respectively, according to the same sort rule, each repeated as many times as its multiplicity, and λ = λ j is the j th eigenvalue of (2.2) with the algebraic multiplicity q , λ = λ j = λ j + 1 = ⋅ ⋅ ⋅ = λ j + q − 1 . Since T h converges to T , q eigenvalues λ j , h , λ j + 1 , h , ⋯ , λ j + q − 1 , h of (2.3) will converge to λ . Let M ( λ ) be the space spanned by all eigenfunctions of (2.2) that is corresponding to λ , and M h ( λ ) be the space spanned by all generalized eigenfunctions of (2.3) corresponding to λ p , h ( p = j , j + 1 , … , j + q − 1 ) .

Let M ^ ( λ ) = { u : u ∈ M ( λ ) , ‖ u ‖ 0 , ∂ Ω = 1 } , and denote

δ h ( λ ) = sup u ∈ M ^ ( λ ) inf v ∈ V h ‖ u − v ‖ 1 , Ω .

It is clear that δ h ( λ ) ≤ ρ h , and from (2.9) and the interpolation error estimates [35], we know that ρ h ≤ C h r , δ h ( λ ) ≤ C h r . From Lemma 3.3 in [36], we have ρ h → 0 ( h → 0 ) .

Theorem 2.5

Let λ and λ h be the j th eigenvalue of (2.2) and (2.3), respectively. Suppose that M ( λ ) ⊂ H 1 + t ( Ω ) ( t ≥ r ) , then

(2.15) ∣ λ − λ h ∣ ≤ C δ h 2 ( λ ) ,

For any eigenfunction u h corresponding to λ h , ‖ u h ‖ 0 , ∂ Ω = 1 , there exists u ∈ M ( λ ) and h 0 > 0 , such that when h ≤ h 0 , it is valid that

(2.16) ‖ u − u h ‖ 1 , Ω ≤ C δ h ( λ ) ,

(2.17) ‖ u − u h ‖ 0 , ∂ Ω ≤ C ρ h δ h ( λ ) .

For any u ∈ M ^ ( λ ) , there exists u h ∈ M h ( λ ) such that when h ≤ h 0 , there holds

(2.18) ‖ u h − u ‖ 1 , Ω ≤ C δ h ( λ ) .

Proof

Since ‖ T h − T ‖ L 2 ( ∂ Ω ) → L 2 ( ∂ Ω ) → 0 , from Theorems 7.3 and 7.4 in [37] and Lemma 2.3 in [38], we know that there exists u ∈ M ( λ ) , such that

‖ u h − u ‖ 0 , ∂ Ω ≤ C ‖ ( A − A h ) u ‖ 0 , ∂ Ω , λ − λ h = ⟨ ( A − A h ) u , u ⟩ + R ,

where ∣ R ∣ ≤ C ‖ ( A − A h ) u ‖ 0 , ∂ Ω 2 .

Then

∣ λ − λ h ∣ ≤ ∣ ⟨ ( A − A h ) u , u ⟩ ∣ + ∣ R ∣ ≤ ∣ ⟨ u , ( A − A h ) u ⟩ ∣ + ∣ R ∣ ≤ ∣ a ( A u , ( A − A h ) u ) ∣ + ∣ R ∣ ≤ ∣ a ( A u − A h u , ( A − A h ) u ) ∣ + ∣ R ∣ ≤ C ∥ ( A − A h ) u ∥ 1 , Ω 2 ≤ C sup u ∈ M ^ ( λ ) ∥ ( A − P h A ) u ∥ 1 , Ω 2 ≤ C sup u ∈ M ^ ( λ ) ∥ u − P h u ∥ 1 , Ω 2 ≤ C sup u ∈ M ^ ( λ ) inf v ∈ V h ∥ u − v ∥ 1 , Ω 2 ≤ C δ h 2 ( λ ) ,

then (2.15) holds.

From (2.12), we can deduce that

∥ u h − u ∥ 0 , ∂ Ω ≤ C ∥ ( A − A h ) u ∥ 0 , ∂ Ω ≤ C sup u ∈ M ^ ( λ ) ∥ ( P h A − A ) u ∥ 0 , ∂ Ω ≤ C sup u ∈ M ^ ( λ ) ρ h ∥ ( P h A − A ) u ∥ 1 , Ω ≤ C ρ h sup u ∈ M ^ ( λ ) ∥ P h u − u ∥ 1 , Ω ≤ C ρ h sup u ∈ M ^ ( λ ) inf v ∈ V h ∥ u − v ∥ 1 , Ω ≤ C ρ h δ h ( λ ) ,

i.e., (2.17) holds.

From the definition of A and A h , it is easy to see that u = − λ A u and u h = − λ h A h u h . A simple calculation shows that

(2.19) u − u h = − λ A u + λ h A h u h = − λ A u + λ A h u − λ A h u + λ h A h u h = ( A h − A ) ( λ u ) + A h ( λ h u h − λ u ) .

From (2.11) and the definition of δ h ( λ ) , we deduce that

(2.20) ∥ ( A − A h ) u ∥ 1 , Ω ≤ C sup u ∈ M ^ ( λ ) ∥ ( A − A h ) u ∥ 1 , Ω ≤ C sup u ∈ M ^ ( λ ) ∥ ( A − P h A ) u ∥ 1 , Ω ≤ C sup u ∈ M ^ ( λ ) ∥ u − P h u ∥ 1 , Ω ≤ C sup u ∈ M ^ ( λ ) inf v ∈ V h ∥ u − v ∥ 1 , Ω ≤ C δ h ( λ ) ;

and from (2.14), (2.15), and (2.17), we have

(2.21) ‖ A h ( λ h u h − λ u ) ‖ 1 , Ω ≤ C ‖ λ h u h − λ u ‖ 0 , ∂ Ω ≤ C ( ∣ λ h − λ ∣ + ‖ u h − u ‖ 0 , ∂ Ω ) ≤ C ρ h δ h ( λ ) .

Hence, from (2.20), (2.21), and (2.19), we derive (2.16).

By using a similar proof method, we can prove (2.18). The proof is completed.□

3 Multigrid discretization

3.1 Multigrid scheme

Let { π h i } 0 l be a family of regular meshes, h i − 1 ≫ h i , V h i be the conforming finite element space defined on π h i satisfying V h 0 = V H , V h i ⊂ V h i + 1 ⊂ H 1 ( Ω ) , i = 0 , 1 , … l − 1 , and ρ h i → 0 ( h i → 0 ) . We give the following multigrid discretization scheme based on the shifted inverse iteration for solving (2.2).

Scheme 1

Given the iteration time l and tolerance ε .

Step 1. Solve (2.2) in V H : Find ( λ H , u H ) ∈ R × V H such that ‖ u H ‖ 0 , ∂ Ω = 1 and

(3.1) a ( u H , v ) = − λ H ⟨ u H , v ⟩ ∀ v ∈ V H .

Step 2. u h 0 ⇐ u H , λ h 0 ⇐ λ H , i ⇐ 1 .

Step 3. Solve a linear system in V h i : Find u ′ ∈ V h i such that

a ( u ′ , v ) + λ h i − 1 ⟨ u ′ , v ⟩ = − ⟨ u h i − 1 , v ⟩ ∀ v ∈ V h i .

Set u h i = u ′ ‖ u ′ ‖ 0 , ∂ Ω .

Step 4. Compute the Rayleigh quotient:

λ h i = − a ( u h i , u h i ) ⟨ u h i , u h i ⟩ .

Step 5. If ∣ λ h i − λ h i − 1 ∣ ≤ ε or i = l , then output ( λ h i , u h i ) or ( λ h l , u h l ) and stop; Else, i ⇐ i + 1 and return to Step 3.

Let ( λ H , u H ) be the j th eigenpair of (3.1), then ( λ h l , u h l ) obtained by Scheme 1 is the j th approximate eigenpair of (2.2).

3.2 Theoretical analysis

In this section, we will conduct an error analysis of multigrid scheme. First, we introduce the following tools that are needed in the sequel.

Lemma 3.1

Let ( λ , u ) be an eigenpair of (2.2), then for any v ∈ H 1 ( Ω ) , ‖ v ‖ 0 , ∂ Ω ≠ 0 , the Rayleigh quotient R ( v ) = − a ( v , v ) ‖ v ‖ 0 , ∂ Ω 2 satisfies

R ( v ) − λ = − a ( v − u , v − u ) ‖ v ‖ 0 , ∂ Ω 2 − λ ⟨ v − u , v − u ⟩ ‖ v ‖ 0 , ∂ Ω 2 .

Proof

See the proof on page 699 in [37].□

Lemma 3.2

For any nonzero elements u , v in any linear normed space ( V , ‖ ⋅ ‖ ) , it is valid that

u ‖ u ‖ − v ‖ v ‖ ≤ 2 ‖ u − v ‖ ‖ u ‖ , u ‖ u ‖ − v ‖ v ‖ ≤ 2 ‖ u − v ‖ ‖ v ‖ .

Proof

See Lemma 3.1 in [22].□

In the coming discussion, let ( λ j , u j ) and ( λ j , h , u j , h ) denote the j th eigenpair of (2.2) and (2.3), respectively, and μ j = − 1 λ j , μ j , h = − 1 λ j , h , M ( μ j ) = M ( λ j ) , M h ( μ j ) = M h ( λ j ) .

Denote d = dim V h , dist ( u , W ) = inf v ∈ W ‖ u − v ‖ 1 , Ω . The following lemma is an extension of Theorem 3 in [22] and Lemma 4.1 in [29].

Lemma 3.3

Suppose that ( μ 0 , w 0 ) is an approximation of the j th eigenpair ( μ , u ) , where μ 0 is not an eigenvalue of A h and w 0 ∈ V h , ‖ w 0 ‖ 0 , ∂ Ω = 1 . Let u 0 = A h w 0 ‖ A h w 0 ‖ 0 , ∂ Ω , and suppose that

inf v ∈ M h ( λ ) ‖ w 0 − v ‖ 0 , ∂ Ω ≤ 1 2 ;
∣ μ 0 − μ ∣ ≤ ρ 4 , ∣ μ p , h − μ p ∣ ≤ ρ 4 , p = j − 1 , j , j + q ( p ≠ 0 ) , where ρ = min p ≠ j ∣ μ p − μ ∣ is the separation constant of the j th eigenvalue μ ;
u ′ ∈ V h , u h ∈ V h satisfy
(3.2) ( μ 0 − A h ) u ′ = u 0 , u h = u ′ ∥ u ′ ∥ 0 , ∂ Ω .

Then

dist ( u h , M h ( λ ) ) ≤ C ρ max j ≤ p ≤ j + q − 1 ∣ μ 0 − μ p , h ∣ dist ( w 0 , M h ( λ ) ) .

Proof

Let { u p , h } p = 1 d be eigenfunctions of A h satisfying ⟨ u p , h , u i , h ⟩ = δ p , i , then

u 0 = ∑ p = 1 d ⟨ u 0 , u p , h ⟩ u p , h .

Since μ 0 is not an eigenvalue of A h , and from (3.2), we can obtain

(3.3) ( μ 0 − μ j , h ) u ′ = ( μ 0 − μ j , h ) ( μ 0 − A h ) − 1 u 0 = ∑ p = 1 d μ 0 − μ j , h μ 0 − μ p , h ⟨ u 0 , u p , h ⟩ u p , h .

By using the triangle inequality and the condition (C2), we have

∣ μ 0 − μ j , h ∣ ≤ ∣ μ 0 − μ ∣ + ∣ μ − μ j , h ∣ ≤ ρ 4 + ρ 4 = ρ 2 , ∣ μ 0 − μ p , h ∣ ≥ ∣ μ − μ p ∣ − ∣ μ 0 − μ ∣ − ∣ μ p − μ p , h ∣ ≥ ρ − ρ 4 − ρ 4 = ρ 2 ,

where p = j − 1 , j + q ( p ≠ 0 ) . Hence, we obtain

(3.4) ∣ μ 0 − μ p , h ∣ ≥ ρ 2 p ≠ j , j + 1 , ⋯ , j + q − 1 .

Because T h is self-adjoint with respect to ⟨ ⋅ , ⋅ ⟩ [4], and A h u h = μ h u h , T h u h = μ h u h , for any p = 1 , 2 , ⋯ , d , it is valid that

(3.5) ⟨ T h w 0 , u p , h ⟩ u p , h = ⟨ w 0 , T h u p , h ⟩ u p , h = ⟨ w 0 , μ p , h u p , h ⟩ u p , h = ⟨ w 0 , u p , h ⟩ μ p , h u p , h = ⟨ w 0 , u p , h ⟩ A h u p , h .

Note that { u p , h } p = j j + q − 1 is a standard orthogonal basis of M h ( λ ) in the sense of L 2 ( ∂ Ω ) -inner product ⟨ ⋅ , ⋅ ⟩ , from u 0 = A h w 0 ‖ A h w 0 ‖ 0 , ∂ Ω , (3.3), (3.5), (2.14), and (3.4), we deduce that

(3.6) ( μ 0 − μ j , h ) u ′ − ∑ p = j j + q − 1 μ 0 − μ j , h μ 0 − μ p , h ⟨ u 0 , u p , h ⟩ u p , h 1 , Ω = ∑ p ≠ j , j + 1 , ⋯ , j + q − 1 μ 0 − μ j , h μ 0 − μ p , h ⟨ u 0 , u p , h ⟩ u p , h 1 , Ω = 1 ∥ A h w 0 ∥ 0 , ∂ Ω ∑ p ≠ j , j + 1 , ⋯ , j + q − 1 μ 0 − μ j , h μ 0 − μ p , h ⟨ A h w 0 , u p , h ⟩ u p , h 1 , Ω = 1 ∥ A h w 0 ∥ 0 , ∂ Ω ∑ p ≠ j , j + 1 , ⋯ , j + q − 1 μ 0 − μ j , h μ 0 − μ p , h ⟨ T h w 0 , u p , h ⟩ u p , h 1 , Ω = 1 ∥ A h w 0 ∥ 0 , ∂ Ω ∑ p ≠ j , j + 1 , ⋯ , j + q − 1 μ 0 − μ j , h μ 0 − μ p , h ⟨ w 0 , u p , h ⟩ A h u p , h 1 , Ω = 1 ∥ A h w 0 ∥ 0 , ∂ Ω A h ∑ p ≠ j , j + 1 , ⋯ , j + q − 1 μ 0 − μ j , h μ 0 − μ p , h ⟨ w 0 , u p , h ⟩ u p , h 1 , Ω ≤ C ∥ A h w 0 ∥ 0 , ∂ Ω ∑ p ≠ j , j + 1 , ⋯ , j + q − 1 μ 0 − μ j , h μ 0 − μ p , h ⟨ w 0 , u p , h ⟩ u p , h 0 , ∂ Ω ≤ 2 C ρ ∥ A h w 0 ∥ 0 , ∂ Ω ∣ μ 0 − μ j , h ∣ ∑ p ≠ j , j + 1 , ⋯ , j + q − 1 ⟨ w 0 , u p , h ⟩ 2 1 2 ≤ 2 C ρ ∥ A h w 0 ∥ 0 , ∂ Ω ∣ μ 0 − μ j , h ∣ w 0 − ∑ p = j j + q − 1 ⟨ w 0 , u p , h ⟩ u p , h 0 , ∂ Ω = 2 C ρ ∥ A h w 0 ∥ 0 , ∂ Ω ∣ μ 0 − μ j , h ∣ inf v ∈ M h ( λ ) ∥ w 0 − v ∥ 0 , ∂ Ω ≤ 2 C ρ ∥ A h w 0 ∥ 0 , ∂ Ω ∣ μ 0 − μ j , h ∣ dist ( w 0 , M h ( λ ) ) .

Taking norm on both sides of (3.3), and noticing that u 0 = A h w 0 ‖ A h w 0 ‖ 0 , ∂ Ω , from the condition (C1) and (3.5), we derive that

(3.7) ∥ ( μ 0 − μ j , h ) u ′ ∥ 0 , ∂ Ω = 1 ∥ A h w 0 ∥ 0 , ∂ Ω ∑ p = 1 d μ 0 − μ j , h μ 0 − μ p , h ⟨ A h w 0 , u p , h ⟩ u p , h 0 , ∂ Ω = 1 ∥ A h w 0 ∥ 0 , ∂ Ω ∑ p = 1 d μ 0 − μ j , h μ 0 − μ p , h ⟨ T h w 0 , u p , h ⟩ u p , h 0 , ∂ Ω = 1 ∥ A h w 0 ∥ 0 , ∂ Ω ∑ p = 1 d μ 0 − μ j , h μ 0 − μ p , h ⟨ w 0 , μ p , h u p , h ⟩ 2 1 2 ≥ C ∥ A h w 0 ∥ 0 , ∂ Ω min j ≤ p ≤ j + q − 1 μ 0 − μ j , h μ 0 − μ p , h ∑ p = j j + q − 1 ⟨ w 0 , u p , h ⟩ 2 1 2 = C ∥ A h w 0 ∥ 0 , ∂ Ω min j ≤ p ≤ j + q − 1 μ 0 − μ j , h μ 0 − μ p , h w 0 − w 0 − ∑ p = j j + q − 1 ⟨ w 0 , u p , h ⟩ u p , h 0 , ∂ Ω ≥ C 2 ∥ A h w 0 ∥ 0 , ∂ Ω min j ≤ p ≤ j + q − 1 μ 0 − μ j , h μ 0 − μ p , h .

From (3.6) and (3.7), we obtain

dist ( u h , M h ( λ ) ) = dist ( sign ( μ 0 − μ j , h ) u h , M h ( λ ) ) ≤ sign ( μ 0 − μ j , h ) u h − 1 ∥ ( μ 0 − μ j , h ) u ′ ∥ 0 , ∂ Ω ∑ p = j j + q − 1 μ 0 − μ j , h μ 0 − μ p , h ⟨ u 0 , u p , h ⟩ u p , h 1 , Ω = ( μ 0 − μ j , h ) u ′ ∥ ( μ 0 − μ j , h ) u ′ ∥ 0 , ∂ Ω − 1 ∥ ( μ 0 − μ j , h ) u ′ ∥ 0 , ∂ Ω ∑ p = j j + q − 1 μ 0 − μ j , h μ 0 − μ p , h ⟨ u 0 , u p , h ⟩ u p , h 1 , Ω ≤ 2 C ∥ A h w 0 ∥ 0 , ∂ Ω max j ≤ p ≤ j + q − 1 μ 0 − μ p , h μ 0 − μ j , h ( μ 0 − μ j , h ) u ′ − ∑ p = j j + q − 1 μ 0 − μ j , h μ 0 − μ p , h ⟨ u 0 , u p , h ⟩ u p , h 1 , Ω ≤ C ρ max j ≤ p ≤ j + q − 1 ∣ μ 0 − μ p , h ∣ dist ( w 0 , M h ( λ ) ) .

The proof is completed.□

Next we will use Lemma 3.3 and Theorem 2.5 to analyze Scheme 1. We first consider the case of l = 1 . Denote H = h 0 , h = h 1 .

Theorem 3.4

Suppose that M ( λ j ) ⊂ H 1 + t ( Ω ) ( t ≥ r ) . Let ( λ j h , u j h ) be an approximate eigenpair obtained by Scheme 1 ( l = 1 ) and H is sufficiently small, then there exists u j ∈ M ( λ j ) such that

(3.8) ∥ u j h − u j ∥ 1 , Ω ≤ C ( δ H 3 ( λ j ) + δ h ( λ j ) ) ,

(3.9) ∥ u j h − u j ∥ 0 , ∂ Ω ≤ C ( δ H 3 ( λ j ) + ρ h δ h ( λ j ) ) ,

(3.10) ∣ λ j h − λ j ∣ ≤ C ( δ H 6 ( λ j ) + δ h 2 ( λ j ) ) .

Proof

Select μ 0 = − 1 λ H , w 0 = u H and u 0 = A h u H ‖ A h u H ‖ 0 , ∂ Ω . From (2.16) and (2.17), we know that there exists u ¯ ∈ M ( λ j ) such that

∥ u H − u ¯ ∥ 1 , Ω ≤ C δ H ( λ j ) , ∥ u H − u ¯ ∥ 0 , ∂ Ω ≤ C ρ H δ H ( λ j ) ,

so, by the triangle inequality and (2.18), we have

(3.11) dist ( u H , M h ( λ j ) ) ≤ ∥ u H − u ¯ ∥ 1 , Ω + dist ( u ¯ , M h ( λ j ) ) ≤ C ( δ H ( λ j ) + δ h ( λ j ) ) ≤ C δ H ( λ j ) ,

thus,

inf v ∈ M h ( λ j ) ∥ u H − v ∥ 0 , ∂ Ω ≤ C δ H ( λ j ) .

When H is sufficiently small, the condition (C1) of Lemma 3.3 is valid.

From (2.15), we have

∣ μ 0 − μ j ∣ = ∣ λ H − λ j ∣ ∣ λ H λ j ∣ ≤ C δ H 2 ( λ j ) ≤ ρ 4 ; ∣ μ p − μ p , h ∣ = ∣ λ p , h − λ p ∣ ∣ λ p , h λ p ∣ ≤ C δ h 2 ( λ j ) ≤ ρ 4 , p = j − 1 , j , … , j + q , p ≠ 0 ,

i.e., the condition (C2) of Lemma 3.3 holds.

According to (2.8), we know that Step 3 of Scheme 1 is equivalent to the following:

a ( u ′ , v ) + λ H a ( A h u ′ , v ) = − a ( A h u H , v ) ∀ v ∈ V h ,

u j h = u ′ ∥ u ′ ∥ 0 , ∂ Ω , that is,

( λ H − 1 + A h ) u ′ = − λ H − 1 A h u H , u j h = u ′ ∥ u ′ ∥ 0 , ∂ Ω .

And − λ H − 1 A h u H and u 0 differ by only one constant, then Step 3 of Scheme 1 is also equivalent to

( μ 0 − A h ) u ′ = u 0 , u j h = u ′ ∥ u ′ ∥ 0 , ∂ Ω .

So far, all the conditions of Lemma 3.3 are valid.

Since M h ( λ j ) is a q -dimensional space, there must exist u ∗ ∈ M h ( λ j ) such that

∥ u j h − u ∗ ∥ 1 , Ω = dist ( u j h , M h ( λ j ) ) .

While for p = j , j + 1 , … , j + q − 1 , from (2.15), we have

(3.12) ∣ μ 0 − μ p , h ∣ = − 1 λ H + 1 λ p , h ≤ ∣ λ H − λ p , h ∣ ∣ λ H λ p , h ∣ ≤ C ( ∣ λ H − λ j ∣ + ∣ λ j − λ p , h ∣ ) ≤ C δ H 2 ( λ j ) .

Hence, from Lemma 3.3, (3.11), and (3.12), we obtain

(3.13) ‖ u j h − u ∗ ‖ 1 , Ω = dist ( u j h , M h ( λ j ) ) ≤ C ρ max j ≤ p ≤ j + q − 1 ∣ μ 0 − μ p , h ∣ dist ( u H , M h ( λ j ) ) ≤ C δ H 3 ( λ j ) .

From (2.16), we know that there is u j ∈ M ( λ j ) such that ‖ u ∗ − u j ‖ 1 , Ω = dist ( u ∗ , M ( λ j ) ) and

‖ u ∗ − u j ‖ 1 , Ω ≤ C δ h ( λ j ) ,

thus,

‖ u j h − u j ‖ 1 , Ω ≤ ‖ u j h − u ∗ ‖ 1 , Ω + ‖ u ∗ − u j ‖ 1 , Ω ≤ C ( δ H 3 ( λ j ) + δ h ( λ j ) ) ,

i.e., (3.8) holds. Next we will prove (3.9).

It follows from (2.17) that

‖ u ∗ − u j ‖ 0 , ∂ Ω ≤ C ρ h δ h ( λ j ) ,

which together with (3.13) yields

‖ u j h − u j ‖ 0 , ∂ Ω ≤ ‖ u j h − u ∗ ‖ 0 , ∂ Ω + ‖ u ∗ − u j ‖ 0 , ∂ Ω ≤ C ( δ H 3 ( λ j ) + ρ h δ h ( λ j ) ) .

Finally, we use Lemma 3.1 to derive (3.10). From Step 4 of Scheme 1, Lemma 3.1, (3.8), and (3.9), we deduce that

∣ λ j h − λ j ∣ = − a ( u j h − u j , u j h − u j ) ‖ u j h ‖ 0 , ∂ Ω 2 − λ j ⟨ u j h − u j , u j h − u j ⟩ ‖ u j h ‖ 0 , ∂ Ω 2 ≤ C ( ‖ u j h − u j ‖ 1 , Ω 2 + ∣ λ j ∣ ‖ u j h − u j ‖ 0 , ∂ Ω 2 ) ≤ C ( δ H 6 ( λ j ) + δ h 2 ( λ j ) ) .

The proof is completed.□

Remark

By using Theorem 3.4 and referring Theorem 4.2 in [29], we can establish the error estimates of Scheme 1. To ensure that the error is independent of the iteration times during the multigrid refinement, we also need to introduce the following condition.

Condition 3.1

For any given number ε ∈ ( 0 , 1 ) , there exists t i ∈ ( 1 , 2 − ε ] ( i = 1 , 2 , ⋯ ) , such that h i = O ( h i − 1 t i ) and h i → 0 ( i → ∞ ) .

Condition 3.1 is easy to be satisfied. For example, for smooth eigenfunctions, using uniform meshes and linear elements and taking ε = 0.1 , h 0 = 2 8 , h 1 = 2 32 , h 2 = 2 128 , … , then t i = log ( h i ) log ( h i − 1 ) = log ( h i − 1 ) − log ( 4 ) log ( h i − 1 ) , thus t 1 ≈ 1.80 , t 2 ≈ 1.44 , t 3 ≈ 1.31 , … , and t i ↘ 1 when i → ∞ .

Theorem 3.5

Suppose that Condition 3.1 holds, and M ( λ j ) ⊂ H 1 + t ( Ω ) ( t ≥ r ) . Let ( λ j h l , u j h l ) be the approximation eigenpair obtained by Scheme 1, then, when h 0 = H is sufficiently small, there exists u j ∈ M ( λ j ) such that

∥ u j h l − u j ∥ 1 , Ω ≤ C δ h l ( λ j ) , ∥ u j h l − u j ∥ 0 , ∂ Ω ≤ C ρ h l δ h l ( λ j ) , ∣ λ j h l − λ j ∣ ≤ C δ h l 2 ( λ j ) , l ≥ 1 .

4 Numerical experiments

In this section, we will report some numerical experiments to show the efficiency of the multigrid scheme. Consider the Steklov eigenvalue problem (2.1) with k = 1 and n ( x ) = 4 in the test domains Ω S = − 2 2 , 2 2 2 and Ω L = ( − 1 , 1 ) 2 \ ( [ 0 , 1 ) × ( − 1 , 0 ] ) . We use MATLAB 2019a to complete our numerical examples on an Lenovo ST558 PC with 16G memory. Our programs are compiled under the package of iFEM [39], and in computation, the solver eigs in MATLAB is used in Step 1 for eigenpairs of the discrete algebraic eigenvalue problem (3.1), and the command “ \ ” in MATLAB is used to solve the linear system in Step 3.

Example 1. We use Scheme 1, the multigrid scheme based on the shifted inverse iteration, by adopting linear element to compute. The results are listed in Tables 1 and 2. In Tables 1 and 2, we use the following notations.

Table 1

The results obtained by Scheme 1 with linear element in Ω S

h 0	h 1	h 2	h 3	λ 1 h 3	t ( s )	λ 1 , h 3	t ( s )
1 8	1 32	1 128	1 2,048	2.20250676799282	257.13	2.20250676813995	294.93
1 16	1 64	1 256	1 2,048	2.20250676799215	263.56	2.20250676813991	300.04
1 32	1 128	1 512	1 2,048	2.20250676799231	267.19	2.20250676813991	301.93

h 0	h 1	h 2	h 3	λ 2 h 3	t ( s )	λ 2 , h 3	t ( s )
1 8	1 32	1 128	1 2,048	‒ 0.21225231703229	268.57	‒ 0.21225235299975	299.39
1 16	1 64	1 256	1 2,048	‒ 0.21225231703244	271.42	‒ 0.21225235299959	299.15
1 32	1 128	1 512	1 2,048	‒ 0.21225231703235	284.58	‒ 0.21225235299983	301.60

h 0	h 1	h 2	h 3	λ 4 h 3	t ( s )	λ 4 , h 3	t ( s )
1 8	1 32	1 128	1 2,048	‒ 0.90805674483814	251.37	‒ 0.90805674481081	296.68
1 16	1 64	1 256	1 2,048	‒ 0.90805674483789	262.83	‒ 0.90805674481087	299.49
1 32	1 128	1 512	1 2,048	‒ 0.90805674483799	269.69	‒ 0.90805674481053	295.96

Table 2

The results obtained by Scheme 1 with linear element in Ω L

h 0	h 1	h 2	h 3	λ 1 h 3	t ( s )	λ 1 , h 3	t ( s )
2 8	2 32	2 128	2 2,048	2.53321335004984	872.58	—	—
2 16	2 64	2 256	2 2,048	2.53321335004904	908.62	—	—
2 32	2 128	2 512	2 2,048	2.53321335004911	956.47	—	—

h 0	h 1	h 2	h 3	λ 2 h 3	t ( s )	λ 2 , h 3	t ( s )
2 8	2 32	2 128	2 2,048	0.85778318273331	848.64	—	—
2 16	2 64	2 256	2 2,048	0.85778318273325	929.15	—	—
2 32	2 128	2 512	2 2,048	0.85778318273359	966.22	—	—

h 0	h 1	h 2	h 3	λ 3 h 3	t ( s )	λ 3 , h 3	t ( s )
2 8	2 32	2 128	2 2,048	0.12452433218469	849.56	—	—
2 16	2 64	2 256	2 2,048	0.12452433218501	844.72	—	—
2 32	2 128	2 512	2 2,048	0.12452433218493	889.48	—	—

h i : the mesh size of π h i ;

λ j , h i : the j th eigenvalue of (3.1) on π h i by directly solving. Here, we first use eigs ( K h 0 , ℳ h 0 , ι , ′ s m ′ ) to obtain λ j , h 0 , then use eigs ( K h i , ℳ h i , 1 , λ j , h 0 ) to obtain λ j , h i , where K h i and ℳ h i are the stiffness matrix and mass matrix, respectively, of generalized matrix eigenvalue problem form associated with (3.1) in V h i ;

λ j h i : the j th eigenvalue obtained by Scheme 1 with π h i being the mesh at iteration stop;

t ( s ) : the CPU time from the program starting to the appearance of computed results;

–: the calculation cannot proceed since the computer runs out of memory.

From Tables 1 and 2, we can see that using the multigrid scheme can obtain approximate eigenvalues with the same accuracy as those obtained by directly solving with less CPU time. Especially from Table 2, it can be seen that in the L-shaped domain, when direct solving cannot proceed, the multigrid scheme can still iterate once to obtain the results.

Example 2. Adaptive computation.

Adaptive algorithm based on the a posteriori error estimate is a popular and efficient numerical method for solving partial differential equations [40–45]. Referring to [10], we combine Scheme 1 and a posteriori error indicators to establish the adaptive multigrid algorithm.

For any element τ ∈ π h , let ℰ τ denote the set of all edges of τ . Denote ℰ ≔ ⋃ τ ∈ π h ℰ τ , ℰ Γ ≔ { ℓ ∈ ℰ : ℓ ⊂ ∂ Ω } , ℰ Ω ≔ ℰ \ ℰ Γ . It is obvious that ℰ Γ and ℰ Ω are disjoint and ℰ = ℰ Γ ⋃ ℰ Ω .

For each ℓ ∈ ℰ Ω , denote the two elements that sharing ℓ as τ in and τ out , and take a unit normal ν ℓ on ℓ that points outwards τ in .

For v h ∈ V h , define the jump of normal derivative of v h on ℓ as follows:

∂ v h ∂ ν ℓ ≕ ∇ ( v h ∣ τ out ) ⋅ ν ℓ − ∇ ( v h ∣ τ in ) ⋅ ν ℓ , ℓ ∈ ℰ Ω .

For each ℓ ∈ ℰ , define

J ℓ ( λ h , u h ) ≔ 1 2 ∂ u h ∂ ν ℓ , ℓ ∈ ℰ Ω , − λ h u h − ∂ u h ∂ ν , ℓ ∈ ℰ Γ .

For each element τ ∈ π h , define the local error indicator η τ ( λ h , u h ) :

η τ ( λ h , u h ) ≔ h τ 2 ‖ Δ u h + k 2 n ( x ) u h ‖ 0 , τ 2 + ∑ ℓ ∈ ℰ τ ∣ ℓ ∣ ‖ J ℓ ‖ 0 , ℓ 2 1 2 ,

and the global error indicator η Ω :

η Ω ≔ ∑ τ ∈ π h η τ 2 1 2 .

On the basis of the aforementioned error indicators and Scheme 1, we design the following adaptive multigrid algorithm based on the shifted inverse iteration.

Algorithm 1.

Given the tolerance ε . Choose the parameter 0 < ε < 1 .

Step 1. Pick any initial mesh π h 1 with mesh size h 1 .

Step 2. Solve (2.2) on π h 1 for discrete solution ( λ j , h 1 , u j , h 1 ) .

Step 3. Let l = 1 , u j h 1 ⇐ u j , h 1 , λ j h 1 ⇐ λ j , h 1 .

Step 4. Compute the local error indicator η τ ( λ j h l , u j h l ) .

Step 5. Construct π ^ h l ⊂ π h l by Marking strategy and the parameter ε .

Step 6. Refine π h l to obtain a new mesh π h l + 1 by procedure REFINE.

Step 7. Solve the following linear system on π h l + 1 for u ′ ∈ V h l + 1 :

a ( u ′ , v ) + λ j h l ⟨ u ′ , v ⟩ = − ⟨ u j h l , v ⟩ ∀ v ∈ V h l + 1 ;

Set u j h l + 1 = u ′ ‖ u ′ ‖ 0 , ∂ Ω , and compute the Rayleigh quotient

λ j h l + 1 = − a ( u j h l + 1 , u j h l + 1 ) ⟨ u j h l + 1 , u j h l + 1 ⟩ .

Step 8. If ∣ λ h l + 1 − λ h l ∣ ≤ ε or l reaches the maximum iterations, output ( λ j h l + 1 , u j h l + 1 ) ; Else, let l = l + 1 , and return to Step 4.

Marking strategy Given the parameter 0 < ε < 1 .

Step 1. Pick some elements in π h l to construct a minimal subset π ^ h l of π h l , such that

∑ τ ∈ π ^ h l η τ 2 ( λ j h l , u j h l ) ≥ ε η Ω 2 .

Step 2. Mark all the elements in π ^ h l .

Moreover, to investigate the efficiency of adaptive multigrid algorithm 1, we give the following Algorithm 2 by referring the standard popular adaptive method [46] for comparison.

Algorithm 2.

Given the tolerance ε . Choose the parameter 0 < ε < 1 .

Step 1. Pick any initial mesh π h 1 with mesh size h 1 .

Step 2. Solve (2.2) on π h 1 for discrete solution ( λ j , h 1 , u j , h 1 ) .

Step 3. Let l = 1 .

Step 4. Compute the local error indicator η τ ( λ j , h l , u j , h l ) .

Step 5. Construct π ^ h l ⊂ π h l by Marking Strategy and the parameter ε .

Step 6. Refine π h l to obtain a new mesh π h l + 1 by procedure REFINE.

Step 7. Solve (2.2) on π h l + 1 for ( λ j , h l + 1 , u j , h l + 1 ) .

Step 8. If ∣ λ h l + 1 − λ h l ∣ ≤ ε or l reaches the maximum iterations, output ( λ j h l + 1 , u j h l + 1 ) ; Else, let l = l + 1 , and return to Step 4.

We use Algorithms 1 and 2 with linear element to compute the first four approximate eigenvalues of (2.1) in Ω S and Ω L , respectively. The diameter of initial mesh is taken as 1 32 for Ω S and 2 32 for Ω L , and the results are listed in Tables 3 and 4. In computation, the solver eigs in MATLAB is used in Steps 2 and 7 in adaptive algorithms for eigenvalues. Besides listing the approximate eigenvalues, we also depict the error curves of approximate eigenvalues and the curves of error indicator to observe the performance of our method directly. To do so, we use the approximations that we can calculate by direct solving as accurate as possible, as the reference values. In Figures 1 and 2 we show the error curves of approximate eigenvalues obtained by Algorithm 1 and the curves of error indicators. In the following tables and figures, we adopt the following notations:

Table 3

The results in Ω S by Algorithms 1 and 2 using linear element with ε = 0.25

j	l	N j , l	λ j h l ˜	CPU j , l ( s )	l	N j , l	λ j , h l ˜	CPU j , l ( s )
1	86	103,345	2.2024645152	70.81	86	103,345	2.2024645152	77.45
1	113	595,527	2.2025015517	469.84	113	595,527	2.2025015517	507.83
1	122	1,034,817	2.2025032754	877.39	122	1,034,817	2.2025032754	941.39
2	40	105,152	‒ 0.2122613227	41.76	45	106,543	‒ 0.2122613290	51.09
2	59	927,775	‒ 0.2122536155	487.43	64	904,241	‒ 0.2122533085	513.45
2	63	1,436,819	‒ 0.2122529469	801.15	68	1,447,751	‒ 0.2122527999	824.00
3	40	100,845	‒ 0.2122616692	41.48	53	107,130	‒ 0.2122613970	67.03
3	54	507,413	‒ 0.2122541916	254.78	76	562,880	‒ 0.2122541948	465.78
3	61	1,126,579	‒ 0.2122533001	612.48	84	1,132,214	‒ 0.2122533375	888.31
4	67	100,205	‒ 0.9080979586	53.11	67	100,205	‒ 0.9080979586	67.05
4	85	417,006	‒ 0.9080665587	278.01	85	417,006	‒ 0.9080665587	309.88
4	97	1,108,435	‒ 0.9080609575	779.72	97	1,108,435	‒ 0.9080609575	846.57

Table 4

The results in Ω L by Algorithms 1 and 2 using linear element with ε = 0.25

j	l	N j , l	λ j h l ˜	CPU j , l ( s )	l	N j , l	λ j , h l ˜	CPU j , l ( s )
1	81	105,764	2.5331030459	66.09	81	105,764	2.5331030459	81.43
1	101	404,798	2.5331816554	304.24	101	404,798	2.5331816554	333.49
1	114	951,656	2.5332007976	784.86	114	951,656	2.5332007977	831.87
2	51	98,242	0.8577393960	45.93	51	98,242	0.8577393960	56.44
2	66	378,630	0.8577763900	219.74	66	378,630	0.8577763901	235.12
2	77	976,497	0.8577845504	635.10	77	976,497	0.8577845504	655.17
3	26	83,171	0.1245088600	27.72	26	83,171	0.1245088600	32.96
3	35	250,423	0.1245197124	106.01	35	250,423	0.1245197124	118.18
3	50	1,556,949	0.1245235543	817.58	50	1,556,949	0.1245235543	865.44
4	73	98,743	‒ 1.0853708906	54.95	73	98,743	‒ 1.0853708906	68.64
4	85	243,700	‒ 1.0853297996	158.45	85	243,700	‒ 1.0853297996	175.09
4	106	1,149,078	‒ 1.0853050268	882.66	106	1,149,078	‒ 1.0853050268	914.02

$Figure 1 The error curves by Algorithm 1 using linear element in Ω S {\Omega }_{S} .$

Figure 1

The error curves by Algorithm 1 using linear element in Ω S .

$Figure 2 The error curves by Algorithm 1 using linear element in Ω L {\Omega }_{L} .$

Figure 2

The error curves by Algorithm 1 using linear element in Ω L .

l : the number of iterations;

λ j h l ˜ : the j th eigenvalue obtained by Algorithm 1 after l iterations;

λ j , h l ˜ : the j th eigenvalue obtained by Algorithm 2 after l iterations;

N j , l : the degrees of freedom at the l th iteration;

CPU j , l ( s ) : the CPU time(s) from the program starting to the calculate results of the l th iteration appearing;

e j : the error of the j th approximate eigenvalue obtained by Algorithm 1;

η j : the error indicator of the j th approximate eigenvalue by Algorithm 1.

From Tables 3 and 4, we can see that, comparing with Algorithm 2, Algorithm 1 takes a little less time to obtain the same accurate approximate eigenvalues. However, the advantage of the adaptive multigrid algorithm in saving time is not significant. We think that this may be because the mesh in Algorithm 1 is refined according to the error indicators and the a posteriori errors, so that the diameters of adjacent two meshes do not match the relationship in Condition 3.1. Another possible reason is that the problem considered is nonpositive definite for a ( ⋅ , ⋅ ) being not coercive. It also can be found in Figure 1 that the error curves e 1 , e 2 , e 3 , and e 4 are all parallel to the line with slope − 1 , which indicate that the first, the second, the third, and the fourth approximate eigenvalue achieve the optimal convergence order; in addition, the curves of η 1 , η 2 , η 3 , and η 4 are parallel to the error curves e 1 , e 2 , e 3 , and e 4 , respectively, which imply that the error indicators are reliable and effective. We can obtain the same conclusion from Figure 2.

Acknowledgments

The authors cordially thank the editor and the referees for their valuable comments and suggestions that lead to the improvement of this article, the authors would like to thank Professor Yang Yidu and Professor Jiayu Han of Guizhou Normal University for their guidance, and Liangkun Xu, a graduate student of Guizhou Normal University, for his assistance in this article.

Funding information: This work was supported by the National Natural Science Foundation of China (Nos. 11761024, 12261024).
Author contributions: Jiali Xie: methodology; investigation; validation; writing – original draft; software. Hai Bi: supervision; validation; investigation; methodology. All authors read and approved the final manuscript.
Conflict of interest: The authors state no conflict of interest.
Ethical approval: The conducted research is not related to either human or animals use.
Data availability statement: All data generated or analysed during this study are included in this published article.

References

[1] J. A. Canavati and A. A. Minzoni, A discontinuous Steklov problem with an application to water waves, J. Math. Anal. Appl. 69 (1979), no. 2, 540–558, DOI: https://doi.org/10.1016/0022-247X(79)90165-3. 10.1016/0022-247X(79)90165-3Search in Google Scholar

[2] L. Cao, L. Zhang, W. Allegretto, and Y. Lin, Multiscale asymptotic method for Steklov eigenvalue equations in composite media, SIAM J. Numer. Anal. 51 (2013), no. 1, 273–296, DOI: https://doi.org/10.1137/110850876. 10.1137/110850876Search in Google Scholar

[3] N. Kuznetsov, T. Kulczycki, M. Kwaśnicki, et al., The legacy of Vladimir Andreevich Steklov, Notices Amer. Math. Soc. 61 (2014), no. 1, 9–22, DOI: https://doi.org/10.1090/noti1073. 10.1090/noti1073Search in Google Scholar

[4] F. Cakoni, D. Colton, S. Meng, and P. Monk, Stekloff eigenvalues in inverse scattering, SIAM J. Appl. Math. 76 (2016), no. 4, 1737–1763, DOI: https://doi.org/10.1137/16M1058704. 10.1137/16M1058704Search in Google Scholar

[5] J. Liu, Y. Liu, and J. Sun, An inverse medium problem using Stekloff eigenvalues and a Bayesian approach, Inverse Problems 35 (2019), no. 9, 094004, DOI: https://doi.org/10.1088/1361-6420/ab1be9. 10.1088/1361-6420/ab1be9Search in Google Scholar

[6] I. Harris, Approximation of the inverse scattering Steklov eigenvalues and the inverse spectral problem, Res. Math. Sci. 8 (2021), 31, DOI: https://doi.org/10.1007/s40687-021-00268-1. 10.1007/s40687-021-00268-1Search in Google Scholar

[7] Y. Yang, Y. Zhang, and H. Bi, Non-conforming Crouzeix-Raviart element approximation for Stekloff eigenvalues in inverse scattering, Adv. Comput. Math. 46 (2020), 81, DOI: https://doi.org/10.1007/s10444-020-09818-7. 10.1007/s10444-020-09818-7Search in Google Scholar

[8] J. Liu, J. Sun, and T. Turner, Spectral indicator method for a non-selfadjoint Steklov eigenvalue problem, J. Sci. Comput. 79 (2019), 1814–1831, DOI: https://doi.org/10.1007/s10915-019-00913-6. 10.1007/s10915-019-00913-6Search in Google Scholar

[9] Y. Zhang, H. Bi, and Y. Yang, A multigrid correction scheme for a new Steklov eigenvalue problem in inverse scattering, Int. J. Comput. Math. 97 (2019), no. 7, 1412–1430, DOI: https://doi.org/10.1080/00207160.2019.1622686. 10.1080/00207160.2019.1622686Search in Google Scholar

[10] Y. Zhang, H. Bi, and Y. Yang, The adaptive finite element method for the Steklov eigenvalue problem in inverse scattering, Open Math. 18 (2020), no. 1, 216–236, DOI: https://doi.org/10.1515/math-2020-0140. 10.1515/math-2020-0140Search in Google Scholar

[11] J. Meng and L. Mei, Discontinuous Galerkin methods of the non-selfadjoint Steklov eigenvalue problem in inverse scattering, Appl. Math. Comput. 381 (2020), 125307, DOI: https://doi.org/10.1016/j.amc.2020.125307. 10.1016/j.amc.2020.125307Search in Google Scholar

[12] S. Ren, Y. Zhang, and Z. Wang, An efficient spectral-Galerkin method for a new Steklov eigenvalue problem in inverse scattering, AIMS Math. 7 (2022), no. 5, 7528–7551, DOI: https://doi.org/10.3934/math.2022423. 10.3934/math.2022423Search in Google Scholar

[13] J. Xu, A new class of iterative methods for nonself adjoint or indefinite problems, SIAM J. Numer. Anal. 29 (1992), no. 2, 303–319, DOI: https://doi.org/10.1137/0729020. 10.1137/0729020Search in Google Scholar

[14] J. Xu, Two-grid discretization techniques for linear and nonlinear PDEs, SIAM J. Numer. Anal. 33 (1996), no. 5, 1759–1777, DOI: https://doi.org/10.1137/S0036142992232949. 10.1137/S0036142992232949Search in Google Scholar

[15] M. R. Racheva and A. B. Andreev, Super convergence postprocessing for eigenvalues, Comput. Methods Appl. Math. 2 (2002), no. 2, 171–185, DOI: https://doi.org/10.2478/cmam-2002-0011. 10.2478/cmam-2002-0011Search in Google Scholar

[16] C. S. Chien and B. W. Jeng, A two-grid discretization scheme for semilinear elliptic eigenvalue problems, SIAM J. Sci. Comput. 27 (2006), no. 4, 1287–1304, DOI: https://doi.org/10.1137/030602447. 10.1137/030602447Search in Google Scholar

[17] X. Dai and A. Zhou, Three-scale finite element discretizations for quantum eigenvalue problems, SIAM J. Numer. Anal. 46 (2008), no. 1, 295–324, DOI: https://doi.org/10.1137/06067780X. 10.1137/06067780XSearch in Google Scholar

[18] H. Chen, S. Jia, and H. Xie, Postprocessing and higher order convergence for the mixed finite element approximations of the Stokes eigenvalue problems, Appl. Math. 54 (2009), no. 3, 237–250, DOI: https://doi.org/10.1007/s10492-009-0015-7. 10.1007/s10492-009-0015-7Search in Google Scholar

[19] H. Xie and X. Yin, Acceleration of stabilized finite element discretizations for the Stokes eigenvalue problem, Adv. Comput. Math. 41 (2015), 799–812, DOI: https://doi.org/10.1007/s10444-014-9386-8. 10.1007/s10444-014-9386-8Search in Google Scholar

[20] J. Chen, Y. Xu, and J. Zou, An adaptive inverse iteration for Maxwell eigenvalue problem based on edge elements, J. Comput. Phys. 229 (2010), no. 7, 2649–2658, DOI: https://doi.org/10.1016/j.jcp.2009.12.013. 10.1016/j.jcp.2009.12.013Search in Google Scholar

[21] A. B. Andreev, R. D. Lazarov, and M. R. Racheva, Postprocessing and higher order convergence of the mixed finite element approximations of biharmonic eigenvalue problems, J. Comput. Appl. Math. 182 (2005), no. 2, 333–349, DOI: https://doi.org/10.1016/j.cam.2004.12.015. 10.1016/j.cam.2004.12.015Search in Google Scholar

[22] Y. Yang and H. Bi, Two-grid finite element discretization schemes based on shifted-inverse power method for elliptic eigenvalue problems, SIAM J. Numer. Anal. 49 (2011), no. 4, 1602–1624, DOI: https://doi.org/10.1137/100810241. 10.1137/100810241Search in Google Scholar

[23] X. Hu, and X. Cheng, Acceleration of a two-grid method for eigenvalue problems, Math. Comp. 80 (2011), no. 275, 1287–1301, DOI: https://doi.org/10.1090/s0025-5718-2011-02458-0. 10.1090/S0025-5718-2011-02458-0Search in Google Scholar

[24] X. Hu, and X. Cheng, Corrigendum to: acceleration of a two-grid method for eigenvalue problems, Math. Comp. 84 (2015), no. 296, 2701–2704, DOI: https://doi.org/10.1090/mcom/2967. 10.1090/mcom/2967Search in Google Scholar

[25] H. Chen, Y. He, Y. Li, and H. Xie, A multigrid method for eigenvalue problems based on shifted-inverse power technique, Eur. J. Math. 1 (2015), 207–228, DOI: https://doi.org/10.1007/s40879-014-0034-0. 10.1007/s40879-014-0034-0Search in Google Scholar

[26] J. Zhou, X. Hu, S. Shu, L. Zhong, and L. Chen, Two-grid methods for Maxwell eigenvalue problems, SIAM J. Number. Anal. 52 (2014), no. 4, 2027–2047, DOI: https://doi.org/10.1137/130919921. 10.1137/130919921Search in Google Scholar PubMed PubMed Central

[27] J. Liu, W. Jiang, F. Lin, N. Liu, and Q. H. Liu, A two-grid vector discretization scheme for the resonant cavity problem with anisotropic media, IEEE Trans. Microw. Theory Tech. 65 (2017), no. 8, 2719–2725, DOI: https://doi.org/10.1109/TMTT.2017.2672545. 10.1109/TMTT.2017.2672545Search in Google Scholar

[28] J. Han, Z. Zhang, and Y. Yang, A new adaptive mixed finite element method based on residual type a posterior error estimates for the Stokes eigenvalue problem, Numer. Methods Partial Differential Equations 31 (2015), no. 1, 31–53, DOI: https://doi.org/10.1002/num.21891. 10.1002/num.21891Search in Google Scholar

[29] Y. Yang, H. Bi, J. Han, and Y. Yu, The shifted-inverse iteration based on the multigrid discretizations for eigenvalue problems, SIAM J. Sci. Comput. 37 (2015), no. 6, DOI: https://doi.org/10.1137/140992011. 10.1137/140992011Search in Google Scholar

[30] J. Xu and A. Zhou, Local and parallel finite element algorithms based on two-grid discretizations, Math. Comp. 69 (2000), no. 231, 881–909, DOI: https://doi.org/10.1090/s0025-5718-99-01149-7. 10.1090/S0025-5718-99-01149-7Search in Google Scholar

[31] G. C. Hsiao and W. L. Wendland, Boundary Integral Equations, Applied Mathematical Sciences, Vol. 164. Springer-Verlag, Berlin, 2008. 10.1007/978-3-540-68545-6Search in Google Scholar

[32] M. Dauge, Elliptic Boundary Value Problems on Corner Domains: Smoothness and Asymptotics of Solutions, Lecture Notes in Mathematics, Vol. 1341. Springer-Verlag, Berlin, 1988. 10.1007/BFb0086682Search in Google Scholar

[33] A. Alonso and A. D. Russo, Spectral approximation of variationally-posed eigenvalue problems by nonconforming methods, J. Comput. Appl. Math. 223 (2009), no. 1, 177–197, DOI: https://doi.org/10.1016/j.cam.2008.01.008. 10.1016/j.cam.2008.01.008Search in Google Scholar

[34] H. Bi, Y. Zhang, and Y. Yang, Two-grid discretizations and a local finite element scheme for a non-selfadjoint Stekloff eigenvalue problem, Comput. Math. Appl. 79 (2020), no. 7, 1895–1913, DOI: https://doi.org/10.1016/j.camwa.2018.08.047. 10.1016/j.camwa.2018.08.047Search in Google Scholar

[35] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, Springer-Verlag, New York, 2008. 10.1007/978-0-387-75934-0Search in Google Scholar

[36] I. Babuska and J. E. Osborn, Finite element-Galerkin approximation of the eigenvalues and eigenvectors of selfadjoint problems, Math. Comp. 52 (1989), no. 186, 275–297, DOI: https://doi.org/10.1090/s0025-5718-1989-0962210-8. 10.1090/S0025-5718-1989-0962210-8Search in Google Scholar

[37] I. Babuska and J. Osborn, Eigenvalue problems, in: Finite Element Methods Part 1, Elsevier, North-Holand, 1991. 10.1016/S1570-8659(05)80042-0Search in Google Scholar

[38] Y. Yang, Z. Zhang, and F. Lin, Eigenvalue approximation from below using non-conforming finite elements, Sci. China Math. 53 (2010), no. 1, 137–150, DOI: https://doi.org/10.1007/s11425-009-0198-0. 10.1007/s11425-009-0198-0Search in Google Scholar

[39] L. Chen, An integrated finite element method package in MATLAB, Technical Report, University of California at Irvine, California, 2009. Search in Google Scholar

[40] R. Becker and R. Rannacher, An optimal control approach to a posteriori error estimation in finite element methods, Acta Numerica 10 (2001), 1–102, DOI: https://doi.org/10.1017/S0962492901000010. 10.1017/S0962492901000010Search in Google Scholar

[41] R. Becker and B. Vexler, A posteriori error estimation for finite element discretization of parameter identification problems, Numer. Math. 96 (2004), 435–459, DOI: https://doi.org/10.1007/s00211-003-0482-9. 10.1007/s00211-003-0482-9Search in Google Scholar

[42] L. Beilina and C. Clason, An adaptive hybrid FEM/FDM method for an inverse scattering problem in scanning acoustic microscopy, SIAM J. Sci. Comput. 28 (2006), no. 1, 382–402, DOI: https://doi.org/10.1137/050631252. 10.1137/050631252Search in Google Scholar

[43] J. B. Malmberg and L. Beilina, An adaptive finite element method in quantitative reconstruction of small inclusions from limited observations, Appl. Math. Inform. Sci. 12 (2018), no. 1, 1–19, DOI: https://doi.org/10.18576/amis/120101. 10.18576/amis/120101Search in Google Scholar

[44] M. J. Grote, M. Kray, and U. Nahum, Adaptive eigenspace method for inverse scattering problems in the frequency domain, Inverse Problems 33 (2017), no. 2, 025006, DOI: https://doi.org/10.1088/1361-6420/aa5250. 10.1088/1361-6420/aa5250Search in Google Scholar

[45] M. J. Grote and U. Nahum, Adaptive eigenspace for multiparameter inverse scattering problems, Comput. Math. Appl. 77 (2019), no. 12, 3264–3280, DOI: https://doi.org/10.1016/j.camwa.2019.02.005. 10.1016/j.camwa.2019.02.005Search in Google Scholar

[46] X. Dai, J. Xu, and A. Zhou, Convergence and optimal complexity of adaptive finite element eigenvalue computations, Numer. Math. 110 (2008), 313–355, DOI: https://doi.org/10.1007/s00211-008-0169-3. 10.1007/s00211-008-0169-3Search in Google Scholar

Received: 2023-01-27

Revised: 2023-06-18

Accepted: 2023-06-25

Published Online: 2023-07-28

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/math-2022-0607

Keywords for this article

Steklov eigenvalue problem in inverse scattering; multigrid discretization; the shifted-inverse iteration; error estimation; adaptive algorithm

Creative Commons

BY 4.0