Finite element analysis for identifying the reaction coefficient in PDE from boundary observations

TL;DR

This paper develops a finite element method with Tikhonov regularization to accurately identify the reaction coefficient in an elliptic PDE from boundary measurements, ensuring stability and convergence.

Contribution

It introduces a novel finite element approach combined with Tikhonov regularization for the inverse problem of reaction coefficient identification in elliptic PDEs.

Findings

Method achieves stable and convergent approximations.

Numerical experiments confirm theoretical results.

Abstract

This work is devoted to the nonlinear inverse problem of identifying the reaction coefficient in an elliptic boundary value problem from single Cauchy data on a part of the boundary. We then examine simultaneously two elliptic boundary value problems generated from the available Cauchy data. The output least squares method with the Tikhonov regularization is applied to find approximations of the sought coefficient. We discretize the PDEs with piecewise linear finite elements. The stability and convergence of this technique are then established. 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} , regularization parameter ρ ℓ subscript 𝜌 ℓ \rho_{\ell} , measurement noise δ ℓ subscript 𝛿 ℓ \delta_{\ell} , and errors L β 2 subscript superscript 𝐿 2 𝛽 L^{2}_{\beta} , L N 2 subscript superscript 𝐿 2 𝑁 L^{2}_{N} , L M 2 subscript superscript 𝐿 2 𝑀 L^{2}_{M} , L D 2 subscript superscript 𝐿 2 𝐷 L^{2}_{D} .

Error history for ${\Gamma }_{\text{observation}}:={\Gamma }_{\text{bottom}}$
$\ell$	$h_{\ell }$	${\rho }_{\ell }$	${\delta }_{\ell }$	$L_{\beta }^2$	$L_N^2$	$L_M^2$	$L_D^2$
4	0.7071	8.4090e-4	0.1116	1.2185	0.3026	0.2674	0.1289
8	0.3536	5.9460e-4	4.0042e-2	0.5989	0.1377	0.1205	9.4845e-2
16	0.1767	4.2045e-4	1.9375e-2	0.2378	0.1014	7.7931e-2	3.6737e-2
32	8.8388e-2	2.9730e-4	7.9021e-3	0.1472	4.6465e-2	3.5284e-2	1.6411e-2
64	4.4194e-2	2.1022e-4	3.2289e-3	7.4264e-2	2.0169e-2	1.6636e-2	6.9882e-3

Table 2. Table 2: Experimental order of convergence between finest and coarsest level for L β 2 subscript superscript 𝐿 2 𝛽 L^{2}_{\beta} , L N 2 subscript superscript 𝐿 2 𝑁 L^{2}_{N} , L M 2 subscript superscript 𝐿 2 𝑀 L^{2}_{M} and L D 2 subscript superscript 𝐿 2 𝐷 L^{2}_{D} .

Experimental order of convergence (${\Gamma }_{\text{observation}}:={\Gamma }_{\text{bottom}}$)
$\ell$	EOC${}_{L_{\beta }^2}$	EOC${}_{L_N^2}$	EOC${}_{L_M^2}$	EOC${}_{L_D^2}$
4	–	–	–	–
8	1.0247	1.1359	1.1500	0.4426
16	1.3326	0.4415	0.6288	1.3683
32	0.6920	1.1258	1.1432	1.1626
64	0.9870	1.2040	1.0847	1.2317
Mean of EOC	1.0091	0.9768	1.0017	1.0513

Table 3. Table 3: Regularization parameter ρ ℓ subscript 𝜌 ℓ \rho_{\ell} , measurement noise δ ℓ subscript 𝛿 ℓ \delta_{\ell} and errors L β 2 subscript superscript 𝐿 2 𝛽 L^{2}_{\beta} , L N 2 subscript superscript 𝐿 2 𝑁 L^{2}_{N} , L M 2 subscript superscript 𝐿 2 𝑀 L^{2}_{M} , L D 2 subscript superscript 𝐿 2 𝐷 L^{2}_{D} .

Regularization parameter	${\delta }_{\ell }$	$L_{\beta }^2$	$L_N^2$	$L_M^2$	$L_D^2$
${\rho }_{\ell }=h_{\ell }^2$e-0 = $1.9531$e-3	3.4240e-5	5.0702	1.5779	1.1618	0.1805
${\rho }_{\ell }=h_{\ell }^2$e-1 = $1.9531$e-4	9.3133e-6	2.4302	1.1983	0.8792	0.1439
${\rho }_{\ell }=h_{\ell }^2$e-2 = $1.9531$e-5	3.3618e-6	0.1028	3.6218e-2	2.9314e-2	1.5681e-2
${\rho }_{\ell }=h_{\ell }^2$e-3 = $1.9531$e-6	1.0207e-6	8.1651e-2	3.4587e-2	2.7745e-2	1.3142e-2

Table 4. Table 4: Refinement level ℓ ℓ \ell , measurement noise δ ℓ subscript 𝛿 ℓ \delta_{\ell} , and errors L β 2 subscript superscript 𝐿 2 𝛽 L^{2}_{\beta} , L N 2 subscript superscript 𝐿 2 𝑁 L^{2}_{N} , L M 2 subscript superscript 𝐿 2 𝑀 L^{2}_{M} , L D 2 subscript superscript 𝐿 2 𝐷 L^{2}_{D} .

Error history for ${\Gamma }_{\text{observation}}:={\Gamma }_{\text{bottom-left}}$
$\ell$	${\delta }_{\ell }$	$L_{\beta }^2$	$L_N^2$	$L_M^2$	$L_D^2$
4	7.6402e-2	1.1379	0.3066	0.2618	0.1911
8	3.2528e-2	0.5102	0.2215	0.1834	0.1125
16	1.4637e-2	0.2056	0.1443	0.1196	5.2079e-2
32	5.4159e-3	0.1192	5.6266e-2	4.2652e-2	2.0737e-2
64	2.3331e-3	6.6896e-2	1.4888e-2	1.0937e-2	6.1714e-3

Table 5. Table 5: Experimental order of convergence between finest and coarsest level for L β 2 subscript superscript 𝐿 2 𝛽 L^{2}_{\beta} , L N 2 subscript superscript 𝐿 2 𝑁 L^{2}_{N} , L M 2 subscript superscript 𝐿 2 𝑀 L^{2}_{M} and L D 2 subscript superscript 𝐿 2 𝐷 L^{2}_{D} .

Experimental order of convergence (${\Gamma }_{\text{observation}}:={\Gamma }_{\text{bottom-left}}$)
$\ell$	EOC${}_{L_{\beta }^2}$	EOC${}_{L_N^2}$	EOC${}_{L_M^2}$	EOC${}_{L_D^2}$
4	–	–	–	–
8	1.1572	0.4691	0.5135	0.7644
16	1.3112	0.6182	0.6168	1.1112
32	0.7865	1.3587	1.4875	1.3285
64	0.8334	1.9181	1.9634	1.7485
Mean of EOC	1.0221	1.0910	1.1453	1.2381

Table 6. Table 6: Numerical result for ℓ = 64 ℓ 64 \ell=64 , ρ = ρ ℓ = 0.001 ​ h ℓ 𝜌 subscript 𝜌 ℓ 0.001 subscript ℎ ℓ \rho=\rho_{\ell}=0.001\sqrt{h_{\ell}} , θ = 0.2 𝜃 0.2 \theta=0.2 , i.e. δ ℓ = 0.4773 subscript 𝛿 ℓ 0.4773 \delta_{\ell}=0.4773 , and with multiple measurements I = 1 , 6 , 16 𝐼 1 6 16 I=1,6,16 .

Numerical result for $\ell =64$, $\theta =0.2$
Number of measurements $I$	Iterate	$L_{\beta }^2$	$L_N^2$	$L_M^2$	$L_D^2$
1	600	0.5889	0.3040	0.2397	0.1338
6	600	0.3846	0.1597	0.1215	6.2670e-2
16	600	0.2746	8.7720e-2	7.1346e-2	4.5228e-2

Equations

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

$$\mathbf{\alpha }\nabla \Phi \cdot \vec{\mathbf{n}}+\sigma \Phi$$

