Variational method for multiple parameter identification in elliptic PDEs

TL;DR

This paper introduces a variational approach with Tikhonov regularization for simultaneously identifying multiple parameters in elliptic PDEs from limited measurement data, ensuring convergence and error analysis.

Contribution

It presents a novel variational method combined with finite element discretization for multi-parameter inverse problems in elliptic PDEs, with proven convergence and error bounds.

Findings

Convergence of the proposed method is established.

Error bounds for the parameter estimates are derived.

The approach effectively handles weaker measurement data.

Abstract

In the present paper we investigate the inverse problem of identifying simultaneously the diffusion matrix, source term and boundary condition as well as the state in the Neumann boundary value problem for an elliptic partial differential equation (PDE) from a measurement data, which is weaker than required of the exact state. A variational method based on energy functions with Tikhonov regularization is here proposed to treat the identification problem. We discretize the PDE with the finite element method and prove the convergence as well as analyse error bounds of this approach.

Tables

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

$\ell$	$h_{\ell }$	${\rho }_{\ell }$	${\delta }_{\ell }$	Iterate
3	0.9428	9.4281e-4	0.1755	800
6	0.4714	4.7140e-4	0.3847	800
12	0.2357	2.3570e-4	0.3334	800
24	0.1179	1.1790e-4	0.1508	800
48	5.8926e-2	5.8926e-5	6.5163e-2	800
96	2.9463e-2	2.9463e-5	2.9896e-2	800

Table 2. Table 2: Errors Δ Δ \Delta , Σ Σ \Sigma and Λ Λ \Lambda and Experimental order of convergence between finest and coarsest level.

$\Delta$	$\Sigma$	$\Lambda$	EOCΔ	EOCΣ	EOCΛ
0.6349	6.2551e-2	0.2789	—	—	—
0.1974	3.7602e-2	0.1847	1.6854	0.7342	0.5946
8.3571e-2	1.7066e-2	0.1382	1.2400	1.1397	0.4184
3.1600e-2	5.4913e-3	6.1769e-2	1.4031	1.6359	1.1618
1.1524e-2	9.4491e-4	2.0742e-2	1.4553	2.5389	1.5743
4.1183e-3	2.2575e-4	8.9372e-3	1.4845	2.0655	1.2147
Mean	of	EOC	1.4537	1.6228	0.9928

Equations

$$-\nabla \cdot (Q\nabla \Phi )$$

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

$${\mathbf{L}}_{\text{mbox}sym}^p(\Omega )$$

$${\mathcal{H}}_{ad}$$

$${\mathcal{Q}}_{ad}$$

$${\mathcal{F}}_{ad}$$

$${\mathcal{G}}_{ad}$$

$$\gamma :H^1(\Omega )\to H^{1/2}(\partial \Omega )$$

$$H_{\diamond }^1(\Omega )$$

$$C_{\Omega }{\int }_{\Omega }{\phi }^{2}dx$$

$$\Vert \phi {\Vert }_{H^1(\Omega )}^2$$

$${\int }_{\Omega }Q\nabla \Phi \cdot \nabla \phi dx$$

$${\Vert \Phi \Vert }_{H^1(\Omega )}$$

$$C_{\mathcal{N}}$$

$$\mathcal{U}:{\mathcal{H}}_{ad}$$

$$\left(Q^{†},f^{†},g^{†}\right)$$

$$\mathcal{R}(Q,f,g)$$

$${\Vert {\Phi }^{†}-z_{\delta }\Vert }_{L^2(\Omega )}$$

$$\underset{(Q,f,g)\in {\mathcal{H}}_{ad}}{\min }{\Vert {\mathcal{U}}_{Q,f,g}^h-{\Pi }^h{z}_{\delta }\Vert }_{L^2(\Omega )}^2+\rho \mathcal{R}(Q,f,g)$$

$$(Q,f,g)\in {\mathcal{H}}_{ad}$$

$$\underset{(Q,f,g)\in {\mathcal{H}}_{ad}}{\min }{\mathcal{J}}_{\delta }^h(Q,f,g)+\rho \mathcal{R}(Q,f,g)$$

$$\mathcal{D}$$

$$\Vert (H,l,s){\Vert }_{{\mathbf{L}}_{\text{mbox}sym}^2(\Omega )\times L^2(\Omega )\times L^2(\partial \Omega )}$$

$$\Vert (H,l,s){\Vert }_{{\mathbf{L}}_{\text{mbox}sym}^{\infty }(\Omega )\times L^2(\Omega )\times L^2(\partial \Omega )}$$

$$\mathcal{U}:{\mathcal{H}}_{ad}\subset {\mathbf{L}}_{\text{mbox}sym}^{\infty }(\Omega )\times L^2(\Omega )\times L^2(\partial \Omega )$$

$$\Gamma :=(Q,f,g)\in {\mathcal{H}}_{ad}$$

$${\int }_{\Omega }Q\nabla {\xi }_{\lambda }\cdot \nabla \phi$$

$$\mathcal{K}$$

$$(A-P_{\mathcal{K}}(A))\cdot (B-P_{\mathcal{K}}(A))$$

$$(\xi \otimes \eta {)}_{1\le i,j\le d}\in {\mathcal{S}}_d\text{mbox}with(\xi \otimes \eta {)}_{ij}:=\frac{1}{2}({\xi }_i{\eta }_j+{\xi }_j{\eta }_i)\text{mbox}foralli,j=1,\cdots ,d.$$

$$\mathbf{D}$$

$$\underset{m\to \infty }{\lim }{\int }_{\Omega }{q}_{n_m}{\theta }_1$$

Full text

Variational method for multiple parameter identification in elliptic PDEs

Tran Nhan Tam Quyen††Email: quyen.tran@uni-hamburg.de,  quyen.tran@uni-goettingen.de††The author gratefully acknowledges support of the Alexander von Humboldt-Foundation and the Lothar Collatz Center for Computing in Science at the University of Hamburg, and the University of Goettingen, State of Lower Saxony, Germany

Department of Mathematics, University of Hamburg, Bundesstraße 55, D-20146 Hamburg, Germany

Institut für Numerische und Angewandte Mathematik, Universität Göttingen, Lotzestraße 16-18, D-37083 Göttingen, Germany

Abstract: In the present paper we investigate the inverse problem of identifying simultaneously the diffusion matrix, source term and boundary condition in the Neumann boundary value problem for an elliptic partial differential equation (PDE) from a measurement data, which is weaker than required of the exact state. A variational method based on energy functions with Tikhonov regularization is here proposed to treat the identification problem. We discretize the PDE with the finite element method and prove the convergence as well as analyse error bounds of this approach.

