Finite element approximation of source term identification with TV-regularization

TL;DR

This paper develops a finite element approach combined with TV-regularization to accurately recover discontinuous source terms in elliptic systems from boundary measurements, ensuring stability and convergence.

Contribution

It introduces a finite element discretization for TV-regularized source identification and proves the convergence and stability of the proposed algorithm.

Findings

The method is stable and convergent.

Numerical experiments validate theoretical results.

Abstract

In this paper we investigate the problem of recovering the source term in an elliptic system from a measurement of the state on a part of the boundary. For the particular interest in reconstructing probably discontinuous sources, we use the standard least squares method with the total variation regularization. The finite element method is then applied to discretize the minimization problem, we show the stability and the convergence of this technique. Furthermore, we have proposed an algorithm to stably solve the minimization problem. We prove the iterate sequence generated by the derived algorithm converging to a minimizer of the regularization problem, and that convergence measurement is also established. Finally, a numerical experiment is presented to illustrate our theoretical 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 tolerances.

$\ell$	$h_{\ell }$	${\rho }_{\ell }$	${\delta }_{\ell }$	Iterate	Tolerance
4	0.7071	8.4090e-4	2.3763e-2	152	-4.0518e-6
8	0.3536	5.9460e-4	8.8717e-3	264	-2.8826e-7
16	0.1766	4.2045e-4	2.5872e-3	343	-1.4460e-8
32	8.8388e-2	2.9730e-4	1.1817e-3	471	-1.1382e-7
64	4.4194e-2	2.1022e-4	5.6112e-4	561	-4.6128e-8

Table 2. Table 2: Errors ‖ f † − f ℓ ‖ Ω subscript norm superscript 𝑓 † subscript 𝑓 ℓ Ω \big{\|}f^{\dagger}-f_{\ell}\big{\|}_{\Omega} , ‖ u † − u ℓ ‖ Ω subscript norm superscript 𝑢 † subscript 𝑢 ℓ Ω \big{\|}u^{\dagger}-u_{\ell}\big{\|}_{\Omega} and ‖ u † − u ℓ ‖ 1 , Ω subscript norm superscript 𝑢 † subscript 𝑢 ℓ 1 Ω \big{\|}u^{\dagger}-u_{\ell}\big{\|}_{1,\Omega} .

$\ell$	${\Vert f^{†}-f_{\ell }\Vert }_{\Omega }$	${\Vert u^{†}-u_{\ell }\Vert }_{\Omega }$	${\Vert u^{†}-u_{\ell }\Vert }_{1,\Omega }$
4	0.6829	1.5022e-2	3.9121e-2
8	0.3506	8.5904e-3	3.3448e-2
16	0.1692	3.2470e-3	2.2919e-2
32	0.1095	1.2926e-3	1.5372e-2
64	0.7403e-2	6.0377e-4	1.2252e-2

Equations

$$\displaystyle-\nabla\cdot\big{(}\alpha\nabla u\big{)}+\beta u$$

$$\alpha \nabla u\cdot \vec{n}+\sigma u$$

$$\underset{f\in F_{ad}}{\min }J(f),J(f):=\frac{1}{2}\Vert u(f)-z{\Vert }_{\Gamma }^2+\rho TV(f)$$

$$\underset{f\in F_{ad}^h}{\min }{J}^h(f),J^h(f):=\frac{1}{2}\Vert u^h(f)-z{\Vert }_{\Gamma }^2+\rho TV(f),$$

$$\Vert f_{n+1}^h-f_n^h{\Vert }^2=\mathcal{O}\left(\frac{1}{n}\right)$$

$$\displaystyle-\nabla\cdot\big{(}\alpha\nabla u\big{)}+\beta u$$

$$\alpha \nabla u\cdot \vec{n}+\sigma u$$

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

$$a(u,v):=(\alpha \nabla u,\nabla v{)}_{\Omega }+(\beta u,v{)}_{\Omega }+(\sigma u,v{)}_{\partial \Omega }\text{mbox}and{l}_j^f(v):=(f,v{)}_{\Omega }+[j,v{]}_{\partial \Omega },$$

$$c_1\Vert u{\Vert }_{1,\Omega }^2\le a(u,u)\le c_2\Vert u{\Vert }_{1,\Omega }^2$$

$$a(u,v)=l_j^f(v)$$

$$\Vert u{\Vert }_{1,\Omega }\le \frac{1}{c_1}\max \left(1,\Vert \gamma {\Vert }_{\mathcal{L}\left(H^1(\Omega ),H^{1/2}(\partial \Omega )\right)}\right)\left(\Vert j{\Vert }_{H^{-1/2}(\partial \Omega )}+\Vert f{\Vert }_{\Omega }\right)$$

$$\Vert \gamma u{\Vert }_{\partial \Omega }\le c_{\gamma }\Vert u{\Vert }_{1,\Omega }$$

$$C_{\Omega }(v,v{)}_{\Omega }\le (\nabla v,\nabla v{)}_{\Omega }\text{mbox}forallv\in H_{\diamond }^1(\Omega ):=\left\{u\in H^1(\Omega )|(u,1{)}_{\Omega }=0\right\}$$

$$\underset{f\in F_{ad}}{\min }J(f),J(f):=\frac{1}{2}\Vert u(f)-z{\Vert }_{\Gamma }^2+\rho TV(f)$$

$$F_{ad}:=\left\{f\in BV(\Omega )|-\infty <\underline{f}\le f(x)\le \overline{f}<\infty \text{mbox}fora.e.in\Omega \right\}$$

$$\underset{f\in F_{ad}^h}{\min }{J}^h(f),J^h(f):=\frac{1}{2}\Vert u^h(f)-z{\Vert }_{\Gamma }^2+\rho TV(f),$$

$$\mathbb{T}:BV(\Omega )\to L^2(\Omega )$$

$$\underset{e\in BV(\Omega )}{\inf }E(e),E(e):=\frac{1}{2}\Vert \mathbb{T}e-y{\Vert }_{\Omega }^2+\rho TV(e),$$

$$\underset{e\in BV(\Omega )}{\inf }E(e)=\underset{e\in BV(\Omega )}{\inf }\underset{p\in L^1(\Omega ;{\mathbb{R}}^d)}{\sup }\mathcal{E}(e,p),$$

$$\mathcal{E}(e,p):=\frac{1}{2}\Vert \mathbb{T}e-y{\Vert }_{\Omega }^2+\rho (\nabla e,p{)}_{\Omega }-I_{{\mathcal{B}}_1(L^1(\Omega ;{\mathbb{R}}^d))}(p),$$

$$\underset{e^h\in {\mathcal{V}}_1^h}{\inf }\underset{p^h\in ({\mathcal{V}}_0^h{)}^d}{\sup }\left(\frac{1}{2}\Vert \mathbb{T}{e}^h-y{\Vert }_{\Omega }^2+\rho (\nabla e^h,p^h{)}_{\Omega }-I_{{\mathcal{B}}_1(L^1(\Omega ;{\mathbb{R}}^d))}(p^h)\right).$$

$$e_{n+1}^h$$

$${\tilde{e}}_{n+1}^h=2e_{n+1}^h-e_n^h$$

$$p_{n+1}^h$$

