A regularization approach for an inverse source problem in elliptic systems from single Cauchy data

TL;DR

This paper develops a Tikhonov-type regularization method for identifying sources in elliptic systems from single boundary measurements, providing convergence analysis, error bounds, and a numerical solution approach.

Contribution

It introduces a novel regularization approach for inverse elliptic problems with theoretical convergence and error analysis, supported by numerical experiments.

Findings

Convergence of finite element approximations to the inverse problem solutions.

Derivation of error bounds and convergence rates under source conditions.

Numerical validation using a conjugate gradient method.

Abstract

In this paper we investigate the problem of identifying the source term in an elliptic system from a single noisy measurement couple of the Neumann and Dirichlet data. A variational method of Tikhonov-type regularization with specific misfit term 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 in infinite dimensional Hilbert spaces, convergence and rate results can be found from the general theory of classical Tikhonov and Lavrentiev regularization. Using the variational discretization concept, where the PDE is discretized with piecewise linear and continuous finite elements, we show the convergence of finite element approximations to solutions of the identification problem. Moreover, we derive an error bound and corresponding convergence rates provided a suitable range-type source condition is satisfied. For the numerical solution we propose a conjugate gradient method. To illustrate the theoretical results, a numerical case study is presented which supports our analytical findings.

Figures

Tables

Table 1. Table 1: Refinement level ℓ ℓ \ell , mesh size h ℓ subscript ℎ ℓ h_{\ell} of the triangulation, regularization parameter ρ ℓ subscript 𝜌 ℓ \rho_{\ell} , measurement noise δ ℓ subscript 𝛿 ℓ \delta_{\ell} , number of iterates and value of Tolerance.

Convergence history
$\ell$	$h_{\ell }$	${\rho }_{\ell }$	${\delta }_{\ell }$	Iterate	Tolerance
4	0.7071	0.7071e-2	0.1916	312	-3.0822e-5
8	0.3536	0.3536e-2	9.3172e-2	387	-1.2739e-6
16	0.1766	0.1766e-2	4.1174e-2	461	-1.4029e-6
32	8.8388e-2	0.8839e-3	2.0932e-2	505	-1.8559e-7
64	4.4194e-2	0.4419e-3	7.2765e-3	579	-7.3540e-9

Table 2. Table 2: Errors L f 2 subscript superscript 𝐿 2 𝑓 L^{2}_{f} , L 𝒩 2 subscript superscript 𝐿 2 𝒩 L^{2}_{\mathcal{N}} , L 𝒟 2 subscript superscript 𝐿 2 𝒟 L^{2}_{\mathcal{D}} , H 𝒩 1 subscript superscript 𝐻 1 𝒩 H^{1}_{\mathcal{N}} and H 𝒟 1 subscript superscript 𝐻 1 𝒟 H^{1}_{\mathcal{D}} .

Convergence history
$\ell$	$L_f^2$	$L_{𝒩}^2$	$L_{𝒟}^2$	$H_{𝒩}^1$	$H_{𝒟}^1$
4	0.5215	2.0441e-2	2.0396e-2	6.9952e-2	6.9713e-2
8	0.3309	6.3175e-3	6.3083e-3	3.1374e-2	3.1311e-2
16	0.1915	2.0132e-3	2.0122e-3	1.7276e-2	1.7243e-2
32	0.1073	5.5434e-4	5.5426e-4	8.9136e-3	8.9130e-3
64	5.2568e-2	1.4669e-4	1.4666e-4	3.9352e-3	3.9347e-3

Table 3. Table 3: Experimental order of convergence between finest and coarsest level for L f 2 subscript superscript 𝐿 2 𝑓 L^{2}_{f} , L 𝒩 2 subscript superscript 𝐿 2 𝒩 L^{2}_{\mathcal{N}} , L 𝒟 2 subscript superscript 𝐿 2 𝒟 L^{2}_{\mathcal{D}} , H 𝒩 1 subscript superscript 𝐻 1 𝒩 H^{1}_{\mathcal{N}} and H 𝒟 1 subscript superscript 𝐻 1 𝒟 H^{1}_{\mathcal{D}} .

Experimental order of convergence
$\ell$	EOC${}_{L_f^2}$	EOC${}_{L_{𝒩}^2}$	EOC${}_{L_{𝒟}^2}$	EOC${}_{H_{𝒩}^1}$	EOC${}_{H_{𝒟}^1}$
4	–	–	–	–	–
8	0.6563	1.6940	1.6930	1.1568	1.1548
16	0.7891	1.6499	1.6485	0.8608	0.8607
32	0.8357	1.8606	1.8601	0.9547	0.9520
64	1.0294	1.9180	1.9181	1.1796	1.1797
Mean of EOC	0.8276	1.7806	1.7799	1.0380	1.0368

Table 4. Table 4: Numerical results for ℓ = 64 ℓ 64 \ell=64 , θ = 0.1 𝜃 0.1 \theta=0.1 , and with multiple measurements I = 1 , 6 , 16 𝐼 1 6 16 I=1,6,16 .

Numerical results for $\ell =64$, $\theta =0.1$ with multiple observations
$I$	Iterate	Tolerance	${\overline{\delta }}_{\ell }$	$L_f^2$	$L_{𝒩}^2$	$L_{𝒟}^2$	$H_{𝒩}^1$	$H_{𝒟}^1$
1	531	-3.2313e-8	0.3292	0.3280	5.9096e-3	5.9090e-3	0.1225	0.1221
6	517	-7.1620e-9	0.3331	0.2583	4.3125e-3	4.3122e-3	7.9322e-2	7.9320e-2
16	536	-6.4706e-8	0.3289	0.1747	2.8465e-3	2.8461e-3	5.2318e-2	5.2314e-2

Equations

$$-\nabla\cdot\big{(}Q\nabla\Phi\big{)}=f\mbox{~{}in~{}}\Omega\subset R^{d},d\in\{2,3\},\,Q\nabla\Phi\cdot\vec{n}=j\mbox{~{}on~{}}\partial\Omega\,\mbox{~{}and~{}}\Phi=g\mbox{~{}on~{}}\partial\Omega$$

$$\displaystyle-\nabla\cdot\big{(}Q\nabla\Phi\big{)}$$

$$Q\nabla \Phi \cdot \vec{n}$$

$$\Phi$$

$$Q(x)\xi \cdot \xi =\sum _{1\le r,s\le d}{q}_{rs}(x){\xi }_r{\xi }_s\ge \underline{q}|\xi {|}^2\text{mbox}a.e.in\Omega$$

$$\displaystyle\big{\|}j_{\delta}-j^{\dagger}\big{\|}_{H^{-1/2}(\partial\Omega)}+\big{\|}g_{\delta}-g^{\dagger}\big{\|}_{H^{1/2}(\partial\Omega)}\leq\delta.$$

$$H^{1}_{\diamond}(\Omega):=\left\{u\in H^{1}(\Omega)~{}\Big{|}~{}\int_{\partial\Omega}\gamma udx=0\right\}\mbox{~{}and~{}}H^{1/2}_{\diamond}(\partial\Omega):=\left\{g\in H^{1/2}(\partial\Omega)~{}\Big{|}~{}\int_{\partial\Omega}g(x)dx=0\right\}$$

$$C_{\Omega }{\int }_{\Omega }{\phi }^2\le {\int }_{\Omega }|\nabla \phi {|}^2\text{mbox}forall\phi \in H_{\diamond }^1(\Omega ).$$

$$\Vert \phi {\Vert }_{H^1(\Omega )}^2\le \frac{1+C_{\Omega }}{C_{\Omega }}{\int }_{\Omega }|\nabla \phi {|}^2\le \frac{1+C_{\Omega }}{C_{\Omega }\underline{q}}{\int }_{\Omega }Q\nabla \phi \cdot \nabla \phi$$

$$-\nabla \cdot (Q\nabla u)=f\text{mbox}in\Omega \text{mbox}andQ\nabla u\cdot \vec{n}=j\text{mbox}on\partial \Omega$$

$$-\nabla \cdot (Q\nabla v)=f\text{mbox}in\Omega \text{mbox}andv=g\text{mbox}on\partial \Omega .$$

$${\int }_{\Omega }Q\nabla u\cdot \nabla \phi =⟨j,\gamma \phi ⟩+\left(f,\phi \right)$$

$$\mathcal{N}:L^2(\Omega )\to H_{\diamond }^1(\Omega )\text{mbox}withf\mapsto {\mathcal{N}}_{f}j,$$

$${\int }_{\Omega }Q\nabla v\cdot \nabla \psi =\left(f,\psi \right)$$

$$\mathcal{D}:L^2(\Omega )\to H_{\diamond }^1(\Omega )\text{mbox}withf\mapsto {\mathcal{D}}_{f}g,$$

$${\Lambda }_f:H^{-1/2}(\partial \Omega )\to H_{\diamond }^{1/2}(\partial \Omega ),j\mapsto {\Lambda }_{f}j:=\gamma {\mathcal{N}}_{f}j.$$

$${\mathcal{N}}_{f}j={\mathcal{N}}_{f}0+{\mathcal{N}}_{0}j\text{mbox}and{\mathcal{D}}_{f}g={\mathcal{D}}_{f}0+{\mathcal{D}}_{0}g$$

$$f_1,f_2\in L^2(\Omega )\text{mbox}with{\Lambda }_{f_1}={\Lambda }_{f_2}\Rightarrow f_1=f_2?$$

$$\underset{f\in L^2(\Omega )}{\min }\Vert {\Lambda }_f-{\Lambda }_{\delta }{\Vert }_{\ast }^2+\rho \Vert f-f^{\ast }{\Vert }_{L^2(\Omega )}^2$$

$$Given\left(j^{†},g^{†}\right)\in H^{-1/2}(\partial \Omega )\times H_{\diamond }^{1/2}(\partial \Omega )with{\Lambda }_f{j}^{†}=g^{†},findf\in L^2(\Omega ).$$

$$\mathcal{I}\left(j^{\dagger},g^{\dagger}\right):=\left\{f\in L^{2}(\Omega)~{}\big{|}~{}\Lambda_{f}j^{\dagger}=g^{\dagger}\right\}=\left\{f\in L^{2}(\Omega)~{}\big{|}~{}\mathcal{N}_{f}j^{\dagger}=\mathcal{D}_{f}g^{\dagger}\right\}$$