$$\Phi$$

$$\mathbf{\alpha }(x)\xi \cdot \xi =\sum _{1\le r,s\le d}{\alpha }_{rs}(x){\xi }_r{\xi }_s\ge \underline{\alpha }|\xi {|}^2$$

$$\beta \in {\mathcal{S}}_{ad}:=\{\beta \in L^{\infty }(\Omega )|0<\underline{\beta }\le \beta (x)\le \overline{\beta }\text{mbox}a.e.in\Omega \}$$

$$\displaystyle\frac{1}{I}\sum_{i=1}^{I}\left(\big{\|}j_{\delta}^{i}-j^{\dagger}\big{\|}_{H^{-1/2}(\Gamma)}+\big{\|}g_{\delta}^{i}-g^{\dagger}\big{\|}_{H^{1/2}(\Gamma)}\right)\leq\delta$$

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

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

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

$$v$$

$$\mathbf{\alpha }\nabla v\cdot \vec{\mathbf{n}}+\sigma v$$

$$\displaystyle\beta^{h}_{\delta,\rho}(x)=\mathcal{P}_{[\underline{\beta},\overline{\beta}]}\left(\frac{1}{\rho}\big{(}N^{h}_{j_{\delta}}(\beta^{h}_{\delta,\rho})(x)A^{h}_{N}(\beta^{h}_{\delta,\rho})(x)-M^{h}_{g_{\delta}}(\beta^{h}_{\delta,\rho})(x)A^{h}_{M}(\beta^{h}_{\delta,\rho})(x)\big{)}+\beta^{*}(x)\right)$$

$${\beta }^{†}=\arg \underset{\left\{\beta \in {\mathcal{S}}_{ad}|N_{j^{†}}(\beta )=M_{g^{†}}(\beta )\right\}}{\min }\Vert \beta -{\beta }^{\ast }{\Vert }_{L^2(\Omega )}.$$

$$[u,v]:=[u,v{]}_{(\mathbf{\alpha },\beta ,\sigma )}:=(\mathbf{\alpha }\nabla u,\nabla v{)}_{\Omega }+(\beta u,v{)}_{\Omega }+(\sigma u,v{)}_{\partial \Omega }$$

$$c_1\Vert u{\Vert }_{1,\Omega }\le [u,u{]}_{(\mathbf{\alpha },\beta ,\sigma )}\le c_2\Vert u{\Vert }_{1,\Omega }$$

$$\displaystyle\big{[}N_{j_{\delta}}(\beta),\phi\big{]}_{(\boldsymbol{\alpha},\beta,\sigma)}=\langle f,\phi\rangle_{\Omega}+\langle j_{\delta},\phi\rangle_{\Gamma}+\langle j_{0},\phi\rangle_{\partial\Omega\setminus\Gamma}$$

$$\displaystyle\big{\|}N_{j_{\delta}}(\beta)\big{\|}_{1,\Omega}\preceq\|j_{\delta}\|_{H^{-1/2}(\Gamma)}+\|j_{0}\|_{H^{-1/2}(\partial\Omega\setminus\Gamma)}+\|f\|_{H^{-1}(\Omega)}.$$

$$\displaystyle\big{[}M_{g_{\delta}}(\beta),\phi\big{]}_{(\boldsymbol{\alpha},\beta,\sigma)}=\langle f,\phi\rangle_{\Omega}+\langle j_{0},\phi\rangle_{\partial\Omega\setminus\Gamma}$$

$$\Vert M_{g_{\delta }}(\beta ){\Vert }_{1,\Omega }⪯\Vert g_{\delta }{\Vert }_{H^{1/2}(\Gamma )}+\Vert j_0{\Vert }_{H^{-1/2}(\partial \Omega \setminus \Gamma )}+\Vert f{\Vert }_{H^{-1}(\Omega )}.$$

$$\Vert N_{j_{\delta }}(\beta ){\Vert }_{2,\Omega }⪯\Vert j_{\delta }{\Vert }_{H^{1/2}(\partial \Omega )}+\Vert j_0{\Vert }_{H^{1/2}(\partial \Omega \setminus \Gamma )}+\Vert f{\Vert }_{L^2(\Omega )}$$

$$\Vert M_{g_{\delta }}(\beta ){\Vert }_{2,\Omega }⪯\Vert g_{\delta }{\Vert }_{H^{3/2}(\partial \Omega )}+\Vert j_0{\Vert }_{H^{1/2}(\partial \Omega \setminus \Gamma )}+\Vert f{\Vert }_{L^2(\Omega )}.$$

$$N_{j_{\delta }},M_{g_{\delta }}:{\mathcal{S}}_{ad}\to H^1(\Omega ).$$

$$\displaystyle\big{[}D_{N}^{(m)},\phi\big{]}_{(\boldsymbol{\alpha},\beta,\sigma)}=-\sum_{i=1}^{m}\left(\kappa_{i}N_{j_{\delta}}^{(m-1)}(\beta)\xi_{i},\phi\right)_{\Omega},\enskip\forall\phi\in H^{1}(\Omega)$$

$$\displaystyle\big{[}D_{M}^{(m)},\phi\big{]}_{(\boldsymbol{\alpha},\beta,\sigma)}=-\sum_{i=1}^{m}\left(\kappa_{i}M_{g_{\delta}}^{(m-1)}(\beta)\xi_{i},\phi\right)_{\Omega},\enskip\forall\phi\in H^{1}_{0}(\Omega\cup\Gamma)$$

$$\displaystyle\max\left(\big{\|}D_{N}^{(m)}\big{\|}_{H^{1}(\Omega)},\big{\|}D_{M}^{(m)}\big{\|}_{H^{1}(\Omega)}\right)\preceq\prod_{i=1}^{m}\|\kappa_{i}\|_{L^{2}(\Omega)}.$$

$${\beta }_n\rightharpoonup \beta \text{mbox}{weakly}^{\ast }\text{mbox}in{L}^{\infty }(\Omega ),\text{mbox}i.e.({\beta }_n,\xi {)}_{\Omega }\to (\beta ,\xi {)}_{\Omega }\text{mbox}forall\xi \in L^1(\Omega ),$$

$$N_n\rightharpoonup {\theta }_N\text{mbox}weaklyin{H}^1(\Omega )\text{mbox}andstronglyin{L}^2(\Omega )\text{mbox}and{L}^2(\partial \Omega )$$

$$⟨f,\varphi {⟩}_{\Omega }+⟨j_{\delta },\varphi {⟩}_{\Gamma }+⟨j_0,\varphi {⟩}_{\partial \Omega \setminus \Gamma }$$

$$=(\mathbf{\alpha }\nabla N_n,\nabla \varphi {)}_{\Omega }+({\beta }_n{N}_n,\varphi {)}_{\Omega }+(\sigma N_n,\varphi {)}_{\partial \Omega }$$

$$=(\mathbf{\alpha }\nabla {\theta }_N,\nabla \varphi {)}_{\Omega }+(\beta {\theta }_N,\varphi {)}_{\Omega }+(\sigma {\theta }_N,\varphi {)}_{\partial \Omega }$$

$$+(\mathbf{\alpha }\nabla (N_n-{\theta }_N),\nabla \varphi {)}_{\Omega }+({\beta }_n-\beta ,{\theta }_N\varphi {)}_{\Omega }$$

$$+({\beta }_n(N_n-{\theta }_N),\varphi {)}_{\Omega }+(\sigma (N_n-{\theta }_N),\varphi {)}_{\partial \Omega }.$$

Full text

{AMS}

35R25, 47A52, 35R30, 65J20, 65J22

1 Introduction

Let $\Omega$ be an open, bounded and connected domain of $\mathbb{R}^{d},~{}d\geq 2$ with the boundary $\partial\Omega$ and $\Gamma\subset\partial\Omega$ be an accessible part of the boundary which is relatively open. In this paper we are related with the following elliptic system

[TABLE]