$$p_{n+1}^h=\frac{p_n^h+\frac{\tau \rho }{\theta }\nabla {\tilde{e}}_{n+1}^h}{\max \{1,|p_n^h+\frac{\tau \rho }{\theta }\nabla {\tilde{e}}_{n+1}^h|\}}.$$

$$f_{n+1}^h$$

$${\tilde{f}}_{n+1}^h=2f_{n+1}^h-f_n^h$$

$$p_{n+1}^h$$

$$f_{n+1}^h={\mathcal{P}}_{F_{ad}^h}^{L^2}\left(f_n^h-\tau (u_a^h(f_n^h)-\rho \text{mbox}div{p}_n^h)\right),$$

$$(-\text{mbox}div{q}^h,g^h)=(q^h,\nabla g^h{)}_{\Omega }:=\sum _{T\in {\mathcal{T}}^h}\sum _{j=1}^d|T|{q_j^h}_{{|}_T}\frac{\partial g^h}{\partial x_j}{|}_T$$

$$\Vert {\mu }_{n+1}^h-{\mu }_n^h{\Vert }_{{\mathcal{V}}_1^h\times ({\mathcal{V}}_0^h{)}^d}^2=\mathcal{O}\left(\frac{1}{n}\right)\text{ for }n\to \infty ,$$

$$\displaystyle TV(f):=\int_{\Omega}\left|\nabla f\right|:=\sup\left\{\int_{\Omega}f\text{div}~{}G~{}\big{|}~{}G\in C^{1}_{c}(\Omega)^{d},~{}\left|G(x)\right|_{\infty}\leq 1,~{}x\in\Omega\right\}<\infty,$$

Full text

Finite element approximation of source term identification with TV-regularization

Michael Hinzea††Email: michael.hinze@uni-hamburg.de  $\cdot$  quyen.tran@uni-goettingen.de  $\cdot$  Tran Nhan Tam Quyenb,∗††∗Corresponding author

aDepartment of Mathematics, University of Hamburg, Bundesstr. 55, 20146 Hamburg, Germany

bInstitute for Numerical and Applied Mathematics, University of Goettingen, Lotzestr. 16-18, 37083 Goettingen, Germany

Abstract: In this paper we investigate the problem of recovering the source $f$ in the elliptic system

[TABLE]

from an observation $z$ of the state $u$ on a part $\Gamma$ of the boundary $\partial\Omega$, where the functionals $\alpha,\beta,\sigma$ and $j$ are given. For the particular interest in reconstructing probably discontinuous sources, we use the standard least squares method with the total variation regularization, i.e. we consider a minimizer of the minimization problem

[TABLE]

as reconstruction. Here $u(f)$ denotes the unique weak solution of the above elliptic system which depends on the source term $f$, $TV(f)$ is the total variation of $f$, $\rho>0$ is the regularization parameter, the admissible set $F_{ad}:=\left\{f\in BV(\Omega)~{}|~{}-\infty<\underline{f}\leq f(x)\leq\overline{f}<\infty\quad\mbox{for a.e. in}\quad\Omega\right\}$ with $\underline{f}$, $\overline{f}$ being given constants, and $BV(\Omega)$ is the Banach space of all bounded total variation functionals. We approximate problem $\left(\mathcal{P}\right)$ with piecewise linear and continuous finite elements, where $u^{h}\in\mathcal{V}^{h}_{1}$ denotes the corresponding finite element approximation of $u$. This leads to the minimization problem

[TABLE]

where $F^{h}_{ad}=F_{ad}\cap\mathcal{V}^{h}_{1}$. We in Theorems 3.1, 3.5 provide the numerical analysis for the discrete solutions $f^{h}$ of \big{(}\mathcal{P}^{h}\big{)}, and also propose an algorithm to stably solve this discrete minimization problem, where we are lead by the algorithmic developments of [8, 49]. In particular we prove that the iteration sequence \big{(}f^{h}_{n}\big{)}_{n} generated by this algorithm converges to a minimizer of \big{(}\mathcal{P}^{h}\big{)}, and that convergence measure of the kind

[TABLE]

is satisfied. Finally, a numerical experiment is presented to illustrate our theoretical findings.

Key words and phrases: Inverse source problem, boundary observation, total variation regularization, ill-posedness, finite element method, stability and convergence, elliptic boundary value problem.

AMS Subject Classifications: 35R25; 47A52; 35R30; 65J20; 65J22.

1 Introduction

Let $\Omega$ be an open bounded connected set of $\mathbb{R}^{d},~{}d\geq 2$ with the polygonal boundary $\partial\Omega$. In this paper we investigate the inverse problem of identifying the source term $f\in L^{2}(\Omega)$ in the elliptic system

[TABLE]

from a boundary measurement $z\in L^{2}(\Gamma)$ of the exact data $g^{\dagger}:=u_{|\Gamma}$, where $\vec{n}$ is the unit outward normal on $\partial\Omega$ and $\Gamma\subset\partial\Omega$ is a relatively open subset of the boundary.

In the system (1.1) the Robin boundary condition j\in H^{-1/2}(\partial\Omega):=\big{(}H^{1/2}(\partial\Omega)\big{)}^{*} and special functionals $\alpha,\beta,\sigma$ are assumed to be given, where $\sigma\in L^{\infty}(\partial\Omega)$ with $\sigma(x)\geq 0$ a.e. on $\partial\Omega$, $\beta\in L^{\infty}(\Omega)$ with $\beta(x)\geq 0$ a.e. in $\Omega$ and $\alpha:=\left(\alpha_{rs}\right)_{1\leq r,s\leq d}\in{L^{\infty}(\Omega)}^{d\times d}$ is a symmetric diffusion matrix satisfying the uniformly elliptic condition

[TABLE]

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

We start with some notations. The expressions $(\cdot,\cdot)_{\Omega}$ and $(\cdot,\cdot)_{\partial\Omega}$ stand respectively for the scalar products of the Lebesgue spaces $L^{2}(\Omega)$ and $L^{2}(\partial\Omega)$, while $[\cdot,\cdot]_{\partial\Omega}$ denotes the dual pair (\cdot,\cdot)_{\big{(}H^{-1/2}(\partial\Omega),H^{1/2}(\partial\Omega)\big{)}}. Let

[TABLE]

where $u,v\in H^{1}(\Omega)$.

We assume that $\beta(x)\geq\beta_{0}>0$ a.e. in $\Omega$ or $\sigma(x)\geq\sigma_{0}>0$ a.e. on $\partial\Omega$. Then, there exist positive constants $c_{1},c_{2}$ such that

[TABLE]

for all $u\in H^{1}(\Omega)$, where $\|u\|_{1,\Omega}:=\left((\nabla u,\nabla u)_{\Omega}+(u,u)_{\Omega}\right)^{1/2}$ denotes the usual norm of the Sobolev space $H^{1}(\Omega)$ (cf. [41]). Thus, the expression $a(u,v)$ generates a scalar product on the space $H^{1}(\Omega)$ equivalent to the usual one. Consequently, for each $f\in L^{2}(\Omega)$ the Robin boundary value problem (1.1) defines a unique weak solution $u$ in the sense that $u:=u_{j}(f):=u(f)\in H^{1}(\Omega)$ such that the variational equation

[TABLE]

is satisfied for all $v\in H^{1}(\Omega)$. Furthermore, the estimate

[TABLE]

