Stable Identiﬁcation of Sources Located on Interface of Nonhomogeneous Media

: This paper presents a stable method for the identiﬁcation of sources located on the separation interface of two homogeneous media (where one of them is contained by the other one), from measurement yielded by those sources on the exterior boundary of the media. This is an ill-posed problem because numerical instability is presented, i.e., minimal errors in the measurement can result in signiﬁcant changes in the solution. To obtain the proposed stable method the identiﬁcation problem is categorized into three subproblems, two of which present numerical instability and regularization methods must be applied to obtain their solution in a stable form. To manage the numerical instability due to the ill-posedness of these subproblems, the Tikhonov regularization and sequential smoothing methods are used. We illustrate this methodology in a circular and irregular region to demonstrate the feasibility of the proposed method, which yields convergent and stable solutions for input data with and without noise.


Introduction
Source identification problems are widely investigated in many research fields [1][2][3][4][5][6][7][8][9], and they are modeled by boundary value problems for which the analysis of the associated forward problem and its corresponding inverse problem must be considered, see, e.g., [10,11]. The latter problem involves determining the source that yields the measurement on the boundary of the region [12][13][14]. The inverse source problem appears in different applications, such as inverse electroencephalography, inverse electrocardiography and inverse geophysics (see, e.g., [15][16][17][18][19][20][21][22][23]), where the problems are modeled using differential equations. Numerical solutions to differential equations are crucial in mathematics and engineering because they appear in many applications, such as population growth, diffusion processes, electromagnetic problems, and elasticity problems. The finite difference method, the finite element method, the boundary element method, and the method of fundamental solutions are used for the numerical solution to boundary value problems, with or without initial conditions on the boundary, governed by certain partial differential equations [24][25][26][27]. The finite element method is one of the most important numerical methods for solving differential equations [28]. It has been used for solving the inverse source problem using the gradient conjugate method and a control approach [19,29,30]. Recently, it has been used to solve differential equations with the Henstock-Kurzweil functions [31]. It is applied to obtain a solution with highly oscillatory functions. Special quadratures have been applied to calculate the integrals that appear in the weak formulation of the finite element method [32]. In this study, we applied the trapezoidal rule to calculate integrals that appear in the applications of the finite element method. The layer potential technique has been used to obtain the solution to the inverse source problem in electrocardiography and electroencephalography (see [15,17,18]), and it was not considered in this study.
In this study, we investigate a case in which the sources are located on the separation interface between two different homogeneous media. The problem of identifying either sources having compact support within a finite number of small subdomains (see [12]), pointwise sources, electrostatic mono and dipole sources, or combinations of these (see [13,14]), are not considered in this work. The sources are identified from measurements on the exterior boundary of the bounded domain, and the relationship between the sources and the measurements is established using an elliptic model, as shown in Section 2. The inverse source problem associated with this elliptic model can be categorized into three problems: The Cauchy problem for the Laplace equation in the homogeneous annular region.

2.
The Dirichlet problem for the Laplace equation in the homogeneous interior region. 3.
The source identification problem from the normal derivatives of the solutions to the Cauchy and Dirichlet problems.
The Cauchy problem is ill-posed because it presents numerical instability, i.e., minimal errors in the measurement can yield significant changes in the solution. To manage such instability, we employed the algorithm proposed in [29], where a penalization method was applied (equivalent to the Tikhonov regularization [33,34]); the conjugate gradient method; and the finite element method, where the regularization parameter was determined by the Tikhonov criterion. To obtain the numerical solution to the forward problem and elliptical problems that appear in the conjugate gradient algorithm and the Dirichlet problem, we performed a finite element approximation [35].
The source was recovered from the normal derivatives of the solutions to the Cauchy and Dirichlet problems. For the stable recovery of the normal derivatives, we applied the sequential smoothing method.
Furthermore, we considered cases involving sources that are typically associated with applications (see [15,17]). This paper is organized as follows. In Section 2, we introduce an elliptic model (a boundary value problem) that relates the sources with the measurements, and the results of existence and uniqueness are included for the classical and weak solutions to the given model. In Section 3, the inverse source problem is uncoupled into three problems that allow us to obtain the stable algorithm. In Section 4, we present the control approach of the Cauchy problem, cost function, and Tikhonov regularization. In Section 5, we present the analytical solution to forward and inverse problems in a circular domain using circular harmonics expansion series. In Section 6, we present the stable algorithm for the inverse source problem. Section 7 presents the numerical results of the proposed method for circular and irregular geometries. In Section 8, we present a discussion about the most relevant results obtained in this work. Finally, conclusions are provided in Section 9.

