An efficient spectral-Galerkin method for a new Steklov eigenvalue problem in inverse scattering

1.   Introduction

Steklov eigenvalue problems have significant physical background and wide applications in many fields of science and engineering^[1,2,3,4]. Its theoretical analysis and numerical calculation have attracted the attention of many scholars, and a variety of finite element methods and spectral methods have been proposed ^{[5,6,7,8,9,10,11,12,13]}, but these numerical methods are mainly based on the selfadjoint Steklov eigenvalue problem.

Recently, a new Steklov eigenvalue problem arising from inverse scattering attracts many researchers' interest. The corresponding weak formulation of the problem is non-selfadjoint and indefinite. We can not use Lax-Milgram theorem to determine the existence and uniqueness of the solution of the corresponding resource problems, which brings some difficulties to deduce the equivalent operator form of the eigenvalue problem. In addition, in order to prove the error estimation of approximated eigenvalues and eigenfunctions, we need not only to introduce conjugate eigenvalue problem and corresponding conjugate operators, but also to prove the error estimation between these operators and their approximated operators. However, there are still some numerical methods to solve the problem. For instance, Cakoni et al. discussed the conforming finite element approximation in ^[14], Liu et al. studied spectral indicator method in ^[15], Bi et al. discussed two-grid discretizations and a local finite element scheme in ^[16], Zhang et al. established a multigrid correction scheme in ^[17], Yang et al. used non-conforming Crouzeix-Raviart element solve it ^[18], Xu et al. discussed an asymptotically exact a posteriori error estimator for non-selfadjoint Steklov eigenvalue problem ^[19], Wang et al. studied a priori and a posteriori error estimates for a virtual element method for the non-self-adjoint Steklov eigenvalue problem ^[20]. To our knowledge, the study of spectral method for solving this problem has not been reported. Since the spectral methods have the characteristics of spectral accuracy ^{[21,22,23,24,25]}, that is to say, we only need to spend less degrees of freedom to obtain higher accurate numerical solutions. Then, it is meaningful to propose an effective spectral method for solving a new Steklov eigenvalue problem in inverse scattering.

Hence, we shall study an effective spectral-Galerkin method for the new Steklov eigenvalue problem. Firstly, we establish the weak formulation and the associated discrete scheme by introducing an appropriate Sobolev space and a corresponding approximation space. Then, according to the Fredholm Alternative, the corresponding operator forms of weak formulation and discrete formulation are derived. After that, the error estimates of approximated eigenvalues and eigenfunctions are proved by using the spectral theory of compact operators and the approximation properties of orthogonal projection operators. We also construct an appropriate set of basis functions in the approximation space and derive the matrix form of the discrete scheme based on tensor product. In addition, we extend the algorithm to the circular domain and transform the original problem into an equivalent form under the polar coordinates. By using orthogonal polynomial approximation in the radial direction and Fourier basis function approximation in the $\theta$ direction, combining with pole conditions, we construct an appropriate approximation space and derive the corresponding matrix form of the discrete scheme. Finally, we provide plenty of numerical examples and compare them with some existing numerical methods, the numerical results confirm the effectiveness and high accuracy of our algorithm.

We derive the weak formulation and the corresponding discrete scheme and prove the error estimations of approximation eigenvalues and eigenfunctions in next section. In §3, we propose an efficient algorithm to solve the Steklov eigenvalue problem in the square domain. In §4, we extend our algorithm to the case of circular domain. In §5, we provide several numerical examples to validate the accuracy and efficiency of our algorithm. In §6, some concluding remarks are presented.

Throughout this article, a notation $a \lesssim b$ is used to mean that $a \leqslant Cb$, where $C$ is a positive constant independent of any function or any discretization parameters.

2.   The weak formulation, discrete scheme and the error estimates

2.1.   The weak formulation and discrete scheme

Denote by $H^{s}(D)$ and $H^{s}(\partial D)$ the usual Sobolev spaces with integer order $s$ in $D$ and on $\partial D$, respectively. In particular, $H^{0}(D) = L^2(D)$ and $H^{0}(\partial D) = L^2(\partial D)$. The norm in $H^{s}(D)$ and $H^{s}(\partial D)$ are expressed as $\|\psi\|_{s, D}$ and $\|\psi\|_{s, \partial D}$, separately.