holds true, where $\|f\|_{\Omega}:=\|f\|_{L^{2}(\Omega)}$ and $\gamma:H^{1}(\Omega)\to H^{1/2}(\partial\Omega)$ denotes the continuous Dirichlet trace operator. Recall that $\gamma:H^{1}(\Omega)\to L^{2}(\partial\Omega)$ is compact (cf., e.g., [37, 40]) and there is a positive constant $c_{\gamma}$ such that

[TABLE]

for all $u\in H^{1}(\Omega)$.

In case $\beta=0$ and $\sigma=0$, by the Poincaré-Friedrichs inequality (cf., e.g., [44])

[TABLE]

for some constant $C_{\Omega}>0$ depending only on $\Omega$, the expression $(\alpha\nabla u,\nabla v)_{\Omega}$ is a scalar product on the space $H^{1}_{\diamond}(\Omega)$ which is equivalent to the usual one. Therefore, in view of Riesz’s representation theorem, for each $f\in L^{2}(\Omega)$ the problem (1.1) also has a unique weak solution $u:=u_{j}(f):=u(f)\in H^{1}_{\diamond}(\Omega)$ which is defined via the variational equation $(\alpha\nabla u,\nabla v)_{\Omega}=l^{f}_{j}(v)$ for all $v\in H^{1}_{\diamond}(\Omega)$ and further satisfied the estimate $\|u\|_{1,\Omega}\leq C\left(\|j\|_{H^{-1/2}(\partial\Omega)}+\|f\|_{\Omega}\right)$. In the sequel we thus consider the case $\inf_{\Omega}\beta(x)+\inf_{\partial\Omega}\sigma(x)>0$ only. All results in present paper are still valid for the case $\beta=\sigma=0$.

The inverse problem is stated as follows:

Given a boundary measurement $z\in L^{2}(\Gamma)$ of the exact $g^{\dagger}:=u_{|\Gamma}$ of (1.1), find $f\in L^{2}(\Omega)$.

For this purpose and with the particular interest in reconstructing probably discontinuous sources, we use the standard least squares method with the total variation regularization, i.e. we consider a minimizer of the minimization problem

[TABLE]

as reconstruction. Here $\rho\in(0,1)$ is the regularization parameter and

[TABLE]

with $\underline{f}$, $\overline{f}$ being given constants. We note that the lower and upper bound imposing upon the admissible set $F_{ad}$ guarantees the existence of a minimizer to $\left(\mathcal{P}\right)$ (cf. [1, 12, 17, 51]).

Let $u^{h}$ be the approximation of $u$ in the finite dimensional space $\mathcal{V}^{h}_{1}:=\left\{v^{h}\in C(\overline{\Omega})~{}|~{}{v^{h}}_{|T}\in\mathcal{P}_{1}(T),~{}~{}\forall T\in\mathcal{T}^{h}\right\}$ of piecewise linear, continuous finite elements. We then consider the discrete total variation regularized problem corresponding to $\left(\mathcal{P}\right)$ which is given by the following minimization problem

[TABLE]

where $F^{h}_{ad}:=F_{ad}\cap\mathcal{V}_{1}^{h}$. Let us point to the reader here that a discrete total variation approach with finite elements has been recently proposed for imaging in [30]. Likewise $\left(\mathcal{P}\right)$, we can show that the problem \big{(}\mathcal{P}^{h}\big{)} admits a minimizer $f^{h}$ for each $h>0$. However, due to the lack of strict convexity of the cost functional, a solution of \big{(}\mathcal{P}^{h}\big{)} may be non-unique.

As $\rho$ is fixed, we show that every solution sequence $(f^{h})_{h}$ of the problems \big{(}\mathcal{P}^{h}\big{)}_{h} has a subsequence which converges to a solution of the problem $\left(\mathcal{P}\right)$, in $L^{s}(\Omega)$-norm as well as in the total variation, where $\forall s\in[1,\infty)$ (see §3, Theorem 3.1). Furthermore, in case $h\to 0$, $\|g^{\dagger}-z\|_{\Gamma}\leq\delta_{z}\to 0$ and $\rho=\rho(h,\delta_{z})$ is chosen in a suitable way, every sequence of solutions to \big{(}\mathcal{P}^{h}\big{)}_{h} can be extracted a subsequence that converges to a sought source with the total variation-minimizing property (see §3, Theorem 3.5). To the best of the author’s knowledge those results are new for problem \big{(}\mathcal{P}\big{)}.

For the numerical solution of problem \big{(}\mathcal{P}^{h}\big{)}, we adopt the linearized primal-dual algorithm proposed in [49]. For the convenience of the reader we in the following sketch this approach. Starting with the equation

[TABLE]

with $\mathbb{T}$ being a bounded linear operator, the authors have considered the unconstrained minimization problem

[TABLE]

where $y$ is the given data. Due to the $TV$-dual representation approach developed by [8, 9, 10, 14], (1.8) is rewritten as a saddle-point problem

[TABLE]

where

[TABLE]

$\mathcal{B}_{1}(L^{1}(\Omega;\mathbb{R}^{d})):=\left\{p=(p_{j})_{j=1}^{d}\in L^{1}(\Omega;\mathbb{R}^{d})~{}|~{}\|p\|_{\infty}:=\max_{j=1,...,d}\|p_{j}\|_{L^{\infty}(\Omega)}\leq 1\right\}$, and $I_{\mathcal{B}_{1}(L^{1}(\Omega;\mathbb{R}^{d}))}$ denotes the indicator function of the set $\mathcal{B}_{1}(L^{1}(\Omega;\mathbb{R}^{d}))$. Using the piecewise constant finite element space $\mathcal{V}_{0}^{h}:=\left\{v^{h}\in L^{1}(\Omega)~{}|~{}{v^{h}}_{|T}=\mbox{~{}constant},~{}\forall T\in\mathcal{T}^{h}\right\}$, problem (1.9) is then approximated in the $\left(\mathcal{V}^{h}_{1},\mathcal{V}^{h}_{0}\right)$-space according to

[TABLE]

For the numerical solution of this problem the scheme

[TABLE]

is proposed in [49, Algorithm 1], where the parameters $\tau,\theta$ are chosen suitably. Furthermore, we note that the solution $p^{h}_{n+1}$ to the sub-problem (1.13) is given by the explicit form (see [8, 49])

[TABLE]

In our framework we consider a constrained minimization problem, i.e. the linear space $\mathcal{V}^{h}_{1}$ in (1.11) is replaced by the bound-constrained set $F^{h}_{ad}$, with the forward operator being affine linear arising from solutions of elliptic PDEs, and data taken on a part of the boundary. With $\mu^{h}_{0}:=(f^{h}_{0},p^{h}_{0})\in F^{h}_{ad}\times(\mathcal{V}^{h}_{0})^{d}$ denoting the initial guess the resulting algorithm in our notation is given by

[TABLE]

with $u^{h}_{a}$ denoting the discrete adjoint state associated to $u^{h}$, see §2 for details. We mention that the sub-problem (1.14) admits a solution (cf. Remark 4.5)

[TABLE]

where $\mathcal{P}^{L^{2}}_{F^{h}_{ad}}:L^{2}(\Omega)\to F^{h}_{ad}$ denotes the $L^{2}$-projection on the set $F^{h}_{ad}$ and $-\mbox{div}$ is the adjoint operator of $\nabla$ defined by

[TABLE]

