The adaptive finite element method for the Steklov eigenvalue problem in
 inverse scattering

Abstract In this study, for the first time, we discuss the posteriori error estimates and adaptive algorithm for the non-self-adjoint Steklov eigenvalue problem in inverse scattering. The differential operator corresponding to this problem is non-self-adjoint and the associated weak formulation is not H 1-elliptic. Based on the study of Armentano et al. [Appl. Numer. Math. 58 (2008), 593–601], we first introduce error indicators for primal eigenfunction, dual eigenfunction, and eigenvalue. Second, we use Gårding’s inequality and duality technique to give the upper and lower bounds for energy norm of error of finite element eigenfunction, which shows that our indicators are reliable and efficient. Finally, we present numerical results to validate our theoretical analysis.

Adaptive finite element methods are favored in the current scientific and engineering computations. For a given tolerance, adaptive finite element methods require less degrees of freedom. So far, many excellent works on the a posteriori error estimates and adaptive algorithm have been summarized in previous studies [23][24][25]. There are some studies on the a posteriori error estimates and adaptive finite element methods for self-adjoint Steklov eigenvalue problems (see, e.g., [11][12][13]17,[26][27][28]). As far as we know, up to now, there does not exist any study on the a posteriori error estimates and adaptive finite methods for the non-self-adjoint Steklov eigenvalue problem considered in this study. Therefore, deriving the a posteriori error estimates for the problem is still a new topic.
In this study, we first propose error indicators for primal eigenfunction, dual eigenfunction, and eigenvalue of the non-self-adjoint Steklov eigenvalue problem in inverse scattering. The theoretical analysis in this study primarily refers to [11]. But different from what Armentano and Patra discussed in [11], the differential operator corresponding to this problem is non-self-adjoint and the associated weak formulation is not H 1 -elliptic, which leads to the difficulty of theoretical analysis. We use Gårding's inequality and duality technique to solve this difficulty. Second, we give upper and lower bounds for the energy norm of error of finite element solution, which shows that our indicators are reliable and efficient. Finally, based on the error indicators, we design an adaptive algorithm and present numerical examples to show the efficiency of our algorithm. Numerical results coincide with our theoretical analysis.
For the basic theory of finite element, we refer to [24,[29][30][31][32][33]. In this study, without special explanation, C denotes the constant independent of mesh size, which may stand for different values at its different occurrences.

Preliminary
Let ⊂ Ω 2 be a bounded polygonal domain with Lipschitz boundary Γ = ∂Ω. We consider the following Steklov eigenvalue problem:

2)
where ν denotes the unit outward normal to Γ, k is the wavenumber, and n(x) is the index of refraction. Assume that n(x) is a bounded complex value function given by ( ) = ( ) + ( ) n x n x i n x k , 1 2 where = − i 1 , n 1 (x) ≥ δ > 0, and n 2 (x) ≥ 0 are bounded and properly smooth. Here, we use the standard notation for Sobolev spaces H s (Ω) with norm ∥⋅∥ ( ) , be defined as follows: Let π h be a family of triangulations of Ω such that any two triangles in π h share at most a vertex or an edge. h E denotes the diameter of element be the diameter of π h . We assume that the family of triangulations π h satisfies a minimal angle condition, i.e., there exists a constant σ > 0 such that where P 1 denotes the space of linear polynomials. Then, the discrete variational formulation associated with (2.3)-(2.4) is given by: Note that the primal and dual eigenvalues are connected via = λ λ ⁎ . The discrete variational formulation associated with (2.7)-(2.8) is given by: Note that the primal and dual eigenvalues are connected via For convenience of reading, we introduce the following notation. Let λ be the jth eigenvalue of (2.3)-(2.4) with the ascent α = 1 and the algebraic multiplicity q, i.e., λ = λ i = ⋯ = λ i+q−1 . Then, there are q eigenvalues λ j,h (j = i, i + 1, … , i + q − 1) of (2.5)-(2.6) converging to λ. Let M(λ) be the space spanned by all eigenfunctions corresponding to the eigenvalue λ. Let M h (λ) be the space spanned by all eigenfunctions of (2.5)-(2.6) corresponding to the eigenvalues λ j,h (j = i, i + 1, … , i + q − 1). For the dual problems (2.7)-(2.8) and (2.9)-(2.10), the definitions of M*(λ*) and ( ) M λ h ⁎ ⁎ are made similar to M(λ) and M h (λ), respectively. Thanks to Theorem 3.2 in [20] and Theorem 2.1 in [21], the following lemmas are valid. Lemma 2.1. Let M(λ) ⊂ H 1+r (Ω). Assume that λ and λ h are the jth eigenvalue of (2.3)-(2.4) and (2.5)-(2.6), respectively. There exists h 0 ∈ (0,1) such that, if h ≤ h 0 , then where r < π/ω (with ω being the largest inner angle of Ω) and r = 1 if Ω is convex.
Similar conclusions hold for the finite element approximation of dual problem (2.7)-(2.8).
Lemma 2.2. Let M*(λ*) ⊂ H 1+r (Ω). Assume that λ* and λ h ⁎ are the jth eigenvalue of (2.7)-(2.8) and (2.9)-(2.10), respectively. There exists h 0 ∈ (0,1) such that, if h ≤ h 0 , then , then there exists u* ∈ M*(λ*) such that where Ê and lˆis the union of all elements sharing at least a vertex with E and l, respectively. It is easy to verify that a(·,·) satisfies Gårding's inequality [31], i.e., there exist constants 0 < M < ∞ (M is large enough) and α 0 > 0 such that where both M and α 0 are dependent on k and ∥ ( )∥ ( ) ∞ n x L Ω . Third, we will use inequality (3.9) and the arguments in [11] to prove the following theorem. The theorem shows that the global error indicator provides upper bound of the error in the energy norm up to two higher order terms, which guarantees the reliability of indicator.  Noting that e = u − u h , and using the triangle inequality and Schwarz's inequality, we have Using Schwarz's inequality, (3.7) and (3.8), and that the triangulation satisfies the minimum angle condition we obtain Using the trace theorems, we get ⁎ and using the similar arguments to (3.10), we can derive (3.11) easily. Hence