where $\vec{\boldsymbol{n}}$ is the unit outward normal on $\partial\Omega$, the boundary conditions $j^{\dagger}\in H^{-1/2}(\Gamma):={H^{1/2}(\Gamma})^{*}$, $j_{0}\in H^{-1/2}(\partial\Omega\setminus\Gamma)$,  $g^{\dagger}\in H^{1/2}(\Gamma)$, the source term $f\in H^{-1}(\Omega):={H^{1}(\Omega)}^{*}$ and the functions $\boldsymbol{\alpha},~{}\sigma$ are assumed to be known. Here, $\sigma\in L^{\infty}(\partial\Omega)$ with $\sigma(x)\geq 0$ a.e. on $\partial\Omega$ and $\boldsymbol{\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]

a.e. in $\Omega$ for all $\xi=\left(\xi_{r}\right)_{1\leq r\leq d}\in\mathbb{R}^{d}$ with some constant $\underline{\alpha}>0$. In case $\boldsymbol{\alpha}=\alpha\cdot\boldsymbol{I}_{d}$, the unit $d\times d$-matrix $\boldsymbol{I}_{d}$ and $\alpha:\Omega\to\mathbb{R}$, then $\alpha$ is called the scalar diffusion.

The system (1)–(4) is overdetermined, i.e. on the part $\Gamma$ of the boundary $\partial\Omega$ the Neumann data and Dirichlet data are supplemented simultaneously. Therefore, if the reaction coefficient

[TABLE]

is given also, there may be no $\Phi$ satisfying the system, where the constants $0<\underline{\beta}\leq\overline{\beta}$ are known. In this paper we assume that the system is consistent and our aim is to reconstruct the coefficient $\beta\in\mathcal{S}_{ad}$ from several sets of observation data \big{(}j_{\delta}^{i},g_{\delta}^{i}\big{)}_{i=1,\ldots I}\subset H^{-1/2}(\Gamma)\times H^{1/2}(\Gamma) of the exact \big{(}j^{\dagger},g^{\dagger}\big{)} obeying the deterministic noise model

[TABLE]

with $\delta>0$ denoting the error level of the observations.

The problem arises from different contexts of applied sciences, e.g., from aquifer analysis, optical tomography which attracted great attention of many scientists in the last 40 years or so. For surveys on the subject, we refer the reader to [5, 41, 42]. Although there have been many papers devoted to the subject, the authors however used the distributed observations, i.e. the measurement data is assumed to be given in the whole domain $\Omega$, see, e.g., Alt [3], Colonius and Kunisch [10], Engl et al. [13], Kaltenbacher and Hofmann [25], Kaltenbacher and Klassen [26], Neubauer [33], Resmerita and Scherzer [36] and [16, 17, 18, 19] and the references therein. We mention that the boundary observation subject, i.e. the measurement data is available only on the part $\Gamma$ of the boundary $\partial\Omega$ as defined in the present paper for this identification problem, which is more realistic from the practical point of view, has not yet been investigated so far.

For simplicity of exposition we below consider one observation pair $\left(j_{\delta},g_{\delta}\right)$ being available, i.e. $I=1$, while the approach described here can be naturally extended to multiple measurements. We from the available observation data $\left(j_{\delta},g_{\delta}\right)$ simultaneously examine the Neumann boundary value problem

[TABLE]

and the mixed boundary value problem

[TABLE]

Let $N_{j_{\delta}}(\beta)$ and $M_{g_{\delta}}(\beta)$ be the unique weak solutions of Eq. 7–Eq. 9 and Eq. 10–Eq. 12, respectively. We then consider a minimizer $\beta_{\delta,\rho}$ of the Tikhonov regularized minimization problem

[TABLE]

as reconstruction, where $\rho>0$ is the regularization parameter and $\beta^{*}$ is an a priori estimate of the true coefficient. The motivation for using the above cost functional is that \big{\|}N_{j_{\delta}}(\beta)-M_{g_{\delta}}(\beta)\big{\|}_{L^{2}(\Omega)}\geq 0 and at the sought coefficient $\beta$ it holds the identity \big{\|}N_{j^{\dagger}}(\beta)-M_{g^{\dagger}}(\beta)\big{\|}_{L^{2}(\Omega)}=0.

Let $N^{h}_{j_{\delta}}(\beta)$ and $M^{h}_{g_{\delta}}(\beta)$ be corresponding approximations of $N_{j_{\delta}}(\beta)$ and $M_{g_{\delta}}(\beta)$ in the finite dimensional space $\mathcal{V}^{h}_{1}$ of piecewise linear, continuous finite elements. Utilizing the variational discretization concept [21] of $\left(\mathcal{P}_{\delta,\rho}\right)$ that avoids explicit discretization of the control variable, we then consider the discrete regularized problem corresponding to $\left(\mathcal{P}_{\delta,\rho}\right)$, i.e. the following minimization problem

[TABLE]

which also attains a minimizer $\beta^{h}_{\delta,\rho}$ satisfying the relation (cf. Section 2)

[TABLE]

a.e. in $\Omega$, where \mathcal{P}_{[\underline{\beta},\overline{\beta}]}(c):=\max\left(\underline{\beta},\min\big{(}c,\overline{\beta}\big{)}\right), the states $A^{h}_{N}$ and $A^{h}_{M}$ are finite $\mathcal{V}^{h}_{1}$-element approximations of solutions to suitably chosen adjoint problems. This identity will be exploited to the gradient projection algorithm presented in Section 4.

In Section 3 we show that the proposed finite element method is stable, i.e. if the regularization parameter and the observation data are both fixed, then the sequence of minimizers \big{(}\beta^{h}_{\delta,\rho}\big{)}_{h>0} to \big{(}\mathcal{P}^{h}_{\delta,\rho}\big{)} can be extracted a subsequence which converges in the $L^{2}(\Omega)$-norm to a solution of \big{(}\mathcal{P}_{\delta,\rho}\big{)} as the mesh size $h$ of the triangulation $\mathcal{T}^{h}$ tends to zero. Furthermore as $h,\delta\to 0$ and with an appropriate a priori regularization parameter choice $\rho=\rho(h,\delta)\to 0$, the whole sequence \big{(}\beta^{h}_{\delta,\rho}\big{)}_{\rho>0} converges in the $L^{2}(\Omega)$-norm to the $\beta^{*}$-minimum-norm solution $\beta^{\dagger}$ of the identification problem defined by

[TABLE]

The corresponding state sequences \big{(}N^{h}_{j_{\delta}}\big{(}\beta^{h}_{\delta,\rho}\big{)}\big{)}_{\rho>0} and \big{(}M^{h}_{g_{\delta}}\big{(}\beta^{h}_{\delta,\rho}\big{)}\big{)}_{\rho>0} then converge in the $H^{1}(\Omega)$-norm to the exact state $\Phi^{\dagger}=\Phi(j^{\dagger},g^{\dagger},\beta^{\dagger})$ of the problem Eq. 1–Eq. 4.

Our numerical implementation will be presented in Section 4. First, for the numerical solution of the discrete regularized problem \big{(}\mathcal{P}^{h}_{\delta,\rho}\big{)} we employ a gradient projection algorithm with Armijo steplength rule. In Section 4.2 we assume that observations are available on the bottom surface of the domain. Section 4.2 is a continuity of the first one, where we investigate the effect of the regularization parameter choice rule and previous iteration processes as well on the final computed numerical result. In case observations taking on the bottom and left surface the computation is given in Section 4.2, while Section 4.2 is devoted to multiple measurements.

To complete this introduction we wish to mention briefly some parameter identification problems in PDEs from boundary observations. The authors Xie and Zou [44], Xu and Zou [46] have used finite element methods to numerically recovered the fluxes on the inaccessible boundary $\Gamma_{i}$ from measurement data of the state on the accessible boundary $\Gamma_{a}$, while the problem of identifying the Robin coefficient on $\Gamma_{i}$ is also investigated by Xu and Zou [45]. Recently, authors of [22, 23] adopted the variational approach of Kohn and Vogelius combined with quadratic stabilizing penalty term and total variation regularization technique to the source term and scalar diffusion coefficient identification, respectively, using observations available on the whole boundary.