for all $q^{h}:=(q^{h}_{j})_{j=1}^{d}\in(\mathcal{V}^{h}_{0})^{d}$ and $g^{h}\in\mathcal{V}^{h}_{1}$.

We then show in §4 that the iteration sequence $\left(\mu^{h}_{n}\right)_{n}:=\left(f^{h}_{n},p^{h}_{n}\right)_{n}$ admits a subsequence $\left(\mu^{h}_{n_{k}}\right)_{k}$ which converges in the finite dimensional space $\mathcal{V}^{h}_{1}\times(\mathcal{V}^{h}_{0})^{d}$ to an element $\mu^{h}_{*}:=(f^{h}_{*},p^{h}_{*})\in F^{h}_{ad}\times\partial TV(f^{h}_{*})$, where $\mu^{h}_{*}$ satisfies the first-order optimality condition for the considered minimization problem and $f^{h}_{*}$ is a minimizer of $(\mathcal{P}^{h})$. Moreover, we derive the identity

[TABLE]

which is to be expected for primal-dual-type algorithms by the results of e.g. [42, 49]. We note that using $L^{1}-$sparsity regularization instead of the $TV-$ regularization the rate $\mathcal{O}\left(\frac{1}{n^{2}}\right)$ is shown in [43].

Source identification appears in many practical applications, like e.g. electroencephalography, geophysical prospecting and pollutant detection. Its PDE model-based treatment attracted great attention during the last decades, where we refer to e.g. [2, 7, 16, 21, 29, 34, 48] and the references given there. Let us briefly refer to publications related to our work. In [23, 38, 39] dual reciprocity boundary element methods are used to treat inverse source problems numerically. The situation where a priori knowledge of the identified source is available, e.g. as a point source, characteristic functional or a harmonic functional, is investigated in [5, 6, 11, 36, 46, 50]. Recently, by using a so-called energy functional method combined with Tikhonov regularization the authors in [31] studied a numerical method for the source identification problem from single Cauchy data. Finite element analysis of the problem of simultaneously identifying the source term and coefficients from distributed observations can be found in [45]. For the treatment of coefficient identification problems employing $TV$-regularization techniques we refer to [15, 18, 26, 27, 28, 47].

In the present work we study the finite element approximation and algorithm development for the identification problem of a source term in an elliptic PDE from partial observations of the state on the boundary. To the best of our knowledge numerical analysis for the problem setting we consider is not yet available, although there are many contributions to the numerical and algorithmic treatment of source identification problems. The main results of our paper are contained in Theorem 3.1, where we prove convergence of the finite element discrete regularized approximations to a solution of the corresponding continuous regularized problem as the mesh size of triangulations tends to zero, in Theorem 3.5, where we establish convergence to a sought source functional when the regularization parameter approaches to zero with a suitable coupling of noise level and mesh size, and in Theorem 4.9, where we prove convergence for the iterates of our discrete algorithm.

To conclude this introduction, we briefly present the space of functionals with bounded total variation, for more details the reader may consult [3, 4, 22, 24]. A scalar functional $f\in L^{1}(\Omega)$ is said to be of bounded total variation if

[TABLE]

where $\left|\cdot\right|_{\infty}$ denotes the $\ell_{\infty}$-norm on $\mathbb{R}^{d}$ defined by $\left|x\right|_{\infty}=\max\limits_{1\leq i\leq d}\left|x_{i}\right|$ and $C^{1}_{c}(\Omega)$ is the space of continuously differentiable functionals with compact support in $\Omega$. The space of all functionals in $L^{1}(\Omega)$ with bounded total variation

[TABLE]

is a non-reflexive Banach space equipped with the norm

[TABLE]

Furthermore, if $\Omega$ is an open bounded set with Lipschitz boundary, then $W^{1,1}(\Omega)\varsubsetneq BV(\Omega)$.

2 Preliminaries

Now we summarize some useful properties of the source-to-solution operator $u=u(f)$. First, we note that the decomposition

[TABLE]

holds, where $\overline{u}(f)$ and $\widehat{u}$ are the solutions to the variational equations

[TABLE]

for all $v\in H^{1}(\Omega)$, respectively. Thus, the operator $u:L^{2}(\Omega)\rightarrow H^{1}(\Omega)$ is continuously Fréchet differentiable on $L^{2}(\Omega)$. For each $f\in L^{2}(\Omega)$ the action of the Fréchet derivative in the direction $\xi\in L^{2}(\Omega)$ satisfies the equation

[TABLE]

Next, along with (1.1) we consider the adjoint problem

[TABLE]

which also admits a unique weak solution $u_{a}(f)\in H^{1}(\Omega)$ defined via the variational equation

[TABLE]

for all $v\in H^{1}(\Omega)$. Furthermore, for $f,\,\xi\in L^{2}(\Omega)$ the Fréchet differential ${u_{a}}^{\prime}(f)\xi$ is the unique weak solution in $H^{1}(\Omega)$ of the variational equation

[TABLE]

Now, for all $f\in L^{2}(\Omega)$ letting

[TABLE]

a computation for all $\xi\in L^{2}(\Omega)$ shows, by (2.5) and (2.3),

[TABLE]

and so, by (2.6),

[TABLE]

We however note that the bilinear form $L^{\prime\prime}(f)$ is in general not positive definite. In fact, considering the particular case $\Omega:=\{x=(x_{1},x_{2})\in\mathbb{R}^{2}~{}|~{}|x|=(x^{2}_{1},x^{2}_{2})^{1/2}<\frac{\sqrt{2}}{2}\}$, $\alpha$ is the unit $d\times d$-matrix, $\beta=0$ and $\sigma=0$, we have that $\xi=16-64|x|^{2}\neq 0$, but it gives $\overline{u}(\xi)=4|x|^{2}-4|x|^{4}-1$ satisfying $\overline{u}(\xi)_{|\Gamma}=0$. Therefore, the cost functional $J$ of the problem $(\mathcal{P})$ is not strictly convex. Furthermore, since the $TV$-regularization term is a semi-norm on the space $BV(\Omega)$ only, it is also not differentiable.

Lemma 2.1** ([24]).**

(i) Let $\left(f_{n}\right)_{n}$ be a bounded sequence in the $BV(\Omega)$-norm. Then a subsequence which is denoted by the same symbol and an element $f\in BV(\Omega)$ exist such that $\left(f_{n}\right)_{n}$ converges to $f$ in the $L^{1}(\Omega)$-norm.

(ii) Let $\left(f_{n}\right)_{n}$ be a sequence in $BV(\Omega)$ converging to $f$ in the $L^{1}(\Omega)$-norm. Then $f\in BV(\Omega)$ and

[TABLE]

Lemma 2.2**.**

Assume that the sequence $\left(f_{n}\right)_{n}\subset F_{ad}$ converges to an element $f$ in the $L^{1}(\Omega)$-norm. Then $f\in F_{ad}$, the sequence $\left(u(f_{n})\right)_{n}$ converges to $u(f)$ in the $H^{1}(\Omega)$-norm and $\left(u(f_{n})_{|\Gamma}\right)_{n}$ converges to $u(f)_{|\Gamma}$ in the $L^{2}(\Gamma)$-norm as well.

Proof.

First, by Lemma 2.1, we have $f\in BV(\Omega)$. Furthermore, a subsequence $\left(f_{n_{m}}\right)_{m}$ exists such that