Consider the following Steklov eigenvalue problem:

where $D \subset \mathbb{R}^{2}$ (or $\mathbb{R}^{3}$) is a bounded polygon with Lipschitz boundary $\partial D$, $\nu$ is the unit outward normal on $\partial D$. Let $k$ be the wavenumber and $\beta(x) = \beta_{1}(x)+ i \frac{\beta_{2}(x)}{k}$ be the index of refraction that is a bounded complex value function with $\beta_{1}(x) > 0$ and $\beta_{2}(x)\geqslant 0$.

The weak formulation of the Eqs (2.1) and (2.2) is to find $(\lambda, \psi)\in \mathbb{C }\times H^{1}(D)$, $\psi \neq 0$ such that

where

Refering to ^[15], we know that $\mathcal{A}(\cdot, \cdot)$ satisfies Garding's inequality, i.e., there exist constants $K < \infty$ and $\alpha_{0} > 0$ such that

Let $K$ be a positive constant which is large enough, and the sesquilinear form is defined as follow:

then it is easy to verify that $\tilde{\mathcal{A}}$ is $H^{1}(D)$-elliptic (see ^[15]).

We first focus on the case of $D = I^{d} (d = 2, 3)$ with $I = (-1, 1)$. Denote by $L_{p}(t)$ the Legendre polynomial of degree $p$. Let $V_{N} = \text{span} \{ L_{0}(t), L_{1}(t), \cdots, L_{N}(t) \}$, then the approximation space $X_{N} = (V_{N})^{d}.$

The spectral-Galerkin approximation for the eigenvalue problem (2.3) reads: Find $(\lambda_{N}, \psi_{N})\in \mathbb{C}\times X_{N}$, $\psi_{N} \neq 0$ such that

2.2.   Error estimation

We first consider the following source problem associated with (2.3): Given $g\in H^{-\frac{1}{2}}(\partial D)$, find $w\in H^{1}(D)$ such that

and the approximation source problem associated with (2.4): Find $w_{N} \in X_{N}$ such that

Further, we introduce Neumann eigenvalue problem as follows:

When $k^{2}$ is not a Neumann eigenvalue of (2.7) and (2.8), from the Fredholm Alternative (see, e.g., Section 5.3 of ^[26] or Lemma 1 in ^[18]), we know that for $g \in H^{-\frac{1}{2}}(\partial D)$, there exists a unique solution $w \in H^{1}(D)$ of (2.5) such that

Define an operator $A: H^{-\frac{1}{2}}(\partial D)\rightarrow H^{1}(D)$ by

and a Neumann-to-Dirichlet mapping $\Gamma:H^{-\frac{1}{2}}(\partial D)\rightarrow H^{\frac{1}{2}}(\partial D)$ by

We can similarly define a discrete operator $A_{N}: H^{-\frac{1}{2}}(\partial D)\rightarrow X_{N}$ such that

and a discrete Neumann-to-Dirichlet operator $\Gamma_{N}: H^{-\frac{1}{2}}(\partial D)\rightarrow X_{N}^{B}$ satisfies

where $X_{N}^{B} = X_{N}|_{\partial D}$. Then from (2.9), we obtain $\|Ag\|_{1, D} \lesssim \|g\|_{-\frac{1}{2}, \partial D}$.

Thus, the equivalent operator forms of (2.3) and (2.4) are:

In this paper, we always assume that $k^{2}$ is not a Neumann eigenvalue of (2.7) and (2.8).

Consider the dual problem of (2.3): Find $(\lambda^{*}, \psi^{*})\in\mathbb{C}\times H^{1}(D)$, $\psi^{*} \neq 0$ such that

It's clearly that the primal and dual eigenvalues satisfy the relation: $\lambda = \overline{\lambda^{*}}$.

The discrete variational formulation associated with (2.14) is given by: Find $(\lambda_{N}^{*}, \psi_{N}^{*})\in\mathbb{C}\times X_{N}$, $\psi_{N}^{*} \neq 0$ such that

Likewise, the primal and dual eigenvalues satisfy the relation: $\lambda_{N} = \overline{\lambda^{*}_{N}}$.

Similarly, define the dual operators $A^{*}: H^{-\frac{1}{2}}(\partial D)\rightarrow H^{1}(D)$ and $A^{*}_{N}: H^{-\frac{1}{2}}(\partial D)\rightarrow X_{N}$ such that