Elliptic Model
In this study, we focus on the following boundary value problem: Find u 1 and u 2 such that where Ω = Ω 1 ∪ Ω 2 ⊂ R n , with n = 2 or 3, is a sufficiently smooth region. The positive constants σ 1 and σ 2 are the conductivities of regions Ω 1 and Ω 2 , respectively (see Figure 1). Function g represents the source defined on the interface S 1 , n 1 and n 2 are the unit outward normal vectors on the boundaries of Ω 1 and Ω 2 , respectively. We denote u as the solution to problem (1)- (5) in Ω and u i = u| Ω i , i = 1, 2. ∆ represents the Laplace operator, which is also denoted by ∇ 2 . Boundary condition (3) corresponds to the continuity of the potential, and condition (4) corresponds to the jump of the normal flux, which is equal to function g. Condition (5) indicates that the conductivity in the exterior region R n Ω is zero (in some applications, it corresponds to the conductivity of the air). From Green's formulas, the following compatibility condition was obtained: Herein, we refer to problem (1)-(5) as the superficial boundary problem (SBP), which has been investigated in [15,17,30], with applications to the inverse problem of identifying bioelectrical sources on the cerebral cortex from electroencephalographic recording on the scalp. In [20,36], similar problems were investigated with applications to capacitance tomography and sources in the brain.

Definition 1.
For g on S 1 , the forward problem entails obtaining the measurement V = u| S 2 , where u is the solution to the SBP. Definition 2. Given a function V defined on S 2 , the inverse problem entails determining a source g defined on S 1 such that the solution u to the forward problem corresponding to g, satisfies u| S 2 = V.
is a classical solution to problem (1)-(5) if the associated relationships are satisfied in the classical sense of differentiation.
The proof of the following theorem is provided in [30]. Theorem 1. For g ∈ C(S 1 ), compatibility condition (6) is necessary and sufficient for the existence and uniqueness of the classical solution to problem (1)-(5).

Weak Solution to SBP
Let L 2 (Ω i ), L 2 (S i ), i = 1, 2, and L 2 (Ω) be the spaces of square integrable functions defined on Ω i , S i , and Ω, respectively. We denote H 1 (Ω i ), i = 1, 2 and H 1 (Ω) as the corresponding Sobolev spaces of the functions in L 2 (Ω i ) and L 2 (Ω), respectively, with first generalized derivatives, which are square integrable functions. We denote H 1 2 (S i ) as the subspaces of L 2 (S i ), which comprised traces to S i of the functions in H 1 (Ω i ), for i = 1, 2.
The following theorem regarding the existence and uniqueness of the weak solution is detailed in [15]. (6) is necessary and sufficient for the existence of a weak solution to the SBP. A unique weak solution u exists that satisfies Ω u(x)dx = 0 and u H 1 (Ω) ≤ C g L 2 (S 1 ) , where the constant C does not depend on g.

Theorem 2. Condition
It is noteworthy that the compatibility condition (6) is necessary and sufficient for the existence and uniqueness of both classical and weak solutions to the SBP. The condition Ω u(x)dx = 0, guarantees the uniqueness of the solution to the SBP and hence the uniqueness of the forward problem. For the inverse problem, the following theorem regarding the uniqueness of its solution is detailed in [30].
Theorem 3. For a measurement V defined on S 2 , if a source g defined on S 1 exists and satisfies condition (6), which is a solution to the inverse problem, then it is unique.

Operational Statement
We consider the following spaces: We can define the following operator T : U → V, by T(g) = u, where u is the solution to problem (1)- (5). The composition of the continuous operator T and the compact trace operator Tr : V → W, defines the compact operator A : U → W, given by A = Tr • T, i.e., A(g) = u| S 2 , which is injective. This operator is used in [30] to study the problem of identifying sources defined on S 1 .
The inverse operator A −1 is not continuous and it leads to numerical instability, which is one of the causes of the ill-posedness of the problem. To handle the numerical instability a regularization method can be applied. In [30], the Tikhonov regularization method is applied for which the following result was necessary: The equation A(g) = V does not have solution for all V ∈ W, which is another reason of the ill-posedness of this inverse source problem. If some conditions of smoothness are imposed on V, like in [33], it is possible to find global conditions of existence of the solution. This operational statement is used in [30] to analyze the inverse problem, and it is compared at the end of Section 7 with the method proposed in this work.
In [15], the inverse problem (nonlinear) is studied to determine the parameters of a dipolar source. These parameters are the dipolar moment and the center of the source. For that, operator A is used to find the center of the source and the dipolar moment is determined by solving a linear inverse problem. Thus, the solution to the nonlinear inverse problem is obtained by solving two linear inverse problems in which the same regularization parameter is used.

Separation of Inverse Problem into Three Subproblems
To solve the inverse problem of identifying the source from measurement V, we considered the following three problems: 1.
The first one is the Cauchy problem for the Laplace equation defined in the annular homogeneous region Ω 2 as follows: where V = u| S 2 is the known measurement.

2.
The second is the Dirichlet problem for the Laplace equation defined in the homogeneous interior region Ω 1 as follows: with the Dirichlet boundary condition ϕ = u 2 | S 1 , where u 2 is the solution to the Cauchy problem (8).

