The inverse electromagnetic scattering problem by a penetrable cylinder at oblique incidence

TL;DR

This paper presents a numerical method using non-linear boundary integral equations to solve the inverse electromagnetic scattering problem for a penetrable cylinder at oblique incidence, demonstrating its feasibility through numerical results.

Contribution

It introduces a novel iterative scheme that linearizes only the far-field equation for efficient inverse scattering reconstruction.

Findings

Numerical results confirm the method's effectiveness.

The approach accurately reconstructs the cylinder's properties.

The scheme is computationally feasible for 3D problems.

Abstract

In this work we consider the method of non-linear boundary integral equation for solving numerically the inverse scattering problem of obliquely incident electromagnetic waves by a penetrable homogeneous cylinder in three dimensions. We consider the indirect method and simple representations for the electric and the magnetic fields in order to derive a system of five integral equations, four on the boundary of the cylinder and one on the unit circle where we measure the far-field pattern of the scattered wave. We solve the system iteratively by linearizing only the far-field equation. Numerical results illustrate the feasibility of the proposed scheme.

Equations

$$\nabla \times {\text{E}}^{ext}+{\mu }_0\frac{\partial {\text{H}}^{ext}}{\partial t}$$

$$\nabla \times {\text{E}}^{int}+{\mu }_1\frac{\partial {\text{H}}^{int}}{\partial t}$$

$$\hat{\mathbf{n}}\times {\text{E}}^{int}=\hat{\mathbf{n}}\times {\text{E}}^{ext},\hat{\mathbf{n}}\times {\text{H}}^{int}=\hat{\mathbf{n}}\times {\text{H}}^{ext},\text{x}\in \Gamma ,$$

$${\text{E}}^{inc}(\text{x};\hat{\mathbf{d}},\hat{\mathbf{p}})={\text{e}}^{inc}(\mathbf{x})e^{-i\beta z},{\text{H}}^{inc}(\text{x};\hat{\mathbf{d}},\hat{\mathbf{p}})={\text{h}}^{inc}(\mathbf{x})e^{-i\beta z},$$

$${\text{e}}^{inc}(\mathbf{x})$$

$$\Delta e_3^{int}+{\kappa }_1^2{e}_3^{int}=0,\Delta h_3^{int}+{\kappa }_1^2{h}_3^{int}=0,\mathbf{x}$$

$$\Delta e_3^{ext}+{\kappa }_0^2{e}_3^{ext}=0,\Delta h_3^{ext}+{\kappa }_0^2{h}_3^{ext}=0,\mathbf{x}$$

$$e_3^{int}$$

$${\tilde{\mu }}_1\omega \frac{\partial h_3^{int}}{\partial n}+{\beta }_1\frac{\partial e_3^{int}}{\partial \tau }$$

$$h_3^{int}$$

$${\widetilde{ϵ}}_1\omega \frac{\partial e_3^{int}}{\partial n}-{\beta }_1\frac{\partial h_3^{int}}{\partial \tau }$$

$$e_3^{inc}(\mathbf{x})=\frac{1}{\sqrt{{ϵ}_0}}\sin \theta e^{i{\kappa }_0(x\cos \varphi +y\sin \varphi )},h_3^{inc}(\mathbf{x})=0.$$

$$\underset{r\to \infty }{\lim }\sqrt{r}\left(\frac{\partial e_3^{sc}}{\partial r}-i{\kappa }_0{e}_3^{sc}\right)=0,\underset{r\to \infty }{\lim }\sqrt{r}\left(\frac{\partial h_3^{sc}}{\partial r}-i{\kappa }_0{h}_3^{sc}\right)=0,$$

$$({\mathcal{S}}_{j}f)(\mathbf{x})$$

$$({\mathcal{D}}_{j}f)(\mathbf{x})$$

$${\Phi }_j(\mathbf{x},\mathbf{y})=\frac{i}{4}{H}_0^{(1)}({\kappa }_j|\mathbf{x}-\mathbf{y}|),\mathbf{x},\mathbf{y}\in {\Omega }_j,\mathbf{x}\ne \mathbf{y},$$

$$(S_{j}f)(\mathbf{x})$$