The corresponding Neumann-to-Dirichlet and discrete Neumann-to-Dirichlet dual operators are defined as follows:

Define the $H^{1}$ projection operators $\Pi_{N}^{1}: H^{1}(D)\rightarrow X_{N}$ and ${\Pi_{N}^{1}}^{*}: H^{1}(D)\rightarrow X_{N}$ by

Then for any $g\in H^{-\frac{1}{2}}(D)$, we have

Since the above equation admits a unique solution, we have $A_{N} = \Pi_{N}^{1}A$. Similarly, we can obtain $A_{N}^{*} = {\Pi_{N}^{1}}^{*}A^{*}$.

There holds the following regular results, which will be used in the following theoretical analysis.

Lemma 1. If $g \in L^{2}(\partial D)$, then $Ag \in H^{1+\frac{\kappa}{2}}(D)$ and

if $g\in H^{\frac{1}{2}}(\partial D)$, then $Ag \in H^{1+\kappa}(D)$ and

where $\kappa = 1$ when the largest inner angle $\theta$ of $D$ satisfying $\theta < \pi$, and $\kappa < \frac{\pi}{\theta}$ which can be arbitrarily close to $\frac{\pi}{\theta}$ when $\theta > \pi$.

Proof. see ^[27].

For the dual problem, we have the same regular results as the corresponding source problem. Defining an $H^{1}$-projection operator $P_{N}^{1}: H^{1}(D)\rightarrow X_{N}$ by

According to Theorem 8.4 in ^[28], we have the following lemma:

Lemma 2. For any $w \in H^{s}(D)$ with $s \geq 1$ and $s \geqslant l \geqslant 0$,

Denote

It is obvious that

Lemma 3. Let $w$ be the solution of $(2.5)$, if $w\in H^{s}(D) (s \geq 2)$, there hold:

Proof. According to the definition of the projection operator $\Pi_{N}^{1}$, we obtain

From Theorem 3.1 in ^[15] and the inequality (2.22), we derive that

Similarly, we can arrive at

From Lemma 2.2 in ^[16], for any $\phi_{N} \in X_{N}$ we derive that

Taking $\phi_{N} = \Pi_{N}^{1} Ag$ and using (2.21) and (2.23) we obtain that

The proof is completed.

Consider the following auxiliary problem: Find $\xi_{f}\in H^{1}(D)$, such that

Referring to ^[16], we obtain the following result.

Lemma 4. If $f \in L^{2}(D)$, then there exists a unique solution $\xi_{f}\in H^{1+\kappa}(D)$ to $(2.26)$ and

where the principle to determine $\kappa$ see Lemma 1.

Remark 1. If the regularity results, Lemmas 1 and 4 hold in the case of $D \in R^{3}$, we can prove the analysis and conclusion in this paper for $D \in R^{3}$. However, there are no such good results in $R^{3}$ (see, e.g., Remark 2.1 in ^[29] and Remark 1 in ^[18]). Our analysis is still valid in the case of $D \in R^{3}$, but the conclusions need minor modifications.

Lemma 5. For $\forall w\in H^{1}(D)$, we have

Proof. According to the definition of the projection ${\Pi_{N}^{1}}^{*}$, for any $\xi \in H^{s}(D) (s \geq 2)$ we have

From Lemma 2.4 in ^[16], we can obtain that

From the inequalities (2.24) and (2.27), we have

Thus, we obtain

This ends our proof.

According to the classical spectral approximation theory: when $\lim\limits_{N \rightarrow \infty} \|\Gamma-\Gamma_{N}\|_{H^{-\frac{1}{2}}(\partial D) \longrightarrow H^{-\frac{1}{2}}(\partial D)} = 0$, we can obtain the error estimates of eigenvalue problem.

Lemma 6. If $s \geq 2$, then $\lim\limits_{N \rightarrow \infty} \|\Gamma-\Gamma_{N}\|_{H^{-\frac{1}{2}}(\partial D) \longrightarrow H^{-\frac{1}{2}}(\partial D)} = 0$ and $\Gamma$ is a compact operator.

Proof. From Lemma 2.6 in ^[16], for any $\phi_{N} \in X_{N}$ we have