3.
The third problem entails obtaining g on S 1 by boundary condition (4), i.e., where u 1 and u 2 are the solutions to (8) and (9), respectively. To compute the normal derivative of u 1 and u 2 in the stable form, we applied the sequential smoothing method, as shown in Section 6.

Control Approach of Cauchy Problem: Cost Function
The Cauchy problem is important and applicable to various fields, e.g., inverse electrocardiography, inverse electroencephalography, and electrical capacitance tomography [16,18,21,36]. To investigate the Cauchy problem, we considered the following linear, injective, and compact operator K : H 1/2 (S 1 ) → L 2 (S 2 ) expressed as where w is the solution to the following state equation: The relationship between problems (8) and (12) can be described using operator K as follows: a solution w to problem (12) is also a solution to problem (8) if we select ϕ on S 1 such that where w is the solution to state Equation (12), and V is the known measurement in problem (8), from where ϕ = K −1 (V). The proof of the following theorem, which allows us to apply the regularization methods, is provided in [37].
Because the operator K is linear, injective, compact, and defined on a space of infinite dimension, its inverse K −1 is not continuous [33]. This results in numerical instability, which causes the ill-posedness of this problem. To manage this numerical instability, we used the Tikhonov regularization method in this study, which will be stated below, in Section 4.1.2. In the following subsection we provide the variational formulation which allows us to relate the control approach and the Tikhonov regularization in a straightforward manner. In [29], the Cauchy problem (8) is solved by minimizing the following cost function: where k is the penalization parameter, ϕ is the boundary control on S 1 in (12), and K is defined by (13). Therefore, we can approximate the controllability problem (14) using one of the following kinds of approximation by penalty: The density property provided in Theorem 5 and the standard convexity arguments guarantee that (15) has a unique minimum (see, for instance [38]), characterized by To minimize the cost function (14), we applied the conjugate gradient method provided in [29], where the first derivative (variation) of J k is required. In this case, the derivative of J k is expressed as where v is the solution to the adjoint problem and u(ϕ) is the solution to state Equation (12).

Tikhonov Regularization Method
This method is used to handle the numerical instability associated to ill-posed problems, which depends on a parameter α called the regularization parameter that is chosen in terms of the error δ that is known [33]. Let V δ be the noisy data of measurement V in Cauchy problem (8) To manage the numerical instability from using equation K(ϕ) = V δ , we used the Tikhonov regularization method. The selection of the penalization parameter k in (14) is equivalent to obtaining a regularization parameter α = 1/k for the Tikhonov functional where α(δ), which depends on δ, must be selected appropriately. Hence, we used the known results of the Tikhonov regularization to select the appropriate α(δ).
Let ϕ α(δ) = R α V δ be an approximation of ϕ, where K * is the adjoint operator of K, and R α = (K * K + αI) −1 K * V δ is the operator associated with the Tikhonov regularization strategy that possesses the properties described in the following theorems [33]: . In this case, R α is known as 'admissible', and Different methods can be used to determine the suitable regularization parameter that balances the two terms of the functional J α , e.g., the discrepancy principle and the L-curve. The term L-curve originates from the shape of the curve which corresponds to a log-log plot of the norm of a regularized solution vs. the norm of the corresponding residual norm. More geometrical and analytical details are provided in [39].

Solution to Forward and Inverse Problems in Circular Domain
To obtain the solution to elliptic problems (9) and (12), we can apply different analytical and numerical methods. For example, the Green-type matrix method (which generalizes the Green function method), Fourier series method, method of fundamental solutions, finite element method, and finite difference method [25][26][27][28]. We considered the finite element and Fourier series methods in a two-dimensional case involving a simple geometry to demonstrate the feasibility of the proposed algorithm. All the results can be extended to a three-dimensional case. The forward problem has a unique solution, except for constants. In the functional spaces defined above, in which the solution is sought, constants are disregarded. For the inverse problem, we considered Theorem 3 to guarantee the uniqueness of the solution.
Subsequently, the analytical solution to the forward and inverse problems were obtained using Fourier series (circular harmonics). The results were used to validate the proposed algorithm.

Analytical Solution to Forward and Inverse Problems in Circular Domain
We considered the case of two concentric circles centered at the origin, i.e., Ω 1 is a circular domain of radius R 1 , and Ω 2 is a circular annular domain with interior radius R 1 and exterior radius R 2 . In this case, we applied the Fourier series method to solve the forward problem.
Let g ∈ L 2 (S 1 ) such that it satisfies condition (6). Subsequently, we considered where g 1 k and g 2 k are the Fourier coefficients of the source g, i.e., for k = 1, 2, . . .
We sought the solution to the SBP in the following form: The functions defined in (21) and (22) are harmonics in Ω 1 and Ω 2 , respectively, i.e., these functions satisfy the Laplace equation in polar coordinates. From conditions (3)-(5), we obtained the Fourier coefficients in (21) and (22): Evaluating based on r = R 2 , the solution to the forward problem is expressed as where k are the Fourier coefficients of V, the solution to the inverse problem (in the case of noiseless data) V is expressed as where g 1