$$(D_{j}f)(\mathbf{x})$$

$$(N{S}_{j}f)(\mathbf{x})$$

$$(N{D}_{j}f)(\mathbf{x})$$

$$(T{S}_{j}f)(\mathbf{x})$$

$$(T{D}_{j}f)(\mathbf{x})$$

$$e_3^{int}(\mathbf{x})$$

$$h_3^{int}(\mathbf{x})$$

$$e_3^{ext}(\mathbf{x})$$

$$h_3^{ext}(\mathbf{x})$$

$$({\text{D}}_0+{\text{K}}_0)\left(\begin{array}{c}{\varphi }_0\\ {\psi }_0\end{array}\right)={\text{b}}_0,$$

$${\text{D}}_0$$

$${\text{b}}_0$$

$$e_3^{sc}(\mathbf{x})=\frac{e^{i{\kappa }_{0}r}}{\sqrt{r}}{e}^{\infty }(\hat{\mathbf{x}})+\mathcal{O}(r^{-3/2}),h_3^{sc}(\mathbf{x})=\frac{e^{i{\kappa }_{0}r}}{\sqrt{r}}{h}^{\infty }(\hat{\mathbf{x}})+\mathcal{O}(r^{-3/2}),$$

$$e^{\infty }(\hat{\mathbf{x}})$$

$$h^{\infty }(\hat{\mathbf{x}})$$

$$(S^{\infty }f)(\hat{\mathbf{x}})$$

$$(D^{\infty }f)(\hat{\mathbf{x}})$$

$${\Phi }^{\infty }(\hat{\mathbf{x}},\mathbf{y})=\frac{e^{i\pi /4}}{\sqrt{8\pi {\kappa }_0}}{e}^{-i{\kappa }_0\hat{\mathbf{x}}\cdot \mathbf{y}}.$$

$$e_3^{int}(\mathbf{x})$$

$$e_3^{sc}(\mathbf{x})$$

$$\frac{1}{{\widetilde{ϵ}}_1}\left(D_1-\frac{1}{2}\right){\varphi }_e-\frac{1}{{\widetilde{ϵ}}_0}{S}_0{\psi }_e$$

$$\omega N{D}_1{\varphi }_h+\frac{{\beta }_1}{{\widetilde{ϵ}}_1}\left(T{D}_1-\frac{1}{2}\frac{\partial }{\partial \tau }\right){\varphi }_e-\omega \left(N{S}_0-\frac{1}{2}\right){\psi }_h-\frac{{\beta }_0}{{\widetilde{ϵ}}_0}T{S}_0{\psi }_e$$

$$\frac{1}{{\tilde{\mu }}_1}\left(D_1-\frac{1}{2}\right){\varphi }_h-\frac{1}{{\tilde{\mu }}_0}{S}_0{\psi }_h$$

Full text

Keywords inverse electromagnetic scattering, oblique incidence, integral equation method

1 Introduction

The inverse obstacle scattering problem is to image the scattering object, i.e. find its shape and location, from the knowledge of the far-field pattern of the scattered wave. The medium is illuminated by light at given direction and polarization. Then, Maxwell’s equations are used to model the propagation of the light through the medium, see [6, 7] for an overview. Due to the complexity of the combined system of equations for the electric and the magnetic fields, it is common to impose additional assumptions on the incident illumination and the nature of the scatterer.

We consider time-harmonic incident electromagnetic plane wave that due to the linearity of the problem will result to a time-independent system of equations. In addition, the penetrable object is considered as an infinitely long homogeneous cylinder. Then, it is characterized by constant permittivity and permeability. The problem is further simplified if we impose oblique incidence for the incident wave.

The three-dimensional scattering problem modeled by Maxwell’s equations is then equivalent to a pair of two-dimensional Helmholtz equations for two scalar fields (the third components of the electric and the magnetic fields). This approach reduces the difficulty of the problem but results to more complicated boundary conditions. The transmission conditions now contain also the tangential derivatives of the electric and magnetic fields. In [10] we showed that the corresponding direct problem is well-posed and we constructed a unique solution using the direct integral equation method. A similar problem has been considered for an impedance cylinder embedded in a homogeneous [24], and in an inhomogeneous medium [22]. A numerical solution of the direct problem has been also proposed using the finite element method [4], the Galerkin method [20], and the method of auxiliary sources [23].