Key words and phrases: Multiple parameter identification, diffusion matrix, source term, boundary condition, ill-posed problem, finite element method.

AMS Subject Classifications: 65N21, 65N12, 35J25, 35R30

1 Introduction

Let $\Omega$ be an open bounded connected domain of $R^{d},1\leq d\leq 3$ with polygonal boundary $\partial\Omega$. In this paper we study the problem of identifying simultaneously the diffusion matrix $Q$, source term $f$ and boundary condition $g$ as well as the state $\Phi$ in the Neumann boundary value problem for the elliptic PDE

[TABLE]

from a measurement $z_{\delta}\in L^{2}(\Omega)$ of the solution $\Phi\in H^{1}(\Omega)$, where $\vec{n}$ is the unit outward normal on $\partial\Omega$.

To formulate precisely our problem, let us first denote by $\mathcal{S}_{d}$ the set of all symmetric, real $d\times d$-matrices equipped with the inner product $M\cdot N:=\mbox{trace}(MN)$ and the corresponding norm $\|M\|_{\mathcal{S}_{d}}=(M\cdot M)^{1/2}=\left(\sum_{i,j=1}^{d}m_{ij}^{2}\right)^{1/2}$, where $M:=(m_{ij})_{i,j=\overline{1,d}}$. Furthermore, for $1\leq p\leq\infty$ we set

[TABLE]

In $\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)$ we use the scalar product $\left(H^{1},H^{2}\right)_{\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)}=\sum_{i,j=1}^{d}(h^{1}_{ij},h^{2}_{ij})_{L^{2}(\Omega)}$ and the corresponding norm $\|H\|_{\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)}:=\left(\sum_{i,j=1}^{d}\|h_{ij}\|^{2}_{L^{2}(\Omega)}\right)^{1/2}=\left(\int_{\Omega}\|H(x)\|^{2}_{\mathcal{S}_{d}}dx\right)^{1/2}$, while the space $\mathbf{L}^{\infty}_{\mbox{\tiny sym}}(\Omega)$ is endowed with the norm $\|H\|_{\mathbf{L}^{\infty}_{\mbox{\tiny sym}}(\Omega)}:=\max_{i,j=\overline{1,d}}\|h_{ij}\|_{L^{\infty}(\Omega)}$.

Let us denote by

[TABLE]

with

[TABLE]

and $\underline{q}$, $\overline{q}$ being given constants satisfying $\overline{q}\geq\underline{q}>0$. Let

[TABLE]

be the continuous Dirichlet trace operator and $H^{1}_{\diamond}(\Omega)$ be the closed subspace of $H^{1}(\Omega)$ consisting all functions with zero-mean on the boundary, i.e.

[TABLE]

while $C_{\Omega}$ stands for the positive constant appearing in the Poincaré-Friedrichs inequality (cf.  [38])

[TABLE]

Then, due to the coervicity condition

[TABLE]

holding for all $\varphi\in H^{1}_{\diamond}(\Omega),Q\in\mathcal{Q}_{ad}$ and the Lax-Milgram lemma, we conclude for each $(Q,f,g)\in\mathcal{H}_{ad}$, there exists a unique weak solution $\Phi$ of (1.1)–(1.2) in the sense that $\Phi\in H^{1}_{\diamond}(\Omega)$ and satisfies the identity

[TABLE]

for all $\varphi\in H^{1}_{\diamond}(\Omega)$. Here the expressions $(\cdot,\cdot)$ and $\langle\cdot,\cdot\rangle$ stand for the scalar product on space $L^{2}(\Omega)$ and $L^{2}(\partial\Omega)$, respectively. Furthermore, there holds the priori estimate

[TABLE]

with

[TABLE]

Then we can define the non-linear coefficient-to-solution operator

[TABLE]

which maps each $(Q,f,g)\in\mathcal{H}_{ad}$ to the unique weak solution $\mathcal{U}_{Q,f,g}:=\Phi$ of the problem (1.1)–(1.2). Here, for convenience in computing numerical solutions of the pure Neumann problem we normalize the solution with vanishing mean on the boundary (cf., e.g., [23, Subsection 5.2], [28, Section 2]); however, all results performed in the present paper are still valid for the normalization of solutions of the Neumann problem with zero-mean over the domain, i.e. \mathcal{U}_{Q,f,g}\in\big{\{}u\in H^{1}(\Omega)~{}\Big{|}~{}\int_{\Omega}udx=0\big{\}}. The identification problem is now stated as follows:

*Given $\Phi^{\dagger}:=\mathcal{U}_{Q,f,g}\in H^{1}_{\diamond}(\Omega)$, find an element $(Q,f,g)\in\mathcal{H}_{ad}$

such that (1.6) is satisfied with $\Phi^{\dagger}$ and $Q,f,g$.*

This inverse problem may have more than one solution and it is highly ill-posed. In fact, assume that the exact $\Phi^{\dagger}\in C^{2}_{c}(\Omega)$, the space of all functions having second-order derivatives with compact support in $\Omega$. Then, for all $Q\in{C^{1}(\Omega)}^{d\times d}\cap\mathcal{Q}_{ad}$ the element $(\bar{Q},\bar{f},\bar{g}):=\left(Q,-\nabla(Q\cdot\nabla\Phi^{\dagger}),0\right)$ is a solution of the above identification problem, i.e. $\mathcal{U}_{\bar{Q},\bar{f},\bar{g}}=\Phi^{\dagger}$. In other words we are considering to solve an equation $\mathcal{U}_{Q,f,g}=\Phi^{\dagger}$, where the forward operator $\mathcal{U}$ is non-linear and non-injective. Without using additional objective a-priori information or without exploiting other observation data as considering here, it is difficult for us to classify sought targets. Following the general convergence theory for ill-posed problems (see, e.g., [9, Chapter 5] and [43, Subsection 3.2.1], or the classical monograph [15, Section 10.1]), in the present paper we are interested in finding exact solutions with penalty minimizing, which is defined as

[TABLE]

where $\mathcal{I}(\Phi^{\dagger}):=\left\{(Q,f,g)\in\mathcal{H}_{ad}~{}|~{}\mathcal{U}_{Q,f,g}=\Phi^{\dagger}\right\}$ and the penalty term

[TABLE]

