Convergence analysis of a Crank-Nicolson Galerkin method for an inverse source problem for parabolic equations with boundary observations

TL;DR

This paper analyzes the convergence of a Crank-Nicolson Galerkin method for an inverse source problem in parabolic equations, providing theoretical proofs and numerical validation for the method's accuracy and stability.

Contribution

It introduces a convergence analysis for a novel Crank-Nicolson Galerkin approach to inverse source problems with boundary data, including error bounds and rates.

Findings

Proved convergence of the regularized approximations as noise and mesh size decrease.

Established error bounds and convergence rates under source conditions.

Numerical experiments confirm theoretical results.

Abstract

This work is devoted to an inverse problem of identifying a source term depending on both spatial and time variables in a parabolic equation from single Cauchy data on a part of the boundary. A Crank-Nicolson Galerkin method is applied to the least squares functional with an quadratic stabilizing penalty term. The convergence of finite dimensional regularized approximations to the sought source as measurement noise levels and mesh sizes approach to zero with an appropriate regularization parameter is proved. Moreover, under a suitable source condition, an error bound and corresponding convergence rates are proved. Finally, several numerical experiments are presented to illustrate the theoretical findings.

Tables

Table 1. Table 1: Space-dependent source: Refinement level l 𝑙 l , mesh size h ℎ h , measurement error δ 𝛿 \delta , regularization parameter ρ 𝜌 \rho , and errors

$l$	$h$	$\delta$	$\rho$	${\Vert u-u_{rec.}\Vert }_{L^2\left({\Omega }_T\right)}$	${\Vert u-u_{rec.}\Vert }_{L^2\left(\Sigma \right)}$	${\Vert f-f_{rec.}\Vert }_{L^2\left({\Omega }_T\right)}$
1	0.8	0.32	0.008	0.2160	0.2762	1.1199
2	0.4	0.08	0.004	0.0534	0.0727	0.3262
3	0.2	0.02	0.002	0.0132	0.0190	0.1083
4	0.1	0.005	0.001	0.0029	0.0049	0.0556

Table 2. Table 2: Space-dependent source: EOC

$l$	${\Vert u-u_{rec.}\Vert }_{L^2\left({\Omega }_T\right)}$	${\Vert u-u_{rec.}\Vert }_{L^2\left(\Sigma \right)}$	${\Vert f-f_{rec.}\Vert }_{L^2\left({\Omega }_T\right)}$
1	–	–	–
2	2.0161	1.9257	1.7795
3	2.0163	1.9360	1.5907
4	2.1864	1.9551	0.9619
Mean of EOC	2.0729	1.9389	1.4440

Table 3. Table 3: Identified source satisfying the condition Eq.   57 : Refinement level l 𝑙 l , mesh size h ℎ h , measurement error δ 𝛿 \delta , regularization parameter ρ 𝜌 \rho , and errors

$l$	$h$	$\delta$	$\rho$	${\Vert u-u_{rec.}\Vert }_{L^2\left({\Omega }_T\right)}$	${\Vert u-u_{rec.}\Vert }_{L^2\left(\Sigma \right)}$	${\Vert f-f_{rec.}\Vert }_{L^2\left({\Omega }_T\right)}$
1	0.8	0.32	0.008	0.2199	0.2709	1.0994
2	0.4	0.08	0.004	0.0543	0.0699	0.3038
3	0.2	0.02	0.002	0.0136	0.0175	0.0842
4	0.1	0.005	0.001	0.0032	0.0042	0.0237
5	0.05	0.00125	0.0005	0.000745	0.000929	0.009616

Table 4. Table 4: Identified source satisfying the condition Eq.   57 : EOC

$l$	${\Vert u-u_{rec.}\Vert }_{L^2\left({\Omega }_T\right)}$	${\Vert u-u_{rec.}\Vert }_{L^2\left(\Sigma \right)}$	${\Vert f-f_{rec.}\Vert }_{L^2\left({\Omega }_T\right)}$
1	–	–	–
2	2.0178	1.9544	1.8555
3	1.9973	1.9979	1.8512
4	2.0875	2.0589	1.8289
5	2.1028	2.1766	1.3014
Mean of EOC	2.0513	2.0470	1.7093

Equations

$$\frac{\partial u}{\partial t}(x,t)+\mathcal{L}u(x,t)=f(x,t)\text{ in }{\Omega }_T:=\Omega \times (0,T],$$

$$\frac{\partial u(x,t)}{\partial \vec{n}}+\sigma (x,t)u(x,t)=g(x,t)\text{ on }\mathcal{S}:=\partial \Omega \times (0,T],$$

$$u(x,0)=q(x)\text{ in }\Omega$$

$$\Vert Z-z_{\delta }{\Vert }_{L^2(\Sigma )}\le \delta ,$$

$$\mathcal{L}u(x,t):=-\sum _{i,j=1}^d\frac{\partial }{\partial x_i}\left(a_{ij}(x,t)\frac{\partial u(x,t)}{\partial x_j}\right)+b(x,t)u(x,t),$$

$$A(x,t)\xi \cdot \xi =\sum _{i,j=1}^d{a}_{ij}(x,t){\xi }_i{\xi }_j\ge \underline{a}\sum _{i,j=1}^d|{\xi }_i{|}^2\text{mbox}in\overline{{\Omega }_T}$$

$$\frac{\partial u(x,t)}{\partial \vec{n}}:=A(x,t)\nabla u(x,t)\cdot \vec{n}$$

$$t^n:=n\tau \text{mbox}withn\in I_0:=\{0,1,\dots ,M\}.$$

$$J_{\rho ,\delta ,h,\tau }(f):=\sum _{n=1}^M{\int }_{t^{n-1}}^{t^n}\Vert U_{h,\tau }^n(f)-z_{\delta }{\Vert }_{L^2(\Gamma )}^{2}dt+\rho \Vert f-f^{\ast }{\Vert }_{L^2({\Omega }_T)}^2$$

$${f_{\rho ,\delta ,h,\tau }}_{|\Omega \times (t^{n-1},t^n]}=f^{\ast }-{\rho }^{-1}{P}_{h,\tau }^{n-1}(f_{\rho ,\delta ,h,\tau })$$

$${\mathcal{V}}_{h,\tau }^{1,0}:=\left\{\Phi \in L^2\left(0,T;{\mathcal{V}}_h^1\right)|{\Phi }_{|\Omega \times (t^{n-1},t^n]}:={\phi }_h^n\in {\mathcal{V}}_h^1\text{mbox}foralln\in I\right\},$$

$$-\frac{\partial F}{\partial t}(x,t)+\mathcal{L}F(x,t)=0\text{ in }{\Omega }_T,$$

$$\frac{\partial F(x,t)}{\partial \vec{n}}+\sigma (x,t)F(x,t)=w{\chi }_{\Sigma }\text{ on }\mathcal{S},$$

$$F(x,T)=0\text{ in }\Omega$$

$$\Vert f_{\rho ,\delta ,h,\tau }-f^{†}{\Vert }_{L^2({\Omega }_T)}^2\le C\left(h^3{\rho }^{-1}+{\tau }^2{h}^{-1}{\rho }^{-1}+\delta +\rho +{\delta }^2{\rho }^{-1}\right)$$