[TABLE]

Since

[TABLE]

sending $m$ to $\infty$, we get $\underline{f}\leq f(x)\leq\overline{f}$ for a.e. in $\Omega$ which implies that $f\in F_{ad}$. Now, for all $v\in H^{1}(\Omega)$ and $m\in\mathbb{N}$ we have from (1.5) that

[TABLE]

which together with (1.4) yield

[TABLE]

that deduces $\lim_{m\to\infty}\left\|u(f_{n_{m}})-u(f)\right\|_{1,\Omega}=0$. Since the limit is unique, the whole sequence $\left(u(f_{n})\right)_{n}$ converges to $u(f)$ in the $H^{1}(\Omega)$-norm also. This completes the proof. ∎

Lemma 2.3**.**

The minimization problem $\left(\mathcal{P}\right)$ admits a minimizer.

Proof.

The assertion follows directly from Lemma 2.1 and Lemma 2.2, therefore omitted here (see [1] for details). ∎

3 Discretization and convergence

Let $\left(\mathcal{T}^{h}\right)_{0<h<1}$ be a quasi-uniform family of regular triangulations of the domain $\overline{\Omega}$ with the mesh size $h$ such that each vertex of the polygonal boundary $\partial\Omega$ is a node of $\mathcal{T}^{h}$ (cf. [13, 19]). For the definition of the discretization space of the state functionals let us denote the piecewise affine, globally continuous finite element space by

[TABLE]

where $\mathcal{P}_{1}$ consists all polynomial functionals of degree less than or equal to 1. The piecewise constant finite element space is denoted by

[TABLE]

The admissible set is discretized as

[TABLE]

Similar to the continuous case, we for each $f\in L^{2}(\Omega)$ have that the equation

[TABLE]

admits a unique solution $u^{h}:=u^{h}_{j}(f):=u^{h}(f)\in\mathcal{V}_{1}^{h}$ which further satisfies the estimate

[TABLE]

where the constant $C$ is independent of the mesh size $h$. Then the problem $\left(\mathcal{P}\right)$ can be discretized by

[TABLE]

which admits a minimizer $f^{h}$ for each $h>0$.

Also, the adjoint state $u_{a}=u_{a}(f)$ in (2.4) has the discrete version $u^{h}_{a}=u^{h}_{a}(f)$ defined via the equation

[TABLE]

for all $v^{h}\in\mathcal{V}^{h}_{1}$.

We now state the following result on the stability of finite element approximations. Here and in the sequel, unless otherwise stated, we indicate by $C$ a generic positive constant which is independent of the mesh size $h$, the regularization parameter $\rho$ and the observation data $z$.

Theorem 3.1**.**

Let $(h_{n})_{n}$ be a sequence with $\lim_{n\to\infty}h_{n}=0$ and $\rho$ be a fixed regularization parameter. Assume that $f^{h_{n}}\in F^{h_{n}}_{ad}$ is an arbitrary minimizer of $\left(\mathcal{P}^{h_{n}}\right)$ for each $n\in\mathbb{N}$. Then a subsequence of \big{(}f^{h_{n}}\big{)}_{n} not relabeled and an element $f\in F_{ad}$ exist such that

[TABLE]

Furthermore, $f$ is a solution to $\left(\mathcal{P}\right)$. If the uniqueness of the solution to $\left(\mathcal{P}\right)$ is satisfied, then convergence (3.4) holds for the whole sequence.

We mention that (3.4) yields the convergence in the $L^{s}(\Omega)$-norm for all $1\leq s<\infty$. In fact, this assertion follows directly from the following estimate

[TABLE]

To prove Theorem 3.1 we need the following crucial result.

Lemma 3.2**.**

[32, Lemma 4.6.]** For any fixed $\widehat{f}\in F_{ad}$ an element $\widehat{f}^{h}\in F^{h}_{ad}$ exists such that

[TABLE]

and

[TABLE]

Lemma 3.3**.**

Assume that $(f_{n})_{n}\subset F_{ad}$ converges to $f\in F_{ad}$ in the $L^{1}(\Omega)$-norm. Then the limit

[TABLE]

holds true.

Proof.

The proof is based on standard arguments, it is therefore omitted here. ∎

Proof of Theorem 3.1.

Let $\widehat{f}\in F_{ad}$ be arbitrary and $\widehat{f}^{h_{n}}$ be generated from $\widehat{f}$ due to Lemma 3.2. The optimality of $f^{h_{n}}$, yields

[TABLE]

By (3.2) and (3.6), we deduce from the last inequality that the sequence \big{(}f^{h_{n}}\big{)}_{n} is bounded in the $BV(\Omega)$-norm. Then, due to Lemma 2.1, a subsequence not relabeled and an element $f\in F_{ad}$ exist such that

[TABLE]

Combining this with Lemma 3.3, we arrive at

[TABLE]

here we used Lemma 3.2 (with noting $\lim_{n\to\infty}h_{n}|\log h_{n}|=0$) and Lemma 3.3 in the last equation, again. Thus, $f$ is a solution to $\left(\mathcal{P}\right)$. Furthermore, replacing $\widehat{f}$ in (3) by $f$, we obtain

[TABLE]

and, together with (3.8), arrive at

[TABLE]

and so that $TV(f)=\lim_{n\to\infty}TV(f^{h_{n}})$. The proof is completed. ∎

In the remaining part of this section we investigate the convergence of discrete regularized approximations to an identified source as $\rho$ is chosen in a suitable way depending on the mesh size and the error level of observations. Before doing so, we introduce the notion of the total variation-minimizing solution of the identification problem.

Lemma 3.4**.**

The problem

[TABLE]

admits a solution, called the total variation-minimizing solution of the identification problem.

Proof.

The assertion follows by stand arguments, therefore omitted here. ∎

Let $f\in L^{2}(\Omega)$ be fixed, we denote by

[TABLE]

Then (cf. [13, 19]),

[TABLE]

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

Theorem 3.5**.**

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

[TABLE]

where $\widehat{f}\in\mathcal{I}\left(g^{\dagger}\right)$ is arbitrary and $\chi^{h_{n}}_{\widehat{f}}$ is defined by (3.10). Furthermore, assume that $\left(z_{n}\right)_{n}\subset L^{2}(\Gamma)$ is a sequence satisfying

[TABLE]

and $f_{n}$ is an arbitrary minimizer of the problem

[TABLE]

for each $n\in N$. Then, a subsequence of $(f_{n})_{n}$ not relabeled and a total variation-minimizing solution $f^{\dagger}$ of the identification problem $\left(\mathcal{IP}\right)$ exist such that

[TABLE]

Furthermore, the discrete state sequence \big{(}u^{h_{n}}(f_{n})\big{)}_{n} converges in the $H^{1}(\Omega)$-norm to the unique weak solution $u(f^{\dagger})$ of the boundary value problem (1.1).

Proof.

By the definition of $f_{n}$, we have

[TABLE]

where $\widehat{f}^{h_{n}}\in F^{h_{n}}_{ad}$ generates from $\widehat{f}$, due to Lemma 3.2. We bound

[TABLE]

Using (1.5) and (3.1), we get for all $v^{h_{n}}\in\mathcal{V}^{h_{n}}_{1}$ that

[TABLE]

and so that

[TABLE]

