TY - JOUR
AU - Milan,, Giacomo
AB - Abstract We investigate the reciprocity gap functional method, which has been developed in the inverse scattering theory, in the context of electrical impedance tomography. In particular, we aim to reconstruct an inclusion contained in a body, whose conductivity is different from the conductivity of the surrounding material. Numerical examples are given, showing the performance of our algorithm. 1. Introduction In this paper, we consider the inverse problem of reconstructing an inclusion contained in an electrical conductor by performing electrostatic boundary measurements. Let |$\varOmega $| be an electrical conductor inside which an anomalous region is located, i.e. a region whose conductivity is different from the conductivity of the surrounding material. For instance, assuming that the conductivity of the background medium is equal to |$1$|⁠, the conductivity inside the inclusion |$D$| will be equal to an unknown constant |$k>0$|⁠, |$k\neq 1$|⁠. Therefore, denoting by |$\gamma $| the conductivity function, |$\gamma $| will be of the form |$\gamma (x)=1+(k-1)\chi _D$|⁠, where |$\chi _D$| is the characteristic function of the set |$D$|⁠. Prescribing a voltage |$f\in H^{1/2}(\partial \varOmega )$| on the boundary of |$\varOmega $|⁠, the induced potential |$u$| will be the solution of the problem $$\begin{align*} &\left\{\begin{array}{rl} \textrm{div}\ \left(\gamma(x)\nabla u\right)=0 & \textrm{in}\ \varOmega,\\ u=f & \textrm{on}\ \partial\varOmega. \end{array}\right. \end{align*}$$ The current density that we measure on the boundary is given by the normal derivative of the potential |$\partial _\nu u|_{\partial \varOmega }$|⁠. The inverse problem we are addressing to is to recover information on the inclusion |$D$| performing infinitely many boundary measurements on |$\partial \varOmega $|⁠. This problem is a special instance of the classical inverse conductivity problem (Calderón, 1980) and the uniqueness issue has been solved by Isakov (1988). The stability analysis has been performed in Alessandrini & Di Cristo (2005) for piecewise constant conductivity and extended later in Di Cristo (2007) for the variable coefficients case. The modulus of continuity obtained is of logarithmic type and it turns out to be optimal as shown by examples in Di Cristo & Rondi (2003). This emphasizes the difficulty in any numerical reconstruction as small errors in the data can lead to large errors in the solution. In this paper, we focus our attention on the reconstruction issue. In particular, we provide a method to detect the inclusion and implement it showing some examples. The method we propose can be traced back to Colton & Haddar (2005) where a new version of the linear sampling method (Colton et al., 2003) based on the reciprocity gap functional (Andrieux & Ben Abda, 1996) has been introduced. We also refer to Cakoni et al. (2006) for a further analysis of this method and a proper comparison with the previous ones. For the boundary measurements case in the inverse scattering framework, this method has been developed in Di Cristo & Sun (2006, 2007) and Monk & Sun (2007) for penetrable and impenetrable obstacles and for cavities in Sun et al. (2016) and Zeng et al. (2017) (see also Cakoni et al., 2012, for multilayered media). We apply this procedure to the inverse inclusion setting. The main difference with respect to the other cases is that the unknown parameter is in the highest-order term of the partial differential equation that models the problem. In particular in the integral equation that we have to solve (see (2.10) below), specific singular functions (see Definition 2.3) have to be used. These have been first introduced in a similar problem by Somersalo (2001) and they take into account the fact that the unknown is in the highest-order term of the partial differential equation (PDE). The paper is organized as follows. In Section 2, we formulate the problem, setting the notations, and we state our main theorem whose proof is shown in Section 3 based on some preliminary properties. In Section 4, we provide several numerical examples based on the algorithm developed in this article. In particular, we compare the algorithm with the linear sampling method, but we would like to point out to the reader other possible approaches like the factorization method directly descending from the linear sampling method (see Kirsch, 2011) or iterative schemes like that contained in Delbary & Kress, 2010. Conclusions are presented in Section 5. 2. Notations and main result To begin with, let us introduce some notations and definitions, which we will use through out the paper. We denote by |$\varOmega $| a bounded domain in |$\mathbb R^n$|⁠, |$n\geq 2$| in which an inclusion |$D$| is located, such that |$D$| is of Lipschitz class and |$\varOmega \setminus \overline D$| is connected. We also denote by |$B$| a ball containing the domain |$\varOmega $|⁠. Denoting by |$\gamma (x)$| the conductivity function, it will be of the form $$\begin{align*} &\gamma(x)=1+(k-1)\chi_D,\end{align*}$$ where |$\chi _D$| is the characteristic function of the set |$D$| and |$k>0$|⁠, |$k\neq 1$|⁠, is a positive unknown constant. Denoting by |$u$| and |$v$| the two solutions of the conductivity problems with the inclusion and without the inclusion respectively, these functions satisfy the equations $$\begin{align} \textrm{div}((1 + (k - 1) \chi_D)\nabla u)&=0\ \textrm{in}\ \varOmega, \end{align}$$(2.1) $$\begin{align}\varDelta v = \textrm{div}(\nabla v)&= 0\ \textrm{in}\ \varOmega \end{align}$$(2.2) and we can use them to define the reciprocity gap functional |$\mathcal{R}$|⁠. Definition 2.1 Let |$u,v$| be solutions to (2.1) and (2.2), respectively, we define the reciprocity gap functional |$\mathcal{R}$| as $$\begin{align}& \mathcal{R}\bigl(u, v\bigr)\coloneqq \int_{\partial \varOmega} \bigl(u(\partial_\nu v) - (\partial_\nu u) v\bigr), \end{align}$$(2.3) where |$\partial _\nu $| denotes the partial derivative with respect to the outer unit normal vector |$\nu $| to |$\partial \varOmega $|⁠. In the sequel, let |$\varOmega \subset B$|⁠. The idea is to test the reciprocity gap functional with singular solutions and move the singularity to detect the inclusion. For this purpose, let us introduce |$\varPhi _D(x,x_0)$| the fundamental solution of the equation (2.1), i.e. $$\begin{align*} &\textrm{div}((1+(k-1)\chi_D)\nabla\varPhi_D(x,x_0))=-\delta(x-x_0),\end{align*}$$ with |$x_0\in \partial B$|⁠. We will denote by $$\begin{align} \mathcal{U}&\coloneqq\bigl\{\varPhi_D(x, x_0): x_0\in \partial B\bigr\}, \end{align}$$(2.4) $$\begin{align} \mathcal{V}(\overline{B})&\coloneqq\bigl\{v\in C^2(B) \cap C(\overline{B}):\varDelta v = 0\bigr\}, \end{align}$$(2.5) $$\begin{align} \mathcal{V}_0(\overline{B})&\coloneqq\left\{v\in\mathcal{V}(\overline{B}):\int_{\partial B} v = 0\right\}=\mathcal{V}(\overline{B})/\textup{span}\bigl(\{1\}\bigr), \end{align}$$(2.6) the spaces of functions we need to set the reciprocity gap equation. The classical definition for |$\mathcal{V}(\overline{B})$| is automatically implied by the use of layer potentials, easy to discretize, but Theorem 2.4 still remains valid for solutions |$\varDelta v=0$| in the weak sense with |$v\in H^1(B)$|⁠. The reciprocity gap functional restricted to |$\mathcal V({\overline B})$| can be seen as an operator, whose properties will be analysed more deeply in Section 3 (see Propositions 3.4 and 3.6 below). More precisely, we have the following definition. Definition 2.2 We denote by |$R: \mathcal{V}(\overline{B})\to L^2(\partial B)$| the ‘reciprocity gap operator’ defined by $$\begin{align}& R(v)(x_0)\coloneqq \mathcal{R}\bigl(u_{x_0},v\bigr) \coloneqq \int_{\partial \varOmega}\bigl(u_{x_0}\partial_\nu v - \partial_\nu u_{x_0}v\bigr), \end{align}$$(2.7) where |$u_{x_0}(x) = \varPhi _D(x, x_0) \in \mathcal{U}$|⁠. We introduce now a set of functions used in Somersalo (2001) that are singular in a ‘sampling point’ |$z\in \varOmega $| and represent the solution for the unperturbed case. Denoting by |$\varPhi $|⁠, the fundamental solution of the Laplace equation (2.2) in |$\mathbb{R}^n$|⁠, we have the following definition. Definition 2.3 We define $$\begin{align} &\vec{\varPsi}_z(x) = \vec{\varPsi}(x,z) \coloneqq \nabla\varPhi(x,z), \end{align}$$(2.8) $$\begin{align} &\varPsi_z(x,\vec{d})\coloneqq \vec{\varPsi}_z\cdot \vec{d} = \partial_{\vec{d}\,\,}\varPhi(x,z), \end{align}$$(2.9) where |$\vec{d}$| is a unit vector. Sometimes, for brevity, we will write |$\varPsi _z$| for |$\varPsi _z(x,\vec{d})$|⁠, dropping the dependence on the unit vector. For a fixed |$z\in \varOmega $|⁠, the core part of our reconstruction algorithm is to look for a solution |$v \in \mathcal{V}(\overline{B})$| of the equation $$\begin{align}& \mathcal{R}(u, v) = \mathcal{R}(u, \varPsi_z)\quad \forall u \in \mathcal{U}, \ \end{align}$$(2.10) where |$\mathcal R$| is defined by (2.3). Theorem 2.4 (Reciprocity gap approximation theorem). Let |$D\subset \varOmega \subset B$| be an inclusion satisfying the aforementioned assumptions, let |$\{(u|_{\partial \varOmega },\partial _\nu u|_{\partial \varOmega })\}$| be a set of measured data with |$u\in \mathcal{U}$| as defined in (2.4) and let consider the set |$\mathcal{V}(\overline{B})$| of harmonic test functions defined in (2.5). Then if |$z \in D$| then there exists a sequence |$\{v_n\} \subset \mathcal{V}(\overline{B})$| such that $$\begin{align}& \lim_{n\to\infty}\mathcal{R}(u,v_n) = \mathcal{R}(u,\varPsi_z)\quad\forall u\in\mathcal{U}, \end{align}$$(2.11) where |$\varPsi _z$| as in Definition 2.3 and |$v_n|_{\partial D}$| converges to some |$g\in L^2(\partial D)$| and consequently it is bounded in the same norm |$\|v_n\|_{L^2(\partial D)}$|⁠; if |$z \in \varOmega \backslash D$| then any sequence |$\{v_n\} \subset \mathcal{V}(\overline{B})$| such that $$\begin{align}& \lim_{n\to\infty}\mathcal{R}(u,v_n) = \mathcal{R}(u,\varPsi_z)\quad\forall u\in\mathcal{U} \end{align}$$(2.12) is unbounded |$\|v_n\|_{L^2(\partial D)}\to \infty $| in |$L^2(\partial D)$|⁠. Remark 2.5 The numerical implementation consists in taking an appropriate set, like the single layer potentials |$\{\mathcal{S}(\partial B, \sigma ): \sigma \in C(\partial B)\}$|⁠, contained in |$\mathcal{V}(\overline{B})$|⁠, (this choice would not affect the final result), and then looking for a solution |$\sigma $| of $$\begin{align}& R_D \sigma = r_z, \end{align}$$(2.13) with |$R_D:C(\partial B)\to C(\partial B)$| and |$r_z\in C(\partial B)$| defined as $$\begin{align} &R_D\sigma(x_0) \coloneqq R(v)(x_0)=\mathcal{R}\bigl(u_{x_0},v\bigr), \end{align}$$(2.14) $$\begin{align} & r_z(x_0)\coloneqq \mathcal{R}\bigl(u_{x_0},\varPsi_z\bigr). \end{align}$$(2.15) 3. Proof of Theorem 2.4 Before proving our main result, we premise some lemmas and propositions that analyse the structure and the properties of the reciprocity gap operator defined in Definition 2.2 and the functional sets where it is defined. For the sake of simplicity, we restrict ourselves to the two-dimensional case. The other dimensions can be treated similarly through minor modifications. We will denote by |$\mathcal S(\partial B,\sigma )$| the single layer potential with density |$\sigma $| defined on |$\partial B$|⁠, i.e. $$\begin{align}& \mathcal S(\partial B,\sigma)(x)\coloneqq\int_{\partial B}\varPhi(x,y)\sigma(y)\,\textrm{d}y,\qquad x\in\mathbb R^2\setminus\partial B, \end{align}$$(3.1) and by |$\mathcal{D}(\partial B, \mu )$| the double layer potential defined as $$\begin{align}& \mathcal{D}(\partial B,\mu)(x)\coloneqq \int_{\partial B} \partial_{\nu(y)}\varPhi(x, y)\mu(y)\, \textrm{d}y,\quad x\in\mathbb{R}^2 \backslash\partial B. \end{align}$$(3.2) We have the following two properties of the space |$\mathcal V(\overline B)$|⁠. Let us mention here that similar results, treated in a more general setting, can be found in Bogomolny (1985). Lemma 3.1 (Density). The set |$\{v|_{\partial D}: v\in \mathcal{V}(\overline{B})\}$| is dense in |$L^2(\partial D)$|⁠. Proof. Let |$\beta \in L^2(\partial D)$| such that the product with all the harmonic functions |$v$| of the form |$v=\mathcal{S}(\partial B,\sigma ) + \mathcal{D}(\partial B, \mu )$| vanishes. Our aim is to show that |$\beta =0$|⁠. Computing the product of |$\beta $| with a single layer potential we get $$\begin{align*} 0=\langle\mathcal{S}({\partial B},\sigma), \beta\rangle_{L^2(\partial D)} &= \int_{\partial D} \int_{\partial B} \sigma(y) \varPhi(x,y)\,\textrm{d}y\,\beta(x)\,\textrm{d}x \\ &= \langle\sigma, \mathcal{S}({\partial D},\beta)\rangle_{L^2(\partial B)}\quad\forall\sigma\in C(\partial B), \end{align*}$$ which implies |$\mathcal{S}(\partial D,\beta ) = 0$| on |$\partial B$|⁠. Testing |$\beta $| with a double layer potential, we also get $$\begin{align*} 0=\langle\mathcal{D}({\partial B},\mu), \beta\rangle_{L^2(\partial D)} &= \int_{\partial D} \int_{\partial B} \mu(y) \partial_{\nu(y)}\varPhi(x,y)\,\textrm{d}y\,\beta(x)\,\textrm{d}x \\ &= \int_{\partial B} \partial_\nu \mathcal{S}(\partial D, \beta)(y) \mu(y)\, \textrm{d}y\quad \forall \mu\in C(\partial B), \end{align*}$$ which implies |$\partial _\nu \mathcal{S}(\partial D, \beta ) = 0$| on |$\partial B$|⁠. Therefore, the total flux, i.e. $$\begin{align*}& 0=\int_{\partial B} \partial_\nu \mathcal{S}(\partial D,\beta)(y)\,\textrm{d}y =\int_{\partial D} \int_{\partial B} \partial_{\nu(y)} \varPhi(x,y)\, \textrm{d}y \, \beta(x) dx = - \int_{\partial D} \beta(x) \, \textrm{d}x, \end{align*}$$ is zero. Noticing that |$\mathcal{S}(\partial D, \beta )$| is harmonic in |$\mathbb{R}^2\backslash \overline{B}$|⁠, and it has homogeneous Cauchy boundary conditions on |$\partial B$|⁠, by uniqueness result for exterior problem, we obtain |$\mathcal{S}(\partial D,\beta ) = 0$| in |$\mathbb{R}^2\backslash \overline{B}$|⁠. Then, by continuation principle, |$\mathcal{S}(\partial D,\beta ) = 0$| in |$\mathbb{R}^2\backslash \overline{D}$|⁠, by maximum principle in |$D$| we have |$\mathcal{S}(\partial D,\beta ) = 0$| in |$\mathbb{R}^2$| and finally by jump relation we conclude |$\beta = 0$|⁠. Lemma 3.2 (Density). The set |$\{(\partial _\nu v )|_{\partial D}:v\in \mathcal{V}(\overline{B})\}$| is dense in |$L^2_0(\partial D)\coloneqq L^2(\partial D)/\textup{span}(\{1\})$|⁠, i.e. the closure of smooth functions with zero mean. Proof. First, we observe that any constant function defined on |$\partial D$| is orthogonal to |$\partial _\nu v$| with |$v\in \mathcal{V}(\overline{B})$|⁠, indeed $$\begin{align*}& \langle \partial_\nu v, 1\rangle_{L^2(\partial D)} = \int_{\partial D}\partial_\nu v\,\textrm{d}y=0. \end{align*}$$ Now we fix a generic |$\beta \in L^2(\partial D)$| such that the product of |$\beta $| with all test functions in |$\mathcal{V}(\overline{B})$| is zero. Consider the set $$\begin{align}& \left\{\mathcal{S}(\partial B, \sigma)(x) = \int_{\partial B} \sigma(y)\varPhi(x,y)\,\textrm{d}y\,:\, \sigma \in C(\partial B)\right\} \subset \mathcal{V}(\overline{B}). \end{align}$$(3.3) For |$\sigma \in C(\partial B)$|⁠, by definition of adjoint operator, we have |$\langle \mathcal{A}\sigma ,\beta \rangle _{L^2(\partial D)} =\langle \sigma ,\mathcal{A}^*\beta \rangle _{L^2(\partial B)}=0$|⁠, i.e. $$\begin{align*} 0 = \int_{\partial D}\partial_\nu \mathcal{S}(\partial B, \sigma)(x) \beta(x) \, \textrm{d}x &= \int_{\partial D}\left( \int_{\partial B}\nabla_x \varPhi(x,y) \sigma(y)\,\textrm{d}y\right)\cdot\nu(x)\beta(x) \, \textrm{d}x \\ &= \int_{\partial B}\left( \int_{\partial D}\partial_{\nu(x)} \varPhi(x,y) \beta(x)\,\textrm{d}x\right)\sigma(y) \, \textrm{d}y. \end{align*}$$ Defining the double layer potential $$\begin{align}& f(y)\coloneqq \mathcal{D}(\partial D, \beta)(y)=\int_{\partial D}\partial_{\nu(x)} \varPhi(x,y) \beta(x)\,\textrm{d}x, \end{align}$$(3.4) we have |$f(y)=0$| on |$\partial B$|⁠. Also, |$f$| is harmonic in |$\mathbb{R}^2\backslash \partial D$| and bounded. By uniqueness of bounded harmonic functions for exterior Dirichlet problem and by continuation principle, we obtain |$f=0$| in |$\mathbb{R}^2\backslash \overline{D}$|⁠. Therefore, the double layer potential |$f$| vanishes outside |$D$|⁠, and the boundary integral operator has nullspace equal to |$\textup{span}(\{1\})$|⁠. This forces |$\beta $| to be a generic constant and this shows that |$\{(\partial _\nu v )|_{\partial D}\}^\perp = \textup{span}(\{1\})$|⁠, which concludes the proof. Remark 3.3 If we consider |$v\in \mathcal{V}_0(\overline{B})$| instead of |$\mathcal{V}(\overline{B})$|⁠, since the set |$\{(\partial _\nu v )|_{\partial D}\}$| is the same in both cases, then the previous density result is still valid. Now let us better characterize the reciprocity gap operator |$R$|⁠. In particular, we have the following injectivity and surjectivity properties. Proposition 3.4 (Injectivity). The operator |$R$| is injective in |$\mathcal{V}_0(\overline{B})$|⁠. The kernel of |$R$| is the one dimensional vector space |$\mathcal{N}(R)=\textup{span}(\{1\})$|⁠, spanned by the constant function. Remark 3.5 Observe, by simple substitution, that for any constant function, |$R(v_c)=0$|⁠, where |$v_c(x)\coloneqq c$|⁠, |$c\in \mathbb{C}$|⁠. Proof. By Remark 3.5, we just need to show that the only harmonic functions such that |$R(v)(x_0)=0$| for |$x_0\in \partial B$| are constant functions. Let |$v\in \mathcal{V}(\overline{B})$| such that |$R(v)=0$|⁠. Applying the Green’s representation formula in the open set |$\varOmega \backslash \overline{D}$| to |$u\in \mathcal{U}$| and |$v$|⁠, since they are both harmonic, with continuous extension up to the boundary, we get $$\begin{align}& \textrm{R}(v)(x_0)\coloneqq\mathcal{R}_{\partial\varOmega}\bigl(u(\cdot,x_0),v(\cdot)\bigr) = \mathcal{R}_{\partial D}\bigl(u^+(\cdot,x_0),v(\cdot)\bigr) = 0, \end{align}$$(3.5) where |$\mathcal R_{\partial \varOmega }$| and |$\mathcal{R}_{\partial D}$| are defined by (2.3) with the boundary integral taken on |$\partial \varOmega $| and |$\partial D$|⁠, respectively, and |$u^+=u\chi _D$|⁠. By transmission conditions $$\begin{align} \bigl[u\bigr]^+_-\bigl(y\bigr)&=u^+(y) - u^-(y) =0, & y \in\partial D, \end{align}$$(3.6) $$\begin{align} \bigl[\partial_\gamma u\bigr]^+_-\bigl(y\bigr) &= \partial_\nu u^+(y) - k\,\partial_\nu u^-(y) =0, & y \in\partial D, \end{align}$$(3.7) where |$\mathcal{R}_{\partial \varOmega }$| and |$\mathcal{R}_{\partial D}$| are defined (2.3) with the boundary integral taken on |$\partial \varOmega $| and |$\partial D$|⁠, respectively, and |$u^+=u\chi _D$|⁠. By Green’s identity, we have $$\begin{align*} \mathcal{R}_{\partial D}\bigl(u^+(\cdot,x_0),v(\cdot)\bigr) &= \int_{\partial D}\bigl(u^-(y,x_0)\partial_\nu v (y) - k \cdot \partial_\nu u^-(y,x_0)v(y)\bigr)\,\textrm{d}y \\ &= \int_{\partial D}u^-(y,x_0)\partial_\nu v (y)\,\textrm{d}y - k\int_{\partial D}u^-(y,x_0)\partial_\nu v (y)\,\textrm{d}y \\ &= (1-k)\int_{\partial D}u^-(y,x_0)\partial_\nu v (y)\,\textrm{d}y = 0. \end{align*}$$ Now define $$\begin{align}& f(x_0)\coloneqq \int_{\partial D} u(y,x_0)\partial_\nu v(y)\,\textrm{d}y\quad x_0 \in \mathbb{R}^2\backslash\partial D. \end{align}$$(3.8) We note that |$f(x_0)=0$| on |$\partial B$| and that |$f(x_0)$| is harmonic in |$\mathbb{R}^2\backslash \overline{D}$|⁠, by reciprocity relation. Also the density |$\partial _\nu v(y)$| has zero mean, since |$v$| is harmonic. This implies that in |$\mathbb{R}^2$|⁠, |$f(x_0)$| is bounded at infinity, while in |$\mathbb{R}^3$| boundedness of |$D$| prevents any kind of this bad growth. By uniqueness of exterior boundary value problems in |$\mathbb{R}^2\backslash \overline{B}$|⁠, and by continuation principle of analytic functions from the open set |$\mathbb{R}^2\backslash \overline{B}$| to the open set |$\mathbb{R}^2\backslash \overline{D}$|⁠, we conclude |$f=0$| in |$\mathbb{R}^2\backslash \overline{D}$|⁠. By continuity on |$\partial D$|⁠, and by uniqueness in |$D$|⁠, we conclude that |$f(x_0)=0$| in |$\mathbb{R}^2$|⁠. By jump the relations of the derivatives of |$f(x_0)$|⁠, we have |$\partial _\nu v = 0$|⁠. Finally, imposing homogeneous Neumann condition to the interior problem in |$D$|⁠, we conclude that |$v$| is constant in |$D$| and in |$\overline{B}$|⁠. Proposition 3.6 (Surjectivity). The range of the operator |$R: \mathcal{V}(\overline{B})\to L^2(\partial B)$| has codimension one, indeed |$\mathcal{R}(R)^\perp =\textup{span}(\{\alpha _c\})$|⁠, where |$\alpha _c$| is the solution of $$\begin{align}& \left(K^{\prime}+\dfrac{1}{2}I\right)\alpha_{c}= 0\quad\textup{ on} \partial B, \end{align}$$(3.9) where |$K^{\prime}$| is the adjoint operator of the double layer operator |$K:L^2(\partial B)\to L^2(\partial B)$|⁠. Proof. Let us consider |$\alpha \in L^2(\partial B)$| such that $$\begin{align*}& 0 =\langle\,R(v), \alpha\rangle_{L^2(\partial B)} = \langle \mathcal{R}\bigl(u(\cdot,x_0),v(\cdot)\bigr), \alpha(x_0)\rangle_{L^2(\partial B)}, \end{align*}$$ for all |$v\in \mathcal{V}(\overline{B})$|⁠. Switching the order of integration between the source boundary |$\partial B$| and the observation boundary |$\partial \varOmega $| and defining a new layer potential with density |$\alpha $|⁠, we have $$\begin{align}& 0 = \mathcal{R}(\mathcal{E},v),\quad\textup{ where} \quad\mathcal{E}(x):=\int_{\partial B}u(x,x_0)\alpha(x_0)\, \textrm{d}x_0. \end{align}$$(3.10) By Green’s identity and the transmission conditions (3.6)–(3.7) for |$\mathcal{E}$|⁠, as in the proof of Proposition 3.4, we have $$\begin{align}& 0 = \mathcal{R}(\mathcal{E},v)=(1-k)\int_{\partial D} \partial_\nu\mathcal{E}^-(y)v(y)\,\textrm{d}y,\quad \forall v \in \mathcal{V}(\overline{B}). \end{align}$$(3.11) By density Lemma 3.2, we obtain |$\partial _\nu \mathcal{E}^- = k\,\partial _\nu \mathcal{E}^+= 0$| on |$\partial D$|⁠. Now since |$u(x,x_0)=\varPhi (x,x_0)+\mathcal{S}(\partial D, \psi _{x_0})$|⁠, we get $$\begin{align} \mathcal{E}(x)&=\mathcal{S}(\partial B, \alpha) + \int_{\partial D} \left(\int_{\partial B}\psi(y,x_0)\alpha(x_0)\,\textrm{d}x_0\right)\varPhi(x,y)\,\textrm{d}y \nonumber\\ &=\mathcal{S}(\partial B, \alpha) + \mathcal{S}(\partial D, \psi_\alpha). \end{align}$$(3.12) By uniqueness of the interior homogeneous Neumann problem, |$\mathcal{E} = c$| in |$D$|⁠, where the constant |$c$| can be computed from the definition of |$\mathcal{E}$|⁠. By the previous jump relation on |$\partial D$|⁠, we get that |$\psi _\alpha =0$| on |$\partial D$|⁠. Now reformulating the expression for |$\mathcal{E}$|⁠, we get $$\begin{align}& \mathcal{E}(x)= \mathcal{S}(\partial B, \alpha) = \int_{\partial B}\varPhi(x,x_0)\alpha(x_0)\, \textrm{d}x_0, \end{align}$$(3.13) and by continuation principle from |$D$| to |$B$|⁠, |$\mathcal{E}(x)=c$| in |$B$|⁠. So, we are forced to look for a single layer potential, if there exists, constant in |$\overline{B}$|⁠. Actually, |$\alpha $| solves the boundary integral equation, which we obtain by imposing the equivalent condition |$\partial _\nu \mathcal{E}^-|_{\partial B}=0$|⁠, i.e. $$\begin{align}& \partial_\nu\mathcal{E}^-=\left(K^{\prime}+\dfrac{1}{2}I\right)\alpha= 0,\quad\textrm{on}\ \partial B. \end{align}$$(3.14) We know that the kernel |$\mathcal{N}(K^{\prime} + 1/2 I)=\textup{span}(\{\alpha _{c}\})$| is one dimensional, where we have denoted by |$\alpha _{c}$| the physical charge layer, which can be measured on the surface of a conductor, with the property of a constant potential inside it. On the other hand, if |$\mathcal{S}(\partial B, \alpha _{c})$| is constant in |$\overline{B}$|⁠, $$\begin{align}& \mathcal{E}_{c}(x)\coloneqq\int_{\partial B}u(x,x_0)\alpha_{c}(x_0)\,\textrm{d}x_0 \end{align}$$(3.15) is also constant in |$\overline{B}$|⁠, since (3.12) solves the equation on |$\partial D$|⁠, by the transmission conditions, $$\begin{align}& \left(K^{\prime} + \frac{\hat{c}(k)}{2} I\right)\psi_{\alpha_c} = -\partial_\nu \mathcal{S}(\partial B, \alpha_c) = 0, \end{align}$$(3.16) where |$\hat{c}(k)\coloneqq ({k+1})/({k-1})$|⁠, yielding to |$\psi _{\alpha _c}=0$|⁠, since we prescribed |$k>0, k\neq 1$| and |$\textrm{spectrum}(K^{\prime})\subset [-1/2, 1/2)$| (see Kress, 1999). Testing the initial product with any |$v\in \mathcal{V}(\overline{B})$|⁠, we conclude that |$\alpha _{c}$| is orthogonal to all of them, i.e. $$\begin{align}& \langle \, R(v),\alpha_{c}\rangle_{L^2(\partial B)}=\mathcal{R}_{\partial \varOmega}(\mathcal{E}_{c}, v)=\int_{\partial\varOmega}\bigl(\mathcal{E}_{c}\,\partial_\nu v - \partial_\nu\mathcal{E}_{c}\, v\bigr)\,\textrm{d}y=0, \end{align}$$(3.17) because |$\mathcal{E}_{c}=c$|⁠, |$\partial _\nu \mathcal{E}_{c}=0$| on |$\partial \varOmega $|⁠, and |$\int _{\partial \varOmega }\partial _\nu v=0$|⁠. We can now proceed with the proof of our approximation theorem. Proof. of Theorem 2.4. We split the proof in two parts. For a fixed |$z \in D$|⁠, and using Neumann-to-Dirichlet operator of Laplace equation in |$D$|⁠, |${N_0}:H^{-1/2}_0(\partial D)\to H^{1/2}(\partial D)$|⁠, we have $$\begin{align} &\mathcal{R}_{\partial \varOmega}(u,\varPsi_z) = \mathcal{R}_{\partial D}(u,\varPsi_z) = \int_{\partial D}(u\partial_\nu\varPsi_z - \partial_\nu u^+ \varPsi_z) \,\textrm{d}y \nonumber\\ =& \int_{\partial D}\bigl(\partial_\nu u^- {N_0}(\partial_\nu\varPsi_z) - k \partial_\nu u^- \varPsi_z\bigr) \,\textrm{d}y = \int_{\partial D}\partial_\nu u^-\Big( {N_0}(\partial_\nu\varPsi_z) - k \varPsi_z)\Big) \,\textrm{d}y. \end{align}$$(3.18) On the other hand, for |$v\in \mathcal{V}(\overline{B})$|⁠, we have $$\begin{align}& \mathcal{R}_{\partial \varOmega}(u,v) = \mathcal{R}_{\partial D}(u,v) = \int_{\partial D}(1-k)(\partial_\nu u^- )v\,\textrm{d}y. \end{align}$$(3.19) By Lemma 3.1, there exists a sequence |$\{v_n\}\subset \mathcal{V}(\overline{B})$| converging strongly and thus weakly in |$L^2(\partial D)$| to the argument of the right integral |$\bigl ( {N_0}(\partial _\nu \varPsi _z) - k \varPsi _z)\bigr )(1-k)^{-1}\in H^{1/2}(\partial D)$|⁠. The convergence holds for all |$u\in \mathcal{U}$|⁠, but in general the convergence is not uniform in |$x_0\in \partial B$|⁠. For a fixed |$z \in \varOmega \backslash D$|⁠, using Green’s representation formula for derivatives of a harmonic function |$u$| in |$E=\varOmega \backslash \overline{D}$|⁠, $$\begin{align}& \partial_{\vec{d}\,}u(z) = \int_{\partial E}(u\, \partial_{\nu} \varPsi_z - \partial_\nu u \, \varPsi_z)\,\textrm{d}y. \end{align}$$(3.20) Then, on the right-hand side of (2.13) by reciprocity relation, we have $$\begin{align}& \mathcal{R}_{\partial \varOmega}(u,\varPsi_z) = \partial_{\vec{d}}\,u(z,x_0) + \mathcal{R}_{\partial D}(u,\varPsi_z) = \partial_{\vec{d}}\,u(z,x_0) + w(x_0), \end{align}$$(3.21) with |$w(x_0)$| harmonic in |$x_0 \in \mathbb{R}^2\backslash \overline{D}$|⁠, while on the left-hand side $$\begin{align}& \mathcal{R}_{\partial \varOmega}(u,v_n) = \mathcal{R}_{\partial D}(u,v_n) = \int_{\partial D}(1-k)(\partial_\nu u^- )v_n\,\textrm{d}y. \end{align}$$(3.22) If we assume the existence of a bounded sequence |$\{v_n\} \subset H^{1/2}(\partial D)$|⁠, then up to subsequences, by weak compactness we can extract |$v_n$| converging weakly to some |$f$| in |$H^{1/2}(\partial D)$|⁠, yielding $$\begin{align}& \mathcal{R}_{\partial \varOmega}(u,v_n) \to \int_{\partial D}(1-k)(\partial_\nu u^- )f\,\textrm{d}y = \widetilde{w}(x_0). \end{align}$$(3.23) Since both |$\partial _{\vec{d}}\,u(z,x_0) + w(x_0)$| and |$\widetilde{w}(x_0)$| are bounded harmonic functions in |$\mathbb{R}^2\backslash \overline{B}$|⁠, coinciding on |$\partial B$|⁠, by uniqueness and continuation principle, we can conclude that they coincide on |$\mathbb{R}^2\backslash (D\cup \{z\})$| arriving at a contradiction, by letting |$x_0\to z$|⁠. 4. Numerical examples In this section, we illustrate some numerical results. As before, we denote by |$v$| the solution of the Laplace equation (2.2) in |$\varOmega $| with Neumann data |$f$|⁠, and by |$u$| the solution of the inclusion equation (2.1), with corresponding boundary maps $$\begin{align} &{N_0}: H^{-1/2}_0(\partial \varOmega) \to H^{1/2}_0(\partial\varOmega) && {N_0} f = v|_{\partial\varOmega}, \end{align}$$(4.1) $$\begin{align} &{N_D}: H^{-1/2}_0(\partial \varOmega) \to H^{1/2}_0(\partial\varOmega) && {N_D}f = u|_{\partial\varOmega}. \end{align}$$(4.2) Both the algorithms, the linear sampling method and the reciprocity gap method, end up with a linear equation. The first method yields an equation of the form $$\begin{align}& (N_D - N_0)f = \psi_{0z}\quad\textup{ on} \partial\varOmega, \end{align}$$(4.3) where |$\psi _{0z} \coloneqq \vec{\varPsi }(x,z)\cdot \vec{d} - m_z - N_0\bigl (\partial _\nu \vec{\varPsi }(x,z) \cdot \vec{d}\bigr )$|⁠, with |$m_z\in \mathbb{R}$| a constant such that |$\psi _{0z}\in H^{1/2}_0(\partial \varOmega )$| has a vanishing mean. Furthermore, we can observe that |$\partial _\nu \psi _{0z} = 0$| on |$\partial \varOmega $|⁠. The method considered in this work yields the equation $$\begin{align}& \mathcal{R}(u_{x_0},v)=Rv=r_z\quad\textup{ for } x_0 \in \partial B, \end{align}$$(4.4) where |$r_z(x_0)\coloneqq \mathcal{R}\bigl (u_{x_0},\varPsi _z\bigr )$| for any |$u_{x_0}\in \mathcal{U}$| defined in (2.4). It is possible to recover this formulation from the first one, computing the duality with the generic function |$\partial _\nu u_{x_0}$|⁠. Since both the linear bounded operators |$N_D - N_0$| and |$R$| are compact, we regularize the linear equations according with the Tikhonov regularization combined with the Morozov discrepancy principle, as done in Kirsch (1998) for the scattering problem. Writing the previous equations in the form |$Kx=y$|⁠, the Tikhonov regularization consists in determining the regularized solution |$x^\alpha $| such that it minimizes the ‘Tikhonov functional’ $$\begin{align}& J_\alpha(x^\alpha)=\min\limits_{x\in X}J_\alpha(x) = \min\limits_{x\in X}\Big\{ \|Kx - y\|^2_Y + \alpha\|x\|^2_X \Big\}. \end{align}$$(4.5) In general, we could expect an error level in the data, which we fix of the order of |$\delta $| $$\begin{align}& \|y - y^\delta\|<\delta. \end{align}$$(4.6) We denote by |$\varDelta (\alpha ), \varDelta _N(\alpha )$|⁠, the ‘discrepancy’ function and its normalized variant $$\begin{align} \varDelta(\alpha)&\coloneqq\|Kx^{\alpha,\delta} - y^\delta\| - \delta, \end{align}$$(4.7) $$\begin{align} \varDelta_N(\alpha)&\coloneqq\|Kx^{\alpha,\delta} - y^\delta\| - \delta\|x^{\alpha,\delta}\|. \end{align}$$(4.8) Proposition 4.1 Let |$\delta $| be fixed, then there exists a unique zero |$\alpha ^*(\delta )$| of the discrepancy function, such that |$\varDelta \bigl (\alpha ^*\bigr )=0$|⁠. The same is still true for |$\varDelta _N(\alpha )$| (see Kirsch, 1998). Therefore, let |$\delta $| be fixed, then there exists a unique couple |$(\alpha (\delta ), x^{\alpha (\delta )})$| such that $$\begin{equation} \begin{cases} \alpha(\delta) x^{\alpha(\delta)} + K^*Kx^{\alpha(\delta)} = K^*y^\delta, \\ \varDelta(\alpha(\delta))=0, \quad\textup{ or} \quad \varDelta_N(\alpha(\delta))=0. \end{cases} \end{equation}$$(4.9) From above, we deduce that |$x^{\alpha (\delta )} = R_{\alpha (\delta )}y^\delta $| is the regularized solution computed from |$y^\delta $|⁠, and the Morozov principle guarantees us the convergence of the regularization strategy (see Kirsch, 2011, for the detailed proof) $$\begin{align}& \lim_{\delta\to0} R_{\alpha(\delta)}y^\delta = x, \quad \textrm{where}\ \|y^\delta-Kx\|<\delta. \end{align}$$(4.10) We will provide a numerical reconstruction of the unknown inclusion by the following function $$\begin{align}& r(z) = \frac{\|r_z\|}{\|v_z\|}, \end{align}$$(4.11) the ratio of the norms of the discrete approximations of the corresponding terms in (4.4). This will be divergent for |$z\notin D$|⁠, according to Theorem 2.4. Throughout these numerical examples, we fix the geometry of regions |$B$| and |$\varOmega $|⁠: a circumference with radius 3 for the observation boundary |$\partial \varOmega $|⁠, and a little larger circle with radius 3.5 for the region |$B$| where there are defined test functions. In this case, we observed that larger values for the radius of |$\partial B$| with respect to |$\partial \varOmega $| can seriously affect the resolution of the final result, since the fast decreasing of the fundamental solutions employed. We fixed the discretization of |$\partial B$|⁠, |$\partial \varOmega $| around about 200 points. The solution of the direct problem |$u_{x_0}(x)=\varPhi _D(x, x_0)\in \mathcal{U}$|⁠, for |$x_0\in \partial B$|⁠, can be represented as |$u_{x_0} = \varPhi _{x_0} + \mathcal{S}(\partial D, \psi ).$| Using this fact, we can obtain a linear integral equation from the inclusion problem (2.1), to solve for the density |$\psi $| the equation $$\begin{align}& \left(K^{\prime} + \frac{\hat{c}(k)}{2} I\right)\psi = -\partial_\nu \varPhi_{x_0}\quad\textrm{on}\ \partial D, \end{align}$$(4.12) with |$\hat{c}(k)\coloneqq (k+1)/(k-1), k\neq 1, k>0$| such that the operator |$K^{\prime} + \hat{c}(k)/2 I$| is injective, as highlighted in the previous section. We report some implementation details about the boundary element method employed for integrals discretization (see Kress, 1999, for more details). Let |$\partial D$| (or any boundaries) be a curve in |$\mathbb{R}^2$| of class |$C^2$| parametrized by |$ c(t):[0,2\pi ]\to \mathbb{R}^2, $| then in the two-dimensional case, the layer operators |$S, K, K^{\prime}$| defined by layer potentials (3.1)–(3.2) with |$x=c(t),\, y=c(\tau )$|⁠, can be rewritten as $$\begin{align}& (S\sigma)(t)\coloneqq\int_{0}^{2\pi}M(\tau, t)\sigma(c(\tau))\,\textrm{d}\tau, \end{align}$$(4.13) $$\begin{align}& (K\mu)(t)\coloneqq\int_{0}^{2\pi}L(\tau, t)\mu(c(\tau))\,\textrm{d}\tau. \end{align}$$(4.14) In our case, for |$C^2$| curves, the function |$L(t, \tau )$| is continuous for |$\tau \to t$|⁠, $$\begin{align*} &L(\tau,t)=\frac{1}{2\pi}\frac{c(t)-c(\tau)}{|c(t) - c(\tau)|^2}\cdot\nu(c(\tau))|c^{\prime}(\tau)|, \quad \lim_{\tau \to t}L(\tau,t) = -\frac{1}{4\pi}\kappa(t)|c^{\prime}(t)|, \end{align*}$$ where |$\kappa (t)$| is the curvature. The integrand function |$M(\tau , t)$| of the single layer has a logarithmic singularity, and the integral can be computed with an appropriate numerical quadrature rule for the second term in the sum, i.e. $$\begin{align}& \begin{cases} M_1(t,\tau)= -1/(4\pi), \\ M_2(t,\tau)=M(t,\tau) - M_1(t,\tau)\ln\Big(4\sin^2\dfrac{t-\tau}{2}\Big), & t\neq \tau. \end{cases} \end{align}$$(4.15) The application of the Nyström method guarantees the convergence of the solution of the discretized equation (4.12). In these cases, the quadrature rule derived from the trapezoidal rule on a uniform mesh, which converges exponentially for analytic density functions, while integrals defined on boundaries with corners can be computed on a graded mesh (a detailed implementation in the scattering context can be found in Colton & Kress, 1998). In Fig. 1, we compare the numerical reconstructions for the two methods. To obtain appreciable results, we need smaller values for the regularization parameter |$\alpha $| in the reciprocity gap method. This consideration is strengthened by the condition numbers of the two discretized maps |$N_r = N_0 - N_D$| and |$R$|⁠, (we report the estimates for the elliptic inclusion, they usually can differ about 1-2 orders of magnitude) $$\begin{align}& \kappa_2(N_r) = 1.9\mathrm{E+}{18}, \quad \kappa_2(R) = 1.1\mathrm{E+}{19}. \end{align}$$(4.16) Fig. 1. Open in new tabDownload slide Reconstruction with the linear sampling method in the upper figures, and the reciprocity gap functional method in the lower ones. The regularization parameter is equal to |$\alpha =\mathrm{1e}{-8}$| in the first two columns (two cases: the elliptic and the kite shape), and then it has been decreased |$\alpha =\mathrm{1e}{-12}$| in the third column, as required by the reciprocal gap method. Fig. 1. Open in new tabDownload slide Reconstruction with the linear sampling method in the upper figures, and the reciprocity gap functional method in the lower ones. The regularization parameter is equal to |$\alpha =\mathrm{1e}{-8}$| in the first two columns (two cases: the elliptic and the kite shape), and then it has been decreased |$\alpha =\mathrm{1e}{-12}$| in the third column, as required by the reciprocal gap method. The choice of |$\alpha $| is critically related to the size of the unknown inclusion, since this parameter can be used to tune the resolution of the reconstruction, and a too large parameter would result in a regularized but not accurate shape, as it can be seen in Fig. 2. Fig. 2. Open in new tabDownload slide Reciprocity gap (RG) method for |$\theta =\textup{Arg}(\vec{d}\,)=\pi /4$|⁠, and two different values for |$\alpha $| to be set according with the inclusion size: |$\alpha =\mathrm{1e}{-12}$| in the three figures above, |$\alpha =\mathrm{1e}{-15}$| below, where we identify the shape of the small inclusion. Fig. 2. Open in new tabDownload slide Reciprocity gap (RG) method for |$\theta =\textup{Arg}(\vec{d}\,)=\pi /4$|⁠, and two different values for |$\alpha $| to be set according with the inclusion size: |$\alpha =\mathrm{1e}{-12}$| in the three figures above, |$\alpha =\mathrm{1e}{-15}$| below, where we identify the shape of the small inclusion. Fig. 3. Open in new tabDownload slide RG method for |$\theta =\textup{Arg}(\vec{d}\,)=0$|⁠, and |$\alpha =\mathrm{1e}{-14}$|⁠. Fig. 3. Open in new tabDownload slide RG method for |$\theta =\textup{Arg}(\vec{d}\,)=0$|⁠, and |$\alpha =\mathrm{1e}{-14}$|⁠. In Figs 3–4, we highlight the dependence on the angle |$\theta $| corresponding to the unit vector |$\vec{d}$|⁠, whose effect is more relevant for large inclusions and not negligible regularization, but the overall region does not change much, and it would be redundant to consider all possible solutions with different |$\theta $|⁠. We used two different values of |$\alpha =0, \pi /4$|⁠, to study the behaviour for symmetric shapes. We provide the reconstructions of different unknown inclusions, to give a wide overview of the possible results, taking into account thin or small objects, corners, non-convex regions and multiple inclusions. We can observe that the method identifies correctly the shape for the single inclusions, and it is quite flexible about the choice of the amount of regularization. We cannot completely distinguish disconnected components, but we can make some guesses trying to decrease the parameter |$\alpha $| in these cases. Fig. 4. Open in new tabDownload slide RG method for |$\theta =\textup{Arg}(\vec{d}\,)=\pi /4$|⁠, and |$\alpha =\mathrm{1e}{-14}$|⁠. Fig. 4. Open in new tabDownload slide RG method for |$\theta =\textup{Arg}(\vec{d}\,)=\pi /4$|⁠, and |$\alpha =\mathrm{1e}{-14}$|⁠. Fig. 5. Open in new tabDownload slide Linear sampling method (LSM) for |$\theta =\textup{Arg}(\vec{d}\,)=\pi /4$|⁠, and |$\alpha =\mathrm{1e}{-12}$|⁠, an equivalent amount of regularization for LSM. Fig. 5. Open in new tabDownload slide Linear sampling method (LSM) for |$\theta =\textup{Arg}(\vec{d}\,)=\pi /4$|⁠, and |$\alpha =\mathrm{1e}{-12}$|⁠, an equivalent amount of regularization for LSM. There is a further aspect in this numeric discretization of |$Rv=r_z$| with the boundary elements method that combined with Tikhonov regularization seems to slightly suffer with respect to the same kind of discretization and Tikhonov regularization of the linear sampling method as we can see in Fig. 5, which show less spiky shapes (e.g. with large triangles). Furthermore, we obtained a positive feedback from the algorithm, testing it on two inclusions with different conductivities, like two ellipses or two kite shapes, as showed in Fig. 6. So far we reported images obtained for a conductivity |$k=2$|⁠, indeed different values lead to very similar results; the only fact we may notice is that for larger values |$k\sim 10,100$| the reconstruction has comparable results with smaller values of |$\alpha $| than one obtained for |$k=2$|⁠. Fig. 6. Open in new tabDownload slide Here we can see the reconstruction of two kete shapes with conductivities |$k_1=2 \,\textrm{and}\, k_2=2$|⁠, |$k_1=2\,\textrm{and}\,k_2=100$|⁠, |$k_1=100\,\textrm{and}\, k_2=100$|⁠, where |$k_2$| is the conductivity of the small inclusion. We can see that the first and the third simulations provide very similar results, and the second simulation is not too much affected (the same results can be seen with other couple of shapes). Fig. 6. Open in new tabDownload slide Here we can see the reconstruction of two kete shapes with conductivities |$k_1=2 \,\textrm{and}\, k_2=2$|⁠, |$k_1=2\,\textrm{and}\,k_2=100$|⁠, |$k_1=100\,\textrm{and}\, k_2=100$|⁠, where |$k_2$| is the conductivity of the small inclusion. We can see that the first and the third simulations provide very similar results, and the second simulation is not too much affected (the same results can be seen with other couple of shapes). Fig. 7. Open in new tabDownload slide In this figure we fix two values for |$\bar{\delta }$|⁠, equal to 0.03 in the first row and 0.05 in the second. In the first row we can see the first two simulations with the same noise and |$\alpha$| equals to |$1\mathrm{e}{-12}$| and |$1\mathrm{e}{-14}$| respectively. The other two pictures are two more examples with |$\alpha =1\mathrm{e}{-12}$|⁠. In the second row the two pairs of reconstructions are done with the same random noise and |$\alpha$| equals to |$1\mathrm{e}{-12}$| and |$1\mathrm{e}{-14}$| respectively. It appears celarly that the recontructions are strongly effected by the noise. Fig. 7. Open in new tabDownload slide In this figure we fix two values for |$\bar{\delta }$|⁠, equal to 0.03 in the first row and 0.05 in the second. In the first row we can see the first two simulations with the same noise and |$\alpha$| equals to |$1\mathrm{e}{-12}$| and |$1\mathrm{e}{-14}$| respectively. The other two pictures are two more examples with |$\alpha =1\mathrm{e}{-12}$|⁠. In the second row the two pairs of reconstructions are done with the same random noise and |$\alpha$| equals to |$1\mathrm{e}{-12}$| and |$1\mathrm{e}{-14}$| respectively. It appears celarly that the recontructions are strongly effected by the noise. It is well known in literature how the results obtained from the inversion of simulated data, like in our case, can lead to better estimates than those we should expect in real cases. Therefore, we studied how this algorithm is numerically stable when we introduce some normal distributed noise in the map that we simulated in the direct problem. We added this noise to the discretized version of the solutions |$u_{x_i}^n\in \mathbb{R}^n$| with |$x_i\in \partial B$| and |$i=1,\dots ,n$|⁠, that represents the physical quantity in real cases, i.e. $$\begin{align}& u_{x_i}^n + \delta^n_{x_i}, \quad \delta^n_{x_i}=\bar{\delta} \max_{\{x_i\in \partial B\}} |u^n_{x_i}|_{\infty} \cdot \mathcal{N}_{n}(0_n,\mathbb{I}_{n}). \end{align}$$(4.17) For the cases with |$\bar{\delta }=0.01, 0.03$|⁠, corresponding to noise equals to 1% and 3%, the method proved not to suffer too much, while for |$\delta =0.05$| the inclusion is correctly localized, and in most cases the size is identified, but whenever we try to increase the resolution diminishing |$\alpha $|⁠, the shape can be seriously affected. With noise |$\bar{\delta }=0.1$| most of the shapes are not correctly recognized, even if the position and partially the size can be deduced. In Fig. 8, we propose the results obtained by the Morozov discrepancy principle to establish a convergence criterion for the computation of the optimal regularization |$\alpha $|⁠, by setting the precision |$\delta $|⁠. The optimal value can be found for each testing point. The results obtained are qualitatively very similar to the previous ones with the simple reciprocal gap method in Fig. 4. In this case, we have the advantage that |$\alpha $| is automatically computed. Fig. 8. Open in new tabDownload slide Morozov discrepancy principle applied for same geometries contained in Fig. 4. In this case, the regularization parameter |$\alpha $| has been automatically computed for each testing point, imposing the vanishing of normalized discrepancy function |$\varDelta _N$| with |$\delta =\mathrm{5e}{-8}$|⁠. In these simulations, |$\alpha $| has been varied into a discrete set of values, which is the reason for the irregular contour lines. Fig. 8. Open in new tabDownload slide Morozov discrepancy principle applied for same geometries contained in Fig. 4. In this case, the regularization parameter |$\alpha $| has been automatically computed for each testing point, imposing the vanishing of normalized discrepancy function |$\varDelta _N$| with |$\delta =\mathrm{5e}{-8}$|⁠. In these simulations, |$\alpha $| has been varied into a discrete set of values, which is the reason for the irregular contour lines. 5. Conclusions From the theoretical perspective, the reciprocity gap method turns out to be equivalent to the linear sampling method, as it is the variational formulation of it. A possible disadvantage could be the additional integration that makes its computational approximation more difficult. In both cases, the main difficulty, which arises in the algorithm is the identification of an approximating sequence for a generic sampling point |$z$| inside the unknown domain |$D$|⁠. The regularization helps the stability of the reciprocity gap method, and the Morozov principle has the advantage of the automatic computation of the regularization parameter. On the other hand, this method is more stable to significant noise introduced in the discretized map |$R$|⁠, compared to other methods present in literature, like the factorization method. Further developments of this work could consider non constant conductivities, where we expect similar results. References Alessandrini , G. & Di Cristo , M. ( 2005 ) Stable determination of an inclusion by boundary measurements . SIAM J. Math. Anal. , 37 , 200 – 217 . Google Scholar Crossref Search ADS WorldCat Andrieux , S. & Ben Abda , A. ( 1996 ) Identification of planar cracks by complete overdetermined data: inversion formulae . Inverse Probl. , 12 , 553 – 563 . Google Scholar Crossref Search ADS WorldCat Bogomolny , A. ( 1985 ) Fundamental solutions method for elliptic boundary value problems . SIAM J. Numer. Anal. , 22 , 644 – 669 . Google Scholar Crossref Search ADS WorldCat Cakoni , F. , Di Cristo , M. & Sun , J. ( 2012 ) A multistep reciprocity gap functional method for the inverse problem in a multilayered medium . Complex Var. Elliptic Equ. , 57 , 261 – 276 . Google Scholar Crossref Search ADS WorldCat Cakoni , F. , Colton , D. & Haddar , H. ( 2001 ) The linear sampling method for anisotropic media . J. Comput. Appl. Math. , 146 , 285 – 299 . Google Scholar Crossref Search ADS WorldCat Cakoni , F. , Fares , M. B. & Haddar , H. ( 2006 ) Analysis of two linear sampling method applied to electromagnetic imagining of buried objects . Inverse Probl. , 22 , 845 – 867 . Google Scholar Crossref Search ADS WorldCat Calderón , A. P. ( 1980 ) On an inverse boundary value problem . Seminar on Numerical Analysis and its Applications to Continuum Physics . Rio de Janeiro : Soc. Brasileira de Matematica . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Colton , C. & Haddar , H. ( 2005 ) An application of the reciprocity gap functional to inverse scattering theory . Inverse Probl. , 21 , 383 – 398 . Google Scholar Crossref Search ADS WorldCat Colton , D. , Haddar , H. & Piana , M. ( 2003 ) The linear sampling method in inverse electromagnetic scattering theory . Inverse Probl. , 19 , S105 – S137 . Google Scholar Crossref Search ADS WorldCat Colton , D. & Kress , R. ( 1998 ) Inverse Acoustic and Electromagnetic Scattering Theory , 2nd edn. Applied Mathematical Sciences Series , vol. 93 . New York : Springer . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Di Cristo , M. ( 2007 ) Stable determination of an inhomogeneous inclusion by local boundary measurements . J. Comput. Appl. Math. , 198 , 414 – 425 . Google Scholar Crossref Search ADS WorldCat Di Cristo , M. & Rondi , L. ( 2003 ) Examples of exponential instability for inverse inclusion and scattering problems . Inverse Probl. , 19 , 685 – 701 . Google Scholar Crossref Search ADS WorldCat Di Cristo , M. & Sun , J. ( 2006 ) An inverse scattering problem for a partially coated buried obstacle . Inverse Probl. , 22 , 2331 – 2350 . Google Scholar Crossref Search ADS WorldCat Di Cristo , M. & Sun , J. ( 2007 ) The determination of the support and surface conductivity of a partially coated buried object . Inverse Probl. , 23 , 1161 – 1179 . Google Scholar Crossref Search ADS WorldCat Delbary , F. & Kress , R. ( 2010 ) Electrical impedance tomography with point electrodes . J. Integral Equations Appl. , 22 , 193 – 216 . Google Scholar Crossref Search ADS WorldCat Isakov , V. ( 1988 ) On uniqueness of recovery of a discontinuous conductivity coefficient . Comm. Pure Appl. Math. , 41 , 865 – 877 . Google Scholar Crossref Search ADS WorldCat Kirsch , A. ( 1998 ) Characterization of the shape of a scattering obstacle using the spectral data of the far field operator . Inverse Probl. , 14 , 1489 – 1512 . Google Scholar Crossref Search ADS WorldCat Kirsch , A. ( 2005 ) The factorization method for a class of inverse elliptic problems . Math. Nachr. , 278 , 258 – 277 . Google Scholar Crossref Search ADS WorldCat Kirsch , A. ( 2011 ) An Introduction to the Mathematical Theory of Inverse Problems , 2nd edn. Applied Mathematical Sciences Series, vol. 120 . New York : Springer . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Kress , R. ( 1999 ) Linear Integral Equations , 2nd edn. Applied Mathematical Sciences Series, vol. 82 . New York : Springer . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Monk , P. & Sun , J. ( 2007 ) Inverse scattering using finite elements and reciprocity gap . Inverse Probl. Imaging , 1 , 643 – 660 . Google Scholar Crossref Search ADS WorldCat Somersalo , E. ( 2001 ) Locating anisotropies in electrical impedance tomography , ArXiv 2001. https://arxiv.org/pdf/math/0110299.pdf . Sun , Y. , Guo , T. & Ma , F. ( 2016 ) The reciprocity gap functional method for the inverse scattering problem for cavities . Appl. Anal. , 96 , 1327 – 1346 . Google Scholar Crossref Search ADS WorldCat Zeng , F. , Liu , X., Sun , J. & Xu , L. ( 2017 ) Reciprocity gap method for an interior inverse scattering problem . J. Inv. Ill. Posed Probl. , 25 , 57 – 68 . Google Scholar OpenURL Placeholder Text WorldCat © The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Reconstruction of inclusions in electrical conductors
JF - IMA Journal of Applied Mathematics
DO - 10.1093/imamat/hxaa030
DA - 2020-11-25
UR - https://www.deepdyve.com/lp/oxford-university-press/reconstruction-of-inclusions-in-electrical-conductors-L72Ufyva9v
SP - 933
EP - 950
VL - 85
IS - 6
DP - DeepDyve
ER -