$$\Vert f_{\rho ,\delta ,h,\tau }-f^{†}{\Vert }_{L^2({\Omega }_T)}=\mathcal{O}(\sqrt{\delta })$$

$$\displaystyle\mathcal{W}(0,T):=\left\{w\in L^{2}\left(0,T;H^{1}(\Omega)\right)~{}\big{|}~{}\frac{\partial w}{\partial t}\in L^{2}\left(0,T;H^{1}(\Omega)^{*}\right)\right\}.$$

$$\Vert w{\Vert }_{\mathcal{W}(0,T)}:={\left(\Vert w{\Vert }_{L^2\left(0,T;H^1(\Omega )\right)}^2+{\Vert \frac{\partial w}{\partial t}\Vert }_{L^2\left(0,T;H^1(\Omega {)}^{\ast }\right)}^2\right)}^{1/2}.$$

$$a(v,w)$$

$$:={\int }_{\Omega }A(x,t)\nabla v(x)\cdot \nabla w(x)dx+{\int }_{\Omega }b(x,t)v(x)w(x)dx+{\int }_{\partial \Omega }\sigma (x,t)v(x)w(x)dx$$

$${⟨\frac{\partial u}{\partial t},\phi ⟩}_{\left(H^1(\Omega {)}^{\ast },H^1(\Omega )\right)}+a(u,\phi )=(f,\phi {)}_{L^2(\Omega )}+(g,\phi {)}_{L^2(\partial \Omega )}$$

$$\Vert u{\Vert }_{\mathcal{W}(0,T)}\le C\left(\Vert f{\Vert }_{L^2({\Omega }_T)}+\Vert g{\Vert }_{L^2(\mathcal{S})}+\Vert q{\Vert }_{L^2(\Omega )}\right)$$

$$u:L^2({\Omega }_T)\to \mathcal{W}(0,T)$$

$${\Vert Z-z_{\delta }\Vert }_{L^2(\Sigma )}\le \delta ,\delta >0$$

$$C_1\Vert \phi {\Vert }_{H^1(\Omega )}^2\le a(\phi ,\phi )\le C_2\Vert \phi {\Vert }_{H^1(\Omega )}^2$$

$$\frac{\partial v}{\partial t}(x,t)-\sum _{i,j=1}^d\frac{\partial }{\partial x_i}\left(a_{ij}(x,t)\frac{\partial v(x,t)}{\partial x_j}\right)+v(x,t)=e^{-t}f(x,t)\text{ in }{\Omega }_T,$$

$$\frac{\partial v(x,t)}{\partial \vec{n}}=e^{-t}g(x,t)\text{ on }\mathcal{S}\text{mbox}andv(x,0)=q(x)\text{ in }\Omega .$$

$$-\frac{\partial p}{\partial t}(x,t)+\mathcal{L}p(x,t)=0\text{ in }{\Omega }_T,$$

$$\frac{\partial p(x,t)}{\partial \vec{n}}+\sigma (x,t)p(x,t)=\left(u(x,t;f)-z_{\delta }(x,t)\right){\chi }_{\Sigma }\text{ on }\mathcal{S},$$

$$p(x,T)=0\text{ in }\Omega ,$$

$$-{⟨\frac{\partial p}{\partial t},\phi ⟩}_{\left(H^1(\Omega {)}^{\ast },H^1(\Omega )\right)}+a(p,\phi )={\int }_{\Gamma }(u(x,t;f)-z_{\delta }(x,t))\phi (x)dx$$

Full text

{AMS}

35R25, 35R30, 65J20, 65J22

1 Introduction

The problem of identifying a source in a heat transfer or diffusion process has got attention of many researchers during last years. This problem leads to determining a term in the right hand side of parabolic equations from some observations of the solution which is well known to be ill-posed. For surveys on the subject, we refer the reader to the books [5, 21, 28, 29, 36], the recent papers [22, 37] and the references therein.

Although there have been many papers devoted to the source identification problems with observations in the whole domain or at the final moment, those with boundary observations are quite few. Furthermore, the sought term depends either on the spatial variable as in [4, 6, 7, 8, 9, 10, 12, 19, 20, 21, 48, 54, 55], or only on the time variable as in [23]. In this paper, we consider the problem of determining the right hand side depending on both spatial and time variables by a variational method. Indeed, let $\Omega$ be an open bounded connected subset of $\mathbb{R}^{d},d\geq 1$ with boundary $\partial\Omega$ and $T>0$ be a given constant. We investigate the problem of identifying the source term $f=f(x,t)$ in the Robin boundary value problem for the parabolic equation

[TABLE]

from a partial boundary measurement $z_{\delta}:=z_{\delta}(x,t)\in L^{2}(\Sigma)$ of the solution $u(x,t)$ on the surface $\Sigma:=\Gamma\times(0,T)\subset\mathcal{S}$ satisfying

[TABLE]

where $Z:=u_{|\Sigma}$, $\Gamma$ is a relatively open subset of $\partial\Omega$ and the positive constant $\delta$ stands for the measurement error.

In Eq. 1 $\mathcal{L}$ is a time-dependent, second order self-adjoint elliptic operator of the form

[TABLE]

where $A:=\left(a_{ij}\right)_{1\leq i,j\leq d}\in{C(\overline{\Omega_{T}})}^{d\times d}$ is a symmetric diffusion matrix satisfying the uniformly elliptic condition

[TABLE]

for all $\xi=\left(\xi_{i}\right)_{1\leq i\leq d}\in\mathbb{R}^{d}$ with some constant $\underline{a}>0$ and $b(x,t)\in C(\overline{\Omega_{T}})$ is a non-negative function. The vector $\vec{n}:=\vec{n}(x,t)$ is the unit outward normal on $\mathcal{S}$ and

[TABLE]

with $\nabla u(x,t):=\nabla_{x}u(x,t)=\left(\frac{\partial u(x,t)}{\partial x_{1}},\ldots,\frac{\partial u(x,t)}{\partial x_{d}}\right)$. In addition, the functions $q\in H^{1}(\Omega)$, $g\in C(\overline{\mathcal{S}})$ and $\sigma\in C(\overline{\mathcal{S}})$ with $\sigma(x,t)\geq 0$ in $\overline{\mathcal{S}}$ are assumed to be given. The source term $f=f(x,t)$ is sought in the space $L^{2}(\Omega_{T})$.

The contents of this paper are as follows: For any fixed $f\in L^{2}(\Omega_{T})$ let $u=u(f)\in\mathcal{W}(0,T)$ denote the unique weak solution of the system Eq. 1, see Section 2 for the definition of related functional spaces. Adopting the output least squares method combined with the Tikhonov regularization, we consider the (unique) minimizer of the minimization problem

[TABLE]

as a reconstruction, where $\rho\in(0,1)$ is the regularization parameter and $f^{*}$ is an a priori estimate of the exact source which is identified. We mention that $f^{*}$ plays the role of a selection criterion. In practice it is not easy to estimate $f^{*}$, however an a priori knowledge of the identified source improves the quality of reconstructions (cf. [13, 32]).