We note that the admissible set $\mathcal{I}(\Phi^{\dagger})$ of the problem (1.8) is non-empty, convex and weakly closed in $\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$, so that the minimizer $(Q^{\dagger},f^{\dagger},g^{\dagger})$ is defined uniquely. Furthermore, the exact data $\Phi^{\dagger}$ may not be known in practice, thus we assume instead of $\Phi^{\dagger}$ to have a measurement $z_{\delta}\in L^{2}(\Omega)$ such that

[TABLE]

Our identification problem is now to reconstruct $\left(Q^{\dagger},f^{\dagger},g^{\dagger}\right)\in\mathcal{H}_{ad}$ from $z_{\delta}$.

Let $\left(\mathcal{T}^{h}\right)_{0<h<1}$ denote a family of triangulations of the domain $\overline{\Omega}$ with the mesh size $h$ and $\mathcal{U}^{h}$ be the approximation of the operator $\mathcal{U}$ on the piecewise linear, continuous finite element space associated with $\mathcal{T}^{h}$. Furthermore, let $\Pi^{h}$ be the Clément’s mollification interpolation operator (cf. §2). The standard method for solving the above mentioned identification problem is the output least squares one with Tikhonov regularization, i.e. one considers a minimizer of the problem

[TABLE]

as a discrete approximation of the identified coefficient $\left(Q^{\dagger},f^{\dagger},g^{\dagger}\right)$, here $\rho>0$ is the regularization parameter. However, due to the non-linearity of the coefficient-to-solution operator, we are faced with certain difficulties in holding the non-convex minimization problem (1.10). Thus, instead of working with the above least squares functional and following the use of energy functions (cf. [37, 35, 48]), in the present work the convex cost function (cf. §2)

[TABLE]

will be taken into account. We then consider a unique minimizer $\left(Q^{h},f^{h},g^{h}\right)$ of the strictly convex problem

[TABLE]

as a discrete regularized solution of the identification problem. Note that, by using variational discretization concept introduced in [22], every solution of the minimization problem (1.11) is proved to automatically belong to finite dimensional spaces. Thus, a discretization of the admissible set $\mathcal{H}_{ad}$ can be avoided. Furthermore, for simplicity of exposition we here restrict ourselves to the case of one set of data $(z_{\delta})_{\delta>0}$. In case with several sets of data $(z_{\delta_{i}})_{i=1}^{I}$ being available, we can replace the misfit term in the problem (1.11) by the term $\dfrac{1}{I}\sum_{i=1}^{I}{\mathcal{J}}^{h}_{\delta_{i}}(Q,f,g)$.

In §3 we will show the convergence of these approximation solutions $\left(Q^{h},f^{h},g^{h}\right)$ to the identification $\left(Q^{\dagger},f^{\dagger},g^{\dagger}\right)$ in the $\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$-norm as well as the convergence of corresponding approximation states \big{(}\mathcal{U}^{h}_{Q^{h},f^{h},g^{h}}\big{)} to the exact $\Phi^{\dagger}$ in the $H^{1}(\Omega)$-norm. Under the structural source condition — but without the smallness requirement — of the general convergence theory for non-linear, ill-posed problems (cf. [15, 16]), we prove in §4 error bounds for these discrete approximations. For the numerical solution of the minimization problem (1.11) we in §5 employ a gradient projection algorithm with Armijo steplength rule. Finally, a numerical implementation will be performed to illustrate the theoretical findings.

The coefficient identification problem in PDEs arises from different contexts of applied sciences, e.g., from aquifer analysis, geophysical prospecting and pollutant detection, and attracted great attention from many scientists in the last 30 years or so. For surveys on the subject one may consult in [3, 9, 29, 43, 45, 46]. The problem of identifying the scalar diffusion coefficient has been extensively studied for both theoretical research and numerical implementation, see e.g., [7, 8, 10, 11, 17, 18, 19, 27, 30, 32, 33, 36, 40, 48]. Some contributions for the case of the simultaneous identification can be found in [2, 20, 21, 34] while some works treated the diffusion matrix case have been obtained in [14, 24, 25, 26, 39].

We conclude this introduction with the following mention. By using the H-convergent concept, the convergence analysis presented in [24] can not be applied directly to the problem of identifying scalar diffusion coefficients. There are two main difficulties for the scalar coefficient identification. First, the set

[TABLE]

is in general not a closed subset of $\mathcal{Q}_{ad}$ under the topology of the H-convergence (cf. [47]), i.e. if the sequence $(q_{n}I_{d})_{n}\subset\mathcal{D}$ is H-convergent to $Q\in\mathcal{Q}_{ad}$, then $Q$ is not necessarily proportional to $I_{d}$ in dimension $d\geq 2$ or $Q\notin\mathcal{D}$. Second, the forward operator $\mathcal{U}$ is not weakly sequentially closed in $L^{2}$, i.e. if $(q_{n},\mathcal{U}(q_{n}))\rightharpoonup(q,\mathcal{Y})$ weakly in $L^{2}(\Omega)\times L^{2}(\Omega)$, it is not guaranteed that $\mathcal{Y}=\mathcal{U}(q)$ (see [14] and the references therein for counterexamples). To overcome these difficulties, a different analysis technique based on the convexity of the cost functional will be taken into counting. Due to the weak∗ closedness of the set $\mathcal{D}$ above in $\mathbf{L}^{\infty}_{\mbox{\tiny sym}}(\Omega)$ (cf. Remark 2.1), the convergence analysis performed in the present paper thus covers the scalar diffusion identification case. On the other hand, in [24] the source term and the boundary condition were assumed to be given. In the present situation they are variables which have to be found simultaneously together with the diffusion from observations.

Throughout the paper we write $\int_{\Omega}\cdots$ instead of $\int_{\Omega}\cdots dx$ for the convenience of relevant notations. We use the standard notion of Sobolev spaces $H^{1}(\Omega)$, $H^{2}(\Omega)$, $W^{k,p}(\Omega)$, etc from, e.g., [1].

2 Finite element discretization

2.1 Preliminaries

In product spaces $\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$ and $\mathbf{L}^{\infty}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$ we use respectively the norm

[TABLE]

We note that the coefficient-to-solution operator

[TABLE]

with

[TABLE]

is Fréchet differentiable on $\mathcal{H}_{ad}$. For each $\Gamma=(Q,f,g)\in\mathcal{H}_{ad}$ the action of its Fréchet derivative in direction $\lambda:=(H,l,s)\in\mathbf{L}^{\infty}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$ denoted by $\xi_{\lambda}:=\mathcal{U}^{\prime}_{\Gamma}(\lambda):=\mathcal{U}^{\prime}(\Gamma)(\lambda)$ is the unique weak solution in $H^{1}_{\diamond}(\Omega)$ to the equation

