# A combination of Kohn-Vogelius and DDM methods for a geometrical inverse problem

Slim Chaabane, Houssem Haddar, Rahma Jerbi

<figure><img src="image_1.png" /></figure>

## To cite this version:

Slim Chaabane, Houssem Haddar, Rahma Jerbi. A combination of Kohn-Vogelius and DDM methods for a geometrical inverse problem. Inverse Problems, 2023, 39 (9), pp.095001. 10.1088/1361-6420/ace64a. hal-03940961v2

HAL Id: hal-03940961

https://inria.hal.science/hal-03940961v2

Submitted on 27 Sep 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A combination of Kohn-Vogelius and DDM methods for a geometrical inverse problem

S. Chaabane¹, H. Haddar²,*, R. Jerbi¹

¹Laboratory LAMHA-LR 11ES52, Department of Mathematics, Faculty of Sciences of Sfax, University of Sfax, BP 1171, Sfax 3000, Tunisia.

²INRIA, UMA, ENSTA Paris, Institut Polytechnique de Paris, 828, Bd des Maréchaux, 91762, Palaiseau Cedex, France.

*Corresponding Author

## Abstract

We consider the inverse geometrical problem of identifying the discontinuity curve of an electrical conductivity from boundary measurements. This standard inverse problem is used as a model to introduce and study a combined inversion algorithm coupling a gradient descent on the Kohn-Vogelius cost functional with a domain decomposition method that includes the unknown curve in the domain partitioning. We prove the local convergence of the method in a simplified case and numerically show its efficiency for some two dimensional experiments.

## 1 Introduction

The goal of this paper is to introduce and study a combined iterative inversion method that uses incomplete solver for the direct problem at each iteration for the inverse problem. This type of approach has the advantage of reducing the iterative cost and speeding up the convergence rate. Many variants of these combinations have been proposed in the literature [19, 8, 18, 21, 16, 29, 6]. We develop here an algorithm in the vein of so-called one-shot iterative methods for optimization problems. We address two novel aspects for this type of methods. The first one is to use a combination of the Kohn-Vogelius energy functional and a non overlapping domain decomposition method as an iterative solver. The second one is that the inverse parameter is the geometry of an unknown object used as a part of the domain partitioning. We study all these aspects in the framework of the inverse conductivity problem, where one would like to identify the discontinuity curve Σ of an electric conductivity from boundary measurements. This is a classic inverse problem that has many applications and has been extensively studied in the literature (we refer to [25, 7, 22, 33, 26, 30]). Our goal here is not to address issues specific to this inverse problem but to rather use it as a toy model to illustrate the feasibility of the combined inversion scheme and study convergence in some simplified configurations.

The Kohn-Vogelius cost functional has been used in the solution of various inverse geometrical problems [31, 12, 4, 23]. It has been specifically applied to inverse conductivity problem in [31, 4, 1, 2]. This functional seems to provide better stability and precision as compared with a classical least squares cost functional [10, 1]. From the theoretical point of view, it also
has the advantage of being differentiable with respect to the discontinuity curve under less
restrictive smoothness assumptions [9, 11]. We shall review some results from the literature on
this functional and study its local convergence for the inverse conductivity problem in the case
of circular domains.

The implementation of this method requires the solution of two direct problems, one as-
sociated with Neumann data and the other one associated with Dirichlet data coming from
the measurements. In order to evaluate these solutions we employ a non overlapping domain
decomposition method. We employ the so-called Optimized Schwarz Method (OSM) where
communications at the interfaces of the domains are enforced through Robin type boundary
conditions [32, 27, 28, 15, 13]. We here study the case where the inverse problem unknown Σ
is part of the partitioning used in the domain decomposition method. We propose a combined
inverse algorithm where the Kohn-Vogelius cost functional is minimized using a gradient de-
scent scheme. At each gradient step, the exact solutions of the Neumann and Dirichlet problems
are approximated using only one or a few iterations of the OSM. Two difficulties arise in this
coupling. The first one is that the domain partitioning changes after a gradient descent step
which requires modifications in the OSM scheme. The second one is that the gradient cannot
be evaluated exactly and therefore a choice has to be made. One should either evaluate the
shape gradient with respect to incomplete exact solutions or propose an approximation of the
gradient of the exact cost functional. The first option would require the introduction of two
adjoint problems and therefore may render the method more costly. This is why we adopt the
second approach, that indeed lead to an incorrect gradient at first iterations, but this gradient
becomes closer to the exact one as the iterations number increases. Concerning the first issue
related to OSM, we solve it by first rewriting the OSM as an iterative scheme on the impedance
interface values (see (3.6)). These interface values are then transported on the updated geom-
etry Σ using the shape gradient (see step 5 of Algorithm 2 ). We explicit this scheme in the
case of star shaped interfaces, but the approach can be easily extended to other type or shape
parametrizations. We study and prove local convergence of the resulting algorithm in the very
simplified case where the geometry is circular and the inverse shape problem is the radius of the
inner circle. We then numerically investigate the effectiveness and the accuracy of this algorithm
in the case of star shaped domains. We show in particular that only one OSM iteration would
achieve a converge rate similar to classical gradient (where the solutions are evaluated exactly
at each iteration). Determining the optimal choice of OSM iteration number in order to have
the best convergence rate is an open issue. Numerical tests suggests that this optimal choice is
among the one using few number of OSM iterations. The current work addresses only the 2D
problem. The extension to 3D problems would not bring any major difficulties with respect to
the theoretical developments in the present work (considering the local convergence in the case
of spheres). Also, the extension of the numerical algorithm can follow a similar procedure but
indeed its implementation would be more involving. In fact, the true efficiency of the proposed
combined inversion method can only be revealed for large scale problems such as 3D problems.
This is the subject of an forthcoming work.

The article is organized as follows. We introduce the direct and inverse problem together with
the Kohn-Vogelius cost functional in Section 2. The combined inversion algorithm is presented
and studied in Section 3. We provide in particular a local convergence result in Section 3.3.
Section 4 is dedicated to some numerical experiments for testing the efficiency of the combined
algorithm and comparing with the classical one.
## 2 A model problem