Stable Algorithm for Solving Inverse Problem
Considering the subproblems in which the inverse problem was separated, we proposed the following algorithm for identifying the source g on S 1 , based on the measurement of S 2 (with noisy or noiseless data).
Step 1. To obtain the stable form of the numerical solution u n 2,h,k to the Cauchy problem (8) minimizing cost function (14) using the conjugate gradient method and the finite element method provided in [29], where h is the size of the triangular mesh of Ω, k is the penalization parameter of function J k , and n denotes the number of conjugate gradient iterations.
Step 2. To use the finite element method to compute the solution u n 1,h,k to problem (9) in region Ω 1 with Dirichlet boundary condition ϕ n 2,h,k = u n 2,h,k | S 1 , where u n 2,h,k is the solution to the Cauchy problem (8) provided in Step 1.
Step 3. To compute the source g h,k,K on boundary S 1 as follows: • First, in each i-th node p i (of the triangular mesh of Ω) on boundary S 1 , we compute the normal derivative of u n 2,h,k as follows: where u n 2,h,k is the solution to Cauchy problem (8) provided in Step 1 (in this case, S 1 is the interior boundary of Ω 2 ). • Second, in each i-th node p i (of the triangular mesh of Ω) on boundary S 1 , we compute the normal derivative of u n 1,h,k as follows: where u n 1,h,k is the solution to the Dirichlet problem (9) provided in Step 2 (in this case, S 1 is the boundary of Ω 1 ). Here, {ϕ i } nnt κ i=1 are functions such that ϕ i (p j ) = δ ij (δ ij is the Kronecker delta), which forms a basis for the discrete space , respectively, where K is the number of sequential smoothing.

Discretization of Elliptic Problems Using Finite Element Method
The computational implementation of the conjugate gradient algorithm requires solutions to problems (12) and (17) at the initialization step, solutions to two other elliptic problems at each iteration (see [29]), and solutions to problem (9). The numerical solutions to these problems can be achieved using any approximation method. In this study, we advocate the use of finite element approximation because it is suitable for managing the variational formulation of the inverse problem. Subsequently, we consider the following elliptic model problem: where Ω is a bounded region in R n (n =1, 2, or 3) with a sufficiently smooth boundary ∂Ω, assuming that ∂Ω comprises two non-intersecting (n − 1)-dimensional boundaries S 1 and S 2 . For problem (12) or (17), σ = σ 2 ; Ω = Ω 2 ; S 1 and S 2 are the interior and exterior boundaries of Ω 2 , respectively. For problem (9), σ = σ 1 , Ω = Ω 1 , S 1 = ∂Ω 1 , and S 2 is the empty set ∅. The variational formulation of problem (28) is as follows: Find where To approximate V ϕ , we considered a finite element triangulation τ h of Ω, where h is the length of the largest edge of τ h , i.e., region Ω is approximated by a polyhedral region Ω h = ∪ T∈τ h T. Space V ϕ is approximated by a discrete space as follows: ∀T ∈ τ h , where C(Ω h ) denotes the space of the continuous functions defined on Ω h . Similarly, V 0,h is the discrete space of V 0 . Subsequently, the variational problem can be approximated using the following discrete variational formulation: where ψ h is a piecewise linear approximation of ψ on the exterior boundary S 2,h of Ω h . The integral (31) satisfies the property Ω h = ∑ T∈T h T , which can be calculated using the trapezoidal rule. The various discrete linear elliptic problems presented in (9), (12) and (17) are associated with the same stiffness matrix, differing only by their right-hand sides. Because the stiffness matrix is symmetrical, positive definite, and sparse, the associated linear systems can be solved using a sparse Cholesky solver, similar to that available in MATLAB. Therefore, we can obtain the solution u h in the nodes of the triangulation τ h (see [35]).

Sequential Smoothing Method
Let f = ( f 1 , f 2 , . . . , f N ) T be a column vector from discrete data of length N, with , and δ i denotes errors. The sequential smoothing method in iterative form using cubic splines can be written as for k = 1, 2, . . . , K, where s i−2,4 is a local basis cubic spline (B-spline), and K > 0 is an integer. The number of sequential smoothing K is selected based on the error in the input data and the discretization parameter N. Hence, f K = S K (x, f ), where K can be selected based on the residual (discrepancy) principle (see [40]). For a uniform grid {x i } N i=1 , the sequential smoothing method is expressed as where M is an N × N matrix, expressed as (for more details, see [40,41]) In our case, because the interior boundary S 1 has a uniform grid, we applied the following formulas: ∂u n 2,h,k,K and ∂u n 1,h,K where M is the N × N matrix provided in Equation (33), N the number of nodes on S 1 , and K the number of sequential smoothing.
To ensure periodicity, we averaged the values of the normal derivative with five points around the endpoints and took this average as the value at the endpoints. This process is repeated for each smoothing step of the normal derivatives.
In this form, we obtained the smoothing of the normal derivatives  (34) and (35), respectively, in Step 3 of the algorithm proposed in Section 6.