[TABLE]

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

In $\mathcal{S}_{d}$ we introduce the convex subset

[TABLE]

together with the orthogonal projection $P_{\mathcal{K}}:\mathcal{S}_{d}\to\mathcal{K}$ that is characterised by

[TABLE]

for all $A\in\mathcal{S}_{d}$ and $B\in\mathcal{K}$. Furthermore, let $\xi:=(\xi_{1},\cdot\cdot\cdot,\xi_{d})$ and $\eta:=(\eta_{1},\cdot\cdot\cdot,\eta_{d})$ be two arbitrary vectors in $R^{d}$, we use the notation

[TABLE]

We close this subsection by the following note.

Remark 2.1**.**

Let

[TABLE]

Then $\mathbf{D}$ is a weakly∗ compact subset of $L^{\infty}(\Omega)$, i.e. for any sequence $(q_{n})_{n}\subset\mathbf{D}$ a subsequence $(q_{n_{m}})_{m}$ and an element $\xi_{\infty}\in\mathbf{D}$ exist such that $(q_{n_{m}})_{m}$ is weakly∗ convergent in $L^{\infty}(\Omega)$ to $\xi_{\infty}$. In other words, for all $\theta_{1}\in L^{1}(\Omega)$ there holds the limit

[TABLE]

We also remark that any $\Psi\in L^{\infty}(\Omega)$ can be considered as an element in ${L^{\infty}(\Omega)}^{*}$ by

[TABLE]

for all $\psi\mbox{~{}in~{}}L^{\infty}(\Omega)$ and $\|\Psi\|_{{L^{\infty}(\Omega)}^{*}}\leq|\Omega|\cdot\|\Psi\|_{L^{\infty}(\Omega)}$. Therefore, due to (2.2), the assertion of Remark 2.1 is a direct consequence of the Banach-Alaoglu theorem.

2.2 Discretization

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

[TABLE]

with $\mathcal{P}_{r}$ consisting all polynomial functions of degree at most $r$. Similar to the continuous case, we have the following result.

Lemma 2.2**.**

Let $(Q,f,g)$ be in $\mathcal{H}_{ad}$. Then the variational equation

[TABLE]

for all $\varphi^{h}\in\mathcal{V}^{h}_{1}$ admits a unique solution $\Phi^{h}\in\mathcal{V}^{h}_{1}$. Furthermore, the priori estimate

[TABLE]

is satisfied.

The map $\mathcal{U}^{h}:\mathcal{H}_{ad}\subset\mathbf{L}^{\infty}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)\rightarrow\mathcal{V}^{h}_{1}$ from each $\Gamma:=(Q,f,g)\in\mathcal{H}_{ad}$ to the unique solution $\mathcal{U}^{h}_{\Gamma}:=\Phi^{h}$ of (2.4) is called the discrete coefficient-to-solution operator. This operator is also Fréchet differentiable on the set $\mathcal{H}_{ad}$. For each $\Gamma=(Q,f,g)\in\mathcal{H}_{ad}$ and $\lambda:=(H,l,s)\in\mathbf{L}^{\infty}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$ the Fréchet differential $\xi^{h}_{\lambda}:={\mathcal{U}^{h}_{\Gamma}}^{\prime}(\lambda)$ is an element of $\mathcal{V}^{h}_{1}$ and satisfies for all $\varphi^{h}$ in $\mathcal{V}^{h}_{1}$ the equation

[TABLE]

Due to the standard theory of the finite element method for elliptic problems (cf. [6, 12]), for any fixed $\Gamma=(Q,f,g)\in\mathcal{H}_{ad}$ there holds the limit

[TABLE]

Let

[TABLE]

be the Clément’s mollification interpolation operator with properties

[TABLE]

and

[TABLE]

for $0\leq k\leq l\leq 2$, where $C$ is independent of $h$ and $\phi$ (cf. [13, 4, 5, 44]). Then, using the discrete operator $\mathcal{U}^{h}$ and the interpolation operator $\Pi^{h}$, we can now introduce the discrete cost functional

[TABLE]

where $(Q,f,g)\in\mathcal{H}_{ad}$.

Lemma 2.3**.**

Assume that the sequence $\left(\Gamma_{n}\right)_{n}:=\left(Q_{n},f_{n},g_{n}\right)_{n}\subset\mathcal{H}_{ad}$ weakly converges to $\Gamma:=(Q,f,g)$ in $\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$. Then for any fixed $h>0$ the sequence $\left(\mathcal{U}^{h}_{\Gamma_{n}}\right)_{n}\subset\mathcal{V}^{h}_{1}$ converges to $\mathcal{U}^{h}_{\Gamma}$ in the $H^{1}(\Omega)$-norm.

Proof.

Due to Remark 2.1, $(Q_{n})_{n}$ has a subsequence denoted by the same symbol which is weakly∗ convergent in $\mathbf{L}^{\infty}_{\mbox{\tiny sym}}(\Omega)$ to $Q$. Furthermore, by (2.5), the corresponding state sequence $\left(\mathcal{U}^{h}_{\Gamma_{n}}\right)_{n}$ is bounded in the finite dimensional space $\mathcal{V}^{h}_{1}$. A subsequence which is not relabelled and an element $\Theta^{h}\in\mathcal{V}^{h}_{1}$ then exist such that $\left(\mathcal{U}^{h}_{\Gamma_{n}}\right)_{n}$ converges to $\Theta^{h}$ in the $H^{1}(\Omega)$-norm. It follows from the equation (2.4) that

[TABLE]

for all $\varphi^{h}\in\mathcal{V}^{h}_{1}$. Taking $\varphi^{h}=\mathcal{U}^{h}_{\Gamma_{n}}-\mathcal{U}^{h}_{\Gamma}$, by (1.5), we obtain that

[TABLE]

Since $Q_{n}\rightharpoonup Q$ weakly∗ in $\mathbf{L}^{\infty}_{\mbox{\tiny sym}}(\Omega)$, we get $\lim_{n\to\infty}\int_{\Omega}\left(Q-Q_{n}\right)\nabla\mathcal{U}^{h}_{\Gamma}\cdot\nabla\left(\Theta^{h}-\mathcal{U}^{h}_{\Gamma}\right)=0.$ Sending $n$ to $\infty$, we thus obtain from the last inequality that $\lim_{n\to\infty}\left\|\mathcal{U}^{h}_{\Gamma_{n}}-\mathcal{U}^{h}_{\Gamma}\right\|_{H^{1}(\Omega)}=0$, which finishes the proof. ∎