On the other hand, the inverse problem is non-linear and ill-posed. The non-linearity is due to the dependence of the solution of the scattering problem on the unknown boundary curve. The smoothness of the mapping from the boundary to the far-field pattern reflects the ill-posedness of the inverse problem. The unique solvability of the inverse problem is still an open problem. The first and only, to our knowledge, uniqueness result was presented recently in [21] for the case of an impedance cylinder using the Lax-Phillips method.

In this work, we solve the inverse problem by formulating an equivalent system of non-linear integral equations that is solved using a regularized iterative scheme. This method was introduced by Kress and Rundell [18] and then considered in many different problems, in acoustic scattering problems [12, 13], in elasticity [5, 9] and in electrical impedance problem [8]. We propose an iterative scheme that is based on the idea of Johansson and Sleeman [14] applied to the inverse problem of recovering a perfect conducting cylinder. See [1, 5] for some recent applications. We assume integral representations for the solutions that results to a system consisting of four integral equations on the unknown boundary (considering the transmission conditions) and one on the unit circle (taking into account the asymptotic expansion of the solutions).

We solve this system in two steps. First, given an initial guess for the boundary curve we solve the well-posed subsystem (equations on the boundary) to obtain the corresponding densities and then we solve the linearized (with respect to the boundary) ill-posed far-field equation to update the initial approximation of the radial function. We consider Tikhonov regularization and the normal equations are solved by the conjugate gradient method.

The paper is organized as follows: in Section 2 we present the direct scattering problem, the elastic potentials and the equivalent system of integral equations that provide us with the far-field data. The inverse problem is stated in Section 3 where we construct an equivalent system of integral equation using the indirect integral equation method. In Section 4 the two-step method for the parametrized form of the system and the necessary Fréchet derivative of the integral operators are presented. The numerical examples give satisfactory results and justify the applicability of the proposed iterative scheme.

2 The direct problem

We consider the scattering of an electromagnetic wave by a penetrable cylinder in $\mathbbm{R}^{3}$. Let $\textbf{x}=(x,y,z)\in\mathbbm{R}^{3}.$ We denote by $\Omega_{int}=\{\textbf{x}:(x,y)\in\Omega,z\in\mathbbm{R}\}$ the cylinder, where $\Omega$ is a bounded domain in $\mathbbm{R}^{2}$ with smooth boundary $\Gamma.$ The cylinder $\Omega_{int}$ is oriented parallel to the $z$-axis and $\Omega$ is its horizontal cross section. We assume constant permittivity $\epsilon_{0}$ and permeability $\mu_{0}$ for the exterior domain $\Omega_{ext}:=\mathbbm{R}^{3}\setminus\overline{\Omega}_{int}.$ The interior domain $\Omega_{int}$ is also characterized by constant parameters $\epsilon_{1}$ and $\mu_{1}.$

We define the exterior magnetic $\textbf{H}^{ext}(\textbf{x},t)$ and electric field $\textbf{E}^{ext}(\textbf{x},t)$ for $\textbf{x}\in\Omega_{ext},\,t\in\mathbbm{R}$ and the interior fields $\textbf{H}^{int}(\textbf{x},t)$ and $\textbf{E}^{int}(\textbf{x},t)$ for $\textbf{x}\in\Omega_{int},\,t\in\mathbbm{R},$ that satisfy the Maxwell’s equations

[TABLE]

and the transmission conditions

[TABLE]

where $\bm{\hat{n}}$ is the outward normal vector, directed into $\Omega_{ext}$.

We illuminate the cylinder with an incident electromagnetic plane wave at oblique incidence, meaning transverse magnetic (TM) polarized wave. We define by $\theta$ the incident angle with respect to the negative $z$ axis and by $\phi$ the polar angle of the incident direction $\bm{\hat{d}}$ (in spherical coordinates), see Figure 1. Then, $\bm{\hat{d}}=(\sin\theta\cos\phi,\sin\theta\sin\phi,-\cos\theta)$ and the polarization vector is given by $\bm{\hat{p}}=(\cos\theta\cos\phi,\cos\theta\sin\phi,\sin\theta),$ satisfying $\bm{\hat{d}}\perp\bm{\hat{p}}$ and assuming that $\theta\in(0,\pi/2)\cup(\pi/2,\pi).$