Numerical Examples
Based on two examples, we present the numerical results obtained using the methodology discussed in Section 6. Each example includes a different type of two-dimensional bounded region: a circular Ω determined by two concentric circles centered in the origin, i.e., Ω 1 is a circular domain of radius R 1 ; Ω 2 , which is a circular annular domain with interior radius R 1 and exterior radius R 2 (see Figure 1a), and the complex region shown in Figure 1b. The forward problem is solved by selecting a source g on S 1 , followed by obtaining the solution u to the SBP (1)-(5) and then computing V = u| S 2 . Subsequently, the obtained potential V is provided as input data to the control problem (15) or cost function (14). Next, we consider examples with noiseless data V and noisy data V δ , where δ is the measurement error, i.e., V δ − V L 2 (S 2 ) ≤ δ. The noisy data V δ were obtained by adding a random error using the rand function in MATLAB. Therefore, we define where Error = δ max |V| [2 rand(1, s) − 1] is a vector of random numbers of length s, where s is the numbers of nodes on S 2 .
Since the measurements are obtained by a device, they contain inherent error generated by different circumstances, such as rounding error (due to approximation or truncation), measurement handling, and device quality. We want to emphasize that the algorithm pro-posed works with different error levels. As an example, the electroencephalogram recorded on the scalp detects the voltage produced by bioelectrical activity and also detects artifacts produced by ocular and muscular movement. Thus, filters are applied to obtain a clean signal. Additionally, the sensitivity of the device generates errors.
However, we want to emphasize that the algorithm proposed works with different error levels. According to the numerical results, the error between the exact and approximate sources obtained by this method is proportional to the measurement error.
The corresponding approximate solution to problem (9) is denoted by u n 1,h,k , and the corresponding approximate solution to state Equation (12) is denoted by u n 2,h,k and V n h,k = u n 2,h,k | S 2 , where n denotes the number of conjugate gradient iterations to achieve a specific tolerance ε, h is the size of the triangular mesh of Ω, and k is the penalization parameter. Moreover, their corresponding approximations to the normal derivative without and with sequential smoothing are denoted by , and ∂u 2,h,k,K ∂n 1 S 1 , respectively, where K is the number of sequential smoothing. Finally, we denote g n h,k,K as the numerical approximation to g from the exact measurement V.
The corresponding approximations to the normal derivative without and with the sequential smoothing of u 1,h,k(δ) and u n 2,h,k(δ) are denoted by and ∂u 2,h,k(α),K ∂n 1 S 1 , respectively. Hence, we recovered an approximate stable source g n h,k(δ),K from the measurement with error V δ . For the numerical implementation of the proposed method in each example, we applied the finite element method to solve the elliptic problems that appeared in the algorithm proposed herein, using three different meshes. These meshes were denoted by M i (nn i , ne i ), i = 1, 2, 3, where nn i and ne i represent the number of nodes (vertices) and elements (triangles), respectively. Meshes M 2 and M 3 were obtained from successive regular refinements of the coarsest mesh M 1 with mesh size h = 0.1. Figure 1a,b show the corresponding mesh M 1 of the circular and irregular region, respectively. Meshes M 1 , M 2 , and M 3 were generated using pdetool in MATLAB.
In the examples below, we considered the following absolute and relative errors for different values of δ: AE(g n h,k(δ),K , g) = g n h,k(δ),K − g L 2 (S 1 ) , RE(g n h,k(δ),K , g) = g n h,k(δ),K − g L 2 (S 1 ) / g L 2 (S 1 ) , where g n h,k(δ),K is the recovered source, and g is the exact source. Furthermore, we considered the relative error between the noisy data V δ and noiseless data V, expressed as Example 1. In this example, we considered the circular region Ω defined above, with σ 1 = 3, σ 2 = 1, R 1 = 1, and R 2 = 1.2. Let g(x, y) = x 2 − y 2 for all (x, y) ∈ S 1 be the exact source; therefore, it belongs to L 2 (S 1 ) and satisfies condition (6). In polar coordinates, g is expressed as where the unique Fourier coefficient different from zero is g 1 2 = R 2 1 . Subsequently, the exact solution u to the SBP for the exact source g, using circular harmonics series, is and where their corresponding exact normal derivatives on the interior boundary S 1 are expressed as The solution to the forward problem for the exact source g is the exact measurement V provided in Equation (23), i.e., Results for circular region with noiseless input data V. For the case of noiseless data V, the numerical results are summarized in Table 1, where the following errors were included: AE(g n h,k,K , g), and RE(g n h,k,K , g). In this case, the penalization parameter, associated to the Cauchy problem, was set to k = 10 20 , and we defined the stopping tolerance (which is the criterion used to stop the conjugate gradient method when the residual norm g n+1 2 L 2 (S 1 ) ≤ ε max{1, g 0 2 L 2 (S 1 ) } for n = 0, 1, 2, . . . (see [29])) as ε = 10 −9 . The relative errors decreased with each refinement of mesh M 1 , showing the convergence of the discrete space. Moreover, the estimated computation times of the algorithm proposed for meshes M 1 , M 2 , and M 3 were 1.7715, 4.1419, and 30.2689 s, respectively. The computation time of mesh M 2 was greater than that of mesh M 1 because the numbers of nodes (vertices) and elements (triangles) were greater than those in mesh M 1 ; furthermore, the number of smoothing increased. The computation time of mesh M 3 was the highest because the numbers of nodes (vertices), elements (triangles), and smoothing increased.  Figure 2 shows the exact source g and recovered source g n h,k,K (at iteration n = 26) on S 1 , for noiseless data V. The approximation is acceptable from the practical perspective. In fact, the highest relative error was 4.0355 × 10 −2 , which was obtained using the coarsest mesh M 1 , where g n h,k,K was calculated as shown in Step 3 of Section 6, i.e., by computing the approximate solutions u n 1,h,k and u n 2,h,k to problems (8) and (9) after applying the sequential smoothing method for computing their normal derivatives in stable form, as denoted by ∂u 1,k,h,K ∂n 1 and ∂u 2,k,h,K ∂n 1 on S 1 , respectively. These are shown in Figure 3, where the number of sequential smoothing was K = 4. A similar behavior was observed when using meshes M 2 and M 3 .  Results for circular region with noisy input data V δ . We selected the noisy data V δ from Equation (36) and the coarsest mesh M 1 . The values differed from the noiseless data V, as RE(V δ , V) defined above. Table 2 shows the results for the circular region for different values of δ, from where we observed that the numerical solutions g n h,k(δ),K converged to the noiseless solutions g n h,k,K as δ approached zero; this demonstrates the stability of the numerical results against perturbations in measurement V for the circular region. In addition, the computation time of the algorithm proposed for each value of δ was 0.2117, 0.2409, 0.2418, and 1.2474 s, with K = 8. The computation time for δ = 0 was the lowest because the number of smoothing was minor (K = 4).  Figure 4 shows the plots of the exact measurement V, recovered data V n h,k(δ) with measurement with error V δ when δ = 0.1, and corresponding recovered source g n h,k(δ),K (at iteration n = 1) on S 1 that was calculated as shown in Step 3 of Section 6. Figure 5 shows the smoothing normal derivatives ∂ ∂n 1 u n 1,h,k(δ),K and ∂ ∂n 1 u n 2,h,k(δ),K from the corresponding recovered derivatives of u n 1,h,k(δ) and u n 2,h,k(δ) , respectively. In this case, the penalization parameter k(δ) was selected as in the case of noisy data V δ , i.e., k(δ) = δ −1 . The stopping tolerance was ε = 10 −3 , and the number of smoothing was K = 8 for δ = 0.001, 0.01, and 0.1. In this case, the highest relative error between the exact source g and the recovered source g n h,k(δ),K was 8.3461 × 10 −2 using the coarsest mesh M 1 . We observed that the recovered solutions g n h,k(δ),K converged to the exact solution g when δ approached zero. (blue), for δ = 0.1. (b) Plots of exact source g (blue) and recovered source g n h,k(δ),K (green) on S 1 , corresponding to data with noise V δ , using mesh M 1 of circular region Ω. Figure 5. (a) Plots of exact normal derivative of u 1 and its approximations with and without sequential smoothing corresponding to input data with error V δ , using mesh M 1 of circular region. (b) Plots of exact normal derivative of u 2 and its approximations with and without sequential smoothing corresponding to input data with error V δ , using mesh M 1 of circular region.
Furthermore, the exact solution u to the SBP and the exact measurement V = u| S 2 (for this source g) were calculated numerically using the finite element method and a very fine mesh, which we denote as M 4 . This mesh was obtained after three regular refinements of the starting mesh, i.e., M 1 shown in Figure 1b. In addition, we computed the exact normal derivatives of u 1 and u 2 on the interior boundary S 1 , as shown in Equations (25) and (26), respectively, using the exact solution u calculated based on the finest mesh M 4 .
Results for irregular region with noiseless data V. The numerical results in Table 3 show the convergence regarding the discrete space of the complex domain Ω for the source g defined in Equation (44), with the penalization parameter k = 10 20 and stopping tolerance ε = 10 −9 .
The computation time in mesh M 3 was the highest because mesh M 3 was the finest, i.e., the numbers of nodes (vertices) and elements (triangles) were greater than those in meshes M 1 and M 2 , and the number of smoothing increased to 30 in mesh M 3 . We observed that the relative numerical errors on the interior boundary, RE(g n h,k,K , g), were slightly higher for the complex region because the interior boundary S 1 was more complex than the corresponding circular region. In addition, the number of conjugate gradient iterations for achieving the desired accuracy was higher for the complex region presented in meshes M 2 and M 3 . Table 3. Convergence of numerical results for irregular region with noiseless data V, where k = 10 20 and ε = 10 −9 .  Figure 6 shows the exact g and recovered sources g n h,k,K (at iteration n = 24) on S 1 , for noiseless data V. In this case, the highest relative error of 4.8550 × 10 −2 was observed in the coarsest mesh M 1 . A similar behavior was observed when using meshes M 2 and M 3 . The plots of the exact normal derivatives of u 1 and u 2 and their corresponding approximations without and with sequential smoothing corresponding to the input data without error V are shown in Figure 7, where the number of sequential smoothing was K = 4.