We now state the following useful result on the convexity of the cost functional.

Lemma 2.4**.**

$\mathcal{J}^{h}_{\delta}$* is convex and continuous on $\mathcal{H}_{ad}$ with respect to the $\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$-norm.*

Proof.

The continuity of $\mathcal{J}^{h}_{\delta}$ follows directly from Lemma 2.3. We show that $\mathcal{J}^{h}_{\delta}$ is convex.

Let $\Gamma:=(Q,f,g)\in\mathcal{H}_{ad}$ and $\lambda:=(H,l,s)\in\mathbf{L}^{\infty}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$. We have that

[TABLE]

We compute for each term in the right hand side of the last equation. First we get

[TABLE]

For the second term we have

[TABLE]

Finally, we have for the third term

[TABLE]

Therefore,

[TABLE]

where

[TABLE]

By (2.6), we infer that

[TABLE]

Therefore, by (2.6) again, we arrive at

[TABLE]

by (1.5), which completes the proof. ∎

Now we are in position to prove the main result of this section.

Theorem 2.5**.**

The strictly convex minimization problem

[TABLE]

attains a unique minimizer. Furthermore, an element $\Gamma:=(Q,f,g)\in\mathcal{H}_{ad}$ is the unique minimizer to $\left(\mathcal{P}^{\rho,h}_{\delta}\right)$ if and only if the system

[TABLE]

holds for a.e. in $\Omega$, where $\bar{\Pi}^{h}$ was generated from $\Pi^{h}$ according to (2.13).

Proof.

Let $(\Gamma_{n})_{n}:=(Q_{n},f_{n},g_{n})_{n}\subset\mathcal{H}_{ad}$ be a minimizing sequence of $\left(\mathcal{P}^{\rho,h}_{\delta}\right)$, i.e.

[TABLE]

The sequence $(\Gamma_{n})_{n}$ is thus bounded in the $\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$-norm. A subsequence not relabelled and an element $\Gamma:=(Q,f,g)\in\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$ exist such that $\Gamma_{n}\rightharpoonup\Gamma$ weakly in $\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$. On the other hand, since $\mathcal{H}_{ad}$ is a convex, closed subset of $\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$, so is weakly closed, it follows that $\Gamma\in\mathcal{H}_{ad}$. By Lemma 2.4, $\mathcal{J}^{h}_{\delta}$ and $\mathcal{R}$ are both weakly lower semi-continuous on $\mathcal{H}_{ad}$ which yields that

[TABLE]

We therefore have that

[TABLE]

and $\Gamma$ is then a minimizer to $\left(\mathcal{P}^{\rho,h}_{\delta}\right)$. Since $\Upsilon^{\rho,h}_{\delta}$ is strictly convex, this minimizer is unique. Next, an element $\Gamma:=(Q,f,g)\in\mathcal{H}_{ad}$ is the minimizer to $\left(\mathcal{P}^{\rho,h}_{\delta}\right)$ if and only if ${\Upsilon^{\rho,h}_{\delta}}^{\prime}(\Gamma)(\overline{\Gamma}-\Gamma)\geq 0$ for all $\overline{\Gamma}=(H,l,s)\in\mathcal{H}_{ad}$. Then, in view of (2.2), we get that

[TABLE]

for all $\overline{\Gamma}=(H,l,s)\in\mathcal{H}_{ad}$. Taking $\overline{\Gamma}_{1}=(H,f,g)$, $\overline{\Gamma}_{2}=(Q,l,g)$ and $\overline{\Gamma}_{3}=(Q,f,s)$ into the above inequality we obtain the system (2.15)–(2.17). The proof is completed. ∎

Remark 2.6**.**

We denote by

[TABLE]

Since $\mathcal{U}^{h}_{\Gamma}\in\mathcal{V}^{h}_{1}$ and $\bar{\Pi}^{h}z_{\delta}\in\mathcal{V}^{h}_{1}$, the system (2.15)–(2.17) shows that every solution of \big{(}\mathcal{P}^{\rho,h}_{\delta}\big{)} automatically belongs to the finite dimensional space ${\mathcal{V}^{h}_{0}}^{d\times d}\times\mathcal{V}^{h}_{1}\times\mathcal{E}^{h}_{1}$.

3 Convergence

For abbreviation in what follows we denote by $C$ a generic positive constant independent of the mesh size $h$, the noise level $\delta$ and the regularization parameter $\rho$. By (2.8) and (2.9), we can introduce for each $\Phi\in H^{1}(\Omega)$

[TABLE]

in case $\Phi\in H^{2}(\Omega)$. Likewise, by (2.7), for all $\Gamma\in\mathcal{H}_{ad}$

[TABLE]

Furthermore, by (2.9), we get

[TABLE]

Thus, it follows from the inverse inequality (cf. [6, 12]):

[TABLE]

that

[TABLE]

The following result shows the convergence of finite element approximations to the unique minimum norm solution $\Gamma^{\dagger}:=\left(Q^{\dagger},f^{\dagger},g^{\dagger}\right)$ of the identification problem, which is defined by (1.8).

Theorem 3.1**.**

Let $\left(h_{n}\right)_{n}$ be a sequence with $\lim_{n\to\infty}h_{n}=0$ and $\left(\delta_{n}\right)_{n}$ and $\left(\rho_{n}\right)_{n}$ are any positive sequences such that

[TABLE]

Assume that $\left(z_{\delta_{n}}\right)_{n}\subset L^{2}(\Omega)$ is a sequence satisfying $\left\|z_{\delta_{n}}-\Phi^{\dagger}\right\|_{L^{2}(\Omega)}\leq\delta_{n}$ and $\Gamma_{n}:=\left(Q_{n},f_{n},g_{n}\right)$ is the unique minimizer of the problem \big{(}\mathcal{P}^{\rho_{n},h_{n}}_{\delta_{n}}\big{)} for each $n\in N$. Then the sequence $\left(\Gamma_{n}\right)_{n}$ converges to $\Gamma^{\dagger}$ in the $\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$-norm as $n\to\infty$. Furthermore, the corresponding discrete state sequence \big{(}\mathcal{U}^{h_{n}}_{\Gamma_{n}}\big{)}_{n} also converges to $\Phi^{\dagger}$ in the $H^{1}(\Omega)$-norm.