$$\underset{f\in \mathcal{I}\left(j^{†},g^{†}\right)}{\min }\Vert f-f^{\ast }{\Vert }_{L^2(\Omega )}^2.$$

$$\Vert {\mathcal{N}}_f{j}_{\delta }-{\mathcal{D}}_f{g}_{\delta }{\Vert }_{L^2(\Omega )}^2+\rho \Vert f-f^{\ast }{\Vert }_{L^2(\Omega )}^2$$

$${\mathcal{J}}_{\delta }(f):={\int }_{\Omega }Q\nabla \left({\mathcal{N}}_f{j}_{\delta }-{\mathcal{D}}_f{g}_{\delta }\right)\cdot \nabla \left({\mathcal{N}}_f{j}_{\delta }-{\mathcal{D}}_f{g}_{\delta }\right)dx,$$

$$\underset{f\in L^2(\Omega )}{\min }{{\rm Y}}_{\rho ,\delta }(f)\text{mbox}with{{\rm Y}}_{\rho ,\delta }(f):={\mathcal{J}}_{\delta }(f)+\rho \Vert f-f^{\ast }{\Vert }_{L^2(\Omega )}^2,$$

$$\frac{1}{2}\nabla {{\rm Y}}_{\rho ,\delta }(f)={\mathcal{N}}_f{j}_{\delta }-{\mathcal{D}}_f{g}_{\delta }+\rho (f-f^{\ast })\text{mbox}forallf\in L^2(\Omega ).$$

$${\mathcal{J}}_0(\xi )$$

$${\mathbf{L}}_{j,g}(f):={\mathcal{N}}_{f}j-{\mathcal{D}}_{f}g+\rho (f-f^{\ast })=0,$$

$$f-f^{\ast }=-\frac{1}{\rho }\left({\mathcal{N}}_f{j}_{\delta }-{\mathcal{D}}_f{g}_{\delta }\right).$$

$$\underset{f\in L^2(\Omega )}{\min }{\int }_{\Omega }Q\nabla \left({\mathcal{N}}_f^h{j}_{\delta }-{\mathcal{D}}_f^h{g}_{\delta }\right)\cdot \nabla \left({\mathcal{N}}_f^h{j}_{\delta }-{\mathcal{D}}_f^h{g}_{\delta }\right)dx+\rho \Vert f-f^{\ast }{\Vert }_{L^2(\Omega )}^2.$$

$$f-f^{\ast }={\mathcal{N}}_w{j}^{†}-{\mathcal{D}}_w{g}^{†},$$

Full text

A regularization approach for an inverse source problem in elliptic systems from single Cauchy data

Michael Hinzea††Email:  michael.hinze@uni-hamburg.de,  bernd.hofmann@mathematik.tu-chemnitz.de,  quyen.tran@uni-goettingen.de, Bernd Hofmannb and Tran Nhan Tam Quyenc

aDepartment of Mathematics, University of Hamburg, 20146 Hamburg, Germany

bFaculty of Mathematics, Chemnitz University of Technology, 09107 Chemnitz, Germany

cInstitute 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

[TABLE]

from a single noisy measurement couple $\left(j_{\delta},g_{\delta}\right)$ of the Neumann and Dirichlet data $(j,g)$ with noise level $\delta>0$. In this context, the diffusion matrix $Q$ is given. A variational method of Tikhonov-type regularization with specific misfit term 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 in infinite dimensional Hilbert spaces, convergence and rate results can be found from the general theory of classical Tikhonov and Lavrentiev regularization. Using the variational discretization concept, where the PDE is discretized with piecewise linear and continuous finite elements, we show the convergence of finite element approximations to solutions of the regularized problem. Moreover, we derive an error bound and corresponding convergence rates provided a suitable range-type source condition is satisfied. For the numerical solution we propose a conjugate gradient method. To illustrate the theoretical results, a numerical case study is presented which supports our analytical 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: 35R25; 47A52; 35R30; 65J20; 65J22.

1 Introduction

Let $\Omega$ be an open, bounded and connected domain of $R^{d},~{}d\in\{2,3\},$ with Lipschitz boundary $\partial\Omega$. We consider the elliptic system

[TABLE]

where $\vec{n}$ is the unit outward normal on $\partial\Omega$ and the diffusion matrix $Q$ is given. Furthermore, we assume that $Q:=\left(q_{rs}\right)_{1\leq r,s\leq d}\in{L^{\infty}(\Omega)}^{d\times d}$ is symmetric and satisfies the uniformly ellipticity condition

[TABLE]

for all $\xi=\left(\xi_{r}\right)_{1\leq r\leq d}\in R^{d}$ with some constant $\underline{q}>0$.

The system (1.1)–(1.3) is overdetermined, i.e. if the Neumann and Dirichlet boundary conditions $j^{\dagger}\in H^{-1/2}(\partial\Omega):={H^{1/2}(\partial\Omega)}^{*},~{}g^{\dagger}\in H^{1/2}(\partial\Omega)$, and the source term $f\in L^{2}(\Omega)$ are given, then there may be no $\Phi$ satisfying this system. In this paper we assume that the system is consistent and our aim is to reconstruct a function $f\in L^{2}(\Omega)$ in the system (1.1)–(1.3) from a noisy measurement couple $\left(j_{\delta},g_{\delta}\right)\in H^{-1/2}(\partial\Omega)\times H^{1/2}(\partial\Omega)$ of the exact Neumann and Dirichlet data \big{(}j^{\dagger},g^{\dagger}\big{)}, where $\delta>0$ stands for the measurement error, i.e. we assume the noise model

[TABLE]

The source identification problem in PDEs arises in many branches of applied science such as electroencephalography, geophysical prospecting and pollutant detection, and attracted great attention from many scientists in the last 30 years or so. For surveys on this subject we may consult in [7, 15, 18, 22, 25, 42] and the references therein. Up to now, only a limited number of works was investigated the general source identification 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 the above 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 the problem have been obtained in [5, 6, 30, 39]. A survey of the problem of simultaneously identifying the source term and coefficients in elliptic systems from distributed observations can be found in [38], where further references can be found.

In the present paper, the general source identification problem in elliptic partial differential equations from a single noisy measurement couple of Neumann and Dirichlet data is studied. So far, we have not yet found investigations on the discretization analysis for this source recovery problem, a fact which also motivated the research presented in the paper. By using a suitable version of the Tikhonov-type regularization with some non-standard misfit term we could outline that the source distribution inside the physical domain $\Omega$ can be reconstructed from a finite number of observations on the boundary $\partial\Omega$, at least by numerical approximations. The specific regularization approach proves to be a version of Lavrentiev regularization with implicit forward operator. One of the main results of the paper is to show convergence of the finite element discretized Tikhonov-regularized solutions to a sought source function. Another main result is the interpretation of an occurring condition of solution smoothness as a range-type source condition of Lavrentiev’s regularization method. This allows us to establish error bounds and corresponding convergence rates for the regularized solutions.

To formulate precisely the problem, we first give some notations. Let us denote by $\gamma:H^{1}(\Omega)\to H^{1/2}(\partial\Omega)$ the continuous Dirichlet trace operator with $\gamma^{-1}:H^{1/2}(\partial\Omega)\to H^{1}(\Omega)$ its continuous right inverse operator, i.e. $(\gamma\circ\gamma^{-1})g=g$ for all $g\in H^{1/2}(\partial\Omega)$. We set

[TABLE]

and denote by $C_{\Omega}$ the positive constant appearing in the Poincaré-Friedrichs inequality (cf. [35])

[TABLE]

Since $H^{1}_{0}(\Omega):=\left\{u\in H^{1}(\Omega)~{}|~{}\gamma u=0\right\}\subset H^{1}_{\diamond}(\Omega)$, the inequality (1.6) is in particular valid for all $\varphi\in H^{1}_{0}(\Omega)$. Furthermore, by (1.4), the coercivity condition

[TABLE]

holds for all $\varphi\in H^{1}_{\diamond}(\Omega)$.

Now, for any fixed $\left(j,g\right)\in H^{-1/2}(\partial\Omega)\times H^{1/2}_{\diamond}(\partial\Omega)$ we can simultaneously consider the Neumann problem

[TABLE]

as well as the Dirichlet problem

[TABLE]

By the aid of (1.7) and the Riesz representation theorem, we conclude that for each $f\in L^{2}(\Omega)$ there exists a unique weak solution $u$ of the problem (1.8) in the sense that $u\in H^{1}_{\diamond}(\Omega)$ and satisfies the identity

[TABLE]

for all $\varphi\in H^{1}_{\diamond}(\Omega)$, where notation $\left\langle j,g\right\rangle$ stands for the value of the function $j\in H^{-1/2}(\partial\Omega)$ at $g\in H^{1/2}(\partial\Omega)$ and the notation $\left(f,\varphi\right)$ is the inner product of $f$ and $\varphi$ in the space $L^{2}(\Omega)$. Then we can define the Neumann operator

[TABLE]

which maps each $f\in L^{2}(\Omega)$ to the unique weak solution $\mathcal{N}_{f}j:=u$ of the problem (1.8). Similarly, the problem (1.9) also attains a unique weak solution $v$ in the sense that $v\in H^{1}(\Omega)$, $\gamma v=g$ and the identity

[TABLE]

holds for all $\psi\in H^{1}_{0}(\Omega)$. The Dirichlet operator is defined as

[TABLE]

which maps each $f\in L^{2}(\Omega)$ to the unique weak solution $\mathcal{D}_{f}g:=v$ of the problem (1.9). Therefore, for any fixed $f\in L^{2}(\Omega)$ we can define the so-called Neumann-to-Dirichlet map

[TABLE]

We mention that since $H^{1}_{0}(\Omega)\subset H^{1}_{\diamond}(\Omega)$, we from (1.10) have that $\int_{\Omega}Q\nabla\mathcal{N}_{f}j\cdot\nabla\psi=\left(f,\psi\right)$ for all $\psi\in H^{1}_{0}(\Omega)$. In view of (1.11) we therefore conclude $\Lambda_{f}j=g\mbox{~{}if and only if~{}}\mathcal{N}_{f}j=\mathcal{D}_{f}g$, where the identities