Mesh
Results for irregular region with noisy input data V δ . Table 4 shows the numerical results for the complex region when the input data V were perturbed, as in Equation (36) for source in Equation (44), using mesh M 1 of the complex region Ω shown in Figure 1b. In each case, one can observe that the qualitative behavior was similar to that of the circular region, and the highest relative error of the approximate solution g n h,k(δ),K was 1.7011 × 10 −1 in the coarsest mesh M 1 . In addition, in this mesh, the computation time for δ = 0 was lower than for δ = 0.001, 0.01, and 0.1 because the number of smoothing was lesser (K = 4). The numerical solutions g n h,k(δ),K converged to g n h,k,K when δ approached zero, demonstrating the stability of the numerical results against perturbations in the irregular region.    Figure 8 shows the plots of the exact measurement V, recovered data V n h,k(δ) with measurement error V δ when δ = 0.1 (left), and corresponding recovered source g n h,k(δ),K (at iteration n = 2) on S 1 of the irregular region (right), which was calculated using the algorithm proposed in Section 6. Figure 9 shows the normal derivatives without and with the sequential smoothing of u n 1,h,k(δ) and u n 2,h,k(δ) . In this case, the penalization parameter k(δ) was selected, as in the case of noisy data V δ , i.e., k(δ) = δ −1 . The stopping tolerance was ε = 10 −3 , and the number of smoothing was K = 10. In this case, the highest relative error between the exact source g and recovered source g n h,k(δ),K was 1.7011 × 10 −1 in the coarsest mesh M 1 .  Therefore, we demonstrated that our approach yielded convergent solutions that were stable against perturbations of the input data V in circular and irregular nonhomogeneous region media Ω.