Remark 3.2**.**

In case $\Phi^{\dagger}=\mathcal{U}_{\Gamma^{\dagger}}\in H^{2}(\Omega)$ we have $0\leq\beta^{h_{n}}_{\mathcal{U}_{\Gamma^{\dagger}}},\chi^{h_{n}}_{\Phi^{\dagger}}\leq Ch_{n}$. Therefore, the convergence of Theorem 3.1 is obtained if $\delta_{n}\sim h_{n}^{2}$ and the sequence $(\rho_{n})_{n}$ is chosen such that

[TABLE]

To prove Theorem 3.1, we need the following auxiliary estimate.

Lemma 3.3**.**

There holds the estimate

[TABLE]

Proof.

We have with $\Phi^{\dagger}=\mathcal{U}_{\Gamma^{\dagger}}$ and (3.2) that

[TABLE]

which finishes the proof. ∎

Proof of Theorem 3.1.

By the optimality of $\Gamma_{n}$ and Lemma 3.3, we have that

[TABLE]

which yields

[TABLE]

and

[TABLE]

A subsequence of the sequence $\left(\Gamma_{n}\right)_{n}$ denoted by the same symbol and an element $\Gamma_{0}:=(Q_{0},f_{0},g_{0})\in\mathcal{H}_{ad}$ then exist such that

[TABLE]

We will show that $\left(\Gamma_{n}\right)_{n}$ converges to $\Gamma_{0}$ in the $\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$-norm and $\Gamma_{0}=\Gamma^{\dagger}$. We have from (3.2) that

[TABLE]

Combining this with $\lim_{n\to\infty}\|\mathcal{U}_{\Gamma_{0}}-\mathcal{U}^{h_{n}}_{\Gamma_{0}}\|_{H^{1}(\Omega)}=0$ from (2.7), we arrive at

[TABLE]

Now for each fixed $n$ we consider an arbitrary subsequence $(\Gamma_{n_{m}})_{m}$ of $(\Gamma_{n})_{n}$. By the weakly l.s.c. property of the functional $\mathcal{J}^{h_{n}}_{\delta_{n}}$ (cf. Lemma 2.4), we obtain that

[TABLE]

Again, using the convexity of $\mathcal{J}^{h_{n}}_{\delta_{n}}$, we get that

[TABLE]

By (1.5), we thus arrive at

[TABLE]

Using (3.4), we infer from the last inequality that

[TABLE]

In view of (2.2) we get that

[TABLE]

Since $Q_{n_{m}}\rightharpoonup Q_{0}$ weakly∗ in $\mathbf{L}^{\infty}_{\mbox{\tiny sym}}(\Omega)$ as $m\to\infty$, we have for the first term that

[TABLE]

since $\lim_{n\to\infty}\int_{\Omega}\left(Q_{0}-Q_{n}\right)\nabla\mathcal{U}_{\Gamma^{\dagger}}\cdot\nabla\mathcal{U}_{\Gamma^{\dagger}}=0$, due to $Q_{n}\rightharpoonup Q_{0}$ weakly∗ in $\mathbf{L}^{\infty}_{\mbox{\tiny sym}}(\Omega)$. Furthermore, by (3.6), we get that

[TABLE]

Therefore,

[TABLE]

On the other hand, we get

[TABLE]

We now have that

[TABLE]

with

[TABLE]

and

[TABLE]

so that

[TABLE]

Next, we rewrite

[TABLE]

By (2.7), likewise as (3.9), we get that

[TABLE]

Furthermore, we have

[TABLE]

By Lemma 2.3, for each fixed $n$ we have that the sequence $\left(\mathcal{U}^{h_{n}}_{\Gamma_{n_{m}}}\right)_{m}\subset\mathcal{V}^{h_{n}}_{1}$ converges to $\mathcal{U}^{h_{n}}_{\Gamma_{0}}$ in the $H^{1}(\Omega)$-norm as $m$ tends to $\infty$. Then we deduce that

[TABLE]

Thus, we obtain

[TABLE]

Finally, we also get that

[TABLE]

and

[TABLE]

Therefore, it follows from the equations (3.8)–(3.14) that

[TABLE]

Combining this with (3.7), we obtain that $\mathcal{U}_{\Gamma_{0}}=\mathcal{U}_{\Gamma^{\dagger}}.$ Then, by the definition of $\Gamma^{\dagger}$, the weakly l.s.c. property of $\mathcal{R}$ and (3.5), we get

[TABLE]

Thus, $\mathcal{R}(\Gamma^{\dagger})=\mathcal{R}(\Gamma_{0})=\lim_{n\to\infty}\mathcal{R}\left(\Gamma_{n}\right).$ By the uniqueness of $\Gamma^{\dagger}$, we have $\Gamma_{0}=\Gamma^{\dagger}$. Furthermore, since $\left(\Gamma_{n}\right)_{n}$ weakly converges in $\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$ to $\Gamma_{0}$, we conclude from the last equation that $\left(\Gamma_{n}\right)_{n}$ converges to $\Gamma_{0}$ in the $\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)\times L^{2}(\Omega)\times L^{2}(\partial\Omega)$-norm.

It remains to show that the sequence \big{(}\mathcal{U}^{h_{n}}_{\Gamma_{n}}\big{)}_{n} converges to $\Phi^{\dagger}=\mathcal{U}_{\Gamma^{\dagger}}$ in the $H^{1}(\Omega)$-norm. We first get from (2.7) that

[TABLE]

Furthermore, in view of (2.2) we also have that

[TABLE]

Since $f_{n}\rightarrow f^{\dagger}$ in the $L^{2}(\Omega)$-norm and $g_{n}\rightarrow g^{\dagger}$ in the $L^{2}(\partial\Omega)$-norm together with the uniform boundedness (2.5), it follows that

[TABLE]

We now rewrite

[TABLE]

We will estimate for two terms in the right hand side of the above equation. For simplicity of notation we here set

[TABLE]

Then, we have

[TABLE]

Similarly, we get

[TABLE]

and arrive at

[TABLE]

Since $Q_{n}\rightarrow Q^{\dagger}$ in the $\mathbf{L}^{2}_{\mbox{\tiny sym}}(\Omega)$-norm, up to a subsequence we assume that $(Q_{n})_{n}$ converges to $Q^{\dagger}$ a.e. in $\Omega$. Then, by the Lebesgue dominated convergence theorem, we deduce that

