1 Introduction↩︎

1.1 Background of inverse interface problems and direct sampling methods↩︎

Interface problems are ubiquitous in various scientific and engineering disciplines. For instance, when a physical domain is composed of two different materials or media, an interface naturally emerges between them, leading to an interface problem. In the context of Partial Differential Equations (PDEs), numerical solutions to interface problems often require partitioning the computational domain into distinct regions or subdomains. These regions may exhibit different material properties or physical conditions and are typically separated by interfaces or boundaries. For example, in electrical conduction problems, resistivity varies across different materials, resulting in irregularities such as discontinuous parameters, jump conditions, and singular source terms in the governing differential equations at the interfaces or boundaries. Such complexities pose significant challenges in solving interface problems.

A flux jump condition describes a discontinuous change in the flux (or flow) of a quantity across an interface. This phenomenon typically arises when there is a sharp transition or inhomogeneity in the physical properties at the interface, and is particularly relevant in problems involving composite materials, phase boundaries, or fractured media.

In this paper, we focus on inverse interface problems in elliptic equations. While addressing two physically distinct scenarios, we unify them mathematically as the problem of recovering the spatial support of a singular perturbation. The first part of our work addresses the identification of the flux jump location at a known interior interface with a flux boundary condition, where the interface represents a surface with cracks, defects, or damages. In the second part, we identify the location of an entirely unknown interface. Both problems are formulated and solved using limited Cauchy measurements on the outer boundary. These problems find applications in diverse fields, including electrical conductivity problems [1], multi-material compressible flow simulations [2], medical diffusion MRI imaging [3], and geological surface reconstruction [4].

Interface problems have been extensively studied through both analytical and computational approaches. For forward problems involving discontinuous coefficients at fixed or moving interfaces, the immersed interface methods (IIM) [5], [6] and certain optimization techniques [7] have emerged as accurate and efficient approaches for problems defined on irregular domains. In inverse interface problems, classical approaches often rely on iterative optimization techniques such as the Newton-type methods [8] and the level-set method [9]. However, these methods face challenges in handling noisy data and require relatively accurate initial guesses, which are crucial for their convergence. Apart from these iterative approaches, the class of sampling methods is also powerful for various inverse problems, including linear sampling methods [10], [11] and orthogonality sampling methods [12], [13]. These sampling methods define an index function on the domain of interest to indicate the inhomogeneous part, yielding satisfactory results. As direct methods, they offer significant computational efficiency and demonstrate promising potential for inverse problems.

We employ the Direct Sampling Method (DSM) to solve these inverse problems. Specifically, we compute the probing index \(I(x)\) of the interface, defined through a weighted inner product between probing functions and boundary measurement data (given as Dirichlet-type measurements). The DSM is highly efficient and provides accurate solutions to our inverse problems. While the DSM has been applied to inverse medium problems [14]–[22], our work presents a novel application for determining the interior interface condition. We also evaluate the index function \(I(x), x\in\Gamma\) by the adjoint method, which results in a PDE-based, highly direct approach.

Another main contribution is the first quantitative proof, to the best of our knowledge, of the almost orthogonal property for a DSM kernel arising from an inverse interface problem. Earlier DSM papers use almost orthogonality to explain why a sampling index is large near the target and small away from it; related notions, such as mutually almost orthogonal probing functions, have also been developed [19]. These arguments do not cover the kernel used here, because both Cauchy data are built into the reference problem through a Robin condition and the probing functions are rescaled in the Sobolev inner product. In the disk case, we derive the Fourier representation of the normalized kernel and prove explicit radial monotonicity estimates: the kernel is maximized at the probing point and decreases quantitatively as the radial variable moves away. This turns the almost orthogonal property from a numerical observation into a proved localization result for the proposed interface DSM.