Throughout the paper the symbol $A\preceq B$ refers to the inequality $A\leq cB$ for some constant $c$ independent of both $A$ and $B$. In the Lebesgue space $L^{2}(Q)$, where $Q$ is either $\Omega$, $\partial\Omega$ or $\Gamma$, we use for all $y,\widehat{y}\in L^{2}(Q)$ the inner product and the corresponding norm as $(y,\widehat{y})_{Q}:=\int_{Q}y(x)\cdot\widehat{y}(x)dx$ and $\|y\|_{Q}:=(y,y)^{1/2}_{Q}.$ We also use the standard notion of Sobolev spaces $H^{k}(Q):=W^{k}_{2}(Q)$ from, e.g., [2] with notations of its inner product $(\cdot,\cdot)_{k,Q}$, the norm $\|\cdot\|_{k,Q}$ and the semi-norm $|\cdot|_{k,Q}$. Note that $\|\cdot\|_{0,Q}=|\cdot|_{0,Q}=\|\cdot\|_{Q}$.

2 Finite element discretization

2.1 Preliminaries

We remark that the expression

[TABLE]

generates an inner product on the space $H^{1}(\Omega)$ which is equivalent to the usual one, i.e. there exist positive constants $c_{1},c_{2}$ such that

[TABLE]

for all $u\in H^{1}(\Omega)$ and $\beta\in\mathcal{S}_{ad}$, where $c_{1}$ and $c_{2}$ are independent of $\beta$. Therefore, for each $\beta\in\mathcal{S}_{ad}$ the Neumann boundary value problem Eq. 7–Eq. 9 defines a unique weak solution $u=u(\beta):=N_{j_{\delta}}(\beta)$ in the sense that $N_{j_{\delta}}(\beta)\in H^{1}(\Omega)$ and the equation

[TABLE]

is satisfied for all $\phi\in H^{1}(\Omega)$, where $\langle\cdot,\cdot\rangle_{\Omega}$ and $\langle\cdot,\cdot\rangle_{\Gamma}$ stand for the dual pairs \langle\cdot,\cdot\rangle_{\big{(}H^{-1}(\Omega),H^{1}(\Omega)\big{)}} and \langle\cdot,\cdot\rangle_{\big{(}H^{-1/2}(\Gamma),H^{1/2}(\Gamma)\big{)}}, respectively. Furthermore, there holds the estimate

[TABLE]

A function $v=v(\beta):=M_{g_{\delta}}(\beta)$ is said to be a (unique) weak solution of the mixed boundary value problem Eq. 10–Eq. 12 if $M_{g_{\delta}}(\beta)\in H^{1}(\Omega)$ with $M_{g_{\delta}}(\beta)_{|\Gamma}=g_{\delta}$ and the equation

[TABLE]

is satisfied for all $\phi\in H^{1}_{0}(\Omega\cup\Gamma)$, where $H^{1}_{0}(\Omega\cup\Gamma):=\overline{C^{\infty}_{c}(\Omega\cup\Gamma)}^{H^{1}(\Omega)}=\{\phi\in H^{1}(\Omega)~{}|~{}\phi_{|\Gamma}=0\}$, the bar denotes the closure in $H^{1}(\Omega)$ and $C^{\infty}_{c}(\Omega\cup\Gamma)$ is the set of all functions $\phi\in C^{\infty}(\overline{\Omega})$ with $\mbox{supp}\phi$ being a compact subset of $\Omega\cup\Gamma$ (see, e.g., [43, pp. 9, 67]). The above weak solution satisfies the estimate

[TABLE]

Remark 2.1**.**

We mention that under additional assumptions $\boldsymbol{\alpha}\in{W^{1,\infty}(\Omega)}^{d\times d}$, $j_{\delta}\in H^{1/2}(\Gamma)$, $j_{0}\in H^{1/2}(\partial\Omega\setminus\Gamma)$, $g_{\delta}\in H^{3/2}(\Gamma)$, $f\in L^{2}(\Omega)$ and either $\partial\Omega$ is smooth of the class $C^{0,1}$ or the domain $\Omega$ is convex, the weak solutions $N_{j_{\delta}}(\beta),~{}M_{g_{\delta}}(\beta)\in H^{2}(\Omega)$ for all $\beta\in\mathcal{S}_{ad}$ (see, e.g., [15, 43]) satisfying

[TABLE]

and

[TABLE]

We now state some properties of the coefficient-to-solution operators

[TABLE]

Lemma 2.2**.**

Assume that the dimension $d\leq 4$. Then the operators $N_{j_{\delta}}$ and $M_{g_{\delta}}$ are infinitely Fréchet differentiable on the set $\mathcal{S}_{ad}$ with respect to the $L^{2}(\Omega)$-norm. For $\beta\in\mathcal{S}_{ad}$ and $(\kappa_{1},\ldots,\kappa_{m})\in{L^{\infty}(\Omega)}^{m}$ the $m$-th order differentials $D_{N}^{(m)}:=N_{j_{\delta}}^{(m)}(\beta)(\kappa_{1},\ldots,\kappa_{m})\in H^{1}(\Omega)$ and $D_{M}^{(m)}:=M_{g_{\delta}}^{(m)}(\beta)(\kappa_{1},\ldots,\kappa_{m})\in H^{1}_{0}(\Omega\cup\Gamma)$ are the unique solutions to the variational equations

[TABLE]

and

[TABLE]

with $\xi_{i}:=(\kappa_{1},\ldots,\kappa_{i-1},\kappa_{i+1},\ldots,\kappa_{m})\in{L^{\infty}(\Omega)}^{m-1}$, respectively. Furthermore,

[TABLE]

Proof 2.3**.**

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

We mention that the restriction on the dimension $d\leq 4$ in the above Lemma 2.2 is removed if the $L^{\infty}(\Omega)$-norm is taken into account instead of the $L^{2}(\Omega)$-norm (see, e.g., [16, 18]).

Lemma 2.4**.**

*Assume that the sequence $\left(\beta_{n}\right)_{n}\subset\mathcal{S}_{ad}$ converges weakly in $L^{2}(\Omega)$ to an element $\beta$. Then the sequences $\left(N_{j_{\delta}}(\beta_{n})\right)_{n}$ and $\left(M_{g_{\delta}}(\beta_{n})\right)_{n}$ converge respectively to $N_{j_{\delta}}(\beta)$ and $M_{g_{\delta}}(\beta)$ weakly in $H^{1}(\Omega)$ and strongly in the $L^{2}(\Omega)$-norm. *

Proof 2.5**.**

We first note that since $\mathcal{S}_{ad}$ is a convex and closed subset of $L^{2}(\Omega)$, it is weakly closed in $L^{2}(\Omega)$ which implies that $\beta\in\mathcal{S}_{ad}$. Furthermore, it is a weakly∗ compact subset of $L^{\infty}(\Omega)$ (see, e.g., [35, Remark 2.1]). Therefore, by the inequality Eq. 15 and the embeddings $H^{1}(\Omega)\hookrightarrow L^{2}(\Omega)$ as well as $H^{1}(\Omega)\hookrightarrow L^{2}(\partial\Omega)$ being compact (see, e.g., [30, 31]), the sequences (N_{n})_{n}:\equiv\big{(}N_{j_{\delta}}(\beta_{n})\big{)}_{n} and $\left(\beta_{n}\right)_{n}$ have subsequences denoted by the same symbol such that

[TABLE]

as $n\to\infty$, where $\theta_{N}$ is an element of $H^{1}(\Omega)$. We can show that $\theta_{N}=N_{j_{\delta}}(\beta)$. In fact, for all $\phi\in H^{1}(\Omega)$, we have from (14) for all $n\in\mathbb{N}$ that

[TABLE]

By (21), we get $(\boldsymbol{\alpha}\nabla(N_{n}-\theta_{N}),\nabla\phi)_{\Omega}\to 0$ as $n\to\infty$. Furthermore, we have

[TABLE]

as $n\to\infty$, here we used (21) again. Since $\theta_{N}\phi\in L^{1}(\Omega)$, it follows from (20) that $(\beta_{n}-\beta,\theta_{N}\phi)_{\Omega}\to 0$ as $n\to\infty$. Therefore, sending $n\to\infty$ in the above equation we arrive at