[TABLE]

Thus, together with (3.15), we have

[TABLE]

It follows from (3)–(3.18) that \lim_{n\to\infty}\big{\|}\mathcal{U}^{h_{n}}_{\Gamma_{n}}-\mathcal{U}^{h_{n}}_{\Gamma^{\dagger}}\big{\|}_{H^{1}(\Omega)}=0. By serving of (3.15) again, we then conclude that \lim_{n\to\infty}\big{\|}\mathcal{U}^{h_{n}}_{\Gamma_{n}}-\mathcal{U}_{\Gamma^{\dagger}}\big{\|}_{H^{1}(\Omega)}=0, which finishes the proof. ∎

4 Error bounds

In this section we investigate error bounds of discrete regularized solutions to the identification problem. For any $\Gamma:=(Q,f,g)\in\mathcal{H}_{ad}$ the mapping

[TABLE]

is linear, continuous with the dual

[TABLE]

Theorem 4.1**.**

Assume that a function $w^{*}\in{H^{1}_{\diamond}(\Omega)}^{*}$ exists such that

[TABLE]

Then

[TABLE]

where \Gamma^{h}:=\big{(}Q^{h},f^{h},g^{h}\big{)} is the unique solution to $\left(\mathcal{P}^{\rho,h}_{\delta}\right)$ and $w\in H^{1}_{\diamond}(\Omega)$ is the unique weak solution of the Neumann problem

[TABLE]

Remark 4.2**.**

Due to Remark 3.2, in case $\mathcal{U}_{\Gamma^{\dagger}},~{}w\in H^{2}(\Omega)$ we have $0\leq\chi^{h}_{\Phi^{\dagger}},~{}\beta^{h}_{\mathcal{U}_{\Gamma^{\dagger}}},~{}\chi^{h}_{w}\leq Ch.$ Therefore, with $\delta\sim h^{2}$ and $\rho\sim h$ we obtain the following error bounds

[TABLE]

Remark 4.3**.**

Let $\bar{\Gamma}:=(\bar{Q},\bar{f},\bar{g})\in\mathcal{I}(\Phi^{\dagger})$ be such that the equation (4.1) satisfying with $\bar{\Gamma}$ for some $w^{*}\in{H^{1}_{\diamond}(\Omega)}^{*}$, i.e. ${\mathcal{U}^{\prime}_{\bar{\Gamma}}}^{*}w^{*}=\bar{\Gamma}.$ Then $\bar{\Gamma}$ is the unique minimum norm solution of the identification, i.e. $\bar{\Gamma}=\Gamma^{\dagger}$.

Indeed, due to (2.2) we have for all $\Gamma:=(Q,f,g)\in\mathcal{I}(\Phi^{\dagger})$ that

[TABLE]

for some $W\in H^{1}_{\diamond}(\Omega)$, since the expression $[u,v]:=\int_{\Omega}\bar{Q}\nabla u\cdot\nabla v$ generates a scalar inner product on the space $H^{1}_{\diamond}(\Omega)$ which is equivalent to the usual one. By (2.1) we then get

[TABLE]

due to (1.6) and the fact $\mathcal{U}_{\bar{\Gamma}}=\mathcal{U}_{\Gamma}=\Phi^{\dagger}$. Therefore, we deduce that

[TABLE]

which completed the proof.

Proof of Theorem 4.1.

Due to the optimality of $\Gamma^{h}$, we get that

[TABLE]

which implies

[TABLE]

by Lemma 3.3. Now, by (2.2) and (4.1), we infer that

[TABLE]

Thus, by the definition of the weak solution to (4.3) and (2.1), we obtain

[TABLE]

which yields

[TABLE]

For $I_{1}$ we have from (3.2) that

[TABLE]

Due to (1.6) and (2.4), we get \int_{\Omega}Q^{h}\nabla\left(\mathcal{U}^{h}_{\Gamma^{h}}-\mathcal{U}_{\Gamma^{h}}\right)\cdot\nabla\Pi^{h}\big{(}w-\mathcal{U}_{\Gamma^{\dagger}}\big{)}=0 and then infer that

[TABLE]

Finally, we have that

[TABLE]

It follows from (4.8)–(4.11) that

[TABLE]

Thus, together with (4)–(4.7), we get

[TABLE]

which finishes the proof. ∎

5 Gradient projection algorithm with Armijo steplength rule

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

We first note that for each $\Gamma:=(Q,f,g)\in\mathcal{H}_{ad}$, in view of (2.2), the $\mathcal{L}^{2}$-gradient of the strictly convex cost function $\Upsilon^{\rho,h}_{\delta}$ of the problem \big{(}\mathcal{P}^{\rho,h}_{\delta}\big{)} is given by \nabla\Upsilon^{\rho,h}_{\delta}(\Gamma):=\big{(}\Upsilon_{Q}(\Gamma),\Upsilon_{f}(\Gamma),\Upsilon_{g}(\Gamma)\big{)} with

[TABLE]

and $\bar{\Pi}^{h}$ generating from $\Pi^{h}$ according to (2.13).

The algorithm is then read as: given a step size control $\beta\in(0,1)$, an initial approximation (cf. Remark 2.6) \Gamma_{0}:=(Q_{0},f_{0},g_{0})\in\mathcal{H}_{ad}\cap\big{(}{\mathcal{V}^{h}_{0}}^{d\times d}\times\mathcal{V}^{h}_{1}\times\mathcal{E}^{h}_{1}\big{)}, number of iteration $N$ and setting $k=0$.

Compute $\mathcal{U}^{h}_{\Gamma_{k}}$ from the variational equation

[TABLE]

as well as

[TABLE] 2. 2.

Compute the gradient \nabla\Upsilon^{\rho,h}_{\delta}(\Gamma_{k}):=\big{(}\Upsilon_{Q_{k}}(\Gamma_{k}),\Upsilon_{f_{k}}(\Gamma_{k}),\Upsilon_{g_{k}}(\Gamma_{k})\big{)} with

[TABLE] 3. 3.

Set \widetilde{\Gamma}_{k}:=\big{(}\widetilde{Q}_{k},\widetilde{f}_{k},\widetilde{g}_{k}\big{)} with \widetilde{Q}_{k}(x):=P_{\mathcal{K}}\big{(}Q_{k}(x)-\beta\Upsilon_{Q_{k}}(\Gamma_{k})(x)\big{)}, $\widetilde{f}_{k}(x):=f_{k}(x)-\beta\Upsilon_{f_{k}}(\Gamma_{k})(x)$ and $\widetilde{g}_{k}(x):=g_{k}(x)-\beta\Upsilon_{g_{k}}(\Gamma_{k})(x)$.