In the following, due to the linearity of the problem we suppress the time-dependence of the fields and because of the cylindrical symmetry of the medium we express the incident fields as separable functions of $\bm{x}:=(x,y)$ and $z.$

Let $\omega>0$ be the frequency and $k_{0}=\omega\sqrt{\mu_{0}\epsilon_{0}}$ the wave number in $\Omega_{ext}$. We define $\beta=k_{0}\cos\theta$ and $\kappa_{0}=\sqrt{k_{0}^{2}-\beta^{2}}=k_{0}\sin\theta$ and it follows that the incident fields can be decomposed to [10]

[TABLE]

where

[TABLE]

After some calculations, we can reformulate Maxwell’s euqations (1) as a system of equations only for the $z$-component of the electric and magnetic fields [10]. The interior fields $e_{3}^{int}(\bm{x})$ and $h_{3}^{int}(\bm{x}),\ \bm{x}\in\Omega_{1}:=\Omega$ and the exterior fields $e_{3}^{ext}(\bm{x})$ and $h_{3}^{ext}(\bm{x}),$ $\bm{x}\in\Omega_{0}:=\mathbbm{R}^{2}\setminus\Omega$ satisfy the Helmholtz equations

[TABLE]

where $\kappa^{2}_{1}=\mu_{1}\epsilon_{1}\omega^{2}-\beta^{2}.$ Here, we assume $\mu_{1}\epsilon_{1}>\mu_{0}\epsilon_{0}\cos^{2}\theta$ in order to have $\kappa^{2}_{1}>0.$

The transmission conditions (2) can also be written only for the $z$-component of the fields. Let $(\bm{\hat{n}},\bm{\hat{\tau}})$ be a local coordinate system, where $\bm{\hat{n}}=(n_{1},n_{2})$ is the outward normal vector and $\bm{\hat{\tau}}=(-n_{2},n_{1})$ the outward tangent vector on $\Gamma.$ We define $\tfrac{\partial}{\partial n}=\bm{\hat{n}}\cdot\nabla_{t},$ $\tfrac{\partial}{\partial\tau}=\bm{\hat{\tau}}\cdot\nabla_{t},$ where $\nabla_{t}=\textbf{e}_{1}\tfrac{\partial}{\partial x}+\textbf{e}_{2}\frac{\partial}{\partial y}$ and $\textbf{e}_{1},\textbf{e}_{2}$ denote the unit vectors in $\mathbbm{R}^{2}.$ Then, we rewrite the boundary conditions as [10]

[TABLE]

where $\tilde{\mu}_{j}=\mu_{j}/\kappa_{j}^{2},\,\tilde{\epsilon}_{j}=\epsilon_{j}/\kappa_{j}^{2},\,\beta_{j}=\beta/\kappa_{j}^{2},$ for $j=0,1.$ The exterior fields are decomposed to $e_{3}^{ext}=e_{3}^{sc}+e_{3}^{inc}$ and $h_{3}^{ext}=h_{3}^{sc}+h_{3}^{inc},$ where $e_{3}^{sc}$ and $h_{3}^{sc}$ denote the scattered electric and magnetic field, respectively. From (3) we see that

[TABLE]

To ensure that the scattered fields are outgoing, we impose in addition the radiation conditions in $\mathbbm{R}^{2}:$

[TABLE]

where $r=|\bm{x}|,$ uniformly over all directions.

Now we are in position to formulate the direct transmission problem for oblique incident wave: Find the fields $h_{3}^{int},h_{3}^{sc},e_{3}^{int}$ and $e_{3}^{sc}$ that satisfy the Helmholtz equations (4), the transmission conditions (5) and the radiation conditions (7).

Theorem \thetheorem.

*If $\kappa_{1}^{2}$ is not an interior Dirichlet eigenvalue and $\kappa_{0}^{2}$ is not an interior Dirichlet and Neumann eigenvalue, then the direct transmission problem (4) – (7) admits a unique solution. *

