A regularization approach for an inverse source problem in elliptic systems from single Cauchy data
AA regularization approach for an inverse source problem in elliptic systemsfrom single Cauchy data
Michael Hinze a , Bernd Hofmann b and Tran Nhan Tam Quyen c a Department of Mathematics, University of Hamburg, 20146 Hamburg, Germany b Faculty of Mathematics, Chemnitz University of Technology, 09107 Chemnitz, Germany c Institute for Numerical and Applied Mathematics, University of Goettingen, 37083 Goettingen, Germany
Abstract:
In this paper we investigate the problem of identifying the source term f in the elliptic system −∇ · (cid:0) Q ∇ Φ (cid:1) = f in Ω ⊂ R d , d ∈ { , } , Q ∇ Φ · (cid:126)n = j on ∂ Ω and Φ = g on ∂ Ωfrom a single noisy measurement couple ( j δ , g δ ) of the Neumann and Dirichlet data ( j, g ) with noise level δ >
0. In thiscontext, the diffusion matrix Q is given. A variational method of Tikhonov-type regularization with specific misfitterm of Kohn-Vogelius-type and quadratic stabilizing penalty term is suggested to tackle this linear inverse problem.The method also appears as a variant of the Lavrentiev regularization. For the occurring linear inverse problem ininfinite dimensional Hilbert spaces, convergence and rate results can be found from the general theory of classicalTikhonov and Lavrentiev regularization. Using the variational discretization concept, where the PDE is discretizedwith piecewise linear and continuous finite elements, we show the convergence of finite element approximations tosolutions of the regularized problem. Moreover, we derive an error bound and corresponding convergence ratesprovided a suitable range-type source condition is satisfied. For the numerical solution we propose a conjugategradient method. To illustrate the theoretical results, a numerical case study is presented which supports ouranalytical findings. Key words and phrases:
Inverse source problem, Tikhonov and Lavrentiev regularization, finite element method,source condition, convergence rates, ill-posedness, conjugate gradient method, Neumann problem, Dirichlet problem.
AMS Subject Classifications:
Let Ω be an open, bounded and connected domain of R d , d ∈ { , } , with Lipschitz boundary ∂ Ω. Weconsider the elliptic system −∇ · (cid:0) Q ∇ Φ (cid:1) = f in Ω , (1.1) Q ∇ Φ · (cid:126)n = j † on ∂ Ω and (1.2)Φ = g † on ∂ Ω , (1.3)where (cid:126)n is the unit outward normal on ∂ Ω and the diffusion matrix Q is given. Furthermore, we assumethat Q := ( q rs ) ≤ r,s ≤ d ∈ L ∞ (Ω) d × d is symmetric and satisfies the uniformly ellipticity condition Q ( x ) ξ · ξ = (cid:88) ≤ r,s ≤ d q rs ( x ) ξ r ξ s ≥ q | ξ | a.e. in Ω (1.4)for all ξ = ( ξ r ) ≤ r ≤ d ∈ R d with some constant q > j † ∈ H − / ( ∂ Ω) := H / ( ∂ Ω) ∗ , g † ∈ H / ( ∂ Ω), and the source term f ∈ L (Ω) are given, then there maybe no Φ satisfying this system. In this paper we assume that the system is consistent and our aim is toreconstruct a function f ∈ L (Ω) in the system (1.1)–(1.3) from a noisy measurement couple ( j δ , g δ ) ∈ H − / ( ∂ Ω) × H / ( ∂ Ω) of the exact Neumann and Dirichlet data (cid:0) j † , g † (cid:1) , where δ > (cid:13)(cid:13) j δ − j † (cid:13)(cid:13) H − / ( ∂ Ω) + (cid:13)(cid:13) g δ − g † (cid:13)(cid:13) H / ( ∂ Ω) ≤ δ. (1.5) Email: [email protected], [email protected], [email protected] a r X i v : . [ m a t h . A P ] M a r he source identification problem in PDEs arises in many branches of applied science such as electroen-cephalography, geophysical prospecting and pollutant detection, and attracted great attention from manyscientists in the last 30 years or so. For surveys on this subject we may consult in [7, 15, 18, 22, 25, 42] andthe references therein. Up to now, only a limited number of works was investigated the general source iden-tification problem and obtained results concentrated on numerical analysis for the identification problem. In[19, 31, 32] authors have used the dual reciprocity boundary element methods to simulate numerically for theabove mentioned identification problem. In case some priori knowledge of the identified source is available,such as a point source, a characteristic function or a harmonic function, numerical methods treating theproblem have been obtained in [5, 6, 30, 39]. A survey of the problem of simultaneously identifying thesource term and coefficients in elliptic systems from distributed observations can be found in [38], wherefurther references can be found.In the present paper, the general source identification problem in elliptic partial differential equations froma single noisy measurement couple of Neumann and Dirichlet data is studied. So far, we have not yet foundinvestigations on the discretization analysis for this source recovery problem, a fact which also motivatedthe research presented in the paper. By using a suitable version of the Tikhonov-type regularization withsome non-standard misfit term we could outline that the source distribution inside the physical domainΩ can be reconstructed from a finite number of observations on the boundary ∂ Ω, at least by numericalapproximations. The specific regularization approach proves to be a version of Lavrentiev regularizationwith implicit forward operator. One of the main results of the paper is to show convergence of the finiteelement discretized Tikhonov-regularized solutions to a sought source function. Another main result isthe interpretation of an occurring condition of solution smoothness as a range-type source condition ofLavrentiev’s regularization method. This allows us to establish error bounds and corresponding convergencerates for the regularized solutions.To formulate precisely the problem, we first give some notations. Let us denote by γ : H (Ω) → H / ( ∂ Ω)the continuous Dirichlet trace operator with γ − : H / ( ∂ Ω) → H (Ω) its continuous right inverse operator,i.e. ( γ ◦ γ − ) g = g for all g ∈ H / ( ∂ Ω). We set H (cid:5) (Ω) := (cid:26) u ∈ H (Ω) (cid:12)(cid:12)(cid:12) (cid:90) ∂ Ω γudx = 0 (cid:27) and H / (cid:5) ( ∂ Ω) := (cid:26) g ∈ H / ( ∂ Ω) (cid:12)(cid:12)(cid:12) (cid:90) ∂ Ω g ( x ) dx = 0 (cid:27) and denote by C Ω the positive constant appearing in the Poincar´e-Friedrichs inequality (cf. [35]) C Ω (cid:90) Ω ϕ ≤ (cid:90) Ω |∇ ϕ | for all ϕ ∈ H (cid:5) (Ω) . (1.6)Since H (Ω) := (cid:8) u ∈ H (Ω) | γu = 0 (cid:9) ⊂ H (cid:5) (Ω), the inequality (1.6) is in particular valid for all ϕ ∈ H (Ω).Furthermore, by (1.4), the coercivity condition (cid:107) ϕ (cid:107) H (Ω) ≤ C Ω C Ω (cid:90) Ω |∇ ϕ | ≤ C Ω C Ω q (cid:90) Ω Q ∇ ϕ · ∇ ϕ (1.7)holds for all ϕ ∈ H (cid:5) (Ω).Now, for any fixed ( j, g ) ∈ H − / ( ∂ Ω) × H / (cid:5) ( ∂ Ω) we can simultaneously consider the Neumann problem −∇ · ( Q ∇ u ) = f in Ω and Q ∇ u · (cid:126)n = j on ∂ Ω (1.8)as well as the Dirichlet problem −∇ · ( Q ∇ v ) = f in Ω and v = g on ∂ Ω . (1.9)By the aid of (1.7) and the Riesz representation theorem, we conclude that for each f ∈ L (Ω) there existsa unique weak solution u of the problem (1.8) in the sense that u ∈ H (cid:5) (Ω) and satisfies the identity (cid:90) Ω Q ∇ u · ∇ ϕ = (cid:104) j, γϕ (cid:105) + ( f, ϕ ) (1.10)for all ϕ ∈ H (cid:5) (Ω), where notation (cid:104) j, g (cid:105) stands for the value of the function j ∈ H − / ( ∂ Ω) at g ∈ H / ( ∂ Ω)and the notation ( f, ϕ ) is the inner product of f and ϕ in the space L (Ω). Then we can define the Neumannoperator N : L (Ω) → H (cid:5) (Ω) with f (cid:55)→ N f j, f ∈ L (Ω) to the unique weak solution N f j := u of the problem (1.8). Similarly, theproblem (1.9) also attains a unique weak solution v in the sense that v ∈ H (Ω), γv = g and the identity (cid:90) Ω Q ∇ v · ∇ ψ = ( f, ψ ) (1.11)holds for all ψ ∈ H (Ω). The Dirichlet operator is defined as D : L (Ω) → H (cid:5) (Ω) with f (cid:55)→ D f g, which maps each f ∈ L (Ω) to the unique weak solution D f g := v of the problem (1.9). Therefore, for anyfixed f ∈ L (Ω) we can define the so-called Neumann-to-Dirichlet map Λ f : H − / ( ∂ Ω) → H / (cid:5) ( ∂ Ω) , j (cid:55)→ Λ f j := γ N f j. We mention that since H (Ω) ⊂ H (cid:5) (Ω), we from (1.10) have that (cid:82) Ω Q ∇N f j ·∇ ψ = ( f, ψ ) for all ψ ∈ H (Ω).In view of (1.11) we therefore conclude Λ f j = g if and only if N f j = D f g , where the identities N f j = N f N j and D f g = D f D g (1.12)are satisfied, and the operators f (cid:55)→ N f f (cid:55)→ D f L (Ω) into itself.Furthermore, Λ f j = γ N f j = γ N j + γ N f j + Λ f
0, where Λ j is linear, self-adjoint, bounded andinvertible, as the diffusion Q is smooth enough (cf. [33]).As in electrical impedance tomography (EIT) or for the Calder´on’s problem [4, 14, 33] one can pose thequestion whether the source distribution f inside a physical domain Ω can be determined from an infinite number of observations on the boundary ∂ Ω, i.e. from the Neumann-to-Dirichlet map Λ f : f , f ∈ L (Ω) with Λ f = Λ f ⇒ f = f ?To the best of our knowledge, the above question is still open so far. In case an observation Λ δ of Λ f beingavailable one can use a certain regularization method to approximate the sought source. For example, onecan consider for operator norms (cid:107) · (cid:107) ∗ a minimizer of the problemmin f ∈ L (Ω) (cid:107) Λ f − Λ δ (cid:107) ∗ + ρ (cid:107) f − f ∗ (cid:107) L (Ω) as a reconstruction along the lines of Tikhonov’s regularization method, where ρ > f ∗ is an a-priori estimate of the sought source.However, in practice we have only a finite number of observations and the task is to reconstruct the identifiedsource, at least by numerical approximations. Furthermore, for simplicity of exposition we below restrictourselves to the case of just one observation pair ( j δ , g δ ) being available, while the approach described herecan be easily extended to multiple measurements (cid:0) j iδ , g iδ (cid:1) i =1 ,...,I , see Section 6, Ex. 6.2. The inverse problemis thus stated as follows. Given (cid:0) j † , g † (cid:1) ∈ H − / ( ∂ Ω) × H / (cid:5) ( ∂ Ω) with Λ f j † = g † , f ind f ∈ L (Ω) . ( IP )In other words, the interested problem is, for given (cid:0) j † , g † (cid:1) ∈ H − / ( ∂ Ω) × H / (cid:5) ( ∂ Ω), to find some f ∈ L (Ω)and consequently Φ ∈ H (cid:5) (Ω) such that the system (1.1)–(1.3) is satisfied in the weak sense. Precisely, wedefine the general solution set I (cid:0) j † , g † (cid:1) := (cid:8) f ∈ L (Ω) (cid:12)(cid:12) Λ f j † = g † (cid:9) = (cid:8) f ∈ L (Ω) (cid:12)(cid:12) N f j † = D f g † (cid:9) (1.13)of the inverse problem ( IP ). The source identification problem as described here is well known to be notuniquely determined from boundary observations (see a counterexample in [3]), i.e. the set I (cid:0) j † , g † (cid:1) fails tobe a singleton. Since not the Neumann-to-Dirichlet map is given, but only one pair (cid:0) j † , g † (cid:1) , the problem iseven highly underdetermined. Thus instead we will search for the uniquely determined f ∗ -minimum-normsolution f † , which is the minimizer of the problemmin f ∈I ( j † ,g † ) (cid:107) f − f ∗ (cid:107) L (Ω) . ( IP −
M N )3s a consequence of item (iii) of Lemma 2.1 below, the set I (cid:0) j † , g † (cid:1) is non-empty, closed and convex, hence f † is uniquely determined. On the other hand, for all f ∈ I (cid:0) j † , g † (cid:1) the equation N f j † = D f g † is fulfilled.However, we have to solve this equation with noise data ( j δ , g δ ) ∈ H − / ( ∂ Ω) × H / (cid:5) ( ∂ Ω) of (cid:0) j † , g † (cid:1) satisfying(1.5). The simplest variety of regularization may be to consider a minimizer of the Tikhonov functional (cid:107)N f j δ − D f g δ (cid:107) L (Ω) + ρ (cid:107) f − f ∗ (cid:107) L (Ω) (1.14)over f ∈ L (Ω) as an approximation solution to f † .In present work we adopt the variational approach of Kohn and Vogelius [27, 28, 29] in using cost functionalcontaining the gradient of forward operators to the above mentioned inverse source problem. More precisely,we use the convex functional J δ ( f ) := (cid:90) Ω Q ∇ ( N f j δ − D f g δ ) · ∇ ( N f j δ − D f g δ ) dx, (1.15)instead of the mapping f (cid:55)→ (cid:107)N f j δ − D f g δ (cid:107) L (Ω) , together with Tikhonov regularization and consider the unique solution f ρ,δ of the strictly convex minimization problemmin f ∈ L (Ω) Υ ρ,δ ( f ) with Υ ρ,δ ( f ) := J δ ( f ) + ρ (cid:107) f − f ∗ (cid:107) L (Ω) , ( P ρ,δ )where the gradient of the functional Υ ρ,δ can be explicitly written as12 ∇ Υ ρ,δ ( f ) = N f j δ − D f g δ + ρ ( f − f ∗ ) for all f ∈ L (Ω) . (1.16)The motivation in using this cost functional J δ as misfit functional is that for all ξ ∈ L (Ω) the inequality J ( ξ ) := (cid:90) Ω Q ∇ (cid:0) N ξ j † − D ξ g † (cid:1) · ∇ (cid:0) N ξ j † − D ξ g † (cid:1) dx ≥ C Ω q C Ω (cid:13)(cid:13) N ξ j † − D ξ g † (cid:13)(cid:13) H (Ω) ≥ J ( f ) = 0 at any f ∈ I (cid:0) j † , g † (cid:1) . The advantage is evident, because the minimizer f ρ,δ ∈ L (Ω)satisfies the equation ∇ Υ ρ,δ ( f ρ,δ ) = 0 such that, for j := j δ , g := g δ and f := f ρ,δ , we have L j,g ( f ) := N f j − D f g + ρ ( f − f ∗ ) = 0 , (1.17)and hence, for f := f ρ,δ , we have f − f ∗ = − ρ ( N f j δ − D f g δ ) . (1.18)Due to formula (1.18), the Tikhonov regularization approach under consideration with specific misfit termalso appears as a variant of the Lavrentiev regularization (see, e.g., [2, 12, 24, 43]). After some operator-theoretic settings and preliminary results in Section 2, concerning also the ill-posedness of the linear inverseproblem under consideration, we apply in Section 3 the general theory of classical Tikhonov and Lavrentievregularization for such problems yielding propositions on convergence and convergence rates for the regular-ized solutions in the infinite dimensional Hilbert spaces. However, for convenience in numerical analysis withthe finite element methods introduced in Section 4 our focus is here on the extremal problem for minimizingthe Tikhonov functional with Kohn-Vogelius misfit term and quadratic penalty. The use of different convexpenalty terms, e.g. total variation, may be a work for us in future.Let N hf j δ and D hf g δ be corresponding approximations of the solution maps N f j δ and D f g δ in the finitedimensional space V h of piecewise linear, continuous finite elements. We then consider the discrete regularizedproblem corresponding to ( P ρ,δ ), i.e., the following strictly convex minimization problemmin f ∈ L (Ω) (cid:90) Ω Q ∇ (cid:0) N hf j δ − D hf g δ (cid:1) · ∇ (cid:0) N hf j δ − D hf g δ (cid:1) dx + ρ (cid:107) f − f ∗ (cid:107) L (Ω) . (cid:16) P hρ,δ (cid:17) Using the variational discretization concept introduced in [23], we show in Section 4 that the unique solution f hρ,δ of the problem (cid:16) P hρ,δ (cid:17) automatically belongs to the finite dimensional space V h . Thus, a discretizationof the admissible set L (Ω) can be avoided. 4s h, δ → ρ = ρ ( h, δ ), also in Section 4,we prove that the sequence (cid:0) f hρ,δ (cid:1) converges to f † in the L (Ω)-norm. Furthermore, the corresponding statesequences (cid:16) N hf hρ,δ j δ (cid:17) and (cid:16) D hf hρ,δ g δ (cid:17) converge in the H (Ω)-norm to Φ † = Φ † ( f † , j † , g † ) solving (1.1)–(1.3).Section 5 is devoted to convergence rates for the discretized problem. In this section we also show that if f ∈ I (cid:0) j † , g † (cid:1) and there is a function w ∈ L (Ω) such that f − f ∗ = L j † ,g † ( w ), or in other notation f − f ∗ = N w j † − D w g † , (1.19)then f = f † , i.e. f is the unique f ∗ -minimum-norm solution of the identification problem. Condition(1.19) appears to be a source condition for both, Tikhonov and Lavrentiev regularization, and allows forcorresponding convergence rates of the continuous setting in infinite dimensional spaces as well as afterincorporating the discretization. In the latter case, precisely for the known matrix Q ∈ C , (Ω) d × d and theexact data (cid:0) j † , g † (cid:1) ∈ H / ( ∂ Ω) × H / ( ∂ Ω), we derive the convergence rates (cid:13)(cid:13)(cid:13) N hf hρ,δ j δ − D hf hρ,δ g δ (cid:13)(cid:13)(cid:13) H (Ω) + ρ (cid:13)(cid:13) f hρ,δ − f † (cid:13)(cid:13) L (Ω) = O (cid:0) δ + h + hρ + δρ + ρ (cid:1) and (cid:13)(cid:13)(cid:13) N hf hρ,δ j δ − Φ † (cid:13)(cid:13)(cid:13) H (Ω) + (cid:13)(cid:13)(cid:13) D hf hρ,δ g δ − Φ † (cid:13)(cid:13)(cid:13) H (Ω) = O (cid:0) δ ρ − + h ρ − + h + δ + ρ (cid:1) . Finally, for the numerical solution of the discrete regularized problem (cid:16) P hρ,δ (cid:17) we employ in Section 6 aconjugate gradient algorithm. Numerical case studies illustrate the analytical results and show the efficiencyof our theoretical findings.We conclude this introduction with a remark that since the main interest is to clearly state our ideas, weonly treat the model elliptic problem (1.1) while the approach described here can be easily extended to moregeneral models, e.g., for the source identification problem in diffusion-reaction equations −∇ · (cid:0) Q ∇ Φ (cid:1) + κ Φ = f in Ω , Q ∇ Φ · (cid:126)n + σ Φ = j † on ∂ Ω and Φ = g † on ∂ Ω (1.20)from a measurement ( j δ , g δ ) of (cid:0) j † , g † (cid:1) , where Q satisfying the condition (1.4), 0 (cid:54) = κ = κ ( x ) ∈ L ∞ (Ω), i.ethe set { x ∈ Ω | κ ( x ) (cid:54) = 0 } has positive Lebesgue measure, and σ = σ ( x ) ∈ L ∞ ( ∂ Ω) with σ ≥ (cid:90) Ω Q ∇ ( R f j δ − D f g δ ) · ∇ ( R f j δ − D f g δ ) dx + (cid:90) Ω κ ( R f j δ − D f g δ ) dx + (cid:90) ∂ Ω σ ( R f j δ − D f g δ ) dx over f ∈ L (Ω), where R and D are the Robin operator and the Dirichlet operator relating with the equation(1.20), respectively.We here would like to mention an inverse problem related closely to the identification in this paper, theproblem of identifying the source term f in the Helmholtz-type equation ∇ · (cid:0) Q ∇ Φ (cid:1) + κ Φ = f in Ωfrom measured Cauchy data ( j δ ( κ ) , g δ ( κ )) which is available for all frequency κ >
0. The uniqueness resultsfor this identification problem can be found in [3, 8], while several effective recovered algorithms have beenpresented in [1, 9].Throughout the paper we use the standard notion of Sobolev spaces H (Ω), H (Ω), W k,p (Ω), etc from, forexample, [44]. If not stated otherwise we write (cid:82) Ω · · · instead of (cid:82) Ω · · · dx . In order to define appropriate operators, we recall the decompositions (1.12) of the corresponding Neumannand Dirichlet problems, where N f D f f ∈ L (Ω). On the other5and, N j and D g depend nonlinearly on j and g , respectively, but both are independent of f . Hence, thedifference of N f j and D f g characterizes, for fixed elements j and g , the affine mapping L j,g of f ∈ L (Ω)defined by formula (1.17). First we introduce the linear operator ˜T : L (Ω) → H (cid:5) (Ω) defined as ˜T ( f ) := N f − D f ∈ H (cid:5) (Ω) . Since the image elements N f −D f ∈ H (cid:5) (Ω) also belong to L (Ω), one can moreover introduce the operator T : L (Ω) → L (Ω) defined by T ( f ) := N f − D f ∈ L (Ω) , (2.1)where T ( f ) = ˜T ( f ) for all f ∈ L (Ω). On the other hand, we remark that the expression[ u, v ] := (cid:90) Ω Q ∇ u · ∇ v (2.2)generates an inner product on the space H (cid:5) (Ω) which is equivalent to the usual one. Now let ˜T ∗ : H (cid:5) (Ω) → L (Ω)be the adjoint operator of ˜T : L (Ω) → H (cid:5) (Ω), where H (cid:5) (Ω) is equipped with the inner product (2.2) above.For all f ∈ L (Ω) and φ ∈ H (cid:5) (Ω) we thus have (cid:104) ˜T f, φ (cid:105) = (cid:90) Ω Q ∇N f · ∇ φ − (cid:90) Ω Q ∇D f · ∇ φ = ( f, φ ) − (cid:90) Ω Q ∇D f · ∇ φ = ( f, ˜T ∗ φ ) , (2.3)by (1.10). We now decompose H (cid:5) (Ω) into the orthogonal direct sum H (cid:5) (Ω) = H (Ω) ⊕ H (Ω) ⊥ with respectto the inner product (2.2). We note for all g ∈ H / (cid:5) ( ∂ Ω) that D g ∈ H (Ω) ⊥ . Furthermore, ∀ g , g ∈ H / (cid:5) ( ∂ Ω) , g (cid:54) = g ⇒ D g (cid:54) = D g which implies dim H (Ω) ⊥ ≥ dim H / (cid:5) ( ∂ Ω) = ∞ . For all f ∈ L (Ω) we deduce from (1.11) and (2.3) that φ ∈ H (Ω) ⇔ (cid:90) Ω Q ∇D f · ∇ φ = ( f, φ ) ⇔ ( f, ˜T ∗ φ ) = 0 ⇔ ˜T ∗ φ = 0 ⇔ φ ∈ ker ˜T ∗ , or in other words ker ˜T ∗ = H (Ω). Furthermore, for all (cid:98) φ ∈ H (Ω) ⊥ we get (cid:82) Ω Q ∇D f · ∇ (cid:98) φ = 0, since D f ∈ H (Ω). Again, the equation (2.3) implies that( f, (cid:98) φ ) = ( f, ˜T ∗ (cid:98) φ ) for all f ∈ L (Ω) , (cid:98) φ ∈ H (Ω) ⊥ . Therefore, ˜T ∗| H (Ω) ⊥ is the compact embedding H (Ω) ⊥ (cid:44) → L (Ω) and ˜T ∗ is the composition of the projectorfrom H (cid:5) (Ω) onto H (Ω) ⊥ and the embedding operator from H (Ω) to L (Ω). Furthermore we have, for all f ∈ L (Ω), T ( f ) = ˜T ∗ [ ˜T ( f )] , (2.4)because range( ˜T ) is orthogonal to ker ˜T ∗ and hence ˜T ∗ acts only as embedding operator. Lemma 2.1. (i) The operator T defined by formula (2.1) is linear, bounded, self-adjoint and non-negative,i.e. we have ( T ( f ) , w ) = ( f, T ( w )) and ( T ( f ) , f ) ≥ for all f, w ∈ L (Ω) . (2.5) Moreover, T is compact and has an infinite dimensional range which is non-closed, i.e. we have range ( T ) (cid:54) = range ( T ) .(ii) For any fixed ( j, g ) ∈ H − / ( ∂ Ω) × H / ( ∂ Ω) the map L j,g : L (Ω) → L (Ω) defined by L j,g ( f ) := N f j − D f g = T ( f ) + N j − D g is affine linear, continuous and monotone, i.e. we have ( L j,g ( f ) − L j,g ( w ) , f − w ) ≥ for all f, w ∈ L (Ω) . (2.6) (iii) The solution set I (cf. (1.13)) is a closed affine subspace of the Hilbert space L (Ω) . roof. (i) It follows from (1.10) that( T ( f ) , w ) = (cid:90) Ω Q ∇N w · ∇ ( N f − D f
0) = (cid:90) Ω Q ∇N w · ∇N f − (cid:90) Ω Q ∇N w · ∇D f , (2.7)and similarly, ( f, T ( w )) = (cid:90) Ω Q ∇N f · ∇N w − (cid:90) Ω Q ∇N f · ∇D w . (2.8)Using (1.10)–(1.11) again, we get (cid:90) Ω Q ∇N w · ∇D f w, D f
0) = (cid:90) Ω Q ∇D w · ∇D f (cid:90) Ω Q ∇N f · ∇D w f, D w
0) = (cid:90) Ω Q ∇D f · ∇D w T now follows directly from (2.7)-(2.9). We further have from (1.10)–(1.11)for all f ∈ L (Ω) that (cid:82) Ω Q ∇D f · ∇ ( N f − D f
0) = 0 . Combining this with the identity ( T ( f ) , f ) = (cid:82) Ω Q ∇N f · ∇ ( N f − D f
0) we arrive at ( T ( f ) , f ) = (cid:82) Ω Q ∇ ( N f − D f · ∇ ( N f − D f ≥ T is compact and has an infinite dimensional range. The operator T : L (Ω) → L (Ω)as a composition of a bounded linear operator ˜T and a compact embedding operator is compact. Next, weshow that dim( ˜T ) = ∞ . For deriving a contradiction we assume that dim range( ˜T ) < ∞ . Then we can write H (cid:5) (Ω) = range( ˜T ) ⊕ range( ˜T ) ⊥ with respect to the inner product (2.2). By (2.3), for all ϕ ∈ range( ˜T ) ⊥ we get (cid:82) Ω Q ∇D f · ∇ ϕ = ( f, ϕ ) holding for all f ∈ L (Ω) which implies that ϕ ∈ H (Ω). Therefore, H (Ω) ⊥ ⊂ (cid:0) range( ˜T ) ⊥ (cid:1) ⊥ = range( ˜T ) = range( ˜T ) and this yields the contradiction ∞ = dim H (Ω) ⊥ ≤ dim range( ˜T ) < ∞ . Consequently, T : L (Ω) → L (Ω) is compact operator and possesses an infinitedimensional range.(ii) The inequality (2.6) follows directly from (2.5).(iii) Since T is a bounded linear operator, the solution set can be written as I ( j † , g † ) = (cid:8) f = f ⊕ f ⊥ ∈ L (Ω) (cid:12)(cid:12) f ∈ ker( T ) , f ⊥ = T † ( D g † − N j † ) (cid:9) with the Moore-Penrose pseudoinverse T † of T . Then the nullspace ker( T ) is a closed subspace and T † ( D g † − N j † ) a well-defined element in L (Ω). Consequently, I ( j † , g † ) is a non-empty closed affinesubspace and hence also a convex set in the Hilbert space L (Ω).Due to item (ii) of Lemma 2.1, we can reformulate the identification problem ( IP ) as an operator equationwith linear bounded self-adjoint non-negative operator T mapping in L (Ω). Finding an element f ∈I (cid:0) j † , g † (cid:1) is then equivalent to solving the linear operator equation T ( f ) = κ ( j † , g † ) , where κ ( j, g ) := D g − N j. (2.10)This makes the inverse problem explicit, but we have to take into account that instead of the exact right-hand side κ ( j † , g † ) only noisy data κ ( j δ , g δ ) satisfying (1.5) for j δ and g δ are available. As a consequence ofitem (i) of Lemma 2.1 we see that the equation (2.10) formulated in the Hilbert space L (Ω) is ill-posed oftype II in the sense of Nashed (cf. [34]). Stable approximate (regularized) solutions f ρ,δ to equation (2.10)satisfy with f = f ρ,δ the auxiliary linear operator equation T ( f ) + ρ ( f − f ∗ ) = κ ( j δ , g δ ) (2.11)in L (Ω) with some regularization parameter ρ > j δ , g δ under consideration and f ∈ L (Ω), anddue to (1.10) – (1.11) T ( f ) ∈ H (Ω) ⊥ ⊂ H (cid:5) (Ω) as well as κ ( j δ , g δ ) ∈ H (Ω) ⊥ , which means that the elements κ ( j δ , g δ ) and ˜T ∗ [ κ ( j δ , g δ )] in L (Ω) coincide. Nevertheless, we have todistinguish the two cases κ ( j δ , g δ ) ∈ H (cid:5) (Ω) and κ ( j δ , g δ ) ∈ L (Ω) in the following lemma.7 emma 2.2. Under the noise model (1.5) there is a constant ≤ ˜ K < ∞ independent of δ such that (cid:107) κ ( j δ , g δ ) − κ ( j † , g † ) (cid:107) H (Ω) ≤ ˜ K δ.
Moreover, there is also a constant ≤ K < ∞ independent of δ such that (cid:107) κ ( j δ , g δ ) − κ ( j † , g † ) (cid:107) L (Ω) ≤ K δ.
Proof.
By Lemma 4.5 below the existence of such constant ˜ K in the first estimate of this lemma follows withthe settings h := 0, f = f := 0 and ˜ K := max {C N , C D } . Then we have under the noise model (1.5) (cid:107) κ ( j δ , g δ ) − κ ( j † , g † ) (cid:107) H (Ω) ≤ (cid:107)N j δ − N j † (cid:107) H (Ω) + (cid:107)D j δ − D j † (cid:107) H (Ω) ≤ ˜ K (cid:0) (cid:107) j δ − j † (cid:107) H − / ( ∂ Ω) + (cid:107) g δ − g † (cid:107) H / ( ∂ Ω) (cid:1) ≤ ˜ Kδ.
By setting K := ˜ K (cid:107) ˜T (cid:107) the second estimate of the lemma gets established. This completes the proof.We conclude this section by mentioning that the functional J δ defined by (1.15) is convex and weaklysequentially lower semi-continuous. In fact, the above defined operator ˜T : L (Ω) → H (cid:5) (Ω) is compact (seethe proof of Lemma 2.1). Using the equivalent inner product (2.2), we therefore conclude that J δ (cid:0) f (cid:1) = [ N f j δ − D f g δ , N f j δ − D f g δ ] = (cid:104) ˜T ( f ) − κ ( j δ , g δ ) , ˜T ( f ) − κ ( j δ , g δ ) (cid:105) , (2.12)which shows the convexity and weak sequential lower semi-continuity of the functional J δ and is the basisfor a classical Tikhonov regularization approach in the subsequent section. Moreover, we have as a basis forLavrentiev regularization in the subsequent section J δ (cid:0) f (cid:1) + ρ (cid:107) f − f ∗ (cid:107) L (Ω) = ( T ( f ) , f ) + ρ (cid:107) f (cid:107) L (Ω) − f, κ ( j δ , g δ ) + ρ f ∗ ) + const ., (2.13)where the constant is independent of f , and it is well-known that, due to the properties of T from Lemma 2.1,the unique minimizer f ρ,δ of this functional coincides with the unique solution of the operator equation (2.11).We close this section by the following note. As discussed in the Introduction section, instead of usingKohn and Vogelius’ function (1.15) we can use the least squares function (cf. (1.14)), and then the jointedminimization problem reads asmin f ∈ L (Ω) Θ( f ) with Θ( f ) := (cid:107) T ( f ) − κ ( j δ , g δ ) (cid:107) L (Ω) + ρ (cid:107) f − f ∗ (cid:107) L (Ω) (2.14)where T ( f ) and κ ( j δ , g δ ) were defined by (2.1) and (2.10), respectively. One can show easily that the problemattains a unique solution ¯ f . Furthermore, for computation this minimizer ¯ f in practice one must derive the L -gradient ∇ Θ( ¯ f ). Let ξ ∈ L (Ω) be arbitrary. We will compute briefly the differential Θ (cid:48) ( ¯ f ) ξ as follows.We have that Θ (cid:48) ( ¯ f ) ξ = (cid:0) T ( ¯ f ) − κ ( j δ , g δ ) , T ( ξ ) (cid:1) + ρ ( ¯ f − f ∗ , ξ ). Consider the adjoint problem −∇ · (cid:0) Q ∇ ¯Φ (cid:1) = T ( ¯ f ) − κ ( j δ , g δ ) in Ω , (2.15) Q ∇ ¯Φ · (cid:126)n = 0 on ∂ Ω . (2.16)(We here do not use the homogeneous Dirichlet boundary condition ¯Φ = 0 on ∂ Ω instead of (2.16), becausein general T ( ξ ) / ∈ H (Ω).) Next, we decompose¯Φ = ¯Φ ⊕ ¯Φ ∈ H (Ω) ⊕ H (Ω) ⊥ (2.17)with respect to the inner product (2.2). Then we obtain that (cid:0) T ( ¯ f ) − κ ( j δ , g δ ) , T ( ξ ) (cid:1) = (cid:90) Ω Q ∇ ¯Φ · ∇ T ( ξ ) = (cid:90) Ω Q ∇N ξ · ∇ ¯Φ − (cid:90) Ω Q ∇D ξ · ∇ ¯Φ − (cid:90) Ω Q ∇D ξ · ∇ ¯Φ (cid:124) (cid:123)(cid:122) (cid:125) =0 = ( ξ, ¯Φ) − ( ξ, ¯Φ ) = ( ¯Φ , ξ ) . As a result, we arrive at ∇ Θ( ¯ f ) = ¯Φ + ρ ( ¯ f − f ∗ ). Compared with (1.16), by utilizing Kohn and Vogelius’function (1.15), we here avoid any computations for the adjoint problem (2.15)–(2.16). Furthermore, weavoid computing numerically for the terms ¯Φ and ¯Φ in the decomposition (2.17) which seems to be stillvery difficult for us. 8 Lavrentiev regularization versus Tikhonov regularization for thecontinuous problem in infinite dimensional Hilbert spaces
In this section, we consider the error analysis of stable approximate solutions to the continuous problem(2.10) in infinite dimensional Hilbert spaces by distinguishing the situations that the right-hand side κ isconsidered as an element in L (Ω) and alternatively as an element in H (Ω).First our focus is on the Lavrentiev regularization approach based on formula (2.13), in which κ ( j δ , g δ ) ∈ L (Ω). The general theory of linear Lavrentiev regularization (see, e.g., [43] and also [2, 12, 24, 36, 37])yields convergence and convergence rates results for the error of regularized solutions f ρ,δ with respect to theuniquely determined f ∗ -minimizing solution f † to problem ( IP −
M N ). Taking into account that Lemma 2.2holds, we immediately derive (see, e.g., [12, Rem. 3.3]) the error estimate (cid:107) f ρ,δ − f † (cid:107) L (Ω) ≤ ρ (cid:107) ( T + ρ I ) − ( f † − f ∗ ) (cid:107) L (Ω) + Kδρ ≤ (cid:107) f † − f ∗ (cid:107) L (Ω) + Kδρ . (3.1)This immediately yields (see [24, section 2]) the following convergence assertion.
Proposition 3.1.
For any a priori parameter choice ρ ( δ ) of the regularization parameter satisfying theconditions ρ ( δ ) → and δρ ( δ ) → as δ → we have for a sequence δ n → , associated data j δ n , g δ n , and associated Lavrentiev-regularized solutions f ρ ( δ n ) ,δ n that lim n →∞ (cid:107) f ρ ( δ n ) ,δ n − f † (cid:107) L (Ω) = 0 , i.e. the regularized solutions are strongly convergent in L (Ω) to the f ∗ -minimum norm solution f † . We can also apply the following well-known result on convergence rates from [43, Theorem 2.2]:
Proposition 3.2.
If there is a source element v ∈ L (Ω) such that the range-type source condition f † − f ∗ = T ( v ) (3.3) is satisfied, we have for an a priori choice ρ ( δ ) ∼ √ δ of the regularization parameter ρ the convergence rate (cid:107) f ρ ( δ ) ,δ − f † (cid:107) L (Ω) = O ( √ δ ) as δ → of Lavrentiev-regularized solutions. Corollary 3.3.
The rate result (3.4) remains true if we have some element w ∈ L (Ω) such that f † − f ∗ = N w j † − D w g † . (3.5) Proof.
Evidently we have N w j † − D w g † = T ( w ) − κ ( j † , g † ) and f † is a solution to the operator equation(2.10) with exact right-hand side κ ( j † , g † ). Hence (3.5) can be rewritten as f † − f ∗ = T ( w − f † ). This yields(3.3) with the new source element v := w − f † . Remark 3.4.
It was shown by a saturation theorem in [36] that (3.4) is the best possible rate for linearLavrentiev regularization, with the exception of singular situations with respect to the forward operator,here T , and with respect to the solution f † .Revisiting the source condition (3.5), we add the following remarkable result. Proposition 3.5.
Assume that f ∈ I (cid:0) j † , g † (cid:1) and that there is a function w ∈ L (Ω) such that f − f ∗ = N w j † −D w g † . Then f is the uniquely determined f ∗ -minimum-norm solution of problem ( IP −
M N ) ,i.e. we have f = f † . roof. We have with ξ ∈ (cid:8) ξ ∈ L (Ω) (cid:12)(cid:12) N ξ j † = D ξ g † (cid:9) that12 (cid:107) ξ − f ∗ (cid:107) L (Ω) − (cid:107) f − f ∗ (cid:107) L (Ω) = 12 (cid:107) ξ − f (cid:107) L (Ω) + ( f − f ∗ , ξ − f ) ≥ ( f − f ∗ , ξ − f ) = (cid:0) N w j † − D w g † , ξ − f (cid:1) = (cid:90) Ω Q ∇N ξ j † · ∇ (cid:0) N w j † − D w g † (cid:1) − (cid:10) j † , γ (cid:0) N w j † − D w g † (cid:1)(cid:11) − (cid:90) Ω Q ∇N f j † · ∇ (cid:0) N w j † − D w g † (cid:1) + (cid:10) j † , γ (cid:0) N w j † − D w g † (cid:1)(cid:11) = (cid:90) Ω Q ∇ (cid:0) N ξ j † − N f j † (cid:1) · ∇ (cid:0) N w j † − D w g † (cid:1) . Since γ N ξ j † = γ N f j † = g † , it follows that N ξ j † − N f j † ∈ H (Ω). We thus obtain from the last inequality12 (cid:107) ξ − f ∗ (cid:107) L (Ω) − (cid:107) f − f ∗ (cid:107) L (Ω) ≥ (cid:90) Ω Q ∇N w j † · ∇ (cid:0) N ξ j † − N f j † (cid:1) − (cid:90) Ω Q ∇D w g † · ∇ (cid:0) N ξ j † − N f j † (cid:1) = (cid:0) w, N ξ j † − N f j † (cid:1) + (cid:10) j † , γ (cid:0) N ξ j † − N f j † (cid:1)(cid:11) − (cid:0) w, N ξ j † − N f j † (cid:1) = 0 , which finishes the proof. Remark 3.6.
Indeed, the statement of Proposition 3.5 is a special case of the general assertion that sourceconditions f − f ∗ = T ( v ) , v ∈ X, for self-adjoint non-negative bounded linear operators T in the Hilbertspace X can only hold if f is a f ∗ -minimum-norm solution (see [37, Remark 6]). However, as shown abovewe have here N w j † − D w g † = T ( v ) with v = w − f † ∈ L (Ω).Our second alternative focus is on the classical Tikhonov regularization approach based on formula (2.12),in which κ ( j δ , g δ ) ∈ H (Ω). Then the regularized solution can be established as f ρ,δ = ( ˜T ∗ ˜T + ρ I ) − (cid:104) ˜T ∗ ( κ ( g δ , j δ )) + ρf ∗ (cid:105) (3.6)(cf., e.g., [18, Sec. 5.1]). From (3.6) and Lemma 2.2 we easily derive the error estimate (cid:107) f ρ,δ − f † (cid:107) L (Ω) ≤ ρ (cid:107) ( ˜T ∗ ˜T + ρ I ) − ( f † − f ∗ ) (cid:107) L (Ω) + ˜ Kδ √ ρ ≤ (cid:107) f † − f ∗ (cid:107) L (Ω) + ˜ Kδ √ ρ , (3.7)on which the following proposition is based. Proposition 3.7.
For any a priori parameter choice ρ ( δ ) of the regularization parameter satisfying theconditions ρ ( δ ) → and δ (cid:112) ρ ( δ ) → as δ → we have for a sequence δ n → , associated data j δ n , g δ n , and associated Tikhonv-regularized solutions f ρ ( δ n ) ,δ n that lim n →∞ (cid:107) f ρ ( δ n ) ,δ n − f † (cid:107) L (Ω) = 0 , i.e. the regularized solutions are strongly convergent in L (Ω) to the f ∗ -minimum norm solution f † . Furthermore, under the source condition (3.3), which is equivalent to f † − f ∗ = [ ˜T ∗ ˜T ]( v ) , (3.9)we find even the bounds (cid:107) f ρ,δ − f † (cid:107) L (Ω) ≤ ρ (cid:107) ( ˜T ∗ ˜T + ρ I ) − [ ˜T ∗ ˜T ]( v ) (cid:107) L (Ω) + ˜ Kδ √ ρ ≤ ρ (cid:107) v (cid:107) L (Ω) + ˜ Kδ √ ρ , (3.10)which as a consequence of T ( v ) = [ ˜T ∗ ˜T ]( v ) for all v ∈ L (Ω) immediately yields the rate assertion of thefollowing proposition. 10 roposition 3.8. If there is a source element v ∈ L (Ω) such that the range-type source condition (3.3)(equivalent to (3.5) due to Corollary 3.3) is satisfied, then we have for an a priori choice ρ ( δ ) ∼ δ / of theregularization parameter ρ the convergence rate (cid:107) f ρ ( δ ) ,δ − f † (cid:107) L (Ω) = O ( δ / ) as δ → of Tikhonov-regularized solutions. If the regularization parameter is chosen as ρ ( δ ) ∼ δ , then we obtain alsohere the convergence rate (3.4) as in Proposition 3.2. Remark 3.9.
It was shown by a saturation theorem of Groetsch 1984 (see, e.g, [18, Proposition 4.20]) that(3.11) is the best possible rate for classical linear Tikhonov regularization, with the exception of singularcases. At first glance, it is amazing that the best possible rate (3.11) in Proposition 3.8 is higher than thebest possible rate (3.4) in Proposition 3.2. However, here the classical Tikhonov regularization makes useof the higher smoothness assumption κ ( j δ , g δ ) ∈ H (Ω), whereas the version of Lavrentiev regularizationemployed here ignores this higher smoothness of the right-hand side of equation (2.10) and considers κ ( j δ , g δ )only as an element in L (Ω), see the two different noise inequalities in Lemma 2.2. We in Section 3 applied the Lavrentiev regularization and the Tikhonov regularization for the continuousidentification problem. In the remain Sections 4 – 6 we will analyze the problem in finite dimensional spaces.So far we have not yet found investigations on the discretization error in a combination of both functionalsfor the fully setting, a fact which motivated the research presented in this paper.Let (cid:0) T h (cid:1) Since N hf j δ and D hf g δ are both in V h , so is f , provided that f ∗ ∈ V h . Thus, taking this intoaccount, a discretization of the set L (Ω) can be avoided. Proof of Theorem 4.2. One can see easily that the problem (cid:16) P hρ,δ (cid:17) has a unique solution. It remains to show(4.5). Let f ∈ L (Ω) be the minimizer to (cid:16) P hρ,δ (cid:17) . The first-order optimality condition yields that Υ hρ,δ (cid:48) ( f ) ξ = J hδ (cid:48) ( f ) ξ + 2 ρ ( ξ, f − f ∗ ) = 0 for all ξ ∈ L (Ω). A short computation shows J hδ (cid:48) ( f )( ξ ) = 2 (cid:16) ξ, N hf j δ − D hf g δ (cid:17) and so obtain (cid:16) ξ, ρ (cid:16) N hf j δ − D hf g δ (cid:17) + f − f ∗ (cid:17) = 0 for all ξ ∈ L (Ω), which finishes the proof.From now on C is a generic positive constant which is independent of the mesh size h of T h , the noise level δ and the regularization parameter ρ . Before presenting the convergence of finite element approximationswe here state some auxiliary results. Lemma 4.4. A projection operator Π h (cid:5) : L (Ω) → V h , (cid:5) exists such that Π h (cid:5) ϕ h = ϕ h for all ϕ h ∈ V h , (cid:5) and Π h (cid:5) (cid:0) H (Ω) (cid:1) ⊂ V h , ⊂ V h , (cid:5) . Furthermore, it satisfies the properties lim h → (cid:13)(cid:13) ϑ − Π h (cid:5) ϑ (cid:13)(cid:13) H (Ω) = 0 for all ϑ ∈ H (cid:5) (Ω) (4.6) and (cid:13)(cid:13) ϑ − Π h (cid:5) ϑ (cid:13)(cid:13) H (Ω) ≤ Ch (cid:107) ϑ (cid:107) H (Ω) for all ϑ ∈ H (cid:5) (Ω) ∩ H (Ω) . (4.7) Proof. Let Π h : L (Ω) → V h be the Clement’s mollification interpolation operator, see [17] and somegeneralizations [10, 11, 41]. We then define the operatorΠ h (cid:5) ϑ := Π h ϑ − | ∂ Ω | (cid:90) ∂ Ω γ Π h ϑ ∈ V h , (cid:5) , ∀ ϑ ∈ L (Ω)which has the properties (4.6) and (4.7). The proof is completed.On the basis of (4.6) and (4.7) we introduce for each Φ ∈ H (cid:5) (Ω) (cid:37) h Φ := (cid:13)(cid:13) Φ − Π h (cid:5) Φ (cid:13)(cid:13) H (Ω) . (4.8)We note that lim h → (cid:37) h Φ = 0 and 0 ≤ (cid:37) h Φ ≤ Ch (4.9)in case Φ ∈ H (Ω). Furthermore, let ( f, j, g ) ∈ L (Ω) × H − / ( ∂ Ω) × H / (cid:5) ( ∂ Ω) be fixed, we denote by α hf,j = (cid:13)(cid:13) N hf j − N f j (cid:13)(cid:13) H (Ω) and β hf,g = (cid:13)(cid:13) D hf g − D f g (cid:13)(cid:13) H (Ω) . (4.10)Then lim h → α hf,j = lim h → β hf,g = 0 . In particular, if N f j ∈ H (Ω) and D f g ∈ H (Ω), the error estimates α hf,j ≤ Ch and β hf,g ≤ Ch (4.11)are satisfied (cf. [13, 16]). 12 emma 4.5. Let ( f , j , g ) and ( f , j , g ) be arbitrary in L (Ω) × H − / ( ∂ Ω) × H / (cid:5) ( ∂ Ω) . Then theestimates (cid:13)(cid:13) N hf j − N hf j (cid:13)(cid:13) H (Ω) ≤ C N (cid:16) (cid:107) f − f (cid:107) L (Ω) + (cid:107) j − j (cid:107) H − / ( ∂ Ω) (cid:17) (4.12) and (cid:13)(cid:13) D hf g − D hf g (cid:13)(cid:13) H (Ω) ≤ C D (cid:16) (cid:107) f − f (cid:107) L (Ω) + (cid:107) g − g (cid:107) H / ( ∂ Ω) (cid:17) (4.13) hold for all h ≥ .Proof. According to the definition of the discrete Neumann operator, we have for all ϕ h ∈ V h , (cid:5) that (cid:82) Ω Q ∇N hf i j i · ∇ ϕ h = (cid:10) j i , γϕ h (cid:11) + (cid:0) f i , ϕ h (cid:1) with i = 1 , 2. Thus, Φ h N := N hf j − N hf j is the unique solu-tion to the variational problem (cid:82) Ω Q ∇ Φ h N · ∇ ϕ h = (cid:10) j − j , γϕ h (cid:11) + (cid:0) f − f , ϕ h (cid:1) for all ϕ h ∈ V h , (cid:5) and sothat (4.12) is satisfied. Likewise, we also obtain (4.13). The proof is completed. Lemma 4.6. Let (cid:0) T h n (cid:1) be a sequence of triangulations with lim n →∞ h n = 0 . Assume that ( j δ n , g δ n ) is asequence in H − / ( ∂ Ω) × H / (cid:5) ( ∂ Ω) convergent to ( j δ , g δ ) in the H − / ( ∂ Ω) × H / ( ∂ Ω) -norm and ( f n ) isa sequence in L (Ω) weakly convergent in L (Ω) to f , then there holds the inequality lim inf n →∞ J h n δ n ( f n ) ≥ J δ ( f ) . (4.14) Proof. The proof is based upon the mollification operator introduced in Lemma 4.4 together with standardarguments, therefore omitted here.We now show the convergence of finite element approximations to the identification problem. Theorem 4.7. Assume that lim n →∞ h n = 0 and ( δ n ) and ( ρ n ) any positive sequences such that ρ n → , δ n √ ρ n → , α h n f † ,j † √ ρ n → and β h n f † ,g † √ ρ n → as n → ∞ , (4.15) where α h n f † ,j † and β h n f † ,g † are defined by (4.10) . Furthermore, assume that ( j δ n , g δ n ) is a sequence in H − / ( ∂ Ω) × H / (cid:5) ( ∂ Ω) satisfying (cid:13)(cid:13) j δ n − j † (cid:13)(cid:13) H − / ( ∂ Ω) + (cid:13)(cid:13) g δ n − g † (cid:13)(cid:13) H / ( ∂ Ω) ≤ δ n and f n := f h n ρ n ,δ n is the unique minimizer of (cid:16) P h n ρ n ,δ n (cid:17) for each n ∈ N . Then:(i) The sequence ( f n ) converges in the L (Ω) -norm to f † .(ii) The corresponding state sequences (cid:16) N h n f n j δ n (cid:17) and (cid:16) D h n f n g δ n (cid:17) converge in the H (Ω) -norm to the uniqueweak solution Φ † = Φ † (cid:0) f † , j † , g † (cid:1) of the boundary value problem (1.1) – (1.3) . Before going to prove the theorem, we make the following short remark. Remark 4.8. In case the weak solution Φ † = Φ † (cid:0) f † , j † , g † (cid:1) of (1.1)–(1.3) belonging to H (Ω), the estimate(4.11) shows that 0 ≤ α h n f † ,j † , β h n f † ,g † ≤ Ch n . Therefore, in view of (4.15), the above convergences (i) and (ii)are obtained if the sequence ( ρ n ) is chosen such that ρ n → , δ n √ ρ n → h n √ ρ n → n → ∞ . By regularity theory for elliptic boundary value problems, the regularity assumption Φ † ∈ H (Ω) is satisfiedif the diffusion matrix Q ∈ C , (Ω) d × d , j † ∈ H / ( ∂ Ω), g † ∈ H / ( ∂ Ω) and either ∂ Ω is smooth of the class C , or the domain Ω is convex (see, for example, [20, 44]).13 roof of Theorem 4.7. We have from the optimality of f n that J h n δ n ( f n ) + ρ n (cid:107) f n − f ∗ (cid:107) L (Ω) ≤ J h n δ n ( f † ) + ρ n (cid:107) f † − f ∗ (cid:107) L (Ω) . (4.16)Since at f † there holds the equation N f † j † = D f † g † , we infer from Lemma 4.5 that J h n δ n ( f † ) ≤ C (cid:18) δ n + (cid:16) α h n f † ,j † (cid:17) + (cid:16) β h n f † ,g † (cid:17) (cid:19) (4.17)which implies from (4.16) that lim n →∞ J h n δ n ( f n ) = 0 and, by the assumption (4.15),lim sup n →∞ (cid:107) f n − f ∗ (cid:107) L (Ω) ≤ (cid:13)(cid:13) f † − f ∗ (cid:13)(cid:13) L (Ω) . (4.18)So that the sequence ( f n ) is bounded in the L (Ω)-norm. A subsequence not relabelled and an element (cid:98) f ∈ L (Ω) exist such that ( f n ) converges weakly in L (Ω) to (cid:98) f and (cid:13)(cid:13) (cid:98) f − f ∗ (cid:13)(cid:13) L (Ω) ≤ lim inf n →∞ (cid:107) f n − f ∗ (cid:107) L (Ω) . (4.19)For any f ∈ L (Ω) we denote by J ( f ) := (cid:82) Ω Q ∇ (cid:0) N f j † − D f g † (cid:1) · ∇ (cid:0) N f j † − D f g † (cid:1) . By (1.7), we have (cid:13)(cid:13)(cid:13) N (cid:98) f j † − D (cid:98) f g † (cid:13)(cid:13)(cid:13) H (Ω) ≤ C Ω C Ω q J (cid:0) (cid:98) f (cid:1) ≤ C Ω C Ω q lim inf n →∞ J h n δ n ( f n ) = 0, here we used Lemma 4.6. Thus, N (cid:98) f j † = D (cid:98) f g † which infers (cid:98) f ∈ I (cid:0) j † , g † (cid:1) . Now we show (cid:98) f = f † and the sequence ( f n ) converges to (cid:98) f in the L (Ω)-norm. By the definition of the f ∗ -minimum-norm solution and (4.18)–(4.19), we get that (cid:13)(cid:13) f † − f ∗ (cid:13)(cid:13) L (Ω) ≤ (cid:13)(cid:13) (cid:98) f − f ∗ (cid:13)(cid:13) L (Ω) ≤ lim inf n →∞ (cid:107) f n − f ∗ (cid:107) L (Ω) ≤ lim sup n →∞ (cid:107) f n − f ∗ (cid:107) L (Ω) ≤ (cid:13)(cid:13) f † − f ∗ (cid:13)(cid:13) L (Ω) and so that (cid:13)(cid:13) f † − f ∗ (cid:13)(cid:13) L (Ω) = (cid:13)(cid:13) (cid:98) f − f ∗ (cid:13)(cid:13) L (Ω) = lim n →∞ (cid:107) f n − f ∗ (cid:107) L (Ω) . By the uniqueness of the minimum-norm solution and the sequence ( f n ) weakly converging in L (Ω) to (cid:98) f , we conclude that (cid:98) f = f † and thesequence ( f n ) in fact converges in the L (Ω)-norm to (cid:98) f . Finally, we show the sequences ( N h n f n j δ n ) and( D h n f n g δ n ) converge to Φ † = N f † j † = D f † g † in the H (Ω)-norm. Indeed, by Lemma 4.5, we obtain that (cid:13)(cid:13)(cid:13) N h n f n j δ n − N f † j † (cid:13)(cid:13)(cid:13) H (Ω) ≤ (cid:13)(cid:13)(cid:13) N h n f n j δ n − N h n f † j † (cid:13)(cid:13)(cid:13) H (Ω) + (cid:13)(cid:13)(cid:13) N h n f † j † − N f † j † (cid:13)(cid:13)(cid:13) H (Ω) ≤ C (cid:16)(cid:13)(cid:13) j δ n − j † (cid:13)(cid:13) H − / ( ∂ Ω) + (cid:13)(cid:13) f n − f † (cid:13)(cid:13) L (Ω) + α h n f † ,j † (cid:17) → n → ∞ . Similarly, we also get (cid:13)(cid:13)(cid:13) D h n f n g δ n − D f † g † (cid:13)(cid:13)(cid:13) H (Ω) ≤ C (cid:16)(cid:13)(cid:13) g δ n − g † (cid:13)(cid:13) H / ( ∂ Ω) + (cid:13)(cid:13) f n − f † (cid:13)(cid:13) L (Ω) + β h n f † ,g † (cid:17) → n tends to ∞ , which finishes the proof. We are now in a position to state the main theorem on convergence rates for the general case of finite elementdiscretized regularized solutions with noise level ( δ > h > Theorem 5.1. Assume that the condition (3.5) is fulfilled. Then, we have the error estimate and convergencerate (cid:13)(cid:13)(cid:13) N hf h j δ − D hf h g δ (cid:13)(cid:13)(cid:13) H (Ω) + ρ (cid:13)(cid:13) f h − f † (cid:13)(cid:13) L (Ω) = O (cid:18) δ + (cid:16) α hf † ,j † (cid:17) + (cid:16) β hf † ,g † (cid:17) + ρ(cid:37) h N w j † + ρ(cid:37) h N f † j † + ρ(cid:37) h D γ N w j † − g † + δρ + ρ (cid:19) , (5.1) where f h := f hρ,δ is the unique minimizer of (cid:16) P hρ,δ (cid:17) and D γ N w j † − g † is the unique weak solution to theDirichlet problem −∇ · ( Q ∇ v ) = 0 in Ω and v = γ N w j † − g † on ∂ Ω and α hf † ,j † , β hf † ,g † , (cid:37) h N w j † , (cid:37) h N f † j † and (cid:37) h D γ N w j † − g † come from (4.8) and (4.10) . emark 5.2. In case (cf. Remark 4.8) N f † j † , N w j † , D γ N w j † − g † ∈ H (Ω) , by (4.9) and (4.11), we have0 ≤ α hf † ,j † , β hf † ,g † , (cid:37) h N w j † , (cid:37) h N f † j † , (cid:37) h D γ N w j † − g † ≤ Ch and so that the following convergence rate is obtained (cid:13)(cid:13)(cid:13) N hf h j δ − D hf h g δ (cid:13)(cid:13)(cid:13) H (Ω) + ρ (cid:13)(cid:13) f h − f † (cid:13)(cid:13) L (Ω) = O (cid:0) δ + h + hρ + δρ + ρ (cid:1) . Remark 5.3. Let Φ † = Φ † (cid:0) f † , j † , g † (cid:1) be the weak solution of (1.1)–(1.3). Then the convergence rate (cid:13)(cid:13)(cid:13) N hf h j δ − Φ † (cid:13)(cid:13)(cid:13) H (Ω) + (cid:13)(cid:13)(cid:13) D hf h g δ − Φ † (cid:13)(cid:13)(cid:13) H (Ω) = O (cid:18) δ ρ − + (cid:16) α hf † ,j † (cid:17) ρ − + (cid:16) β hf † ,g † (cid:17) ρ − + (cid:37) h N w j † + (cid:37) h N f † j † + (cid:37) h D γ N w j † − g † + δ + ρ + α hf † ,j † + β hf † ,g † (cid:19) is also established. Indeed, the desired equation directly follows from (5.1) and the following inequalities (cid:13)(cid:13)(cid:13) N hf h j δ − N f † j † (cid:13)(cid:13)(cid:13) H (Ω) ≤ C (cid:16)(cid:13)(cid:13) j δ − j † (cid:13)(cid:13) H − / ( ∂ Ω) + (cid:13)(cid:13) f h − f † (cid:13)(cid:13) L (Ω) + α hf † ,j † (cid:17) ≤ C (cid:16) δ + (cid:13)(cid:13) f h − f † (cid:13)(cid:13) L (Ω) + α hf † ,j † (cid:17) and (cid:13)(cid:13)(cid:13) D hf h g δ − D f † g † (cid:13)(cid:13)(cid:13) H (Ω) ≤ C (cid:16) δ + (cid:13)(cid:13) f h − f † (cid:13)(cid:13) L (Ω) + β hf † ,g † (cid:17) , here we used Lemma 4.5. Proof of Theorem 5.1. In view of (4.17) we first have that J hδ (cid:0) f † (cid:1) ≤ C (cid:18) δ + (cid:16) α hf † ,j † (cid:17) + (cid:16) β hf † ,g † (cid:17) (cid:19) . Theoptimality of f h yields J hδ (cid:0) f h (cid:1) + ρ (cid:13)(cid:13) f h − f ∗ (cid:13)(cid:13) L (Ω) ≤ J hδ (cid:0) f † (cid:1) + ρ (cid:13)(cid:13) f † − f ∗ (cid:13)(cid:13) L (Ω) . This gives J hδ (cid:0) f h (cid:1) + ρ (cid:13)(cid:13) f h − f † (cid:13)(cid:13) L (Ω) ≤ J hδ (cid:0) f † (cid:1) + ρ (cid:16)(cid:13)(cid:13) f † − f ∗ (cid:13)(cid:13) L (Ω) − (cid:13)(cid:13) f h − f ∗ (cid:13)(cid:13) L (Ω) + (cid:13)(cid:13) f h − f † (cid:13)(cid:13) L (Ω) (cid:17) ≤ C (cid:18) δ + (cid:16) α hf † ,j † (cid:17) + (cid:16) β hf † ,g † (cid:17) (cid:19) + 2 ρ (cid:0) f † − f ∗ , f † − f h (cid:1) . (5.2)Since N f † j † = D f † g † , it follows that (cid:0) f † − f ∗ , f † − f h (cid:1) = (cid:0) f † − f h , N w j † − N f † j † (cid:1) + (cid:0) f † − f h , D f † g † − D w g † (cid:1) . (5.3)From (1.10), we infer (cid:0) f † , N w j † − N f † j † (cid:1) = (cid:90) Ω Q ∇N f † j † · ∇ (cid:0) N w j † − N f † j † (cid:1) − (cid:10) j † , γ (cid:0) N w j † − N f † j † (cid:1)(cid:11) , and (cid:0) f h , N w j † − N f † j † (cid:1) = (cid:90) Ω Q ∇N f h j † · ∇ (cid:0) N w j † − N f † j † (cid:1) − (cid:10) j † , γ (cid:0) N w j † − N f † j † (cid:1)(cid:11) . This in turn implies (cid:0) f † − f h , N w j † − N f † j † (cid:1) = (cid:90) Ω Q ∇ (cid:0) N f † j † − N f h j † (cid:1) · ∇ (cid:0) N w j † − N f † j † (cid:1) = (cid:90) Ω Q ∇ (cid:0) N f † j † − D f h g † (cid:1) · ∇ (cid:0) N w j † − N f † j † (cid:1) + (cid:90) Ω Q ∇ (cid:0) D f h g † − N f h j † (cid:1) · ∇ (cid:0) N w j † − N f † j † (cid:1) . (5.4)Since γ (cid:0) D f † g † − D w g † (cid:1) = 0, it follows from (1.11) that (cid:0) f † − f h , D f † g † − D w g † (cid:1) = (cid:90) Ω Q ∇ (cid:0) D f † g † − D f h g † (cid:1) · ∇ (cid:0) D f † g † − D w g † (cid:1) (5.5)15olds. We thus infer from (5.3)–(5.5) the identity (cid:0) f † − f ∗ , f † − f h (cid:1) = (cid:90) Ω Q ∇ (cid:0) D f h g † − N f h j † (cid:1) · ∇ (cid:0) N w j † − N f † j † (cid:1) + (cid:90) Ω Q ∇ (cid:0) N f † j † − D f h g † (cid:1) · ∇ (cid:0) N w j † − N f † j † (cid:1) + (cid:90) Ω Q ∇ (cid:0) D f † g † − D f h g † (cid:1) · ∇ (cid:0) D f † g † − D w g † (cid:1) . (5.6)We note again that N f † j † = D f † g † and γ (cid:0) D f † g † − D f h g † (cid:1) = 0. Then, together with (1.10) and (1.11), thelast two terms on the right hand side of (5.6) satisfy (cid:90) Ω Q ∇ (cid:0) N f † j † − D f h g † (cid:1) · ∇ (cid:0) N w j † − N f † j † (cid:1) + (cid:90) Ω Q ∇ (cid:0) D f † g † − D f h g † (cid:1) · ∇ (cid:0) D f † g † − D w g † (cid:1) = 0 . Thus, we obtain from (5.6) (cid:0) f † − f ∗ , f † − f h (cid:1) = (cid:90) Ω Q ∇ (cid:0) D f h g † − N f h j † (cid:1) · ∇ (cid:0) N w j † − N f † j † (cid:1) . Next, we abbreviate W = N w j † − N f † j † and note γW = γ N w j † − g † . (5.7)Then we get (cid:0) f † − f ∗ , f † − f h (cid:1) = (cid:90) Ω Q ∇ (cid:16) D f h g † − D hf h g † (cid:17) · ∇ W − (cid:90) Ω Q ∇ (cid:16) N f h j † − N hf h j † (cid:17) · ∇ W + (cid:90) Ω Q ∇ (cid:16) D hf h g † − D hf h g δ (cid:17) · ∇ W − (cid:90) Ω Q ∇ (cid:16) N hf h j † − N hf h j δ (cid:17) · ∇ W + (cid:90) Ω Q ∇ (cid:16) D hf h g δ − N hf h j δ (cid:17) · ∇ W := I + I + I . (5.8)To prepare the estimation of those three addends we start with writing (cid:90) Ω Q ∇ (cid:16) D f h g † − D hf h g † (cid:17) · ∇ W = (cid:90) Ω Q ∇ (cid:16) D f h g † − D hf h g † (cid:17) · ∇D γW + (cid:90) Ω Q ∇ (cid:16) D f h g † − D hf h g † (cid:17) · ∇ Π h (cid:5) ( W − D γW )+ (cid:90) Ω Q ∇ (cid:16) D f h g † − D hf h g † (cid:17) · ∇ (cid:0) W − D γW − Π h (cid:5) ( W − D γW ) (cid:1) . Since D f h g † − D hf h g † ∈ H (Ω), we then get (cid:90) Ω Q ∇ (cid:16) D f h g † − D hf h g † (cid:17) · ∇D γW = (cid:90) Ω Q ∇D γW · ∇ (cid:16) D f h g † − D hf h g † (cid:17) = 0 . Since γ ( W − D γW ) = γW − γ D γW = γW − γW = 0, we infer Π h (cid:5) ( W − D γW ) ∈ V h , = V h ∩ H (Ω)and then obtain from (1.11) and (4.3) that (cid:82) Ω Q ∇ (cid:16) D f h g † − D hf h g † (cid:17) · ∇ Π h (cid:5) ( W − D γW ) = 0 holds. Hencewe have (cid:12)(cid:12)(cid:12)(cid:12) (cid:90) Ω Q ∇ (cid:16) D f h g † − D hf h g † (cid:17) · ∇ W (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω Q ∇ (cid:16) D f h g † − D hf h g † (cid:17) · ∇ (cid:0) W − D γW − Π h (cid:5) ( W − D γW ) (cid:1)(cid:12)(cid:12)(cid:12)(cid:12) ≤ C (cid:16)(cid:13)(cid:13) f h (cid:13)(cid:13) L (Ω) + (cid:13)(cid:13) g † (cid:13)(cid:13) H / ( ∂ Ω) (cid:17) (cid:16) (cid:37) h N w j † + (cid:37) h N f † j † + (cid:37) h D γ N w j † − g † (cid:17) , (5.9)where we use (5.7). Similarly, since Π h (cid:5) W ∈ V h (cid:5) and by (1.10) and (4.1), we get (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω Q ∇ (cid:16) N f h j † − N hf h j † (cid:17) · ∇ W (cid:12)(cid:12)(cid:12)(cid:12) ≤ C (cid:16)(cid:13)(cid:13) f h (cid:13)(cid:13) L (Ω) + (cid:13)(cid:13) j † (cid:13)(cid:13) H − / ( ∂ Ω) (cid:17) (cid:16) (cid:37) h N w j † + (cid:37) h N f † j † (cid:17) . (5.10)16ow we are in the position to estimate I − I . Combining (5.9) with (5.10), we obtain ρ | I | ≤ C (cid:18) δ + (cid:16) α hf † ,j † (cid:17) + (cid:16) β hf † ,g † (cid:17) + ρ(cid:37) h N w j † + ρ(cid:37) h N f † j † + ρ(cid:37) h D γ N w j † − g † (cid:19) . (5.11)Now, using Lemma 4.5, we arrive at ρ | I | ≤ Cρ (cid:16)(cid:13)(cid:13) g δ − g † (cid:13)(cid:13) H / ( ∂ Ω) + (cid:13)(cid:13) j δ − j † (cid:13)(cid:13) H − / ( ∂ Ω) (cid:17) ≤ Cδρ. (5.12)Since for a.e. in Ω the matrix Q ( x ) is positive definite, the root Q ( x ) / is then well defined. Thus, usingthe Cauchy-Schwarz inequality and Young’s inequality, we estimate I as ρ | I | ≤ Cρ (cid:0) J hδ (cid:0) f h (cid:1)(cid:1) / ≤ C ρ + 14 J hδ (cid:0) f h (cid:1) ≤ Cρ + 14 J hδ (cid:0) f h (cid:1) . (5.13)It follows from (5.8) and (5.11)–(5.13) that2 ρ (cid:0) f † − f ∗ , f † − f h (cid:1) ≤ C (cid:18) δ + (cid:16) α hf † ,j † (cid:17) + (cid:16) β hf † ,g † (cid:17) + ρ(cid:37) h N w j † + ρ(cid:37) h N f † j † + ρ(cid:37) h D γ N w j † − g † + ρδ + ρ (cid:19) + 12 J hδ (cid:0) f h (cid:1) holds, which together with (5.2) implies12 J hδ (cid:0) f h (cid:1) + ρ (cid:13)(cid:13) f h − f † (cid:13)(cid:13) L (Ω) ≤ C (cid:18) δ + (cid:16) α hf † ,j † (cid:17) + (cid:16) β hf † ,g † (cid:17) + ρ(cid:37) h N w j † + ρ(cid:37) h N f † j † + ρ(cid:37) h D γ N w j † − g † + ρδ + ρ (cid:19) . (5.14)Since (cid:13)(cid:13)(cid:13) D hf h g δ − N hf h j δ (cid:13)(cid:13)(cid:13) H (Ω) ≤ C J hδ (cid:0) f h (cid:1) , (5.1) now directly follows from (5.14), which finishes the proof. In this section we will utilize the conjugate gradient (CG) method (see, for example, [21, 26]) to find theminimizes of the strictly convex, discrete regularized problem (cid:16) P hρ,δ (cid:17) . Let ∇ Υ hρ,δ ( f ) = 2 (cid:16) N hf j δ − D hf g δ (cid:17) +2 ρ ( f − f ∗ ) be the L -gradient of the cost function Υ hρ,δ at f (see Proof of Theorem 4.2), where f ∗ ∈ V h .Then the sequence of iterates via this algorithm is generated by f ∈ L (Ω) ∩ V h and f k +1 := f k + t k d k for k ≥ 0, where d k := (cid:40) −∇ Υ hρ,δ ( f k ) if k = 0 , −∇ Υ hρ,δ ( f k ) + β k d k − if k > β k := (cid:107)∇ Υ hρ,δ ( f k ) (cid:107) (cid:107)∇ Υ hρ,δ ( f k − ) (cid:107) and t k := arg min t ≥ Υ hρ,δ ( f k + td k ) . A short computation shows that t k = − (cid:82) Ω Q ∇ (cid:0) N hd k − D hd k (cid:1) · ∇ (cid:16) N hf k j δ − D hf k g δ (cid:17) + ρ (cid:0) d k , f k − f ∗ (cid:1)(cid:82) Ω Q ∇ (cid:0) N hd k − D hd k (cid:1) · ∇ (cid:0) N hd k − D hd k (cid:1) + ρ (cid:107) d k (cid:107) L (Ω) = − (cid:16) d k , ∇ Υ hρ,δ ( f k ) (cid:17)(cid:0) d k , N hd k − D hd k (cid:1) + ρ (cid:107) d k (cid:107) L (Ω) . Consequently, the CG method then reads as follows: giving an initial approximation f ∈ V h , number ofiterations N and a positive constants τ , τ . Computing ∇ Υ hρ,δ ( f ) = 2 (cid:0) N hf j δ − D hf g δ (cid:1) + 2 ρ ( f − f ∗ ) , d = −∇ Υ hρ,δ ( f ) , t = 12 (cid:13)(cid:13) d (cid:13)(cid:13) L (Ω) (cid:0) d , N hd − D hd (cid:1) + ρ (cid:107) d (cid:107) L (Ω) and setting f = f + t d and k = 1 , Tolerance := (cid:13)(cid:13) ∇ Υ hρ,δ (cid:0) f k (cid:1)(cid:13)(cid:13) L (Ω) − τ − τ (cid:13)(cid:13) ∇ Υ hρ,δ (cid:0) f (cid:1)(cid:13)(cid:13) L (Ω) . hile ( Tolerance > 0) & ( k ≤ N ) do r = (cid:13)(cid:13) ∇ Υ hρ,δ ( f k − ) (cid:13)(cid:13) L (Ω) , r = (cid:13)(cid:13) ∇ Υ hρ,δ ( f k ) (cid:13)(cid:13) L (Ω) , β k = rr ,d k = −∇ Υ hρ,δ ( f k ) + β k d k − , t k = − (cid:16) d k , ∇ Υ hρ,δ ( f k ) (cid:17)(cid:0) d k , N hd k − D hd k (cid:1) + ρ (cid:107) d k (cid:107) L (Ω) ,f k +1 = f k + t k d k ,k := k + 1 , Tolerance := (cid:13)(cid:13) ∇ Υ hρ,δ (cid:0) f k (cid:1)(cid:13)(cid:13) L (Ω) − τ − τ (cid:13)(cid:13) ∇ Υ hρ,δ (cid:0) f (cid:1)(cid:13)(cid:13) L (Ω) . end Algorithm 1: CG iterationBelow we illustrate the theoretical result with numerical examples. For this purpose we consider the theboundary value problem −∇ · (cid:0) Q ∇ Φ (cid:1) = f † in Ω := ( − , × ( − , , (6.1) Q ∇ Φ · (cid:126)n = j † on ∂ Ω and Φ = g † on ∂ Ω . (6.2)We assume that entries of the known symmetric diffusion matrix Q are discontinuous which are defined as q = 3 χ Ω + χ Ω \ Ω , q = χ Ω , q = 4 χ Ω + 2 χ Ω \ Ω , where χ D is the characteristic function of theLebesgue measurable set D andΩ := (cid:8) ( x , x ) ∈ Ω (cid:12)(cid:12) | x | ≤ / | x | ≤ / (cid:9) , Ω := (cid:8) ( x , x ) ∈ Ω (cid:12)(cid:12) | x | + | x | ≤ / (cid:9) andΩ := (cid:8) ( x , x ) ∈ Ω (cid:12)(cid:12) x + x ≤ / (cid:9) . The identified source function f † ∈ L (Ω) in (6.1) is assumed to be discontinuous and defined as f † = 2 χ Ω − χ Ω + 5 π π − χ Ω \ (Ω ∪ Ω ) , where Ω := (cid:8) ( x , x ) ∈ Ω (cid:12)(cid:12) x + 1 / + 16( x − / ≤ (cid:9) andΩ := (cid:8) ( x , x ) ∈ Ω (cid:12)(cid:12) ( x − / + ( x + 1 / ≤ / (cid:9) . For the discretization we divide the interval ( − , 1) into (cid:96) equal segments and so that the domain Ω = ( − , is divided into 2 (cid:96) triangles, where the diameter of each triangle is h (cid:96) = √ (cid:96) . In the minimization problem (cid:16) P hρ,δ (cid:17) we take h = h (cid:96) and ρ = ρ (cid:96) = 0 . h (cid:96) . We use Algorithm 1 which is described above for computingthe numerical solution of the problem (cid:16) P h (cid:96) ρ (cid:96) ,δ (cid:96) (cid:17) . As an a-priori estimate and the initial approximation wechoose f ∗ := 0 and f ( x ) := χ (0 , × [ − , − χ [ − , × [ − , . Example 6.1. In this first example j † ∈ H − / ( ∂ Ω) is chosen to be the piecewise constant function definedby j † = χ (0 , ×{− } − χ [ − , ×{ } + 2 χ (0 , ×{ } − χ [ − , ×{− } + 3 χ {− }× ( − , − χ { }× (0 , + 4 χ { }× ( − , − χ {− }× (0 , . (6.3)Then g † ∈ H / (cid:5) ( ∂ Ω) is defined as g † = γ N f † j † . We mention that, to avoid a so-called inverse crime, wegenerate the given data on a finer grid than those used in the computations. For this purpose we firstsolve the problem (6.1) supplemented with the Neumann boundary condition in (6.2) on the very fine gridwith (cid:96) = 128, and then use this numerical approximation as substitute for ( j † , g † ) in our computationalconsiderations below.For observations with noise we assume that( j δ (cid:96) , g δ (cid:96) ) = (cid:0) j † + θ (cid:96) · R j † , g † + θ (cid:96) · R g † (cid:1) for some θ (cid:96) > (cid:96), (6.4)18here R j † and R g † are ∂M h (cid:96) × − , 1) which are generated bythe MATLAB function “rand”, and ∂M h (cid:96) is the number of boundary nodes of the triangulation T h (cid:96) . Themeasurement error is then computed as δ (cid:96) = (cid:13)(cid:13) j δ (cid:96) − j † (cid:13)(cid:13) L ( ∂ Ω) + (cid:13)(cid:13) g δ (cid:96) − g † (cid:13)(cid:13) L ( ∂ Ω) . To satisfy the condition δ (cid:96) · ρ − / (cid:96) → (cid:96) → ∞ in Theorem 4.7 we below take θ (cid:96) = h (cid:96) √ ρ (cid:96) .We start with the coarsest level (cid:96) = 4. At each iteration k we computeTolerance := (cid:13)(cid:13) ∇ Υ h (cid:96) ρ (cid:96) ,δ (cid:96) (cid:0) f k(cid:96) (cid:1)(cid:13)(cid:13) L (Ω) − τ − τ (cid:13)(cid:13) ∇ Υ h (cid:96) ρ (cid:96) ,δ (cid:96) (cid:0) f (cid:96) (cid:1)(cid:13)(cid:13) L (Ω) , where τ := 10 − h / (cid:96) and τ := 10 − h / (cid:96) . Then the iteration was stopped if Tolerance ≤ (cid:96) = 4, we use its interpolation on the next finer mesh (cid:96) = 8as an initial approximation f for the algorithm on this finer mesh, and proceed similarly on the precedingrefinement levels.Let f (cid:96) be the function which is obtained at the final iterate of Algorithm 1 corresponding to the refinementlevel (cid:96) . Furthermore, let N h (cid:96) f (cid:96) j δ (cid:96) and D h (cid:96) f (cid:96) g δ (cid:96) denote the computed numerical solution to the Neumann andDirichlet problem −∇ · ( Q ∇ u ) = f (cid:96) in Ω and Q ∇ u · (cid:126)n = j δ (cid:96) on ∂ Ω and − ∇ · ( Q ∇ v ) = f (cid:96) in Ω and v = g (cid:96) on ∂ Ω , respectively. The notations N h (cid:96) f † j † and D h (cid:96) f † g † of the exact numerical solutions are to be understood similarly.We use the following abbreviations for the errors L f = (cid:13)(cid:13) f (cid:96) − f † (cid:13)(cid:13) L (Ω) , L N = (cid:13)(cid:13) N h (cid:96) f (cid:96) j δ (cid:96) − N h (cid:96) f † j † (cid:13)(cid:13) L (Ω) , H N = (cid:13)(cid:13) N h (cid:96) f (cid:96) j δ (cid:96) − N h (cid:96) f † j † (cid:13)(cid:13) H (Ω) and L D = (cid:13)(cid:13) D h (cid:96) f (cid:96) g δ (cid:96) − D h (cid:96) f † g † (cid:13)(cid:13) L (Ω) , H D = (cid:13)(cid:13) D h (cid:96) f (cid:96) g δ (cid:96) − D h (cid:96) f † g † (cid:13)(cid:13) H (Ω) . The numerical results are summarized in Table 1 and Table 2, where we present the refinement level (cid:96) , meshsize h (cid:96) of the triangulation, regularization parameter ρ (cid:96) , measured noise δ (cid:96) , number of iterations, value oftolerances and errors L f , L N , L D , H N and H D . Their experimental order of convergence (EOC) is presentedin Table 3, where EOC Θ := ln Θ( h ) − ln Θ( h )ln h − ln h and Θ( h ) is an error function of the mesh size h .All figures presented correspond to (cid:96) = 64. Figure 1 from left to right shows the computed numerical solution f (cid:96) of the algorithm at the final 579th-iteration, and the differences N h (cid:96) f † j † − N h (cid:96) f (cid:96) j δ (cid:96) , D h (cid:96) f † g † − D h (cid:96) f (cid:96) g δ (cid:96) and D h (cid:96) f (cid:96) g δ (cid:96) − N h (cid:96) f (cid:96) j δ (cid:96) . Convergence history (cid:96) h (cid:96) ρ (cid:96) δ (cid:96) Iterate Tolerance (cid:96) , mesh size h (cid:96) of the triangulation, regularization parameter ρ (cid:96) , measurementnoise δ (cid:96) , number of iterates and value of Tolerance. 19 onvergence history (cid:96) L f L N L D H N H D L f , L N , L D , H N and H D . Experimental order of convergence (cid:96) EOC L f EOC L N EOC L D EOC H N EOC H D Mean of EOC L f , L N , L D , H N and H D .Figure 1: Computed numerical solution f (cid:96) of the algorithm at the final iteration, and the differences N h (cid:96) f † j † −N h (cid:96) f (cid:96) j δ (cid:96) , D h (cid:96) f † g † − D h (cid:96) f (cid:96) g δ (cid:96) and D h (cid:96) f (cid:96) g δ (cid:96) − N h (cid:96) f (cid:96) j δ (cid:96) . Example 6.2. In present example we assume that multiple measurements are available, say (cid:0) j iδ , g iδ (cid:1) i =1 ,...,I .Then, problem (cid:16) P hρ,δ (cid:17) in Section 4 is given bymin f ∈ L (Ω) ¯Υ hρ,δ ( f ) := min f ∈ L (Ω) I I (cid:88) i =1 (cid:90) Ω Q ∇ (cid:0) N hf j iδ − D hf g iδ (cid:1) · ∇ (cid:0) N hf j iδ − D hf g iδ (cid:1)(cid:124) (cid:123)(cid:122) (cid:125) := ¯ J hδ ( q ) + ρ (cid:107) f − f ∗ (cid:107) L (Ω) (cid:16) ¯ P hρ,δ (cid:17) , which also attains a solution ¯ f hρ,δ . The Neumann boundary condition in the equation (6.2) is chosen in thesame form as (6.3), i.e. j † ( A,B,C,D ) = A · χ (0 , ×{− } − A · χ [ − , ×{ } + B · χ (0 , ×{ } − B · χ [ − , ×{− } + C · χ {− }× ( − , − C · χ { }× (0 , + D · χ { }× ( − , − D · χ {− }× (0 , , (6.5)and depends on the constants A, B, C and D . Let g † ( A,B,C,D ) := γ N f † j † ( A,B,C,D ) and assume that noisyobservations are given by (cid:16) j ( A,B,C,D ) δ (cid:96) , g ( A,B,C,D ) δ (cid:96) (cid:17) = (cid:16) j † ( A,B,C,D ) + θ · R j † ( A,B,C,D ) , g † ( A,B,C,D ) + θ · R g † ( A,B,C,D ) (cid:17) , (6.6)20here R j † ( A,B,C,D ) and R g † ( A,B,C,D ) denote ∂M h (cid:96) × − , θ appeared in the equation (6.6) is now independent of the grid level (cid:96) .In the case ( A, B, C, D ) = (1 , , , 4) we have a single noisy measurement couple, i.e. I = 1. We nowfix D = 4, and let ( A, B, C ) take all permutations S of the set { , , } . Then, the equations (6.5)–(6.6)generate I = 6 measurements. Similarly, if ( A, B, C, D ) takes all permutations S of { , , , } we get I = 16measurements. With θ = 0 . (cid:96) = 64 we compute the noise level¯ δ (cid:96) = (cid:13)(cid:13) j (1 , , , δ (cid:96) − j † (1 , , , (cid:13)(cid:13) L ( ∂ Ω) + (cid:13)(cid:13) g (1 , , , δ (cid:96) − g † (1 , , , (cid:13)(cid:13) L ( ∂ Ω) if ( A, B, C, D ) = (1 , , , , (cid:80) ( A,B,C ) ∈S (cid:13)(cid:13) j ( A,B,C, δ (cid:96) − j † ( A,B,C, (cid:13)(cid:13) L ( ∂ Ω) + (cid:13)(cid:13) g ( A,B,C, δ (cid:96) − g † ( A,B,C, (cid:13)(cid:13) L ( ∂ Ω) if D = 4 , (cid:80) ( A,B,C,D ) ∈S (cid:13)(cid:13) j ( A,B,C,D ) δ (cid:96) − j † ( A,B,C,D ) (cid:13)(cid:13) L ( ∂ Ω) + (cid:13)(cid:13) g ( A,B,C,D ) δ (cid:96) − g † ( A,B,C,D ) (cid:13)(cid:13) L ( ∂ Ω) . The corresponding numerical results for the multiple measurement case are presented in the Table 4. Numerical results for (cid:96) = 64 , θ = 0 . with multiple observations I Iterate Tolerance ¯ δ (cid:96) L f L N L D H N H D (cid:96) = 64, θ = 0 . 1, and with multiple measurements I = 1 , , I h (cid:96) f † of the exact source and the computednumerical solution q (cid:96) of the algorithm at the final iteration for (cid:96) = 64, θ = 0 . 1, and I = 16 , , 1, respectively.Figure 2: Interpolation I h (cid:96) f † , computed numerical solution f (cid:96) of the algorithm at the final iteration for (cid:96) = 64, θ = 0 . 1, and with multiple measurements I = 16 , , 1, respectively. Acknowledgments The authors thank the Referee and Editor for their valuable comments and suggestions.M. Hinze gratefully acknowledges support of the Lothar Collatz Center for Computing in Science at theUniversity of Hamburg.B. Hofmann gratefully acknowledges support by the German Research Foundation (DFG) under grantHO 1454/12-1.T. N. T. Quyen gratefully acknowledges support of the Alexander von Humboldt Foundation. References [1] S. Acosta, S. Chow, J. Taylor and V. Villamizar, On the multi–frequency inverse source problem inheterogeneous media, Inverse Problems Nonlinear Ill-posed Problems of Monotone Type , Dordrecht, Springer,2006.[3] C. J. S. Alves, N. F. M. Martins and N. C. Roberty, Full identification of acoustic sources with multiplefrequencies and boundary measurements, Inverse Probl. Imaging Ann. Math. InverseProblems Inverse Problems Estimation Techniques for Distributed Parameter Systems , Systems &Control: Foundations & Applications, Boston: Birkh¨auser, 1989.[8] G. Bao, J. Lin and F. Triki, A multi-frequency inverse source problem, J. Differential Equations Inverse Problems SIAM J. Numer. Anal. SIAM J. Numer. Anal. Inverse Problems The Mathematical Theory of Finite Element Methods, New York: Springer,2008.[14] A. P. Calder´on, On an inverse boundary value problem, in Seminar on Numerical Analysis and itsApplications to Continuum Physics (Rio de Janeiro, 1980), 65-73, Soc. Brasil. Mat., Rio de Janeiro,1980.[15] G. Chavent, Nonlinear Least Squares for Inverse Problems. Theoretical Foundations and Step-by-StepGuide for Applications , New York: Springer, 2009.[16] P. G. Ciarlet, Basis Error Estimates for Elliptic Problems, Handbook of Numerical Analisis, Vol. II, P.G. Ciarlet and J. -L. Lions, eds., North-Holland Amsterdam: Elsevier, 1991.[17] P. Cl´ement, Approximation by finite element functions using local regularization, RAIRO Anal. Num´er. Regularization of Inverse Problems , Dortdrecht: Kluwer, 1996.[19] A. Farcas, L. Elliott, DB. Ingham, D. Lesnic and NS. Mera, A dual reciprocity boundary elementmethod for the regularized numerical solution of the inverse source problem associated to the Poissonequation, Inverse Problems in Engineering Elliptic Problems in Nonsmooth Domains , Boston: Pitman, 1985.[21] W. W. Hager and H. Zhang, A new conjugate gradient method with guaranteed descent and an efficientline search, SIAM J. Optim. Rev. Mod.Phys. Comput. Optim. Appl. Inverse Probl. Imaging Inverse Source Problems , Rhode-Island: American Mathematical Society, 1989.[26] C. T. Kelley, Iterative Methods for Optimization , Philadelphia: SIAM, 1999.[27] R. V. Kohn and M. Vogelius, Determining conductivity by boundary measurements, Comm. Pure Appl.Math. Comm. Pure Appl. Math. Comm. Pure Appl. Math. SIAM J. Control Optim. Advances in Boundary Element Techniques IV , R. Gallego, MH. Aliabadi,eds., University of London: Queen Mary, 2003, 177–182.[32] T. Matsumoto, M. Tanaka M and T. Tsukamoto, Identifications of source distributions using BEM withdual reciprocity method, Inverse Problems in Engineering Mechanics IV , M. Tanaka, ed., Amsterdam,New York: Elsevier, 2003, 127–135.[33] A. I. Nachman, Reconstructions from boundary measurements, Ann. Math. Inverse and Ill-posed Problems (Sankt Wolfgang, 1986),volume 4 of Notes Rep. Math. Sci. Engrg. , pages 53–75. Academic Press, Boston, MA, 1987.[35] C. Pechstein, Finite and Boundary Element Tearing and Interconnecting Solvers for Multiscale Prob-lems , Heidelberg New York Dordrecht London: Springer, 2010.[36] R. Plato, Converse results, saturation and quasi-optimality for Lavrentiev regularization of accretiveproblems, SIAM J. Numer. Anal. Mathematics of Computation J. Math.Anal. Appl. SIAM J. Appl. Math Regularization Methods in BanachSpaces , Berlin Boston: Walter de Gruyter GmbH & Co. KG, 2012.[41] R. Scott and S. Y. Zhang, Finite element interpolation of nonsmooth function satisfying boundaryconditions, Math. Comp. Inverse Problem Theory and Methods for Model Parameter Estimation , Philadelphia:SIAM, 2005.[43] U. Tautenhahn, On the method of Lavrentiev regularization for nonlinear ill-posed problems, InverseProblems