Taking $\phi_{N} = \Pi_{N}^{1} Ag$, from (2.23) and (2.21) we derive that

From the above inequality, together with (2.23) and (2.9) we obtain that

Since $s \geqslant 2$, then we obtain $\lim\limits_{N \rightarrow \infty} \|\Gamma-\Gamma_{N}\|_{H^{-\frac{1}{2}}(\partial D) \longrightarrow H^{-\frac{1}{2}}(\partial D)} = 0$. Note that $\Gamma_{N}$ is a finite rank operator. Then, $\Gamma$ is a compact operator. This ends our proof.

Let $\lambda$ be the $p$-th eigenvalue of (2.3) with the algebraic multiplicity $h$ and the ascent $\alpha$. Since $\Gamma_{N}$ converges to $\Gamma$, there are $h$ eigenvalues $\lambda_{q, N}(q = \delta, \delta+1, \delta+2, \ldots, \delta+h-1)$ of (2.4) converging to $\lambda$. Let $M(\lambda)$ be the generalized eigenfunction space of (2.3) associated with the eigenvalue $\lambda$, and $M_{N}(\lambda)$ be the generalized eigenfunction space of (2.4) associated with the eigenvalue $\lambda_{q, N}\; (q = \delta, \delta+1, \delta+2, \ldots, \delta+h-1)$. As for the dual problems (2.14) and (2.15), we can also define $M^{*}(\lambda^{*})$ and $M_{N}^{*}(\lambda^{*})$.

Theorem 1. Assume that $\psi_{N}$ is the eigenfunction approximation of $(2.4)$, then there exists an eigenfunction $\psi \in M(\lambda)$ of $(2.3)$ corresponding $\lambda$ such that

and

Proof. From the Theorem 7.3 and Theorem 7.4 in ^[30] we know that there exists an eigenfunction $\psi$ corresponding to $\lambda$ and

where $\varphi_{\delta}, \ldots, \varphi_{\delta+h-1}$ are the any basis for $M(\lambda)$ and $\varphi_{\delta}^{*}, \ldots, \varphi_{\delta+h-1}^{*}$ are the dual basis in $M^{*}(\lambda^{*})$.

From Theorem 2.1 in ^[16], combining (2.23) and (2.24) we derive that

We derive from Lemma 6, (2.22) and (2.23) that

Similarly, we have

Combining the inequalities (2.34) and (2.36), we obtain (2.29). From (2.33), (2.35), (2.36) and (2.37), we obtain (2.32).

From Theorem 2.1 in ^[16], and the inequalities (2.9), (2.32), (2.29) and (2.23), we deduce

From Theorem 2.1 in ^[16], and the inequalities (2.28), (2.9), (2.32), (2.29) and (2.23), we obtain

This ends our proof.

For the dual problems (2.14) and (2.15), we can also obtain the corresponding conclusion as Theorem 1.

3.   Efficient implementation of the algorithm

In this section, we shall present an efficient algorithm to solve the discrete scheme (2.4). We first construct a group of basis functions of the approximation space $X_{N}$. Denote by $L_p(t)$ the Legendre polynomial of degree $p$. Let $\hat{\varphi}_{l}(t) = L_l(t)-L_{l+2}(t), l = 0, 1, \cdots, N-2$, $\hat{\varphi}_{N-1}(t) = L_{0}(t), \hat{\varphi}_{N}(t) = L_{1}(t)$. Then the approximation space $X_{N} = \text{span} \{\hat{\varphi}_{l}(x_{1}) \hat{\varphi}_{j}(x_{2}): l, j = 0, 1, \ldots, N \}$.

Let $a_{l j} = \int_{-1}^{1}\hat{\varphi}_j'\hat{\varphi}_{l}'dt$, $b_{{l} j} = \int_{-1}^{1}\hat{\varphi}_j\hat{\varphi}_{l}dt$, $c_{{l} j} = \hat{\varphi}_j(-1)\hat{\varphi}_{l}(-1)$, $d_{l j} = \hat{\varphi}_j(1)\hat{\varphi}_l (1)$. By utilizing the orthogonal properties of the Legendre polynomials, we obtain that

(1)

where $l, j = 0, 1, \ldots, N-2$.

(2)

where $l = N-1, N, j = 0, 1, \ldots, N-2$.