[TABLE]

are satisfied, and the operators $f\mapsto\mathcal{N}_{f}0$ and $f\mapsto\mathcal{D}_{f}0$ are linear and bounded from $L^{2}(\Omega)$ into itself. Furthermore, $\Lambda_{f}j=\gamma\mathcal{N}_{f}j=\gamma\mathcal{N}_{0}j+\gamma\mathcal{N}_{f}0=\Lambda_{0}j+\Lambda_{f}0$, where $\Lambda_{0}j$ is linear, self-adjoint, bounded and invertible, as the diffusion $Q$ is smooth enough (cf. [33]).

As in electrical impedance tomography (EIT) or for the Calderón’s problem [4, 14, 33] one can pose the question whether the source distribution $f$ inside a physical domain $\Omega$ can be determined from an infinite number of observations on the boundary $\partial\Omega$, i.e. from the Neumann-to-Dirichlet map $\Lambda_{f}$:

[TABLE]

To the best of our knowledge, the above question is still open so far. In case an observation $\Lambda_{\delta}$ of $\Lambda_{f}$ being available one can use a certain regularization method to approximate the sought source. For example, one can consider for operator norms $\|\cdot\|_{*}$ a minimizer of the problem

[TABLE]

as a reconstruction along the lines of Tikhonov’s regularization method, where $\rho>0$ is the regularization parameter and $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 identified source, at least by numerical approximations. Furthermore, for simplicity of exposition we below restrict ourselves to the case of just one observation pair $\left(j_{\delta},g_{\delta}\right)$ being available, while the approach described here can be easily extended to multiple measurements $\left(j_{\delta}^{i},g_{\delta}^{i}\right)_{i=1,\ldots,I}$, see Section 6, Ex. 6.2. The inverse problem is thus stated as follows.

[TABLE]

In other words, the interested problem is, for given $\left(j^{\dagger},g^{\dagger}\right)\in H^{-1/2}(\partial\Omega)\times H^{1/2}_{\diamond}(\partial\Omega)$, to find some $f\in L^{2}(\Omega)$ and consequently $\Phi\in H^{1}_{\diamond}(\Omega)$ such that the system (1.1)–(1.3) is satisfied in the weak sense. Precisely, we define the general solution set

[TABLE]

of the inverse problem $\left(\mathcal{IP}\right)$. The source identification problem as described here is well known to be not uniquely determined from boundary observations (see a counterexample in [3]), i.e. the set $\mathcal{I}\left(j^{\dagger},g^{\dagger}\right)$ fails to be a singleton. Since not the Neumann-to-Dirichlet map is given, but only one pair $\left(j^{\dagger},g^{\dagger}\right)$, the problem is even highly underdetermined. Thus instead we will search for the uniquely determined $f^{*}$-minimum-norm solution $f^{\dagger}$, which is the minimizer of the problem

[TABLE]

As a consequence of item (iii) of Lemma 2.1 below, the set $\mathcal{I}\left(j^{\dagger},g^{\dagger}\right)$ is non-empty, closed and convex, hence $f^{\dagger}$ is uniquely determined. On the other hand, for all $f\in\mathcal{I}\left(j^{\dagger},g^{\dagger}\right)$ the equation $\mathcal{N}_{f}j^{\dagger}=\mathcal{D}_{f}g^{\dagger}$ is fulfilled. However, we have to solve this equation with noise data $\left(j_{\delta},g_{\delta}\right)\in H^{-1/2}(\partial\Omega)\times H^{1/2}_{\diamond}(\partial\Omega)$ of $\left(j^{\dagger},g^{\dagger}\right)$ satisfying (1.5). The simplest variety of regularization may be to consider a minimizer of the Tikhonov functional

[TABLE]

over $f\in L^{2}(\Omega)$ as an approximation solution to $f^{\dagger}$.

In present work we adopt the variational approach of Kohn and Vogelius [27, 28, 29] in using cost functional containing the gradient of forward operators to the above mentioned inverse source problem. More precisely, we use the convex functional

[TABLE]

instead of the mapping $f\mapsto\|\mathcal{N}_{f}j_{\delta}-\mathcal{D}_{f}g_{\delta}\|^{2}_{L^{2}(\Omega)}$, together with Tikhonov regularization and consider the unique solution $f_{\rho,\delta}$ of the strictly convex minimization problem

[TABLE]

where the gradient of the functional $\Upsilon_{\rho,\delta}$ can be explicitly written as

[TABLE]

The motivation in using this cost functional $\mathcal{J}_{\delta}$ as misfit functional is that for all $\xi\in L^{2}(\Omega)$ the inequality

[TABLE]

holds true and $\mathcal{J}_{0}(f)=0$ at any $f\in\mathcal{I}\left(j^{\dagger},g^{\dagger}\right)$. The advantage is evident, because the minimizer $f_{\rho,\delta}\in L^{2}(\Omega)$ satisfies the equation $\nabla\Upsilon_{\rho,\delta}(f_{\rho,\delta})=0$ such that, for $j:=j_{\delta}$, $g:=g_{\delta}$ and $f:=f_{\rho,\delta}$, we have

[TABLE]

and hence, for $f:=f_{\rho,\delta}$, we have

[TABLE]

Due to formula (1.18), the Tikhonov regularization approach under consideration with specific misfit term also 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 inverse problem under consideration, we apply in Section 3 the general theory of classical Tikhonov and Lavrentiev regularization for such problems yielding propositions on convergence and convergence rates for the regularized solutions in the infinite dimensional Hilbert spaces. However, for convenience in numerical analysis with the finite element methods introduced in Section 4 our focus is here on the extremal problem for minimizing the Tikhonov functional with Kohn-Vogelius misfit term and quadratic penalty. The use of different convex penalty terms, e.g. total variation, may be a work for us in future.

Let $\mathcal{N}^{h}_{f}j_{\delta}$ and $\mathcal{D}^{h}_{f}g_{\delta}$ be corresponding approximations of the solution maps $\mathcal{N}_{f}j_{\delta}$ and $\mathcal{D}_{f}g_{\delta}$ in the finite dimensional space $\mathcal{V}^{h}_{1}$ of piecewise linear, continuous finite elements. We then consider the discrete regularized problem corresponding to $\left(\mathcal{P}_{\rho,\delta}\right)$, i.e., the following strictly convex minimization problem

[TABLE]

Using the variational discretization concept introduced in [23], we show in Section 4 that the unique solution $f^{h}_{\rho,\delta}$ of the problem $\left(\mathcal{P}^{h}_{\rho,\delta}\right)$ automatically belongs to the finite dimensional space $\mathcal{V}^{h}_{1}$. Thus, a discretization of the admissible set $L^{2}(\Omega)$ can be avoided.

As $h,\delta\to 0$ and with an appropriate a-priori regularization parameter choice $\rho=\rho(h,\delta)$, also in Section 4, we prove that the sequence \big{(}f^{h}_{\rho,\delta}\big{)} converges to $f^{\dagger}$ in the $L^{2}(\Omega)$-norm. Furthermore, the corresponding state sequences \Big{(}\mathcal{N}^{h}_{f^{h}_{\rho,\delta}}j_{\delta}\Big{)} and \Big{(}\mathcal{D}^{h}_{f^{h}_{\rho,\delta}}g_{\delta}\Big{)} converge in the $H^{1}(\Omega)$-norm to $\Phi^{\dagger}=\Phi^{\dagger}(f^{\dagger},j^{\dagger},g^{\dagger})$ 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\in\mathcal{I}\left(j^{\dagger},g^{\dagger}\right)$ and there is a function $w\in L^{2}(\Omega)$ such that $f-f^{*}=\mathbf{L}_{j^{\dagger},g^{\dagger}}(w)$, or in other notation

[TABLE]

then $f=f^{\dagger}$, 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 for corresponding convergence rates of the continuous setting in infinite dimensional spaces as well as after incorporating the discretization. In the latter case, precisely for the known matrix $Q\in{C^{0,1}(\Omega)}^{d\times d}$ and the exact data $\left(j^{\dagger},g^{\dagger}\right)\in H^{1/2}(\partial\Omega)\times H^{3/2}(\partial\Omega)$, we derive the convergence rates

[TABLE]

and

[TABLE]

Finally, for the numerical solution of the discrete regularized problem $\left(\mathcal{P}^{h}_{\rho,\delta}\right)$ we employ in Section 6 a conjugate gradient algorithm. Numerical case studies illustrate the analytical results and show the efficiency of our theoretical findings.

We conclude this introduction with a remark that since the main interest is to clearly state our ideas, we only treat the model elliptic problem (1.1) while the approach described here can be easily extended to more general models, e.g., for the source identification problem in diffusion-reaction equations

[TABLE]

from a measurement $\left(j_{\delta},g_{\delta}\right)$ of \big{(}j^{\dagger},g^{\dagger}\big{)}, where $Q$ satisfying the condition (1.4), $0\neq\kappa=\kappa(x)\in L^{\infty}(\Omega)$, i.e the set $\{x\in\Omega|\kappa(x)\neq 0\}$ has positive Lebesgue measure, and $\sigma=\sigma(x)\in L^{\infty}(\partial\Omega)$ with $\sigma\geq 0$ are given. The variational approach is now formulated as the minimizing problem with the misfit

[TABLE]

over $f\in L^{2}(\Omega)$, 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, the problem of identifying the source term $f$ in the Helmholtz-type equation

[TABLE]

from measured Cauchy data $\left(j_{\delta}(\kappa),g_{\delta}(\kappa)\right)$ which is available for all frequency $\kappa>0$. The uniqueness results for this identification problem can be found in [3, 8], while several effective recovered algorithms have been presented in [1, 9].

Throughout the paper we use the standard notion of Sobolev spaces $H^{1}(\Omega)$, $H^{1}_{0}(\Omega)$, $W^{k,p}(\Omega)$, etc from, for example, [44]. If not stated otherwise we write $\int_{\Omega}\cdots$ instead of $\int_{\Omega}\cdots dx$.

2 Preliminaries and operator-theoretic settings