For discretization we employ the Crank-Nicolson Galerkin method, where the finite dimensional space \mathcal{V}_{h}^{1}:=\big{\{}\varphi_{h}\in C\left(\overline{\Omega}\right)~{}|~{}{\varphi_{h}}_{|T}\in\mathcal{P}_{1}(T)\mbox{~{}for all~{}}T\in\mathcal{T}_{h}\big{\}} of piecewise linear, continuous finite elements is used to discretize the state with respect to the spatial variable. Further, to discretize the state with respect to the time variable, we divide the time interval $(0,T)$ into $M$ equal subintervals and introduce a time step $\tau:=T/M$ together with time levels

[TABLE]

As a result, the state $u(f)$ is then approximated by the finite sequence $(U^{n}_{h,\tau}(f))_{n=0}^{M}$ in which for each $n\in I_{0}$ the element $U^{n}_{h,\tau}(f)\in\mathcal{V}_{h}^{1}.$ With these notions at hand, we examine the discrete regularized problem corresponding to $\left(\mathcal{P}_{\rho,\delta}\right)$ i.e. the following strictly convex minimization problem

[TABLE]

with

[TABLE]

which admits a unique solution $f_{\rho,\delta,h,\tau}$ obeying the relation (see Section 3)

[TABLE]

for any $n\in I:=\{1,\ldots,M\}$, where $(P^{n}_{h,\tau}(f))_{n=0}^{M}$ is the approximation of the adjoint state $p(f)$. Using the variational discretization concept introduced in [24], the minimizer automatically belongs to the finite dimensional space

[TABLE]

thus a discretization of the admissible set $L^{2}(\Omega_{T})$ can be avoided.

As $\delta,h,\tau\to 0$ and with an appropriate a priori regularization parameter choice $\rho=\rho(\delta,h,\tau)\to 0$, we show in Section 4 that the whole sequence \big{(}f_{\rho,\delta,h,\tau}\big{)}_{\rho>0} converges in the $L^{2}(\Omega_{T})$-norm to the unique $f^{*}$-minimum-norm solution $f^{\dagger}$ of the identification problem defined by

[TABLE]

The corresponding state sequence then converges in the $L^{2}(0,T;H^{1}(\Omega))$-norm to the state $u(f^{\dagger})$.

Section 5 is devoted to convergence rates for the discretized problem, where we first show that if $f\in\left\{f\in L^{2}(\Omega_{T})~{}|~{}u(f)_{|\Sigma}=Z\right\}$ and there exists a function $w\in L^{2}(\Sigma)$ such that $f=F(w)+f^{*}$, where $F(w)$ is the unique weak solution of the parabolic equation

[TABLE]

with $\chi_{\Sigma}$ being the characteristic function of $\Sigma\subset\mathcal{S}$, then $f=f^{\dagger}$. Furthermore, if the data appearing in the system Eq. 1 are regular enough the error bound

[TABLE]

is established, that yields the convergence rate

[TABLE]

with the suitable choice of the parameters $h$, $\tau$ and $\rho$.

For the numerical solution of the discrete regularized problem \big{(}\mathcal{P}_{\rho,\delta,h,\tau}\big{)} we in Section 6 utilize a conjugate gradient algorithm. Numerical studies are presented for several cases where the sought source is smooth and discontinuous as well, that illustrates the efficiency of the proposed approach.

To conclude this introduction we wish to mention that to the best of our knowledge, although there have been many papers devoted to source identification problems for parabolic equations, we however have not yet found investigations on the discretization analysis for those with boundary observations — which is more realistic from the practical point of view, a fact that motivated the research presented in the paper.

Regarding the identification problem in elliptic equations utilizing boundary measurements, we here would like to comment briefly some previously published works. In [49, 51] the authors used finite element methods to numerically recover 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 considered in [50]. In [25, 26, 27, 40] some source and coefficient identification problems have been investigated, and in [39] the problem of simultaneously identifying the source term and coefficients from distributed observations. A survey of the optimal control problems for parabolic equations can be found in the book [47], where distributed data is also taken into account.

Throughout the paper we use the standard notion of Sobolev spaces $H^{1}(\Omega)$, $H^{2}(\Omega)$, $W^{k,p}(\Omega)$, etc. from, for example, [1].

2 Problem setting and preliminaries

To formulate the identification problem, we first give some notations [53]. Let $H^{1}(\Omega)^{*}$ be the dual space of $H^{1}(\Omega)$, we use the notation

[TABLE]

It is a Banach space equipped with the norm

[TABLE]

We note that, since $\mathcal{W}(0,T)$ with respect to the norm Eq. 7 is a closed subspace of the reflexive space $L^{2}\left(0,T;H^{1}(\Omega)\right)\times L^{2}\left(0,T;H^{1}(\Omega)^{*}\right)$, it is itself reflexive.

2.1 Direct and inverse problems

For considering Equation 1, we set

[TABLE]

where $t\in(0,T]$ and $v,w\in H^{1}(\Omega)$. Then, for each $f\in L^{2}(\Omega_{T})$ Equation 1 defines a unique weak solution $u:=u(x,t;f):=u(f)$ in the sense that $u(f)\in\mathcal{W}(0,T)$ with $u(x,0)=q(x)$ for a.e. $x\in\Omega$ and the following variational equation is satisfied (cf. [47, 53])

[TABLE]

for all $\varphi\in H^{1}(\Omega)$ and a.e. $t\in(0,T]$. Furthermore, the estimate

[TABLE]

holds, where $C$ is a positive constant independent of $f$, $g$ and $q$. Therefore, we can define the source-to-state operator

[TABLE]

which maps each $f\in L^{2}(\Omega_{T})$ to the unique weak solution $u:=u(f)$ of Equation 1. The inverse problem is stated as follows:

Given the boundary data $Z={u}_{|\Sigma}$ of the exact solution $u$, find an element $f\in L^{2}(\Omega_{T})$ such that $u(f)_{|\Sigma}=Z$.

2.2 Variational method

In practice only the observation $z_{\delta}\in L^{2}(\Sigma)$ of the exact $Z$ with an error level

[TABLE]

is available. Hence, our problem is to reconstruct an element $f\in L^{2}(\Omega_{T})$ in Eq. 1 from noisy observation $z_{\delta}$ of $Z$. For this purpose we use the standard least squares method with Tikhonov regularization, i.e. we consider a minimizer of the minimization problem $\left(\mathcal{P}_{\rho,\delta}\right)$ as a reconstruction.

Remark 2.1**.**

In case $\inf_{(x,t)\in\Omega_{T}}b(x,t)>0$ or $\inf_{(x,t)\in\mathcal{S}}\sigma(x,t)>0$ the expression $a(v,w)$ generates a scalar product on the space $H^{1}(\Omega)$ equivalent to the usual one, i.e. there exist positive constants $C_{1},C_{2}$ such that (cf. [34, 38])

[TABLE]

for all $\varphi\in H^{1}(\Omega)$ and $t\in(0,T]$.

Now we assume that $b=\sigma=0$. A change of the variable $u=e^{t}v$, Equation 1 has the form

[TABLE]

*In the sequel we thus consider the case $\inf_{(x,t)\in\Omega_{T}}b(x,t)>0$ or $\inf_{(x,t)\in\mathcal{S}}\sigma(x,t)>0$ only. All results in present paper are still valid for the case $b=\sigma=0$. *