Taking $v^{h_{n}}=u^{h_{n}}(\widehat{f})-u^{h_{n}}(\widehat{f}^{h_{n}})$, by (1.4), (3.5) and (3.10), we deduce

[TABLE]

and then

[TABLE]

Combining this with (3)–(3.14), we get

[TABLE]

It follows from the last inequality and (3.12), (3.6) that

[TABLE]

and

[TABLE]

Due to Lemma 2.1 and Lemma 3.3, a subsequence of $\left(f_{n}\right)_{n}$ not relabeled and an element $f^{\dagger}\in F_{ad}$ exist such that

[TABLE]

Thus, it follows from (3.17) that

[TABLE]

and so that $f^{\dagger}\in\mathcal{I}\left(g^{\dagger}\right)$. Furthermore, by (3.18) and (3.19), we also get

[TABLE]

for any $\widehat{f}\in\mathcal{I}\left(g^{\dagger}\right)$. Therefore, $f^{\dagger}$ is a total variation-minimizing solution of the identification problem. Now, by setting $\widehat{f}=f^{\dagger}$, it implies that

[TABLE]

Finally, in view of (3.1) and (1.4) we arrive at

[TABLE]

and therefore conclude that

[TABLE]

which finishes the proof. ∎

4 An iterated total variation algorithm

The aim of this section is to propose an algorithm to reach minimizers of the problem \big{(}\mathcal{P}^{h}\big{)}. We start with the following note.

Remark 4.1**.**

(i) Any $\phi\in\mathcal{V}^{h}_{1}$ can be considered as an element of ${(\mathcal{V}^{h}_{1})}^{*}$, the dual space of $\mathcal{V}^{h}_{1}$, by

[TABLE]

(ii) The inclusion $(\mathcal{V}^{h}_{0})^{d}\subset{(\mathcal{V}^{h}_{1})}^{*}$ holds via the identity

[TABLE]

Lemma 4.2**.**

For each $f^{h}\in\mathcal{V}^{h}_{1}$ the relation

[TABLE]

holds, where

[TABLE]

and for $p^{h}:=(p^{h}_{j})_{j=1}^{d}\in(\mathcal{V}^{h}_{0})^{d}$, $\|p^{h}\|_{\infty}:=\max_{j=1,...,d}\|p^{h}_{j}\|_{L^{\infty}(\Omega)}$.

Proof.

Let $f^{h}\in\mathcal{V}^{h}_{1}$ be fixed. As shown in [8] that

[TABLE]

Assume that $\mathcal{T}^{h}:=\{T_{1},T_{2},...,T_{n^{h}}\}$ with $n^{h}\in\mathbb{N}$ and $\overline{\Omega}=\cup_{i=1}^{n^{h}}\overline{T_{i}}$. We then consider $p^{*}:=(p^{*}_{ij})_{n^{h}\times d}\in\mathcal{B}_{1}\left((\mathcal{V}^{h}_{0})^{d}\right)$ with, for all $i=1,...,n^{h}$ and $j=1,...,d$,

[TABLE]

and deduce from (4.2), (4.4) that

[TABLE]

The subdifferential $\partial TV(f^{h})$ is then given by (cf. [52, Theorem 2.4.18])

[TABLE]

and so that (4.3) follows, which finishes the proof. ∎

We mention that, similar to the decomposition (2.1), for all $f\in L^{2}(\Omega)$ we have

[TABLE]

where

[TABLE]

for all $v^{h}\in\mathcal{V}^{h}_{1}$. Furthermore, the identity

[TABLE]

holds for all $f,\xi\in L^{2}(\Omega)$.

Lemma 4.3**.**

The functional $f^{h}\in F^{h}_{ad}$ is a solution of $(\mathcal{P}^{h})$ if and only if there exists $p^{h}\in\partial TV(f^{h})$ such that

[TABLE]

where $u^{h}_{a}$ is the discrete adjoint state defined by (3.3).

Proof.

The functionals $\|u^{h}(\cdot)-z\|_{\Gamma}$ and $TV(\cdot)$ are both convex on the convex set $F^{h}_{ad}$, an element $f^{h}\in F^{h}_{ad}$ is thus a solution of $(\mathcal{P}^{h})$ if and only if

[TABLE]

for all $g^{h}\in F^{h}_{ad}$. Using (3.3) and (4.7)–(4.6), we get

[TABLE]

and thus deduce from (4.10) that

[TABLE]

We here consider the convex functional

[TABLE]

where $I_{F^{h}_{ad}}$ is the indicator functional of the set $F^{h}_{ad}$. Then, the inequality (4.11) is equivalent to the relation

[TABLE]

Since, for all $\xi^{h}\in\mathcal{V}^{h}_{1}$

[TABLE]

the relation (4.12) means that there exists $p^{h}\in\partial TV(f^{h})$ such that

[TABLE]

for all $g^{h}\in F^{h}_{ad}$. This together with (4.1)–(4.2) yields (4.8). Now, for all $q^{h}\in(\mathcal{V}^{h}_{0})^{d}$ with $\|q^{h}\|_{\infty}\leq 1$ we have from the fact $p^{h}\in\partial TV(f^{h})$ defined by the equation (4.3) that

[TABLE]

which finishes the proof. ∎

Remark 4.4**.**

The system (4.8)–(4.9) is equivalent to the following inequality

[TABLE]

for all $(g^{h},q^{h})\in F^{h}_{ad}\times(\mathcal{V}^{h}_{0})^{d}$ with $\|q^{h}\|_{\infty}\leq 1$.

In the sequel, we make use the following notation (cf. [8, 10])

[TABLE]

By the inverse inequality (cf. [13, 19]) $\|v^{h}\|_{1,\Omega}\leq Ch^{-1}\|v^{h}\|_{\Omega}$ for all $v^{h}\in\mathcal{V}^{h}_{1}$, we obtain that $\|\nabla\|\leq Ch^{-1}$.

Remark 4.5**.**

Due to the optimality of $f^{h}_{n+1}$ in (1.14), we have that

[TABLE]

for all $g^{h}\in F^{h}_{ad}$ which can be rewritten as

[TABLE]

or

[TABLE]

where $\mathcal{P}^{L^{2}}_{F^{h}_{ad}}:L^{2}(\Omega)\to F^{h}_{ad}$ denotes the $L^{2}$-projection on the set $F^{h}_{ad}$ characterized by

[TABLE]

for each $\phi\in L^{2}(\Omega)$.

Remark 4.6**.**

For all $f\in L^{2}(\Omega)$ let $\overline{u^{h}_{a}}(f)\in\mathcal{V}^{h}_{1}$ be the unique weak solution of the following variational equation

[TABLE]

for all $v^{h}\in\mathcal{V}^{h}_{1}$, where $\overline{u^{h}}(f)$ is defined by (4.6). Then we have for all $v^{h}\in\mathcal{V}^{h}_{1}$ and $f,g\in L^{2}(\Omega)$ that

[TABLE]

Thus, by (1.4), taking $v^{h}=u^{h}_{a}(f)-u^{h}_{a}(g)-\overline{u^{h}_{a}}(f-g)\in\mathcal{V}^{h}_{1}$, we obtain

[TABLE]

for all $f,g\in L^{2}(\Omega)$. Furthermore, by (1.4), (1.7) and (4.6), it holds the estimate

[TABLE]

Remark 4.7**.**

Let

[TABLE]