In order to define appropriate operators, we recall the decompositions (1.12) of the corresponding Neumann and Dirichlet problems, where $\mathcal{N}_{f}0$ and $\mathcal{D}_{f}0$ characterize linear mappings of $f\in L^{2}(\Omega)$. On the other hand, $\mathcal{N}_{0}j$ and $\mathcal{D}_{0}g$ depend nonlinearly on $j$ and $g$, respectively, but both are independent of $f$. Hence, the difference of $\mathcal{N}_{f}j$ and $\mathcal{D}_{f}g$ characterizes, for fixed elements $j$ and $g$, the affine mapping $\mathbf{L}_{j,g}$ of $f\in L^{2}(\Omega)$ defined by formula (1.17). First we introduce the linear operator $\mathbf{\tilde{T}}:L^{2}(\Omega)\to H^{1}_{\diamond}(\Omega)$ defined as

[TABLE]

Since the image elements $\mathcal{N}_{f}0-\mathcal{D}_{f}0\in H^{1}_{\diamond}(\Omega)$ also belong to $L^{2}(\Omega)$, one can moreover introduce the operator $\mathbf{T}:L^{2}(\Omega)\to L^{2}(\Omega)$ defined by

[TABLE]

where $\mathbf{T}(f)=\mathbf{\tilde{T}}(f)$ for all $f\in L^{2}(\Omega)$. On the other hand, we remark that the expression

[TABLE]

generates an inner product on the space $H^{1}_{\diamond}(\Omega)$ which is equivalent to the usual one. Now let

[TABLE]

be the adjoint operator of $\mathbf{\tilde{T}}:L^{2}(\Omega)\to H^{1}_{\diamond}(\Omega)$, where $H^{1}_{\diamond}(\Omega)$ is equipped with the inner product (2.2) above. For all $f\in L^{2}(\Omega)$ and $\phi\in H^{1}_{\diamond}(\Omega)$ we thus have

[TABLE]

by (1.10). We now decompose $H^{1}_{\diamond}(\Omega)$ into the orthogonal direct sum $H^{1}_{\diamond}(\Omega)=H^{1}_{0}(\Omega)\oplus{H^{1}_{0}(\Omega)}^{\perp}$ with respect to the inner product (2.2). We note for all $g\in H^{1/2}_{\diamond}(\partial\Omega)$ that $\mathcal{D}_{0}g\in{H^{1}_{0}(\Omega)}^{\perp}$. Furthermore,

[TABLE]

which implies $\dim{H^{1}_{0}(\Omega)}^{\perp}\geq\dim H^{1/2}_{\diamond}(\partial\Omega)=\infty$. For all $f\in L^{2}(\Omega)$ we deduce from (1.11) and (2.3) that

[TABLE]

or in other words $\ker\mathbf{\tilde{T}}^{*}=H^{1}_{0}(\Omega)$. Furthermore, for all $\widehat{\phi}\in{H^{1}_{0}(\Omega)}^{\perp}$ we get $\int_{\Omega}Q\nabla\mathcal{D}_{f}0\cdot\nabla\widehat{\phi}=0$, since $\mathcal{D}_{f}0\in H^{1}_{0}(\Omega)$. Again, the equation (2.3) implies that

[TABLE]

Therefore, $\mathbf{\tilde{T}}^{*}_{|{H^{1}_{0}(\Omega)}^{\perp}}$ is the compact embedding ${H^{1}_{0}(\Omega)}^{\perp}\hookrightarrow L^{2}(\Omega)$ and $\mathbf{\tilde{T}}^{*}$ is the composition of the projector from $H^{1}_{\diamond}(\Omega)$ onto ${H^{1}_{0}(\Omega)}^{\perp}$ and the embedding operator from $H^{1}(\Omega)$ to $L^{2}(\Omega)$. Furthermore we have, for all $f\in L^{2}(\Omega)$,

[TABLE]

because ${\rm range}(\mathbf{\tilde{T}})$ is orthogonal to $\ker\mathbf{\tilde{T}}^{*}$ and hence $\mathbf{\tilde{T}}^{*}$ acts only as embedding operator.

Lemma 2.1**.**

(i) The operator $\mathbf{T}$ defined by formula (2.1) is linear, bounded, self-adjoint and non-negative, i.e. we have

[TABLE]

Moreover, $\mathbf{T}$ is compact and has an infinite dimensional range which is non-closed, i.e. we have $\overline{\text{range}(\mathbf{T})}\not=\text{range}(\mathbf{T})$.

(ii) For any fixed $(j,g)\in H^{-1/2}(\partial\Omega)\times H^{1/2}(\partial\Omega)$ the map $\mathbf{L}_{j,g}:L^{2}(\Omega)\to L^{2}(\Omega)$ defined by

[TABLE]

is affine linear, continuous and monotone, i.e. we have

[TABLE]

(iii) The solution set $\mathcal{I}$ (cf. (1.13)) is a closed affine subspace of the Hilbert space $L^{2}(\Omega)$.

Proof.

(i) It follows from (1.10) that

[TABLE]

and similarly,

[TABLE]

Using (1.10)–(1.11) again, we get

[TABLE]

and the self-adjoint property of $\mathbf{T}$ now follows directly from (2.7)-(2.9). We further have from (1.10)–(1.11) for all $f\in L^{2}(\Omega)$ that $\int_{\Omega}Q\nabla\mathcal{D}_{f}0\cdot\nabla(\mathcal{N}_{f}0-\mathcal{D}_{f}0)=0.$ Combining this with the identity $\left(\mathbf{T}(f),f\right)=\int_{\Omega}Q\nabla\mathcal{N}_{f}0\cdot\nabla(\mathcal{N}_{f}0-\mathcal{D}_{f}0)$ we arrive at $\left(\mathbf{T}(f),f\right)=\int_{\Omega}Q\nabla(\mathcal{N}_{f}0-\mathcal{D}_{f}0)\cdot\nabla(\mathcal{N}_{f}0-\mathcal{D}_{f}0)\geq 0$.

We now show that $\mathbf{T}$ is compact and has an infinite dimensional range. The operator $\mathbf{T}:L^{2}(\Omega)\to L^{2}(\Omega)$ as a composition of a bounded linear operator $\mathbf{\tilde{T}}$ and a compact embedding operator is compact. Next, we show that $\dim(\mathbf{\tilde{T}})=\infty$. For deriving a contradiction we assume that $\dim\text{range}(\mathbf{\tilde{T}})<\infty$. Then we can write $H^{1}_{\diamond}(\Omega)=\text{range}(\mathbf{\tilde{T}})\oplus{\text{range}(\mathbf{\tilde{T}})}^{\perp}$ with respect to the inner product (2.2). By (2.3), for all $\varphi\in{\text{range}(\mathbf{\tilde{T}})}^{\perp}$ we get $\int_{\Omega}Q\nabla\mathcal{D}_{f}0\cdot\nabla\varphi=(f,\varphi)$ holding for all $f\in L^{2}(\Omega)$ which implies that $\varphi\in H^{1}_{0}(\Omega)$. Therefore, {H^{1}_{0}(\Omega)}^{\perp}\subset{\big{(}{\text{range}(\mathbf{\tilde{T}})}^{\perp}\big{)}}^{\perp}=\overline{\text{range}(\mathbf{\tilde{T}})}=\text{range}(\mathbf{\tilde{T}}) and this yields the contradiction $\infty=\dim{H^{1}_{0}(\Omega)}^{\perp}\leq\dim\text{range}(\mathbf{\tilde{T}})<\infty$. Consequently, $\mathbf{T}:L^{2}(\Omega)\to L^{2}(\Omega)$ is compact operator and possesses an infinite dimensional range.

(ii) The inequality (2.6) follows directly from (2.5).

(iii) Since $\mathbf{T}$ is a bounded linear operator, the solution set can be written as

[TABLE]

with the Moore-Penrose pseudoinverse $\mathbf{T}^{\dagger}$ of $\mathbf{T}$. Then the nullspace $\ker(\mathbf{T})$ is a closed subspace and $\mathbf{T}^{\dagger}(\mathcal{D}_{0}g^{\dagger}-\mathcal{N}_{0}j^{\dagger})$ a well-defined element in $L^{2}(\Omega)$. Consequently, $\mathcal{I}(j^{\dagger},g^{\dagger})$ is a non-empty closed affine subspace and hence also a convex set in the Hilbert space $L^{2}(\Omega)$. ∎

Due to item (ii) of Lemma 2.1, we can reformulate the identification problem $\left(\mathcal{IP}\right)$ as an operator equation with linear bounded self-adjoint non-negative operator $\mathbf{T}$ mapping in $L^{2}(\Omega)$. Finding an element $f\in\mathcal{I}\left(j^{\dagger},g^{\dagger}\right)$ is then equivalent to solving the linear operator equation

[TABLE]

This makes the inverse problem explicit, but we have to take into account that instead of the exact right-hand side $\kappa(j^{\dagger},g^{\dagger})$ only noisy data $\kappa(j_{\delta},g_{\delta})$ satisfying (1.5) for $j_{\delta}$ and $g_{\delta}$ are available. As a consequence of item (i) of Lemma 2.1 we see that the equation (2.10) formulated in the Hilbert space $L^{2}(\Omega)$ is ill-posed of type II in the sense of Nashed (cf. [34]). Stable approximate (regularized) solutions $f_{\rho,\delta}$ to equation (2.10) satisfy with $f=f_{\rho,\delta}$ the auxiliary linear operator equation

[TABLE]

in $L^{2}(\Omega)$ with some regularization parameter $\rho>0$.

At this point we should recall and note that we have, for all $j_{\delta}$, $g_{\delta}$ under consideration and $f\in L^{2}(\Omega)$, and due to (1.10) – (1.11)

[TABLE]

which means that the elements $\kappa(j_{\delta},g_{\delta})$ and $\mathbf{\tilde{T}}^{*}[\kappa(j_{\delta},g_{\delta})]$ in $L^{2}(\Omega)$ coincide. Nevertheless, we have to distinguish the two cases $\kappa(j_{\delta},g_{\delta})\in H^{1}_{\diamond}(\Omega)$ and $\kappa(j_{\delta},g_{\delta})\in L^{2}(\Omega)$ in the following lemma.

Lemma 2.2**.**

Under the noise model (1.5) there is a constant $0\leq\tilde{K}<\infty$ independent of $\delta$ such that

[TABLE]

Moreover, there is also a constant $0\leq K<\infty$ independent of $\delta$ such that