Proof*:*

The proof is based on the integral representation of the solution resulting to a Fredholm type system of boundary integral equations. For more details see [10, Theorem 3.2]. $\square$

In the following, $j=0,1$ counts for the exterior ($\bm{x}\in\Omega_{0}$) and interior domain ($\bm{x}\in\Omega_{1}$), respectively. We introduce the single- and double-layer potentials defined by

[TABLE]

where $\Phi_{j}$ is the fundamental solution of the Helmholtz equation in $\mathbbm{R}^{2}:$

[TABLE]

and $H_{0}^{(1)}$ is the Hankel function of the first kind and zero order. We define also the integral operators

[TABLE]

The following theorem was proven in [10].

Theorem \thetheorem.

Let the assumptions of Theorem 2 still hold. Then, the potentials

[TABLE]

solve the direct transmission problem (4) – (7) provided that the densities $\phi_{0}\in H^{1/2}(\Gamma)$ and $\psi_{0}\in H^{1/2}(\Gamma)$ satisfy the system of integral equations

[TABLE]

where

[TABLE]

and $K_{j}:=\left(NS_{j}\pm\tfrac{1}{2}I\right)^{-1}ND_{j},\,L_{j}:=2(TD_{j}-TS_{j}K_{j}).$ The rest of the densities satisfy $\phi_{1}=\phi_{0}+e^{inc}_{3},\,\psi_{1}=\psi_{0},\,\eta_{j}=K_{j}\phi_{j}$ and $\xi_{j}=K_{j}\psi_{j}.$

The solutions $e_{3}^{sc}$ and $h_{3}^{sc}$ of (4) – (7) have the asymptotic behavior

[TABLE]

where $\bm{\hat{x}}=\bm{x}/|\bm{x}|.$ The pair $(e^{\infty},h^{\infty})$ is called the far-field pattern corresponding to the scattering problem (4) – (7). Its knowledge is essential for the inverse problem and using (11) we can compute it by

[TABLE]

where $\mathbbm{S}$ is the unit ball. The far-field operators are given by

[TABLE]

where $\Phi^{\infty}$ is the far-field of the Green function $\Phi,$ given by

[TABLE]

3 The inverse problem

The inverse scattering problem, we address here, reads: Find the shape and the position of the inclusion $\Omega,$ meaning reconstruct its boundary $\Gamma$, given the far-field patterns $(e^{\infty}(\bm{\hat{x}}),h^{\infty}(\bm{\hat{x}})),$ for all $\bm{\hat{x}}\in\mathbbm{S},$ for one or few incident fields (6).

3.1 The integral equation method

To solve the inverse problem we apply the method of nonlinear boundary integral equations, which in our case results to a system of four integral equations on the unknown boundary and one on the unit circle where the far-field data are defined. This method was first introduced in [18] and further considered in various inverse problems, see for instance [2, 3, 5, 9, 13]. Since the direct problem was solved with the direct method (Green’s formulas), in order to obtain our numerical data, here we adopt a different approach based on the indirect integral equation method, using simple representations for the fields.

We assume a double-layer representation for the interior fields and a single-layer representation for the exterior fields. Thus, we set

[TABLE]

Substituting the above representations in the transmission conditions (5) and considering the well-known jump relations, we get the system of integral equations

[TABLE]

In addition, given the far-field operators (15) and the representations (16) of the exterior fields we see that the unknown boundary $\Gamma$ and the densities $\psi_{e}$ and $\psi_{h}$ satisfy also the far-field equations

[TABLE]

where the right-hand sides are the known far-field patterns from the direct problem. The equation (17) in matrix form reads

[TABLE]

where

[TABLE]

The matrix T due to its special form and the boundness of $\partial_{\tau}:H^{1/2}(\Gamma)\rightarrow H^{-1/2}(\Gamma)$ has a bounded inverse given by

[TABLE]

Then, equation (19) takes the form

[TABLE]

where now I is the identity matrix and

[TABLE]