[TABLE]

*for all $\phi\in H^{1}(\Omega)$, that means $\theta_{N}=N_{j_{\delta}}(\beta)$. With similar arguments we also obtain that $\left(M_{g_{\delta}}(\beta_{n})\right)_{n}$ converges to $M_{g_{\delta}}(\beta)$ weakly in $H^{1}(\Omega)$, which finishes the proof. *

Together with Eq. 7–Eq. 12, we consider two adjoint problems

[TABLE]

and

[TABLE]

that attain unique weak solutions $A_{N}=A_{N}(\beta)$ and $A_{M}=A_{M}(\beta)$ in the sense that $A_{N}(\beta)\in H^{1}(\Omega)$ and $A_{M}(\beta)\in H^{1}_{0}(\Omega\cup\Gamma)$ satisfy the variational equations

[TABLE]

Furthermore,

[TABLE]

Theorem 2.6**.**

The minimization problem

[TABLE]

attains a minimizer $\beta_{\delta,\rho}$ which satisfies the identity

[TABLE]

*a.e. in $\Omega$, where $A_{N},~{}A_{M}$ come from Eq. 22–Eq. 26. *

Proof 2.7**.**

The existence of a minimizer $\beta$ follows directly from Lemma 2.4. It remains to show the above identity. Due to the first order optimality condition for the minimizer $\beta$, we get for all $\gamma\in\mathcal{S}_{ad}$ that $J^{\prime}_{\delta,\rho}(\beta)(\gamma-\beta)\geq 0$. Setting $\kappa:=\gamma-\beta$, we then obtain

[TABLE]

It follows from Lemma 2.2 that

[TABLE]

which yields

[TABLE]

*for all $\gamma\in\mathcal{S}_{ad}$. This completes the proof. *

2.2 Finite element discretization

Let $\left(\mathcal{T}^{h}\right)_{0<h<1}$ be a family of 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}_{k}$ consists of all polynomial functions of degree less than or equal to $k$. For each $\beta\in\mathcal{S}_{ad}$ the variational equations

[TABLE]

admit unique solutions $u^{h}:={N}^{h}_{j_{\delta}}(\beta)\in\mathcal{V}_{1}^{h}$ and $v^{h}:={M}^{h}_{g_{\delta}}(\beta)\in\mathcal{V}_{1}^{h}$, respectively. Furthermore, the estimates

[TABLE]

hold true.

The solutions $A_{N}=A_{N}(\beta)$ and $A_{M}=A_{M}(\beta)$ of the adjoint problems Eq. 22–Eq. 26 are approximated by $A^{h}_{N}=A^{h}_{N}(\beta)\in\mathcal{V}_{1}^{h}$ and $A^{h}_{M}=A^{h}_{M}(\beta)\in\mathcal{V}_{1,0}^{h}$ satisfying

[TABLE]

The Sobolev number of $W^{m}_{q}(\Omega)$ is defined by $\mbox{sob}(W^{m}_{q}(\Omega)):=m-d/q.$

Lemma 2.8** (Quasi-interpolation operator, see, e.g., [9, 34, 39]).**

There exists an operator $\Pi^{h}:L^{1}(\Omega)\rightarrow\mathcal{V}^{h}_{1}$ such that $\Pi^{h}\varphi^{h}=\varphi^{h}$ for all $\varphi^{h}\in\mathcal{V}^{h}_{1}$ and the limit

[TABLE]

Furthermore, for all $T\in\mathcal{T}^{h}$ we have the local estimate

[TABLE]

*where $0\leq k\leq m\leq 2$ and $1\leq p,q\leq\infty$ such that $\mbox{sob}(W^{m}_{q}(\Omega))>\mbox{sob}(W^{k}_{p}(\Omega))$. If $\phi\in W^{1}_{1}(\Omega)$ has a vanishing trace on a part $\Gamma\subset\partial\Omega$, then so does $\Pi^{h}\phi$. *

With the above notations at hand, the continuous regularized problem \big{(}\mathcal{P}_{\delta,\rho}\big{)} can be discretized by

[TABLE]

Theorem 2.9**.**

The discrete problem \big{(}\mathcal{P}^{h}_{\delta,\rho}\big{)} has a minimizer $\beta^{h}_{\delta,\rho}$ which satisfies for a.e. in $\Omega$ the relation

[TABLE]

Proof 2.10**.**

*The proof follows exactly as in the continuous case, we therefore omit here. *

3 Stability and Convergence

We are in position to prove the stability of the proposed finite element method and the convergence of the regularized finite element approximations to the $\beta^{*}$-minimum-norm solution of the identification problem.

Theorem 3.1** (Stability).**

*Assume that the regularization parameter $\rho$ and the observation data \big{(}j_{\delta},g_{\delta}\big{)}\in H^{-1/2}(\Gamma)\times H^{1/2}(\Gamma) are fixed. For each $n\in\mathbb{N}$ let $\beta_{n}:=\beta^{h_{n}}_{\delta,\rho}$ be an arbitrary minimizer of \big{(}\mathcal{P}^{h_{n}}_{\delta,\rho}\big{)}. Then the sequence $(\beta_{n})_{n}$ has a subsequence converging in the $L^{2}(\Omega)$-norm to an element $\beta_{\delta,\rho}\in\mathcal{S}_{ad}$. Furthermore, $\beta_{\delta,\rho}$ is a minimizer of \big{(}\mathcal{P}_{\delta,\rho}\big{)}. *

Proof 3.2**.**

In view of the proof of Lemma 2.4 we deduce that a subsequence of $(\beta_{n})_{n}$ which is not relabeled and an element $(\beta_{\infty},\theta_{N},\theta_{M})\in\mathcal{S}_{ad}\times H^{1}(\Omega)\times H^{1}(\Omega)$ exist such that

[TABLE]

weakly in $H^{1}(\Omega)$ as $n\to\infty$. We first show that $\theta_{N}=N_{j_{\delta}}(\beta_{\infty})$, i.e. the limit \lim_{n\to\infty}\big{[}N_{j_{\delta}}^{h_{n}}(\beta_{n})-N_{j_{\delta}}(\beta_{\infty}),\phi\big{]}_{\boldsymbol{\alpha},\beta_{\infty},\sigma}=0 holds true for all $\phi\in H^{1}(\Omega)$. In fact, we can rewrite

[TABLE]

It follows from Eq. 29 and Eq. 14 that

[TABLE]

and so

[TABLE]

where we used Eq. 31 and Eq. 35. Furthermore, we have that

[TABLE]

We thus derive from 3.2–3.2 \lim_{n\to\infty}\big{[}N_{j_{\delta}}^{h_{n}}(\beta_{n}),\phi\big{]}_{\boldsymbol{\alpha},\beta_{\infty},\sigma}=\big{[}N_{j_{\delta}}(\beta_{\infty}),\phi\big{]}_{\boldsymbol{\alpha},\beta_{\infty},\sigma}. This also yields

[TABLE]

Likewise, we can show $\theta_{M}=M_{g_{\delta}}(\beta_{\infty})$ and

[TABLE]

Consequently, for all $\beta\in\mathcal{S}_{ad}$ we arrive at

[TABLE]

This means that $\beta_{\infty}$ is a minimizer of \big{(}\mathcal{P}_{\delta,\rho}\big{)}. It remains to show that $(\beta_{n})_{n}$ converges to $\beta_{\infty}$ in the $L^{2}(\Omega)$-norm. For this purpose we take $\beta=\beta_{\infty}$ in the last equation to get

[TABLE]

and then write

[TABLE]

By the equations Eq. 40–3.2, we obtain

[TABLE]

*which finishes the proof. *

To go further, we remark that, due to the assumption on consistency of the system Eq. 1–Eq. 4, the set

[TABLE]

is nonempty, convex, bounded and closed in the $L^{2}(\Omega)$-norm. As a result, there exists a unique solution $\beta^{\dagger}$ of the problem

[TABLE]

which is called by $\beta^{*}$-minimum-norm solution to the identification problem. Let

[TABLE]