1. (a)

Compute $\mathcal{U}^{h}_{\widetilde{\Gamma}_{k}}$ according to (5.1), \Upsilon_{\rho,\delta}^{h}\big{(}\widetilde{\Gamma}_{k}\big{)} according to (5.2), and with $\tau=10^{-4}$

[TABLE] 2. (b)

If $L\leq 0$

go to the next step (c) below

else

set $\beta:=\frac{\beta}{2}$ and then go back (a) 3. (c)

Update $\Gamma_{k}=\widetilde{\Gamma}_{k}$, set $k=k+1$. 4. 4.

Compute

[TABLE]

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

6 Numerical implementation

For illustrating the theoretical result we consider the Neumann problem

[TABLE]

with $\Omega=\{x=(x_{1},x_{2})\in{R}^{2}~{}|~{}-1<x_{1},x_{2}<1\}$.

The special constants in the equation (1.3) are chosen as $\underline{q}=0.05$ and $\overline{q}=10$. 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}$.

We assume that entries of the symmetric diffusion matrix $Q^{\dagger}$ are discontinuous which are defined as

[TABLE]

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

[TABLE]

The source functional $f^{\dagger}$ is assumed to be also discontinuous and defined as

[TABLE]

where

[TABLE]

The Neumann boundary condition $g^{\dagger}$ is chosen with

[TABLE]

The exact state $\Phi^{\dagger}$ is then computed from the finite element equation $KU=F$, where $K$ and $F$ are the stiffness matrix and the load vector associated with the problem (6.1)–(6.2), respectively.

We mention that in the above example the sought functions are chosen to be discontinuous. To reconstruct such discontinuous functions one usually employs the total variation regularization which was originally introduced in image denoising by authors of [41]. This regularization method was proved to be very effective and analyzed by many authors over the last decades for several ill-posed and inverse problems. We also note that the space of all functions with bounded total variation is a non-reflexive Banach space and the Tikhonov-function of the total variation regularization is non-differentiable, which cause some certain difficulties in numerically treating for non-linear, ill-posed inverse problems. In the present work the cost function is convex and differentiable, the convergence history given in Table 1 and Table 2 below shows that the algorithm presented in Section 5 performs well for the identification problem with the discontinuous coefficents.

We start the computation with the coarsest level $\ell=3$. To this end, for constructing observations with noise of the exact state $\Phi^{\dagger}$ on this coarsest grid we use

[TABLE]

where $\overline{\delta_{\ell}}=10\rho_{\ell}^{1/2}h_{\ell}^{3/2}$, $\rho_{\ell}=10^{-3}h_{\ell}$ and $\mathcal{N}_{\overline{\delta_{\ell}}}$ is a $M^{h_{\ell}}\times 1$-matrix of random numbers in the interval \big{(}-\overline{\delta_{\ell}},\overline{\delta_{\ell}}\big{)}, $M^{h_{\ell}}=(\ell+1)^{2}$ is the number of nodes of the triangulation $\mathcal{T}^{h_{\ell}}$. Therefore, the exact state $\Phi^{\dagger}$ is only measured at 16 nodes of $\mathcal{T}^{h_{\ell}}$.

We use the algorithm described in §5 for computing the numerical solution of the problem \big{(}\mathcal{P}_{\rho_{\ell},\delta_{\ell}}^{h_{\ell}}\big{)}. The step size control is chosen with $\beta=0.75$. As the initial approximation we choose

[TABLE]

At each iteration $k$ we compute Tolerance defined by (5.3). Then the iteration was stopped if $\mbox{Tolerance}\leq 0$ or the number of iterations reached the maximum iteration count of 800.

After obtaining the numerical solution $\Gamma_{\ell}=(Q_{\ell},f_{\ell},g_{\ell})$ and the computed numerical state $\mathcal{U}_{\ell}=\mathcal{U}^{h_{\ell}}_{\Gamma_{\ell}}$ of the first iteration process with respect to the coarsest level $\ell=3$, we use their interpolations on the next finer mesh $\ell=6$ as an initial approximation and an observation of the exact state for the algorithm on this finer mesh, i.e. for the next iteration process with respect to the level $\ell=6$ we employ

[TABLE]

and $I^{h_{\ell}}_{1}$ being the usual node value interpolation operator on $\mathcal{T}^{h_{\ell}}$, and so on $\ell=12,24,\ldots$. We note that the computation process only requires the measurement data of the exact data for the coarsest level $\ell=3$.

The numerical results are summarized in Table 1 and Table 2, where we present the refinement level $\ell$, mesh size $h_{\ell}$ of the triangulation, regularization parameter $\rho_{\ell}$, noise $\delta_{\ell}$ and number of iterates as well as the final $L^{2}$-error in the coefficients, the final $L^{2}$ and $H^{1}$-error in the states, and their experimental order of convergence (EOC), where $\mbox{EOC}_{\Phi}:=\frac{\ln\Phi(h_{1})-\ln\Phi(h_{2})}{\ln h_{1}-\ln h_{2}}$ and $\Phi(h)$ is an error function with respect to $h$.

All figures are here presented corresponding to $\ell=96$. Figure 1 from left to right shows the graphs of $\Phi^{\dagger}$, computed numerical state $\mathcal{U}_{\ell}$ of the algorithm at the last iteration, and the difference to $\Phi^{\dagger}$. In Figure 2 we display the computed numerical source term and boundary condition $f_{\ell}$, $g_{\ell}$ at the last iteration as well as the differences $f_{\ell}-f^{\dagger}$, $g_{\ell}-g^{\dagger}$. We write the computed numerical diffusion matrix at the last iteration as

[TABLE]

Figure 3 then shows $q_{\ell,11}$, $q_{\ell,12}$ and $q_{\ell,22}$ while Figure 4 shows differences $q_{\ell,11}-q^{\dagger}_{11}$, $q_{\ell,12}-q^{\dagger}_{12}$ and $q_{\ell,22}-q^{\dagger}_{22}$. For abbreviation we denote by \Gamma^{\dagger}:=\big{(}Q^{\dagger},f^{\dagger},g^{\dagger}\big{)} and errors

[TABLE]

Acknowledgments

The author would like to thank the Referees and the Editor for their valuable comments and suggestions which helped to improve the present paper.