Using the mapping properties of the integral operators [17], we see that the operator $\textbf{C}:(H^{-1/2}(\Gamma)\times H^{1/2}(\Gamma))^{2}\rightarrow(H^{-3/2}(\Gamma)\times H^{-1/2}(\Gamma))^{2}$ is compact.

We observe that we have six equations (21) and (18) for the five unknowns: $\Gamma$ and the four densities. Thus, we consider the linear combination $\tilde{\epsilon}_{0}\cdot$(18a) + $\tilde{\mu}_{0}\cdot$(18b) as a replacement for the far-field equations in order to state the following theorem as a formulation of the inverse problem.

Theorem \thetheorem.

Given the incident field (6) and the far-field patterns $(e^{\infty}(\bm{\hat{x}}),h^{\infty}(\bm{\hat{x}})),$ for all $\bm{\hat{x}}\in\mathbbm{S},$ if the boundary $\Gamma$ and the densities $\psi_{e},\,\phi_{h},\,\psi_{h}$ and $\phi_{e}$ satisfy the system of equations

[TABLE]

*then, $\Gamma$ solves the inverse problem. *

The integral operators in (22) are linear with respect to the densities but non-linear with respect to the unknown boundary $\Gamma.$ The smoothness of the kernels in the far-field equation (22e) reflects the ill-posedness of the inverse problem.

To solve the above system of equations, we consider the method first introduced in [14] and then applied in different problems, see for instance [1, 5, 19]. More precisely, given an initial approximation for the boundary $\Gamma$, we solve the subsystem (22a) - (22d) for the densities $\psi_{e},\,\phi_{h},\,\psi_{h}$ and $\phi_{e}.$ Then, keeping the densities $\psi_{e}$ and $\psi_{h}$ fixed we linearize the far-field equation (22e) with respect to the boundary. The linearized equation is solved to obtain the update for the boundary. The linearization is performed using Fréchet derivatives of the operators and we also regularize the ill-posed last equation.

To present the proposed method in details, we consider the following parametrization for the boundary

[TABLE]

where $\bm{z}:\mathbbm{R}\rightarrow\mathbbm{R}^{2}$ is a $C^{2}$-smooth, $2\pi$-periodic, injective in $[0,2\pi),$ meaning that $\bm{z}^{\prime}(t)\neq 0,$ for all $t\in[0,2\pi].$ The non-negative function $r$ represents the radial distance of $\Gamma$ from the origin. Then, we define

[TABLE]

and the parametrized form of (22) is given by

[TABLE]

with the parametrized operators

[TABLE]

and the right-hand side

[TABLE]

In addition, we set $\bm{\mathcal{A}}_{2}=\bm{\mathcal{B}}_{1}=\bm{\mathcal{B}}_{4}=\bm{\mathcal{B}}_{5}=\bm{\mathcal{C}}_{4}=\bm{\mathcal{D}}_{2}=\bm{\mathcal{D}}_{3}=\bm{\mathcal{D}}_{5}=\bm{\mathcal{F}}_{2}=0.$ The matrix $M^{K_{j}}$ denotes the discretized kernel of the operator $K_{j}.$ The explicit forms of the kernels can be found for example in [10, Equation 4.3]. The operators $\bm{\mathcal{A}}_{k},\,\bm{\mathcal{B}}_{k},\,\bm{\mathcal{C}}_{k},\,\bm{\mathcal{D}}_{k},$ $k=1,2,3,4,5$ act on the densities and the first variable $r$ shows the dependence on the unknown parametrization of the boundary. Only $\bm{\mathcal{F}}_{5}$ is independent of the radial function.

Let the function $q$ stand for the radial function of the perturbed boundary

[TABLE]

Then the iterative method reads:

Iterative Scheme \theiterative.

Let $r^{(0)}$ be an initial approximation of the radial function. Then, in the $k$th iteration step:

1. i.

We assume that we know $r^{(k-1)}$ and we solve the subsystem

[TABLE]

to obtain the densities $\zeta_{e}^{(k)},\,\xi_{h}^{(k)},\,\zeta_{h}^{(k)}$ and $\xi_{e}^{(k)}.$ 2. ii.

Keeping the densities $\zeta_{e}$ and $\zeta_{h}$ fixed, we linearize the fifth equation of (23), namely