Due to the standard theory of the finite element method for elliptic problems (see, e.g., [6]), we get

[TABLE]

in case $N_{j_{\delta}}(\beta),~{}M_{g_{\delta}}(\beta)\in H^{2}(\Omega)$.

Theorem 3.3** (Convergence).**

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

[TABLE]

Furthermore, assume that \big{(}j_{\delta_{n}},g_{\delta_{n}}\big{)}\in H^{-1/2}(\Gamma)\times H^{1/2}(\Gamma) is a sequence satisfying the inequality

[TABLE]

and $\beta_{n}:=\beta^{h_{n}}_{\delta_{n},\rho_{n}}$ denotes an arbitrary minimizer of \big{(}\mathcal{P}^{h_{n}}_{\delta_{n},\rho_{n}}\big{)} for each $n\in\mathbb{N}$. Then:

(i) The whole sequence $(\beta_{n})_{n}$ converges in the $L^{2}(\Omega)$-norm to $\beta^{\dagger}$.

*(ii) The corresponding state sequences \big{(}N^{h_{n}}_{j_{\delta_{n}}}(\beta_{n})\big{)}_{n} and \big{(}M^{h_{n}}_{g_{\delta_{n}}}(\beta_{n})\big{)}_{n} converge in the $H^{1}(\Omega)$-norm to the exact state $\Phi^{\dagger}=\Phi(j^{\dagger},g^{\dagger},\beta^{\dagger})$ of the problem Eq. 1–Eq. 4. *

Note that in case the exact solution $\Phi^{\dagger}\in H^{2}(\Omega)$ (cf. Remark 2.1) then the convergences (i) and (ii) are obtained if the regularization parameter is chosen such that $\rho_{n}\to 0$, $\delta_{n}/\sqrt{\rho_{n}}\to 0$ and $h^{2}_{n}/\sqrt{\rho_{n}}\to 0$ as $n\to\infty$.

Proof 3.4** (Proof of Theorem 3.3).**

We have from the optimality of $\beta_{n}$ for each $n$ that

[TABLE]

Note that $N_{j^{\dagger}}(\beta^{\dagger})=M_{g^{\dagger}}(\beta^{\dagger})$, we thus have

[TABLE]

By the identities Eq. 29–Eq. 30 and the inequality Eq. 45, we get

[TABLE]

which yields

[TABLE]

It follows from 3.4–Eq. 48 that

[TABLE]

and

[TABLE]

Next, we get that

[TABLE]

as $n\to\infty$, by Eq. 49. Furthermore, since $\beta_{n}\in\mathcal{S}_{ad}$ for all $n\in\mathbb{N}$, in view of the proof of Theorem 3.1, a subsequence of $\left(\beta_{n}\right)_{n}$ not relabeled and an element $\widehat{\beta}\in\mathcal{S}_{ad}$ exist such that

[TABLE]

which also implies that

[TABLE]

and so that $\widehat{\beta}\in\Pi_{\mathcal{S}_{ad}}(j^{\dagger},g^{\dagger})$. Then, combining Eq. 51 with Eq. 50 and using the uniqueness of the $\beta^{*}$-minimum-norm solution, we obtain

[TABLE]

the assertion (i) is thus proved. For (ii) we have for all $\phi^{h_{n}}\in\mathcal{V}^{h_{n}}_{1}$ that

[TABLE]

which together with the inequality Eq. 13 and the continuous embedding $H^{1}(\Omega)\hookrightarrow L^{4}(\Omega)$ imply

[TABLE]

Therefore, we arrive at

[TABLE]

with the aid of the limit Eq. 52. Consequently, we obtain from 3.4 and Eq. 53 that

[TABLE]

*as $n\to\infty$. By the similar arguments, we also get \lim_{n\to\infty}\big{\|}M^{h_{n}}_{g_{\delta_{n}}}(\beta_{n})-M_{g^{\dagger}}(\beta^{\dagger})\big{\|}_{1,\Omega}=0, which finishes the proof. *

4 Gradient projection algorithm and numerical implementation

4.1 Algorithm

In this section we present the gradient projection algorithm with Armijo steplength rule (cf. [27, 38]) for numerical solution of the minimization problem \big{(}\mathcal{P}^{h}_{\delta,\rho}\big{)}.

In view of the proof of Theorem 2.6, we first note that for each $\beta\in\mathcal{S}_{ad}$ the $L^{2}$-gradient of the cost functional $J^{h}_{\delta,\rho}$ of the problem \big{(}\mathcal{P}^{h}_{\delta,\rho}\big{)} at $\beta$ is given by

[TABLE]

The algorithm is then read as: given a step size control $\mu\in(0,1)$, an initial approximation $\beta_{0}$, number of iteration $N$ and setting $k=0$.

Compute $N^{h}_{j_{\delta}}(\beta_{k})$ and $M^{h}_{g_{\delta}}(\beta_{k})$ from the equations

[TABLE]

and then the solutions $A^{h}_{N}(\beta_{k})$ and $A^{h}_{M}(\beta_{k})$ of the adjoint problems

[TABLE] 2. 2.

Compute the corresponding value of the cost functional

[TABLE]

as well as the gradient

[TABLE] 3. 3.

Compute

[TABLE]

and then the corresponding states $N^{h}_{j_{\delta}}(\widehat{\beta}_{k})$ and $M^{h}_{g_{\delta}}(\widehat{\beta}_{k})$, the value of the cost functional $J^{h}_{\delta,\rho}(\widehat{\beta}_{k})$ as well. 4. 4.

Compute the quantity

[TABLE]

for a small positive constant $\tau=10^{-4}$.

1. (a)

If $Q\geq 0$

go to the next step (b) below

else

set $\mu:=\frac{\mu}{2}$ and then go back to the step 3. 2. (b)

Update $\beta_{k}=\widehat{\beta}_{k}$, set $k=k+1$. 5. 5.

Compute

[TABLE]

for $\tau_{1}:=10^{-3}h$ and $\tau_{2}:=10^{-2}h$. If $\mbox{Tolerance}\leq 0$ or $k>N$, then stop; otherwise go back to the step 1.

4.2 Numerical implementation

For illustrating the theoretical result we consider the system Eq. 1–Eq. 4 with $\Omega=\{x=(x_{1},x_{2})\in\mathbb{R}^{2}~{}|~{}-1<x_{1},x_{2}<1\}$. For discretization we divide the interval $(-1,1)$ into $\ell$ equal segments, and so 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}$.

The source function $f$ is assumed to be discontinuous and defined as

[TABLE]

where $\chi_{D}$ is the characteristic function of the Lebesgue measurable set

[TABLE]

We assume that entries of the symmetric diffusion matrix $\boldsymbol{\alpha}$ are discontinuous which are defined as

[TABLE]

with

[TABLE]

Furthermore, the special function $\sigma$ is chosen to be zero while the constants appearing in the admissible set $\mathcal{S}_{ad}$ defined by (5) are chosen as $\underline{\beta}=0.05$ and $\overline{\beta}=10$.

The sought reaction coefficient $\beta^{\dagger}$ is assumed to be discontinuous and given by

[TABLE]

with

[TABLE]

The Neumann boundary condition on the bottom and left surface is given by

[TABLE]

and on the right and top surface

[TABLE]

with the constants $A,B,C$ and $D$ discussed in details later. The exact state $\Phi^{\dagger}$ is then computed from the finite element equation $K\Phi^{\dagger}=F$, where $K$ and $F$ are the stiffness matrix and the load vector associated with the problem Eq. 1–Eq. 3, respectively. The Dirichlet boundary condition $g^{\dagger}$ in Eq. 4 is then defined as $g^{\dagger}=\gamma_{\Gamma}\Phi^{\dagger}$, the Dirichlet trace of $\Phi^{\dagger}$ on the boundary $\Gamma$.