Comparison with Other Method
We present the comparison of the solutions to the inverse source problem using the method proposed in this work that we called Method 1 (see Section 6) and the one presented in reference [30] that we called Method 2, which uses operator A shown in Section 2.3 (Operational statement).
We consider the circular and irregular region given in the examples from this section with the same values for R 1 = 1, R 2 = 1.2, σ 1 = 3, and σ 2 = 1. In the case of the circular region Ω, we consider the exact source g(x, y) = x 2 − y 2 , for all (x, y) ∈ S 1 , that in polar coordinates is given by (38), and the exact source g given by (44) with β 2 = 0.1 and (a 1 , a 2 ) = (−1, 0) ∈ S 1 , that expressed in polar coordinates takes the form where , for all θ ∈ [0, 2π], θ 0 = π, and β 2 = 0.1. In this case, g belongs to L 2 (S 1 ) and satisfies condition (6). For this exact source g, we approximated function g ∈ L 2 (S 1 ) by its truncated Fourier series g N using the first N terms. In this case, both the exact solution u to the forward problem (1)-(5) and the exact measurement are generated with the first N = 30 terms of the Fourier series (21)-(23), and g − g N L 2 (S 1 ) = 1.4011 × 10 −8 . The Fourier coefficients g 1 k , g 2 k , k = 1, 2, . . . , N are obtained numerically using the function quadl of MATLAB. The plots of g and its approximation g N are shown in Figure 10. In the following, we denote the approximation g N by g. In the case of the irregular region Ω shown in Figure 1(b), we considered the exact g Figure 10. Plots of exact source g (blue) and its approximation g N by its truncated Fourier series (magenta) for N = 30 terms.
In the case of the irregular region Ω shown in Figure 1b, we considered the exact source g given by (44), with the same values for β 2 = 0.1 and (a 1 , a 2 ) = (−1.0067, 0.0093) ∈ S 1 , and the exact source g expressed as for all (x, y) ∈ S 1 , where f 2 (x, y) = x 2 − y 2 . In this case, g belongs to L 2 (S 1 ) and satisfies condition (6) too. Furthermore, the exact solution u to the forward problem (1)-(5) and the exact measurement V = u| S 2 (for this source g) were calculated numerically using the finite element method with the finest mesh M 4 of the irregular region Ω. In addition, we computed the exact normal derivatives of u 1 and u 2 on the interior boundary S 1 , as shown in Equations (25) and (26), respectively, using the exact solution u calculated based on the finest mesh M 4 . Tables 5 and 6 show the numerical results for the circular and irregular region for input noisy data V δ , for δ = 0.1 in mesh M 1 of each region, respectively, and where the computation time and the following relative errors are included: RE(g n h,k(δ),K , g), RE(g n h,α(δ) , g) of each method, where g n h,α(δ) is the corresponding recovered source given by Method 2. In this case, the regularization parameter α(δ) was chosen by the L-curve method.
In order to compare the efficiency of Methods 1 and 2, we considered the same measurement with error V δ of each source g in each region Ω, and the same stopping tolerance ε = 10 −3 for the conjugate gradient method proposed in Methods 1 and 2. In this case, the regularization parameter α(δ) = 10 −3 for the circular region and α(δ) = 10 −4 for the irregular region. From the numerical results, we observed the following important properties (for δ = 0.1): (a) All relative errors corresponding to the sources recovered by each method are proportional to the relative perturbation ER S 2 (V δ , V), showing the stability of the numerical results regarding perturbations for the circular and irregular region, respectively; (b) All relative errors corresponding to the sources recovered by each method are of the same order in each region Ω; (c) The computation time of Method 1 was lower than that of Method 2, because Method 1 solves the Cauchy problem for the Laplace Equation (using the conjugate gradient method provided in [29]), in the smaller sub-region Ω 2 contained in Ω and with a smaller number of iterations than Method 2 that solves the inverse problem applying the conjugate gradient method proposed in [30], wherein each iteration two elliptical problems are solved in the largest region Ω. In addition, Method 1 also includes the computation time to solve the Dirichlet problem for the Laplace equation in the smaller sub-region Ω 1 contained in Ω (which is calculated only once), the calculation of normal derivatives and the applying of the sequential smoothing method for the calculation of the stable derivatives, and the calculation of the recovered source g n h,k(δ),K using formula (10).
The plots of the sources recovered from the measurement with error V δ (δ = 0.1) by each method are shown in Figures 11 and 12, using mesh M 1 of the circular and irregular region, respectively.
From Tables 5 and 6, we observed that the method proposed in this work is better, regarding the computational cost, than the one presented in [30]. Furthermore, the relative errors of the sources recovered by Method 1 are of the same order as the ones recovered by Method 2, which are proportional to the relative perturbation RE(V δ , V), showing the stability of the numerical results regarding perturbations for both the circular and irregular region. Table 5. Numerical results applying methods 1 and 2 for circular region with noisy data V δ for δ = 0.1, using mesh M 1 (774, 1470), and with the same measurement with error V δ for each source. NA means not applicable.