[TABLE]

Proof.

By Lemma 4.5 below the existence of such constant $\tilde{K}$ in the first estimate of this lemma follows with the settings $h:=0$, $f_{1}=f_{2}:=0$ and $\tilde{K}:=\max\{\mathcal{C_{N}},\mathcal{C_{D}}\}$. Then we have under the noise model (1.5)

[TABLE]

By setting $K:=\tilde{K}\|\mathbf{\tilde{T}}\|$ the second estimate of the lemma gets established. This completes the proof. ∎

We conclude this section by mentioning that the functional $\mathcal{J}_{\delta}$ defined by (1.15) is convex and weakly sequentially lower semi-continuous. In fact, the above defined operator $\mathbf{\tilde{T}}:L^{2}(\Omega)\to H^{1}_{\diamond}(\Omega)$ is compact (see the proof of Lemma 2.1). Using the equivalent inner product (2.2), we therefore conclude that

[TABLE]

which shows the convexity and weak sequential lower semi-continuity of the functional $\mathcal{J}_{\delta}$ and is the basis for a classical Tikhonov regularization approach in the subsequent section. Moreover, we have as a basis for Lavrentiev regularization in the subsequent section

[TABLE]

where the constant is independent of $f$, and it is well-known that, due to the properties of $\mathbf{T}$ from Lemma 2.1, the unique minimizer $f_{\rho,\delta}$ 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 using Kohn and Vogelius’ function (1.15) we can use the least squares function (cf. (1.14)), and then the jointed minimization problem reads as

[TABLE]

where $\mathbf{T}(f)$ and $\kappa(j_{\delta},g_{\delta})$ were defined by (2.1) and (2.10), respectively. One can show easily that the problem attains a unique solution $\bar{f}$. Furthermore, for computation this minimizer $\bar{f}$ in practice one must derive the $L^{2}$-gradient $\nabla\Theta(\bar{f})$. Let $\xi\in L^{2}(\Omega)$ be arbitrary. We will compute briefly the differential $\Theta^{\prime}(\bar{f})\xi$ as follows.

We have that $\frac{1}{2}\Theta^{\prime}(\bar{f})\xi=\left(\mathbf{T}(\bar{f})-\kappa(j_{\delta},g_{\delta}),\mathbf{T}(\xi)\right)+\rho(\bar{f}-f^{*},\xi)$. Consider the adjoint problem

[TABLE]

(We here do not use the homogeneous Dirichlet boundary condition $\bar{\Phi}=0$ on $\partial\Omega$ instead of (2.16), because in general $\mathbf{T}(\xi)\notin H^{1}_{0}(\Omega)$.) Next, we decompose

[TABLE]

with respect to the inner product (2.2). Then we obtain that

[TABLE]