Let $\Omega$ be a simply connected bounded domain of $\mathbb{R}^2$ with $C^{1,\beta}$ boundary $\Gamma := \partial\Omega$, $\beta \in ]0, 1[$. We denote by $\sigma: \Omega \to \mathbb{R}$ the electric conductivity of $\Omega$ assumed to be a piecewise constant function with a regular discontinuity surface $\Sigma$. More specifically, $\sigma \in S_{ad}$ where

$$ S_{ad} := \{ \sigma = \sigma_1\chi_{\Omega_1} + \sigma_2\chi_{\Omega_2}; \sigma_1 > 0; \sigma_2 > 0; \Omega_1 \subset \Omega; \Sigma = \partial\Omega_1 \text{ a } C^{4,\beta} \text{ Jordan curve}; \Omega_2 = \Omega \setminus \bar{\Omega}_1 \}, \\ \chi_{\Omega_i} \text{ denotes the characteristics function of a domain } \Omega_i. $$

We denote by $u \in H^1(\Omega)$ the electric potential which satisfies the following Neumann boundary value problem

$$ (N_{\sigma}) \left\{ \begin{array}{ll} -\mathrm{div}(\sigma\nabla u) = 0 & \text{in } \Omega, \\ \sigma \frac{\partial u}{\partial \nu} = \phi & \text{on } \Gamma, \end{array} \right. $$

where $\nu$ denotes the outward unit normal on $\Gamma$ (see Figure 1) and $\phi \in L^2(\Gamma)$ the current flux through $\Gamma$ that satisfies the compatibility condition: $\int_{\Gamma} \phi \, ds = 0$. To ensure uniqueness of the solution to problem $(N_{\sigma})$, we impose that

$$ \int_{\Sigma} u \, ds = 0. \qquad (2.1) $$

Indeed one can choose other conditions that ensures uniqueness of the solutions. The advantage of using (2.1) is that it can be naturally encoded in the domain decomposition scheme introduced later.

<figure><img src="image_3.png" /><figcaption>Figure 1: The domain $\Omega$ divided into two subdomains $\Omega_1$ and $\Omega_2$.</figcaption></figure>

**Remark 2.1.** All the following still hold if we assume that $\sigma_1$ and $\sigma_2$ are known bounded regular functions on all of $\Omega$ and are positive definite.

### 2.1 The inverse problem and the Kohn-Vogelius functional

We consider in this section the inverse shape problem that consists in recovering the discontinuity interface $\bar{\Sigma}$ of an electrical conductivity $\bar{\sigma} \in S_{ad}$ from the knowledge of the flux $\phi$ together with the potential $f = u_{\bar{\sigma}|_{\Gamma}}$, where $u_{\bar{\sigma}}$ is the solution of $(N_{\bar{\sigma}})$.

**Remark 2.2.** Notice that, in general, a single couple $(\phi, f)$ is not sufficient to uniquely determine the unknown parameter $\bar{\sigma}$ [20]. However, two linearly independent pairs of measurements $(\phi_1, f_1)$, $(\phi_2, f_2)$ are sufficient to define the piecewise constant conductivity $\sigma$. Lipschitz stability results for the inverse geometrical problem using Dirichlet to Neumann maps have been
established in [5] when the interface is polygonal and extended to polyhedral interfaces for 3D problems in [3].

Our adopted approach consists in transforming the inverse problem into an optimization one by constructing a cost function $J$ modeling the energy gap between the solution of the direct problem and the solution of the following Dirichlet problem: $v \in H^1(\Omega)$, such that

$$ (\mathcal{D}_\sigma) \begin{cases} -\mathrm{div}(\sigma\nabla v) = 0 & \text{in } \Omega, \\ v = f & \text{on } \Gamma. \end{cases} $$

More precisely, we define the Kohn-Vogelius cost function as

$$ J(\sigma) := \int_{\Omega} \sigma \nabla(u_{\sigma} - v_{\sigma}) \cdot \nabla(u_{\sigma} - v_{\sigma}) dx $$

where $u_{\sigma} \in H^1(\Omega)$ the solution of the Neumann problem ($\mathcal{N}_{\sigma}$) and $v_{\sigma} \in H^1(\Omega)$, the solution of the Dirichlet problem ($\mathcal{D}_{\sigma}$). Indeed, the solution $\bar{\sigma}$ of the inverse problem is a zero of the non negative functional $J$ (since $u_{\bar{\sigma}} = v_{\bar{\sigma}}$) and therefore is a minimizer of $J$.

To numerically minimize the function $J$, we shall use a gradient descent algorithm based on the shape derivative of $J$ with respect to the singularity surface $\Sigma$ of $\sigma$. The existence and expression of the shape derivative has been studied in [1] and we outline in the following the main related results.

Let $\sigma \in S_{ad}$ with $\Sigma = \bar{\Omega}_1 \cap \bar{\Omega}_2$, $h > 0$ and $\zeta : \mathbb{R}^2 \to \mathbb{R}^2$ a $C^{4,\beta}$ vector field such that $\zeta = 0$ in a neighborhood of the boundary $\Gamma$. Then, there exists $h_0 > 0$, such that for all $h < h_0$, the mapping $\mathcal{F}_h = \mathrm{Id} + h\zeta$ is a $C^{4,\beta}$ diffeomorphism transforming the domain $\Omega$ into itself. We denote by $\Omega_{1,h} := \mathcal{F}_h(\Omega_1)$, $\Omega_{2,h} := \mathcal{F}_h(\Omega_2)$ and $\sigma_h = \sigma_1\chi_{\Omega_{1,h}} + \sigma_2\chi_{\Omega_{2,h}}$. Referring to [1, 2], the shape derivative of the cost function $J$ is given by the following theorem.

**Theorem 2.3.** The Kohn-Vogelius cost function $J$ is differentiable with respect to the shape $\Sigma$ and its derivative is given by:

$$ DJ(\sigma) \cdot \zeta = [\sigma] \int_{\Sigma} \left[ \frac{1}{\sigma_1 \sigma_2} \left( \left| \sigma \frac{\partial v_{\sigma}}{\partial \nu} \right|^2 - \left| \sigma \frac{\partial u_{\sigma}}{\partial \nu} \right|^2 \right) + \left( |\nabla_{\tau} v_{\sigma}|^2 - |\nabla_{\tau} u_{\sigma}|^2 \right) \right] \zeta \cdot \nu \, ds, $$

where $\zeta$ is as defined above, $[\sigma] := \sigma_1 - \sigma_2$, $\nabla_{\tau}$ denotes the tangential gradient and $\nu$ is the normal on $\Sigma$ oriented to the exterior of $\Omega_1$. The shape derivative is to be understood in the sense that

$$ \frac{J(\sigma_h) - J(\sigma)}{h} = DJ(\sigma) \cdot \zeta + \epsilon(h) \quad \text{where} \quad \lim_{h \to 0} |\epsilon(h)| = 0. $$

## 2.2 The gradient algorithm in the case of starlike domains

As a preparatory step to the combined algorithm we explicit the gradient descent algorithm in the case of starlike interfaces $\Sigma$. Let $\mathcal{C}$ be the set of $\mathcal{C}^1$ piecewise Jordan curves of $\mathbb{R}^2$, $n \in \mathbb{N}^*$ and $R = (R_0, \dots, R_{n-1}) \in (\mathbb{R}_+^*)^n$. For $i = 0, 1, \dots, n$, we denote by $\theta_i := \frac{2\pi}{n}i$ and $M_i := (R_i \cos(\theta_i), R_i \sin(\theta_i))$. Let $\Sigma := \Sigma_R$ be the interface defined as the union of the $n$ following arcs, $i = 0, \dots, n-1$

$$ S_i := \{ M = \hat{M}_i(t) := (\hat{R}_i(t) \cos(\hat{\theta}_i(t)), \hat{R}_i(t) \sin(\hat{\theta}_i(t)), t \in [0, 1] \}, \quad (2.2) $$
where $\hat{R}_i(t) := tR_{i+1} + (1-t)R_i$ and $\hat{\theta}_i(t) := t\theta_{i+1} + (1-t)\theta_i$ (see Figure 2 for an illustration) and where for the notation convenience we have set $R_n = R_0$. We also set $S_n = S_0$. We define the interface operator $\mathcal{T}_n$ by:

$$ \mathcal{T}_n: (\mathbb{R}_+^*)^n \to \mathcal{C} \\ R \mapsto \Sigma_R = \mathcal{T}_n(R) := \bigcup_{i=0}^{n-1} S_i. \qquad (2.3) $$

For $R \in (\mathbb{R}_+^*)^n$ sufficiently small, we denote by $\Omega_1 := \Omega_{1,R}$ the interior domain limited by the interface $\Sigma_R$ and by $\Omega_2 := \Omega_{2,R} = \Omega \setminus \overline{\Omega_{1,R}}$. The unknown of the inverse problem is $\bar{R} \in (\mathbb{R}_+^*)^n$ that corresponds with $\bar{\Sigma} = \Sigma_{\bar{R}}$.

<figure><img src="image_4.png" /><figcaption>Figure 2: The subdomain $\Omega_{1,R}$.</figcaption></figure>

Let us set $\sigma(R) := \sigma_1\chi_{\Omega_{1,R}} + \sigma_2\chi_{\Omega_{2,R}}$ and define the function $\mathcal{J}$ by

$$ \mathcal{J}: (\mathbb{R}_+^*)^n \to \mathbb{R} \\ R \mapsto \mathcal{J}(R) = J(\sigma(R)). \qquad (2.4) $$

The partial derivative $\frac{\partial \mathcal{J}}{\partial R_i}$ can be evaluated by applying Theorem 2.3 to a deformation field $\zeta = \zeta_i\nu$ on $\Sigma_R$ where $\zeta_i$ is hat function defined by, for $i = 0, \dots, n-1$,

$$ \zeta_i(M) = t\chi_{\{M=\hat{M}_{i-1}(t)\in S_{i-1}\}} + (1-t)\chi_{\{M=\hat{M}_i(t)\in S_i\}} $$

with $S_{-1} = S_{n-1}$, $\hat{M}_i(t)$ is defined in (2.2) and $\nu$ is the outward normal vector to $\Omega_1$. We then get, by using Theorem 2.3 that the derivative of the cost function $\mathcal{J}$ with respect to $R_i$ is given by the following formula:

$$ \frac{\partial \mathcal{J}}{\partial R_i}(R) = [\sigma(R)] \int_{\Sigma_R} \left[ \frac{1}{\sigma_1\sigma_2} \left( \left|\sigma(R)\frac{\partial v_{\sigma(R)}}{\partial \nu}\right|^2 - \left|\sigma(R)\frac{\partial u_{\sigma(R)}}{\partial \nu}\right|^2 \right) + \left( |\nabla_{\tau}v_{\sigma(R)}|^2 - |\nabla_{\tau}u_{\sigma(R)}|^2 \right) \right] \zeta_i ds. \qquad (2.5) $$

A gradient descent scheme to solve the inverse problem is summarized in Algorithm 1.
**Algorithm 1:** Gradient descent algorithm with exact direct solver

* Fix the number of parameters $n \in \mathbb{N}^*$ that serve to define the starlike interface.

* Consider an initial guess $R^0 \in (\mathbb{R}_+^*)^n$, the initial interface $\Sigma = \Sigma_{R^0} = \mathcal{T}_n(R^0)$ and the corresponding conductivity $\sigma(R^0) = \sigma_1\chi_{\Omega_{1,R^0}} + \sigma_2\chi_{\Omega_{2,R^0}}$ as defined above.

* $k = 0$.

repeat until $k \le \text{maximum number of iterations}$

* Use a direct solver to calculate $u_{\sigma(R^k)}$ and $v_{\sigma(R^k)}$, the respective solutions of $(\mathcal{N}_{\sigma(R^k)})$ and $(\mathcal{D}_{\sigma(R^k)})$.

* Calculate $\frac{\partial \mathcal{J}}{\partial R_i}(R^k)$, for $i = 0, \dots, n-1$ using formula (2.5).

* Update $\Sigma = \mathcal{T}_n(R^{k+1})$ with

$$R_i^{k+1} := R_i^k - \tau \frac{\partial \mathcal{J}}{\partial R_i}(R^k), \quad i = 0, \dots, n-1,$$

where $\tau > 0$ is chosen sufficiently small (a step adaptation can be incorporated here).

* $R^k = R^{k+1}$.

* $k = k + 1$.

end

As a first step for the analysis of the combined inverse scheme introduced later, we study the convergence of this algorithm in the particular case where both of the domains $\Omega$ and $\Omega_1$ are circular.

## 2.3 Local convergence analysis of the gradient descent algorithm

We study here the convergence of Algorithm 1 in the case $n = 1$, i.e $\Omega_{1,R}$ is an open disk of center $(0,0)$ and radius $R > 0$. We also choose the domain $\Omega$ to be the open disk of center $(0,0)$ and radius $R_2 > 0$. In this case the interface $\Sigma = \Sigma_R$ coincides with the circle of center $(0,0)$ and radius $R$ (note that $S_{-1} = S_0 = S_1$). The unknown of the inverse problem is $\bar{R}$ that corresponds with $\bar{\Sigma} = \Sigma_{\bar{R}}$. We impose the current flux $\phi(\theta) = m \cos(m\theta)$ or $\phi(\theta) = m \sin(m\theta)$, $\theta \in [0, 2\pi]$ and $m \in \mathbb{N}^*$. The solution of the direct problem $(\mathcal{N}_{\sigma(R)})$ then can be explicitly expressed as

$$u_{\sigma(R)}(r, \theta) = \begin{cases} u_1(r, \theta) = \alpha_N r^m \frac{\phi(\theta)}{m} & \text{in } \Omega_{1,R}, \\ u_2(r, \theta) = \left(\beta_N r^m + \frac{\gamma_N}{r^m}\right) \frac{\phi(\theta)}{m} & \text{in } \Omega_{2,R}, \end{cases}$$

where:

$$\left\{
\begin{aligned}
\alpha_N &:= \alpha_N(R) = \frac{2R_2^{m+1}}{\sigma_1(R_2^{2m} + R^{2m}) + \sigma_2(R_2^{2m} - R^{2m})}, \\
\beta_N &:= \beta_N(R) = \frac{(\sigma_2 + \sigma_1)R_2^{m+1}}{\sigma_2[\sigma_1(R_2^{2m} + R^{2m}) + \sigma_2(R_2^{2m} - R^{2m})]}, \\
\gamma_N &:= \gamma_N(R) = \frac{(\sigma_2 - \sigma_1)R_2^{m+1}R^{2m}}{\sigma_2[\sigma_1(R_2^{2m} + R^{2m}) + \sigma_2(R_2^{2m} - R^{2m})]}.
\end{aligned}
\right.
\quad (2.6)$$
The measurement f is then given by:

$$f(\theta) = C_f \frac{\phi(\theta)}{m}$$

where

$$C_f := \beta_N(\bar{R}) R_2^m + \frac{\gamma_N(\bar{R})}{R_2^m} = \frac{R_2 [\sigma_1(R_2^{2m} - \bar{R}^{2m}) + \sigma_2(R_2^{2m} + \bar{R}^{2m})]}{\sigma_2 [\sigma_1(R_2^{2m} + \bar{R}^{2m}) + \sigma_2(R_2^{2m} - \bar{R}^{2m})]} \quad (2.7)$$

Consequently, the solution of the Dirichlet problem ($D_{\sigma(R)}$) is

$$v_{\sigma(R)}(r, \theta) = \begin{cases} v_1(r, \theta) = \alpha_D r^m \frac{\phi(\theta)}{m} & \text{in } \Omega_{1,R}, \\ v_2(r, \theta) = \left(\beta_D r^m + \frac{\gamma_D}{r^m}\right) \frac{\phi(\theta)}{m} & \text{in } \Omega_{2,R}, \end{cases}$$

where:

$$\left\{
\begin{aligned}
\alpha_D &:= \alpha_D(R) = \frac{2\sigma_2 R_2^m C_f}{\sigma_1(R_2^{2m} - R^{2m}) + \sigma_2(R_2^{2m} + R^{2m})}, \\
\beta_D &:= \beta_D(R) = \frac{(\sigma_2 + \sigma_1) R_2^m C_f}{\sigma_1(R_2^{2m} - R^{2m}) + \sigma_2(R_2^{2m} + R^{2m})}, \\
\gamma_D &:= \gamma_D(R) = \frac{(\sigma_2 - \sigma_1) R_2^m R^{2m} C_f}{\sigma_1(R_2^{2m} - R^{2m}) + \sigma_2(R_2^{2m} + R^{2m})}.
\end{aligned}
\right.
\quad (2.8)$$

The function $\mathcal{J}$ depends here only on one variable $R$ and its derivative is given by (2.5) for $i=0$ with $\zeta_0 = 1$, namely,

$$\begin{aligned}
\mathcal{J}'(R) &= [\sigma(R)] \int_0^{2\pi} \left[ \frac{1}{\sigma_1 \sigma_2} \left( \left|\sigma(R) \frac{\partial v_{\sigma(R)}}{\partial \nu}\right|^2 - \left|\sigma(R) \frac{\partial u_{\sigma(R)}}{\partial \nu}\right|^2 \right) + \left( |\nabla_{\tau} v_{\sigma(R)}|^2 - |\nabla_{\tau} u_{\sigma(R)}|^2 \right) \right] R \, d\theta \\
&= \frac{\pi m^2 (\sigma_1^2 - \sigma_2^2)}{\sigma_2} R^{2m-1} (\alpha_D^2(R) - \alpha_N^2(R)).
\end{aligned}$$

Let us consider now the iterative sequence $R^k$ obtained by the gradient algorithm:

$$R^{k+1} = R^k - \tau \mathcal{J}'(R^k)$$

where $\tau > 0$ denotes the descent step of the gradient algorithm.

**Definition 2.4.** We say that the sequence $R^k$ is locally convergent, if there exists $\varepsilon > 0$, such that, for every $R^0 \in [\bar{R} - \varepsilon, \bar{R} + \varepsilon]$, the sequence $R^k$ converges to $\bar{R}$.

We then have the following proposition.

**Proposition 2.1.** Assume that $\sigma_1 \neq \sigma_2$ and let $\phi(\theta) = m \cos(m\theta)$ or $\phi(\theta) = m \sin(m\theta)$, $m \in \mathbb{N}^*$. Then $\mathcal{J}''(\bar{R}) > 0$ and the sequence $R^k$ is locally convergent if and only if the descent step $\tau < 2/\mathcal{J}''(\bar{R})$.
Proof. The function $\mathcal{J}$ is twice differentiable on $]0, R_2[$ and we have

$$ \mathcal{J}''(R) = F_1'(R)F_2(R) + F_1(R)F_2'(R) $$

where

$$ F_1(R) := \frac{\pi m^2 (\sigma_1^2 - \sigma_2^2)}{\sigma_2} R^{2m-1} \quad \text{and} \quad F_2(R) := \alpha_D^2(R) - \alpha_N^2(R). $$

From the fact $\alpha_N(\bar{R}) = \alpha_D(\bar{R})$, we deduce that

$$ \mathcal{J}''(\bar{R}) = F_1(\bar{R})F_2'(\bar{R}) = 2F_1(\bar{R})\alpha_N(\bar{R}) (\alpha'_D(\bar{R}) - \alpha'_N(\bar{R})), $$

where:

$$ \alpha'_N(\bar{R}) = \frac{4(\sigma_2 - \sigma_1)mR_2^{m+1}\bar{R}^{2m-1}}{[\sigma_1(R_2^{2m} + \bar{R}^{2m}) + \sigma_2(R_2^{2m} - \bar{R}^{2m})]^2} $$

and

$$ \alpha'_D(\bar{R}) = \frac{4(\sigma_1 - \sigma_2)m\sigma_2 C_f R_2^m \bar{R}^{2m-1}}{[\sigma_1(R_2^{2m} - \bar{R}^{2m}) + \sigma_2(R_2^{2m} + \bar{R}^{2m})]^2} $$

By using (2.7), we obtain

$$ \alpha'_D(\bar{R}) = \frac{4(\sigma_1 - \sigma_2)mR_2^{m+1}\bar{R}^{2m-1}}{[\sigma_1(R_2^{2m} - \bar{R}^{2m}) + \sigma_2(R_2^{2m} + \bar{R}^{2m})][\sigma_1(R_2^{2m} + \bar{R}^{2m}) + \sigma_2(R_2^{2m} - \bar{R}^{2m})]} $$

consequently

$$ F_2'(\bar{R}) = 2 \alpha_N(\bar{R}) (\alpha'_D(\bar{R}) - \alpha'_N(\bar{R})) = \frac{32m(\sigma_1^2 - \sigma_2^2)R_2^{4m+2}\bar{R}^{2m-1}}{[\sigma_1(R_2^{2m} - \bar{R}^{2m}) + \sigma_2(R_2^{2m} + \bar{R}^{2m})][\sigma_1(R_2^{2m} + \bar{R}^{2m}) + \sigma_2(R_2^{2m} - \bar{R}^{2m})]^3}, $$

and then

$$ \mathcal{J}''(\bar{R}) = \frac{32\pi m^3 (\sigma_1^2 - \sigma_2^2)^2 R_2^{4m+2} \bar{R}^{4m-2}}{\sigma_2 [\sigma_1(R_2^{2m} - \bar{R}^{2m}) + \sigma_2(R_2^{2m} + \bar{R}^{2m})][\sigma_1(R_2^{2m} + \bar{R}^{2m}) + \sigma_2(R_2^{2m} - \bar{R}^{2m})]^3} $$

The iterations for $R^k$ can be written as $R^{k+1} = g(R^k)$, where

$$ g(R) = R - \tau \mathcal{J}'(R). $$

Then, the sequence $R^k$ converges locally if and only if $|g'(\bar{R})| = |1 - \mathcal{J}''(\bar{R})| < 1$. From the condition: $\sigma_1 \neq \sigma_2$, we deduce that $\mathcal{J}''(\bar{R}) > 0$, and the sequence $R^k$ converges locally if and only if $0 < \tau < \frac{2}{\mathcal{J}''(\bar{R})}$. $\square$
# 3 A combined inversion method

We develop in this section some inversion algorithms combining the previous gradient algorithm with a Domain Decomposition Method (DDM) that respects the partitioning of the domain $\Omega$ into $\Omega_1 \cup \Omega_2 \cup \Sigma$. The main idea consists in approximating the direct problems ($\mathcal{N}_{\sigma(R^k)}$) and ($\mathcal{D}_{\sigma(R^k)}$) (at each iteration $k$ of Algorithm 1) using only one or a few DDM steps. Prior to giving the details of the resulting algorithm, we shall present first the so-called Optimized Schwarz Method (OSM) that we use as DDM.

## 3.1 Multidomain formulation using OSM

### 3.1.1 The case of the direct problem ($\mathcal{N}_{\sigma}$)

Let us present here the OSM for solving the direct problem ($\mathcal{N}_{\sigma}$). Let $u_i := u_{|\Omega_i}$ be the restriction of $u$ to $\Omega_i$, $i = 1, 2$. Then problem ($\mathcal{N}_{\sigma}$) can be reformulated as an equivalent multidomain problem consisting of the following subdomain problems

$$-\sigma_i \Delta u_i = 0 \quad \text{in } \Omega_i,$$

$$\sigma_2 \frac{\partial u_2}{\partial \nu} = \phi \quad \text{on } \Gamma,$$

together with the transmission conditions on the interface $\Sigma$

$$u_1 = u_2 \quad \text{on } \Sigma, \qquad (3.1)$$

$$\sigma_1 \frac{\partial u_1}{\partial \nu} = \sigma_2 \frac{\partial u_2}{\partial \nu} \quad \text{on } \Sigma. \qquad (3.2)$$

Alternatively and equivalently, one may impose the Robin transmission conditions,

$$\sigma_i \frac{\partial u_i}{\partial \nu} \pm \alpha u_i = \sigma_j \frac{\partial u_j}{\partial \nu} \pm \alpha u_j \quad \text{on } \Sigma \quad \text{for all } i = 1, 2 \text{ and } j \neq i, \qquad (3.3)$$

where $\alpha > 0$ is a fixed parameter that may be optimized to improve the convergence rate of the iterative domain decomposition method (see [32, 27, 28, 15, 13]). The resulting method is referred to as the Optimized Schwarz Method, which can be described as follows. Given arbitrary initial guesses ($u_i^0$)$_{1 \le i \le 2} \in H^2(\Omega_i)$, we inductively build the sequences $u_i^\ell \in H^1(\Omega_i)$, $i = 1, 2$, by solving (in parallel) for all $\ell \ge 0$

$$-\sigma_i \Delta u_i^{\ell+1} = 0 \quad \text{in } \Omega_i,$$

$$\sigma_2 \frac{\partial u_2^{\ell+1}}{\partial \nu} = \phi \quad \text{on } \Gamma,$$

$$\sigma_1 \frac{\partial u_1^{\ell+1}}{\partial \nu} + \alpha u_1^{\ell+1} = \sigma_2 \frac{\partial u_2^{\ell}}{\partial \nu} + \alpha u_2^{\ell} \quad \text{on } \Sigma, \qquad (3.4)$$

$$\sigma_2 \frac{\partial u_2^{\ell+1}}{\partial \nu} - \alpha u_2^{\ell+1} = \sigma_1 \frac{\partial u_1^{\ell}}{\partial \nu} - \alpha u_1^{\ell} \quad \text{on } \Sigma.$$

A direct use of (3.3) would require a numerical evaluation of the normal derivatives along the interfaces $\Sigma$ in order to compute the right-hand sides in the transmission conditions of (3.4). This can be avoided by renaming the problematic quantities

$$\lambda_{i,N}^{\ell} := \sigma_j \frac{\partial u_j^{\ell}}{\partial \nu} \pm \alpha u_j^{\ell} \quad \text{respectively for } i = 1, 2 \text{ and } j = 1, 2, j \neq i. \qquad (3.5)$$
where $\lambda_{i,N}^\ell$ is the information coming from the neighboring subdomain $\Omega_j$ ($j \neq i$) at step $\ell$ of the algorithm. One can easily verify that

$$ egin{cases} 
\lambda_{1,N}^{\ell+1} = \sigma_2 \frac{\partial u_2^{\ell+1}}{\partial \nu} + \alpha u_2^{\ell+1} = \sigma_2 \frac{\partial u_2^{\ell+1}}{\partial \nu} - \alpha u_2^{\ell+1} + 2\alpha u_2^{\ell+1}, \\ 
\lambda_{2,N}^{\ell+1} = \sigma_1 \frac{\partial u_1^{\ell+1}}{\partial \nu} - \alpha u_1^{\ell+1} = \sigma_1 \frac{\partial u_1^{\ell+1}}{\partial \nu} + \alpha u_1^{\ell+1} - 2\alpha u_1^{\ell+1}. 
\end{cases} $$

Therefore, the parameters $\lambda_{i,N}^\ell$ satisfy the induction

$$ \left\{ \begin{array}{ll} \lambda_{1,N}^{\ell+1} = \lambda_{2,N}^{\ell} + 2\alpha u_2^{\ell+1}, \\ \lambda_{2,N}^{\ell+1} = \lambda_{1,N}^{\ell} - 2\alpha u_1^{\ell+1} \end{array} \right. \qquad (3.6) $$

which can replace (3.5). We now can rewrite the iterative algorithm as follows. Given an initial guess $(\lambda_{i,N}^0)_{1 \le i \le 2} \in L^2(\Sigma)$ solve for each iteration $\ell \in \mathbb{N}$ the two following problems:

$$ -\sigma_i \Delta u_i^{\ell+1} = 0 \quad \text{in } \Omega_i, $$

$$ \sigma_2 \frac{\partial u_2^{\ell+1}}{\partial \nu} = \phi \quad \text{on } \Gamma, $$

$$ \sigma_1 \frac{\partial u_1^{\ell+1}}{\partial \nu} + \alpha u_1^{\ell+1} = \lambda_{1,N}^{\ell} \quad \text{on } \Sigma, \qquad (3.7) $$

$$ \sigma_2 \frac{\partial u_2^{\ell+1}}{\partial \nu} - \alpha u_2^{\ell+1} = \lambda_{2,N}^{\ell} \quad \text{on } \Sigma, $$

where the boundary values $\lambda_{i,N}^\ell$, $i = 1, 2$ are updated using (3.6). For the convergence analysis of this algorithm, we refer to [32, 17, 14]. Let us observe that condition (2.1) is ensured at convergence as soon as it is verified by the initial guess $(\lambda_{i,N}^0)_{1 \le i \le 2}$. This is what we summarize in the following lemma.

**Lemma 3.1.** Let $(\lambda_{i,N}^0)_{1 \le i \le 2} \in L^2(\Sigma)$, such that $\int_{\Sigma} \lambda_{i,N}^0 ds = 0$; $i = 1, 2$, then for $\ell \in \mathbb{N}^*$, we have: $\int_{\Sigma} \lambda_{i,N}^\ell ds = 0$, and $\int_{\Sigma} u_i^\ell ds = 0$.

*Proof.* We prove this Lemma by induction. By integrating equations (3.7) for $\ell = 0$ in $\Omega_i$ against a constant function in $\Omega_i$ we obtain

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Sigma} u_1^1 ds = -\frac{1}{\alpha} \left[ \int_{\Sigma} \lambda_{1,N}^0 ds \right] = 0, \\ \displaystyle \int_{\Sigma} u_2^1 ds = \frac{1}{\alpha} \left[ -\int_{\Sigma} \lambda_{2,N}^0 ds + \int_{\Gamma} \phi ds \right] = 0, \end{array} \right. $$

which proves the statement for $\ell = 1$. Assume that

$$ \int_{\Sigma} \lambda_{i,N}^{\ell} ds = \int_{\Sigma} u_{i}^{\ell} ds = 0, \text{ for } i = 1, 2. $$

By integrating equations (3.7) in $\Omega_i$ against a constant function in $\Omega_i$ we obtain

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Sigma} u_1^{\ell+1} ds = -\frac{1}{\alpha} \left[ \int_{\Sigma} \lambda_{1,N}^{\ell} ds \right] = 0, \\ \displaystyle \int_{\Sigma} u_2^{\ell+1} ds = \frac{1}{\alpha} \left[ -\int_{\Sigma} \lambda_{2,N}^{\ell} ds + \int_{\Gamma} \phi ds \right] = 0 \end{array} \right. $$
and by using (3.6), we obtain

$$ \int_{\Sigma} \lambda_{1,N}^{\ell+1} ds = \int_{\Sigma} \lambda_{2,N}^{\ell+1} ds = 0, $$

which prove the lemma.

<figure><img src="image_5.png" /></figure>

We study in the sequel the convergence rate of OSM in the case of circular interfaces, which will be useful for the convergence analysis of the combined inversion algorithm.

### 3.1.2 Convergence rate for OSM in circular domains

We consider again the particular case where the domain $\Omega_1$ is the open disk of center $(0, 0)$ and radius $R > 0$ of $\mathbb{R}^2$ and an $\Omega_2$ is the annulus domain $R < |x| < R_2$. The interface $\Sigma$ coincides with the circle of center $(0, 0)$ and radius $R$. We set for $m > 0$,

$$ k_m^{\alpha}(R) := \frac{\sigma_1 m - \alpha R}{\sigma_1 m + \alpha R}, \qquad (3.8) $$

$$ p_m^{\alpha}(R) := \frac{\sigma_2 m(R_2^{2m} - R^{2m}) - \alpha R(R_2^{2m} + R^{2m})}{\sigma_2 m(R_2^{2m} - R^{2m}) + \alpha R(R_2^{2m} + R^{2m})}. \qquad (3.9) $$

Then, we have the following proposition:

**Proposition 3.2.** For $\phi(\theta) = m \cos(m\theta)$ or $\phi(\theta) = m \sin(m\theta)$, $\theta \in [0, 2\pi]$ and $m \in \mathbb{N}^*$, the OSM (3.6)-(3.7), initialized with $\lambda_{i,N}^0 = 0$, $i = 1, 2$, geometrically converges for all $\alpha > 0$ and the convergence rate is given by the spectral radius

$$ \rho_m^{\alpha}(R) := \sqrt{|k_m^{\alpha}(R)p_m^{\alpha}(R)|} < 1. $$

*Proof.* For $\phi(\theta) = m \cos(m\theta)$ or $\phi(\theta) = m \sin(m\theta)$, $\theta \in [0, 2\pi]$ and $m \in \mathbb{N}^*$, the solutions of the problem (3.7) can be written as:

$$ u^{\ell+1}(r, \theta) = \begin{cases} u_1^{\ell+1}(r, \theta) = \alpha_N^{\ell+1} r^m \frac{\phi(\theta)}{m} & \text{in } \Omega_1, \\ u_2^{\ell+1}(r, \theta) = \left( \beta_N^{\ell+1} r^m + \frac{\gamma_N^{\ell+1}}{r^m} \right) \frac{\phi(\theta)}{m} & \text{in } \Omega_2, \end{cases} $$

where the constants $\alpha_N^{\ell+1}$, $\beta_N^{\ell+1}$ and $\gamma_N^{\ell+1}$ are determined by the following equations:

$$ \left\{ \begin{aligned} \sigma_2 \frac{\partial u_2^{\ell+1}}{\partial \nu} &= \phi(\theta) & \text{on } \Gamma, \\ \sigma_1 \frac{\partial u_1^{\ell+1}}{\partial \nu} + \alpha u_1^{\ell+1} &= \lambda_{1,N}^{\ell} & \text{on } \Sigma, \\ \sigma_2 \frac{\partial u_2^{\ell+1}}{\partial \nu} - \alpha u_2^{\ell+1} &= \lambda_{2,N}^{\ell} & \text{on } \Sigma. \end{aligned} \right. \qquad (3.10) $$

For $i = 1, 2$, we denote by $\hat{\lambda}_{i,N}^{\ell} = \lambda_{i,N}^{\ell} \frac{m}{\phi(\theta)}$. Then, from the second equation of (3.10), we obtain:

$$ u_1^{\ell+1}(r, \theta) = \frac{\hat{\lambda}_{1,N}^{\ell}}{R^{m-1}(m\sigma_1 + \alpha R)} r^m \frac{\phi(\theta)}{m} \qquad (3.11) $$
and from the first and the third equations of (3.10), we obtain:

$$u_2^{\ell+1}(r, \theta) = \left[ \frac{R_2^{m+1} R^{m+1} \left( (\sigma_2 m - \alpha R) R^{m-1} - \sigma_2 \hat{\lambda}_{2,N} R_2^{m-1} \right)}{\sigma_2 (\sigma_2 m (R_2^{2m} - R^{2m}) + \alpha R (R_2^{2m} + R^{2m}))} \left( \frac{r^m}{R_2^{2m}} + \frac{1}{r^m} \right) + \frac{r^m}{\sigma_2 R_2^{m-1}} \right] \frac{\phi(\theta)}{m}.$$

From (3.6) and using the fact that $\lambda_{i,N}^{\ell+1} = \hat{\lambda}_{i,N}^{\ell+1} \frac{\phi(\theta)}{m}$, we obtain:

$$\hat{\lambda}_{2,N}^{\ell+1} = k_m^{\alpha}(R) \hat{\lambda}_{1,N}^{\ell}.$$

In the same way we get:

$$\hat{\lambda}_{1,N}^{\ell+1} = p_m^{\alpha}(R) \hat{\lambda}_{2,N}^{\ell} + \eta_m^{\alpha}(R),$$

where

$$\eta_m^\alpha(R) := \frac{4m\alpha R^m R_2^{m+1}}{\sigma_2 m(R_2^{2m} - R^{2m}) + \alpha R(R_2^{2m} + R^{2m})} \quad (3.12)$$

Therefore, the DDM iterations can be written as:

$$\begin{pmatrix} \hat{\lambda}_{1,N}^{\ell+1} \\ \hat{\lambda}_{2,N}^{\ell+1} \end{pmatrix} = \begin{pmatrix} 0 & p_m^\alpha(R) \\ k_m^\alpha(R) & 0 \end{pmatrix} \begin{pmatrix} \hat{\lambda}_{1,N}^\ell \\ \hat{\lambda}_{2,N}^\ell \end{pmatrix} + \begin{pmatrix} \eta_m^\alpha(R) \\ 0 \end{pmatrix}. \quad (3.13)$$

This induction converges if and only if $\rho_m^\alpha(R) = \sqrt{|k_m^\alpha(R)p_m^\alpha(R)|} < 1$, which is always true for $\alpha > 0, m > 0$ and $0 < R < R_2$. $\square$

### 3.1.3 The case of the Dirichlet problem ($D_\sigma$)

Similarly to above, we also apply the OSM for solving the Dirichlet problem ($D_\sigma$). Let $v_i := v_{|\Omega_i}$ be the restriction of $v$ to $\Omega_i$, $i = 1, 2$. Then problem ($D_\sigma$) can be reformulated as an equivalent multidomain problem consisting of the following subdomain problems

$$\begin{array}{l} -\sigma_i \Delta v_i = 0 \quad \text{in } \Omega_i, \\ v_2 = f \quad \text{on } \Gamma, \end{array}$$

together with the Robin transmission conditions:

$$\sigma_i \frac{\partial v_i}{\partial \nu} \pm \alpha v_i = \sigma_j \frac{\partial v_j}{\partial \nu} \pm \alpha v_j \quad \text{on } \Sigma \text{ for all } i = 1, 2 \text{ and } j \neq i. \quad (3.14)$$

The OSM for this problem can be then formulated as: Given an initial guess $(\lambda_{i,D}^0)_{1\le i\le 2} \in L^2(\Sigma)$ solve for each iteration $\ell \in \mathbb{N}$ the two following problems:

$$\begin{array}{l} -\sigma_i \Delta v_i^{\ell+1} = 0 \quad \text{in } \Omega_i, \\ v_2^{\ell+1} = f \quad \text{on } \Gamma, \\ \\ \sigma_1 \frac{\partial v_1^{\ell+1}}{\partial \nu} + \alpha v_1^{\ell+1} = \lambda_{1,D}^{\ell} \quad \text{on } \Sigma, \\ \\ \sigma_2 \frac{\partial v_2^{\ell+1}}{\partial \nu} - \alpha v_2^{\ell+1} = \lambda_{2,D}^{\ell} \quad \text{on } \Sigma, \end{array} \quad (3.15)$$
where the parameters $\lambda_{i,D}^{\ell+1}$, $i=1,2$, verify the induction

$$ \left\{ \begin{array}{l} \lambda_{1,D}^{\ell+1} = \lambda_{2,D}^{\ell} + 2\alpha v_{2}^{\ell+1}, \\ \lambda_{2,D}^{\ell+1} = \lambda_{1,D}^{\ell} - 2\alpha v_{1}^{\ell+1}. \end{array} \right. \qquad (3.16) $$

We can state and prove similar convergence results as in Proposition 3.2. We set

$$ q_m^\alpha(R) := \frac{\sigma_2 m (R_2^{2m} + R^{2m}) - \alpha R (R_2^{2m} - R^{2m})}{\sigma_2 m (R_2^{2m} + R^{2m}) + \alpha R (R_2^{2m} - R^{2m})}. \qquad (3.17) $$

**Proposition 3.3.** For $f(\theta) = C_f \frac{\phi(\theta)}{m}$, $C_f \in \mathbb{R}$, $m \in \mathbb{N}^*$, the OSM (3.15)-(3.16) geometrically converges for all $\alpha > 0$ and the convergence rate is given by the spectral radius

$$ \hat{\rho}_m^\alpha(R) = \sqrt{|k_m^\alpha(R)q_m^\alpha(R)|} < 1. $$

*Proof.* For $f(\theta) = C_f \frac{\phi(\theta)}{m}$, $C_f \in \mathbb{R}$, the approximate solutions of the problem (3.15) can be written as:

$$ v^{\ell+1}(r, \theta) = \begin{cases} v_1^{\ell+1}(r, \theta) = \alpha_D^{\ell+1} r^m \frac{\phi(\theta)}{m} & \text{in } \Omega_1, \\ v_2^{\ell+1}(r, \theta) = \left( \beta_D^{\ell+1} r^m + \frac{\gamma_D^{\ell+1}}{r^m} \right) \frac{\phi(\theta)}{m} & \text{in } \Omega_2. \end{cases} $$

where the constants $\alpha_D^{\ell+1}$, $\beta_D^{\ell+1}$ and $\gamma_D^{\ell+1}$ are calculated by the following equations:

$$ \left\{ \begin{array}{ll} v_2^{\ell+1} = f(\theta) & \text{on } \Gamma, \\ \sigma_1 \frac{\partial v_1^{\ell+1}}{\partial \nu} + \alpha v_1^{\ell+1} = \lambda_{1,D}^{\ell} & \text{on } \Sigma, \\ \sigma_2 \frac{\partial v_2^{\ell+1}}{\partial \nu} - \alpha v_2^{\ell+1} = \lambda_{2,D}^{\ell} & \text{on } \Sigma. \end{array} \right. \qquad (3.18) $$

For $i=1,2$, we denote by $\hat{\lambda}_{i,D}^{\ell} = \lambda_{i,D}^{\ell} \frac{m}{\phi(\theta)}$. Then, from the second equation of (3.18), we obtain:

$$ v_1^{\ell+1}(r, \theta) = \frac{\hat{\lambda}_{1,D}^{\ell}}{R^{m-1}(m\sigma_1 + \alpha R)} r^m \frac{\phi(\theta)}{m}, \qquad (3.19) $$

and from the first and the third equation of (3.18), we get

$$ v_2^{\ell+1}(r, \theta) = \left[ \frac{R_2^m R^{m+1} \left( (\sigma_2 m - \alpha R) C_f R^{m-1} - \hat{\lambda}_{2,D}^{\ell} R_2^m \right)}{\sigma_2 m (R_2^{2m} + R^{2m}) + \alpha R (R_2^{2m} - R^{2m})} \left( \frac{1}{r^m} - \frac{r^m}{R_2^{2m}} \right) + \frac{C_f}{R_2^m} r^m \right] \frac{\phi(\theta)}{m}. $$

From (3.16) and using the fact that $\lambda_{i,D}^{\ell+1} = \hat{\lambda}_{i,D}^{\ell+1} \frac{\phi(\theta)}{m}$, we obtain

$$ \hat{\lambda}_{2,D}^{\ell+1} = k_m^\alpha(R) \hat{\lambda}_{1,D}^{\ell}. $$
In the same way we get

$$ \hat{\lambda}_{1,D}^{\ell+1} = q_m^\alpha(R) \hat{\lambda}_{2,D}^\ell + \omega_m^\alpha(R), $$

where

$$ \omega_m^\alpha(R) = \frac{4\alpha\sigma_2mC_fR^mR_2^m}{\sigma_2m(R_2^{2m} + R^{2m}) + \alpha R(R_2^{2m} - R^{2m})}. \quad (3.20) $$

Therefore, the DDM iterations can be written as:

$$ \begin{pmatrix} \hat{\lambda}_{1,D}^{\ell+1} \\ \hat{\lambda}_{2,D}^{\ell+1} \end{pmatrix} = \begin{pmatrix} 0 & q_m^\alpha(R) \\ k_m^\alpha(R) & 0 \end{pmatrix} \begin{pmatrix} \hat{\lambda}_{1,D}^\ell \\ \hat{\lambda}_{2,D}^\ell \end{pmatrix} + \begin{pmatrix} \omega_m^\alpha(R) \\ 0 \end{pmatrix}. \quad (3.21) $$

This induction converges if and only if $\hat{\rho}_m^\alpha(R) = \sqrt{|k_m^\alpha(R)q_m^\alpha(R)|} < 1$, which is always true for $\alpha > 0, m > 0$ and $0 < R < R_2$. $\square$

## 3.2 Combined inversion method

We now present the combined inversion algorithm that couples Algorithm 1 with OSM. Given an integer $L > 0$, an interface $\Sigma$ and an initial guess $\lambda_{i,N}^0 \in L^2(\Sigma)$, $i = 1, 2$ we denote by $N^L(\Sigma, \lambda_{1,N}, \lambda_{2,N})$ the $L^{th}$ iterate of (3.7). More precisely, we set

$$ \left\{ \begin{array}{ll} N^L(\Sigma, \lambda_{1,N}^0, \lambda_{2,N}^0) := u_1^L & \text{in } \Omega_1, \\ N^L(\Sigma, \lambda_{1,N}^0, \lambda_{2,N}^0) := u_2^L & \text{in } \Omega_2, \end{array} \right. $$

where $(u_1^\ell)$ and $(u_2^\ell)$, $\ell = 1, \dots, L$ verify the induction (3.7). We also set

$$ \Lambda_{i,N}^L(\Sigma, \lambda_{1,N}^0, \lambda_{2,N}^0) := \lambda_{i,N}^L \text{ on } \Sigma, \quad i = 1, 2, $$

where $\lambda_{i,N}^L$ is the iterate number $L$ of (3.7).

Similarly, for some initial guess $\lambda_{i,D}^0 \in L^2(\Sigma)$, $i = 1, 2$ we define:

$$ \left\{ \begin{array}{ll} D^L(\Sigma, \lambda_{1,D}^0, \lambda_{2,D}^0) := v_1^L & \text{in } \Omega_1, \\ D^L(\Sigma, \lambda_{1,D}^0, \lambda_{2,D}^0) := v_2^L & \text{in } \Omega_2, \end{array} \right. $$

where $(v_1^\ell)$ and $(v_2^\ell)$, $\ell = 1, \dots, L$ verify the induction (3.15) and by

$$ \Lambda_{i,D}^L(\Sigma, \lambda_{1,D}^0, \lambda_{2,D}^0) := \lambda_{i,D}^L \text{ on } \Sigma, \quad i = 1, 2, $$

where $\lambda_{i,D}^L$ is the iterate number $L$ of (3.15).

Roughly speaking, the combined algorithm consists in replacing $u_{\sigma(R^k)}$ and $v_{\sigma(R^k)}$ respectively with $N^L(\Sigma, \lambda_{1,N}^0, \lambda_{2,N}^0)$ and $D^L(\Sigma, \lambda_{1,D}^0, \lambda_{2,D}^0)$ where $\Sigma = \mathcal{T}_n(R^k)$ as in Algorithm 1. The main ambiguity in the construction of the induction is the update for $\lambda_{i,N}$ and $\lambda_{i,D}$ on $\Sigma = \mathcal{T}_n(R^{k+1})$ using the boundary values of the solutions at previous iterate. We choose to transport these values on the new update of the interface (i.e $\Sigma = \mathcal{T}_n(R^{k+1})$) using the same mapping that transforms $\mathcal{T}_n(R^k)$ into $\mathcal{T}_n(R^{k+1})$. We explicit this for starlike interfaces as in Algorithm 1 where we suppose that the number of subdivisions $n$ of the interface is kept fixed during iterations. A function $\lambda$ defined on $\Sigma = \Sigma_R$ can be parameterized on $S_j, j = 0, \dots, n-1$ as a function of the variable $t \in [0, 1]$ as follows

$$ \lambda(t) = \lambda(\hat{M}_j(t)) \quad (3.22) $$
where $\hat{M}_j(t)$ is given by (2.2). The gradient flow as defined by Algorithm 1 preserves the parametrization of $\Sigma$ in terms of $S_j$. Therefore the update for the boundary values on $\Sigma$ can be written easily in terms of the variable $t \in [0, 1]$ since the latter is independent from the iteration index.

There is also an ambiguity in defining the partial shape derivative given by (2.5) since the OSM iterates are not continuous across $\Sigma_R$ nor are the fluxes in general. One can either evaluate the shape gradient with respect to $N^L(\Sigma, \lambda_{1,N}^0, \lambda_{2,N}^0)$ and $D^L(\Sigma, \lambda_{1,D}^0, \lambda_{2,D}^0)$ or approximate the gradient of the exact cost functional using these approximate solutions of the direct problems. The first option would require the introduction of two adjoint problems and therefore may render the method more costly. This is why we adopt the second approach, that indeed lead to an incorrect gradient at first iterations, but this gradient becomes close to the exact one as the iterations number increases. More precisely, the gradient expression given by (2.5) requires the trace and the normal trace values of the solutions to the Dirichlet and Neumann problems which are continuous for the exact solutions. However, these values are no longer the same if we use the restriction of solutions to $\Omega_1$ or $\Omega_2$ coming from incomplete OSM iterations. In the following algorithm, we choose to apply (2.5) using the traces and the normal traces of $u_{\sigma(R)|\Omega_1}$ and $v_{\sigma(R)|\Omega_2}$.
Algorithm 2: A combined inversion algorithm

• Fix the number of parameters $n \in \mathbb{N}^*$ that serve to define the starlike interface.

• Consider an initial guess $R^0 \in (\mathbb{R}_+^*)^n$, the initial interface $\Sigma = \Sigma_{R^0} = \mathcal{T}_n(R^0)$ and the corresponding conductivity $\sigma(R^0) = \sigma_1\chi_{\Omega_{1,R^0}} + \sigma_2\chi_{\Omega_{2,R^0}}$ as defined above.

• Choose as initial boundary values:

$$ \lambda_{i,N}(\Sigma) = 0, \lambda_{i,D}(\Sigma) = 0 \text{ on } \Sigma, i = 1, 2. $$

• k=0;

repeat until $k \le \text{maximum number of iterations}$

1. Set

$$ \lambda_{i,N}^0 = \lambda_{i,N}(\Sigma), \lambda_{i,D}^0 = \lambda_{i,D}(\Sigma) \text{ on } \Sigma, i = 1, 2. $$

2. Use $L$ iterations of OSM to evaluate

$$ u_{\sigma(R^k)} = N^L(\Sigma, \lambda_{1,N}^0, \lambda_{2,N}^0), \lambda_{i,N}^L(R^k) = \Lambda_{i,N}^L(\Sigma, \lambda_{1,N}^0, \lambda_{2,N}^0). $$

$$ v_{\sigma(R^k)} = D^L(\Sigma, \lambda_{1,D}^0, \lambda_{2,D}^0), \lambda_{i,D}^L(R^k) = \Lambda_{i,D}^L(\Sigma, \lambda_{1,D}^0, \lambda_{2,D}^0). $$

3. Evaluate $\frac{\partial \mathcal{J}}{\partial R_j}(R^k)$, for $j = 0, \dots, n-1$ using formula (2.5) where the boundary values are calculated using

$$ u_{\sigma(R^k)|_{\Omega_1}} \text{ and } v_{\sigma(R^k)|_{\Omega_1}}. $$

4. Update $\Sigma = \mathcal{T}_n(R^{k+1})$ with

$$ R_j^{k+1} := R_j^k - \tau \frac{\partial \mathcal{J}}{\partial R_j}(R^k), \quad j = 0, \dots, n-1, $$

where $\tau > 0$ is chosen sufficiently small. (A step adaptation can be incorporated here but only after a few iterations, when the gradient becomes sufficiently accurate and provides a descent direction).

5. Update the interface values on $S_j$ as

$$ \lambda_{i,N}(\Sigma)(\hat{M}_j(t)) = \lambda_{i,N}^L(R^k)(t), \quad t \in [0, 1], $$

$$ \lambda_{i,D}(\Sigma)(\hat{M}_j(t)) = \lambda_{i,D}^L(R^k)(t), \quad t \in [0, 1], $$

following the identification (3.22).

6. $R^k = R^{k+1}$.

7. $k = k + 1$.

end

As indicated earlier, in step 3 of Algorithm 2 we could have made the choice to compute the derivative using $u_{\sigma(R^k)|_{\Omega_2}}$ and $v_{\sigma(R^k)|_{\Omega_2}}$. This would not affect the theoretical or numerical results below. We keep the choice made in Algorithm 2 in the remainder of this paper for the sake of conciseness.
### 3.3 Local convergence analysis of Algorithm 2 for circular domains

This section is dedicated to the convergence analysis of Algorithm 2 in the case of circular domains as in Section (2.3) and Section (3.1.2) and for $L = 1$, i.e we perform only one DDM step per iteration. With the notation of Section 2.3 and Algorithm 2, the iterative scheme leads to the following. The boundary values at iteration $k$ are of the form

$$ \lambda_{i,N}^{L}(R^{k}) = \hat{\lambda}_{i,N}^{L}(R^{k}) \frac{\phi(\theta)}{m}, \quad \lambda_{i,D}^{L}(R^{k}) = \hat{\lambda}_{i,D}^{L}(R^{k}) \frac{\phi(\theta)}{m}, $$

for some constants $\hat{\lambda}_{i,N}^{L}(R^{k})$ and $\hat{\lambda}_{i,D}^{L}(R^{k})$ and for $i = 1, 2$. Therefore according to equations (3.11) and (3.19), the solutions after one step DDM iterate are given in $\Omega_1$ by

$$ u_{\sigma(R^k)}(r, \theta) = \alpha_N^L(R^k) r^m \frac{\phi(\theta)}{m}, \quad v_{\sigma(R^k)}(r, \theta) = \alpha_D^L(R^k) r^m \frac{\phi(\theta)}{m}, $$

for $R^k < r < R_2$ with

$$ \left\{ \begin{array}{l} \alpha_N^L(R) := \frac{\hat{\lambda}_{1,N}^L(R)}{R^{m-1}(m\sigma_1 + \alpha R)}, \\ \alpha_D^L(R) := \frac{\hat{\lambda}_{1,D}^L(R)}{R^{m-1}(m\sigma_1 + \alpha R)}. \end{array} \right. $$

Consequently, by step 3 of Algorithm 2, the gradient is approximated by

$$ \frac{\partial \mathcal{J}}{\partial R}(R^k) \approx \frac{\pi m^2 (\sigma_1^2 - \sigma_2^2) R^k}{\sigma_2 (\sigma_1 m + \alpha R^k)^2} (|\hat{\lambda}_{1,D}^L(R^k)|^2 - |\hat{\lambda}_{1,N}^L(R^k)|^2). $$

The iterative scheme of Algorithm 2 now reads

$$ R^{k+1} = R^k - \tau \frac{\pi m^2 (\sigma_1^2 - \sigma_2^2) R^k}{\sigma_2 (\sigma_1 m + \alpha R^k)^2} (|\hat{\lambda}_{1,D}^L(R^k)|^2 - |\hat{\lambda}_{1,N}^L(R^k)|^2) \quad (3.23) $$

and the updates for the boundary values are given by

$$ \hat{\lambda}_{1,N}^{L}(R^{k+1}) = p_m^{\alpha}(R^k) \hat{\lambda}_{2,N}^{L}(R^k) + \eta_m^{\alpha}(R^k), \quad (3.24) $$

$$ \hat{\lambda}_{2,N}^{L}(R^{k+1}) = k_m^{\alpha}(R^k) \hat{\lambda}_{1,N}^{L}(R^k), \quad (3.25) $$

$$ \hat{\lambda}_{1,D}^{L}(R^{k+1}) = q_m^{\alpha}(R^k) \hat{\lambda}_{2,D}^{L}(R^k) + \omega_m^{\alpha}(R^k), \quad (3.26) $$

$$ \hat{\lambda}_{2,D}^{L}(R^{k+1}) = k_m^{\alpha}(R^k) \hat{\lambda}_{1,D}^{L}(R^k), \quad (3.27) $$

where $p_m^{\alpha}(R^k)$, $\eta_m^{\alpha}(R^k)$, $k_m^{\alpha}(R^k)$, $q_m^{\alpha}(R^k)$ and $\omega_m^{\alpha}(R^k)$ are respectively given by (3.9), (3.12), (3.8), (3.17) and (3.20).

If we set $X^k := (R^k, \hat{\lambda}_{1,N}^L(R^k), \hat{\lambda}_{2,N}^L(R^k), \hat{\lambda}_{1,D}^L(R^k), \hat{\lambda}_{2,D}^L(R^k))$, then the iterative scheme formed by (3.23)-(3.27) can be synthetically written as

$$ X^{k+1} = G(X^k), $$
where $G:]0, R_2[\times\mathbb{R}^4 \to \mathbb{R}^5$, $X = (x, y, z, t, h) \mapsto (g_1(X), g_2(X), g_2(X), g_2(X), g_2(X))$ is given by

$$g_1(X) := x - \tau \frac{\pi m^2 (\sigma_1^2 - \sigma_2^2)x}{\sigma_2 (\sigma_1 m + \alpha x)^2} (t^2 - y^2),$$

$$g_2(X) := p_m^\alpha(x)z + \eta_m^\alpha(x),$$

$$g_3(X) := k_m^\alpha(x)y,$$

$$g_4(X) := q_m^\alpha(x)h + \omega_m^\alpha(x),$$

$$g_5(X) := k_m^\alpha(x)t.$$

From the analysis of section (3.1.2) we observe that the sequence $\hat{\lambda}_{i,N}^\ell(\overline{R})$, $i = 1,2$ converges to $\hat{\lambda}_{i,N}^\infty$, $i = 1,2$, where

$$\hat{\lambda}_{i,N}^\infty = \frac{2R_2^{m+1}\overline{R}^{m-1}}{\sigma_1(R_2^{2m} + \overline{R}^{2m}) + \sigma_2(R_2^{2m} - \overline{R}^{2m})} \psi_i$$

with

$$\psi_1 = (\sigma_1 m + \alpha \overline{R}) \quad \text{and} \quad \psi_2 = (\sigma_1 m - \alpha \overline{R}).$$

Similarly, the sequence $\hat{\lambda}_{i,D}^\ell(\overline{R})$, $i = 1,2$ converges to $\hat{\lambda}_{i,D}^\infty$, $i = 1,2$, where

$$\hat{\lambda}_{i,D}^\infty = \frac{2\sigma_2 C_f R_2^m \overline{R}^{m-1}}{\sigma_1(R_2^{2m} - \overline{R}^{2m}) + \sigma_2(R_2^{2m} + \overline{R}^{2m})} \psi_i.$$

Moreover, we have

$$\left\{ \begin{aligned} \hat{\lambda}_{1,N}^\infty &= p_m^\alpha(\overline{R})\hat{\lambda}_{2,N}^\infty + \eta_m^\alpha(\overline{R}) \\ \hat{\lambda}_{2,N}^\infty &= k_m^\alpha(\overline{R})\hat{\lambda}_{1,N}^\infty \\ \hat{\lambda}_{1,D}^\infty &= q_m^\alpha(\overline{R})\hat{\lambda}_{2,D}^\infty + w_m^\alpha(\overline{R}) \\ \hat{\lambda}_{2,D}^\infty &= k_m^\alpha(\overline{R})\hat{\lambda}_{1,D}^\infty \end{aligned} \right.$$

Using the expression of $C_f$ given by (2.7) one can easily check that $\hat{\lambda}_{i,N}^\infty = \hat{\lambda}_{i,D}^\infty$ and then $\overline{X} = (\overline{R}, \hat{\lambda}_{1,N}^\infty, \hat{\lambda}_{2,N}^\infty, \hat{\lambda}_{1,D}^\infty, \hat{\lambda}_{2,D}^\infty) = (\overline{R}, \hat{\lambda}_{1,N}^\infty, \hat{\lambda}_{2,N}^\infty, \hat{\lambda}_{1,N}^\infty, \hat{\lambda}_{2,N}^\infty)$ is a fixed point of $G$. We first establish in the following lemma that $\overline{X}$ is the unique fixed point of $G$.

**Lemma 3.4.** The point $\overline{X}$ defined above is the unique fixed point of $G$.

*Proof.* Let $X = (x, y, z, t, h)$, $x > 0$, be a fixed point of $G$. Then we have from $g_1(X) = x$ that $t = y$ and the other equalities imply

$$\left\{ \begin{aligned} y &= p_m^\alpha(x)z + \eta_m^\alpha(x), \\ z &= k_m^\alpha(x)y, \\ t &= q_m^\alpha(x)h + w_m^\alpha(x) \\ h &= k_m^\alpha(x)t. \end{aligned} \right.$$

Consequently

$$\begin{pmatrix} y \\ z \end{pmatrix} = \begin{pmatrix} 0 & p_m^\alpha(x) \\ k_m^\alpha(x) & 0 \end{pmatrix} \begin{pmatrix} y \\ z \end{pmatrix} + \begin{pmatrix} \eta_m^\alpha(x) \\ 0 \end{pmatrix}$$
and

$$ \begin{pmatrix} t \\ h \end{pmatrix} = \begin{pmatrix} 0 & q_m^\alpha(x) \\ k_m^\alpha(x) & 0 \end{pmatrix} \begin{pmatrix} t \\ h \end{pmatrix} + \begin{pmatrix} w_m^\alpha(x) \\ 0 \end{pmatrix}. $$

Then we get

$$ \left\{ \begin{aligned} y &= \frac{2R_2^{m+1}x^{m-1}(\sigma_1 m + \alpha x)}{\sigma_1(R_2^{2m} + x^{2m}) + \sigma_2(R_2^{2m} - x^{2m})} \\ t &= \frac{2\sigma_2 C_f R_2^m x^{m-1}(\sigma_1 m + \alpha x)}{\sigma_1(R_2^{2m} - x^{2m}) + \sigma_2(R_2^{2m} + x^{2m})}. \end{aligned} \right. $$

The equality $y = t$ then implies

$$ C_f = \frac{R_2 [\sigma_1(R_2^{2m} - x^{2m}) + \sigma_2(R_2^{2m} + x^{2m})]}{\sigma_2 [\sigma_1(R_2^{2m} + x^{2m}) + \sigma_2(R_2^{2m} - x^{2m})]}. $$

Since the function $g$ defined on $]0, R_2[$ by

$$ g(x) = \frac{R_2 [\sigma_1(R_2^{2m} - x^{2m}) + \sigma_2(R_2^{2m} + x^{2m})]}{\sigma_2 [\sigma_1(R_2^{2m} + x^{2m}) + \sigma_2(R_2^{2m} - x^{2m})]} $$

is one to one, we deduce from the expression of $C_f$ given by (2.7) that $x = \bar{R}$. This finally gives $X = (\bar{R}, \hat{\lambda}_{1,N}^\infty, \hat{\lambda}_{2,N}^\infty, \hat{\lambda}_{1,D}^\infty, \hat{\lambda}_{2,D}^\infty) = \bar{X}$. $\square$

**Theorem 3.1.** There exists $\alpha_0 > 0$ such that for all $\alpha \in ]0, \alpha_0[$ Algorithm 2 with $L = 1$ is locally convergent for all $\tau \in ]0, \delta_\alpha[$ for some $\delta_\alpha > 0$.

*Proof.* For short notation, let us set $\bar{x} = \bar{R}, \bar{y} = \hat{\lambda}_{1,N}^\infty$, and

$$ k := k_m^\alpha(\bar{x}) = \frac{\sigma_1 m - \alpha \bar{x}}{\sigma_1 m + \alpha \bar{x}} = \frac{\psi_2}{\psi_1}. $$

Then from the previous Lemma, $\bar{X} = (\bar{x}, \bar{y}, k\bar{y}, \bar{y}, k\bar{y})$ is the unique fixed point of $G$. In order to prove the local convergence we shall establish that $G$ is a contraction in a neighborhood of $\bar{X}$. This requires the study of the Jacobian matrix $DG(\bar{X})$. The latter is given by

$$ DG(\bar{X}) = \begin{pmatrix} 1 & d\tau & 0 & -d\tau & 0 \\ a & 0 & p & 0 & 0 \\ b & k & 0 & 0 & 0 \\ c & 0 & 0 & 0 & q \\ b & 0 & 0 & k & 0 \end{pmatrix} $$
where we have set for short notation

$$d := \frac{1}{\tau} \frac{\partial g_1}{\partial y}(\bar{X}) = \frac{2\pi m^2 (\sigma_1^2 - \sigma_2^2) \bar{x}}{\sigma_2 (\sigma_1 m + \alpha \bar{x})^2} \bar{y}$$

$$a := \frac{\partial g_2}{\partial x}(\bar{X}) = \frac{-2m\alpha\sigma_2(4mR_2^{2m}\bar{x}^{2m} - \bar{x}^{4m} + R_2^{4m})}{[\sigma_2 m(R_2^{2m} - \bar{x}^{2m}) + \alpha\bar{x}(R_2^{2m} + \bar{x}^{2m})]^2} k\bar{y} \\ + \frac{4m\alpha R_2^{m+1}\bar{x}^{m-1}[(\sigma_2 m^2 - \alpha\bar{x})(R_2^{2m} + \bar{x}^{2m}) + \alpha\bar{x}m(R_2^{2m} - \bar{x}^{2m})]}{[\sigma_2 m(R_2^{2m} - \bar{x}^{2m}) + \alpha\bar{x}(R_2^{2m} + \bar{x}^{2m})]^2}$$

$$p := \frac{\partial g_2}{\partial z}(\bar{X}) = p_m^\alpha(\bar{x})$$

$$b := \frac{\partial g_3}{\partial x}(\bar{X}) = \frac{-2\sigma_1 m \alpha}{(\sigma_1 m + \alpha \bar{x})^2} \bar{y}$$

$$c := \frac{\partial g_4}{\partial x}(\bar{X}) = \frac{2m\alpha\sigma_2(4mR_2^{2m}\bar{x}^{2m} + \bar{x}^{4m} - R_2^{4m})}{[\sigma_2 m(R_2^{2m} + \bar{x}^{2m}) + \alpha\bar{x}(R_2^{2m} - \bar{x}^{2m})]^2} k\bar{y} \\ + \frac{4m\sigma_2\alpha C_f R_2^m \bar{x}^{m-1} [(\sigma_2 m^2 - \alpha\bar{x})(R_2^{2m} - \bar{x}^{2m}) + \alpha\bar{x}m(R_2^{2m} + \bar{x}^{2m})]}{[\sigma_2 m(R_2^{2m} + \bar{x}^{2m}) + \alpha\bar{x}(R_2^{2m} - \bar{x}^{2m})]^2},$$

$$q := \frac{\partial g_4}{\partial h}(\bar{X}) = q_m^\alpha(\bar{x}).$$

For every $(\lambda, \tau) \in \mathbb{R} \times \mathbb{R}$, we define the characteristic polynomial of $DG(\bar{X})$

$$P(\lambda, \tau) := \det(DG(\bar{X}) - \lambda I_5) = \begin{vmatrix} 1 - \lambda & d\tau & 0 & -d\tau & 0 \\ a & -\lambda & p & 0 & 0 \\ b & k & -\lambda & 0 & 0 \\ c & 0 & 0 & -\lambda & q \\ b & 0 & 0 & k & -\lambda \end{vmatrix}$$

and we denote by $\lambda_i(\tau) \in \mathbb{C}, i = 1, \dots, 5$, the eigenvalues of $DG(\bar{X})$. Then, we have

$$P(\lambda, \tau) = (1 - \lambda)(\lambda^2 - kp)(\lambda^2 - kq) - \tau d [(c - a)\lambda^3 + b(q - p)\lambda^2 + (aq - cp)k\lambda],$$

and $P(\lambda_i(\tau), \tau) = 0$ for $i = 1, \dots, 5$. Let $\bar{\rho}(\tau) := \max\{|\lambda_i(\tau)|; i = 1, \dots, 5\}$ be the spectral radius of $DG(\bar{X})$. We remark that for $\tau = 0$, we have

$$P(\lambda, 0) = (1 - \lambda)(\lambda^2 - kp)(\lambda^2 - kq),$$

and therefore, the eigenvalues of the $DG(\bar{X})$ for $\tau = 0$ are given by

$$(\lambda_1(0), \lambda_2(0), \lambda_3(0), \lambda_4(0), \lambda_5(0)) = (1, \sqrt{kp}, -\sqrt{kp}, \sqrt{kq}, -\sqrt{kq})$$

where $|kp| < 1$, $|kq| < 1$, and $kq \neq kp$ if $\alpha$ sufficiently small. Then, for $\tau = 0$ and $\alpha$ sufficiently small, all eigenvalues of $DG(\bar{X})$ are simple and verify: $|\lambda_1(0)| = 1$, and $|\lambda_i(0)| < 1$ for $i \in \{2, 3, 4, 5\}$. Consequently, there exists a small enough number $\delta > 0$ such that all eigenvalues of $DG(\bar{X})$ are simple for all $\tau \in [0, \delta]$ and

$$\lim_{\tau \to 0} \lambda_i(\tau) = \lambda_i(0) \quad \forall i \in \{1, 2, 3, 4, 5\}.$$
Thus one can choose $\delta > 0$ sufficiently small such that $|\lambda_i(\tau)| < 1$ for $i \in \{2, 3, 4, 5\}$, $\bar{\rho}(\tau) = |\lambda_0(\tau)|$ for all $\tau \in ]0, \delta[$. Moreover $\lim_{\tau \to 0} \bar{\rho}(\tau) = 1$. To show that if $\bar{\rho}(\tau) < 1$ for small enough $\tau > 0$, we use a first order Taylor expansion of $(\lambda, \tau) \mapsto P(\lambda, \tau)$ in a neighborhood of the point $(1, 0)$. Indeed

$$ 
\begin{aligned}
P(\lambda, \tau) &= P(1, 0) + \frac{\partial P}{\partial \lambda}(1, 0) (\lambda - 1) + \frac{\partial P}{\partial \tau}(1, 0) \tau + \epsilon(\lambda, \tau) \\
&= C_1(1 - \lambda) - C_2\tau + \epsilon(\lambda, \tau)
\end{aligned}
$$

with $\lim_{(\lambda, \tau) \to (1, 0)} \frac{|\epsilon(\lambda, \tau)|}{\sqrt{(1 - \lambda)^2 + \tau^2}} = 0$, and

$$ 
\left\{ 
\begin{aligned}
C_1 &:= (1 - kp)(1 - kq), \\
C_2 &:= d \left[ (c - a)(1 - kp) + (q - p)(b + ak) \right].
\end{aligned} 
\right. 
\qquad (3.28) 
$$

This shows that for $\tau$ sufficiently small

$$ \lambda_0(\tau) = 1 - \frac{C_2}{C_1}\tau + o(\tau). \qquad (3.29) $$

Clearly $C_1 > 0$ from the conditions $|kp| < 1$ and $|kq| < 1$. Let $\alpha \in ]0, \alpha_0[$ where $\alpha_0$ is as in Lemma 3.5. Then we also have that $C_2 > 0$. We then infer from equation (3.29) the existence of $\delta_\alpha > 0$ such that $\bar{\rho}(\tau) = |\lambda_0(\tau)| < 1$ for all $\tau \in ]0, \delta_\alpha[$. This proves that $G$ is a contraction for $\tau \in ]0, \delta_\alpha[$, which gives the desired result. $\square$

**Lemma 3.5.** There exists $\alpha_0 > 0$ such that $C_2 > 0$ for all $0 < \alpha < \alpha_0$ where $C_2 = C_2(\alpha)$ is given by (3.28).

*Proof.* We first observe that

$$ \bar{y} \sim k\bar{y} \sim \frac{2m\sigma_1 R_2^{m+1} \bar{R}^{m-1}}{\sigma_1(R_2^{2m} + \bar{R}^{2m}) + \sigma_2(R_2^{2m} - \bar{R}^{2m})} \quad \text{as } \alpha \to 0. $$

Then, using a first order Taylor expansion of $a, b, c, d, k, p, q$, in a neighborhood of $\alpha = 0$ we get the following equivalences as $\alpha \to 0$

$$ d \sim \frac{2\pi(\sigma_1^2 - \sigma_2^2)\bar{x}}{\sigma_1^2\sigma_2}\bar{y} \sim \frac{4m\pi(\sigma_1^2 - \sigma_2^2)R_2^{m+1}\bar{R}^m}{\sigma_1\sigma_2[\sigma_1(R_2^{2m} + \bar{R}^{2m}) + \sigma_2(R_2^{2m} - \bar{R}^{2m})]} \qquad (3.30) $$
$$b \sim \frac{-2\alpha}{m\sigma_1}\bar{y} \sim \frac{-4\alpha R_2^{m+1}\bar{R}^{m-1}}{\sigma_1(R_2^{2m} + \bar{R}^{2m}) + \sigma_2(R_2^{2m} - \bar{R}^{2m})}$$

$$k \sim 1 - \frac{2\alpha\bar{R}}{m\sigma_1}$$

$$p \sim 1 - \frac{2\alpha\bar{R}(R_2^{2m} + \bar{R}^{2m})}{m\sigma_2(R_2^{2m} - \bar{R}^{2m})}$$

$$q \sim 1 - \frac{2\alpha\bar{R}(R_2^{2m} - \bar{R}^{2m})}{m\sigma_2(R_2^{2m} + \bar{R}^{2m})}$$

$$a \sim \frac{-2\alpha(4mR_2^{2m}\bar{R}^{2m} - \bar{R}^{4m} + \bar{R}_2^{4m})}{m\sigma_2(R_2^{2m} - \bar{R}^{2m})^2}\bar{y} + \frac{4m\alpha R_2^{m+1}\bar{R}^{m-1}(R_2^{2m} + \bar{R}^{2m})}{\sigma_2(R_2^{2m} - \bar{R}^{2m})^2}$$

$$\sim \frac{4\alpha R_2^{m+1}\bar{R}^{m-1}[m\sigma_1(R_2^{2m} - \bar{R}^{2m}) + (m\sigma_2 - \sigma_1)(R_2^{2m} + \bar{R}^{2m})]}{\sigma_2(R_2^{2m} - \bar{R}^{2m})[\sigma_1(R_2^{2m} + \bar{R}^{2m}) + \sigma_2(R_2^{2m} - \bar{R}^{2m})]}$$

$$c \sim \frac{2\alpha(4mR_2^{2m}\bar{R}^{2m} + \bar{R}^{4m} - \bar{R}_2^{4m})}{m\sigma_2(R_2^{2m} + \bar{R}^{2m})^2}\bar{y} + \frac{4m\alpha C_f R_2^m \bar{R}^{m-1}(R_2^{2m} - \bar{R}^{2m})}{(R_2^{2m} + \bar{R}^{2m})^2}$$

$$\sim \frac{4\alpha R_2^{m+1}\bar{R}^{m-1}[m\sigma_1(R_2^{2m} + \bar{R}^{2m}) + (m\sigma_2 - \sigma_1)(R_2^{2m} - \bar{R}^{2m})]}{\sigma_2(R_2^{2m} + \bar{R}^{2m})[\sigma_1(R_2^{2m} + \bar{R}^{2m}) + \sigma_2(R_2^{2m} - \bar{R}^{2m})]}$$

Then we infer that

$$c - a \sim \frac{16\alpha(\sigma_1 - m\sigma_2)R_2^{3m+1}\bar{R}^{3m-1}}{\sigma_2(R_2^{4m} - \bar{R}^{4m})[\sigma_1(R_2^{2m} + \bar{R}^{2m}) + \sigma_2(R_2^{2m} - \bar{R}^{2m})]}$$
$$1 - kp \sim \frac{2\alpha\bar{R}[\sigma_1(R_2^{2m} + \bar{R}^{2m}) + \sigma_2(R_2^{2m} - \bar{R}^{2m})]}{m\sigma_1\sigma_2(R_2^{2m} - \bar{R}^{2m})} \quad (3.31)$$

$$1 - kq \sim \frac{2\alpha\bar{R}[\sigma_1(R_2^{2m} - \bar{R}^{2m}) + \sigma_2(R_2^{2m} + \bar{R}^{2m})]}{m\sigma_1\sigma_2(R_2^{2m} + \bar{R}^{2m})} \quad (3.32)$$

$$q - p \sim \frac{8\alpha R_2^{2m}\bar{R}^{2m+1}}{m\sigma_2(R_2^{4m} - \bar{R}^{4m})}$$

$$b + ak \sim b + a \sim \frac{4\alpha R_2^{m+1}\bar{R}^{m-1}[m\sigma_1 - \sigma_2](R_2^{2m} - \bar{R}^{2m}) + (m\sigma_2 - \sigma_1)(R_2^{2m} + \bar{R}^{2m})}{\sigma_2(R_2^{2m} - \bar{R}^{2m})[\sigma_1(R_2^{2m} + \bar{R}^{2m}) + \sigma_2(R_2^{2m} - \bar{R}^{2m})]}$$
Consequently,

$$ (c-a)(1-kp) \sim \frac{32\alpha^2(\sigma_1 - m\sigma_2)R_2^{3m+1}\bar{R}^{3m}}{m\sigma_1\sigma_2(R_2^{4m} - \bar{R}^{4m})(R_2^{2m} - \bar{R}^{2m})} $$

$$ (q-p)(b+ak) \sim \frac{32\alpha^2 R_2^{3m+1} \bar{R}^{3m} \left[ (m\sigma_1 - \sigma_2)(R_2^{2m} - \bar{R}^{2m}) + (m\sigma_2 - \sigma_1)(R_2^{2m} + \bar{R}^{2m}) \right]}{m\sigma_2^2 (R_2^{4m} - \bar{R}^{4m})(R_2^{2m} - \bar{R}^{2m}) \left[ \sigma_1(R_2^{2m} + \bar{R}^{2m}) + \sigma_2(R_2^{2m} - \bar{R}^{2m}) \right]}, $$

and from equations (3.28) and (3.30), we obtain

$$ C_2 \sim \frac{128\pi\alpha^2 m (\sigma_1^2 - \sigma_2^2)^2 R_2^{4m+2} \bar{R}^{4m}}{\sigma_1^2 \sigma_2^3 (R_2^{4m} - \bar{R}^{4m}) \left[ \sigma_1 (R_2^{2m} + \bar{R}^{2m}) + \sigma_2 (R_2^{2m} - \bar{R}^{2m}) \right]^2} > 0. $$

□

**Remark 3.2.** Using equations (3.28), (3.31) and (3.32), we obtain

$$ C_1 \sim \frac{4\alpha^2 \bar{R}^2 \left[ \sigma_1 (R_2^{2m} + \bar{R}^{2m}) + \sigma_2 (R_2^{2m} - \bar{R}^{2m}) \right] \left[ \sigma_1 (R_2^{2m} - \bar{R}^{2m}) + \sigma_2 (R_2^{2m} + \bar{R}^{2m}) \right]}{m^2 \sigma_1^2 \sigma_2^2 (R_2^{4m} - \bar{R}^{4m})} $$

as $\alpha \to 0$. Therefore,

$$ \frac{C_2}{C_1} \sim \mathcal{J}''(\bar{R}) \text{ as } \alpha \to 0. \qquad (3.33) $$

This shows that the local convergence of Algorithm 2 is equivalent to the local convergence of Algorithm 1 in this specific case. In fact, some calculations with Maple for $m = 1$ and $m = 2$ show that we have exact equality $C_2/C_1 = \mathcal{J}''(\bar{R})$ for all $\alpha > 0$ and $\bar{R} > 0$. We conjecture that this also holds true for all $m > 0$.

# 4 Numerical experiments and validation

The goal of this section is to test the efficiency of Algorithm 2 in comparison with the standard descent gradient described in Algorithm 1. We shall employ synthetic data numerically generated using a finite elements solver designed with the help of FreeFem++ [24]. Indeed, we use a direct solver to generate the data, while the OSM is used in the inversion for Algorithm 2. We further avoid any inverse crime by making sure that the meshes used for generating the data have no connections with the ones used in the inversion. Actually, the latter vary during iterations since the interface $\Sigma$ changes. Moreover, in most of the examples below the geometry of $\Sigma$ can not be exactly represented by the parametrization (2.3) used in the inversion algorithms.

For all the examples below, the domain $\Omega$ is the open disk of center $(0,0)$ and radius $R_2 = 2$, the current flux $\phi(\theta) = \cos(\theta)$, $\theta \in [0, 2\pi]$, the conductivity parameters $\sigma_1 = 1$, $\sigma_2 = 2$ and the OSM parameter $\alpha = 1$. The measured data $f$ is represented by the values $f_i$, $i = 1, \dots, N_\Gamma$ of the numerical solution $u_{\text{num}}$ at the nodes belonging to $\Gamma$. In order to simulate noise in the data $f$ we artificially corrupt the computed values $f_i$ with random noise as follows

$$ f_i^\epsilon := f_i + \epsilon(1 - 2r_i)f_i, \quad i = 1, \dots, N_\Gamma, $$
where $r_i$ are randomly chosen between 0 and 1 and $\epsilon$ denotes the noise level. In addition to representing the obtained geometrical reconstructions we shall also give the evolution of the dimensionless square-root cost functional $\sqrt{\mathcal{J}/\mathcal{J}_0}$ during iterations where $\mathcal{J}$ is given by (2.4) and

$$ \mathcal{J}_0 := \int_{\Omega} \sigma |\nabla u_{\text{num}}|^2 dx. $$

## 4.1 Discussion of Algorithm 2 for a kite shape

In the first example we choose $\Sigma$ to be a kite defined by

$$ x(t) = \cos(t) + 0.5 \cos(2t) - 0.4 \text{ and } y(t) = 1.2 \sin(t), \ t \in [0, 2\pi]. \quad (4.1) $$

For the inversion algorithms we use the parametrization (2.3) with $n = 19$. The initial guess is $R_j^0 = 1.8$, for $j = 0, \dots, n-1$. The results of the inversion are given in Figure 3(b) for Algorithm 1 and Figure 3(c) for Algorithm 2 with $L = 1$, i.e only one OSM step is used at each gradient descent iteration. We observe that we qualitatively obtain the same accuracy for both algorithms. The evolution of the cost functional is depicted in Figure 4(left). We remark that the cost functional increases in the first iterations for Algorithm 2, which means that the approximated OSM solution is not yet sufficiently close to the exact one and therefore the approximate gradient is not yet a descent direction. This is corrected as the iteration number increases. Whence the iteration number is sufficiently large we notice that the speed of convergence represented by the slope of the curves is roughly the same between the two algorithms. This indeed shows the superiority and relevance of Algorithm 2 which achieves comparative performances with potentially much cheaper numerical cost per iteration. Figure 4(middle) and (right) show the effect of increasing the number $L$ of OSM steps. We naturally observe that as this number increases, Algorithm 2 becomes closer to Algorithm 1. Let us also mention that, other numerical tests not reported here suggest that this observation also holds when we fix $L$, for instance $L = 1$, and decrease the descent step $\tau$ (see also Figure 9(right)).

<figure><img src="image_6.png" /><figcaption>Figure 3: Comparison between Algorithm 1 and Algorithm 2 with $L = 1$ for the case of the kite parameterized by (4.1) and for noise free data. The exact shape and initial guess are shown in Figure (a). The gradient descent parameter $\tau = 0.05$ for both algorithms.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 4: Comparison of the evolution of $log_{10}(\sqrt{\mathcal{J}/\mathcal{J}_0})$ between Algorithm 1 and Algorithm 2 with $L = 1$ (left), $L = 2$ (middle) and $L = 5$ (right) for the example shown in Figure 3. The gradient descent parameter is $\tau = 0.05$.</figcaption></figure>

Figure 5 shows the reconstructions obtained by Algorithm 2 for the example discussed in Figure 3 but for noisy data with noise level $\epsilon = 1\%$ (left), $\epsilon = 3\%$ (middle) and $\epsilon = 5\%$ (right). We observe robustness of the obtained results and a good accuracy which is very similar to the one obtained by Algorithm 1 (not reported here).

<figure><img src="image_8.png" /><figcaption>Figure 5: Reconstructions obtained by Algorithm 2 with $L = 1$ for the example discussed in Figure 3 but for noisy data with noise level $\epsilon = 1\%$ (left), $\epsilon = 3\%$ (middle) and $\epsilon = 5\%$ (right).</figcaption></figure>

## 4.2 The example of other geometries

We report in Figure 6 and Figure 7 the reconstructions obtained by Algorithm 2, $L = 1$ for other geometries, keeping the same parameters as in Example 4.1 and with noise levels $\epsilon = 3\%$ and $\epsilon = 5\%$. The parametrizations of these geometries are given by Table 1.

<table><caption>Table 1: The boundary parametrization of the geometries.</caption><thead><tr><th>Geometry type</th><th>Parametrization $(x(t), y(t))$</th></tr></thead><tbody><tr><td>Peanut shaped</td><td>1.5√(cos<sup>2</sup>(t) + 0.23 sin<sup>2</sup>(t))(cos(t), sin(t))</td></tr><tr><td>Pear</td><td>5.5 + sin(3t)/5 (cos(t), sin(t))</td></tr><tr><td>Star</td><td>(1.2 + 0.4 sin(5t))(cos(t), sin(t))</td></tr></tbody></table>
We clearly observe good performances in terms of robustness and accuracy.

<figure><img src="image_9.png" /><figcaption>Figure 6: Reconstructions obtained by Algorithm 2 for the geometries discussed in Table 1 but for noisy data with noise level $\epsilon = 3\%$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7: Reconstructions obtained by Algorithm 2 for the geometries discussed in Table 1 but for noisy data with noise level $\epsilon = 5\%$.</figcaption></figure>

## 4.3 Discussion on the choice of $L$

As mentioned earlier for a given $\tau$, it is not clear which value of $L$ would lead to potentially best performances of Algorithm 2. We provide here some elementary numerical investigations in the case where the exact geometry is given by the parametrization $\mathcal{T}_n$ (2.3) for given $n$. In that case it is possible to use a stopping rule related to the accuracy of the reconstruction, $e^k$ defined as $e^k = \frac{\|\bar{R}-R^k\|_2}{\|\bar{R}\|_2}$ where $\Sigma_{\bar{R}}$ is the exact interface.

In the example below we choose $n = 9$, $\bar{R} = (1, 0.8, 0.7, 0.9, 1, 1.7, 1.6, 1.5, 1.4)$ and the initial guess $R^0 = (1.7, 1.6, 1.55, 1.4, 1.3, 1.3, 1.2, 1.4, 1.5)$ (see Figure 8(a)).
<figure><img src="image_11.png" /><figcaption>Figure 8: Comparison between Algorithm 1 and Algorithm 2 with $L = 1$ to achieve 5% relative accuracy in the case of the geometry given by Figure (a) and for noise free data. The gradient descent parameter $\tau = 0.1$ for both algorithms.</figcaption></figure>

Table 2 gives the number of iterations $K$ need to achieve an accuracy $e^k \simeq 5\%$ for noise free data. We also indicate the virtually total cost $2KL$ that represents the total number of iterations for direct problems. We observe that $L = 1$ provides the lowest values for $KL$ and $L = 2$ or $L = 3$ comparative performances. This shows again the relevance of Algorithm 2 and that a few DDM iterations are sufficient to provide good solution for the inverse problem.

<table><caption>Table 2: Total number of iterations for Algorithm 2 to achieve 5% relative accuracy in the case of the geometry given by Figure 8 and for noise free data.</caption><thead><tr><th rowspan="2">L</th><th colspan="2">τ = 0.1</th><th colspan="2">τ = 0.05</th><th colspan="2">τ = 0.005</th></tr><tr><th>K</th><th>2KL</th><th>K</th><th>2KL</th><th>K</th><th>2KL</th></tr></thead><tbody><tr><td>1</td><td>53</td><td>106</td><td>55</td><td>110</td><td>359</td><td>718</td></tr><tr><td>2</td><td>29</td><td>116</td><td>41</td><td>164</td><td>355</td><td>1420</td></tr><tr><td>3</td><td>29</td><td>174</td><td>39</td><td>234</td><td>355</td><td>2130</td></tr><tr><td>4</td><td>30</td><td>240</td><td>41</td><td>328</td><td>354</td><td>2832</td></tr><tr><td>5</td><td>30</td><td>300</td><td>39</td><td>390</td><td>354</td><td>3540</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>Figure 9: Comparison of the evolution of $\log_{10}(\sqrt{\mathcal{J}/\mathcal{J}_0})$ between Algorithm 1 and Algorithm 2 applied to the example of Figure 8 with different gradient descent parameter $\tau = 0.1$ (left), $\tau = 0.05$ (middle) and $\tau = 0.005$ (right).</figcaption></figure>
References

[1] L. Afraites, M. Dambrine, and D. Kateb. *Shape Methods for the Transmission Problem with a Single Measurement*. Numer. Funct. Anal. Optim., 28(5-6), 519-551, 2007.

[2] L. Afraites, M. Dambrine, and D. Kateb. *On second order shape optimization methods for electrical impedance tomography*. SIAM J. Control Optim., 47(3), 1556-1590, 2008.

[3] A. Aspri, E. Beretta, E. Francini and S. Vessella, *Lipschitz Stable Determination of Polyhedral Conductivity Inclusions from Local Boundary Measurements*, SIAM Journal on Mathematical Analysis, 54, 5, 5182-5222, 2022.

[4] K. Astala and L. Päivärinta. *Calderón's inverse conductivity problem in the plane*. Ann. of Math., 163, 265-299, 2006.

[5] E. Beretta and E. Francini. *Global Lipschitz stability estimates for polygonal conductivity inclusions from boundary measurements*. Appl. Anal., 01(10), 35363549, 2022.

[6] M. Bonazzoli, H. Haddar, and T. A. Vu. *Convergence analysis of multi-step one-shot methods for linear inverse problems*. [Research Report] RR-9477, Inria Saclay; ENSTA ParisTech, 2022.

[7] L. BORCEA. *Electrical impedance tomography*. Inverse Problems, 18(6), R99R136, 2002.

[8] M. Burger and W. Mhlhuber. *Iterative regularization of parameter identification problems by sequential quadratic programming methods*. Inverse Problems, 18(4), 943969, 2002.

[9] S. Chaabane, B. Charfi, and H. Haddar. *Reconstruction of discontinuous parameters in a second order impedance boundary operator*. Inverse Problems, 32(10), 105004, 2016.

[10] S. Chaabane, C. Elhechmi, and M. Jaoua. *A stable recovery algorithm for the Robin inverse problem*. IMACS, Math. Comput. Simul., 66(4-5), 367-383, 2004.

[11] S. Chaabane, I. Feki, and N. Mars. *Numerical reconstruction of a piecewise constant Robin parameter in the two-or three-dimensional case*. Inverse Problems, 28(6), 065016, 2012.

[12] S. Chaabane and M. Jaoua. *Un algorithme d'identification de frontières soumises à des conditions aux limites de Signorini*. M2AN, 34(3), 707-722, 2000.

[13] V. Dolean, P. Jolivet, and F. Nataf. *An Introduction to Domain Decomposition Methods*. SIAM, Philadelphia, 2015.

[14] O. Dubois. *Optimized Schwarz methods for the advection-diffusion equation and for problems with discontinuous coefficients*. PhD thesis. McGill University, 2007.

[15] M. J. Gander. *Optimized Schwarz methods*. SIAM J. Numer. Anal., 44(2), 699-731, 2006.

[16] N. Gauger, A. Griewank, A. Hamdi, C. Kratzenstein, E. zkaya, and T. Slawig. *Automated extension of fixed point PDE solvers for optimal design with bounded retardation*. In: Constrained Optimization and Optimal Control for Partial Differential Equations. International Series of Numerical Mathematics, 160, 99122, Springer Basel, Basel, 2012.
[17] L. Gerardo-Giorda and F. Nataf. *Optimized Schwarz methods for unsymmetric layered problems with strongly discontinuous and anisotropic coefficients*. Technical report, 561, submitted. Ecole Polytechnique, CMAP, CNRS UMR 7641, Paris, France, 2004.

[18] A. Griewank. *Projected Hessians for Preconditioning in One-Step One-Shot Design Optimization*. In: Large-Scale Nonlinear Optimization. Nonconvex Optimization and Its Applications, 83, 151171, Springer US, Boston, MA, 2006.

[19] E. Haber and U. M. Ascher. *Preconditioned all-at-once methods for large, sparse parameter estimation problems*. Inverse Problems, 17(6), 18471864, 2001.

[20] H. Haddar and R. Kress. *Conformal mapping and impedance tomography*. Inverse Problems, 26(7), 074002, 2010.

[21] A. Hamdi and A. Griewank. *Reduced quasi-Newton method for simultaneous design and optimization*. Computational Optimization and Applications, 49(3), 521548, 2009.

[22] M. Hanke, B. Harrach, and N. Hyvnen. *Justification of point electrode models in electrical impedance tomography*. Math. Models Methods Appl. Sci., 21(6), 1395-1413, 2011.

[23] M. Hassine and I. Kallel. *Kohn-Vogelius formulation and topological sensitivity analysis based method for solving geometric inverse problems*. Arab. J. Math. Sci., 24(1), 4362, 2018.

[24] F. Hecht. *New development in FreeFem++*. J. Numer. Math., 20(3-4), 251265, 2012.

[25] F. Hettlich and W. Rundell. *The determination of a discontinuity in a conductivity from a single boundary measurement*. Inverse Problems, 14(1), 67-82, 1998.

[26] V. Isakov. *Inverse Problems for Partial Differential Equations*. Applied Mathematical Sciences, 127, Springer, 2017.

[27] C. Japhet and F. Nataf. *The best interface conditions for domain decomposition methods: Absorbing boundary conditions*. In: Absorbing Boundaries and Layers, Domain Decomposition Methods: Applications to Large Scale Computations, 348-373, Nova Science, New York, 2001.

[28] C. Japhet, F. Nataf, and F. Rogier. *The optimized order 2 method. Application to convection-diffusion problems*. Future Gener. Comput. Sys., 18(1), 17-30, 2001.

[29] B. Kaltenbacher, A. Kirchner, and B. Vexler. *Goal oriented adaptivity in the IRGNM for parameter identification in PDEs II: all-at-once formulations*. Inverse Problems, 30(4), 045002, 2014.

[30] T. A. Khan. *Review on Electrical Impedance Tomography: Artificial Intelligence Methods and its Applications*. Algorithms, 12(5), 88, 2019.

[31] R. V. Kohn and M. Vogelius. *Relaxation of a variational method for impedance computed tomography*. Communications on Pure and Applied Mathematics, 40(6), 745-777, 1987.

[32] P.-L. Lions. *On the Schwarz alternating method. III. A variant for nonoverlapping subdomains*. In: Third International Symposium on Domain Decomposition Methods for Partial Differential Equations. SIAM, 202-223, 1990.
[33] J. L. Mueller and S. Siltanen. *Linear and Nonlinear Inverse Problems with Practical Applications*. SIAM, 2012.