Discussion
In this work, we presented a stable algorithm for recovering sources located on the interface that separates two homogeneous media (which compound the nonhomogeneous media studied in this work). This algorithm works for simple and complex geometries and consists of three steps given in Section 6. Each problem is solved appropriately. Two of them present numerical instability, namely, the Cauchy problem for the Laplace equation defined in an annular region and the calculation of the normal derivatives, which involves the derivative operator that is discontinuous. We consider the results given in [29], in which an operational statement of the problem with a linear, injective and compact operator defined between Hilbert spaces is used. The ill-posedness due to the numerical instability of the Cauchy problem is proved using this operational statement. A control approach and the Tikhonov regularization are considered to handle the numerical instability. We used the sequential smoothing method as a regularization technique for the stable calculations of the derivatives. In some inverse problems, such as inverse electroencephalography and capacitance tomography it is necessary to obtain the normal derivatives to find the solution to these inverse problems [20,36]. This work shows the feasibility of applying the sequential smoothing method to find the normal derivatives in a stable form. The sequential smoothing method is simple to program and apply. In [30], the authors presented another method to find the source in which it is not necessary to find the normal derivatives that represent a disadvantage regarding the method presented in this work. Furthermore, the computational cost is lower in the method presented here. In this work, we numerically chose the number of sequential smoothing. One interesting point is to find a theoretical result to choose the number of sequential smoothing in terms of the measurement error (taking into account the results for the Cauchy problem). We consider this point in future work. Finally, we want to emphasize that the algorithm presented in this work can be extended to a region Ω in the 3D space.

Conclusions
The problem of identifying sources defined on separation interfaces in nonhomogeneous media can be solved in a stable form using the proposed algorithm, which can be applied to regions with simple and complex geometries. Its efficiency was demonstrated based on two examples involving nonhomogeneous circular and irregular regions. The numerical results demonstrated the convergence and stability of the algorithm regarding the input data with and without noise. The algorithm employs the penalization and sequential smoothing methods, where the penalization parameter can be selected as the inverse of the Tikhonov regularization parameter, whereas the number of sequential smoothing K can be selected based on the residual principle. The sequential smoothing method allows us to find the normal derivatives of the potentials on the separation interface in a stable form.