Remark 2.2**.**

*If the sequence $\left(f_{k}\right)\subset L^{2}(\Omega_{T})$ weakly converges in $L^{2}(\Omega_{T})$ to an element $f$, then the sequence $\left(u(f_{k})\right)$ weakly converges in $\mathcal{W}(0,T)$ to $u(f)$. *

For each $f\in L^{2}(\Omega_{T})$ we now consider the adjoint problem

[TABLE]

where $\chi_{\Sigma}$ is the characteristic function of $\Sigma\subset\mathcal{S}$, i.e. $\chi_{\Sigma}(x,t)=1$ if $(x,t)\in\Sigma$ and equals to zero otherwise. A function $p\in\mathcal{W}(0,T)$ is said to be a weak solution to this problem, if $p(x,T)=0$ for a.e. $x\in\Omega$ and

[TABLE]

for all $\varphi\in H^{1}(\Omega)$ and a.e. $t\in(0,T]$. Since $u\in\mathcal{W}(0,T)$, the boundary value $\left(u(x,t;f)-z_{\delta}(x,t)\right)\chi_{\Sigma}$ belongs to $L^{2}(\Sigma)$ and by changing the time direction we see that Eq. 11 attains a unique weak solution $p(x,t;f):=p(f)\in\mathcal{W}(0,T)$.

We close this section with the following result.

Theorem 2.3**.**

The minimization problem $\left(\mathcal{P}_{\rho,\delta}\right)$ attains a unique minimizer $f_{\rho,\delta}$ which satisfies the equation

[TABLE]

Proof 2.4**.**

Let $(f_{k})$ be a minimizing sequence of the problem $\left(\mathcal{P}_{\rho,\delta}\right)$, i.e.

[TABLE]

*The sequence $(f_{k})$ is then bounded in the $L^{2}(\Omega_{T})$-norm. A subsequence not relabeled and an element $f\in L^{2}(\Omega_{T})$ exist such that $(f_{k})$ weakly converges to $f$ in $L^{2}(\Omega_{T})$. Due to Remark 2.2 and the continuity of the mapping $L^{2}(0,T;H^{1}(\Omega))\to L^{2}(0,T;L^{2}(\Gamma))=L^{2}(\Sigma)$, it follows that $u(f_{k})$ weakly converges to $u(f)$ in $L^{2}(\Sigma)$. And since the $L^{2}$-norm is weakly lower semi-continuous, is so the functional $J_{\rho,\delta}(\cdot)$. Furthermore, it is clear that $J_{\rho,\delta}(f)$ is strictly convex, hence it attains a unique minimizer. Meanwhile, the equation Eq. 13 follows directly from the first order optimality condition, which finishes the proof. *

3 Crank-Nicolson Galerkin discretization

In this section we present the Crank-Nicolson Galerkin method (see, e.g., [46]) to discretize the regularized minimization problem in finite dimensional spaces.

Let $\left(\mathcal{T}_{h}\right)_{0<h<1}$ be a regular family of quasi-uniform triangulations of the domain $\overline{\Omega}$ with the mesh size $h$ (cf. e.g., [3]). We first consider the finite element space $\mathcal{V}_{h}^{1}$ and recall the Clément interpolation operator $\Pi_{h}:L^{1}(\Omega)\rightarrow\mathcal{V}_{h}^{1}$ which satisfies the following properties

[TABLE]

and

[TABLE]

for $0\leq k<l\leq 2$ (see [11, 2, 45]). We also mention that for all $v\in\mathcal{V}^{1}_{h}$, the inverse inequality

[TABLE]

holds true. Further, for all $v\in H^{1}(\Omega)$ and $T\in\mathcal{T}_{h}$ there holds the local estimate ([35])

[TABLE]

To discretize the state functions with respect to the time variable, we consider the partition Eq. 4 of the interval $(0,T)$. For a continuous function $\zeta:[0,T]\to\mathbb{R}$ and $t^{n}\in[0,T]$, we denote by

[TABLE]

Then, we set

[TABLE]

for all $v,w\in H^{1}(\Omega)$.

Linking to the above partition, we introduce the constant piecewise, discontinuous interpolation operator

[TABLE]

which is defined for each $\eta\in L^{2}(0,T)$ by

[TABLE]

By Proposition 9 of [41, pp. 129], we get the limit

[TABLE]

and the estimate (cf. [33, Proposition 5.1])

[TABLE]

in case $\eta\in H^{1}(0,T)$.

For a sequence $(w^{k})\in L^{2}(\Omega)$ we respectively introduce the backward difference quotient and the mean as follows

[TABLE]

We are in the position to present the Crank-Nicolson Galerkin method applied to Eq. 1. Let $q_{h}:=\Pi_{h}q\in\mathcal{V}_{h}^{1}$, find $U^{n}(f):=U^{n}_{h,\tau}(f)\in\mathcal{V}^{1}_{h}$ such that $U^{0}(f)=q_{h}$ and

[TABLE]

for all $\varphi_{h}\in\mathcal{V}^{1}_{h}$, $n\in I$.

Definition 3.1**.**

The mapping $u_{h,\tau}:L^{2}(\Omega_{T})\to\left(\mathcal{V}_{h}^{1}\right)^{M}$ defined for each $f\in L^{2}(\Omega_{T})$ by

[TABLE]

*is called the discrete source-to-state operator. *

The operator $u_{h,\tau}$ is Fréchet differentiable on $L^{2}(\Omega_{T})$. For each $f\in L^{2}(\Omega_{T})$ in the direction $\xi\in L^{2}(\Omega_{T})$ its differential $u_{h,\tau}^{\prime}(f)\xi$ is $\left(\widehat{U}^{1}(\xi),\ldots,\widehat{U}^{M}(\xi)\right)$, where $\widehat{U}^{n}(\xi):=\widehat{U}^{n}_{h,\tau}(\xi)\in\mathcal{V}^{1}_{h}$ is defined by the variational equation

[TABLE]

The adjoint problem Eq. 11 is discretized via the process that for $n\in I_{0}$ the element $P^{n}(f):=P^{n}_{h,\tau}(f)\in\mathcal{V}^{1}_{h}$ satisfies the system

[TABLE]

for all $\varphi_{h}\in\mathcal{V}^{1}_{h}$ with $P^{M}(f)=0$. Note that from $P^{M}(f)=0$ we can compute $P^{M-1}(f)$ due to Eq. 23, and so $P^{M-2}(f)$, …, $P^{1}(f),P^{0}(f)$.

To begin, we recall the discrete Gronwall inequality.

Lemma 3.2**.**

Assume that $(Y_{n})_{n}$,  $(X_{n})_{n}$,  and $(\alpha_{n})_{n}$ are non-negative sequences such that

[TABLE]

Then

[TABLE]

*for all $n\geq 0$. *

Lemma 3.3**.**

(i) There holds the estimate

[TABLE]

(ii) The inequalities

[TABLE]

for all $\theta\in L^{2}\left(0,T;H^{1}(\Omega)\right)$ and

[TABLE]

hold true.

(iii) The limit

[TABLE]

*is satisfied for all $\theta\in L^{2}\left(0,T;H^{1}(\Omega)\right)$. *

Proof 3.4**.**