[TABLE]

We solve this equation for $q$ and we update the radial function $r^{(k)}=r^{(k-1)}+q.$

*The iteration stops when a suitable stopping criterion is satisfied. *

Remark \theremark.

In order to take advantage of the available measurement data, we can also keep the overdetermined system (17) and (18) instead of (22e) and replace equation (25) with

[TABLE]

*where now $\bm{\mathcal{F}}_{e}=\tilde{\epsilon}_{0}\,e^{\infty}$ and $\bm{\mathcal{F}}_{h}=\tilde{\mu}_{0}\,h^{\infty}.$ *

The Fréchet derivatives of the operators are calculated by formally differentiating their kernels with respect to $r$

[TABLE]

Recall that $\bm{\mathcal{A}}_{5}=\bm{\mathcal{C}}_{5}=S^{\infty}.$ If $\kappa_{0}^{2}$ is not an interior Neumann eigenvalue, then the operator $\bm{\mathcal{A}}^{\prime}_{5}$ is injective [13]. Using similar arguments as in [1, 11] we can relate the above iterative scheme to the classical Newton’s method.

The iterative scheme 3.1 can also be generalized to the case of multiple illuminations $\textbf{e}^{inc}_{l},\,l=1,...,L.$

Iterative Scheme \theiterative (Multiple illuminations).

Let $r^{(0)}$ be an initial approximation of the radial function. Then, in the $k$th iteration step:

1. i.

We assume that we know $r^{(k-1)}$ and we solve the $L$ subsystems

[TABLE]

to obtain the densities $\zeta_{e,l}^{(k)},\,\xi_{h,l}^{(k)},\,\zeta_{h,l}^{(k)}$ and $\xi_{e,l}^{(k)}.$ 2. ii.

Then, keeping the densities fixed, we solve the overdetermined version of the linearized fifth equation of (23)

[TABLE]

for $q$ and we update the radial function $r^{(k)}=r^{(k-1)}+q.$

*The iteration stops when a suitable stopping criterion is satisfied. *

4 Numerical implementation

In this section, we present numerical examples that illustrate the applicability of the proposed method. We use quadrature rules for integrating the singularities considering trigonometric interpolation. The convergence and error analyses are given in [15, 16]. Then, the system of integral equations is solved using the Nyström method. The parametrized forms of the integral operators are presented in [10, Section 4]. We approximate the smooth kernels with the trapezoidal rule and the singular ones with the well-known quadratures rules [16].

In the following examples, we consider two different boundary curves. A peanut-shaped and an apple-shaped boundary with radial function

[TABLE]

and

[TABLE]

respectively.

To avoid an inverse crime, we construct the simulated far-field data using the numerical scheme (12) and considering double amount of quadrature points compared to the inverse problem. We approximate the radial function $q$ by a trigonometric polynomial of the form

[TABLE]

and we consider $2n$ equidistant points $t_{j}=j\pi/n,\,j=0,...,2n-1.$ The well-posed subsystem (24) does not require any special treatment. The ill-posed linearized far-field equation (25) is solved by Tikhonov regularization. We rewrite (25) as

[TABLE]

for $\zeta^{(k)}:=\zeta_{e}^{(k)}+\zeta_{h}^{(k)},$ and we decompose (27) as

[TABLE]

where

[TABLE]

We replace the derivative of $q$ by the derivative of the trigonometric interpolation polynomial

[TABLE]

with weight

[TABLE]

Then, at the $k$th step we minimize the Tikhonov functional of the discretized equation

[TABLE]

where $\textbf{x}\in\mathbbm{R}^{(2m+1)\times 1}$ is the vector with the unknowns coefficients $a_{0},...,a_{m},$ $b_{1},...,b_{m}$ of the radial function, and $\textbf{A}\in\mathbbm{C}^{2n\times 2n},\,\textbf{b}\in\mathbbm{C}^{2n\times 1}$ are given by

[TABLE]

for $k,j=0,...,2n-1.$ The multiplication matrix $\textbf{T}\in\mathbbm{R}^{2n\times(2m+1)}$ stands for the trigonometric functions of the approximated radial function and is given by