(3)

where $l = N-1, N, j = N-1, N$.

Next, we shall derive the matrix form based on the tensor product of the discrete scheme (2.4). For simplicity, we only consider the case of $d = 2$. It can be derived similarly for the case of $d = 3$.

Let $\psi_{N}(x_{1}, x_{2}) = \sum\limits_{l, j = 0}^{N} \psi_{lj} \hat{\varphi}_{l}(x_{1})\hat{\varphi}_{j}(x_{2})$ and

We denote by $\overline{\Psi}$ a column vectors with $(N+1)^{2}$ elements, which is consist of the $N+1$ columns of $\Psi$. Let $\phi_{N}(x_{1}, x_{2}) = \hat{\varphi}_{m}(x_{1})\hat{\varphi}_{n}(x_{2}), m, n = 0, 1, \ldots, N$, then we derive that

where $\hat{A} = (a_{lj})_{l, j = 0}^{N}$, $\hat{B} = (b_{lj})_{l, j = 0}^{N}$, $\hat{C} = (c_{lj})_{l, j = 0}^{N}$, $\hat{D} = (d_{lj})_{l, j = 0}^{N}$. $\hat{A}(m, :)$ indicates the m-th row of the matrix $\hat{A}$, $\hat{B}(m, :)$, $\hat{C}(m, :)$, $\hat{D}(m, :)$ are similar to $\hat{A}(m, :)$. $\otimes$ is a notation of tensor product of matrix, then $\hat{A} \otimes \hat{B} = (a_{lj} \hat{B})_{l, j = 0}^{N}$.

Let $\xi_{\mu}, w_{\mu} (\mu = 0, 1, \ldots, N_{1})$ be the Labatto-Gauss points and weights, respectively. Then we have

where

From the above deduce, we obtain the matrix form of the discrete scheme (2.4) as follows:

Note that the matrix $\hat{P}$ can also be written as matrix form base on the tensor product and the stiffness matrix and mass matrix in (3.1) are all sparse when $\beta$ is a constant, so we can solve (3.1) efficiently.

4.   Extension to circular domain

We extend the above algorithm to a circular domain. Utilizing polar coordinate transformation $x = (x_{1}, x_{2}) = (r \text{cos}\theta, r \text{sin}\theta)$, the functions $\psi(x)$ and $\beta(x)$ are represented as $\widetilde{\psi}(r, \theta) = \psi (r \text{cos}\theta, r \text{cos}\theta)$ and $\tilde{\beta}(r, \theta) = \beta(r \text{cos}\theta, r \text{cos}\theta)$. Let $Lv = \frac{1}{r} \frac{\partial}{\partial r}(r\frac{\partial v}{\partial r})+\frac{1}{r^{2}} \frac{\partial^2 v}{\partial \theta^2}$, then the Eqs (2.1) and (2.2) can be rewritten as follows:

Then from (2.3) we derive that the weak form of (4.1) and (4.2) is as follows:

where

Since $\widetilde{\psi}$ is $2\pi$ periodic in $\theta$, by using the expansion of the Fourier basis function we have

Substituting the expansion (4.4) into (4.3), and taking $\tilde{\phi}(r, \theta) = v_{\tilde{n}}(r)e^{i\tilde{n}\theta}$, we derive that

In order to make (4.5) well posed, we need to introduce the following polar condition (see ^[22])

Let $r = \frac{R}{2}(t+1)$, $t\in (-1, 1)$, $w_{\tilde{m}}(t) = u_{\tilde{m}}(r)$, $z_{\tilde{n}}(t) = v_{\tilde{n}}(r)$, $\gamma(t, \theta) = \widetilde{\beta}(r, \theta)$, then (4.5) can be rewritten as follows:

Define the weighted Sobolev space:

which is endowed with the norm as follows:

Then the weak formulation of (4.5) is to find $w_{\tilde{m}}(t)e^{i\tilde{m} \theta}\in \mathcal{H}_{\tilde{m}}^{1}(\Omega)$, $\lambda \in \mathbb{C}$ such that

for $z_{\tilde{n}}(t)e^{-i\tilde{n} \theta}\in \mathcal{H}_{\tilde{n}}^{1}(\Omega)$.