(i) Taking $\varphi_{h}=2\tau^{2}\partial U^{n}(f)$ in Eq. 21, we have

[TABLE]

By Eq. 10, we thus get

[TABLE]

For an arbitrary $\epsilon>0$, an application of Young’s inequality yields that

[TABLE]

Meanwhile, we have

[TABLE]

We thus arrive at

[TABLE]

as $\epsilon$ is small enough with $C$ independent of $n$. Since $\|U^{0}(f)\|_{H^{1}(\Omega)}=\|\Pi_{h}q\|_{H^{1}(\Omega)}\leq C\|q\|_{H^{1}(\Omega)}$, we deduce from 3.4 that

[TABLE]

for all $n\geq 1$, where $\alpha_{0}=\ldots=\alpha_{n-2}=0$ and $\alpha_{n-1}=C$. Therefore, an application of the discrete Gronwall inequality implies Eq. 24 for all $n\geq 0$.

(ii) By Eq. 21, for all $\theta\in L^{2}\left(0,T;H^{1}(\Omega)\right)$ we rewrite

[TABLE]

We have that

[TABLE]

and

[TABLE]

Using Eq. 24, we further get that

[TABLE]

Next, we deduce from 3.4 that

[TABLE]

Utilizing Eq. 15, we have that

[TABLE]

by 3.4 as $h^{2}\tau^{-1}\leq C$. Therefore, Lemma 3.3 follows from 3.4 and the above estimates. Furthermore, in the same manner we also get Lemma 3.3.

(iii) We consider the piecewise constant function with respect to $t$ defined as follows

[TABLE]

Due to Eq. 24, the sequence $\left(\Phi_{h,\tau}\right)_{h,\tau}$ is bounded in the $L^{2}\left(0,T;H^{1}(\Omega)\right)$-norm. Therefore, there exist a subsequence of it denoted the same symbol and an element $\Phi\in L^{2}\left(0,T;H^{1}(\Omega)\right)$ such that for all $\theta\in L^{2}\left(0,T;H^{1}(\Omega)\right)$

[TABLE]

and $\|\Phi\|_{L^{2}(\Omega_{T})}\leq\liminf_{h,\tau\to 0}\left\|\Phi_{h,\tau}\right\|_{L^{2}(\Omega_{T})}$, since the mapping $L^{2}\left(0,T;H^{1}(\Omega)\right)\hookrightarrow L^{2}(\Omega_{T})$ is continuous. On the other hand, it follows from 3.4 that

[TABLE]

*as $h,\tau\to 0$, which finishes the proof. *

Lemma 3.5**.**

*Assume that the sequence $\left(f_{k}\right)$ weakly converges in $L^{2}(\Omega_{T})$ to an element $f$. Then for any fixed $n\in I_{0}$, the sequence $\left(U^{n}(f_{k})\right)$ converges to $U^{n}(f)$ in the $H^{1}(\Omega)$-norm. *

Proof 3.6**.**

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

Theorem 3.7**.**

The problem $\left(\mathcal{P}_{\rho,\delta,h,\tau}\right)$ attains a unique solution $f:=f_{\rho,\delta,h,\tau}$ satisfying the equation

[TABLE]

*for any $n\in I$. *

Proof 3.8**.**

In virtue of Lemma 3.5, it is straightforward to show the uniqueness existence of a solution $f$ to $\left(\mathcal{P}_{\rho,\delta,h,\tau}\right)$. We now show the equation Eq. 31. We have $J^{\prime}_{\rho,\delta,h,\tau}(f)\xi=0$ for all $\xi\in L^{2}(\Omega_{T})$. By Eq. 22–Eq. 23, we get

[TABLE]

Using the identities

[TABLE]

together with $P^{M}(f)=\widehat{U}^{0}(\xi)=0$, we obtain

[TABLE]

For any fixed $n\in I$ and for all $\varphi\in L^{2}(\Omega\times(t^{n-1},t^{n}])$ we consider $\xi:=\varphi\chi_{\Omega\times(t^{n-1},t^{n}]}\in L^{2}(\Omega_{T})$ and then have $\overline{\xi}^{k}=0$ as $k\neq n$ and

[TABLE]

as well as

[TABLE]

Thus, we arrive at

[TABLE]

*where $\varphi\in L^{2}(\Omega\times(t^{n-1},t^{n}])$ is arbitrary. This implies Eq. 31. The proof is finished. *

Remark 3.9**.**

For any fixed $f\in L^{2}(\Omega_{T})$, denote by

[TABLE]

In view of the identity 3.8, the $L^{2}$-gradient of the cost functional at $f$ is given by

[TABLE]

i.e. the equality

[TABLE]

*holds true for all $\xi\in L^{2}(\Omega_{T})$. *

4 Convergence of finite dimensional approximations

The aim of this section is to show finite dimensional approximations, i.e. solutions of $\left(\mathcal{P}_{\rho,\delta,h,\tau}\right)$, converge to the sought source. To do so, we state some auxiliary results.

Lemma 4.1**.**

(i) For all $\phi\in L^{2}\left(0,T;H^{k}(\Omega)\right)$ with $k\in\{0,1\}$ there hold

[TABLE]

(ii) Assume that $\phi\in L^{\infty}\left(\Omega_{T}\right)$, then

[TABLE]

*for all $u\in L^{2}(\Omega_{T})$ and any bounded sequence $(v_{\tau})_{\tau}$ in $L^{2}(\Omega_{T})$. *

Proof 4.2**.**

(i) The first and second statements of Eq. 35 follow directly from Eq. 14 and Eq. 19, resectively, and the Lebesgue’s dominated convergence theorem. Meanwhile, for a.e in $t\in(0,T)$, since

[TABLE]

the third assertion follows from the first and second ones.

For (ii) we take an arbitrary $\epsilon>0$ and $u^{\epsilon}\in C(\overline{\Omega_{T}})$ such that $\|u-u^{\epsilon}\|_{L^{2}(\Omega_{T})}<\epsilon$. Then we have

[TABLE]

*Sending $\tau$ to zero, we thus have that $\lim_{\tau\to 0}\left|\int_{0}^{T}\int_{\Omega}(\phi_{\tau}-\phi)uv_{\tau}dxdt\right|\leq C\epsilon$ for all $\epsilon>0$ with the constant $C$ independent of $\epsilon$, which yields Eq. 36. The proof is completed. *

Lemma 4.3**.**

Assume that $\lim_{k\to\infty}h_{k}=\lim_{k\to\infty}\tau_{k}=0$ and the sequence $\left(z_{\delta_{k}}\right)\subset L^{2}(\Sigma)$ converge to $z_{\delta}$ in the $L^{2}(\Sigma)$-norm. Let the sequence $(f_{k})\subset L^{2}(\Omega_{T})$ weakly converge in $L^{2}(\Omega_{T})$ to $f$. Then,

[TABLE]

*where $M_{k}=T/\tau_{k}$ and $U^{n}_{h_{k},\tau_{k}}(f_{k})$ defined by Eq. 21. *

Proof 4.4**.**

For convenience of exposition we denote by $U^{n}_{k}:=U^{n}_{h_{k},\tau_{k}}(f_{k})$. Let $\Phi_{k}=\Phi_{k}(x,t)$ be the piecewise linear, continuous interpolation of $\left(U^{n}_{k}\right)_{n=0,\ldots,M_{k}}$ with respect to $t$, i.e.