[TABLE]

Here $p\geq 0$ defines the corresponding Sobolev norm. Since $q$ is real valued we solve the following regularized equation

[TABLE]

on the $k$th step, where the matrix $\textbf{I}_{p}\in\mathbbm{R}^{(2m+1)\times(2m+1)}$ corresponds to the Sobolev $H^{p}$ penalty term. We solve (31) using the conjugate gradient method. We update the regularization parameter in each iteration step $k$ by

[TABLE]

for some given initial parameter $\lambda_{0}>0.$ To test the stability of the iterative method against noisy data, we add also noise to the far-field patterns with respect to the $L^{2}-$norm

[TABLE]

for some given noise levels $\delta_{1},\,\delta_{2}$ where $u=u_{1}+\mathrm{i}u_{2},\,v=v_{1}+\mathrm{i}v_{2},$ for $u_{1},u_{2},v_{1},v_{2}\in\mathbbm{R}$ normally distributed random variables.

Already in simpler cases [1], the knowledge of the far-field patterns for one incident wave is not enough to produce satisfactory reconstructions. Thus, we will also use multiple incident directions. To do so, we have to consider different values of the polar angle $\phi$ since in $\mathbbm{R}^{2},$ as we see from (6), corresponds to the incident direction. We set

[TABLE]

4.1 Numerical results

In all examples we set the exterior parameters $(\epsilon_{0},\mu_{0})=(1,1)$ and the interior $(\epsilon_{1},\mu_{1})=(2,2).$ We use $n=64$ collocation points for the direct problem and $n=32$ for the inverse. We set $\theta=\pi/3$ and $\lambda_{0}\in[0.5,\,0.8]$ as the initial regularization parameter.

We present reconstructions for different boundary curves, different number of incident directions and initial guesses for exact and perturbed far-field data. In all figures the initial guess is a circle with radius $r_{0},$ a green solid line, the exact curve is represented by a dashed red line and the reconstructed by a solid blue line. The arrows denote the directions of the incoming incident fields.

In the first three examples we consider the peanut-shaped boundary. In the first example, the regularized equation (31) is solved with $L^{2}$ penalty term, meaning $p=0$ and $m=3$ coefficients. We solve equation (26) for different incident directions. The reconstructions for $\omega=2.5$ and $r_{0}=0.6$ are presented in Figure 2 for two incident fields with directions $\textbf{d}_{l+1/2}.$ On the left picture, we see the reconstructed curve for exact data and 9 iterations and on the right picture for noisy data with $\delta_{1}=\delta_{2}=5\%$ and 14 iterations. In the second example, we consider equation (25), four incident fields, noisy data $\delta_{1}=\delta_{2}=5\%$ and we keep all the parameters as before. The reconstructions for $r_{0}=0.6$ and 14 iterations are shown in the left picture of Figure 4, and for $r_{0}=1$ and 20 iterations in the right one. We set $m=5$ and $p=1$ ($H^{1}$ penalty term) in the third example. The results for $r_{0}=1$ and four incident fields are shown in Figure 4. Here $\omega=2$ and we use equation (26). We need 26 iterations for the exact data and 30 iterations for the noisy data ($\delta_{1}=\delta_{2}=3\%$).

In the last two examples we consider the apple-shaped boundary, $H^{1}$ penalty term, $\omega=3$ and $m=3$ coefficients. In the fourth example, we consider equation (25), noise-free data and four incident fields in order to examine the dependence of the iterative scheme on the initial radial guess. On the left picture of Figure 6, we see the reconstructed curve for $r_{0}=0.5$ after 13 iterations and on the right picture for $r_{0}=1$ after 20 iterations. In the last example we consider $\delta_{1}=\delta_{2}=3\%$ noise and $r_{0}=0.6.$ Figure 6 shows the improvement of the reconstruction for more incident fields. On the left picture we see the results for three incident fields, equation (26) and 7 iterations and the reconstructed curve for 4 incident fields, equation (25) and 15 iterations is shown on the right picture.

Our examples show the feasibility of the proposed iterative scheme and the stability against noisy data. Considering more than one incident field improves considerably the reconstructions. The choice of the initial guess is also crucial.