A NEW KOHN-VOGELIUS TYPE FORMULATION FOR INVERSE SOURCE PROBLEMS

In this paper we propose a Kohn-Vogelius type formulation for an inverse source problem of partial differential equations. The unknown source term is to be determined from both Dirichlet and Neumann boundary conditions. We introduce two different boundary value problems, which depend on two different positive real numbers α and β, and both of them incorporate the Dirichlet and Neumann data into a single Robin boundary condition. This allows noise in both boundary data. By using the Kohn-Vogelius type Tikhonov regularization, data to be fitted is transferred from boundary into the whole domain, making the problem resolution more robust. More importantly, with the formulation proposed here, satisfactory reconstruction could be achieved for rather small regularization parameter through choosing properly the values of α and β. This is a desirable property to have since a smaller regularization parameter implies a more accurate approximation of the regularized problem to the original one. The proposed method is studied theoretically. Two numerical examples are provided to show the usefulness of the proposed method.


(Communicated by Fioralba Cakoni)
Abstract. In this paper we propose a Kohn-Vogelius type formulation for an inverse source problem of partial differential equations. The unknown source term is to be determined from both Dirichlet and Neumann boundary conditions. We introduce two different boundary value problems, which depend on two different positive real numbers α and β, and both of them incorporate the Dirichlet and Neumann data into a single Robin boundary condition. This allows noise in both boundary data. By using the Kohn-Vogelius type Tikhonov regularization, data to be fitted is transferred from boundary into the whole domain, making the problem resolution more robust. More importantly, with the formulation proposed here, satisfactory reconstruction could be achieved for rather small regularization parameter through choosing properly the values of α and β. This is a desirable property to have since a smaller regularization parameter implies a more accurate approximation of the regularized problem to the original one. The proposed method is studied theoretically. Two numerical examples are provided to show the usefulness of the proposed method.