[TABLE]

with $\partial U^{n}_{k}=\tau_{k}^{-1}(U^{n}_{k}-U^{n-1}_{k})$ and $(x,t)\in\Omega\times(t^{n-1},t^{n}]$, $n=1,\ldots,M_{k}$. We first note that for all $t\in(t^{n-1},t^{n})$

[TABLE]

Thus, $\left(\Phi_{k}\right)\subset H^{1}\left(0,T;\mathcal{V}^{1}_{h}\right)\subset H^{1}\left(0,T;H^{1}(\Omega)\right)\subset\mathcal{W}(0,T)$. Further, the inequalities Eq. 24 and Lemma 3.3 yield that the sequence $\left(\Phi_{k}\right)$ is bounded in the reflexive space $\mathcal{W}(0,T)$. There exists a subsequence of $\left(\Phi_{k}\right)$ denoted again by $\left(\Phi_{k}\right)$ and an element $u\in\mathcal{W}(0,T)$ such that $\left(\Phi_{k}\right)$ weakly converges in $\mathcal{W}(0,T)$ to $u$.

We show that $u=u(f)$. In fact, for all $\phi\in L^{2}\left(0,T,H^{1}(\Omega)\right)$ we have

[TABLE]

For all $t\in(t^{n-1},t^{n}]$, we have

[TABLE]

and, by Lemma 3.3 and Eq. 35,

[TABLE]

as $k\to\infty$. Using the estimates Eq. 24 and Lemma 3.3 as well as the equalities Eq. 35, Eq. 36, we decompose the remainder in the right hand side of 4.4 as follows

[TABLE]

by Eq. 38, where

[TABLE]

We remark that due to the continuity of data and Eq. 35, the relation $\lim_{k\to\infty}\mathcal{R}_{k}=0$ holds true. Therefore, we obtain from 4.4–4.4 that

[TABLE]

here we used Eq. 21. Further, with standard arguments it can be shown that $u(x,0)=q(x)$. Thus, we get that the sequence $\left(\Phi_{k}\right)$ weakly converges in $\mathcal{W}(0,T)$ to $u(f)$.

Next, for each $\theta\in C^{1}\left([0,T];H^{1}(\Omega)\right)$ we will show that

[TABLE]

which then holds true for each $\theta\in L^{2}\left(0,T;H^{1}(\Omega)\right)$, by the density argumentation. We have that

[TABLE]

and so

[TABLE]

by Eq. 27 and Eq. 24.

*Finally, using the continuity of the mapping $L^{2}\left(0,T;H^{1}(\Omega)\right)\to L^{2}\left(0,T;L^{2}(\Gamma)\right)$ and the fact $\|z_{\delta_{k}}-z_{\delta}\|_{L^{2}(\Sigma)}\to 0$ as $k\to\infty$, we arrive at Eq. 37, which finishes the proof. *

Next we introduce the notion of the $f^{*}$-minimum-norm solution of the identification problem.

Lemma 4.5**.**

The problem

[TABLE]

attains a unique solution, which is called the $f^{*}$-minimum-norm solution of the identification problem, where

[TABLE]

Proof 4.6**.**

*Due to Remark 2.2, $\mathcal{I}\left(Z\right)\neq\emptyset$ is a close subset of $L^{2}(\Omega_{T})$. Furthermore, it is a convex set. Therefore, the minimization problem has a unique solution, which finishes the proof. *

We now show the main result of this section on the convergence of finite dimensional approximations $f_{\rho,\delta,h,\tau}$ of $\left(\mathcal{P}_{\rho,\delta,h,\tau}\right)$ to the $f^{*}$-minimizing-norm solution of the idendification problem $\left(\mathcal{IP}\right)$. For any fixed $f\in L^{2}(\Omega_{T})$, let $u(f)$ and $U^{n}(f)$ define by Eq. 8 and Eq. 21, respectively. We recall the convergence of the Crank-Nicolson Galerkin method for linear parabolic problems

[TABLE]

Theorem 4.7**.**

Let $f^{\dagger}$ be the unique $f^{*}$-minimum-norm solution of the identification problem $\left(\mathcal{IP}\right)$. Let $(h_{k}),(\tau_{k})$, $(\delta_{k})$ and $(\rho_{k})$ be any positive sequences such that

[TABLE]

as $k\to\infty$. Furthermore, assume that $\left(z_{\delta_{k}}\right)\subset L^{2}(\Sigma)$ is a sequence satisfying

[TABLE]

and $f_{k}$ denotes the unique minimizer of $(\mathcal{P}_{\rho_{k},\delta_{k},h_{k},\tau_{k}})$ for each $k\in\mathbb{N}$. Then:

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

(ii) The following equality holds true

[TABLE]

Proof 4.8**.**

For $n=1,\ldots,M_{k}$ we write $U^{n}_{k}$ and ${U^{n}_{k}}^{\dagger}$ instead $U^{n}_{h_{k},\tau_{k}}(f_{k})$ and $U^{n}_{h_{k},\tau_{k}}(f^{\dagger})$ for short, respectively. Since $f_{k}$ is the solution of $(\mathcal{P}_{\rho_{k},\delta_{k},h_{k},\tau_{k}})$, we have

[TABLE]

We get

[TABLE]

It follows from 4.8 and 4.8 that

[TABLE]

Therefore, by Eq. 46, we have

[TABLE]

and

[TABLE]

Applying Lemma 4.3, we deduce from the boundedness of $(f_{k})$ due to Eq. 51 that there are a subsequence of it denoted by the same symbol and an element $\widehat{f}\in L^{2}(\Omega)$ such that

[TABLE]

We thus obtain from Eq. 50 and Eq. 52 that

[TABLE]

Further, combining Eq. 51 with Eq. 52, we also obtain

[TABLE]

and so, by the uniqueness of the $f^{*}$-minimum-norm solution $f^{\dagger}$,

[TABLE]

Next, by Eq. 21, for all $\varphi_{h_{k}}\in\mathcal{V}^{1}_{h_{k}}$ and $n\in I$ we have

[TABLE]

Denoting by $e^{n}_{k}:=U^{n}_{k}-{U^{n}_{k}}^{\dagger}$ and taking $\varphi_{h_{k}}=2\tau^{2}_{k}\partial e^{n}_{k}$ in the above equation, we can estimate the left hand side of Eq. 53 by

[TABLE]

and the right hand side by

[TABLE]

for any $\epsilon>0$. We then get from Eq. 53–4.8

[TABLE]

that yields

[TABLE]

where $\alpha_{0}=\ldots=\alpha_{n-2}=0$ and $\alpha_{n-1}=C$. Since $e^{0}_{k}=0$, applying Gronwall’s inequality, we obtain

[TABLE]

and so that

[TABLE]

as $k\to\infty$. Therefore, we in view of Eq. 45 conclude that

[TABLE]

*as $k\to\infty$, which finishes the proof. *

5 Convergence rates

To obtain convergence rates for the Tikhonov regularization, some additional assumptions are required. We here assume that the data appearing in the system Eq. 1 are regular enough such that the following error bound of the Crank-Nicolson Galerkin method for linear parabolic problems is fulfilled, see, e.g., [46].