We use the algorithm which is described in Section 4.1 for computing the numerical solution of the problem \big{(}\mathcal{P}_{\delta,\rho}^{h}\big{)}. The step size control is chosen with $\mu=0.75$. The initial approximation and an a priori estimate are the constant functions defined by $\beta_{0}=\beta^{*}=1.5$. Our computational process will be started with the coarsest level $\ell=4$. In each iteration $k$ we compute Tolerance defined by Eq. 54. Then the iteration is stopped if $\mbox{Tolerance}\leq 0$ or the number of iterations reaches the maximum iteration counted 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 the initial approximation and an a priori estimate as well for the algorithm on this finer mesh, and so on for $\ell=16,32,64$.

We mention that in our numerical implementation the sought reaction coefficient is chosen to be discontinuous. To reconstruct such a discontinuous function one usually employs the total variation regularization. We will discuss the details in the last section §5.

{exam}

In this example we assume $(A,B,C,D)=(1,-2,3,-4)$ while observations are taken on the bottom surface $\Gamma_{\mbox{observation}}:=\Gamma_{\mbox{bottom}}:=[-1,1]\times\{-1\}$ only. We assume that noisy observations are available in the form

[TABLE]

for some $\theta_{\ell}>0$ depending on $\ell$, 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 set of boundary nodes of the triangulation $\mathcal{T}^{h_{\ell}}$ which belong to $\Gamma_{\mbox{observation}}$. The measurement error is then computed as \delta_{\ell}=\big{\|}j_{\delta_{\ell}}-j^{\dagger}\big{\|}_{L^{2}(\Gamma)}+\big{\|}g_{\delta_{\ell}}-g^{\dagger}\big{\|}_{L^{2}(\Gamma)}. To satisfy the condition Eq. 44 in Theorem 3.3 we below take $\theta_{\ell}=h_{\ell}\sqrt{10\cdot\rho_{\ell}}$ and the regularization parameter $\rho=\rho_{\ell}=0.001\sqrt{h_{\ell}}$.

Let $\beta_{\ell}$ denote the reaction coefficient obtained at the final iteration of the algorithm corresponding to the refinement level $\ell$. We then use the following abbreviations for the errors

[TABLE]

where

[TABLE]

$D^{h_{\ell}}_{g^{\dagger}_{\delta_{\ell}}}({\beta_{\ell}})$ and $D^{h_{\ell}}_{\widehat{g}^{\dagger}}(\beta^{\dagger})$ are numerical solutions of the problem (1) with $\beta=\beta_{\ell}$ and $\beta=\beta^{\dagger}$, respectively supplemented with the Dirichlet boundary condition $\Phi_{|\partial\Omega}=g^{\dagger}_{\delta_{\ell}}$ and $\Phi_{|\partial\Omega}=\widehat{g}^{\dagger}$.

The numerical result is summarized in Table 1 and Table 2, where we present the refinement level $\ell$, the mesh size $h_{\ell}$ of the triangulation, the regularization parameter $\rho_{\ell}$, the measurement noise $\delta_{\ell}$, and the errors $L^{2}_{\beta}$, $L^{2}_{N}$, $L^{2}_{M}$, $L^{2}_{D}$.

The experimental order of convergence (EOC) is presented in Table 2, where

[TABLE]

and $\Theta(h)$ is an error function with respect to the mesh size $h$.

The error history given in Table 1 shows that the algorithm performs well for our identification problem. We also observe a decrease of all errors as the noise level and the mesh size gets smaller as expected from our convergence result in Theorem 3.3.

All figures presented hereafter correspond to the finest level $\ell=64$. Figure 1 from left to right shows the numerical solution $\beta_{\ell}$ computed by the algorithm at the final 523${}^{\mbox{\tiny th}}$-iteration and the differences $\beta_{\ell}-I^{h_{\ell}}_{1}\beta^{\dagger}$, $N^{h_{\ell}}_{j^{\dagger}}(\beta^{\dagger})-N^{h_{\ell}}_{j_{\delta_{\ell}}}({\beta_{\ell}})$, $D^{h_{\ell}}_{{\widehat{g}^{\dagger}}}(\beta^{\dagger})-D^{h_{\ell}}_{g^{\dagger}_{\delta_{\ell}}}({\beta_{\ell}})$, where $I^{h_{\ell}}_{1}$ is the usual Lagrange node value interpolation operator.

{exam}

In this example we consider regularization parameter in the forms

[TABLE]

while $\theta_{\ell}=10^{-2}h_{\ell}\sqrt{\rho_{\ell}}$. Using the computational process which was described as in Section 4.2 starting with $\ell=4$, in Table 3 we perform the numerical result for the finest grid $\ell=64$ and with different values of $\rho_{\ell}$ above.

We note that the numerical result is not so good in first two cases, where $\rho_{\ell}=h_{\ell}^{2}$ and $\rho_{\ell}=0.1h_{\ell}^{2}$. This indicates that previous iteration processes corresponding to the coarse grid levels $\ell=4,8,...$ strongly affect the final obtained numerical result with respect to the finest grid level $\ell=64$.

A crucial problem of Tikhonov regularization as well as other regularization methods is the choice of the regularization parameter. If the regularization parameter is too large, one obtains only a poor approximation of the exact solution even for exact data (i.e. the error level of observations is zero), while if it is too small, the reconstruction becomes slow or unstable. The discussion on this subject is still ongoing and for a survey the reader may consult [12, 28, 32]. We mention that by the relation (44), if $\delta_{\ell}\rho_{\ell}^{-1/2}\to 0$ the regularized approximation $\beta_{\ell}$ is convergent in the $L^{2}(\Omega)$-norm to the sought coefficient $\beta^{\dagger}$ if $\rho_{\ell}$ is chosen to be $h^{2}_{\ell}$. However, the computation given in Table 3 shows in case $\rho_{\ell}=h^{2}_{\ell}$ and $\rho_{\ell}=0.1h^{2}_{\ell}$ we get only poor approximations of the sought coefficient, though the data noise levels $\delta_{\ell}$ are quite small.

With $\rho_{\ell}=h_{\ell}^{2}$, Figure 2 from left to right shows the computed numerical solution $\beta_{\ell}$ at the final iteration, the differences $N^{h_{\ell}}_{j^{\dagger}}(\beta^{\dagger})-N^{h_{\ell}}_{j_{\delta_{\ell}}}({\beta_{\ell}})$ and $D^{h_{\ell}}_{\widehat{g}^{\dagger}}(\beta^{\dagger})-D^{h_{\ell}}_{g^{\dagger}_{\delta_{\ell}}}({\beta_{\ell}})$, while the differences $\beta_{\ell}-I^{h_{\ell}}_{1}\beta^{\dagger}$ and $\beta_{\ell}(x_{1},0)-I^{h_{\ell}}_{1}\beta^{\dagger}(x_{1},0)$ as the second variable $x_{2}=0$ are shown respectively in left two figures of Figure 3. The right one of Figure 3 performs the graphs $\beta_{\ell}(x_{1},0)$ (blue one) and $I^{h_{\ell}}_{1}\beta^{\dagger}(x_{1},0)$.

Figure 4 and Figure 5 perform analogous numerical computations for $\rho=0.1h_{\ell}^{2}$.

{exam}

We now consider the case $\Gamma_{\mbox{observation}}:=\Gamma_{\mbox{bottom-left}}$ includes the bottom surface and the left surface of the domain $\Omega$, i.e. $\Gamma_{\mbox{bottom-left}}=\{x=(x_{1},x_{2})\in\overline{\Omega}~{}|~{}x_{2}=-1\}\cup\{x=(x_{1},x_{2})\in\overline{\Omega}~{}|~{}x_{1}=-1\}$. In this case $\theta_{\ell}$ in Eq. 56 is changed to $\theta_{\ell}=\frac{1}{2}h_{\ell}\sqrt{10\cdot\rho_{\ell}}$, while the regularization parameter $\rho=\rho_{\ell}=0.001\sqrt{h_{\ell}}$ as in Section 4.2. The computational result shows in Table 4 and Table 5 below.