Then the system (4.8)–(4.9) can be rewritten in the abbreviation form

[TABLE]

Furthermore, the optimality of $f^{h}_{n+1}$ and $p^{h}_{n+1}$ due to Algorithm 1 yields that

[TABLE]

for all $\nu^{h}=(g^{h},q^{h})\in F^{h}_{ad}\times\mathcal{B}_{1}\left((\mathcal{V}^{h}_{0})^{d}\right)$.

We have the following auxiliary result.

Lemma 4.8**.**

Let $\mu^{h}_{n}:=\left(f^{h}_{n},p^{h}_{n}\right)$ and $\mu^{h}_{n+1}:=\left(f^{h}_{n+1},p^{h}_{n+1}\right)$. Then the system (4.21) has the form

[TABLE]

where the matrix

[TABLE]

is symmetric, positive definite.

Proof.

Using the identity (4.18), we rewrite the system (4.21) in the form

[TABLE]

which verifies the inequality (4.22), where the matrix $B$ by the equation (1.17) is symmetric. Now we show that $B$ is positive definite. For all $\mu^{h}=(f^{h},p^{h})\in F^{h}_{ad}\times\mathcal{B}_{1}\left((\mathcal{V}^{h}_{0})^{d}\right)$, by the inequality (4.19), we get

[TABLE]

This together with (4.14) yields that

[TABLE]

with

[TABLE]

satisfying the relation $PR>Q^{2}$ due to (4.15), which finishes the proof. ∎

Theorem 4.9**.**

For any fixed $h>0$ let $\left(\mu^{h}_{n}\right)_{n}:=\left(f^{h}_{n},p^{h}_{n}\right)_{n}$ be the sequence generated by Algorithm 1.

(i) Then a subsequence $\left(\mu^{h}_{n_{k}}\right)_{k}$ and an element $\mu^{h}_{*}:=(f^{h}_{*},p^{h}_{*})\in F^{h}_{ad}\times\partial TV(f^{h}_{*})$ exist such that $\left(\mu^{h}_{n_{k}}\right)_{k}$ converges to $\mu^{h}_{*}$ in the finite dimensional space $\mathcal{V}^{h}_{1}\times(\mathcal{V}^{h}_{0})^{d}$. Furthermore, $\mu^{h}_{*}$ satisfies the condition (4.20), i.e. $\left(A(\mu^{h}_{*}),\nu^{h}-\mu^{h}_{*}\right)\geq 0$ for all $\nu^{h}\in F^{h}_{ad}\times\mathcal{B}_{1}\left((\mathcal{V}^{h}_{0})^{d}\right)$ and $f^{h}_{*}$ is thus a minimizer of $(\mathcal{P}^{h})$.

(ii) The equation

[TABLE]

holds true, where $\|\cdot\|$ is some norm of the finite dimensional space $\mathcal{V}^{h}_{1}\times(\mathcal{V}^{h}_{0})^{d}$.

Proof.

(i) Denote $\mu^{h}:=(f^{h},p^{h})$, where $f^{h}\in F^{h}_{ad}$ is an arbitrary solution to $(\mathcal{P}^{h})$ and $p^{h}\in\partial TV(f^{h})$ generates from $f^{h}$, due to Lemma 4.3. Using the notations (4.20) and (4.23), we have that

[TABLE]

Talking $\nu^{h}=\mu^{h}$ in (4.22), we get

[TABLE]

and so that (4) yields

[TABLE]

Talking $\nu^{h}=\mu^{h}_{n+1}$ in (4.20), we have

[TABLE]

On the other hand, by (1.17), we get

[TABLE]

Combining this with (4.27)–(4), we arrive at

[TABLE]

Then for all $n\in\mathbb{N}$ it follows from (4.29) that

[TABLE]

Therefore, the sequence ${(\mu^{h}_{n})}_{n}$ is bounded while the series $\sum_{n=0}^{\infty}\|\mu^{h}_{n+1}-\mu^{h}_{n}\|^{2}_{B}$ is convergent. A subsequence of ${(\mu^{h}_{n})}_{n}$ not relabeled and an element $\mu^{h}_{*}:=(f^{h}_{*},p^{h}_{*})\in F^{h}_{ad}\times\mathcal{B}_{1}\left((\mathcal{V}^{h}_{0})^{d}\right)$ then exist such that

[TABLE]

Thus sending $n$ to $\infty$ in (4.22), we obtain

[TABLE]

for all $\nu^{h}\in F^{h}_{ad}\times\mathcal{B}_{1}\left((\mathcal{V}^{h}_{0})^{d}\right)$. As a direct consequence of (4.31) (cf. Remark 4.4) we get

[TABLE]

and so that

[TABLE]

This together with Lemma 4.2 yields $p^{h}_{*}\in\partial TV(f^{h}_{*})$. Consequently, $f^{h}_{*}$ is then a minimizer to $(\mathcal{P}^{h})$, due to (4.31) and Lemma 4.3.

(ii) It remains to show (4.24). For all $n\in\mathbb{N}$ we get from (4.30) that

[TABLE]

On the other hand, in view of (4.22) we have

[TABLE]

which implies

[TABLE]

or

[TABLE]

In view of (4) we get

[TABLE]

which together with

[TABLE]

imply that

[TABLE]

and so that

[TABLE]

for all $n\in\mathbb{N}$. Therefore, we arrive at

[TABLE]

by (4.32). The proof is completed. ∎

5 Numerical test

We now illustrate the theoretical result with numerical examples. For simplicity of exposition we restrict ourselves to the case where $\beta=\sigma=0$ in the system (1.1), and thus consider the Neumann problem

[TABLE]

with $\Omega=\{x=(x_{1},x_{2})\in{R}^{2}~{}|~{}-1<x_{1},x_{2}<1\}$. Let $\chi_{D}$ denote the characteristic functional of a Lebesgue measurable set $D\subset\mathbb{R}^{2}$. We assume that entries of the known symmetric diffusion matrix $\alpha$ are defined as

[TABLE]

where

[TABLE]

see Figure 1 (bottom row, right). The sought source functional $f^{\dagger}$ in (5.1) is assumed to be discontinuous and defined as

[TABLE]

We adopt the mesh data structure of [25] to our numerical implementation. For this purpose we divide the interval $(-1,1)$ into $\ell$ equal segments, so that the domain $\Omega=(-1,1)^{2}$ is divided into $2\ell^{2}$ triangles, where the diameter of each triangle is $h=h_{\ell}=\frac{\sqrt{8}}{\ell}$. The Neumann boundary data $j$ in the equation (5.2) is chosen as

[TABLE]

We note that (5.1)–(5.2) is the pure Neumann boundary value problem, so $f^{\dagger}$ and $j$ were taken such that the compatibility condition $\int_{\Omega}f^{\dagger}+\int_{\partial\Omega}j=0$ is satisfied. Let $u_{j}(f^{\dagger})$ denote the unique weak solution of (5.1)–(5.2). The Dirichlet boundary data on $\Gamma$ of the problem (5.1)–(5.2) is computed by

[TABLE]

We use Algorithm 1 which is described in Section 4 for computing the numerical solution of the problem $\left(\mathcal{P}^{h}\right)$. Before doing so, we discuss the constants $c_{\gamma}$ and $c_{1}$ appearing in (4.15). First, we can show the inequality (see, e.g., [33, §7])