Denote by $P_{N}$ the space of polynomials of degree less than or equal to $N$, and set $X_{N}^{M} = \bigoplus\limits_{|\tilde{m}| = 0}^{M} \mathcal{H}_{\tilde{m}}^{1}(\Omega)\bigcap P_{\tilde{m}N},$ where $P_{\tilde{m}N} = \{p_{N}e^{i\tilde{m}\theta}: p_{N}\in P_{N}\}$.

Then the discrete form of (4.8) is to find $w_{\tilde{m}N}e^{i\tilde{m} \theta}\in X_{N}^{M}$, $\lambda_{N} \in \mathbb{C}$ such that

for $\forall\; z_{\tilde{n}N}e^{-i\tilde{n} \theta}\in X_{N}^{M}$.

4.1.   Implementation of the algorithm

Let

It is clear that

Setting

Let $w_{\tilde{m}N} = \sum\limits_{\tilde{i} = 0}^{N-sign(|\tilde{m}|)} w_{\tilde{i}}^{\tilde{m}}\varphi_{\tilde{i}}^{\tilde{m}}(|\tilde{m}| = 0, \ldots, M)$, $z_{\tilde{n}N} = \varphi_{\tilde{j}}^{\tilde{n}}(|\tilde{n}| = 0, \ldots, M, \tilde{j} = 0, \ldots, N-sign(|\tilde{n}|))$, and inserting the expressions into (4.9), we obtain that

Then the corresponding matrix form of (4.10) can be derived as follows:

where

$\overline{W}$ is a column vector with $N(2M+1)+1$ elements, which is consist of the column of $W$. By using the orthogonality of the Fourier system ${e^{i\tilde{m} \theta}}$, we obtain that

It's obvious that the matrices F, G, and H are diagonal block matrices.

Let $t_{\tilde{s}}, \; w_{\tilde{s}} (\tilde{s} = 0, 1, \ldots, N_{1})$ be the Gauss-Lobatto points and the weights and $\theta_{\tilde{l}}, \; \hat{\omega}_{\tilde{l}} (\tilde{l} = 0, 1, \ldots, \; M_{1}-1)$ be the Fourier points and the weights, respectively. Then we have

5.   Numerical experiments

We shall present some numerical experiments to validate the effectiveness of the algorithm. We carry out our programs in Matlab 2018a.

5.1.   Square domain

Example 1. Consider the problem (2.3) with $k = 1$ and $\beta(x) = 4$ on the domain $D = (-\frac{\sqrt{2}}{2}, \frac{\sqrt{2}}{2})^{2}$. The approximation eigenvalue $\lambda_{N}^{i} (i = 1, 2, 3, 4)$ for different $N$ are shown in Table 1.

We observe from Table 1 that the approximation eigenvalues reach at least fourteen-digit accuracy with $N \geq 20$. For comparison, we list in Table 2 the numerical results obtained by Multigrid Correction Scheme in ^[17]. However, the numerical results reported in Table 2 have at most six-digit accuracy despite utilizing a great quantity of degrees of freedom.

Next, we choose the numerical solutions of $N = 40$ as reference solutions, the corresponding error figures of the approximate eigenvalues $\lambda^{i}_{N} (i = 1, 2, 3, 4)$ with different $N$ are listed in Figure 1. We know from Figure 1 that the approximation eigenvalues converge gradually with the increase of $N$.

Example 2. We take $k = 1$, $\beta(x) = 4+4i$ and $D = (-\frac{\sqrt{2}}{2}, \frac{\sqrt{2}}{2})^{2}$ as our second example. The approximation eigenvalue $\lambda_{N}^{i} (i = 1, 2, 3, 4)$ of complex Steklov eigenvalues with largest imaginary parts on the square domain are shown in Table 3.

We can see from Table 3 that the complex Steklov eigenvalues reach at least thirteen-digit accuracy with $N \geq 20$. The numerical solutions obtained by Multigrid Correction Scheme of ^[17] in Table 4 have at most seven-digit accuracy despite utilizing a great quantity of degrees of freedom.

Similarly, we also choose the numerical solutions of $N = 40$ as reference solutions, the corresponding error figures of the approximate eigenvalues $\lambda^{i}_{N} (i = 1, 2, 3, 4)$ with different $N$ are listed in Figure 2. We observe from Figure 2 that the approximation eigenvalues also converge gradually with the increase of $N$.