Introduction.
Let Ω ⊂ R d (d ≤ 3: space dimension) be an open bounded set with a Lipschitz boundary Γ, and let Ω 0 ⊂ Ω be a Lipschitz subset. We consider the following inverse source problem: Given g 1 and g 2 on Γ, find p such that the solution of boundary value problem (BVP) (1) −∆u + u = pχ Ω0 in Ω, ∂u ∂n = g 2 on Γ satisfies (2) u = g 1 on Γ.
Here ∂ ∂n is the outward normal derivative, and χ Ω0 is the characteristic function of Ω 0 , i.e., its value is 1 in Ω 0 while 0 outside Ω 0 . In the context of applications in bioluminescence tomography, Ω 0 is known as a permissible region of the source function p.
It is well known that inverse source problems are under-determined and ill-posed. In this paper, we consider using Tikhonov regularization [9,13] to obtain a stable approximate solution for the inverse source problem associated with the BVP (1). We note that the inverse source problem (1)- (2) has been studied extensively through Tikhonov methods in the field of bioluminescence tomography, see e.g. [7,8,11] and references therein. With the Tikhonov regularization, the original inverse source problem (1)-(2) is converted to the following minimization problem: (3) p ε = arg min p∈Q ad where Q ad ⊂ Q := L 2 (Ω 0 ) is an admissible set, incorporating a priori information about the source function p; for each p ∈ Q, u(p) is the weak solution of the BVP (1) in the Sobolev space V := H 1 (Ω); · 0,Γ and · 0,Ω0 are standard norms of the spaces L 2 (Γ) and L 2 (Ω 0 ). Under some assumptions, Problem (3) admits a unique solution p ε , and as ε → 0, p ε converges to p * , a solution of (1)-(2) with minimal L 2 -norm ( [7]). In [12], a Kohn-Vogelius type functional is used for the data fitting in solving the inverse source problem. In general, Kohn-Vogelius functionals are expected to lead to more robust optimization procedures [1]. For the inverse problem (1)-(2), the formulation studied in [12] is (4) p ε = arg min p∈Q ad where u 1 (p), u 2 (p) ∈ V are the weak solutions of BVPs −∆u 1 + u 1 = pχ Ω0 in Ω, u 1 = g 1 on Γ, and −∆u 2 + u 2 = pχ Ω0 in Ω, ∂u2 ∂n = g 2 on Γ, respectively; · 0,Ω is the standard L 2 (Ω) norm. Under certain assumptions, Problem (4) admits a unique stable solution p ε which converges to p * ( [12]). The regularization parameter ε plays an important role in solution accuracy and stability. In this paper, we propose a new Kohn-Vogelius type functional to solve the inverse source problem which allows the use of very small values of the parameter ε. Consequently, we can compute accurate approximations of the solution corresponding to rather small ε. This is a desirable property since the smaller the regularization parameter, the better the approximation of the regularized problem to the original one. Moreover, with introduction of two positive real parameters, boundary condition data and boundary measurement data can be incorporated into a single boundary value condition, and this improves the stability of the solution with respect to the noise in the data. Moreover, data to be fitted is transferred from boundary into the whole domain, making the problem resolution more robust.
The structure for the rest of this paper is as follows. The new Kohn-Vogelius type method is proposed for the problem (1)-(2) in section 2. Some theoretical properties are shown in this section. We discuss in section 3 the model when the data contains random noise, and a convergence result is proved. In section 4, we provide a direct study of a system that defines the source function from the regularized formulations, with or without the noise. In section 5, an iterative algorithm is stated for practical reconstructions. Its convergence is ensured when the two parameters α and β are sufficiently close. Several numerical examples are presented in section 6 to demonstrate the feasibility and efficiency of the proposed method. Finally, concluding remarks are given in section 7.
The source function is sought within an admissible set Q ad , assumed to be a closed convex subset of Q. We introduce the following regularized problem. Problem 1. Find p ε ∈ Q ad such that Regarding Problem 1, we have the following well-posedness results.
Proposition 1. For ε > 0, Problem 1 admits a unique solution p ε ∈ Q ad , which depends continuously on g 1 , g 2 ∈ H −1/2 (Γ) and ε. The following first order necessary and sufficient condition for p ε holds: where Φ 1 , Φ 2 ∈ V are weak solutions of the adjoint BVPs Proof. It is not difficult to verify that for ε > 0, J ε (·) is strictly convex on Q. Because Q ad is a convex and closed subset of space Q, it is a standard result ( [2,10]) that Problem 1 admits a unique solution p ε ∈ Q ad , depending continuously on g 1 , g 2 ∈ H −1/2 (Γ) and ε > 0, and is characterized by the inequality Next we derive an equivalent formulation for (11). For any q ∈ Q, we consider the auxiliary BVPs: Then for any t ∈ R, the solutions u 1 (p ε + tq), u 2 (p ε + tq) of (5), (6) (p replaced by p ε + tq) equal to u 1 (p ε ) + tw 1 (q), u 2 (p ε ) + tw 2 (q), respectively. Therefore, the Gateaux derivative of J ε at p ε in the direction q ∈ Q is By using integration by parts together with the definitions of w 1 (q), w 2 (q), Φ 1 and Φ 2 , the above expression reduces to Thus, (8) follows from (11). The proof is completed.
Throughout this paper, we assume that there is at least one solution to the original noise-free inverse source problem (1)-(2), and denote by S ⊂ Q the set of all the solutions, which is also the solution set of Problem 1 with ε = 0. It is not difficult to verify that S is closed and convex. Denote by p * the unique element in S with minimal L 2 -norm. Since we have p ε 0,Ω0 ≤ p * 0,Ω0 . Then with arguments similar to those in [7, Proposition 3.5], we have the following limiting property.
The optimality condition (8) is equivalent to a nonlinear equation where P Q ad is a orthogonal projection from Q onto Q ad . On the function 1 ε (αΦ 1 − βΦ 2 ), we have the following result.
Proposition 3 indicates that a reasonable reconstruction could be achieved for rather small regularization parameter with properly selected α and β. More pre- , and this provides a guidance on how to choose α and β in numerical simulations, cf. Section 6.
3. Noise model and solution convergence. In this section, the discussion focuses on the case where the Dirichlet data g 1 and Neumann data g 2 are allowed to contain random noise. Thus, assume noise data g δ 1 and g δ 2 satisfy (23) g δ 1 − g 1 0,Γ ≤ δ and g δ 2 − g 2 0,Γ ≤ δ, where δ > 0 is known as the noise level. Here and below, we make the natural assumption that g 1 , g 2 , g δ 1 and g δ 2 all belong to the space L 2 (Γ). The problem we need to solve is: find a source function p in Q ad so that the weak solution u δ ,Ω ≤ c δ. Similar to the noise-free model, we define a Kohn-Vogelius type objective functional , and introduce the following regularized problem.
Remark 1. In the ordinary Tikhonov regularization, noisy measurements are only used in objective functional for data fitting. In our framework here, noisy data, including both Dirichlet and Neumann type, are used in defining the forward problems.
Similar to Problem 1, for ε > 0, Problem 2 admits a unique solution p δ ε ∈ Q ad , which depends continuously on boundary data g δ 1 , g δ 2 and parameter ε. Thus p δ ε → p ε in Q as δ → 0. Moreover, similar to (8) and (14), both the first order necessary and sufficient condition of p δ ε , and the nonlinear projection equation, Recall that we always assume exact measurement g 1 is attainable. Therefore the solution set S of the original noise-free problem is non-empty and p * , the exact solution with minimum norm, exists. Next we discuss the convergence of p δ ε to p * as the noise level δ → 0.
Therefore, from (30)-(32), {(p n , u n 1 , u n 2 )} is uniformly bounded with respect to n in Q × V × V , and there is a subsequence {n } of the sequence {n}, and some elements We verify that u ∞ 1 = u 1 (p ∞ ). From the definition of u n 1 , we have Let n → 0, and use (23) and convergence relations (33) above to get which shows u ∞ 1 = u 1 (p ∞ ). Similarly, we have u ∞ 2 = u 2 (p ∞ ). Subtracting (35) from (34), taking v = u n 1 − u ∞ 1 and using convergence relations in (33) again give to u n Therefore, as n → ∞, where we use the uniform boundness of p n (cf. (30)). Since is a solution of original noise-free inverse source problem (1)- (2). Hence, p ∞ ∈ S.
Next we prove p ∞ = p * . From the lower semi-continuity of norm · 0,Ω0 and the weak convergence of p n to p ∞ , we have Therefore, for any fixed η > 0, there exists a positive integer N such that ∀ n > N , the following relation holds: (36) p n 2 0,Ω0 ≥ p ∞ 2 0,Ω0 − η. We note that (30) also holds when p * is replaced by p ∞ . Therefore, together with (36), −η ≤ p n 2 0,Ω0 − p ∞ 2 0,Ω0 ≤ c δ 2 n ε n holds for big enough n > N . Letting first n → ∞ and then η → 0 in the relation above, we have By the definition of p * , p * 0,Ω0 ≤ p ∞ 0,Ω0 . Combining it with (30) and (36), for big enough n > N , the following relation holds: p n 2 0,Ω0 − p ∞ 2 0,Ω0 ≤ p n 2 0,Ω0 − p * 2 0,Ω0 ≤ c δ 2 n ε n Letting n → ∞ again in the relation above and using (37), we arrive at (38) p ∞ 0,Ω0 = p * 0,Ω0 . Using the definition of p * again, (38) means p ∞ = p * and p n p * in Q, as n → ∞. Thus the limit does not depend on the subsequence selected, and then the entire solution sequence p n 0 in Q, as n → ∞. The strong convergence holds due to lim n→∞ p n 0,Ω0 = p * 0,Ω0 and weak convergence. 4. Study of a mapping. From the discussion in sections 2 and 3, we see that the source function from the regularized problem can be equivalently defined by a system: (5), (6), (9), (10) and (14) for the noise-free case, and (24), (25), (28), (29) and (27) for the noise model. The study of the system is the same for either case, and here for simplicity in notation, we consider the noise-free case. To further simplify the notation, we write g α := g 2 + α g 1 , g β := g 2 + β g 1 .
Convergence of Algorithm 5.1 can be studied similar to the fixed-point iteration, as stated in Theorem 4.1. Indeed, instead of (43), we define Instead of (54), we have and then The inequality (53) still holds and we use it in (57): (58) Π(p 1 ) − Π(p 2 ) 0,Ω0 ≤ (1 − ρ) + c ρ ε −1 (β − α) 2 p 1 − p 2 0,Ω0 . Thus, for ρ ∈ (0, 1), if c 0 > 0 in (44) is small enough, the mapping defined by (39)-(42) and (56) is a contraction, and the statement of Theorem 4.1 holds for this mapping. In particular, we have the convergence of Algorithm 5.1. Note that the the same convergence statement is valid also for the discrete version of the system. We comment that in numerical simulations, convergence occurs even for moderately large values of c 0 .
In practice, Q ad often has one of three forms: (1) Q ad = Q; (2) Q ad = {p ∈ Q | p ≥ 0 a.e. in Ω 0 }; (3) Q ad = {p ∈ Q | p ≤ p ≤ p a.e. in Ω 0 } with p, p ∈ C(Ω 0 ). In these cases, P Q ad has simple forms. Particularly, in the case Q ad = Q, the iteration in Algorithm 5.1 can be avoided. In fact, in this case, P Q ad = I, the identical operator in Q, and p ε = 1 ε (αΦ 1 − βΦ 2 )χ Ω0 . Substituting 1 ε (αΦ 1 − βΦ 2 )χ Ω0 for p back into (5) and (6), and together with (9) and (10), we arrive at a forward system of partial differential equations. Once Φ 1 and Φ 2 are computed from the reduced system, we obtain p ε = 1 ε (αΦ 1 − βΦ 2 )χ Ω0 . We note that all statements above hold for noise model in Section 3, just with g 1 and g 2 being replaced by g δ 1 and g δ 2 . In the next section, finite element methods will be used for obtaining discrete approximate solutions based on Algorithm 5.1 and some numerical experiments are presented. 6. Numerical experiments. In this section, some numerical results are reported based on our new reconstruction framework for the inverse source problem (1)-(2). Standard linear finite element methods are applied to obtain discrete approximate solutions of the BVPs. Let Ω be the problem domain and h be the mesh parameter, i.e., the maximal diameter of the elements in a partition T h of Ω. These results are computed in a Matlab environment. For a two dimensional domain, triangular meshes are produced and refined by using of the PDE Toolbox. For a three dimensional domain, tetrahedral meshes are produced and refined through a COMSOL Multiphysics software. In all examples, biconjugate gradient (BICG) methods are used to solve the resulting linear algebraic system. The Dirichlet data g 1 , which may or may not include noise, are obtained by solving the forward BVP (1) on a rather fine mesh for a given true source p and Neumann data g 2 . Indicated by Proposition 3, without loss of generality, we set with m to be chosen. For better assessing the accuracy of approximate solutions, we define the L 2 -norm relative error in an approximate solution: where p is the true source function, p h ε stands for an approximate solution of p ε (for noise-free data) or p δ ε (for noise data). Example 1. In this example, Ω = {(x, y) ∈ Ω | x 2 +y 2 < 1}, Ω 0 is a circle, centered at (0.55, 0.45) with radius 0.2. The true source function p(x, y) = 1 + x + y in Ω 0 . The Neumann data g 2 = 0.2 on Γ. Then Dirichlet data g 1 is computed through (1) on a mesh with h = 0.0003278, 351232 elements and 176177 nodes. We consider first the case where the source function p is searched in the whole space Q, i.e., Q ad = Q. For fixed m = 10 5 , reconstruction is repeated on a mesh with h = 0.02134, 1372 elements and 722 nodes for different ε, and L 2 -norm relative errors in approximate source functions are reported in fourth column of Table 1. For comparison with methods in [7] and [12], approximate source functions are also computed with these two methods. The associated L 2 -norm relative errors are listed in the second and third columns of Table 1 respectively. In Table 1, symbols 'HCW', 'SH' and 'CGH' mean that results in the corresponding column are obtained with reconstruction models (3), (4) and (7), respectively. Table 1 shows that all the three methods perform well for properly selected parameters. We conclude from Table 1  that, although Proposition 3 is discussed for accurate p ε , it also holds for p h ε , that is, p h ε is uniform with respect to ε when ε is small enough. In other words, with our new method, a reasonable source function can be reconstructed for rather small regularization parameter ε. It is meaningful by noticing the fact that the smaller ε is, the better the regularized problem approaches to the original one is. As a result, without the need for concern of possible instability associated with a very small ε, in our new method, we can choose a small value of ε based on the consideration of the solution accuracy, and the method will work as long as β − α is small enough. To measure the stability of the new K-V type reconstruction model, 5% uniformly distributed noise is added to g 1 . Set m = 10 4 . The above experiments are repeated, and the resulting approximate solution errors are shown in the fourth column of Table 2. Counterpart with methods in [7] and [12] are also computed, and shown in the second and third columns. Table 2 shows that the results for model with noise are similar to those for noise-free case: 1 ε (αΦ δ 1 − βΦ δ 2 ) is uniform with respect to ε for small enough ε; with properly selectly α and β, reasonable approximate source functions could be obtained for rather small ε. More importantly, we can see from it that the reconstruction model is stable. For clarity, Tables 1 and 2 are plotted in   with the best accuracy. Specifically, Figure 2 is p h ε for δ = 0, ε = 10 −7 and m = 10 5 ; Figure 3 is p h ε for δ = 0.05, ε = 10 −11 and m = 10 4 . Example 2. In this example, a 3D source is reconstructed. For this purpose, set the problem domain Ω: Assume Ω 0 is sphere centered at (5,5,6) with unit radius and we place a true source p = 5 in Ω 0 . Then for given g 2 = | sin(x + y + z)|, (1) is used to compute measurement g 1 , on a mesh with h = 0.1355, 162759 tetrahedrons and 29043 nodes. A sketch of computational mesh is shown in Figure 4. With m = 10 3 , Algorithm (5.1) is used again to obtain approximate source functions for different ε on a mesh with h = 0.2210, 8463 tetrahedrons and 1656 nodes. We list the L 2 -norm relative errors, computed with methods in [7], [12] and this paper respectively, in Table 3. As 2D example above, for the 3D problem, similar conclusions can be drawn from Table 3. Table 3. Err(p h ε ) for noise-free data. Again, for verifying the stability of our formulation, 5% random noise is added to Dirichlet data g 1 . Set m = 10 4 , reconstruction above is repeated for this 3D problem. Approximate solution errors are shown in Table 4. Again, we plot Tables 3 and 4 in Figure 5. Finally, the approximate source function p h ε with the best accuracy for noise-free model and 5% noise model are displayed in Figures 6 and  7, respectively. Specifically, Figure 6 is p h ε for δ = 0, ε = 10 −3 and m = 10 4 while Figure 7 is p h ε for δ = 0.05, ε = 10 −3 and m = 10 4 . Conclusion could be drawn from these results that the model explored here is stable.
We note that all numerical results above are obtained for Case 1: Q ad = Q. Using Algorithm 5.1, we repeat the experiments for Case 2: Q ad = {p ∈ Q | p ≥ 0 a.e. in Ω 0 } and Case 3: Q ad = {p ∈ Q | p ≤ p ≤p a.e. in Ω 0 }, and the reconstructed results are similar to those for Case 1. In summary, we conclude from them that the method proposed here is feasible and effective.   7. Concluding remarks. In this paper, a new Kohn-Vogelius type formulation is proposed for the inverse source problem of elliptic PDEs. Different from existing methods in the literature, we introduce two positive parameter to form two different but related BVPs, so that the original boundary fitting is transferred to domain fitting. With two new BVPs, the boundary measurement is incorporated into Robin boundary condition. From the practical point of view, this is meaningful since the effects of noise in forward boundary condition and noise in measurement on solution accuracy are different. The most important feature of the new method is that an accurate approximate source function can be reconstructed for rather small regularization parameter, even in the case where the data is polluted by random noise. Therefore, although in standard regularization methods for the inverse source problem, choosing a proper regularization parameter is delicate issue, our new method allows the use of a very small regularization parameter based mostly on the consideration of the solution accuracy.