[TABLE]

for all $u\in H^{1}(\Omega)$. Therefore, the constant $c_{1}$ can be chosen by $\frac{\underline{\alpha}}{1+\left(\sqrt{\frac{3}{2}}\right)^{(d+2)/2}|\Omega|^{1/d}}$. Furthermore, in the case $\Omega=(-1,1)^{2}$, the following estimate provides an upper bound for $c_{\gamma}$.

Lemma 5.1**.**

If $\Omega:=(a_{1},b_{1})\times(a_{2},b_{2})\times\ldots\times(a_{d},b_{d})\ni 0$, then $c_{\gamma}$ can be chosen by $\sqrt{\frac{d+|\overline{\hbar}|^{2}}{|\underline{\hbar}|}}$, i.e. for all $u\in H^{1}(\Omega)$, the estimate

[TABLE]

holds true, where $\underline{\hbar}=\min\{|a_{1}|,|b_{1}|,\ldots,|a_{d}|,|b_{d}|\}$ and $\overline{\hbar}=\max\{|a_{1}|,|b_{1}|,\ldots,|a_{d}|,|b_{d}|\}$.

Proof.

Denote by \underline{S}_{d}:=\left\{\big{(}x^{\prime},\underline{\varphi}(x^{\prime})\big{)}\in\overline{\Omega}~{}\big{|}~{}x^{\prime}:=(x_{1},\ldots,x_{d-1})~{}\mbox{and}~{}\underline{\varphi}(x^{\prime})\equiv a_{d}\right\}. First, we assume that $u\in C^{1}(\overline{\Omega})$. Then, we get

[TABLE]

which implies that

[TABLE]

Multiplying this inequality by $\sqrt{1+\left(\frac{\partial\underline{\varphi}}{x_{1}}\right)^{2}+\ldots+\left(\frac{\partial\underline{\varphi}}{x_{d-1}}\right)^{2}}\equiv 1$, integrating over $Q:=(a_{1},b_{1})\times(a_{2},b_{2})\times\ldots\times(a_{d-1},b_{d-1})$, we obtain

[TABLE]

Likewise, denoting by \overline{S}_{d}:=\left\{\big{(}x^{\prime},\overline{\varphi}(x^{\prime})\big{)}\in\overline{\Omega}~{}\big{|}~{}x^{\prime}:=(x_{1},\ldots,x_{d-1})~{}\mbox{and}~{}\overline{\varphi}(x^{\prime})\equiv b_{d}\right\}, we also get

[TABLE]

and arrive at

[TABLE]

This implies that (5.5) is satisfied for all $u\in C^{1}(\overline{\Omega})$. By the everywhere dense property of the set $C^{1}(\overline{\Omega})$ in $H^{1}(\Omega)$, we conclude that (5.5) is also fulfilled for all $u\in H^{1}(\Omega)$, which finishes the proof. ∎

The constant $\underline{\alpha}$ in (1.2) is taken by $0.1$. With $d=2$, $\Omega=(-1,1)^{2}$ we then can take $c_{1}=0.025$ and $c_{\gamma}=\sqrt{3}$. As an initial approximation we choose $f^{h_{\ell}}_{0}(x)\equiv 1\in\mathcal{V}^{h_{\ell}}_{1}$ and $p^{h_{\ell}}_{0}(x)\equiv(\frac{1}{2},\ldots,\frac{1}{2})\in(\mathcal{V}^{h_{\ell}}_{0})^{d}$. In the computation below we choose $\tau=2\cdot 10^{-4}$ and $\theta=5\cdot 10^{-2}$.

For observations with noise we assume that

[TABLE]

where $R_{g^{\dagger}}$ is the matrix of random numbers on the interval $(-1,1)$. We here refer to the discussion in [20] which states that the $L^{\infty}$-misfit term is well suited for this uniformly distributed noise model. However, since $L^{\infty}$ is a non-reflexive Banach space and the $L^{\infty}$-norm is not differentiable the use of this misfit term would cause difficulties in the numerical treatment of ill-posed inverse problems.

The measurement error is then computed as \delta_{\ell}=\big{\|}z_{\delta_{\ell}}-g^{\dagger}\big{\|}_{L^{2}(\Gamma)}. To satisfy the condition (3.12) in Theorem 3.5 we take $h=h_{\ell}$, $\rho=\rho_{\ell}=10^{-3}h^{1/2}_{\ell}$ and $\theta_{\ell}=h_{\ell}\rho^{1/2}_{\ell}$.

Starting with the coarsest level $\ell=4$, we follow [35, pp. 16] in the use of the stopping criterion. Here at each iteration $k$ we compute

[TABLE]

with $\tau_{1}:=10^{-5}h_{\ell}^{1/2}$ and $\tau_{2}:=10^{-4}h_{\ell}^{1/2}$ and the iteration was stopped if $\mbox{Tolerance}\leq 0$ or the number of iterations reached the maximum iteration count of $N=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^{h_{\ell}}_{0},p^{h_{\ell}}_{0})$ for the algorithm on this finer mesh, and proceed similarly on the refinement levels.

Let $f_{\ell}$ denote the functional which is obtained at the final iterate of Algorithm 1 corresponding to the refinement level $\ell$. Furthermore, we denote by $u_{\ell}$ the computed numerical solution to the Dirichlet problem

[TABLE]

while by $u^{\dagger}$ the solution of the equation (5.1) supplemented with the Dirichlet boundary condition \big{(}u_{j}(f^{\dagger})\big{)}_{|\partial\Omega}.

First, we consider the case where $\Gamma$ is the bottom surface of the domain $\Omega$, i.e. $\Gamma=\{x=(x_{1},x_{2})\in\overline{\Omega}~{}|~{}x_{2}=-1\}$. The numerical results of this case 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 \big{\|}f^{\dagger}-f_{\ell}\big{\|}_{\Omega}, \big{\|}u^{\dagger}-u_{\ell}\big{\|}_{\Omega} and \big{\|}u^{\dagger}-u_{\ell}\big{\|}_{1,\Omega}.

We now consider the case where $\Gamma$ includes the bottom surface and the left surface of the domain $\Omega$, i.e. $\Gamma=\{x=(x_{1},x_{2})\in\overline{\Omega}~{}|~{}x_{2}=-1\}\cup\{x=(x_{1},x_{2})\in\overline{\Omega}~{}|~{}x_{1}=-1\}$. To try obtaining an error level similar to the case where observations were taken on the bottom surface only, $\theta_{\ell}$ in (5.6) is changed to $\theta_{\ell}=\frac{1}{2}h_{\ell}\rho^{1/2}_{\ell}$.

Computations show that the measurement noise $\delta_{\ell}=4.6039.10^{-4}$, $\mbox{Tolerance}=-1.8493.10^{-8}$, the iteration stops at $n=477$. Errors \big{\|}f^{\dagger}-f_{\ell}\big{\|}_{\Omega}=6.1605.10^{-2}, \big{\|}u^{\dagger}-u_{\ell}\big{\|}_{\Omega}=5.8110.10^{-4} and \big{\|}u^{\dagger}-u_{\ell}\big{\|}_{1,\Omega}=1.1164.10^{-2}.

Acknowledgements

The authors M. Hinze and T.N.T. Quyen would like to thank the Referee and the Editor for their valuable comments and suggestions which helped to improve our paper.