Compared with Table 1, we do not see the difference clearly in the obtained numerical result between $\Gamma_{\mbox{observation}}:=\Gamma_{\mbox{bottom}}$ and $\Gamma_{\mbox{observation}}:=\Gamma_{\mbox{bottom-left}}$. Table 2 and Table 5 show that the experimental order of convergence for $L^{2}$-error of the sought coefficient is of the first order. To the best of our knowledge, there have been no papers providing an error bound for the identification problem mentioned in the present work so far. Nevertheless, we would like to distinguish that utilizing distributed observations for the diffusion coefficient identification problem, under the suitable source conditions authors of [14, 20, 29] also proved the first order convergence rate $\mathcal{O}(h)$ of the piecewise linear finite element approximations to the identified solution in the $L^{2}$-norm.

Figure 6 below performs the graphs of the computation in this case.

{exam}

In the last example we assume that $I$ multiple measurements on the bottom surface and the left surface $\left(j_{\delta}^{i},g_{\delta}^{i}\right)_{i=1,\ldots,I}$ are available. Then, the cost functional $J^{h}_{\delta,\rho}$ and the problem \big{(}\mathcal{P}^{h}_{\rho,\delta}\big{)} can be rewritten as

[TABLE]

which also attains a solution $\bar{\beta}^{h}_{\delta,\rho}$. Let $j^{\dagger}_{(A,B,C,D)}$ be given by Eq. 55 which depends on the constants $A,B,C,D$ and $g^{\dagger}_{(A,B,C,D)}:=\gamma_{\Gamma}N_{j^{\dagger}_{(A,B,C,D)}}(\beta^{\dagger})$. The noisy observations are assumed to give by

[TABLE]

where $\theta>0$ is independent of the mesh size, regularization parameter and noise level. $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)$, and $\partial M^{h_{\ell}}$ is the set of boundary nodes of the triangulation $\mathcal{T}^{h_{\ell}}$ which belong to $\Gamma_{\mbox{observation}}:=\Gamma_{\mbox{bottom-left}}$.

The numerical result in Table 6 presents for $\theta=0.2$ and with respect to

$\bullet$ $I=1$ measurement: $(A,B,C,D)=(1,-2,3,-4)$,

$\bullet$ $I=6$ measurements: fixing $D=-4$ and taking $(A,B,C)$ equals to all permutations of the set $\{1,-2,3\}$,

$\bullet$ $I=16$ measurements: taking $(A,B,C,D)$ equals to all permutations of the set $\{1,-2,3,-4\}$.

We observe that the use of multiple measurements improves the solution to yield an acceptable result even in the presence of relatively large noise as Table 6 below.

Finally, in Figure 7–Figure 9 we perform the graphs of the computation for the multiple measurements $I=1,6,16$, respectively, which include the computed numerical solution $\beta_{\ell}$, and the differences $\beta_{\ell}-I^{h_{\ell}}_{1}\beta^{\dagger}$, $N^{h_{\ell}}_{j^{\dagger}}(\beta^{\dagger})-N^{h_{\ell}}_{j_{\delta_{\ell}}}({\beta_{\ell}})$, $D^{h_{\ell}}_{\widehat{g}^{\dagger}}(\beta^{\dagger})-D^{h_{\ell}}_{g^{\dagger}_{\delta_{\ell}}}({\beta_{\ell}})$ for $\ell=64$, $\rho=\rho_{\ell}=0.001\sqrt{h_{\ell}}$ and $\delta_{\ell}=0.4773$.

5 Conclusions

The paper is devoted to the problem of identifying the reaction coefficient $\beta=\beta(x)$ in the equation -\nabla\cdot\big{(}\boldsymbol{\alpha}(x)\nabla\Phi(x)\big{)}+\beta(x)\Phi(x)=f(x),~{}x\in\Omega supplemented with the boundary conditions $\boldsymbol{\alpha}(x)\nabla\Phi(x)\cdot\vec{\boldsymbol{n}}(x)+\sigma(x)\Phi(x)=j^{\dagger}(x)$,  $x\in\Gamma$,  $\boldsymbol{\alpha}(x)\nabla\Phi(x)\cdot\vec{\boldsymbol{n}}(x)+\sigma(x)\Phi(x)=j_{0}(x),~{}x\in\partial\Omega\setminus\Gamma$ and $\Phi(x)=g^{\dagger}(x),~{}x\in\Gamma$ from the measurement data \big{(}j_{\delta},g_{\delta}\big{)} of the exact \big{(}j^{\dagger},g^{\dagger}\big{)}. In this context, the special functions $\boldsymbol{\alpha},f$ and $\sigma$ are given, where $\vec{\boldsymbol{n}}$ is the unit outward normal on $\partial\Omega$.

With the available measurement data \big{(}j_{\delta},g_{\delta}\big{)} at hand, we for each $\beta\in\mathcal{S}_{ad}$ consider simultaneously two problems

[TABLE]

and denote respectively by $N_{j_{\delta}}(\beta)$ and $M_{g_{\delta}}(\beta)$ their unique weak solutions. A minimizer $\beta_{\delta,\rho}$ of the Tikhonov regularized minimization problem

[TABLE]

is considered as reconstruction, where $\rho>0$ is the regularization parameter and $\beta^{*}$ is an a priori estimate of the true coefficient, with the regularization term $R(\beta,\beta^{*}):=\|\beta-\beta^{*}\|^{2}_{L^{2}(\Omega)}$. Let $N^{h}_{j_{\delta}}(\beta)$ and $M^{h}_{g_{\delta}}(\beta)$ be corresponding approximations of $N_{j_{\delta}}(\beta)$ and $M_{g_{\delta}}(\beta)$ in the finite dimensional space $\mathcal{V}^{h}_{1}$ of piecewise linear, continuous finite elements. We then examine the discrete regularized problem

[TABLE]

which also attains a minimizer $\beta^{h}_{\delta,\rho}$.

We show that when $\delta$ and $\rho$ are fixed the sequence of minimizers \big{(}\beta^{h}_{\delta,\rho}\big{)}_{h>0} to \big{(}\mathcal{P}^{h}_{\delta,\rho}\big{)} can be extracted a subsequence which converges in the $L^{2}(\Omega)$-norm to a solution $\beta_{\delta,\rho}$ of \big{(}\mathcal{P}_{\delta,\rho}\big{)} as the mesh size $h\to 0$. Furthermore as $h,\delta\to 0$ and with an appropriate a priori regularization parameter choice $\rho=\rho(h,\delta)\to 0$, the whole sequence \big{(}\beta^{h}_{\delta,\rho}\big{)}_{\rho>0} converges in the $L^{2}(\Omega)$-norm to the $\beta^{*}$-minimum-$R$ solution $\beta^{\dagger}$ of the identification problem. The corresponding state sequences \big{(}N^{h}_{j_{\delta}}\big{(}\beta^{h}_{\delta,\rho}\big{)}\big{)}_{\rho>0} and \big{(}M^{h}_{g_{\delta}}\big{(}\beta^{h}_{\delta,\rho}\big{)}\big{)}_{\rho>0} then converge in the $H^{1}(\Omega)$-norm to the exact state $\Phi^{\dagger}=\Phi(j^{\dagger},g^{\dagger},\beta^{\dagger})$ of the problem.

For the particular interest in estimating probably discontinuous reaction coefficients one can employ the total variation regularization $R(\beta)=\int_{\Omega}|\nabla\beta|$ which was originally introduced in image denoising by authors of [37] and was also applied to several ill-posed and inverse problems, see, e.g., [1, 4, 7, 8].

Alternatively, starting with [11], the sparsity regularization has been applied to the parameter identification problem for PDEs, see the review paper [24]. In this situation the penalty term is defined by $R(\beta)=\sum_{i}w_{i}|\langle\beta,\varphi_{i}\rangle|^{p}$ for $1\leq p\leq 2$ and with weights $w_{i}\geq\underline{w}>0$. Here $\{\varphi_{i}\}_{i}$ is an orthonormal basis or overcomplete frame of the space $L^{2}(\Omega)$. In practice the system $\{\varphi_{i}\}_{i}$ is chosen to be highly dependent on the structure of the reconstructed solution, that could be a Fourier representation for oscillatory features, wavelets for pointwise singularities [40]. Together with the total variation penalty term, adopting the sparsity regularization method for the identification problem may be a work for us in future.

Acknowledgments

The author would like to thank the Editor and Referees for their valuable comments and suggestions.