In this work, we develop a direct sampling method to solve a two-dimensional inverse interface problem where the flux jump across the interface is proportional to the flow. Specifically, we consider the following interface Laplace problem in a domain \(\Omega\subset\mathbb{R}^{2}\): \[\left\{ \begin{align} -\Delta u=0 \quad \text{in} \quad &\Omega\setminus\Gamma, \\ \left[\frac{\partial u}{\partial \nu}\right]=\gamma u\quad \text{on} \quad &\Gamma,\\ [u] = 0 \quad \text{on} \quad &\Gamma,\\ \end{align} \right. \label{con:model95problem}\tag{1}\] where \(\Gamma\subset\Omega\) is a one-dimensional interface. The jump condition is defined as \[\left[ \frac{\partial u(\hat{x})}{\partial \nu} \right] := \left( \lim_{x \to \hat{x}, x \in \Omega_1} \nabla u(x) - \lim_{x \to \hat{x}, x \in \Omega_2} \nabla u(x) \right) \cdot \nu \text{ for } \hat{x} \in \Gamma,\] where \(\nu\) denotes the unit normal vector on \(\Gamma\). A pair of Cauchy data \((f, g)\) is measured on the outer boundary \(\partial\Omega\) as follows: \[f = u\bigg|_{L}\quad \text{and}\quad g = \frac{\partial u}{\partial n}\bigg|_{\partial\Omega},\label{Cauchy95def}\tag{2}\] where \(n\) represents the outward unit normal vector on \(\partial\Omega\), and \(L\subset\partial\Omega\) is the portion of the outer boundary where the measurement data \(f\) is available.

Our work addresses two distinct inverse problems. In Problem I, we consider a known interface \(\Gamma\) where the jump coefficient \(\gamma(x)\) consists of a known background value \(\gamma_{0}\) and a sparse perturbation. Mathematically, this sparsity is characterized by \(D = \text{supp}(\gamma-\gamma_{0})\) being a small subset of \(\Gamma\). The objective is to construct a probing index whose large values indicate the defect support \(D=\operatorname{supp}(\gamma-\gamma_0)\subset\Gamma\) from a single pair of Cauchy data \((f,g)\). In Problem II, the interface \(\Gamma\) and the jump coefficient \(\gamma\) are both unknown, and we aim to reconstruct a probing index that highlights the location of the unknown interface \(\Gamma\) from a single pair of Cauchy data \((f,g)\). Since the response of a defect may depend on the chosen boundary excitation, a single pair of data does not necessarily illuminate all connected components of \(D\). When several Cauchy pairs are available, the corresponding single pair probing indices can be combined to obtain a more complete reconstruction. The analysis and results for the general case that \([u]\neq 0\) will be addressed in future work.

The remainder of the paper is structured as follows. In Section 1.2, we state a generic well-posedness result for the forward interface problem, with the proof deferred to the Appendix. Section 2 introduces the general philosophy of the DSM and details the construction of probing functions for both problems. In this section, we propose an innovative approach by utilizing the two Cauchy data together, reformulated as a Robin boundary condition, to solve the forward reference problem \(u_0\), while retaining the Dirichlet boundary condition for measurement data fitting. It significantly enhances the robustness of the results in the presence of noise. Section 3 derives the Fourier representation of the DSM kernel and states its quantitative almost orthogonal property, whose rigorous radial monotonicity proof is provided in the Appendix. Furthermore, it employs a frequency-domain analysis to establish a systematic criterion for selecting the Robin parameter \(\alpha\). In Section 4, we present a novel refinement approach to sharpen the performance of the DSM. Section 5 showcases our numerical results. Finally, Section 6 summarizes the main findings and discusses potential future research directions.

1.2 Well-posedness of the interface problem↩︎

In this section, we state the well-posedness result for the forward problem that will be used throughout the paper. It is classical that the problem is uniquely solvable if \(\gamma<0\) almost everywhere on \(\Gamma\), due to the coercivity of the associated bilinear form. The following lemma records a generic solvability statement for the more delicate case in which \(\gamma\) may take positive values.

The proof, based on compactness and the analytic Fredholm theorem, is deferred to Appendix 7.1.

2 Direct sampling methods for inverse interface problems↩︎

This section presents the Direct Sampling Method (DSM) for solving the two inverse interface problems introduced in Section 1. Before introducing the DSM construction, we establish a unified mathematical framework that characterizes the essence of their inverse identification. The key observation is that both problems can be viewed as the detection of a singular perturbation relative to a properly chosen reference or background model. This viewpoint separates the known part of the physics, encoded in a background operator, from the unknown object to be recovered, encoded in the support of a perturbation term. More precisely, we denote the corresponding background operator by \(\mathcal{L}_0\) and write, at a formal distributional level, \[(\mathcal{L}_0 + \mathcal{S})u = 0 \quad \text{in } \mathcal{D}'(\Omega),\] where \(\mathcal{S}\) represents the unknown singular perturbation relative to the chosen background operator and is supported on a subset \(\Sigma \subset \Omega\). The primary objective is to determine the spatial support \(\Sigma = \text{supp}(\mathcal{S})\) using a limited set of Cauchy data pairs \((f, g)\). In this way, the known-interface coefficient problem and the unknown-interface geometry problem are placed in the same form: the DSM only has to identify where the discrepancy between the true model and the reference model is located. Within this unified framework, we consider the following two cases.

1. Problem I (Local Anomaly Detection): In the interface model 1 , the interface \(\Gamma\subset\Omega\) is known a priori and may be non-closed, while the Robin jump coefficient \(\gamma\in L^\infty(\Gamma)\) is unknown. We assume that \(\gamma\) has a sparse deviation from a known background value \(\gamma_0\). The background operator \(\mathcal{L}_0\) includes the known background jump condition \([\partial_\nu u]=\gamma_0u\) on \(\Gamma\), and the perturbation operator \(\mathcal{S}\) represents only the deviation of \(\gamma\) from \(\gamma_0\), acting distributionally as \(\mathcal{S}u = -(\gamma - \gamma_0)u\delta_{\Gamma}\), where \(\delta_{\Gamma}\) is the Dirac measure on \(\Gamma\). Thus the inverse task is to locate the defect set \(\Sigma = D = \text{supp}(\gamma - \gamma_0) \subseteq \Gamma\), rather than to reconstruct the whole known carrier \(\Gamma\).

2. Problem II (Full Topology Reconstruction): In the same interface model 1 , the interface \(\Gamma\) may be non-closed and is itself unknown, as is the coefficient \(\gamma\in L^\infty(\Gamma)\). In this case the background operator is the interface-free Laplacian, \(\mathcal{L}_0=-\Delta\), and \(\mathcal{S}\) represents the entire interface jump condition, acting as \(\mathcal{S}u = -\gamma u \delta_{\Gamma}\). The target support is therefore the unknown interface itself, \(\Sigma = \Gamma \subset \Omega\), and the sampling points must be searched over the interior of the domain.

When \(\Gamma\) is not closed, we use the following convention. Choose a closed Lipschitz extension \(\widetilde{\Gamma}\Subset\Omega\) containing \(\Gamma\) and partitioning \(\Omega\), and extend the jump coefficient by \(\widetilde{\gamma}=\gamma\) on \(\Gamma\) and \(\widetilde{\gamma}=0\) on \(\widetilde{\Gamma}\setminus\Gamma\) (and similarly for the background coefficient when used). This zero extension leaves the weak formulation unchanged, since the added part contributes no interface term; equivalently, the artificial part carries zero flux jump. Therefore, the closed-interface well-posedness result in Lemma 1 applies to this extended formulation. Under this formulation, the DSM constructs an index function that is expected to take relatively large values near \(\Sigma\) by exploiting the discrepancy between the observed data and a reference “unperturbed” state \(u_0\). The specific implementations are detailed in the following subsections.

2.1 The DSM Approach↩︎

The Direct Sampling Method (DSM) is a mathematical approach designed to efficiently and directly recover the location of inhomogeneities in a given equation. These inhomogeneities may appear as discontinuities in a coefficient or source term, or as the presence of an interface within the domain. The method operates by considering a suitably chosen reference or background model of the original equation, in which the unknown singular perturbation is removed. Because the problem is overdetermined, the reference equation cannot simultaneously satisfy both sets of Cauchy boundary data. Therefore, the DSM enforces the reference equation to satisfy one prescribed boundary condition, which allows us to capture the discrepancy between the second boundary measurement and the solution to the homogeneous equation.

To quantify this discrepancy, we employ Green’s representation theorem to express the difference between the observed data and the homogeneous solution in terms of Green’s functions. Subsequently, we construct a family of probing functions that is almost orthogonal to the family of Green’s functions under a specific inner product with proper scaling. The index function is then defined as the scaled inner product of the probing functions with the boundary data discrepancy. This almost orthogonal property ensures that the index function attains significantly larger values near the inhomogeneities, effectively identifying their locations.

We now adapt this DSM framework to our specific interface problems. Recall that the boundary flux \(g\) is prescribed on the whole outer boundary \(\partial\Omega\), while the Dirichlet measurement \(f\) is available only on a subset \(L\subset\partial\Omega\). We define a reference state \(u_0\) through a mixed Neumann–Robin background system. On \(\partial\Omega\setminus L\), where there is no Dirichlet measurement, we impose the measured flux \(g\). On the measurement part \(L\), we combine the two Cauchy data through the Robin condition \[\frac{\partial u_0}{\partial n}+\alpha u_0=g+\alpha f, \qquad \alpha>0.\] The parameter \(\alpha\) controls the relative weight assigned to the Dirichlet measurement. In particular, larger values of \(\alpha\) enforce the trace of \(u_0\) to be closer to the measured data \(f\), while smaller values place more emphasis on the Neumann datum.

For Problem I, the known background interface condition is retained in the reference system. Thus \(u_0\) is defined as the solution of 2 below. For Problem II, where the interface itself is unknown, the reference state is taken to be the solution of the interface-free background problem 3 below.

Throughout the paper, we assume that the corresponding reference systems are uniquely solvable for the chosen parameters. In the interface case, this amounts to excluding the same type of exceptional resonant values as in Lemma 1.

After obtaining the reference solution, we extract its trace \(f_0 = u_0|_L\) and exclusively use the Dirichlet discrepancy \(f - f_0\) for data fitting. We construct a family of probing functions \(\eta_{x}:L\rightarrow\mathbb{R}\) such that the inner product between \(f-f_{0}\) and \(\eta_{x}\) serves as an effective probing index. For Problem I, this index locates the defect \(D = \text{supp}(\gamma - \gamma_0) \subset \Gamma\). For Problem II, it identifies the entirely unknown interface \(\Gamma\). Specifically, we expect the index function \[I(x)=\frac{(f-f_{0},\eta_{x})}{S(x)}\] to correctly locate the inhomogeneities for both problems. Here, \((\cdot, \cdot)\) denotes a specific inner product and \(S(x)\) is a scaling function, both of which will be appropriately chosen in the subsequent sections along with the probing functions \(\eta_x\).

2.2 Representations of boundary measurement↩︎

In this subsection, we derive a representation formula for the boundary data discrepancy \(f-f_0\) via the Green’s function under noise-free assumptions. Let \(G_{x}\) be the Green’s function satisfying the following boundary value problem: \[\left\{ \begin{align} -\Delta G_{x}=\delta_{x} \quad in \quad &\Omega,\\ \frac{\partial G_{x}}{\partial n} = 0 \quad on \quad & \partial\Omega\setminus L,\\ \frac{\partial G_{x}}{\partial n}+\alpha G_{x}=0 \quad on \quad & L.\\ \end{align} \right. \label{Gx-def}\tag{3}\] It is well known that the Green’s function \(G_x\) possesses the symmetric property \(G_{x}(y) = G_{y}(x)\) for all \(x, y\in\Omega\).
To apply Green’s second identity, we extend the interface \(\Gamma\) to be a closed surface \(\widetilde{\Gamma}\subset\Omega\) which partitions the domain such that \(\Omega = \Omega_1 \cup \Omega_2\) and \(\Gamma \subseteq \partial\Omega_1\) (as illustrated in Figure [fig:extend95surface]). Define \(\widetilde{\gamma}: \widetilde{\Gamma}\to\mathbb{R}\) by \(\widetilde{\gamma}(x) = \gamma(x)\) on \(\Gamma\), and \(\widetilde{\gamma}(x) = 0\) otherwise. By applying Green’s second identity on the partitioned subdomains \(\Omega_1\) and \(\Omega_2\), we have for \(x\in\text{int}(\Omega_{1})\cup\text{int}(\Omega_{2})\), \[\begin{align} (u-u_{0})(x) &= \int_{\Omega} (G_x \Delta(u-u_0) - (u-u_0)\Delta G_x) dy \\ &= \int_{\Omega_1}(G_x \Delta(u-u_0) - (u-u_0)\Delta G_x) dy + \int_{\Omega_2}(G_x \Delta(u-u_0) - (u-u_0)\Delta G_x) dy\\ &= \int_{\Gamma} \left( G_x \left[ \frac{\partial (u-u_0)}{\partial \nu} \right] - (u-u_0) \left[ \frac{\partial G_x}{\partial \nu} \right] \right) ds_y + \int_{\partial\Omega} \left( G_{x}\frac{\partial (u-u_0)}{\partial n}-(u-u_{0})\frac{\partial G_{x}}{\partial n} \right) ds_y. \end{align}\] Since \(u\), \(u_0\), and the Green’s function \(G_x\) are continuous across the interface \(\Gamma\), the jump term \(\left[ \frac{\partial G_x}{\partial \nu} \right]\) vanishes. Utilizing the transmission conditions \(\left[ \frac{\partial u}{\partial \nu} \right] = \gamma u\) and \(\left[ \frac{\partial u_0}{\partial \nu} \right] = \gamma_0 u_0\), the integral over \(\Gamma\) simplifies directly to the flux jump discrepancy, yielding: \[(u-u_{0})(x) = \int_{\Gamma}(\gamma u-\gamma_{0}u_{0})G_{x} ds_y + \int_{\partial\Omega} \left( G_{x}\frac{\partial (u-u_0)}{\partial n}-(u-u_{0})\frac{\partial G_{x}}{\partial n} \right) ds_y.\] We evaluate the boundary integral over \(\partial\Omega\) by splitting it into \(\partial\Omega \setminus L\) and \(L\). On \(\partial\Omega \setminus L\), both \(u\) and \(u_0\) satisfy the Neumann condition \(\frac{\partial u}{\partial n} = \frac{\partial u_0}{\partial n} = g\), which causes the integral on this portion to vanish. On the measurement boundary \(L\), by substituting \(\frac{\partial u}{\partial n} = g\), \(u = f\), and \(\frac{\partial G_x}{\partial n} = -\alpha G_x\), we obtain: \[\int_{L} \left( G_{x} \left( g - \frac{\partial u_0}{\partial n} \right) - (f-u_{0})(-\alpha G_{x}) \right) ds_y = \int_{L} G_{x} \left( g - \frac{\partial u_{0}}{\partial n} + \alpha(f - u_{0}) \right) ds_y.\] From the Robin boundary condition 2 for the reference system on \(L\), we have \(\frac{\partial u_0}{\partial n} + \alpha u_0 = g + \alpha f\), which implies \(g - \frac{\partial u_0}{\partial n} = -\alpha(f - u_0)\). Substituting this relationship into the equation above, the terms inside the parenthesis exactly cancel out. Thus, the entire boundary integral over \(\partial\Omega\) is zero, and we arrive at: \[(u-u_{0})(x) = \int_{\Gamma}(\gamma(y) u(y)-\gamma_{0}u_{0}(y))G_{x}(y)ds_y.\] By the symmetry of Green’s function \(G_{x}\), we have for \(x\in\text{int}(\Omega_{1})\cup\text{int}(\Omega_{2})\), \[(u-u_{0})(x) = \int_{\Gamma}(\gamma(y) u(y)-\gamma_{0}u_{0}(y))G_{y}(x)ds_y.\] Therefore, by the continuity of \(u-u_{0}\), taking the limit as \(x\) approaches the boundary \(\partial\Omega\) yields the representation for the discrepancy data: \[(f-f_{0})(x) = \int_{\Gamma}(\gamma(y) u(y)-\gamma_{0}u_{0}(y))G_{y}(x)ds_y. \label{step1repre}\tag{4}\] For later use, we denote the effective density in the above representation by \[J_I(y):=\gamma(y)u(y)-\gamma_0(y)u_0(y). \label{eq:effective-density-I}\tag{5}\] Since \(\gamma=\gamma_0\) on \(\Gamma\setminus D\), this density can be decomposed as \[J_I = (\gamma-\gamma_0)u+\gamma_0(u-u_0). \label{eq:effective-density-decomposition-I}\tag{6}\]

Similarly, for Problem II, by applying the exact same procedure and noting that the reference system employs an identical boundary condition matching strategy, the boundary integrals vanish identically. For \(x \in \text{int}(\Omega_1) \cup \text{int}(\Omega_2)\), we directly obtain: \[(u-u_{0})(x) = \int_{\Gamma}\gamma(y) u(y)G_{y}(x)ds_y.\] Again, by the continuity of \(u - u_0\), we obtain the representation formula for \(x \in \partial\Omega\): \[(f-f_{0})(x) = \int_{\Gamma}\gamma(y) u(y)G_{y}(x)ds_y. \label{step2repre}\tag{7}\] At this stage, we have derived the representation formulae for the boundary data discrepancy \(f-f_0\) in both Problem I (4 ) and Problem II (7 ).

2.3 Choices of inner products and probing functions↩︎

This subsection introduces the specific inner products and probing functions, which are crucial for the successful implementation of the DSM. Following the notation in [22], let \(\nabla_{\Gamma}\) denote the surface gradient operator, and let \(-\Delta_{\Gamma}\) denote the Laplace-Beltrami operator (surface Laplacian) on a given curve. We define the inner product \((\cdot,\cdot)\) as the \(H^{s}\) semi-inner product. For a smooth curve \(L\) and \(s\in\mathbb{N}\), this \(H^{s}\) semi-inner product is expressed as: \[(\phi, \psi)_s = \int_{L} (-\Delta_{\Gamma})^s \phi \cdot \psi d\ell, \label{semi-inner-def}\tag{8}\] where \(d\ell\) denotes the arc length element. The subscript \(s\) in (8 ) denotes the regularity order of the Sobolev space \(H^{s}(L)\). We modify the family of probing functions in [17] slightly to suit our boundary condition matching strategy. For a sampling point \(x \in \text{int}(\Omega)\), let \(\omega_x\) solve the following mixed boundary value problem: \[\left\{ \begin{align} -\Delta \omega_{x}=\delta_{x} \quad in \quad &\Omega,\\ \frac{\partial \omega_{x}}{\partial n} = 0\quad on \quad & \partial\Omega\setminus L,\\ \omega_{x}=0 \quad on \quad & L,\\ \end{align} \right. \label{probe-def}\tag{9}\] where \(\delta_x\) is the Dirac delta distribution centered at \(x\). Although Problem I and II represent distinct inverse problems, they share a fundamental mathematical structure as both involve recovering the spatial location of inhomogeneities within the Laplace equation framework. This unified characterization motivates our adoption of a common probing function approach. To maintain notational clarity, we employ identical symbols for the probing formulations in both problems. We define the probing function \(\eta_x\) and the corresponding index function \(I_s\) as: \[\eta_x(y) = \frac{\partial \omega_x}{\partial n}(y) \quad \text{for } y \in L, \quad \text{and} \quad I_s(x) = \frac{(f - f_0, \eta_x)_s}{S(x)}. \label{index-def}\tag{10}\] The spatial domain of interest differs between the two problems. For Problem I (defect detection), the probing is performed over the known interface, meaning \(x \in \Gamma\). For Problem II (full interface reconstruction), the search domain extends to the entire interior region, meaning \(x \in \text{int}(\Omega)\). The scaling function \(S(x)\) and the measurement boundary \(L \subset \partial\Omega\) remain problem-independent. This unified formulation effectively concentrates the sensitivity of the index function on subdomains containing the target inhomogeneities: Either the support \(D = \text{supp}(\gamma - \gamma_0)\) in Problem I, or the unknown interface \(\Gamma\) in Problem II.

2.4 Kernels of the DSM↩︎

This subsection introduces the kernel functions, which constitute the theoretical core of the DSM. Based on the representation formulae (4 ) and (7 ), the scattered potential \(f - f_0\) can be approximated numerically via a quadrature rule as: \[(f-f_{0})(x) \approx \sum_{k} a_{k}G_{x_{k}}(x), \quad x \in \partial\Omega,\] where \(x_k \in \Gamma\) are the discrete sampling points on the interface, and \(a_k\) denote the corresponding quadrature weights that represent the localized intensity of the flux jump discrepancy. Consequently, the probing index \(I_s(x)\) can be approximated as: \[I_{s}(x) \approx \frac{\sum_{k} a_{k}(G_{x_{k}}, \eta_{x})_{s}}{S(x)}.\]

The probing index \(I_s(x)\) is expected to attain significantly larger values when the sampling point \(x\) approaches the inhomogeneities and remains small when \(x\) is distant from them. For Problem I, the inhomogeneities correspond to the defect support \(D = \text{supp}(\gamma - \gamma_0) \subset \Gamma\), while for Problem II, they represent the unknown interface \(\Gamma \subset \Omega\).

The interpretation of the quadrature weights depends on the problem under consideration. For Problem I, the weights approximate the effective density \(J_I\) in 5 ; more precisely, \[a_k \approx J_I(x_k)\Delta s_k, \qquad x_k\in\Gamma,\] where \(\Delta s_k\) denotes the local quadrature weight. By 6 , this density has two parts. The term \((\gamma-\gamma_0)u\) is the direct contribution of the defect and is supported on \(D=\operatorname{supp}(\gamma-\gamma_0)\). The second term \(\gamma_0(u-u_0)\) is caused by the difference between the true solution and the reference solution. Since a local defect may change the solution along the whole interface, this term is not necessarily supported only on \(D\). Therefore, the DSM index for Problem I should be understood as detecting the parts of \(D\) that produce a sufficiently strong signal in the effective density \(J_I\). For a given boundary flux \(g\), some components of \(D\) may generate strong peaks in the index, while others may be weak or partially masked by the background term \(\gamma_0(u-u_0)\).

For Problem II, the corresponding effective density is \[J_{II}(y):=\gamma(y)u(y),\] which is supported on the unknown interface \(\Gamma\). Thus, provided that \(J_{II}\) does not vanish on the relevant interface component for the chosen boundary excitation, the kernel localization property makes the probing index large near that component.

This observation also explains the use of multiple boundary excitations in the numerical experiments. A single Cauchy pair produces one legitimate probing index, but different choices of \(g\) may illuminate different components of \(D\) or \(\Gamma\). When several Cauchy pairs are available, we therefore combine the independently computed single-pair indices, for example through a composite indicator of the form \[I_{\rm comp}(x) = \sum_{m=1}^{N_{\rm data}} \left|\widetilde{I}^{(m)}(x)\right|^2,\] where \(\widetilde{I}^{(m)}\) denotes the normalized probing index obtained from the \(m\)-th Cauchy pair.

Therefore, the effectiveness of the DSM relies on the property that the probing functions \(\{\eta_x\}_{x \in \Omega}\) are almost orthogonal to the family of Green’s functions \(\{G_y|_L\}_{y \in \Omega}\) in the sense of the \(H^s\) semi-inner product. That is, the kernel function \[K(x,y) = \frac{(\eta_{x}, G_{y})_{s}}{S(x)} \label{kernel}\tag{11}\] should have large magnitude when \(x \approx y\) and decay rapidly in magnitude as the distance between \(x\) and \(y\) increases. The analysis of this almost-orthogonal property provides the theoretical foundation for the DSM’s performance, which we further verify through both rigorous proof for a circular domain and numerical visualization in Section 3.

2.5 Alternative characterization of the index function↩︎

In this subsection, we introduce the adjoint method, which is an efficient and accurate technique to make the calculation of index functions less computationally expensive. The implementation of the index function \(I_s(x)\) defined in Section 2.4 requires solving the mixed boundary value problem 9 for the probing function \(\eta_x\) at every sampling point \(x \in \Omega\). In practice, when the sampling grid is dense, this point-by-point approach becomes computationally expensive. To enhance the efficiency of the DSM, we introduce an alternative characterization of the index function using the adjoint method, which allows for the simultaneous evaluation of \(I_s(x)\) across the entire domain.

For a fixed regularity order \(s \ge 0\), we define the adjoint state \(\lambda_s \in H^1(\Omega)\) as the solution to the following adjoint problem: \[\left\{ \begin{align} -\Delta \lambda_{s}=0 \quad &\text{in} \quad \Omega,\\ \frac{\partial\lambda_{s}}{\partial n}=0 \quad &\text{on} \quad \partial\Omega\setminus L,\\ \lambda_{s}=\mathbb{P}(f-f_{0}) \quad &\text{on} \quad L.\\ \end{align} \right. \label{adjoint-def}\tag{12}\] Here, \(\mathbb{P}\) denotes a projection operator with respect to the chosen inner product. In the case of the \(H^s\) semi-inner product, it is given by \[\mathbb{P} = -\frac{(-\Delta_{\Gamma})^{s}}{\delta}, \label{proj-def}\tag{13}\] where \(\delta\) is a constant parameter, and \(-\Delta_{\Gamma}\) is the Laplace-Beltrami operator. Then, applying Green’s second identity, we obtain: \[\lambda_s(x) = \int_\Omega (\omega_x\Delta\lambda_s - \lambda_s\Delta\omega_x) dy = \int_{\partial\Omega} \left(\omega_x\frac{\partial\lambda_s}{\partial n} - \lambda_s\frac{\partial\omega_x}{\partial n}\right) ds_y = \frac{1}{\delta}\int_L (-\Delta_\Gamma)^s(f-f_0)\frac{\partial\omega_x}{\partial n} ds_y = \frac{(f-f_0,\eta_x)_s}{\delta}.\] Therefore, by the definition of the index function in (10 ), we have \[I_{s}(x)=\delta\frac{\lambda_{s}(x)}{S(x)}.\] The constant \(\delta>0\) is introduced only as a scaling parameter, so that the adjoint formulation is consistent with the refinement scheme in Section 4. For the unrefined DSM index, its value has no effect on the peak locations. In the following, for the \(H^1\) semi-inner product, a detailed analysis and calculations regarding the kernel function (11 ) with the specific scaling function \(S(x)\) are presented in Section 3.

3 Almost orthogonal property of kernel functions↩︎

In this section, we provide a detailed theoretical and numerical analysis of the DSM kernel functions. We first verify the almost orthogonal property for the case where \(L = \partial\Omega\), meaning that full boundary measurement data \(f\) are available on the entire boundary. Since numerical experiments indicate that the Sobolev regularity order \(s=1\) is the most suitable choice for the DSM to maintain both accuracy and robustness against noise, we present the corresponding analytical calculations and numerical verifications for \(s=1\) in Section 3.1. Furthermore, in Section 3.2, we conduct a frequency-domain analysis to establish a systematic criterion for selecting the Robin parameter \(\alpha\).

3.1 Theoretical analysis on a circular domain↩︎

We consider the unit disk \(\Omega = B_1(0)\). Let \((r_x, \theta_x)\) and \((r_y, \theta_y)\) be the polar coordinates of \(x\) and \(y\), respectively. For a harmonic function \(v\) in \(B_1(0)\) satisfying: \[\left\{ \begin{align} -\Delta v=0\quad\text{in}\quad B_{1}(0),\\ \frac{\partial v}{\partial n} + \alpha v =h\quad\text{on}\quad\mathbb{S}^{1}, \end{align} \right. \label{circle-harmonic}\tag{14}\] we have: \[v(x) = \int_{\Omega}G_{x}\Delta v - v\Delta G_{x} = \int_{\partial\Omega}G_{x}\frac{\partial v}{\partial n}-v\frac{\partial G_{x}}{\partial n} = \int_{\partial\Omega}G_{x}(h-\alpha v)-v(-\alpha G_{x}) = \int_{\partial\Omega}G_{x}h. \label{compare1}\tag{15}\] Since \(\Delta v = 0\), by using the polar coordinates, we obtain: \[\frac{\partial^{2} v}{\partial r^{2}} + \frac{1}{r}\frac{\partial v}{\partial r}+\frac{1}{r^{2}}\frac{\partial^{2}v}{\partial\theta^{2}}=0.\] Suppose \(v(r,\theta) = R(r)\Theta(\theta)\). By using the separation of variables, we obtain: \[r^{2}R^{''}(r) + rR^{'}(r)-cR(r) = 0,\qquad \Theta^{''}(\theta) + c\Theta(\theta) = 0.\] By the periodical property of \(\Theta\), we find that \(\{e^{in\theta}\}_{n\in\mathbb{Z}}\) is the fundamental solution set of \(\Theta(\theta)\), and \(\{r^{n}\}_{n\in\mathbb{Z}}\cup\{\log r\}\) is the fundamental solution set of \(R(r)\). Hence, the solution can be shown as the following series: \[v(r,\theta) = c+d\log(r) + \sum_{n=1}^{\infty}(d_{n}r^{n}+d_{-n}r^{-n})(c_{n}e^{in\theta}+c_{-n}e^{-in\theta}). \label{general-sol}\tag{16}\] Since \(\{r=0\}\) (the origin) is contained in \(B_1(0)\), singular terms at the origin must be omitted. Hence, \(v\) can be represented by \[v(r,\theta) = \sum_{n\in\mathbb{Z}}a_{n}r^{|n|}e^{in\theta} \label{sol-0include}\tag{17}\] with coefficients \(\{a_n\}_{n\in\mathbb{Z}}\). Using (14 ) and (17 ), we have \[h(\theta) = \sum_{n\in\mathbb{Z}}(|n|+\alpha)a_{n}e^{in\theta}.\] Therefore, the Fourier coefficients of \(h\in L^{2}(\mathbb{S}^{1})\) can be given as \[(|n|+\alpha)a_{n} = \frac{1}{2\pi}\int_{0}^{2\pi}h(\theta)e^{-in\theta}d\theta,\quad\text{hence}\quad a_{n} = \frac{1}{2\pi(|n|+\alpha)}\int_{0}^{2\pi}h(\theta)e^{-in\theta}d\theta.\] Consequently, we have \[v(r_{x}, \theta_{x}) = \frac{1}{2\pi}\sum_{n\in\mathbb{Z}}\frac{1}{|n|+\alpha}r_{x}^{|n|}e^{in\theta_{x}}\int_{0}^{2\pi}h(\theta)e^{-in\theta}d\theta,\] which is equivalent to \[v(x) = \int_{\partial\Omega}h(y)\sum_{n\in\mathbb{Z}}\frac{1}{2\pi(|n|+\alpha)}r_{x}^{|n|}e^{in(\theta_{x}-\theta_{y})}ds_{y}. \label{compare2}\tag{18}\] Comparing (15 ) and (18 ), the Green’s function \(G_{x}(y)\) for \(y\in\mathbb{S}^{1}\) can be written as \[G_{x}(y) = \sum_{n\in\mathbb{Z}}\frac{1}{2\pi(|n|+\alpha)}r_{x}^{|n|}e^{-in\theta_{x}}e^{in\theta_{y}}. \label{Gx-Fourier}\tag{19}\] Using a similar technique, the probing function \(\eta_{x}\) on \(\mathbb{S}^{1}\) can be evaluated as \[\eta_{x}(y) = \frac{\partial w_{x}}{\partial n} = -\frac{1}{2\pi}\sum_{n\in\mathbb{Z}}r_{x}^{|n|}e^{-in\theta_{x}}e^{in\theta_{y}}. \label{etax-Fourier}\tag{20}\] Thus, the inner product of \(\eta_{x}\) and \(G_{z}\) is given by \[(\eta_{x}, G_{z})_{1} = -\frac{1}{\pi}\text{Re}\left(\sum_{n=1}^{\infty}\frac{n^{2}}{n+\alpha}r_{x}^{n}r_{z}^{n}e^{in(\theta_{x}-\theta_{z})}\right).\] Furthermore, the \(H^{0}\) and \(H^{1}\) norms of \(\eta_{x}\) and \(G_{x}\) are computed as \[|\eta_{x}|_{H^{0}}^{2} = \frac{1}{2\pi}\sum_{n\in\mathbb{Z}}r_{x}^{2|n|} =\frac{1}{2\pi}\frac{1+r_{x}^{2}}{1-r_{x}^{2}},\quad |\eta_{x}|_{H^{1}}^{2} = \frac{1}{2\pi}\sum_{n\in\mathbb{Z}}n^{2}r_{x}^{2|n|} = \frac{1}{\pi}\frac{r_{x}^{2}(1+r_{x}^{2})}{(1-r_{x}^{2})^{3}}.\] \[|G_{x}|_{H^{0}}^{2} = \frac{1}{2\pi}\sum_{n\in\mathbb{Z}}\frac{1}{(|n|+\alpha)^{2}}r_{x}^{2|n|}, \quad |G_{x}|_{H^{1}}^{2} = \frac{1}{2\pi}\sum_{n\in\mathbb{Z}}\frac{n^{2}}{(|n|+\alpha)^{2}}r_{x}^{2|n|}.\] Intuitively, we take the scaling function as \[S(x) = |\eta_{x}|_{H^{1}}^{1/2}|G_{x}|_{H^{1}}^{1/2}.\] Then, for a fixed \(z\in\Omega\), the kernel function is expressed as \[K(x,z) = \frac{(\eta_{x}, G_{z})_{1}}{|\eta_{x}|_{H^{1}}^{1/2}|G_{x}|_{H^{1}}^{1/2}} = -c\frac{\sum_{n=1}^{\infty}\frac{n^{2}}{n+\alpha}r_{x}^{n}r_{z}^{n}\cos(n(\theta_{x}-\theta_{z}))}{\left(\frac{r_{x}^{2}(1+r_{x}^{2})}{(1-r_{x}^{2})^{3}}\right)^{1/4}\left(\sum_{n\in\mathbb{Z}}\frac{n^{2}}{(|n|+\alpha)^{2}}r_{x}^{2|n|}\right)^{1/4}}. \label{kernel95Fourier}\tag{21}\] Here, \(c>0\) is a constant independent of \(r_x\), \(r_z\) and \(\theta_x\), \(\theta_z\). The minus sign in (21 ) is a global sign and does not affect peak locations; in the following numerical plots and in Appendix 7.2, we use the sign-normalized kernel and still denote it by \(K\) for simplicity. We can see from (21 ) that for fixed \(r_x\) and \(r_z\), this kernel \(K\) attains its maximum when \(\theta_x = \theta_z\). Therefore, for each fixed \(r_{z}\), we only need to verify whether \(K(x,z)\) reaches its maximum at \(r_{x}\approx r_{z}\) under the assumption that \(\theta_x = \theta_z\). A rigorous proof establishing this almost orthogonal property for the kernel (21 ) with \(\alpha=1\) is provided in Appendix 7.2. In this section, we first verify it numerically.

We compute the partial sum of the first 100 terms of the Fourier series for \(r_{z}\in(0,1)\)and plot the results, where the y-axis represents \(\text{argmax}_{r_{x}}K(x,z)\) and the x-axis represents \(r_z\). Figure 4 shows the argmax graphs of the kernel function \(K\) for different \(\alpha\), and Figure ¿fig:circle? displays the profiles of \(K(x,z)\) for two different \(z\) in \(B_{1}(0)\).

For the square domain \(\Omega = [-1,1]\times[-1,1]\), we adopt the same scaling function. To enhance computational efficiency, we avoid directly evaluating the Green’s functions \(w_{x}\) and \(G_{x}\) at each point \(x\). We follow [19] by utilizing the fundamental solution \(\Phi_{x}(y) = \frac{1}{2\pi}\log |x-y|\) to define the modified scaling function as \[S(x) = \left|\frac{\partial\Phi_{x}}{\partial n}\right|_{H^{0}}^{1/2}\left|\frac{\partial\Phi_{x}}{\partial n}\right|_{H^{1}}^{1/2}, \label{scaling}\tag{22}\] which can be computed directly. Here, we employ the following approximations: \[|G_{x}|_{H^{1}}\approx\left|\frac{\partial\Phi_{x}}{\partial n}\right|_{H^{0}},\quad |\eta_{x}|_{H^{1}}\approx\left|\frac{\partial\Phi_{x}}{\partial n}\right|_{H^{1}}.\] To verify the almost orthogonal property of the modified kernel function in the square domain \(\Omega = [-1,1]\times[-1,1]\), we plot a 3D surface \((z_0,D(z_0))\) for \(z_0 \in \Omega\), where \(D(z_0)\) is defined as \[D(z_0) = \left|\text{argmax}_{z} \frac{(\eta_{z}, G_{z_{0}})}{|\frac{\partial\Phi_{z}}{\partial n}|_{H^{0}}^{1/2}|\frac{\partial\Phi_{z}}{\partial n}|_{H^{1}}^{1/2}} - z_{0}\right|. \label{distance}\tag{23}\] The surface graph is shown in Figure 5. We can see from Figure 5 that when \(K(z,z_{0})\) attains its maximum in \(\Omega=[-1,1]\times[-1,1]\), the distance between \(z\) and \(z_{0}\) remains strictly less than 0.1. This consistently tight bound numerically verifies the almost orthogonal property of the kernel \(K\) utilizing the modified scaling function (22 ).

Figure [square-kernel1] illustrates the kernel function \(K(x,z)\) for two different source points \(z\) (marked by dark blue dots). The heatmap demonstrates that the maximum values (represented by red regions) are sharply localized around each \(z\) while exhibiting rapid decay away from them. This concentration behavior numerically confirms the almost orthogonal property of the kernel, with \(K(x,z)\) effectively serving as a spatial indicator that peaks near \(x=z\) and decreases with increasing distance from the source point. This concludes the verification of our chosen scaling function \(S(x)\) for the Sobolev regularity order \(s=1\).

3.2 Frequency-domain analysis and the selection of the Robin parameter \(\alpha\)↩︎

The Robin parameter \(\alpha\) is pivotal in the reference system, as it governs the balance between high-frequency spatial resolution and the preservation of low-frequency macroscopic structures. Because the scaling function \(S(x)\) in the DSM acts as a data-independent spatial normalizer, the spatial profile and the spectral energy distribution of the index function \(I_{1}(x)\) are uniquely determined by its unscaled numerator \((f^{\epsilon}-f_{0}^{\epsilon},\eta_{x})_{1}\). To establish a systematic criterion for selecting \(\alpha\), we analyze the spectral behavior of this numerator and how \(\alpha\) weights the signal and noise across different frequencies using Fourier series.

In this subsection, we distinguish between the exact and noisy reference states. Let \(u_0^{\rm ex}\) be the reference solution generated from the exact boundary value \(f_{\rm exact}\), and let \(u_0^{\epsilon}\) be the reference solution generated from the noisy boundary value \(f^{\epsilon}\). We write \(f_0^{\rm ex}=u_0^{\rm ex}|_L\) and \(f_0^{\epsilon}=u_0^{\epsilon}|_L\).

This inequality highlights the ill-posedness of the inverse problem: high-frequency structural information of the source is subjected to exponential decay \(R_{max}^{|n|}\) before reaching the boundary. In practice, the exact data is polluted by noise \(\epsilon\), yielding the noisy Fourier coefficients \(\hat{f}^\epsilon_n = \hat{f}_{{\rm exact},n} + \epsilon_n\). The exact and noisy reference solutions satisfy \[\widehat{u_0^{\rm ex}}_n = \frac{\hat{g}_n + \alpha \hat{f}_{{\rm exact},n}}{|n|+\alpha}, \qquad \widehat{u_0^{\epsilon}}_n = \frac{\hat{g}_n + \alpha \hat{f}^\epsilon_n}{|n|+\alpha}.\] Hence the Fourier coefficients of the noisy input discrepancy \(d^{\epsilon}:=f^{\epsilon}-u_0^{\epsilon}\) satisfy \[\hat{d}^{\epsilon}_n = \widehat{(f_{\rm exact} - u_0^{\rm ex})}_n + \frac{|n|}{|n| + \alpha} \epsilon_n.\] The \(H^1\) inner product in the index function introduces a multiplier \(|n|^2\). Following (20 ), the numerator of the probing index becomes, up to the same global sign normalization as in Section 3.1, \[\text{Numerator}^{\epsilon}(x) = \frac{1}{2\pi} \sum_{n \in \mathbb{Z}} |n|^2 \widehat{(f_{\rm exact} - u_0^{\rm ex})}_n r_x^{|n|} e^{in\theta_x} + \frac{1}{2\pi} \sum_{n \in \mathbb{Z}} \frac{|n|^3}{|n| + \alpha} \epsilon_n r_x^{|n|} e^{in\theta_x}.\] Based on the physical decay bound (?? ), we define the effective signal weight \(SW(n)\) and noise weight \(NW(n)\) as: \[SW(n) = \frac{|n|^2}{|n| + \alpha} R_{max}^{|n|}, \quad NW(n) = \frac{|n|^3}{|n| + \alpha}. \label{eq:spectral95weights}\tag{24}\]

Regarding the selection of \(\alpha\), numerical experiments indicate that utilizing a small parameter such as \(\alpha = 0.1\) yields an overly flat index profile. When \(\alpha \to 0\), the signal weight degenerates to \(SW(n) \approx |n| R_{max}^{|n|}\). Without sufficient penalization of fundamental modes, the signal energy is excessively concentrated at the lowest frequencies (e.g., \(n=1,2\)), which spatially manifests as broad, smooth waves. Moreover, a small \(\alpha\) leads to an abnormally large scaling factor \(S(x, \alpha)\), uniformly squashing the amplitude of the index function.

To systematically evaluate the optimal \(\alpha\), Table 1 compares these weights at low-frequency (\(n=2\)) and mid-frequency (\(n=5\)) modes for \(R_{max} = 0.5\).

For Problem I, localizing sharp flux jumps necessitates isolating higher-frequency features. Table 1 shows that \(\alpha=10\) shifts the signal energy centroid toward higher frequencies. Specifically, the high-to-low ratio \(SW(5)/SW(2)\) increases from \(0.39\) (when \(\alpha=1\)) to \(0.63\). By damping the absolute high-frequency noise \(NW(5)\), this spectral tilt achieves high spatial resolution for local anomalies. Conversely, for Problem II, reconstructing an unknown macroscopic topology relies on the survival of low-frequency modes. In contrast to the high-pass effect of \(\alpha=10\), a moderate \(\alpha=1\) yields a smaller ratio \(SW(5)/SW(2) = 0.39\). This ensures the spectral centroid remains anchored at lower frequencies, which is necessary to recover the macroscopic geometry against the \(R_{max}^{|n|}\) decay.

4 A refinement approach for DSM↩︎

In this section, the refinement approach is presented and tested only for the full Dirichlet measurement case \(L=\partial\Omega\). Experiments show that the numerical results of the DSM with the \(H^{0}\) semi-inner product are more robust to noise but less precise than those with the \(H^{1}\) semi-inner product. We then present a refinement approach to enhance the performance of the DSM with the \(H^{0}\) semi-inner product. Using Green’s second identity and equations (1 ) (2)(3)(9 )(12 ), we have for Problem I: \[\lambda_{0}(x)=(f-f_{0}, \eta_{x})_{0} = -(u-u_{0})(x)+\int_{\Gamma}(\gamma u-\gamma_{0}u_{0})\omega_{x},\quad x\in\text{int}(\Omega). \label{lam-step1-rep}\tag{25}\] Similarly, for Problem II, we obtain: \[\lambda_{0}(x)=(f-f_{0}, \eta_{x})_{0} = -(u-u_{0})(x)+\int_{\Gamma}\gamma u\omega_{x},\quad x\in\text{int}(\Omega). \label{lam-step2-rep}\tag{26}\] Note that in numerical tests, these two equations are valuable tools for verifying the correct solution of the forward model problem (1 ). If a better approximation \(\hat{u}\) of the true solution \(u\) is available compared to the reference solution \(u_{0}\), we can sharpen the numerical results of the DSM using the \(H^{0}\) semi-inner product. To find a better approximation \(\hat{u}\) to replace \(u_0\) in (25 ) and (26 ), we propose the following iterative algorithm.

1. Initialize \(\lambda^{0}(x)=0\), and set a constant \(\delta > 0\).

2. For \(k = 0, 1\):
   For Problem I, solve \[\left\{ \begin{align} -\Delta u_{k}=-\lambda^{k} \quad \text{in} \quad &\Omega\setminus\Gamma,\\ \left[\frac{\partial u_{k}}{\partial \nu}\right] = \gamma_{0}u_{k} \quad \text{on}\quad &\Gamma, \\ \frac{\partial u_{k}}{\partial n}+\alpha u_{k}=g+\alpha f \quad \text{on} \quad & L,\\ \frac{\partial u_{k}}{\partial n} = g\quad \text{on} \quad &\partial\Omega\setminus L. \end{align} \right.\] for Problem II, solve \[\left\{ \begin{align} -\Delta u_{k}=-\lambda^{k} \quad \text{in} \quad &\Omega,\\ \frac{\partial u_{k}}{\partial n}+\alpha u_{k}=g+\alpha f \quad \text{on} \quad & L,\\ \frac{\partial u_{k}}{\partial n} = g\quad \text{on} \quad &\partial\Omega\setminus L, \end{align} \right. \label{u95k95equation}\tag{27}\] to obtain \(f_{k} = u_{k}|_{\partial\Omega}\). Then, solve \[\left\{ \begin{align} -\Delta \lambda^{k+1}=0 \quad \text{in} \quad &\Omega,\\ \lambda^{k+1}=-\frac{f-f_{k}}{\delta}\quad \text{on} \quad &\partial\Omega,\\ \end{align} \right.\] to obtain \(\lambda^{k+1}\), and update \(k \leftarrow k+1\).

3. Compute the modified probing index \[\hat{I}_{0}(x) = \frac{\lambda^{2}(x)}{S(x)} .\label{index-refine-def}\tag{28}\] To satisfy the almost orthogonal property, we define the scaling function for the \(H^0\) semi-inner product as \[S(x) = \left|\frac{\partial\Phi_{x}}{\partial n}\right|_{L^{2}}.\]

From this iterative scheme, we can utilize Green’s second identity to show that for Problem I, \[\lambda^{2}(x) = \frac{1}{\delta}\left(-(u-u_{1})(x) + \int_{\Omega}\lambda^{1}\omega_{x} + \int_{\Gamma}(\gamma u-\gamma_{0}u_{1})\omega_{x}\right). \label{lamk-step1-rep}\tag{29}\] Similarly, for Problem II, we have: \[\lambda^{2}(x) = \frac{1}{\delta}\left(-(u-u_{1})(x) + \int_{\Omega}\lambda^{1}\omega_{x} + \int_{\Gamma}\gamma u\omega_{x}\right). \label{lamk-step2-rep}\tag{30}\] Comparing (29 )(30 ) with (25 )(26 ), the \(u_{0}\) terms are replaced by \(u_{1}\), and an additional term \(\int_{\Omega}\lambda^{1}\omega_{x} dy\) is introduced. This new term is continuously differentiable due to the elliptic regularity of Poisson’s equation. The numerical results presented in Section 5.3 demonstrate that this refinement enhances the performance of the DSM with the \(H^0\) semi-inner product while maintaining its excellent robustness against noise.

To theoretically justify the proposed refinement, we perform the following calculations for the case where \(L = \partial\Omega\), meaning that the full measurement data \(f\) is available on the entire outer boundary. For Problem I, using the calculations in Section 2.6 and equation (4 ), we obtain: \[\lambda^{1}(x) = \frac{1}{\delta}\int_{\partial\Omega}(f-f_0)(y)\eta_x(y)ds_y = \frac{1}{\delta}\int_{\partial\Omega}\left(\int_{\Gamma}(\gamma(z)u(z)-\gamma_0u_0(z))G_{z}(y)ds_z\right)\eta_{x}(y)ds_y.\] Then, by Fubini’s theorem and the definition of the kernel function (11 ), the original index function is given by: \[I_0(x) = \frac{\lambda^{1}(x)}{S(x)} = \frac{1}{\delta S(x)}\int_{\Gamma}(\gamma(z)u(z)-\gamma_0u_0(z))K(x,z)S(x)ds_z = \frac{1}{\delta}\int_{\Gamma}(\gamma(z)u(z)-\gamma_0u_0(z))K(x,z)ds_z . \label{I95step195repre}\tag{31}\] Using a technique similar to that in Section 2.3 along with the refinement scheme, we expand the difference: \[\begin{align} (u-u_{1})(x) &= \int_{\Omega}-(u-u_{1})\Delta G_{x} dy + \int_{\Omega_{1}}G_{x}\Delta(u-u_1) dy + \int_{\Omega_{2}}G_{x}\Delta(u-u_1) dy + \int_{\Omega}G_x\Delta u_1 dy \\ &=\int_{\Gamma}(\gamma(y) u(y)-\gamma_0 u_1(y))G_x(y)ds_y + \int_{\Omega}G_x(y)\lambda^{1}(y)dy. \end{align}\] Subsequently, according to the refinement scheme, we have: \[\begin{align} \lambda^{2}(x) &= \frac{1}{\delta}\int_{\partial\Omega}(f-f_1)(y)\eta_x(y)ds_y = \frac{1}{\delta}\int_{\partial\Omega}(u-u_1)(y)\eta_x(y)ds_y \\ &= \frac{1}{\delta}\int_{\Gamma}(\gamma(y) u(y)-\gamma_0 u_1(y))S(x)K(x,y)ds_y + \frac{1}{\delta}\int_{\Omega}\lambda^{1}(y)S(x)K(x,y)dy. \end{align}\] Therefore, the refined index function can be expressed as: \[\begin{align} \hat{I}_0(x) &= \frac{\lambda^{2}(x)}{S(x)} = \frac{1}{\delta}\int_{\Gamma}(\gamma(y) u(y)-\gamma_0 u_1(y))K(x,y)ds_y + \frac{1}{\delta}\int_{\Omega}\lambda^{1}(y)K(x,y)dy\\ &=\frac{1}{\delta}\int_{\Gamma}(\gamma(y) u(y)-\gamma_0 u_1(y))K(x,y)ds_y + \frac{1}{\delta^2}\int_{\Gamma}(\gamma(y) u(y)-\gamma_0 u_0(y))\left(\int_{\Omega}S(z)K(y,z)K(x,z)dz\right)ds_y. \label{hat95I95step195repre} \end{align}\tag{32}\] Similarly, for Problem II, the original index function is represented as: \[I_0(x) = \frac{\lambda^{1}(x)}{S(x)} = \frac{1}{\delta S(x)}\int_{\Gamma}\gamma(z)u(z)K(x,z)S(x)ds_z = \frac{1}{\delta}\int_{\Gamma}\gamma(z)u(z)K(x,z)ds_z . \label{I95step295repre}\tag{33}\] Furthermore, the refined index function is expressed as: \[\begin{align} \hat{I}_0(x) &= \frac{\lambda^{2}(x)}{S(x)} = \frac{1}{\delta}\int_{\Gamma}\gamma(y)u(y)K(x,y)ds_y + \frac{1}{\delta}\int_{\Omega}\lambda^{1}(y)K(x,y)dy\\ &=\frac{1}{\delta}\int_{\Gamma}\gamma(y)u(y)K(x,y)ds_y + \frac{1}{\delta^2}\int_{\Gamma}\gamma(y)u(y)\left(\int_{\Omega}S(z)K(y,z)K(x,z)dz\right)ds_y .\label{hat95I95step295repre} \end{align}\tag{34}\] Comparing (31 )(32 ) and (33 )(34 ), we observe that the refinement term \(\int_{\Omega}S(z)K(y,z)K(x,z)dz\) enhances the original index function through a fundamental mechanism: The product \(K(y,z)K(x,z)\) effectively reinforces the almost orthogonal property of the kernel function. This reinforcement occurs because the kernel product amplifies coherence when both \(x\) and \(y\) approach the true inhomogeneity location \(z\), and accelerates decay when either \(x\) or \(y\) moves away from \(z\). Consequently, the modified index function \(\hat{I}_0(x)\) exhibits superior localization properties while maintaining the inherent robustness of the \(H^0\) formulation.

5 Numerical results↩︎

In this section, we first present numerical results using the \(H^{1}\) semi-inner product for Problem I and Problem II without the refinement approach. The corresponding figures are shown in Sections 5.1 and 5.2. Subsequently, the numerical results of the refinement approach utilizing the \(H^{0}\) semi-inner product are presented in Section 5.3. We first outline the computation procedure without refinement.

First, we designed a specific interface \(\Gamma\subset\Omega = [-1,1]\times[-1,1]\) and a function \(g:\partial\Omega\rightarrow\mathbb{R}\) as the boundary flux. Second, we applied the finite difference method to solve 1 , supplemented with the Neumann boundary condition \(\partial u/\partial n=g\) on \(\partial\Omega\), to obtain the exact boundary measurement \(f=u|_L\). To test the robustness of the DSM, we added random noise to \(f\) and obtained \(f^{\epsilon}\), so that \(f^{\epsilon} = (1+\epsilon\xi)f\). Here, \(\xi\) is a uniformly distributed random variable in \([-1,1]\), and \(\epsilon\) is the noise level, which is specified in each figure caption. In our numerical tests, we took \(\gamma\) as a piecewise continuous function. Third, we solved the reference system (2)(3) using the boundary flux \(g\) and the noisy measurement data \(f^{\epsilon}\) to obtain the reference solution \(u_{0}\). Fourth, using the reference data \(f_{0} = u_{0}|_{L}\), we solved the adjoint system (12 ) to obtain the adjoint function \(\lambda_{s}\). Finally, we computed our desired probing index \(I_{s}(x)\): \[I_{s}(x) = \frac{\lambda_{s}(x)}{S(x)}.\] Here, the scaling function \(S\) was computed by (22 ) for \(s=1\). To normalize our probing index, we plotted the graph of \(\widetilde{I}_{1}(x)\) for our numerical results, where \[\widetilde{I}_{1}(x)= \left({I_{1}(x)}/{\max(\sup I_{1}(x),10^{-6})}\right)^2 \quad \text{if } I_{1}(x)\geq0,\qquad \left({I_{1}(x)}/{\min(\inf I_{1}(x),-10^{-6})}\right)^2 \quad \text{if } I_{1}(x)<0.\]

For Problem I, the probing index \(I_1\) is defined on \(\Gamma\) and is set to 0 otherwise. We set the background value \(\gamma_{0}=1\), where \(\gamma\) exhibited inhomogeneity on \(\Gamma\). On the graph of the normalized probing index \(\widetilde{I}_{1}\), we plot small red marks at the location of \(D=\text{supp}(\gamma - \gamma_{0})\subset \Gamma\). Ideally, these marks should coincide precisely with the peaks of the probing index. For each interface \(\Gamma\), 3 pairs of Cauchy data \((f,g)\) were provided, yielding three corresponding normalized probing indices \(\widetilde{I}_{1}\). We emphasize that each probing index is generated from one single Cauchy pair. The use of three boundary excitations in the numerical experiments is a fusion strategy: each index provides partial information, and the composite index aggregates the visible components from different measurements. We conducted numerical tests using both full boundary data \(f^{\epsilon}\) and \(75\%\) partial data (utilizing measurements from only three sides of \(\partial\Omega\)). We added \(4\%\) noise to the full data and \(2\%\) noise to the partial data. We observe that the probing indices \(\widetilde{I}_{1}\) can reliably capture at least one location of \(D\subset\Gamma\) for each boundary flux \(g\). By combining results from three distinct fluxes \(g\), all inhomogeneous regions \(D\) were successfully identified. The results for Problem I are presented in Section 5.1.

For Problem II, we assumed \(\gamma(x)\) to be piecewise continuous on \(\Gamma\). As in Problem I, three pairs of Cauchy data \((f,g)\) were provided for each interface, and we investigated both the full and \(75\%\) partial data cases under identical noise levels (\(4\%\) and \(2\%\), respectively). Remarkably, despite the inherent ill-posedness of Problem II characterized by multiple unknown interfaces and limited noisy data, our DSM approach demonstrated robust performance. Each probing index \(\widetilde{I}_{1}\) successfully captured at least one interface segment, indicated by larger index values near the true locations. By integrating the three individual results, a comprehensive and accurate reconstruction of \(\Gamma\) was achieved. The results for Problem II are presented in Section 5.2.

For both Problem I and Problem II, we made use of a \(181\times 181\) uniform mesh to perform our finite difference PDE solver to solve for the true solution \(u\). To avoid inverse crime, we used a \(121\times 121\) uniform mesh alternatively to solve \(u_{0}\) and \(\lambda\) using the interpolated value of \(f^{\epsilon}\). The values of \(\alpha\) were set be 10 and 1 for Problem I and Problem II, respectively. The settings for the examples are listed as follows.

In Section 5.1, we detail the specific numerical examples for Problem I. Throughout these examples, we set the Robin parameter to \(\alpha=10\).
Example 1: Non-closed interface. \(\Gamma = [-0.766, 0.5]\times\{-0.5\}\cup \{-0.766\}\times[-0.5,0.233]\). \[g_{1}(x,y) = |x+y+1|,\quad g_{2}(x,y) = 2,\quad g_{3}(x,y) = \cos(2\pi y)\] \[\gamma(x,y) = \left\{ \begin{align} 10+10x,& \quad y = -0.5, |x-0.2|\leq 0.05\\ 1,& \quad\text{else} \end{align} \right.\] Example 2: Closed interface. \(\Gamma = \partial\Omega_{1}\), where \(\Omega_{1} = [-0.766, 0.5]\times[-0.766, 0.733]\). \[g_{1}(x,y) = |x+y+1|,\quad g_{2}(x,y) = \sin(\pi x),\quad g_{3}(x,y) = \cos(2\pi y) - x\] \[\gamma(x,y) = \left\{ \begin{align} 5,& \quad y = -0.766, |x-0.2|\leq 0.05 \\ 5,& \quad y = 0.733, |x-0.3|\leq 0.05\\ 5,& \quad x = -0.766, |y-0.1|\leq 0.05\\ 1,& \quad\text{else} \end{align} \right.\] Example 3: Non-closed interface. \(\Gamma = [-0.733, 0.767]\times\{0.5\}\). \[g_{1}(x,y) = \sin(2\pi(x+y)),\quad g_{2}(x,y) = \sin(4\pi(x+y)),\quad g_{3}(x,y) = \sin(\pi(x+y))\] \[\gamma(x,y) = \left\{ \begin{align} 10,& \quad y = 0.5, |x-0.5|\leq 0.05 \\ 10,& \quad y = 0.5, |x+0.5|\leq 0.05 \\ 1,& \quad\text{else} \end{align} \right.\] Example 4: Non-closed interface. \(\Gamma = [-0.266, 0.733]\times\{0.5\}\cup\{-0.5\}\times[-0.766, 0.233]\). \[g_{1}(x,y) = \sin(2\pi(x+y)),\quad g_{2}(x,y) = \sin(4\pi(x+y)),\quad g_{3}(x,y) = \sin(\pi(x+y))\] \[\gamma(x,y) = \left\{ \begin{align} 10,& \quad y = 0.5, |x-0.3|\leq 0.05 \\ 10,& \quad x = -0.5, |y+0.2|\leq 0.05 \\ 1,& \quad\text{else} \end{align} \right.\] Example 5: \(75\%\) data available. The same setting as in Example 1. The unavailable data part \(\partial\Omega\setminus L=\{1\}\times[-1,1]\).
Example 6: \(75\%\) data available. The same setting as in Example 3. The unavailable data part \(\partial\Omega\setminus L=\{1\}\times[-1,1]\).
Example 7: \(75\%\) data available. The same setting as in Example 4. The unavailable data part \(\partial\Omega\setminus L=\{1\}\times[-1,1]\).
In Section 5.2, we detail the specific numerical examples for Problem II. Throughout these examples, we set the Robin parameter to \(\alpha=1\). For Examples 1 and 2, the interfaces consisted of the boundaries of small rectangles. In Example 1, \(\gamma=2\) on all four edges of each rectangle. In Example 2, for each rectangle, \(\gamma=2\) on the left and right edges, \(\gamma=-1\) on the top edge, and \(\gamma=3\) on the bottom edge; corner values are immaterial. For Examples 3 and 4, the interfaces were composed of several separate line segments.
Example 1: Closed interface. \(\Gamma = \partial\Omega_{1}\cup \partial\Omega_{2}\cup \partial\Omega_{3}\), where \(\Omega_{1} = [-0.083, 0.167]\times[-0.667, -0.417]\), \(\Omega_{2} = [-0.667, -0.333]\times[0.333, 0.583]\), \(\Omega_{3} = [0.333, 0.617]\times[0.333, 0.667]\); \(\gamma=2\) on each \(\partial\Omega_j\). \[g_{1}(x,y) = 2,\quad g_{2}(x,y) = |x+y+1|,\quad g_{3}(x,y) = |\cos(x+y)|\] Example 2: Closed interface. \(\Gamma = \partial\Omega_{1}\cup \partial\Omega_{2}\), where \(\Omega_{1} = [-0.667, -0.333]\times[-0.667, -0.333]\), \(\Omega_{2} = [0.333, 0.667]\times[0, 0.333]\). For each \(\Omega_j=[a_j,b_j]\times[c_j,d_j]\), \(\gamma=2\) on \(\{a_j,b_j\}\times[c_j,d_j]\), \(\gamma=-1\) on \([a_j,b_j]\times\{d_j\}\), and \(\gamma=3\) on \([a_j,b_j]\times\{c_j\}\). \[g_{1}(x,y) = \sin(\pi(x+y)),\quad g_{2}(x,y) = |x+y+1|,\quad g_{3}(x,y) = \sin(2\pi(x+y))\] Example 3: Non-closed interface. \(\Gamma = [-0.75, -0.167]\times\{-0.667\}\cup [0.2, 0.667]\times\{0.333\}\), \(\gamma(x,y) = 5\) on \(\Gamma\). \[g_{1}(x,y) = x+y,\quad g_{2}(x,y) = \sin(\pi(x+y)),\quad g_{3}(x,y) = \sin(2\pi(x+y))\] Example 4: Non-closed interface. \(\Gamma = \{-0.167\}\times[-0.667, 0.667]\). \[\gamma(x,y) = \left\{ \begin{align} -10,& \quad |y|\leq 0.4 \\ 10,& \quad\text{else} \end{align} \right.\] \[g_{1}(x,y) = 2,\quad g_{2}(x,y) = |\sin(\pi(x+y))|,\quad g_{3}(x,y) = |\cos(\pi(x+y))|\] Example 5: \(75\%\) data available. The same setting as in Example 2. The unavailable data part \(\partial\Omega\setminus L=[-1,1]\times\{1\}\).
Example 6: \(75\%\) data available. The same setting as in Example 3. The unavailable data part \(\partial\Omega\setminus L=[-1,1]\times\{1\}\).
Example 7: \(75\%\) data available. The same setting as in Example 4. The unavailable data part \(\partial\Omega\setminus L=\{1\}\times[-1,1]\).

5.1 Problem I↩︎

For Problem I with full data available on the outer boundary, we present four numerical examples to demonstrate the capability of the proposed DSM in reconstructing flux jump locations on a known interface under four scenarios: (i) a single jump on a non-closed interface; (ii) three jumps on a closed interface; (iii) two jumps on a non-closed interface; and (iv) two jumps distributed across two non-closed interfaces. For the partial data case (\(75\%\) available), we investigate scenarios (i), (iii), and (iv). As illustrated in Figure 6, a single pair of Cauchy data was sufficient for the \(H^1\) probing index to accurately capture at least one jump location. More importantly, by combining the probing indices derived from multiple distinct pairs of Cauchy data, we successfully identified every jump location, thereby achieving a complete and highly accurate reconstruction of the defect distribution.

To quantitatively assess the performance of the proposed DSM and to facilitate an objective evaluation of the probing indices, we introduce two quantitative measures: Mean Localization Error (MLE) and Contrast-to-Noise Ratio (CNR).

To comprehensively evaluate the reconstruction performance across all available measurements, we first define the composite evaluation index \(I_{\text{comp}}(x)\) as the sum of the squared normalized probing indices obtained from the three distinct pairs of Cauchy data. That is, \[I_{\text{comp}}(x) = \sum_{i=1}^3 |\widetilde{I}^{(i)}(x)|^2,\] where \(\widetilde{I}^{(i)}(x)\) represents either the refined or unrefined normalized index for the \(i\)-th measurement. Both the MLE and the CNR are evaluated based on this composite index.

The Mean Localization Error (MLE) evaluates the spatial accuracy of the overall reconstruction. Let \(x_j^{\dagger}\) be the true geometric center of the \(j\)-th flux jump. The reconstructed peak \(\hat{x}_j\) is identified by searching for the local maximum of the composite index \(I_{\text{comp}}(x)\) within a neighborhood of radius \(R_{\text{search}} = 0.2\): \[\hat{x}_j = \arg\max_{\substack{x \in \Gamma_h \\ \|x - x_j^{\dagger}\| \le 0.2}} I_{\text{comp}}(x).\] With the defect radius fixed at \(0.05\), the single-defect localization error and the mean localization error are defined by \[LE_j=\max(0,\|\hat{x}_j-x_j^{\dagger}\|-0.05),\qquad MLE=\frac{1}{M}\sum_{j=1}^{M}LE_j .\]

The Contrast-to-Noise Ratio (CNR) quantifies the sharpness of the composite discrete probing index against background artifacts. Let \(\Gamma_h\) denote the set of discrete grid points on the known interface. We define \(S_{\text{target}}\) as the subset of grid points on the interface located within a distance of \(0.05\) from the true defect centers \(x^{\dagger}\), i.e., \(S_{\text{target}} = \{ x \in \Gamma_h \mid \|x - x^{\dagger}\| \le 0.05 \}\). The background set is defined as the remainder of the interface, \(S_{\text{bg}} = \Gamma_h \setminus S_{\text{target}}\). We then define \[CNR=\frac{| \bar{I}_{\text{target}}-\bar{I}_{\text{bg}} |}{\sigma_{\text{bg}}}, \quad \text{where}\quad \bar{I}_{\text{target}}=\frac{1}{|S_{\text{target}}|} \sum_{x \in S_{\text{target}}}I_{\text{comp}}(x),\quad \bar{I}_{\text{bg}}=\frac{1}{|S_{\text{bg}}|} \sum_{x \in S_{\text{bg}}}I_{\text{comp}}(x).\] Here \(\sigma_{\text{bg}}\) is the sample standard deviation of the discrete background index values: \[\sigma_{\text{bg}} =\sqrt{\frac{1}{|S_{\text{bg}}|-1} \sum_{x \in S_{\text{bg}}}(I_{\text{comp}}(x)-\bar{I}_{\text{bg}})^2}.\]

Table 2 summarized the quantitative metrics for Problem I under various configurations. The results confirmed that the \(H^1\) probing index provides highly precise localizations, with MLE values achieving exactly zero in most full-data scenarios. Even under the severe condition of 25% missing boundary data (Examples 5-7), the spatial deviation remained minimal (e.g., bounded by \(0.0167\), which is the size of one mesh grid), while maintaining a sufficiently high CNR to effectively distinguish the defects.

5.2 Problem II↩︎

For Problem II with full data available on the outer boundary, we present four numerical examples to demonstrate the capability of our proposed DSM for reconstructing the interface locations for four scenarios: (i) Three closed interfaces; (ii) Two closed interfaces; (iii) Two non-closed interfaces with relatively small size; (iv) One non-closed interface with relatively large size. For the partial data case (\(75\%\) available), we investigated scenarios (i)(iii)(iv). We can see from Figure 7 that the \(H^1\) probing index function accurately captured the exact interface locations for scenarios (ii)(iii)(iv) with various boundary flux \(g\). For the more complex scenario with three interfaces, the figure similarly demonstrated that for each boundary flux \(g\), the \(H^1\) probing index captured at least two interfaces. Combining these results yielded the desirable full identification of all jump locations.

To quantitatively evaluate Problem II, we adapted the CNR and MLE from Problem I. We utilized the same composite evaluation index \(I_{\text{comp}}(x) = \sum_{i=1}^3 |\widetilde{I}^{(i)}(x)|^2\) (or \(I_{\text{comp}}(x) = \sum_{i=1}^3 |\hat{I}^{(i)}(x)|^2\) for the refined case), which fuses the indices from the three distinct pairs of Cauchy data. We then extended the distance metric to 2D shapes by replacing the point-to-point distance \(\|x - x_j^*\|\) with the shortest Euclidean distance to the true interface, \(d(x, \Gamma^*) = \inf_{y \in \Gamma^*} \|x - y\|\). The CNR target region becomes \(S_{\text{target}} = \{ x \in \Omega_h \mid d(x, \Gamma^*) \le \Delta x \}\), where the mesh size \(\Delta x\) captures the numerical ridge thickness, leaving the remainder of the domain as the background \(S_{\text{bg}} = \Omega_h \setminus S_{\text{target}}\).

For closed rectangular interfaces (Examples 1 and 2), the localization error for the \(j\)-th interface is \(LE_j = d(\hat{x}_j, \Gamma_j^*)\), where the local peak \(\hat{x}_j\) is identified within a \(0.3\)-neighborhood of the rectangle’s geometric center \(x_j^c\), i.e., \(\hat{x}_j = \arg\max_{\substack{x \in \Omega_h \\ \|x - x_j^c\| \le 0.3}} I_{\text{comp}}(x)\). For open line segments (Examples 3 and 4), a single global peak cannot adequately capture trajectory distortions. Instead, we employed a cross-sectional ridge evaluation: the true line is uniformly divided into \(10\) orthogonal bins. Within each bin, we searched for the local peak strictly within a perpendicular distance of \(0.25\) from the true line. The Mean LE was then calculated by averaging the perpendicular distances of these local peaks to the true line, ensuring a rigorous assessment of the overall geometric alignment.

Table 3 summarized the quantitative metrics for Problem II, corresponding to the seven configurations illustrated in Figure 7. The results confirmed that the \(H^1\) probing index provides highly precise localizations for unknown interfaces. Under the full data scenario with \(4\%\) noise (Examples 1-4), the spatial deviation remained exceptionally small, with several sub-interfaces achieving perfectly zero LE or being bounded tightly within \(O(\Delta x)\) (where the mesh size \(\Delta x \approx 0.0167\)). Even under the severe condition of \(25\%\) missing boundary data (Examples 5-7), the algorithm exhibited strong robustness; the maximum localization error slightly shifted but remained highly constrained, while maintaining a robust CNR (consistently above 2.0) across all scenarios to effectively distinguish the interfaces.

The numerical results in this paper demonstrated that the DSM delivers robust performance across a variety of problem settings, highlighting its strong versatility. Moreover, its computational cost was significantly lower than that of the iterative methods discussed in the literature. The experiments further indicated that the DSM can generate accurate initial guesses for iterative schemes, effectively addressing a major challenge faced by most iterative methods. These findings underscored the considerable potential of combining DSM with iterative approaches to enhance computational efficiency and accuracy.

5.3 The Refinement approach↩︎

We applied the same procedure, solving 1 with the Neumann boundary condition \(\partial u/\partial n=g\) on \(\partial\Omega\), to obtain the noisy boundary measurement \(f^{\epsilon}\) on \(\partial\Omega\). We used a uniform mesh of size \(181\times 181\) in the finite difference method to compute the true solution \(u\). To avoid the inverse crime, we then utilized a different uniform mesh of size \(121\times 121\) to execute the iterative algorithm presented in Section 4, sequentially solving for \(u_{1}\), \(\lambda^{1}\), and \(\lambda^{2}\). The noise level was set to \(15\%\) for these tests. The numerical results demonstrate that the refined probing index \(\hat{I}_{0}\) defined in (28 ) outperforms the original probing index \(I_{0}\). This refinement approach was applied to both Problem I and Problem II. For each example, three pairs of Cauchy data \((f, g)\) were provided, yielding three corresponding pairs of normalized probing index functions \(I_0\) and \(\hat{I}_0\). The graphs of \(I_0\) and \(\hat{I}_0\) are displayed in the first and second columns, respectively. The specific settings for these examples are detailed below.
 
Problem I: Known interface. We set \(\delta = 0.1\), \(\alpha=10\).
Example 1: The same setting as in Example 2, Section 5.1.
Example 2: \(\Gamma = \partial\Omega_{1}\), where \(\Omega_{1} = [-0.5, 0.5]\times[-0.5, 0.5]\). \[\gamma(x,y) = \left\{ \begin{align} 100,& \quad y = 0.5, |x+0.1|\leq 0.05, \\ 100,& \quad x = -0.5, |y+0.25|\leq 0.05,\\ 100,& \quad x = 0.5, |y+0.25|\leq 0.05,\\ 1,& \quad\text{else.} \end{align} \right.\] \[g_1(x,y) = 2, g_2(x,y) = 3\cos(x)-1, g_3(x,y) = |\cos(\pi(x+y))|.\] Example 3: \(\Gamma = \partial\Omega_{1}\), where \(\Omega_{1} = [-0.5, 0.5]\times[-0.5, 0.5]\). \[\gamma(x,y) = \left\{ \begin{align} 100,& \quad y = -0.5, |x+0.1|\leq 0.05, \\ 100,& \quad x = -0.5, |y-0.2|\leq 0.05,\\ 100,& \quad x = 0.5, |y-0.3|\leq 0.05,\\ 1,& \quad\text{else.} \end{align} \right.\] \[g_1(x,y) = 2, g_2(x,y) = 3\cos(x)-1, g_3(x,y) = |\cos(\pi(x+y))|.\] Problem II: Unknown interface.
Example 1: The same interface \(\Gamma\) as in Example 2, Section 5.2, \(\gamma = 5\) on \(\Gamma\). We set \(\delta = 0.1\), \(\alpha=10\). \[g_1(x,y) = 2,\quad g_2(x,y) = 3\cos(x)-1,\quad g_3(x,y) = |\cos(\pi(x+y))|.\] Example 2: The same interface \(\Gamma\) and \(\gamma\) as in Example 4, Section 5.2. We set \(\delta = 1\), \(\alpha=1\). \[g_1(x,y) = 2,\quad g_2(x,y) = |\sin(\pi(x+y))|,\quad g_3(x,y) = |\cos(\pi(x+y))|.\] Example 3: \(\Gamma = \partial\Omega_{1}\cup\partial\Omega_{2}\), where \(\Omega_{1} = [0.5, 0.75]\times[-0.75, -0.5]\), \(\Omega_{2} = [0.5, 0.75]\times[0.5, 0.75]\), \(\gamma=3\) on \(\Gamma\). We set \(\delta = 0.1\), \(\alpha=10\). \[g_1(x,y) = 2,\quad g_2(x,y) = 3\cos(x)-1,\quad g_3(x,y) = |\cos(\pi(x+y))|.\]

For Problem I with full data available on the outer boundary, we present 3 numerical examples to demonstrate the efficacy of our refinement approach. Figure 8 shows that in Example 1, both the original \(H^0\) index \(I_0\) and the refined index \(\hat{I}_0\) successfully identified the exact jump locations (marked in red). However, in the next two examples, while \(I_0\) produced a spurious jump that does not exist, the refined \(\hat{I}_0\) remained accurate and correctly detected all true jump locations.

The quantitative improvements are detailed in Table 4. The refinement approach not only consistently enhanced the CNR across all examples but also significantly reduced the MLE when dealing with complex multi-jump scenarios.

For Problem II with full boundary data, we present three numerical examples to validate the effectiveness of our refinement approach, covering scenarios involving multiple closed interfaces and a relatively large non-closed interface. As shown in Figure 9, the original \(H^0\) probing index \(I_0\) failed to distinguish between different interfaces in the case of multiple small closed interfaces. In contrast, our refined \(\hat{I}_0\) significantly improved performance by clearly resolving individual interfaces. Moreover, for the large non-closed interface case, \(\hat{I}_0\) achieved substantially more accurate reconstruction results.

These visual improvements were quantitatively detailed in Table 5. The refinement approach consistently enhanced the CNR across all examples, indicating a sharper index profile and stronger suppression of background artifacts. Notably, in Example 3, the refined index \(\hat{I}_0\) achieved perfect localization with the MLE dropping to 0.0000, while in Example 2, the MLE was significantly reduced from 0.0167 to 0.0033. Even in cases where the overall MLE remained numerically similar (such as Example 1), the elevated CNR confirmed that the refined method yields a more robust and discernible reconstruction against high noise levels.

The most significant advantage of our refinement approach lies in its ability to dramatically enhance the noise robustness of the DSM for both Problem I and Problem II, while maintaining the same level of accuracy.

6 Conclusion↩︎

In this work, we developed a novel Direct Sampling Method (DSM) to address two distinct yet highly ill-posed inverse interface problems for the Laplace equation. The key innovations of this paper lie in two aspects: (i) an optimized boundary condition selection strategy for the reference system, and (ii) a newly proposed refinement approach. Specifically, our numerical evaluations revealed that the \(H^{0}\) probing index function is significantly more robust to noise than the \(H^{1}\) index function, albeit at the cost of reduced localization accuracy. This crucial observation motivated the development of our iterative refinement scheme, which successfully enhances the accuracy of the \(H^{0}\) index while preserving its inherent robustness. Theoretical analysis and numerical experiments demonstrate that this approach maintains the fundamental advantages of the DSM, namely, direct implementation, computational efficiency, and high precision. Remarkably, our method achieved reliable reconstructions using very limited pairs of Cauchy data, thereby bypassing the need for extensive boundary measurements.

This framework not only resolves the targeted problems but also provides a methodological reference for addressing broader classes of complex inverse interface problems. Beyond its standalone efficacy, this framework serves as a valuable pre-processing tool for iterative inversion algorithms. By providing high-quality initial guesses for unknown interfaces or flux jump distributions, our DSM-based approach can significantly accelerate convergence in gradient-based optimization methods. Looking ahead, the integration of our refinement-enhanced DSM with deep learning techniques holds promising potential for tackling more complex PDE-based inverse problems. Such hybrid methodologies could substantially reduce computational costs in real-world applications where data scarcity and noise resilience remain critical challenges.

7 Appendix: Proofs of auxiliary results↩︎

7.1 Proof of Lemma 1↩︎

7.2 Verification of the almost orthogonal property↩︎

The Fourier representation of the probing kernel has been derived in 21 . Here we give the expanded algebraic verification of the radial monotonicity estimates used to justify the almost orthogonal property for \(\alpha=1\). The proof uses the standard series for \(\log(1-t)\) and \(\operatorname{Li}_2(t)=\sum_{n=1}^{\infty}t^n/n^2\), and Sturm’s theorem is used to confirm the signs of the resulting polynomials; see, for example, [25] and [26].

When \(\theta_x = \theta_z\), the sign-normalized kernel used in the localization analysis is given by \[K(x,z) = d\frac{\sum_{n=1}^{\infty}\frac{n^{2}}{n+\alpha}r_{x}^{n}r_{z}^{n}}{\left(\frac{r_{x}^{2}(1+r_{x}^{2})}{(1-r_{x}^{2})^{3}}\right)^{1/4}\left(\sum_{n\in\mathbb{Z}}\frac{n^{2}}{(|n|+\alpha)^{2}}r_{x}^{2|n|}\right)^{1/4}}.\] where \(d>0\). This sign normalization does not affect peak locations, and \(K(x,z)^4\) is proportional to the following function.

References↩︎

---

1. Department of Mathematics and Center for Research in Scientific Computation, North Carolina State University, Raleigh, North Carolina (kito@ncsu.edu).↩︎

2. Department of Mathematics, The Chinese University of Hong Kong, Shatin, N.T., Hong Kong (twu@math.cuhk.edu.hk).↩︎

3. Department of Mathematics, The Chinese University of Hong Kong, Shatin, N.T., Hong Kong. The work of J. Zou was substantially supported by the Hong Kong RGC General Research Fund (projects 14310324 and 14306623) and NSFC / Hong Kong RGC Joint Research Scheme 2022/23 (project N CUHK465/22). (zou@math.cuhk.edu.hk).↩︎