Example 3. When $\beta(x)$ is a variable coefficient, we take $k = 1$, $\beta(x_{1}, x_{2}) = [(x_{1}+x_{2})^{2}+1]+(x_{1}-x_{2})^{2}i$ and $D = (-\frac{\sqrt{2}}{2}, \frac{\sqrt{2}}{2})^{2}$. The approximation eigenvalue $\lambda_{N}^{i} (i = 1, 2, 3, 4)$ of complex Steklov eigenvalues with largest imaginary parts on the square domain are shown in Table 5.

We can see from Table 5 that the first four of complex Steklov eigenvalues reach at least thirteen-digit accuracy with $N \geq 20$. Again, we choose the numerical solutions of $N = 40$ as reference solutions, the corresponding error figures of the approximate eigenvalues $\lambda^{i}_{N} (i = 1, 2, 3, 4)$ with different $N$ are listed in Figure 3. From Figure 3, we can see that the approximation eigenvalues also converge gradually with the increase of $N$.

Example 4. We consider the problem (2.3) in a three-dimensional case, where we take $k = 2$, $\beta(x) = 1+i$ and $\overline{D} = [0, 1]^{3}$. The approximation eigenvalue $\lambda_{N}^{i} (i = 1, 2, 3, 4)$ of complex Steklov eigenvalues with largest imaginary parts on $\overline{D} = [0, 1]^{3}$ are shown in Table 6.

We can see from Table 6 that the first four of complex Steklov eigenvalues reach at least thirteen-digit accuracy with $N \geq 15$.

5.2.   Circular disk

Example 5. Consider the problem (2.3) with $k = 1$ in the unit disk centered at $(0, 0)$ with radius $R = 1$. We choose the index of refraction $\beta(x) = 4$. The three largest Steklov eigenvalues for different $N$ and different $M$ are listed in Tables 7–9, respectively.

We observe from Tables 7–9 that approximation eigenvalues to $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$ reach at least fourteen-digit accuracy with $N \geq 15$ and $M \geq 4$. For comparison, we list the numerical results of ^[15] in Table 10. The numerical results to $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$ reported in Table 10 have at most three-digit accuracy despite utilizing a great quantity of degrees of freedom.

Example 6. When $\beta(x)$ is complex, we take $\beta(x) = 4+4i$, $k = 1$ and $\Omega$ is a unit disk centered at $(0, 0)$ with radius $R = 1$. The numerical results of the first three complex Steklov eigenvalues with largest imaginary parts in a unit circle are shown in Tables 11–13, respectively.

We observe from Tables 11–13 that approximation eigenvalues to $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$ achieve at least thirteen-digit accuracy with $N \geq 15$ and $M \geq 4$. The results in Table 14 are obtained in ^[15], which are comparison to the results in Tables 11–13. The numerical eigenvalues to $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$ reported in Table 14 have at most four-digit accuracy despite utilizing a great quantity of degrees of freedom.

6.   Conclusions

In this paper, we propose an efficient spectral method for solving a new Steklov eigenvalue problem. First, we give an efficient spectral Galerkin approximation for the new Steklov eigenvalue problem in rectangular domain, and prove the error estimates of approximation eigenvalues and eigenfunctions. Secondly, we derive the matrix form of discrete scheme based on tensor product, and analyze the sparsity of mass matrix and stiffness matrix. In addition, we give an efficient spectral Galerkin approximation for the new Steklov eigenvalue problem in circular domain. Finally, we present ample numerical examples which validate the effectiveness and high accuracy of the algorithm.

The method proposed in this paper can be extended to some more complex problems, such as transmission eigenvalue problem, electromagnetic eigenvalue problem, nonlinear eigenvalue problem and so on, which will be our goal in the future.

Acknowledgments

The authors would like to thank the editor and the referees for helpful comments and suggestions. This work is supported by the National Natural Science Foundation of China (No. 11961009), Guizhou Provincial Graduate Education Innovation Program (No. YJSCXJH [2020] 097) and the Scientific Research Foundation of Guizhou University of Finance and Economics(No. 2020XYB10), the Project for Young Talents Growth of Guizhou Provincial Department of Education under (Grant No. KY[2022]179).

Conflict of interest

The authors declare that they have no competing interests.