proved. □
Finally, we use the arguments in [11] to prove the efficiency of our indicator for practical adaptive refinement. To achieve our goal, we must prove that the local indicator is bounded by the error in the energy norm and higher order terms. Now, we will prove that the following Lemmas 3.1 and 3.2 are true, which will be used to estimate the first term of η E .
Proof. Let δ 1,E , δ 2,E , and δ 3,E denote the barycentric coordinates of E ∈ π h . Define the cubic bubble b E by and noting that Thus, The adaptive finite element method  223 We estimate the first term on the right-hand side of (3.24) by using Schwarz's inequality and (3.23).
We estimate the second term on the right-hand side of (3.24) by using Schwarz's inequality and (3.22).
The following Lemma 3.2 will be used to estimate the second term of η E .
Lemma 3.2. The following conclusions are valid: be the two triangles sharing l. Then, there exists a constant C such that Proof. Using similar arguments to Lemma 3.2 in [11], we complete the proof. For ∈ ∩ l ε ε E Ω , we denote by E 1 and E 2 the two triangles sharing l and we enumerate the vertices of E 1 and E 2 such that the vertices of l are numbered first. Then, we introduce the edge-bubble function b l Then, such v l exists and is unique.
where ψ i = δ i,E b l and δ 1,E , δ 2,E are the barycentric coordinates associated with the vertices of l. Then, using that ∫ = || Using the similar arguments to (3.22), we get Analogously, combining (3.29) we obtain (b) For l ∈ ε E ∩ ε Ω , let E 1 , E 2 ∈ π h be the two triangles sharing l. Let ∈ ( ∪ ) v H E E l 0 1 1 2 be an edge-bubble function such that v l | E i ∈ P 2 , i = 1,2, and ; . l l l l l L l l L l Using similar arguments to (a), we deduce From (3.34) and Lemma 3.1, we conclude (3.28). □ The following theorem provides the lower bound for the error and high-order terms.
where E u denotes the union of E and the triangles sharing an edge with E.
Proof. It follows from Lemmas 3.1 and 3.

□
There are similar conclusions for the dual problem.
Referring to Lemma 9.1 in [29], we deduce the following lemma.

Edge residual error estimator
In this section, for primal eigenfunction and dual eigenfunction, we propose simpler global error indicators which are equivalent to the errors up to higher order terms, respectively.
Actually, they are obtained by omitting the first term in the residual error indicator given in (3.1) and (3.3), respectively. Next, we prove that they are reliable and efficient up to higher order terms. For any P ∈ N Ω , we define Ω P = ∪{E ∈ π h : P ∈ E}.
Proof. Using similar arguments to Lemma 4.1 in [11], we complete the proof. Let I 0 (u h ) be the L 2 (Ω P ) projection of u h onto the constants, i.e., Next, we estimate the first term on the right-hand side of (4.3). Let Φ P be the corresponding Lagrange basis function with supp Φ P = Ω P , then The relationships (4.4) and (4.5) imply that Inequality (4.1) can be concluded by using the standard estimate for the L 2 (Ω P ) projection. □ The following theorem indicates that the indicator is globally reliable and locally efficient up to higher order terms.
and if ∂ ∩ = ∅ E Γ for E ∈ π h , then Proof. It follows immediately from Theorem 3. (4.10) The adaptive finite element method  229 Proof. It follows immediately from Theorem 3.1, Theorem 3.3, and Lemma 4.1. □ On adaptive mesh with local refinement, the convergence rate of eigenvalue is measured by the total number of degrees of freedom. There are some literature studies that have discussed the relationship between the mesh diameter and the degrees of freedom (e.g., see Remark 2.1 in [35]).