Lemma 5.1**.**

Let $u(f^{\dagger})$ and $U^{n}(f^{\dagger})$ be the solutions of Eq. 8 and Eq. 21, respectively. Then, the estimate

[TABLE]

*holds true for all $n\in I$, where $s\in\{0,1\}$. *

We note that, due to the regularity theory for the parabolic equations (see, e.g., [14, Theorem 6, pp. 365]), there exist the high order derivatives $\Delta u_{tt}(f^{\dagger})$ and $u_{ttt}(f^{\dagger})$ if $f^{\dagger}_{tt}\in L^{2}(0,T;L^{2}(\Omega))$ and $f^{\dagger}_{ttt}\in L^{2}(0,T;H^{-2}(\Omega))$, respectively.

To obtain the convergence rates of the regularized solutions to the identification, some smooth assumptions on the sought source $f^{\dagger}$ and the exact data $u(f^{\dagger})$ should be addressed. We first assume they are such that the estimate Lemma 5.1 is fulfilled, and have the following result.

Theorem 5.2**.**

Assume that $\tilde{f}\in\mathcal{I}\left(Z\right)$ and there exists a function $w\in L^{2}(\Sigma)$ such that

[TABLE]

where $F(w)$ is the weak solution of the equation Eq. 6. Then:

(i) $\tilde{f}=f^{\dagger}$, i.e. it is the unique solution of the identification problem $\left(\mathcal{IP}\right)$.

(ii) The estimate

[TABLE]

holds, where $f$ denotes the unique minimizer of \big{(}\mathcal{P}_{\rho,\delta,h,\tau}\big{)}.

(iii) With the noise level $\delta$, the choice $h\sim\delta$, $\tau\sim\delta^{3/2}$ and $\rho\sim\delta$ leads to the convergence rate

[TABLE]

Assume that the coefficients appearing in the elliptic differentiable operator $\mathcal{L}$ and the domain $\Omega$ are smooth enough. Then the solution $F(w)$ to Eq. 6 has the regularity property $F(w)\in H^{k}(0,T;L^{2}(\Omega))$ with $\frac{\partial^{k+1}F(w)}{\partial t^{k+1}}\in L^{2}(0,T;H^{-1}(\Omega))$, $k\geq 0$, provided $w$ has the similar regularity (cf. [14, 53]). Therefore, with an a priori smooth estimate $f^{*}$ the identification $f^{\dagger}=F(w)+f^{*}$ and the exact $u(f^{\dagger})$ satisfy the inequality Lemma 5.1. We also wish to mention that the adjoint approach in the present paper yields the convergence rate $\mathcal{O}(\sqrt{\delta})$ for the Tikhonov regularization under the source condition Eq. 57. Other source conditions may provide different convergence rates or even the optimal-order one $\mathcal{O}(\delta^{2/3})$, which are unfortunately still open to us.

Proof 5.3**.**

(i) For all $\theta\in\mathcal{I}\left(Z\right)$ we rewrite

[TABLE]

and so need to show that $(\tilde{f}-f^{*},\theta-\tilde{f})_{L^{2}(\Omega_{T})}=0$. In fact, by Eq. 8, we have

[TABLE]

Since $F(x,T;w)=u(x,0;\theta)-u(x,0;\tilde{f})=0$, we thus have that

[TABLE]

by $\tilde{f},\theta\in\mathcal{I}\left(Z\right)$.

(ii) By the optimality of $f$, we get that

[TABLE]

with

[TABLE]

Furthermore, by Eq. 17, we get that

[TABLE]

Due to Lemma 5.1, the first term in the right hand side of 5.3 is bounded by

[TABLE]

where we used the estimate

[TABLE]

Likewise, we have that

[TABLE]

We deduce from 5.3–5.3 that

[TABLE]

Since $f^{\dagger}-f^{*}=F(w)$, we have

[TABLE]

We bound for each term in the right hand side of the above equation. First, we get

[TABLE]

and

[TABLE]

Further, as the above estimates 5.3–5.3, we get

[TABLE]

It follows from 5.3–Eq. 67 that

[TABLE]

We therefore conclude from 5.3 and 5.3

[TABLE]

*which completes the proof. *

6 Numerical examples

In this section we employ the conjugate gradient (CG) method (cf. [31]) for reaching the unique minimizer of the finite dimensional minimization problem $(\mathcal{P}_{\rho,\delta,h,\tau})$. Starting with an initial guess $f_{0}$, the iterative process is given by

[TABLE]

where $d_{k}$ is the update direction and $\alpha_{k}$ is the update step size until it meets a stopping criterion of the form

[TABLE]

The update direction is computed as

[TABLE]

where coefficient $\beta_{k}$ is computed using the Polak-Ribière formula

[TABLE]

To compute $\alpha_{k}$, we consider the quadratic minimization problem

[TABLE]

It is shown that the unique minimizer of Eq. 73 is given by

[TABLE]

where for each $n\in I$, $U^{n}_{0}(d_{k})\in\mathcal{V}^{1}_{h}$ is the unique solution of the equation

[TABLE]

in which $\left(U^{n}(f_{k})\right)_{n=1}^{M}$ is defined by the system Eq. 21. Practical steps can be summarized in Algorithm 1.

For numerical implementations below, a common setting is used, unless stated otherwise. We consider a spatial two-dimensional case of Eq. 1, where $\Omega_{T}=(-1,1)\times(-1,1)\times(0,1]$ and

[TABLE]

We denote by $f_{mean}$ the mean of the exact source $f$. For initial iterate and predicted source, we respectively choose $f_{0}=0$ and $f^{*}=f+0.2(f-f_{mean})$. To mimic the noise measurement, we set $z_{\delta}=u(f)_{|\Sigma}+c_{w}z_{rand}$ where $z_{rand}$ is produced by MATLAB command ‘rand’ and $c_{w}$ is chosen such that the measurement error $\|u(f)-z_{\delta}\|_{L^{2}(\Sigma)}$ is equal to the prescribed noise level $\delta$.

Except for the first example where we choose different mesh sizes to show the experimental order of convergence (EOC), we will take $h=0.05$ and generate the mesh using MATLAB function ’generateMesh’. Choice for other constant parameters conceptually follows from Theorem 4.7

[TABLE]

6.1 Time dependent source $f(t)$

In this subsection, we specify the right hand side $f(t)$ independent of the spatial variables, the initial value $u(\cdot,\cdot,0)=0.4$, and boundary condition $g=0.4$. Then the state is numerically computed by solving the corresponding forward problem. We will refer to this source function and the resulting numerical state as the exact source and the exact state, respectively. To generate the noise data $z_{\delta}$, we simulate the forward problem with the aforementioned data but on a finer grid with $h_{fine}=0.025$ to prevent the so called inverse crime (cf. [52]). The derived solution is then interpolated on the grid with mesh size $h$ and perturbed. The observation boundary part is $\Sigma=[-1,1]\times\{-1\}\times(0,1]$ and the result is checked at $P_{1}$ and $P_{2}$, defined as the closest nodal point to $(-0.1,-0.5)$ and $(0.5,0.6)$, respectively. In Figure 1 one can see that at both $P_{1}$ and $P_{2}$, the recovered source matches the exact source quite well in various cases.