As a result, we arrive at $\frac{1}{2}\nabla\Theta(\bar{f})=\bar{\Phi}^{2}+\rho(\bar{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, we avoid computing numerically for the terms $\bar{\Phi}^{1}$ and $\bar{\Phi}^{2}$ in the decomposition (2.17) which seems to be still very difficult for us.

3 Lavrentiev regularization versus Tikhonov regularization for the continuous 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 $\kappa$ is considered as an element in $L^{2}(\Omega)$ and alternatively as an element in $H^{1}(\Omega)$.

First our focus is on the Lavrentiev regularization approach based on formula (2.13), in which $\kappa(j_{\delta},g_{\delta})\in L^{2}(\Omega)$. 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_{\rho,\delta}$ with respect to the uniquely determined $f^{*}$-minimizing solution $f^{\dagger}$ to problem ($\mathcal{IP}-MN$). Taking into account that Lemma 2.2 holds, we immediately derive (see, e.g., [12, Rem. 3.3]) the error estimate

[TABLE]

This immediately yields (see [24, section 2]) the following convergence assertion.

Proposition 3.1**.**

For any a priori parameter choice $\rho(\delta)$ of the regularization parameter satisfying the conditions

[TABLE]

we have for a sequence $\delta_{n}\to 0$, associated data $j_{\delta_{n}},g_{\delta_{n}}$, and associated Lavrentiev-regularized solutions $f_{\rho(\delta_{n}),\delta_{n}}$ that

[TABLE]

i.e. the regularized solutions are strongly convergent in $L^{2}(\Omega)$ to the $f^{*}$-minimum norm solution $f^{\dagger}$.

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\in L^{2}(\Omega)$ such that the range-type source condition

[TABLE]

is satisfied, we have for an a priori choice $\rho(\delta)\sim\sqrt{\delta}$ of the regularization parameter $\rho$ the convergence rate

[TABLE]

of Lavrentiev-regularized solutions.

Corollary 3.3**.**

The rate result (3.4) remains true if we have some element $w\in L^{2}(\Omega)$ such that

[TABLE]

Proof.

Evidently we have $\mathcal{N}_{w}j^{\dagger}-\mathcal{D}_{w}g^{\dagger}=\mathbf{T}(w)-\kappa(j^{\dagger},g^{\dagger})$ and $f^{\dagger}$ is a solution to the operator equation (2.10) with exact right-hand side $\kappa(j^{\dagger},g^{\dagger})$. Hence (3.5) can be rewritten as $f^{\dagger}-f^{*}=\mathbf{T}(w-f^{\dagger})$. This yields (3.3) with the new source element $v:=w-f^{\dagger}$. ∎

Remark 3.4**.**

It was shown by a saturation theorem in [36] that (3.4) is the best possible rate for linear Lavrentiev regularization, with the exception of singular situations with respect to the forward operator, here $\mathbf{T}$, and with respect to the solution $f^{\dagger}$.

Revisiting the source condition (3.5), we add the following remarkable result.

Proposition 3.5**.**

Assume that $f\in\mathcal{I}\left(j^{\dagger},g^{\dagger}\right)$ and that there is a function $w\in L^{2}(\Omega)$ such that $f-f^{*}=\mathcal{N}_{w}j^{\dagger}-\mathcal{D}_{w}g^{\dagger}$. Then $f$ is the uniquely determined $f^{*}$-minimum-norm solution of problem $\left(\mathcal{IP}-MN\right)$, i.e. we have $f=f^{\dagger}$.

Proof.

We have with \xi\in\left\{\xi\in L^{2}(\Omega)~{}\big{|}~{}\mathcal{N}_{\xi}j^{\dagger}=\mathcal{D}_{\xi}g^{\dagger}\right\} that

[TABLE]

Since $\gamma\mathcal{N}_{\xi}j^{\dagger}=\gamma\mathcal{N}_{f}j^{\dagger}=g^{\dagger}$, it follows that $\mathcal{N}_{\xi}j^{\dagger}-\mathcal{N}_{f}j^{\dagger}\in H^{1}_{0}(\Omega)$. We thus obtain from the last inequality

[TABLE]

which finishes the proof. ∎

Remark 3.6**.**

Indeed, the statement of Proposition 3.5 is a special case of the general assertion that source conditions $f-f^{*}=\mathbf{T}(v),\;v\in X,$ for self-adjoint non-negative bounded linear operators $\mathbf{T}$ in the Hilbert space $X$ can only hold if $f$ is a $f^{*}$-minimum-norm solution (see [37, Remark 6]). However, as shown above we have here $\mathcal{N}_{w}j^{\dagger}-\mathcal{D}_{w}g^{\dagger}=\mathbf{T}(v)$ with $v=w-f^{\dagger}\in L^{2}(\Omega)$.

Our second alternative focus is on the classical Tikhonov regularization approach based on formula (2.12), in which $\kappa(j_{\delta},g_{\delta})\in H^{1}(\Omega)$. Then the regularized solution can be established as

[TABLE]

(cf., e.g., [18, Sec. 5.1]). From (3.6) and Lemma 2.2 we easily derive the error estimate

[TABLE]

on which the following proposition is based.

Proposition 3.7**.**

For any a priori parameter choice $\rho(\delta)$ of the regularization parameter satisfying the conditions

[TABLE]

we have for a sequence $\delta_{n}\to 0$, associated data $j_{\delta_{n}},g_{\delta_{n}}$, and associated Tikhonv-regularized solutions $f_{\rho(\delta_{n}),\delta_{n}}$ that

[TABLE]

i.e. the regularized solutions are strongly convergent in $L^{2}(\Omega)$ to the $f^{*}$-minimum norm solution $f^{\dagger}$.

Furthermore, under the source condition (3.3), which is equivalent to

[TABLE]

we find even the bounds

[TABLE]

which as a consequence of $\mathbf{T}(v)=[\mathbf{\tilde{T}}^{*}\mathbf{\tilde{T}}](v)$ for all $v\in L^{2}(\Omega)$ immediately yields the rate assertion of the following proposition.

Proposition 3.8**.**

If there is a source element $v\in L^{2}(\Omega)$ 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 $\rho(\delta)\sim\delta^{2/3}$ of the regularization parameter $\rho$ the convergence rate

[TABLE]

of Tikhonov-regularized solutions. If the regularization parameter is chosen as $\rho(\delta)\sim\delta$, then we obtain also here 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 singular cases. At first glance, it is amazing that the best possible rate (3.11) in Proposition 3.8 is higher than the best possible rate (3.4) in Proposition 3.2. However, here the classical Tikhonov regularization makes use of the higher smoothness assumption $\kappa(j_{\delta},g_{\delta})\in H^{1}(\Omega)$, whereas the version of Lavrentiev regularization employed here ignores this higher smoothness of the right-hand side of equation (2.10) and considers $\kappa(j_{\delta},g_{\delta})$ only as an element in $L^{2}(\Omega)$, see the two different noise inequalities in Lemma 2.2.

4 Finite element approximation and convergence for the discretized problem

We in Section 3 applied the Lavrentiev regularization and the Tikhonov regularization for the continuous identification 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 functionals for the fully setting, a fact which motivated the research presented in this paper.

Let $\left(\mathcal{T}^{h}\right)_{0<h<1}$ be a family of regular and quasi-uniform triangulations of the domain $\overline{\Omega}$ with the mesh size $h$. For the definition of the discretization space of the state functions let us denote

[TABLE]

where $\mathcal{P}_{1}$ consists all polynomial functions of degree less than or equal to 1.

Proposition 4.1**.**

(i) Let $f$ be in $L^{2}(\Omega)$ and $j$ be in $H^{-1/2}(\partial\Omega)$. Then the variational equation

[TABLE]

admits a unique solution $u^{h}\in\mathcal{V}_{1,\diamond}^{h}$. Furthermore, the estimate

[TABLE]

is satisfied. The map $\mathcal{N}^{h}:L^{2}(\Omega)\rightarrow\mathcal{V}_{1,\diamond}^{h}$ from each $f\in L^{2}(\Omega)$ to the unique solution $u^{h}=:\mathcal{N}^{h}_{f}j$ of (4.1) is then called the discrete Neumann operator.

(ii) Let $f$ be in $L^{2}(\Omega)$ and $g$ be in $H^{1/2}_{\diamond}(\partial\Omega)$. The equation

[TABLE]

has a unique solution $v^{h}\in\mathcal{V}_{1,\diamond}^{h}$ with $\gamma v^{h}=g$. Furthermore, the inequality

[TABLE]

is satisfied. The map $\mathcal{D}^{h}:L^{2}(\Omega)\rightarrow\mathcal{V}_{1,\diamond}^{h}$ from each $f\in L^{2}(\Omega)$ to the unique solution $v^{h}=:\mathcal{D}^{h}_{f}g$ of (4.3) is called the discrete Dirichlet operator.

We now can introduce the strictly convex, discrete cost function

[TABLE]

Theorem 4.2**.**

The problem

[TABLE]

attains a unique minimizer $f$ which satisfies the equation

[TABLE]

Remark 4.3**.**

Since $\mathcal{N}^{h}_{f}j_{\delta}$ and $\mathcal{D}^{h}_{f}g_{\delta}$ are both in $\mathcal{V}^{1}_{h}$, so is $f$, provided that $f^{*}\in\mathcal{V}^{1}_{h}$. Thus, taking this into account, a discretization of the set $L^{2}(\Omega)$ can be avoided.

Proof of Theorem 4.2.

One can see easily that the problem $\left(\mathcal{P}_{\rho,\delta}^{h}\right)$ has a unique solution. It remains to show (4.5). Let $f\in L^{2}(\Omega)$ be the minimizer to $\left(\mathcal{P}_{\rho,\delta}^{h}\right)$. The first-order optimality condition yields that ${\Upsilon_{\rho,\delta}^{h}}^{\prime}(f)\xi={\mathcal{J}_{\delta}^{h}}^{\prime}(f)\xi+2\rho(\xi,f-f^{*})=0$ for all $\xi\in L^{2}(\Omega)$. A short computation shows ${\mathcal{J}_{\delta}^{h}}^{\prime}(f)(\xi)=2\left(\xi,\mathcal{N}^{h}_{f}j_{\delta}-\mathcal{D}^{h}_{f}g_{\delta}\right)$ and so obtain $\left(\xi,\frac{1}{\rho}\left(\mathcal{N}^{h}_{f}j_{\delta}-\mathcal{D}^{h}_{f}g_{\delta}\right)+f-f^{*}\right)=0$ for all $\xi\in L^{2}(\Omega)$, which finishes the proof. ∎

From now on $C$ is a generic positive constant which is independent of the mesh size $h$ of $\mathcal{T}^{h}$, the noise level $\delta$ and the regularization parameter $\rho$. Before presenting the convergence of finite element approximations we here state some auxiliary results.

Lemma 4.4**.**

A projection operator $\Pi^{h}_{\diamond}:L^{1}(\Omega)\rightarrow\mathcal{V}^{h}_{1,\diamond}$ exists such that

[TABLE]

Furthermore, it satisfies the properties

[TABLE]

and

[TABLE]

Proof.

Let $\Pi^{h}:L^{1}(\Omega)\rightarrow\mathcal{V}^{h}_{1}$ be the Clement’s mollification interpolation operator, see [17] and some generalizations [10, 11, 41]. We then define the operator

[TABLE]

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 $\Phi\in H^{1}_{\diamond}(\Omega)$

[TABLE]

We note that $\lim_{h\to 0}\varrho^{h}_{\Phi}=0$ and

[TABLE]

in case $\Phi\in H^{2}(\Omega)$. Furthermore, let $(f,j,g)\in L^{2}(\Omega)\times H^{-1/2}(\partial\Omega)\times H^{1/2}_{\diamond}(\partial\Omega)$ be fixed, we denote by

[TABLE]

Then $\lim_{h\to 0}\alpha^{h}_{f,j}=\lim_{h\to 0}\beta^{h}_{f,g}=0.$ In particular, if $\mathcal{N}_{f}j\in H^{2}(\Omega)$ and $\mathcal{D}_{f}g\in H^{2}(\Omega)$, the error estimates

[TABLE]

are satisfied (cf. [13, 16]).

Lemma 4.5**.**

Let $(f_{1},j_{1},g_{1})$ and $(f_{2},j_{2},g_{2})$ be arbitrary in $L^{2}(\Omega)\times H^{-1/2}(\partial\Omega)\times H^{1/2}_{\diamond}(\partial\Omega)$. Then the estimates

[TABLE]

and

[TABLE]

hold for all $h\geq 0$.

Proof.

According to the definition of the discrete Neumann operator, we have for all $\varphi^{h}\in\mathcal{V}_{1,\diamond}^{h}$ that $\int_{\Omega}Q\nabla\mathcal{N}^{h}_{f_{i}}j_{i}\cdot\nabla\varphi^{h}=\left\langle j_{i},\gamma\varphi^{h}\right\rangle+\left(f_{i},\varphi^{h}\right)\mbox{~{}with~{}}i=1,2$. Thus, $\Phi^{h}_{\mathcal{N}}:=\mathcal{N}^{h}_{f_{1}}j_{1}-\mathcal{N}^{h}_{f_{2}}j_{2}$ is the unique solution to the variational problem $\int_{\Omega}Q\nabla\Phi^{h}_{\mathcal{N}}\cdot\nabla\varphi^{h}=\left\langle j_{1}-j_{2},\gamma\varphi^{h}\right\rangle+\left(f_{1}-f_{2},\varphi^{h}\right)$ for all $\varphi^{h}\in\mathcal{V}_{1,\diamond}^{h}$ and so that (4.12) is satisfied. Likewise, we also obtain (4.13). The proof is completed. ∎

Lemma 4.6**.**

Let $\left(\mathcal{T}^{h_{n}}\right)$ be a sequence of triangulations with $\lim_{n\to\infty}h_{n}=0$. Assume that $\left(j_{\delta_{n}},g_{\delta_{n}}\right)$ is a sequence in $H^{-1/2}(\partial\Omega)\times H^{1/2}_{\diamond}(\partial\Omega)$ convergent to $\left(j_{\delta},g_{\delta}\right)$ in the $H^{-1/2}(\partial\Omega)\times H^{1/2}(\partial\Omega)$-norm and $(f_{n})$ is a sequence in $L^{2}(\Omega)$ weakly convergent in $L^{2}(\Omega)$ to $f$, then there holds the inequality

[TABLE]

Proof.

The proof is based upon the mollification operator introduced in Lemma 4.4 together with standard arguments, therefore omitted here. ∎

We now show the convergence of finite element approximations to the identification problem.

Theorem 4.7**.**

Assume that $\lim_{n\to\infty}h_{n}=0$ and $(\delta_{n})$ and $(\rho_{n})$ any positive sequences such that

[TABLE]

where $\alpha^{h_{n}}_{f^{\dagger},j^{\dagger}}$ and $\beta^{h_{n}}_{f^{\dagger},g^{\dagger}}$ are defined by (4.10). Furthermore, assume that $\left(j_{\delta_{n}},g_{\delta_{n}}\right)$ is a sequence in $H^{-1/2}(\partial\Omega)\times H^{1/2}_{\diamond}(\partial\Omega)$ satisfying

[TABLE]

and $f_{n}:=f^{h_{n}}_{\rho_{n},\delta_{n}}$ is the unique minimizer of $\left(\mathcal{P}^{h_{n}}_{\rho_{n},\delta_{n}}\right)$ for each $n\in N$. Then:

(i) The sequence $(f_{n})$ converges in the $L^{2}(\Omega)$-norm to $f^{\dagger}$.

(ii) The corresponding state sequences $\left(\mathcal{N}^{h_{n}}_{f_{n}}j_{\delta_{n}}\right)$ and $\left(\mathcal{D}^{h_{n}}_{f_{n}}g_{\delta_{n}}\right)$ converge in the $H^{1}(\Omega)$-norm to the unique weak solution \Phi^{\dagger}=\Phi^{\dagger}\big{(}f^{\dagger},j^{\dagger},g^{\dagger}\big{)} 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 \Phi^{\dagger}=\Phi^{\dagger}\big{(}f^{\dagger},j^{\dagger},g^{\dagger}\big{)} of (1.1)–(1.3) belonging to $H^{2}(\Omega)$, the estimate (4.11) shows that $0\leq\alpha^{h_{n}}_{f^{\dagger},j^{\dagger}},~{}\beta^{h_{n}}_{f^{\dagger},g^{\dagger}}\leq Ch_{n}$. Therefore, in view of (4.15), the above convergences (i) and (ii) are obtained if the sequence $(\rho_{n})$ is chosen such that

[TABLE]

By regularity theory for elliptic boundary value problems, the regularity assumption $\Phi^{\dagger}\in H^{2}(\Omega)$ is satisfied if the diffusion matrix $Q\in{C^{0,1}(\Omega)}^{d\times d}$, $j^{\dagger}\in H^{1/2}(\partial\Omega)$, $g^{\dagger}\in H^{3/2}(\partial\Omega)$ and either $\partial\Omega$ is smooth of the class $C^{0,1}$ or the domain $\Omega$ is convex (see, for example, [20, 44]).

Proof of Theorem 4.7.

We have from the optimality of $f_{n}$ that

[TABLE]

Since at $f^{\dagger}$ there holds the equation $\mathcal{N}_{f^{\dagger}}j^{\dagger}=\mathcal{D}_{f^{\dagger}}g^{\dagger}$, we infer from Lemma 4.5 that

[TABLE]

which implies from (4.16) that $\lim_{n\to\infty}\mathcal{J}^{h_{n}}_{\delta_{n}}\left(f_{n}\right)=0$ and, by the assumption (4.15),

[TABLE]

So that the sequence $\left(f_{n}\right)$ is bounded in the $L^{2}(\Omega)$-norm. A subsequence not relabelled and an element $\widehat{f}\in L^{2}(\Omega)$ exist such that $\left(f_{n}\right)$ converges weakly in $L^{2}(\Omega)$ to $\widehat{f}$ and

[TABLE]

For any $f\in L^{2}(\Omega)$ we denote by $\mathcal{J}_{0}(f):=\int_{\Omega}Q\nabla\left(\mathcal{N}_{f}j^{\dagger}-\mathcal{D}_{f}g^{\dagger}\right)\cdot\nabla\left(\mathcal{N}_{f}j^{\dagger}-\mathcal{D}_{f}g^{\dagger}\right).$ By (1.7), we have \left\|\mathcal{N}_{\widehat{f}}j^{\dagger}-\mathcal{D}_{\widehat{f}}g^{\dagger}\right\|^{2}_{H^{1}(\Omega)}\leq\frac{1+C_{\Omega}}{C_{\Omega}\underline{q}}\mathcal{J}_{0}\big{(}\widehat{f}\big{)}\leq\frac{1+C_{\Omega}}{C_{\Omega}\underline{q}}\liminf_{n\to\infty}\mathcal{J}^{h_{n}}_{\delta_{n}}\left(f_{n}\right)=0, here we used Lemma 4.6. Thus, $\mathcal{N}_{\widehat{f}}j^{\dagger}=\mathcal{D}_{\widehat{f}}g^{\dagger}$ which infers $\widehat{f}\in\mathcal{I}\left(j^{\dagger},g^{\dagger}\right).$ Now we show $\widehat{f}=f^{\dagger}$ and the sequence $\left(f_{n}\right)$ converges to $\widehat{f}$ in the $L^{2}(\Omega)$-norm. By the definition of the $f^{*}$-minimum-norm solution and (4.18)–(4.19), we get that

[TABLE]

and so that \big{\|}f^{\dagger}-f^{*}\big{\|}^{2}_{L^{2}(\Omega)}=\big{\|}\widehat{f}-f^{*}\big{\|}^{2}_{L^{2}(\Omega)}=\lim_{n\to\infty}\|f_{n}-f^{*}\|^{2}_{L^{2}(\Omega)}. By the uniqueness of the minimum-norm solution and the sequence $\left(f_{n}\right)$ weakly converging in $L^{2}(\Omega)$ to $\widehat{f}$, we conclude that $\widehat{f}=f^{\dagger}$ and the sequence $\left(f_{n}\right)$ in fact converges in the $L^{2}(\Omega)$-norm to $\widehat{f}$. Finally, we show the sequences $(\mathcal{N}^{h_{n}}_{f_{n}}j_{\delta_{n}})$ and $(\mathcal{D}^{h_{n}}_{f_{n}}g_{\delta_{n}})$ converge to $\Phi^{\dagger}=\mathcal{N}_{f^{\dagger}}j^{\dagger}=\mathcal{D}_{f^{\dagger}}g^{\dagger}$ in the $H^{1}(\Omega)$-norm. Indeed, by Lemma 4.5, we obtain that

[TABLE]

Similarly, we also get $\left\|\mathcal{D}^{h_{n}}_{f_{n}}g_{\delta_{n}}-\mathcal{D}_{f^{\dagger}}g^{\dagger}\right\|_{H^{1}(\Omega)}\leq C\left(\left\|g_{\delta_{n}}-g^{\dagger}\right\|_{H^{1/2}(\partial\Omega)}+\left\|f_{n}-f^{\dagger}\right\|_{L^{2}(\Omega)}+\beta^{h_{n}}_{f^{\dagger},g^{\dagger}}\right)\rightarrow 0$ as $n$ tends to $\infty$, which finishes the proof. ∎

5 Convergence rates for the discretized problem

We are now in a position to state the main theorem on convergence rates for the general case of finite element discretized regularized solutions with noise level ($\delta>0$ and $h>0$). The source condition (3.5) will play a prominent role in this context.

Theorem 5.1**.**

Assume that the condition (3.5) is fulfilled. Then, we have the error estimate and convergence rate

[TABLE]

where $f^{h}:=f_{\rho,\delta}^{h}$ is the unique minimizer of $\left(\mathcal{P}_{\rho,\delta}^{h}\right)$ and $\mathcal{D}_{0}\gamma\mathcal{N}_{w}j^{\dagger}-g^{\dagger}$ is the unique weak solution to the Dirichlet problem

[TABLE]

and $\alpha^{h}_{f^{\dagger},j^{\dagger}}$, $\beta^{h}_{f^{\dagger},g^{\dagger}}$, $\varrho^{h}_{\mathcal{N}_{w}j^{\dagger}}$, $\varrho^{h}_{\mathcal{N}_{f^{\dagger}}j^{\dagger}}$ and $\varrho^{h}_{\mathcal{D}_{0}\gamma\mathcal{N}_{w}j^{\dagger}-g^{\dagger}}$ come from (4.8) and (4.10).

Remark 5.2**.**

In case (cf. Remark 4.8) $\mathcal{N}_{f^{\dagger}}j^{\dagger},\mathcal{N}_{w}j^{\dagger},\mathcal{D}_{0}\gamma\mathcal{N}_{w}j^{\dagger}-g^{\dagger}\in H^{2}(\Omega),$ by (4.9) and (4.11), we have

[TABLE]

and so that the following convergence rate is obtained

[TABLE]

Remark 5.3**.**

Let \Phi^{\dagger}=\Phi^{\dagger}\big{(}f^{\dagger},j^{\dagger},g^{\dagger}\big{)} be the weak solution of (1.1)–(1.3). Then the convergence rate

[TABLE]

is also established. Indeed, the desired equation directly follows from (5.1) and the following inequalities

[TABLE]

and \left\|\mathcal{D}^{h}_{f^{h}}g_{\delta}-\mathcal{D}_{f^{\dagger}}g^{\dagger}\right\|_{H^{1}(\Omega)}\leq C\left(\delta+\big{\|}f^{h}-f^{\dagger}\big{\|}_{L^{2}(\Omega)}+\beta^{h}_{f^{\dagger},g^{\dagger}}\right), here we used Lemma 4.5.

Proof of Theorem 5.1.

In view of (4.17) we first have that \mathcal{J}^{h}_{\delta}\big{(}f^{\dagger}\big{)}\leq C\left(\delta^{2}+\left(\alpha^{h}_{f^{\dagger},j^{\dagger}}\right)^{2}+\left(\beta^{h}_{f^{\dagger},g^{\dagger}}\right)^{2}\right). The optimality of $f^{h}$ yields \mathcal{J}^{h}_{\delta}\big{(}f^{h}\big{)}+\rho\big{\|}f^{h}-f^{*}\big{\|}^{2}_{L^{2}(\Omega)}\leq\mathcal{J}^{h}_{\delta}\big{(}f^{\dagger}\big{)}+\rho\big{\|}f^{\dagger}-f^{*}\big{\|}^{2}_{L^{2}(\Omega)}. This gives

[TABLE]

Since $\mathcal{N}_{f^{\dagger}}j^{\dagger}=\mathcal{D}_{f^{\dagger}}g^{\dagger}$, it follows that

[TABLE]

From (1.10), we infer

[TABLE]

and

[TABLE]

This in turn implies

[TABLE]

Since $\gamma\left(\mathcal{D}_{f^{\dagger}}g^{\dagger}-\mathcal{D}_{w}g^{\dagger}\right)=0$, it follows from (1.11) that

[TABLE]

holds. We thus infer from (5.3)–(5.5) the identity

[TABLE]

We note again that $\mathcal{N}_{f^{\dagger}}j^{\dagger}=\mathcal{D}_{f^{\dagger}}g^{\dagger}$ and $\gamma\left(\mathcal{D}_{f^{\dagger}}g^{\dagger}-\mathcal{D}_{f^{h}}g^{\dagger}\right)=0$. Then, together with (1.10) and (1.11), the last two terms on the right hand side of (5) satisfy

[TABLE]

Thus, we obtain from (5)

[TABLE]

Next, we abbreviate $W=\mathcal{N}_{w}j^{\dagger}-\mathcal{N}_{f^{\dagger}}j^{\dagger}$ and note

[TABLE]

Then we get

[TABLE]

To prepare the estimation of those three addends we start with writing

[TABLE]

Since $\mathcal{D}_{f^{h}}g^{\dagger}-\mathcal{D}^{h}_{f^{h}}g^{\dagger}\in H^{1}_{0}(\Omega)$, we then get

[TABLE]

Since $\gamma\left(W-\mathcal{D}_{0}\gamma W\right)=\gamma W-\gamma\mathcal{D}_{0}\gamma W=\gamma W-\gamma W=0$, we infer $\Pi^{h}_{\diamond}\left(W-\mathcal{D}_{0}\gamma W\right)\in\mathcal{V}_{1,0}^{h}=\mathcal{V}_{1}^{h}\cap H^{1}_{0}(\Omega)$ and then obtain from (1.11) and (4.3) that $\int_{\Omega}Q\nabla\left(\mathcal{D}_{f^{h}}g^{\dagger}-\mathcal{D}^{h}_{f^{h}}g^{\dagger}\right)\cdot\nabla\Pi^{h}_{\diamond}\left(W-\mathcal{D}_{0}\gamma W\right)=0$ holds. Hence we have

[TABLE]

where we use (5.7). Similarly, since $\Pi^{h}_{\diamond}W\in\mathcal{V}^{h}_{\diamond}$ and by (1.10) and (4.1), we get

[TABLE]

Now we are in the position to estimate $I_{1}-I_{3}$. Combining (5) with (5.10), we obtain

[TABLE]

Now, using Lemma 4.5, we arrive at

[TABLE]

Since for a.e. in $\Omega$ the matrix $Q(x)$ is positive definite, the root $Q(x)^{1/2}$ is then well defined. Thus, using the Cauchy-Schwarz inequality and Young’s inequality, we estimate $I_{3}$ as

[TABLE]

It follows from (5) and (5.11)–(5.13) that

[TABLE]

holds, which together with (5) implies

[TABLE]

Since \left\|\mathcal{D}^{h}_{f^{h}}g_{\delta}-\mathcal{N}^{h}_{f^{h}}j_{\delta}\right\|^{2}_{H^{1}(\Omega)}\leq C\mathcal{J}^{h}_{\delta}\big{(}f^{h}\big{)}, (5.1) now directly follows from (5), which finishes the proof. ∎

6 Conjugate gradient method and numerical test

In this section we will utilize the conjugate gradient (CG) method (see, for example, [21, 26]) to find the minimizes of the strictly convex, discrete regularized problem $\left(\mathcal{P}_{\rho,\delta}^{h}\right)$. Let $\nabla\Upsilon_{\rho,\delta}^{h}(f)=2\left(\mathcal{N}^{h}_{f}j_{\delta}-\mathcal{D}^{h}_{f}g_{\delta}\right)+2\rho(f-f^{*})$ be the $L^{2}$-gradient of the cost function $\Upsilon_{\rho,\delta}^{h}$ at $f$ (see Proof of Theorem 4.2), where $f^{*}\in\mathcal{V}^{h}_{1}$. Then the sequence of iterates via this algorithm is generated by $f^{0}\in L^{2}(\Omega)\cap\mathcal{V}^{h}_{1}$ and $f^{k+1}:=f^{k}+t^{k}d^{k}$ for $k\geq 0$, where

[TABLE]

A short computation shows that

[TABLE]

Consequently, the CG method then reads as follows: giving an initial approximation $f^{0}\in\mathcal{V}^{h}_{1}$, number of iterations $N$ and a positive constants $\tau_{1},\tau_{2}$. Computing

[TABLE]

and setting

[TABLE]

Below we illustrate the theoretical result with numerical examples. For this purpose we consider the the boundary value problem

[TABLE]

We assume that entries of the known symmetric diffusion matrix $Q$ are discontinuous which are defined as $q_{11}=3\chi_{\Omega_{11}}+\chi_{\Omega\setminus\Omega_{11}},~{}q_{12}=\chi_{\Omega_{12}},~{}q_{22}=4\chi_{\Omega_{22}}+2\chi_{\Omega\setminus\Omega_{22}}$, where $\chi_{D}$ is the characteristic function of the Lebesgue measurable set $D$ and

[TABLE]

The identified source function $f^{\dagger}\in L^{2}(\Omega)$ in (6.1) is assumed to be discontinuous and defined as

[TABLE]

where

[TABLE]

For the discretization we divide the interval $(-1,1)$ into $\ell$ equal segments and so that the domain $\Omega=(-1,1)^{2}$ is divided into $2\ell^{2}$ triangles, where the diameter of each triangle is $h_{\ell}=\frac{\sqrt{8}}{\ell}$. In the minimization problem $\left(\mathcal{P}_{\rho,\delta}^{h}\right)$ we take $h=h_{\ell}$ and $\rho=\rho_{\ell}=0.01h_{\ell}$. We use Algorithm 1 which is described above for computing the numerical solution of the problem $\left(\mathcal{P}_{\rho_{\ell},\delta_{\ell}}^{h_{\ell}}\right)$. As an a-priori estimate and the initial approximation we choose $f^{*}:=0$ and $f^{0}(x):=\chi_{(0,1]\times[-1,1]}-\chi_{[-1,0]\times[-1,1]}$.

Example 6.1**.**

In this first example $j^{\dagger}\in H^{-1/2}(\partial\Omega)$ is chosen to be the piecewise constant function defined by

[TABLE]

Then $g^{\dagger}\in H^{1/2}_{\diamond}(\partial\Omega)$ is defined as $g^{\dagger}=\gamma\mathcal{N}_{f^{\dagger}}j^{\dagger}$. We mention that, to avoid a so-called inverse crime, we generate the given data on a finer grid than those used in the computations. For this purpose we first solve the problem (6.1) supplemented with the Neumann boundary condition in (6.2) on the very fine grid with $\ell=128$, and then use this numerical approximation as substitute for $(j^{\dagger},g^{\dagger})$ in our computational considerations below.

For observations with noise we assume that

[TABLE]

where $R_{j^{\dagger}}$ and $R_{g^{\dagger}}$ are $\partial M^{h_{\ell}}\times 1$-matrices of random numbers on the interval $(-1,1)$ which are generated by the MATLAB function “rand”, and $\partial M^{h_{\ell}}$ is the number of boundary nodes of the triangulation $\mathcal{T}^{h_{\ell}}$. The measurement error is then computed as \delta_{\ell}=\big{\|}j_{\delta_{\ell}}-j^{\dagger}\big{\|}_{L^{2}(\partial\Omega)}+\big{\|}g_{\delta_{\ell}}-g^{\dagger}\big{\|}_{L^{2}(\partial\Omega)}. To satisfy the condition $\delta_{\ell}\cdot\rho^{-1/2}_{\ell}\to 0$ as $\ell\to\infty$ in Theorem 4.7 we below take $\theta_{\ell}=h_{\ell}\sqrt{\rho_{\ell}}$.

We start with the coarsest level $\ell=4$. At each iteration $k$ we compute

[TABLE]

where $\tau_{1}:=10^{-6}h_{\ell}^{1/2}$ and $\tau_{2}:=10^{-4}h_{\ell}^{1/2}$. Then the iteration was stopped if $\mbox{Tolerance}\leq 0$ or the number of iterations reached the maximum iteration count of 600. After obtaining the numerical solution of the first iteration process with respect to the coarsest level $\ell=4$, we use its interpolation on the next finer mesh $\ell=8$ as an initial approximation $f^{0}$ for the algorithm on this finer mesh, and proceed similarly on the preceding refinement levels.

Let $f_{\ell}$ be the function which is obtained at the final iterate of Algorithm 1 corresponding to the refinement level $\ell$. Furthermore, let $\mathcal{N}^{h_{\ell}}_{f_{\ell}}j_{\delta_{\ell}}$ and $\mathcal{D}^{h_{\ell}}_{f_{\ell}}g_{\delta_{\ell}}$ denote the computed numerical solution to the Neumann and Dirichlet problem

[TABLE]

respectively. The notations $\mathcal{N}^{h_{\ell}}_{f^{\dagger}}j^{\dagger}$ and $\mathcal{D}^{h_{\ell}}_{f^{\dagger}}g^{\dagger}$ of the exact numerical solutions are to be understood similarly. We use the following abbreviations for the errors

[TABLE]

The numerical results are summarized in Table 1 and Table 2, where we present the refinement level $\ell$, mesh size $h_{\ell}$ of the triangulation, regularization parameter $\rho_{\ell}$, measured noise $\delta_{\ell}$, number of iterations, value of tolerances and errors $L^{2}_{f}$, $L^{2}_{\mathcal{N}}$, $L^{2}_{\mathcal{D}}$, $H^{1}_{\mathcal{N}}$ and $H^{1}_{\mathcal{D}}$. Their experimental order of convergence (EOC) is presented in Table 3, where $\mbox{EOC}_{\Theta}:=\dfrac{\ln\Theta(h_{1})-\ln\Theta(h_{2})}{\ln h_{1}-\ln h_{2}}$ and $\Theta(h)$ is an error function of the mesh size $h$.

All figures presented correspond to $\ell=64$. Figure 1 from left to right shows the computed numerical solution $f_{\ell}$ of the algorithm at the final $579^{\mbox{th}}$-iteration, and the differences $\mathcal{N}^{h_{\ell}}_{f^{\dagger}}j^{\dagger}-\mathcal{N}^{h_{\ell}}_{f_{\ell}}j_{\delta_{\ell}}$,  $\mathcal{D}^{h_{\ell}}_{f^{\dagger}}g^{\dagger}-\mathcal{D}^{h_{\ell}}_{f_{\ell}}g_{\delta_{\ell}}$ and $\mathcal{D}^{h_{\ell}}_{f_{\ell}}g_{\delta_{\ell}}-\mathcal{N}^{h_{\ell}}_{f_{\ell}}j_{\delta_{\ell}}$.

Example 6.2**.**

In present example we assume that multiple measurements are available, say $\left(j_{\delta}^{i},g_{\delta}^{i}\right)_{i=1,\ldots,I}$. Then, problem $\left(\mathcal{P}^{h}_{\rho,\delta}\right)$ in Section 4 is given by

[TABLE]

which also attains a solution $\bar{f}^{h}_{\rho,\delta}$. The Neumann boundary condition in the equation (6.2) is chosen in the same form as (6.3), i.e.

[TABLE]

and depends on the constants $A,B,C$ and $D$. Let $g^{\dagger}_{(A,B,C,D)}:=\gamma\mathcal{N}_{f^{\dagger}}j^{\dagger}_{(A,B,C,D)}$ and assume that noisy observations are given by

[TABLE]

where $R_{j^{\dagger}_{(A,B,C,D)}}$ and $R_{g^{\dagger}_{(A,B,C,D)}}$ denote $\partial M^{h_{\ell}}\times 1$-matrices of random numbers on the interval $(-1,1)$. Different from (6.4), the constant $\theta$ appeared in the equation (6.6) is now independent of the grid level $\ell$.

In the case $(A,B,C,D)=(1,2,3,4)$ we have a single noisy measurement couple, i.e. $I=1$. We now fix $D=4$, and let $(A,B,C)$ take all permutations $\mathcal{S}_{3}$ of the set $\{1,2,3\}$. Then, the equations (6.5)–(6.6) generate $I=6$ measurements. Similarly, if $(A,B,C,D)$ takes all permutations $\mathcal{S}_{4}$ of $\{1,2,3,4\}$ we get $I=16$ measurements. With $\theta=0.1$ and $\ell=64$ we compute the noise level

[TABLE]

The corresponding numerical results for the multiple measurement case are presented in the Table 4.

Finally, in Figure 2 from left to right we show the interpolation $I_{1}^{h_{\ell}}f^{\dagger}$ of the exact source and the computed numerical solution $q_{\ell}$ of the algorithm at the final iteration for $\ell=64$, $\theta=0.1$, and $I=16,6,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 the University of Hamburg.

B. Hofmann gratefully acknowledges support by the German Research Foundation (DFG) under grant HO 1454/12-1.

T. N. T. Quyen gratefully acknowledges support of the Alexander von Humboldt Foundation.