Adaptive algorithm and numerical experiments
First of all, we give the following Algorithm 5.1 which is fundamental and important and can be found in [36].
Step 1. Pick any initial mesh π h0 with mesh size h 0 .
Step 6. Refine π hi to get a new mesh + π hi 1 by Procedure REFINE.
Step 1. Construct a minimal subset π hi of π hi by selecting some elements in π hi such that Step 2. Mark all the elements in π hi .
Marking strategy M was first proposed by [37].
The procedure REFINE is some iterative or recursive bisection (see, e.g., [38,39]) of elements with the minimal refinement condition that marked elements are bisected at least once.
Next, we will report some numerical examples for Algorithm 5.1 to validate our theoretical results. For comparison, using linear finite element, we also solve the problem on uniform meshes. The discrete eigenvalue problems are solved in MATLAB 2016b on an Lenovo IdeaPad 500-14ACZ PC with 1.8 GHZ CPU and 8 GB RAM. Our program is compiled under the package of iFEM [40] and we take θ = 0.25 in our program. In step 2 and step 7, we use the sparse solver "eigs" in MATLAB 2016b to solve the discrete eigenvalue problems. Numerical experiments are considered on the L-shaped domain (−1,1) 2 \([0,1) × (−1,0]) and the square with a slit . The initial mesh is made up of congruent triangles with the mesh size = h 0 2 32 . We select the index of refraction n(x) = 4 or n(x) = 4 + 4i when k = 1 (refer to [7]). When n(x) is the real number, the first 30 eigenvalues are sorted in the descending order. When n(x) is the complex number, the first 30 eigenvalues are sorted in the descending order of the absolute value of imaginary part. We select n(x) = 4 + 2i, 4 + i when k = 2 and k = 4, respectively. For convenience and simplicity, the following notations are introduced in tables and figures.
• i: the ith iteration of Algorithm 5.1.
• λ j h A , i : the jth eigenvalue from the ith iteration of the Adaptive Algorithm 5.1. • λ j,h : the jth eigenvalue in uniform meshes.
• dof: number of degrees of freedom. • η 2 : the a posteriori error estimator of approximate eigenvalue obtained by Algorithm 5.1.
For the L-shaped domain, using the linear finite element, we solve the first three eigenvalues on uniform meshes, the error curves are depicted in Figure 1. The error curves of eigenvalues obtained by Algorithm 5.1 are depicted in Figures 2 and 3. The numerical results are listed in Tables 1 and 2.
We know that the optimal convergence order of numerical eigenvalues on adaptive meshes can achieve O(dof −1 ), and error curve of numerical eigenvalue should be parallel to the line with slope −1. In the case of k = 1, in Figure 1, we see that the error curves of the second eigenvalues are not parallel to the line with slope −1, which implies that the eigenfunction corresponding to the eigenvalue is of singularity and the convergence order of numerical eigenvalue computed on uniform meshes cannot achieve the optimal order. In Figure 2, the error curve of eigenvalue computed on adaptive meshes is parallel to the line with slope −1. Therefore, the numerical eigenvalue obtained by Algorithm 5.1 can achieve the optimal order. Using adaptive meshes is more effective than using uniform meshes when the eigenfunction corresponding to eigenvalue is of singularity. In addition, we know that the a posteriori error indicator of numerical eigenvalue is reliable and efficient as the error curve of numerical eigenvalue is basically parallel to the curve of error indicator. In the case of k ≠ 1, the eigenfunctions corresponding to these eigenvalues are of singularity. In Figure 3, we also see that the numerical eigenvalues computed in adaptive meshes can achieve the optimal convergence order O(dof −1 ).
For the square with a slit, in the case of k = 1, the error curves of the first three eigenvalues are depicted in Figure 4. The error curves of eigenvalues obtained by Algorithm 5.1 are depicted in Figures 5 and 6. The numerical results are listed in Tables 3 and 4. From the figures and tables, we obtain similar conclusion to the L-shaped domain, which coincides in theory results. We also see that the numerical eigenvalues computed in adaptive meshes can achieve the optimal convergence order O(dof −1 ).
The adaptive finite element method  231 Error of i-th eigenvalues The line with slope -1       Finally, adaptively refined meshes obtained by Algorithm 5.1 are displayed in Figure 7 as k = 1 and n(x) = 4 + 4i. In theory, the errors of eigenfunctions on the elements around concave vertices should be greater than other elements. In Figure 7, we see that, on each domain, these elements around concave vertices are more refined than other elements, which coincides in theoretical results.