We would like to mention that for this special case, the considered inverse problem has a unique solution (cf. [18, 44]). That is, wherever we choose an a priori estimate $f^{*}$, the algorithm should always yield a good approximation of the exact source. In practice, a poorly predicted source may however reduce the quality of the result. In Figure 2 we perform the same test as in Figure 1 but with $f^{*}=0$ for all cases. One can see that the obtaining numerical solutions are not so good, compared with the previous simulations (cf. Figure 1).

6.2 Space-dependent source $f(x,y)$

This implementation is constructed similar to that of Section 6.1, where the source function is independent of the time variable and chosen as

[TABLE]

for all $t$. We note that if the source term depends on the spatial variables, certain space observations are required to guarantee the uniqueness of the identification problem, for example, the additional final data measurement. For this subject we mention to [15, 16, 17, 30, 43] and references given there for detailed discussions.

We assume that the measurement $z_{\delta}$ is available on the whole boundary $\Sigma=\partial\Omega\times[0,T]$. In Figure 3 we compare the recovered state/source with exact state/source with respect to the spatial variables at the time level $t=0.5$.

We also present the test on different mesh sizes and the numerical result is shown in Table 1 with its corresponding EOC being given in Table 2.

The convergence history given in Table 1 and Table 2 shows that the algorithm performs well for our identification problem.

6.3 General source $f(x,y,t)$

In the example, we consider a general case, where $f$ depends on both space and time variables but with a somewhat simplified setting. We set $A=I_{2}$ the $2\times 2$ identity matrix, $b=0,\sigma=0$ and $\Sigma=\partial\Omega\times(0,1]$. Moreover, we choose the exact state as $u(x,y,t)=t(x^{2}-1)^{2}(y^{2}-1)^{2}$. With this setting, it is straightforward to verify that the initial condition $q$, boundary condition $g$ and exact data $Z$ are all zeros. The exact source is then

[TABLE]

We in Figure 4 compare $f_{1}$ and $f_{rec.}$ at nodal point $P_{1}$, while in Figure 5(a) and Figure 5(b) respectively perform them at $t=0.5$. The numerical result shows that our reconstruction method produces a good approximation of the sought source in the general case.

We would like to discuss the role of the predicted source $f^{*}$ in the numerical solution to the identification problem. In principle, the CG method converges to the $f^{*}$-minimum-norm solution $f^{\dagger}$. Therefore, for those inverse problems having more than one solution, the obtaining numerical solution approximates the one that is nearest to $f^{*}$. For the case considered here, we can easily observe that it accepts $f_{2}=0$ (corresponding the state $u_{2}=0$) as another solution. With the setting $f^{*}=0.5$ to imitate the situation where we have no information about the sought source, the algorithm converges to an approximation of $f_{2}$ as shown in Figure 5(c).

6.4 Identified source satisfying the condition Eq. 57

The last example aims at illustrating the convergence rate of the regularized approximations to the identification given in Theorem 5.2. To do so, the setting is similar to the previous subsections with the observations taking on the whole boundary instead. For generating the exact source, we start with the constant function $w=0.2$ on $\Sigma$ and then solve (6) for the solution $F(w)$. As $w$ and the coefficients given in Eq. 75 are constant functions, it deduces that $F(w)\in C^{\infty}(\overline{\Omega}\times[0,T])$ (cf. [14, 53]). Let $f^{*}=(x^{2}+y)t$, we then have $f^{\dagger}=F(w)+f^{*}\in C^{\infty}(\overline{\Omega}\times[0,T])$.

In Table 3 and Table 4, we show the errors corresponding to different hierarchical mesh sizes and the resulting EOC, respectively, where all errors and noisy levels get together smaller. Figure 6 shows the comparison between the exact source and recovered one which match each other so well, as expected from our convergence result.

To close this section we wish to discuss about the Gibbs phenomenon that possibly appears when numerically recovering of discontinuous functions concerning in this section, say Figures 1(c), 2(c), and 3(d). We observe that Gibbs phenomenon slightly happens in Figures 1(c) and 3(d), but it seems to be not in 2(c). These facts might be explained as follows. First, the use of differentiable regularization terms, e.g., the quadratic stabilizing penalty term, may smoothen the recovered solutions. We mention that to reconstruct such discontinuous functions one usually employs the total variation regularization, which was originally introduced in image denoising [42]. And second, by setting, in Figure 2(c) the a priori estimate $f^{*}=0$, while in the other cases, it is discontinuous. Also, obtaining numerical solutions performing in Figure 2(b) seem to be smooth, where the a priori estimate is differentiable; meanwhile we in 1(b) utilize a non-differentiable one.

7 Conclusion

In this paper we investigate the inverse problem of identifying the source function in the parabolic equation

[TABLE]

from a partial boundary measurement $z_{\delta}\in L^{2}(\Sigma)$ of the solution $u(x,t)$ on the surface $\Sigma:=\Gamma\times(0,T)\subset\mathcal{S}$, where $\mathcal{L}$ is a time-dependent, second order self-adjoint elliptic operator, $\Gamma$ is a relatively open subset of $\partial\Omega$ and $\delta>0$ is the error level of the observation.

The Crank-Nicolson Galerkin method is employed to fully discretize the parabolic equation. As a result, the state $u(f)$ is then approximated by the finite sequence $(U^{n}_{h,\tau}(f))_{n=0}^{M}$ in which for each $n\in I_{0}$ the element $U^{n}_{h,\tau}(f)\in\mathcal{V}_{h}^{1}$ — the space of piecewise linear, continuous finite elements, where $h$ and $\tau$ is respectively the mesh size of the space discretization and the time step.

The least squares method combining with the quadratic stabilizing penalty term is utilized to tackle the identification problem, we then consider the unique minimizer $f_{\rho,\delta,h,\tau}$ of the minimization problem

[TABLE]

as a reconstruction, where $\rho>0$ is the regularization parameter and $f^{*}$ is an a priori estimate of the identified source.

We show that with $\delta,h,\tau$ approaching zero and an appropriate a priori regularization parameter choice $\rho=\rho(\delta,h,\tau)$ the whole sequence \big{(}f_{\rho,\delta,h,\tau}\big{)}_{\rho>0} converges in the $L^{2}(\Omega_{T})$-norm to the unique $f^{*}$-minimum-norm solution $f^{\dagger}$ of the identification problem as $\rho$ tends to zero. The corresponding state sequence then converges in the $L^{2}(0,T;H^{1}(\Omega))$-norm to the state $u(f^{\dagger})$. Furthermore, the convergence rate

[TABLE]

is established for an additional suitable source condition and an appropriate choice of the parameters $h$, $\tau$ and $\rho$ coupling with $\delta$. The numerical experiments are presented to illustrate the efficiency of the theoretical findings.

Acknowledgments

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

T.N.T. Quyen gratefully acknowledges support of the University of Goettingen, Germany. N.T. Son is supported in part by the Vietnam National Foundation for Science and Technology Development (NAFOSTED), grant 101.01-2017.319. The paper is completed when he is with the Université catholique de Louvain under the support of EOS Project no. 30468160 funded by FNRS and FWO.