# Numerical identification of initial temperatures in heat equation with
dynamic boundary conditions

S. E. Chorfi  ,  G. El Guermai  ,  L. Maniar  and  W. Zouhair  S. E. Chorfi,
G. El Guermai, L. Maniar and W. Zouhair, Cadi Ayyad University, Faculty of
Sciences Semlalia, LMDP, UMMISCO (IRD-UPMC), B.P. 2390, Marrakesh, Morocco
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

We investigate the inverse problem of numerically identifying unknown initial
temperatures in a heat equation with dynamic boundary conditions whenever some
overdetermination data is provided after a final time. This is a backward
parabolic problem which is severely ill-posed. As a first step, the problem is
reformulated as an optimization problem with an associated cost functional.
Using the weak solution approach, an explicit formula for the Fréchet gradient
of the cost functional is derived from the corresponding sensitivity and
adjoint problems. Then the Lipschitz continuity of the gradient is proved.
Next, further spectral properties of the input-output operator are
established. Finally, the numerical results for noisy measured data are
performed using the regularization framework and the conjugate gradient
method. We consider both one- and two-dimensional numerical experiments using
finite difference discretization to illustrate the efficiency of the designed
algorithm. Aside from dealing with a time derivative on the boundary, the
presence of a boundary diffusion makes the analysis more complicated. This
issue is handled in the 2-D case by considering the polar coordinate system.
The presented method implies fast numerical results.

###### Key words and phrases:

inverse problem, backward parabolic problem, dynamic boundary condition,
adjoint problem, conjugate gradient.

###### 2020 Mathematics Subject Classification:

Primary: 35R30, 49N45; Secondary: 35K05, 47A05.

##  1\. Introduction

Inverse problems arise in various fields of science and engineering, including
the modeling of physical phenomena such as heat transfer problems and
diffusion processes, among others. Most of these problems are modeled by PDEs
where some parameters are unknown and need to be identified from some measured
data. This is exactly where the inverse problems theory comes in. More
precisely, inverse problems theory consists of identifying some input
parameters in a given system from a partial measurement on the solution, e.g.,
forcing terms [14, 13, 17], thermal conductivity and radiative coefficient
[15], potential, damping coefficient and source terms [5, 6, 8], or initial
temperature [21]. Up to our knowledge, the most practical and realistic
measurements considered in literature use the state at a fixed final time.

In this paper, we focus on identifying unknown initial temperatures in a heat
equation with dynamic boundary conditions, from a noisy measurement of the
final temperature. More precisely, let $T>0$ be a fixed final time and let
$\Omega\subset\mathbb{R}^{N}$ be a bounded domain ($N\geq 1$ is an integer)
with boundary $\Gamma=\partial\Omega$ of class $C^{2}$. We denote by
$\Omega_{T}=(0,T)\times\Omega$ and $\Gamma_{T}=(0,T)\times\Gamma$. We consider
the following heat equation with dynamic boundary conditions of surface
diffusion type

$\displaystyle\begin{aligned} &\partial_{t}y-d\Delta
y+a(x)y=0,&&\qquad\text{in }\Omega_{T},\\\
&\partial_{t}y_{\Gamma}-\gamma\Delta_{\Gamma}y_{\Gamma}+d\partial_{\nu}y+b(x)y_{\Gamma}=0,&&\qquad\text{on
}\Gamma_{T},\\\ &y_{\Gamma}(t,x)=y_{|\Gamma}(t,x),&&\qquad\text{on
}\Gamma_{T},\\\
&(y,y_{\Gamma})\rvert_{t=0}=(f,g),&&\qquad\Omega\times\Gamma,\end{aligned}$
(1.1)

where $\mathcal{G}:=(f,g)\in L^{2}(\Omega)\times L^{2}(\Gamma)$ is the unknown
initial temperature. The diffusive constants are such that $d,\gamma>0$ if
$N\geq 2$, with the convention $\gamma=0$ if $N=1$. The radiative potentials
are assumed to be bounded $a\in L^{\infty}(\Omega)$, $b\in
L^{\infty}(\Gamma)$, and $D=\max(\|a\|_{\infty},\|b\|_{\infty})$. The Laplace
operator is denoted by $\Delta=\Delta_{x}$. The trace of $y$ is
$y_{\mid\Gamma}$, and the normal derivative is $\partial_{\nu}y:=(\nabla
y\cdot\nu)_{|\Gamma}$, where $\nu(x)$ is the unit outward normal at
$x\in\Gamma$. Let $\mathrm{g}$ be the standard Riemannian metric on $\Gamma$
induced by $\mathbb{R}^{N}$. We fix a coordinate system $x=(x^{j})$ and we
denote by $\left(\partial_{x^{j}}\right)$ the corresponding tangent vector
field. The Laplace-Beltrami operator $\Delta_{\Gamma}$ is given locally by the
formula

$\Delta_{\Gamma}=\frac{1}{\sqrt{|\mathrm{g}|}}\sum_{i,j=1}^{N-1}\partial_{x^{i}}\left(\sqrt{|\mathrm{g}|}\,\mathrm{g}^{ij}\partial_{x^{j}}\right),$
(1.2)

where $\mathrm{g}=\left(\mathrm{g}_{ij}\right)$ still denotes the tensor
matrix of $\mathrm{g}$, $\mathrm{g}^{-1}=\left(\mathrm{g}^{ij}\right)$ its
inverse with determinant $|\mathrm{g}|=\det\left(\mathrm{g}_{ij}\right)$.

We refer to [22, 20, 11] for a detailed exposition of dynamic boundary
conditions. We also refer to the seminal paper [12], where the physical
derivation and interpretation of dynamic boundary conditions
$\eqref{eq1to4}_{2}$ are discussed. Throughout the paper, we will often use
the surface divergence formula

$\int_{\Gamma}\Delta_{\Gamma}y\,z\,\mathrm{d}S=-\int_{\Gamma}\langle\nabla_{\Gamma}y,\nabla_{\Gamma}z\rangle_{\Gamma}\,\mathrm{d}S,\quad
y\in H^{2}(\Gamma),\,z\in H^{1}(\Gamma),$ (1.3)

where $\nabla_{\Gamma}y$ is the tangential gradient given by
$\nabla_{\Gamma}y=\nabla y-\left(\partial_{\nu}y\right)\nu$, and
$\langle\cdot,\cdot\rangle_{\Gamma}$ is the Riemannian inner product of
tangential vectors on $\Gamma$.

The object of this paper is to numerically identify some unknown initial
temperatures from the final ones. This is the so-called inverse initial data
problem (or the backward parabolic problem) known as a typical problem which
is severely ill-posed in the sense of Hadamard, that is, any arbitrarily small
change in the initial datum lead to drastically huge change in the output
solution [18]. Physically, this can be seen as a consequence of the time
irreversibility and the smoothing effect of the heat process.

Many contributions have been devoted to the enrichment of this field. For
instance, the authors of [27] present a method based on reproducing kernels
applied to a two-dimensional parabolic inverse source problem. The paper [28]
investigate the simultaneous determination of an unknown initial temperature
and a heat source from a given observation at a fixed internal location in
heat equation with Dirichlet boundary conditions. The authors of [23] develop
some numerical schemes based on positivity-preserving Padé approximations to
solve an inverse initial data problem for heat equation with Dirichlet
boundary conditions. In [26], the authors have been interested in recovering
the initial temperature from the measurement on a part of the boundary. The
paper [24] has considered a similar problem with a regularizing parameter that
approximates and regularizes the heat conduction model. More recently, the
authors of [4] have analyzed numerically an inverse source problem in a
degenerate parabolic equation. They have adapted an optimal control approach
using finite element discretization and the augmented Lagrangian method. For
more recent developments within this framework, we refer to the excellent book
[16].

The approach presented in this work can be seen as a continuation of the
approach given in [3] for the identification of source terms in heat equation
with dynamic boundary conditions from the final time measured data. Compared
to the aforementioned inverse source problem, the ill-posedness of the inverse
initial temperatures problem is more severe. In view of the stability results
in [1, 2], one can expect a weak convergence rate and slow numerical results
compared to the fast numerical results for the inverse source problem obtained
in [3] by the Landweber iteration. In fact, in [2], the authors have obtained
a Lipschitz stability estimate for the source terms, while in [1] only a
logarithmic stability estimate is proved for the initial temperatures. It
should be pointed out that such conditional stability results serve to
suitably select the regularizing parameter and determine the convergence rate
of the regularized solution to an exact solution [7]. As a remedy to this
issue, we propose a modified algorithm based on the conjugate gradient (CG)
method to speed up the convergence process. We refer to [9] for a detailed
exposition of such algorithms.

To the best of our knowledge, the problem considered in this paper has only
been shown to be solvable in a particular form of static boundary conditions
(Dirichlet, Neumann or Robin); see, for instance, [24, 26, 23]. However, for
more general boundary conditions as dynamic boundary conditions, the problem
has not been studied, and the numerical reconstruction needs to be
investigated. Therefore, the main objective of the present work is to extend
the weak solution strategy for the determination of initial temperatures to
undertake the numerical aspect of such a problem.

The paper is structured as follows: in Section 2, we recall some well-
posedness and regularity results related to system (1.1) that are needed in
the sequel. In Section 3, we study the minimization problem associated to the
cost functional. Next, we prove an explicit formula for the gradient of the
cost functional. This is done by introducing suitable adjoint and sensitivity
systems. Then the Lipschitz continuity of the gradient is obtained. In Section
5, the gradient formula is algorithmically implemented via the conjugate
gradient method for the numerical reconstruction of initial temperatures in
one- and two-dimensional systems. Finally, Section 6 is devoted to some
conclusions.

##  2\. Wellposedness of the Cauchy problem

We start by recalling some relevant results concerning the wellposedness and
the regularity of the solution of system (1.1). For convenience, we refer to
[22] for more details.

We consider the Lebesgue measure $\mathrm{d}x$ on $\Omega$ and the surface
measure $\mathrm{d}S$ on $\Gamma$. Let us introduce the functional spaces

$\mathbb{L}^{2}:=L^{2}(\Omega,\mathrm{d}x)\times
L^{2}(\Gamma,\mathrm{d}S)\qquad\text{and}\qquad\mathbb{L}_{T}^{2}:=L^{2}(\Omega_{T})\times
L^{2}(\Gamma_{T}).$

Note that $\mathbb{L}^{2}$ and $\mathbb{L}_{T}^{2}$ are Hilbert spaces
respectively endowed with the scalar products defined by

$\displaystyle\langle(y,y_{\Gamma}),(z,z_{\Gamma})\rangle$
$\displaystyle=\langle y,z\rangle_{L^{2}(\Omega)}+\langle
y_{\Gamma},z_{\Gamma}\rangle_{L^{2}(\Gamma)},$
$\displaystyle\langle(y,y_{\Gamma}),(z,z_{\Gamma})\rangle_{\mathbb{L}_{T}^{2}}$
$\displaystyle=\langle y,z\rangle_{L^{2}(\Omega_{T})}+\langle
y_{\Gamma},z_{\Gamma}\rangle_{L^{2}(\Gamma_{T})}.$

We further introduce the functional space

$\mathbb{H}^{k}:=\left\\{(y,y_{\Gamma})\in H^{k}(\Omega)\times
H^{k}(\Gamma)\colon y_{|\Gamma}=y_{\Gamma}\right\\}\text{ for }k=1,2,$

endowed with the standard product norm.

The system (1.1) can be written as an abstract Cauchy problem as follows

$\text{(ACP) }\;\begin{cases}\partial_{t}Y=\mathcal{A}Y,\quad 0<t\leq T,\\\
Y(0)=\mathcal{G}:=(f,g),\end{cases}$

where $Y:=(y,y_{\Gamma})$, and $\mathcal{A}\colon
D(\mathcal{A})\subset\mathbb{L}^{2}\longrightarrow\mathbb{L}^{2}$ is a linear
operator defined by

$\mathcal{A}=\begin{pmatrix}d\Delta-a&0\\\
-d\partial_{\nu}&\gamma\Delta_{\Gamma}-b\end{pmatrix},\qquad\qquad
D(\mathcal{A})=\mathbb{H}^{2}.$

The operator $\mathcal{A}$ is a bounded perturbation of the linear operator
$\mathcal{A}_{0}\colon
D(\mathcal{A}_{0})\subset\mathbb{L}^{2}\longrightarrow\mathbb{L}^{2}$ defined
by

$\mathcal{A}_{0}=\begin{pmatrix}d\Delta&0\\\
-d\partial_{\nu}&\gamma\Delta_{\Gamma}\end{pmatrix},\qquad\qquad
D(\mathcal{A}_{0})=\mathbb{H}^{2}.$

It has been shown in [22] that the operator $\mathcal{A}_{0}$ is the
infinitesimal generator of an analytic $C_{0}$-semigroup on $\mathbb{L}^{2}$.
Thus, the system (1.1) is well-posed and enjoys the maximal regularity
property by a standard perturbation argument.

The following spectral lemma will be useful in the sequel.

######  Lemma 2.1.

There exists an increasing sequence of eigenvalues $(\lambda_{n})_{n\geq 1}$
that constitutes the spectrum of $-\mathcal{A}_{0}$ such that
$0=\lambda_{1}<\lambda_{2}<\cdots$, and $\lambda_{n}\to\infty$ as
$n\to\infty$. Moreover, the sequence of the corresponding eigenfunctions
$(\phi_{n})_{n\geq 1}$ forms an orthonormal basis of $\mathbb{L}^{2}$.

###### Proof.

By [22, Proposition 2.1], the operator $-\mathcal{A}_{0}$ is self-adjoint and
nonnegative. Moreover, since $\mathbb{H}^{2}\hookrightarrow
H^{2}(\Omega)\times H^{2}(\Gamma)$ is continuous, then
$\mathbb{H}^{2}\hookrightarrow\mathbb{L}^{2}$ is compact. By [10, Proposition
4.25], we deduce that $-\mathcal{A}_{0}$ has a compact resolvent. Now, the
existence of $(\lambda_{n})_{n\geq 1}$ and $(\phi_{n})_{n\geq 1}$ follows from
the standard spectral theory. ∎

##  3\. Gradient formula and Lipschitz continuity of the gradient

In this section, we consider the following inverse initial data problem.

Backward Parabolic Problem (BPP). The initial temperature $\mathcal{G}:=(f,g)$
in the system (1.1) is unknown and to be restored from the measured final
temperature

$Y_{T}:=\left(y(T,\cdot),y_{\Gamma}(T,\cdot)\right)\in\mathbb{L}^{2}.$

The measured data is assumed to contain a random noise.

Let $Y(t,\cdot,\mathcal{G})$ be the mild solution of (1.1) corresponding to
the initial datum $\mathcal{G}=(f,g)$. We introduce the input-output operator
$\Psi\colon\mathbb{L}^{2}\longrightarrow\mathbb{L}^{2}$ as follows

$(\Psi\mathcal{G})(\cdot)=Y_{T}(\cdot).$

Then the BPP can be reformulated as the following operator equation

$\Psi\mathcal{G}=Y_{T},\quad Y_{T}\in\mathbb{L}^{2}.$ (3.1)

The following lemma shows that the operator
$\Psi\colon\mathbb{L}^{2}\longrightarrow\mathbb{L}^{2}$ is compact.

######  Lemma 3.1.

The input-output operator

$\Psi\colon\mathbb{L}^{2}\ni\mathcal{G}\longmapsto
Y(T,\cdot,\mathcal{G})\in\mathbb{L}^{2}$

is a linear compact operator.

###### Proof.

The proof is similar to the proof of [3, Lemma 2]. ∎

Consequently, the inverse problem (3.1) is ill-posed (unstable). Therefore, we
define a quasi-solution $\mathcal{G}_{*}$ of the BPP as a minimizer of the
problem

$\displaystyle\mathcal{J}(\mathcal{G}_{*})=\inf_{\mathcal{G}\in\mathbb{L}^{2}}\mathcal{J}(\mathcal{G}),$
(3.2) $\displaystyle\mathcal{J}(\mathcal{G})$
$\displaystyle=\frac{1}{2}\|Y(T,\cdot,\mathcal{G})-Y_{T}^{\delta}\|_{\mathbb{L}^{2}}^{2},\qquad\mathcal{G}\in\mathbb{L}^{2},$
(3.3)

where $Y_{T}^{\delta}=(y_{T}^{\delta},y_{T,\Gamma}^{\delta})$ is a noisy data
of $Y_{T}$ such that $\|Y_{T}-Y_{T}^{\delta}\|\leq\delta$ for $\delta\geq 0$.

As a remedy to the ill-posedness of the BPP, the regularized Tikhonov
functional is defined as

$\displaystyle\mathcal{J}_{\varepsilon}(\mathcal{G})$
$\displaystyle=\frac{1}{2}\|Y(T,\cdot,\mathcal{G})-Y_{T}^{\delta}\|_{\mathbb{L}^{2}}^{2}+\frac{\varepsilon}{2}\|\mathcal{G}\|_{\mathbb{L}^{2}}^{2},\qquad\mathcal{G}\in\mathbb{L}^{2},$
(3.4)

where $\varepsilon>0$ is the regularizing parameter.

Although the problem is unstable, one can obtain some conditional stability
provided that a priori bound on unknown initial data is known. This is shown
in the following lemma using the logarithmic convexity method [25].

######  Lemma 3.2.

The solution $Y(t)$ of the system (1.1) satisfies the following estimate

$\|Y(t)\|_{\mathbb{L}^{2}}\leq\|Y(0)\|_{\mathbb{L}^{2}}^{1-\frac{t}{T}}\|Y(T)\|_{\mathbb{L}^{2}}^{\frac{t}{T}},\qquad
0\leq t\leq T.$ (3.5)

###### Proof.

Let $\mathcal{G}:=(f,g)\neq(0,0)$. Let us define the energy of the system
(1.1) as

$E(t)=\frac{1}{2}\|Y(t)\|_{\mathbb{L}^{2}}^{2},\qquad t\in[0,T].$

Using the equations of system (1.1) and integrating by parts in $\Omega$ and
on $\Gamma$, we obtain

$\displaystyle E^{\prime}(t)$ $\displaystyle=\left\langle
Y(t),\partial_{t}Y(t)\right\rangle_{\mathbb{L}^{2}}$
$\displaystyle=\int_{\Omega}y\partial_{t}y\,\mathrm{d}x+\int_{\Gamma}y_{\Gamma}\partial_{t}y_{\Gamma}\,\mathrm{d}S$
$\displaystyle=\int_{\Omega}y(d\Delta
y-ay)\,\mathrm{d}x+\int_{\Gamma}y_{\Gamma}(\gamma\Delta_{\Gamma}y_{\Gamma}-d\partial_{\nu}y-by_{\Gamma})\,\mathrm{d}S$
$\displaystyle=-d\int_{\Omega}|\nabla
y|^{2}\,\mathrm{d}x-\int_{\Omega}ay^{2}\,\mathrm{d}x-\gamma\int_{\Gamma}|\nabla_{\Gamma}y_{\Gamma}|^{2}\,\mathrm{d}S-\int_{\Gamma}by_{\Gamma}^{2}\,\mathrm{d}S.$

Differentiating in time $t$ and integrating by parts again using (1.1), we
obtain

$\displaystyle E^{\prime\prime}(t)$
$\displaystyle=2\int_{\Omega}\left(-d\nabla
y\cdot\nabla(\partial_{t}y)-ay\partial_{t}y\right)\,\mathrm{d}x+2\int_{\Gamma}\left(-\gamma\nabla_{\Gamma}y_{\Gamma}\cdot\nabla_{\Gamma}(\partial_{t}y_{\Gamma})-by_{\Gamma}\partial_{t}y_{\Gamma}\right)\,\mathrm{d}S$
$\displaystyle=2\|\partial_{t}Y(t)\|_{\mathbb{L}^{2}}^{2}.$

By Cauchy-Schwarz inequality, we infer that

$E(t)E^{\prime\prime}(t)-(E^{\prime}(t))^{2}=\|Y(t)\|_{\mathbb{L}^{2}}^{2}\|\partial_{t}Y(t)\|_{\mathbb{L}^{2}}^{2}-\left\langle
Y(t),\partial_{t}Y(t)\right\rangle^{2}_{\mathbb{L}^{2}}\geq 0.$

It follows then that

$(\log
E(t))^{\prime\prime}=\frac{E(t)E^{\prime\prime}(t)-(E^{\prime}(t))^{2}}{E^{2}(t)}\geq
0.$

This implies that $\log E(t)$ is a convex function of $t$, which yields

$\log E(t)\leq\left(1-\frac{t}{T}\right)\log E(0)+\frac{t}{T}\log E(T).$

This completes the proof. ∎

Now if we assume that $\|Y(0)\|_{\mathbb{L}^{2}}\leq M$ in (3.5), where $M$ is
a given positive constant, the following conditional stability holds

$\|Y(t)\|_{\mathbb{L}^{2}}\leq
M^{1-\frac{t}{T}}\|Y(T)\|_{\mathbb{L}^{2}}^{\frac{t}{T}},\qquad 0\leq t\leq
T.$

Such a result would be of practical interest if one can compute the prior
bound $M$ from the accessible data of the physical model.

Next, we derive a gradient formula for $\mathcal{J}$ via the mild solution
$\Phi=(\varphi,\varphi_{\Gamma})$ of an adjoint system.

######  Proposition 3.3.

The cost functional $\mathcal{J}$ is Fréchet differentiable and its gradient
at each $\mathcal{G}\in\mathbb{L}^{2}$ is given by

$\mathcal{J}^{\prime}(\mathcal{G})(x)=\Phi(x,0,\mathcal{G})\qquad\text{for a.e
}\,x\in\Omega\times\Gamma,$ (3.6)

where $\Phi(t,\cdot,\mathcal{G})=(\varphi,\varphi_{\Gamma})$ is the mild
solution of the following adjoint system

$\displaystyle\begin{aligned}
&-\partial_{t}\varphi-d\Delta\varphi+a(x)\varphi=0,&&\qquad\text{in
}\Omega_{T},\\\
&-\partial_{t}\varphi_{\Gamma}-\gamma\Delta_{\Gamma}\varphi_{\Gamma}+d\partial_{\nu}\varphi+b(x)\varphi_{\Gamma}=0,&&\qquad\text{on
}\Gamma_{T},\\\ &\varphi_{\Gamma}(t,x)=\varphi_{|\Gamma}(t,x),&&\qquad\text{on
}\Gamma_{T},\\\
&\varphi\rvert_{t=T}=y(T,\cdot,\mathcal{G})-y_{T}^{\delta},&&\qquad\text{in
}\Omega,\\\
&\varphi_{\Gamma}\rvert_{t=T}=y_{\Gamma}(T,\cdot,\mathcal{G})-y_{T,\Gamma}^{\delta},&&\qquad\text{on
}\Gamma.\end{aligned}$ (3.7)

###### Proof.

Let $\mathcal{G},\mathcal{G}+\delta\mathcal{G}\in\mathbb{L}^{2}$ be given
initial data and $Y(x,t,\mathcal{G})$, $Y(x,t,\mathcal{G}+\delta\mathcal{G})$
be corresponding solutions of (1.1). Hence, By linearity of the systems,

$\delta
Y(x,t,\mathcal{G}):=Y(x,t,\mathcal{G}+\delta\mathcal{G})-Y(x,t,\mathcal{G}),$

is the mild solution of the following system

$\displaystyle\begin{aligned} &\partial_{t}\delta y-d\Delta\delta y+a(x)\delta
y=0,&&\qquad\text{in }\Omega_{T},\\\ &\partial_{t}\delta
y_{\Gamma}-\gamma\Delta_{\Gamma}\delta y_{\Gamma}+d\partial_{\nu}\delta
y+b(x)\delta y_{\Gamma}=0,&&\qquad\text{on }\Gamma_{T},\\\ &\delta
y_{\Gamma}(t,x)=\delta y_{|\Gamma}(t,x),&&\qquad\text{on }\Gamma_{T},\\\
&(\delta y,\delta y_{\Gamma})\rvert_{t=0}=(\delta f,\delta
g),&&\qquad\Omega\times\Gamma\end{aligned}$ (3.8)

Let us compute the difference

$\delta\mathcal{J}(\mathcal{G}):=\mathcal{J}(\mathcal{G}+\delta\mathcal{G})-\mathcal{J}(\mathcal{G}).$

We have

$\displaystyle\delta\mathcal{J}(\mathcal{G})$
$\displaystyle=\frac{1}{2}\|Y(T,\cdot,\mathcal{G}+\delta\mathcal{G})-Y_{T}^{\delta}\|_{\mathbb{L}^{2}}^{2}-\frac{1}{2}\|Y(T,\cdot,\mathcal{G})-Y_{T}^{\delta}\|_{\mathbb{L}^{2}}^{2},$
$\displaystyle=\frac{1}{2}\left(\|y(T,\cdot,\mathcal{G}+\delta\mathcal{G})-y_{T}^{\delta}\|_{L^{2}(\Omega)}^{2}+\|y_{\Gamma}(T,\cdot,\mathcal{G}+\delta\mathcal{G})-y_{T,\Gamma}^{\delta}\|_{L^{2}(\Gamma)}^{2}\right)$
$\displaystyle\qquad-\frac{1}{2}\left(\|y(T,\cdot,\mathcal{G})-y_{T}^{\delta}\|_{L^{2}(\Omega)}^{2}+\|y_{\Gamma}(T,\cdot,\mathcal{G})-y_{T,\Gamma}^{\delta}\|_{L^{2}(\Gamma)}^{2}\right)$
$\displaystyle=\frac{1}{2}\int_{\Omega}\left[(y(T,x,\mathcal{G}+\delta\mathcal{G})-y_{T}^{\delta}(x))^{2}-(y(T,x,\mathcal{G})-y_{T}^{\delta}(x))^{2}\right]\,\mathrm{d}x$
$\displaystyle\quad+\frac{1}{2}\int_{\Gamma}\left[(y_{\Gamma}(T,x,\mathcal{G}+\delta\mathcal{G})-y_{T,\Gamma}^{\delta}(x))^{2}-(y_{\Gamma}(T,x,\mathcal{G})-y_{T,\Gamma}^{\delta})^{2}\right]\,\mathrm{d}S.$

Using the identity

$\frac{1}{2}\left[(x-z)^{2}-(y-z)^{2}\right]=(y-z)(x-y)+\frac{1}{2}(x-y)^{2},\;x,y\in\mathbb{R},$

in the last two terms, we obtain

$\displaystyle\delta\mathcal{J}(\mathcal{G})$
$\displaystyle=\int_{\Omega}(y(T,x,\mathcal{G})-y_{T}^{\delta}(x))\delta
y(T,x,\mathcal{G})\,\mathrm{d}x+\frac{1}{2}\int_{\Omega}[\delta
y(T,x,\mathcal{G})]^{2}\,\mathrm{d}x$ (3.9)
$\displaystyle+\int_{\Gamma}(y_{\Gamma}(T,x,\mathcal{G})-y_{T,\Gamma}^{\delta}(x))\delta
y_{\Gamma}(T,x,\mathcal{G})\,\mathrm{d}S+\frac{1}{2}\int_{\Gamma}[\delta
y_{\Gamma}(T,x,\mathcal{G})]^{2}\,\mathrm{d}S,$ (3.10)

We rewrite the first integral in the right-hand side of (3.9) using
$\Phi(t,\cdot,\mathcal{G})$ and $\delta Y(t,\cdot,\mathcal{G})$, the mild
solutions of (3.7) and (3.8) respectively. We have

$\displaystyle\int_{\Omega}(y(T,x,\mathcal{G})-y_{T}^{\delta}(x))\delta
y(T,x,\mathcal{G})\,\mathrm{d}x=\int_{\Omega}\varphi(T,x,\mathcal{G})\delta
y(T,x,\mathcal{G})\,\mathrm{d}x$
$\displaystyle=\bigintss_{\Omega}\left[\int_{0}^{T}\partial_{t}(\varphi(t,x,\mathcal{G})\delta
y(t,x,\mathcal{G}))\,\mathrm{d}t\right]\,\mathrm{d}x+\int_{\Omega}\varphi(0,x)\,\delta
f(x)\mathrm{d}x$
$\displaystyle=\int_{\Omega_{T}}\left[(\partial_{t}\varphi)\delta
y+\varphi\partial_{t}(\delta
y)\right]\,\mathrm{d}x\,\mathrm{d}t+\int_{\Omega}\varphi(0,x)\,\delta
f(x)\mathrm{d}x$
$\displaystyle=\int_{\Omega_{T}}\left[(-d\Delta\varphi+a(x)\varphi)\delta
y+\varphi(d\Delta(\delta y)-a(x)(\delta y))\right]\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle\qquad+\int_{\Omega}\varphi(0,x)\,\delta f(x)\mathrm{d}x$
$\displaystyle=\int_{\Omega_{T}}-d[(\Delta\varphi)\delta y-\Delta(\delta
y)\varphi]\,\mathrm{d}x\,\mathrm{d}t+\int_{\Omega}\varphi(0,x)\,\delta
f(x)\mathrm{d}x$
$\displaystyle=\int_{\Gamma_{T}}-d[(\partial_{\nu}\varphi)\delta
y_{\Gamma}-\varphi_{\Gamma}\partial_{\nu}(\delta
y)]\,\mathrm{d}S\,\mathrm{d}t+\int_{\Omega}\varphi(0,x)\,\delta
f(x)\mathrm{d}x.$ (3.11)

Similarly, for the first integral in the right-hand side of (3.10), we obtain

$\displaystyle\int_{\Gamma}(y_{\Gamma}(T,x,\mathcal{G})-y_{T,\Gamma}^{\delta}(x))\delta
y_{\Gamma}(T,x,\mathcal{G})\,\mathrm{d}S$
$\displaystyle=\int_{\Gamma_{T}}-\gamma[(\Delta_{\Gamma}\varphi_{\Gamma})\delta
y_{\Gamma}-\Delta_{\Gamma}(\delta
y_{\Gamma})\varphi_{\Gamma}]\,\mathrm{d}S\,\mathrm{d}t$ (3.12)
$\displaystyle\qquad+\int_{\Gamma_{T}}d[(\partial_{\nu}\varphi)\delta
y_{\Gamma}-\varphi_{\Gamma}\partial_{\nu}(\delta
y)]\,\mathrm{d}S\,\mathrm{d}t+\int_{\Gamma}\varphi_{\Gamma}(0,x)\,\delta
g(x)\,\mathrm{d}S$
$\displaystyle=\int_{\Gamma_{T}}d[(\partial_{\nu}\varphi)\delta
y_{\Gamma}-\varphi_{\Gamma}\partial_{\nu}(\delta
y)]\,\mathrm{d}S\,\mathrm{d}t+\int_{\Gamma}\varphi_{\Gamma}(0,x)\,\delta
g(x)\,\mathrm{d}S$ (3.13)

The first integral in the right-hand side of (3.12) is null by the surface
divergence formula (1.3). Adding up the two integrals (3.11) and (3.13), we
obtain

$\displaystyle\int_{\Omega}(y(T,x,\mathcal{G})-y_{T}^{\delta}(x))\delta
y(T,x,\mathcal{G})\,\mathrm{d}x+\int_{\Gamma}(y_{\Gamma}(T,x,\mathcal{G})-y_{T,\Gamma}^{\delta}(x))\delta
y_{\Gamma}(T,x,\mathcal{G})\,\mathrm{d}S$
$\displaystyle=\int_{\Omega_{T}}\varphi(0,x)\,\delta
f(x)\,\mathrm{d}x+\int_{\Gamma_{T}}\varphi_{\Gamma}(0,x)\,\delta
g(x)\,\mathrm{d}S.$

The second integrals in the right-hand sides of (3.9) and (3.10) are of order
$\mathcal{O}\left(\|\delta\mathcal{G}\|^{2}_{\mathbb{L}^{2}}\right)$. In fact,
multiplying (3.8)1 by $\delta y$, (3.8)2 by $\delta y_{\Gamma}$ and using
Green and surface divergence formula, we obtain

$\dfrac{1}{2}\partial_{t}\left(\int_{\Omega}|\delta
y|^{2}\,\mathrm{d}x\right)+d\int_{\Omega}|\nabla(\delta
y)|^{2}\,\mathrm{d}x+\dfrac{d}{2}\int_{\Gamma}\partial_{\nu}(|\delta
y|^{2})\,\mathrm{d}S+\int_{\Omega}a|\delta y|^{2}\,\mathrm{d}x=0,$ (3.14)
$\dfrac{1}{2}\partial_{t}\left(\int_{\Gamma}|\delta
y_{\Gamma}|^{2}\,\mathrm{d}S\right)+\gamma\int_{\Gamma}|\nabla_{\Gamma}(\delta
y_{\Gamma})|^{2}\,\mathrm{d}S-\dfrac{d}{2}\int_{\Gamma}\partial_{\nu}(|\delta
y|^{2})\,\mathrm{d}S+\int_{\Gamma}b|\delta y_{\Gamma}|^{2}\,\mathrm{d}S=0.$
(3.15)

Combining (3.14) and (3.15), we arrive at

$\displaystyle\dfrac{1}{2}\dfrac{\mathrm{d}}{\mathrm{d}t}\|\delta
Y(t)\|_{\mathbb{L}^{2}}^{2}+d\int_{\Omega}|\nabla(\delta
y)|^{2}\,\mathrm{d}x+\gamma\int_{\Gamma}|\nabla_{\Gamma}(\delta
y_{\Gamma})|^{2}\,\mathrm{d}S+\int_{\Omega}a|\delta y|^{2}\,\mathrm{d}x$
$\displaystyle\qquad+\int_{\Gamma}b|\delta y_{\Gamma}|^{2}\,\mathrm{d}S=0.$

By Cauchy-Schwarz inequality, we obtain

$\dfrac{1}{2}\dfrac{\mathrm{d}}{\mathrm{d}t}\|\delta
Y(t)\|_{\mathbb{L}^{2}}^{2}+\min(d,\gamma)\|(\nabla(\delta
y),\nabla_{\Gamma}(\delta y_{\Gamma}))\|_{\mathbb{L}^{2}}^{2}\leq D\|\delta
Y(t)\|_{\mathbb{L}^{2}}^{2}.$

Therefore,

$\dfrac{1}{2}\dfrac{\mathrm{d}}{\mathrm{d}t}\|\delta
Y(t)\|_{\mathbb{L}^{2}}^{2}\leq D\|\delta Y(t)\|_{\mathbb{L}^{2}}^{2}.$

Then Gronwall inequality imply that

$\|\delta Y(t)\|_{\mathbb{L}^{2}}^{2}\leq\mathrm{e}^{2DT}\|\delta
Y(0)\|_{\mathbb{L}^{2}}^{2}=\mathrm{e}^{2DT}\|\delta\mathcal{G}\|_{\mathbb{L}^{2}}^{2},$
(3.16)

for every $t\in[0,T]$. Consequently,

$\frac{1}{2}\int_{\Omega}[\delta
y(T,x,\mathcal{G})]^{2}\,\mathrm{d}x+\frac{1}{2}\int_{\Gamma}[\delta
y_{\Gamma}(T,x,\mathcal{G})]^{2}\,\mathrm{d}S=\mathcal{O}\left(\|\delta\mathcal{G}\|^{2}_{\mathbb{L}^{2}}\right).$

This completes the proof. ∎

Now, we show the Lipschitz continuity of the Fréchet gradient
$\mathcal{J}^{\prime}$.

######  Lemma 3.4.

Let $\mathcal{G},\delta\mathcal{G}\in\mathbb{L}^{2}$. The Fréchet gradient
$\mathcal{J}^{\prime}$ is Lipschitz continuous,

$\|\mathcal{J}^{\prime}(\mathcal{G}+\delta\mathcal{G})-\mathcal{J}^{\prime}(\mathcal{G})\|_{\mathbb{L}^{2}}\leq
L\|\delta\mathcal{G}\|_{\mathbb{L}^{2}},$

where the Lipschitz constant $L>0$ depends on $T$ and $D$ as follows

$L=\sqrt{\mathrm{e}^{2DT}\left(1+2TD\mathrm{e}^{2TD}\right)}.$ (3.17)

###### Proof.

By Proposition (3.3), we have

$\|\mathcal{J}^{\prime}(\mathcal{G}+\delta\mathcal{G})-\mathcal{J}^{\prime}(\mathcal{G})\|^{2}_{\mathbb{L}^{2}}:=\|\delta\Phi(0,\cdot,\mathcal{G})\|_{\mathbb{L}^{2}}^{2},$
(3.18)

where
$\delta\Phi(t,\cdot,\mathcal{G}):=\Phi(t,\cdot,\mathcal{G}+\delta\mathcal{G})-\Phi(t,\cdot,\mathcal{G})$
be the strong solution of the adjoint system

$\displaystyle\begin{aligned}
&-\partial_{t}(\delta\varphi)-d\Delta(\delta\varphi)+a(x)(\delta\varphi)=0,&\qquad\text{in
}\Omega_{T},\\\
&-\partial_{t}(\delta\varphi_{\Gamma})-\gamma\Delta_{\Gamma}(\delta\varphi_{\Gamma})+d\partial_{\nu}(\delta\varphi)+b(x)(\delta\varphi_{\Gamma})=0,&\qquad\text{on
}\Gamma_{T},\\\
&\delta\varphi_{\Gamma}(t,x)=(\delta\varphi)_{|\Gamma}(t,x),&\qquad\text{on
}\Gamma_{T},\\\ &(\delta\varphi,\delta\varphi_{\Gamma})\rvert_{t=T}=(\delta
y,\delta y_{\Gamma})\rvert_{t=T},&\qquad\Omega\times\Gamma.\end{aligned}$
(3.19)

Next, we estimate the norm
$\|\delta\Phi(0,\cdot,\mathcal{G})\|_{\mathbb{L}^{2}}^{2}$. In the adjoint
system (3.19), multiplying the first equation by $\delta\varphi$ and the
second by $\delta\varphi_{\Gamma}$, we obtain the following identities

$\displaystyle-\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left(\int_{\Omega}|\delta\varphi|^{2}\,\mathrm{d}x\right)+d\int_{\Omega}|\nabla(\delta\varphi)|^{2}\,\mathrm{d}x-\frac{1}{2}\int_{\Gamma}d\partial_{\nu}(|\delta\varphi|^{2})\,\mathrm{d}S+\int_{\Omega}a|\delta\varphi|^{2}\,\mathrm{d}x=0,$
$\displaystyle-\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left(\int_{\Gamma}|\delta\varphi_{\Gamma}|^{2}\,\mathrm{d}S\right)+\gamma\int_{\Gamma}|\nabla_{\Gamma}(\delta\varphi_{\Gamma})|^{2}\,\mathrm{d}S+\frac{1}{2}\int_{\Gamma}d\partial_{\nu}(|\delta\varphi|^{2})\,\mathrm{d}S+\int_{\Gamma}b|\delta\varphi_{\Gamma}|^{2}\,\mathrm{d}S=0,$

where we employed Green and the surface divergence formulae. Combining the
last identities, we obtain

$\displaystyle\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left(\int_{\Omega}|\delta\varphi|^{2}\,\mathrm{d}x+\int_{\Gamma}|\delta\varphi_{\Gamma}|^{2}\,\mathrm{d}S\right)-\int_{\Omega}a|\delta\varphi|^{2}\,\mathrm{d}x-\int_{\Gamma}b|\delta\varphi_{\Gamma}|^{2}\,\mathrm{d}S$
$\displaystyle=d\int_{\Omega}|\nabla(\delta\varphi)|^{2}\,\mathrm{d}x+\gamma\int_{\Gamma}|\nabla_{\Gamma}(\delta\varphi_{\Gamma})|^{2}\,\mathrm{d}S.$
(3.20)

Then

$\displaystyle\dfrac{\mathrm{d}}{\mathrm{d}t}\|\delta\Phi(t,\cdot,\mathcal{G})\|_{\mathbb{L}^{2}}^{2}-2\int_{\Omega}a|\delta\varphi|^{2}\,\mathrm{d}x-2\int_{\Gamma}b|\delta\varphi_{\Gamma}|^{2}\,\mathrm{d}S\geq
0.$

Integrating between $0$ and $T$ and using the fact that
$\delta\Phi(T,\cdot,\mathcal{G})=\delta Y(T,\cdot,\mathcal{G})$, we obtain

$\|\delta\Phi(0,\cdot,\mathcal{G})\|_{\mathbb{L}^{2}}^{2}\leq\|\delta
Y(T,\cdot,\mathcal{G})\|_{\mathbb{L}^{2}}^{2}+2D\|\Phi\|_{\mathbb{L}_{T}^{2}}^{2}.$
(3.21)

Now, we estimate $\|\Phi\|_{\mathbb{L}_{T}^{2}}^{2}$, by (3.20), we obtain

$\dfrac{\mathrm{d}}{\mathrm{d}t}\|\delta\Phi(t,\cdot,\mathcal{G})\|_{\mathbb{L}^{2}}^{2}+2D\|\delta\Phi(t,\cdot,\mathcal{G})\|_{\mathbb{L}^{2}}^{2}\geq
0.$

This inequality implies that the function $h$ define by

$h(t)=\mathrm{e}^{2Dt}\|\delta\Phi(t,\cdot,\mathcal{G})\|_{\mathbb{L}^{2}}^{2},$

is non-decreasing on $[0,T]$. Therefore,

$\displaystyle\|\delta\Phi(t,\cdot,\mathcal{G})\|_{\mathbb{L}_{T}^{2}}^{2}$
$\displaystyle=\int_{0}^{T}\|\delta\Phi(t,\cdot,\mathcal{G})\|_{\mathbb{L}^{2}}^{2}\,\mathrm{d}t$
$\displaystyle\leq\int_{0}^{T}\mathrm{e}^{2Dt}\|\delta\Phi(t,\cdot,\mathcal{G})\|_{\mathbb{L}^{2}}^{2}\,\mathrm{d}t$
$\displaystyle\leq
T\mathrm{e}^{2DT}\|\delta\Phi(T,\cdot,\mathcal{G})\|_{\mathbb{L}^{2}}^{2}.$

Using the last inequality and (3.19), we arrive at

$\|\delta\Phi(t,\cdot,\mathcal{G})\|_{\mathbb{L}_{T}^{2}}^{2}\leq
T\mathrm{e}^{2DT}\|\delta Y(T,\cdot,\mathcal{G})\|_{\mathbb{L}^{2}}^{2}.$
(3.22)

From (3.21) and (3.22), we obtain

$\|\delta\Phi(0,\cdot,\mathcal{G})\|_{\mathbb{L}^{2}}^{2}\leq\left(1+2TD\mathrm{e}^{2DT}\right)\|Y(T,\cdot,\mathcal{G})\|_{\mathbb{L}^{2}}^{2}.$

Using inequality (3.16), we obtain

$\displaystyle\|\delta\Phi(0,\cdot,\mathcal{G})\|_{\mathbb{L}^{2}}^{2}\leq\mathrm{e}^{2DT}\left(1+2TD\mathrm{e}^{2DT}\right)\|\delta\mathcal{G}\|_{\mathbb{L}^{2}}^{2}.$

With the equality (3.18), this implies that

$\displaystyle\|\mathcal{J}^{\prime}(\mathcal{G}+\delta\mathcal{G})-\mathcal{J}^{\prime}(\mathcal{G})\|^{2}_{\mathbb{L}^{2}_{T}}\leq
2\mathrm{e}^{2DT}\left(1+2TD\mathrm{e}^{2DT}\right)\|\delta\mathcal{G}\|_{\mathbb{L}^{2}}^{2}.$

This yields the desired result. ∎

##  4\. Spectral analysis of the BPP

Here, we will give a necessary and sufficient condition for the existence of
the solution to BPP. We also provide a series representation of the unique
solution to (3.4) for the system (1.1) when $a=0$ and $b=0$. For this purpose,
we need the following lemma.

######  Lemma 4.1.

The input-output operator
$\Psi\colon\mathbb{L}^{2}\longrightarrow\mathbb{L}^{2}$ defined in (3.1) is
self-adjoint. Moreover,

$\left(\Psi\phi_{k}\right)(x)=\mathrm{e}^{-\lambda_{k}T}\phi_{k}(x)\qquad\text{
in }\Omega\times\Gamma,$ (4.1)

where $(\phi_{k})_{k\geq 1}$ is the sequence pf eigenfunctions of the operator
$-\mathcal{A}_{0}$ corresponding to the eigenvalues $(\lambda_{k})_{k\geq 1}$
given by Lemma 2.1.

###### Proof.

By virtue of the orthonormal basis $(\phi_{k})_{k}\subset\mathbb{L}^{2}$, the
solution of the system (ACP) corresponding to the initial data $\mathcal{G}$
and the operator $\mathcal{A}_{0}$ can be written as

$Y(t,x)=\sum_{k=1}^{\infty}\left\langle\mathcal{G},\phi_{k}\right\rangle\mathrm{e}^{-\lambda_{k}t}\phi_{k}\qquad\text{
in }\Omega\times\Gamma.$

Then, the Fourier series representation of the input-output operator is given
by

$(\Psi\mathcal{G})(x)=Y(T,x,\mathcal{G})=\sum_{k=1}^{\infty}\left\langle\mathcal{G},\Phi_{k}\right\rangle\mathrm{e}^{-\lambda_{k}T}\phi_{k}(x),$
(4.2)

Thus,

$(\Psi\phi_{n})(x)=\mathrm{e}^{-\lambda_{n}T}\phi_{n}(x)$

Moreover, the operator $\Psi$ is self adjoint. Indeed,

$\displaystyle\left\langle\Psi\mathcal{G}_{1},\mathcal{G}_{2}\right\rangle$
$\displaystyle=\left\langle\sum_{k=1}^{\infty}\left\langle\mathcal{G}_{1},\phi_{k}\right\rangle\mathrm{e}^{-\lambda_{k}T}\phi_{k},\sum_{k=1}^{\infty}\left\langle\mathcal{G}_{2},\phi_{k}\right\rangle\phi_{k}\right\rangle$
$\displaystyle=\sum_{k=1}^{\infty}\left\langle\mathcal{G}_{1},\phi_{k}\right\rangle\left\langle\mathcal{G}_{2},\phi_{k}\right\rangle\mathrm{e}^{-\lambda_{k}T}$
$\displaystyle=\left\langle\mathcal{G}_{1},\Psi\mathcal{G}_{2}\right\rangle.$

∎

Next, we discuss the solvability of the BPP.

######  Theorem 4.2.

Let $Y_{T}:=\left(y(T,\cdot),y_{\Gamma}(T,\cdot)\right)\in\mathbb{L}^{2}$ be a
noise free measured data. Then the BPP in (3.1) admits a solution if and only
if

$\sum_{k=1}^{\infty}\mathrm{e}^{\lambda_{k}T}Y_{T,k}^{2}<\infty,\qquad\text{
where }\,Y_{T,k}:=\left\langle Y_{T},\phi_{k}\right\rangle,$

in that case the solution is given by

$\mathcal{G}=\sum_{k=1}^{\infty}\mathrm{e}^{\lambda_{k}T}Y_{T,n}\phi_{k}.$

###### Proof.

Let $\mathcal{G}\in\mathbb{L}^{2}$. Using (4.2), $\Psi\mathcal{G}$ is given by

$Y(T,x,\mathcal{G})=\sum_{k=1}^{\infty}\left\langle\mathcal{G},\phi_{k}\right\rangle\mathrm{e}^{-\lambda_{k}T}\phi_{k}(x)\qquad\text{
in }\Omega\times\Gamma.$

Moreover, the output data can be expressed as follows

$Y_{T}=\sum_{k=1}^{\infty}Y_{T,k}\phi_{k}.$

Then, the Tikhonov functional (3.4) is represented by

$\displaystyle\mathcal{J}_{\varepsilon}(\mathcal{G})$
$\displaystyle=\frac{1}{2}\left\|\sum_{k=1}^{\infty}\left(\left\langle\mathcal{G},\phi_{k}\right\rangle\mathrm{e}^{-\lambda_{k}T}-\left\langle
Y_{T},\phi_{k}\right\rangle\right)\phi_{k}\right\|^{2}+\frac{\varepsilon}{2}\left\|\sum_{k=1}^{\infty}\left\langle\mathcal{G},\phi_{k}\right\rangle\phi_{k}\right\|^{2}$
$\displaystyle=\frac{1}{2}\sum_{k=1}^{\infty}\left(\left\langle\mathcal{G},\phi_{k}\right>\mathrm{e}^{-\lambda_{k}T}-Y_{T,k}\right)^{2}+\varepsilon\left\langle\mathcal{G},\phi_{k}\right\rangle^{2}.$

Clearly, this functional achieves its minimum value if

$\left\langle\mathcal{G},\phi_{k}\right\rangle\mathrm{e}^{-\lambda_{k}T}-Y_{T,k}=0.$

Then, the $k$th Fourier coefficient of the unique minimum of the Tikhonov
functional (3.4) is given by

$\left\langle\mathcal{G},\Phi_{k}\right\rangle=\mathrm{e}^{\lambda_{k}T}Y_{T,k}.$

As a result, we obtain

$\mathcal{G}=\sum_{k=1}^{\infty}\left\langle\mathcal{G},\phi_{k}\right\rangle\phi_{k}=\sum_{k=1}^{\infty}\mathrm{e}^{\lambda_{k}T}Y_{T,k}\phi_{k}.$

∎

##  5\. Numerical reconstruction of initial temperatures

Now we are in position to present the CG algorithm for reconstructing initial
temperatures. We will consider both 1-D and 2-D initial temperatures
reconstruction.

1 Set $n=0$ and choose an initial temperature $\mathcal{G}_{0}$;

2 Solve the direct problem (1.1) to obtain $Y(t,x,\mathcal{G}_{0})$;

3 Knowing the computed $Y(T,x,\mathcal{G}_{0})$ and the measured
$Y_{T}^{\delta}$, solve the adjoint problem (3.7) to obtain
$\Phi(t,x,\mathcal{G}_{0})$;

4 Compute the gradient
$p_{0}=\mathcal{J}_{\varepsilon}^{\prime}\left(\mathcal{G}_{0}\right)$ using
(3.6);

5 Solve the direct problem (1.1) with initial datum $p_{n}$ to obtain the
solution $\Psi p_{k}$;

6 Compute the relaxation parameter
$\displaystyle\alpha_{n}=\frac{\|\mathcal{J}_{\varepsilon}^{\prime}(\mathcal{G}_{n})\|_{\mathbb{L}^{2}}^{2}}{\|\Psi
p_{n}\|_{\mathbb{L}^{2}}^{2}+\varepsilon\|p_{n}\|_{\mathbb{L}^{2}}^{2}}$ ;

7 Find the next iteration $\mathcal{G}_{n+1}=\mathcal{G}_{n}-\alpha_{n}p_{n}$;

Stop the iteration process if the stopping criterion
$\mathcal{J}_{\varepsilon}(\mathcal{G}_{n+1})<e_{\mathcal{J}}$ holds.
Otherwise, set $n:=n+1$ and compute
$\gamma_{n}=\frac{\|\mathcal{J}_{\varepsilon}^{\prime}(\mathcal{G}_{n})\|_{\mathbb{L}^{2}}^{2}}{\|\mathcal{J}_{\varepsilon}^{\prime}(\mathcal{G}_{n-1})\|_{\mathbb{L}^{2}}^{2}}\qquad\text{and}\qquad
p_{n}=\mathcal{J}_{\varepsilon}^{\prime}(\mathcal{G}_{n})+\gamma_{n}p_{n-1},$
and go to Step 6;

Algorithm 1 CG algorithm

In all numerical simulations, we assume that the measured output data contain
a random noise. Instead of using the exact data $Y_{T}$, we generate the noisy
data by the following formula

$Y_{T}^{\delta}(x)=Y_{T}(x)+p\times\|Y_{T}\|_{\mathbb{L}^{2}}\times\mathrm{RandomReal[]},$

where $p$ stands for the percentage of the noise level, and the function
RandomReal[] produces random real numbers.

Comparing to static boundary conditions, the numerical implementation of
dynamic boundary conditions is quite complicated. For this aim, the numerical
solutions for the direct problems and the corresponding adjoint problems will
be obtained by using the method of lines.

###  5.1. 1-D reconstruction of initial temperatures

Let $\Omega=(0,\ell),\;\ell>0$ and $T>0$. We seek to reconstruct the unknown
initial temperature $g(x)\in L^{2}(0,\ell)$ in the 1-D heat equation with
dynamic boundary conditions

$\displaystyle\begin{aligned}
&y_{t}(t,x)-y_{xx}(t,x)=0,&&\qquad(t,x)\in(0,T)\times(0,\ell),\\\
&y_{t}(t,0)-y_{x}(t,0)=0,&&\qquad t\in(0,T),\\\
&y_{t}(t,\ell)+y_{x}(t,\ell)=0,&&\qquad t\in(0,T),\\\
&y(0,x)=g(x),\quad(y(0,0),y(0,\ell))=(a,b),&&\qquad
x\in(0,\ell).\end{aligned}$ (5.1)

The problem then is to find the initial temperature $(g(x),a,b)$ by using the
final time overdetermination data.

As above, $Y(t,x,\mathcal{G}):=\left(y(t,x),y(t,0),y(t,\ell)\right)$ denote
the solution of the system (5.1), where $\mathcal{G}:=(g,a,b)$. The Tikhonov
functional is given by

$\displaystyle\mathcal{J}_{\varepsilon}(\mathcal{G})$
$\displaystyle=\frac{1}{2}\left\|Y(T,\cdot,\mathcal{G})-Y_{T}^{\delta}\right\|_{\mathbb{L}^{2}}^{2}+\frac{\varepsilon}{2}\|\mathcal{G}\|_{\mathbb{L}^{2}}^{2}$
(5.2)
$\displaystyle=\frac{1}{2}\left(\left\|y(T,\cdot)-y_{T}^{\delta}\right\|_{L^{2}(0,\ell)}^{2}+\left|y(T,0)-y_{T}^{0,\delta}\right|^{2}+\left|y(T,\ell)-y_{T}^{\ell,\delta}\right|^{2}+\varepsilon\|\mathcal{G}\|_{\mathbb{L}^{2}}^{2}\right),$

where
$Y_{T}^{\delta}:=\left(y_{T}^{\delta},y_{T}^{0,\delta},y_{T}^{\ell,\delta}\right)\in\mathbb{L}^{2}$
is the observation data.

In this 1-D case, the governing operator takes the form

$\mathcal{A}_{0}:=\left(\begin{array}[]{ccc}\partial_{xx}&0&0\\\
\left.\partial_{x}\right|_{x=0}&0&0\\\
-\left.\partial_{x}\right|_{x=\ell}&0&0\end{array}\right),\qquad
D\left(\mathcal{A}_{0}\right)=\mathbb{H}^{2}=\left\\{(f,f(0),f(\ell)):f\in
H^{2}(0,\ell)\right\\}.$

######  Remark 5.1.

It is well-known that the decay to zero of singular values $\sigma_{k}$ of the
input-output operator $\Psi$ indicates the degree of ill-posedness of the
corresponding inverse problem. It has been shown in [19, Proposition 5.3.1]
that the eigenvalues of $-\mathcal{A}_{0}$ are of order
$\lambda_{k}=\mathcal{O}\left(k^{2}\right)$. It follows that the singular
values $\sigma_{k}=\mathrm{e}^{-\lambda_{k}T}$ of $\Psi$ are of order
$\mathcal{O}\left(\mathrm{e}^{-k}\right)$. This means that the BPP we consider
is severely exponentially ill-posed.

The adjoint system corresponding to (5.1) is given by

$\displaystyle\begin{aligned}
&\varphi_{t}(t,x)+\varphi_{xx}(t,x)=0,&&(t,x)\in(0,T)\times(0,\ell),\\\
&\varphi_{t}(t,0)+\varphi_{x}(t,0)=0,&&t\in(0,T),\\\
&\varphi_{t}(t,\ell)-\varphi_{x}(t,\ell)=0,&&t\in(0,T),\\\
&\varphi(T,x)=y(T,x)-y_{T}^{\delta}(x),&&x\in(0,\ell),\\\
&(\varphi(T,0),\varphi(T,\ell))=\left(y(T,0)-y_{T}^{0,\delta},y(T,\ell)-y_{T}^{\ell,\delta}\right).\end{aligned}$
(5.3)

Similarly to Proposition 3.3, we obtain the following formula for the gradient
of $\mathcal{J}_{\varepsilon}$

$\mathcal{J}_{\varepsilon}^{\prime}(\mathcal{G})(x)=\left(\varphi(0,x,\mathcal{G}),\varphi(0,0,\mathcal{G}),\varphi(0,\ell,\mathcal{G})\right)+\varepsilon\mathcal{G}.$

Next, we discretize the space domain $[0,\ell]$ into a uniform grid
$(x_{j})_{j=0}^{N_{x}}$ of step $\Delta x=\frac{\ell}{N_{x}}$ using finite
differences. Denoting by $y_{j}(t)=y(t,x_{j})$, the the second-order
derivative of $y$ in system (5.1) is approximated as follows

$y_{xx}(t,x_{j})\approx\frac{y_{j-1}(t)-2y_{j}(t)+y_{j+1}(t)}{(\Delta
x)^{2}},\qquad j=\overline{1,N_{x}-1}.$

The first-order derivatives on the boundary can be approximated by

$\displaystyle y_{x}(t,0)$
$\displaystyle\approx\frac{y_{1}(t)-y_{0}(t)}{\Delta x}$ $\displaystyle
y_{x}(t,\ell)$ $\displaystyle\approx\frac{y_{N_{x}}(t)-y_{N_{x}-1}(t)}{\Delta
x}.$

Now, it suffices to numerically solve the resulting ordinary differential
system. This is done with help of the Wolfram language.

In the following, we will test the stability of our algorithm using some
numerical experiments. For convenience, in all experiments, the system
parameters are taken as follows

$T=0.03,\quad\ell=1,\quad N_{x}=25,\quad a=b=0.$

For each example, an illustration is provided for the convergence error

$e(n,g_{n}):=\|\Psi(g_{n},0,0)-Y_{T}\|_{L^{2}(0,1)\times\mathbb{R}^{2}}^{2},$

and the accuracy error

$E(n,g_{n}):=\|g-g_{n}\|_{L^{2}(0,1)}$

in terms of the iteration number $n$.

### Example 1

We take the exact initial temperature as

$g(x)=\frac{7}{24}\left(\sin(\pi x)+1-x\right)\left(\sin(\pi
x)+\sqrt{x}\right),\quad x\in(0,1).$

Let us set the starting iteration as $g_{0}=0$, the regularization parameter
$\varepsilon=10^{-8}$, and stopping parameter $e_{\mathcal{J}}=10^{-6}$. The
algorithm stops at iteration $n=9$ for $p\in\\{1\%,3\%,5\%\\}$.

Figure 1. Exact and recovered $g(x)$ from noisy data (left) and corresponding
errors for first iterations for noise free data (right)

### Example 2

We choose the exact initial temperature as

$g(x)=6(1-x)\log\left(1+x^{2}\right),\;x\in(0,1).$

The initial iteration is $f_{0}=0$. The regularization parameter is
$\varepsilon=10^{-8}$ and the stopping parameter is $e_{\mathcal{J}}=10^{-6}$.
The algorithm stops at iterations $n\in\\{39,42,43\\}$ for
$p\in\\{1\%,3\%,5\%\\}$, respectively.

Figure 2. Exact and recovered $g(x)$ from noisy data (left) and corresponding
errors for first iterations for noise free data (right)

######  Remark 5.2.

The above numerical tests demonstrate the stability, the accuracy and the
speed of the designed CG algorithm for the reconstruction of 1-D unknown
initial temperatures in heat equation with dynamic boundary conditions. It is
clear that the restoration of the initial temperature $g(x)$ is not much
affected by random noise $p$. Furthermore, Figs. 1 and 2 show the decay of the
convergence error and the accuracy error as the iteration number increases.
This fact validates the convergence of the regularized solutions to the exact
initial temperatures.

###  5.2. 2-D reconstruction of initial temperatures

Let $\Omega=\left\\{(x_{1},x_{2})\in\mathbb{R}^{2}\colon
x_{1}^{2}+x_{2}^{2}<1\right\\}$ be the unit disk. Its boundary
$\Gamma=\left\\{(x_{1},x_{2})\in\mathbb{R}^{2}\colon
x_{1}^{2}+x_{2}^{2}=1\right\\}=\mathbb{S}^{1}$ is the unit circle. In polar
coordinates,

$\Omega=\\{(r\cos\theta,r\sin\theta)\colon
r\in[0,1),\theta\in[0,2\pi)\\}\text{ and
}\Gamma=\\{(\cos\theta,\sin\theta)\colon\theta\in[0,2\pi)\\}.$

The gradient of any smooth function $y$ on $\Omega$ is given locally by

$\nabla y=\partial_{r}y\,e_{r}+\frac{1}{r}\partial_{\theta}y\,e_{\theta},$

where $e_{r}$ and $e_{\theta}$ are the unit vectors in $r$ and in $\theta$
directions respectively. At any point on $\Gamma$, the outer unit normal
vector is $\nu=e_{r}$. Then, by the previous formula the normal derivative of
$y$ is given by

$\partial_{\nu}y=\nabla y\cdot\nu\rvert_{\Gamma}=\partial_{r}y\,e_{r}\cdot
e_{r}\rvert_{\Gamma}=\partial_{r}y\rvert_{r=1}.$

The standard Laplacian $\Delta$ on $\Omega$ takes the form

$\Delta
y=\partial^{2}_{r}y+\frac{1}{r}\partial_{r}y+\frac{1}{r^{2}}\partial^{2}_{\theta}y.$

For any $(x_{1},x_{2})=(\cos\theta,\sin\theta)\in\Gamma$, the tangent space at
$(x_{1},x_{2})$ is
$T_{(x_{1},x_{2})}\Gamma=\text{span}{(-\sin\theta,\cos\theta)}$. Denoting
$\partial_{\theta}:=(-\sin\theta,\cos\theta)$, we have
$|\partial_{\theta}|=1$. Then the metric tensor here is a scalar
$\mathrm{g}=1$. Using the local formula (1.2), the Laplace-Beltrami operator
on $\Gamma$ is given by

$\Delta_{\Gamma}y_{\Gamma}=\partial^{2}_{\theta}y_{\Gamma}.$

Consequently, a simplified version of system (1.1) in polar coordinates reads
as follows

$\displaystyle\begin{aligned}
&\partial_{t}y-\left[\partial^{2}_{r}y+\frac{1}{r}\partial_{r}y+\frac{1}{r^{2}}\partial^{2}_{\theta}y\right]=0,&&\qquad
0<r<1,\;0\leq\theta\leq 2\pi\\\
&\partial_{t}y_{\Gamma}-\partial^{2}_{\theta}y_{\Gamma}+\partial_{r}y\rvert_{r=1}=0,&&\qquad
0\leq\theta\leq 2\pi,\\\
&y_{\Gamma}(t,\theta)=y(t,r,\theta)\rvert_{r=1},&&\qquad 0\leq\theta\leq
2\pi,\\\ &(y,y_{\Gamma})\rvert_{t=0}=(g,g_{\Gamma}),&&\qquad
0<r<1,\;0\leq\theta\leq 2\pi.\end{aligned}$ (5.4)

Our purpose is to determine the initial temperature $g(r,\theta)$ using the
final time overdetermination data. The corresponding adjoint system is given
by

$\displaystyle\begin{aligned}
&-\partial_{t}\varphi-\left[\partial^{2}_{r}\varphi+\frac{1}{r}\partial_{r}\varphi+\frac{1}{r^{2}}\partial^{2}_{\theta}\varphi\right]=0,&&\qquad
0<r<1,\;0\leq\theta\leq 2\pi\\\
&-\partial_{t}\varphi_{\Gamma}-\partial^{2}_{\theta}\varphi_{\Gamma}+\partial_{r}\varphi\rvert_{r=1}=0,&&\qquad
0\leq\theta\leq 2\pi,\\\
&\varphi_{\Gamma}(t,\theta)=\varphi(t,r,\theta)\rvert_{r=1},&&\qquad
0\leq\theta\leq 2\pi,\\\
&\varphi\rvert_{t=T}=y(T,\cdot,\mathcal{G})-y_{T}^{\delta},&&\qquad
0<r<1,\;0\leq\theta\leq 2\pi,\\\
&\varphi_{\Gamma}=y_{\Gamma}(T,\cdot,\mathcal{G})-y_{T,\Gamma}^{\delta}),&&\qquad
0<r<1,\;0\leq\theta\leq 2\pi.\end{aligned}$

Then the gradient formula for $\mathcal{J}_{\varepsilon}$ at
$\mathcal{G}=(g,g_{\Gamma})$ is given by

$\mathcal{J}_{\varepsilon}^{\prime}(\mathcal{G})(r,\theta)=\left(\varphi(0,r,\theta,\mathcal{G}),\varphi(0,r,\theta,\mathcal{G})\rvert_{r=1}\right)+\varepsilon\mathcal{G},\qquad
0<r<1,\;0\leq\theta\leq 2\pi.$

To solve direct and adjoint problems, we discretize the space domain $\Omega$
into a uniform grid $\left(r_{j},\theta_{n}\right),0\leq j\leq N_{r},\,0\leq
n\leq N_{\theta}$, of steps $\Delta r=\frac{1}{N_{r}}$ and
$\Delta\theta=\frac{2\pi}{N_{\theta}}$ using finite differences (Fig. 3).

Figure 3. Polar grid for the disk.

Multiplying both sides of $\eqref{2deq1to4}_{1}$ by $r^{2}$, we obtain

$r^{2}\partial_{t}y-\left[r^{2}\partial^{2}_{r}y+r\partial_{r}y+\partial^{2}_{\theta}y\right]=0.$
(5.5)

Denoting by $y_{j,n}(t)=y(t,r_{j},\theta_{j})$, the the second-order partial
derivatives of $y$ in system (5.4) are approximated as follows

$\displaystyle\partial_{r}^{2}y(t,r_{j},\theta_{n})$
$\displaystyle\approx\frac{y_{j-1,n}(t)-2y_{j,n}(t)+y_{j+1,n}(t)}{(\Delta
r)^{2}},\qquad j=\overline{1,N_{r}-1}$
$\displaystyle\partial_{\theta}^{2}y(t,r_{j},\theta_{n})$
$\displaystyle\approx\frac{y_{j,n-1}(t)-2y_{j,n}(t)+y_{j,n+1}(t)}{(\Delta\theta)^{2}},\qquad
n=\overline{1,N_{\theta}-1}.$

The partial derivatives on the boundary are approximated by

$\displaystyle\partial_{r}y\rvert_{r=1}$
$\displaystyle\approx\frac{y_{N_{r},n}(t)-y_{N_{r}-1,n}(t)}{\Delta r},\qquad
n=\overline{0,N_{\theta}}$
$\displaystyle\partial_{\theta}^{2}y_{\Gamma}(t,\theta_{n})$
$\displaystyle\approx\frac{y_{N_{r},n-1}(t)-2y_{N_{r},n}(t)+y_{N_{r},n+1}(t)}{(\Delta\theta)^{2}},\qquad
n=\overline{1,N_{\theta}-1}.$

Note that the boundary condition for $\theta$ is given by periodicity
$y(t,r,0)=y(t,r,2\pi)$ which corresponds to $y_{j,0}=y_{j,N_{\theta}}$. Since
the system (5.4) is singular at $r=0$, some boundedness condition should be
added to avoid the blow up. We refer to [29] for more details.

Henceforth, we take the final time $T=0.01$ and the mesh parameters
$N_{r}=N_{\theta}=25$. We define the error surface as the difference between
the exact initial temperature and the recovered one.

### Example 1

We take the exact initial temperature as the separate variables function

$g(r,\theta)=\sin(\pi r)\sin\left(\frac{\theta}{2}\right),\quad
r\in(0,1),\;\theta\in[0,2\pi].$

We set the initial iteration as $g_{0}=0$, the regularization parameter
$\varepsilon=10^{-8}$, and the stopping parameter $e_{\mathcal{J}}=5.82\times
10^{-5}$. The algorithm stops at iteration $n=32$ for a noise level $p=1\%$.

Figure 4. Exact $g(r,\theta)$ (left) and recovered one from noisy data (right)
for Example 1 Figure 5. Error surface for the noise level $p=1\%$.

### Example 2

We take the exact initial temperature

$g(r,\theta)=[r(1-r)]^{2}\cos\left(\frac{r\theta}{4}\right),\quad
r\in(0,1),\;\theta\in[0,2\pi].$

We set the initial iteration as $g_{0}=0$, the regularization parameter
$\varepsilon=10^{-8}$, and the stopping parameter $e_{\mathcal{J}}=2.92\times
10^{-7}$. The algorithm stops at iteration $n=60$ for a noise level $p=1\%$.

Figure 6. Exact $g(r,\theta)$ (left) and recovered one from noisy data (right)
for Example 2. Figure 7. Error surface for the noise level $p=1\%$.

######  Remark 5.3.

In the previous experiments, the optimal regularization parameter
$\varepsilon^{*}$ can be chosen in terms of the noise $\delta$ by the rule
$\varepsilon<1$ and $\frac{\delta^{2}}{\varepsilon}<1$, while the optimal
stopping parameter $e^{*}_{\mathcal{J}}$ can be estimated by the iteration
number $n$ from the convergence the accuracy errors. See [16] for more
details.

##  6\. Conclusion

The inverse problem of identifying unknown initial temperatures in heat
conduction models from the final time data is of interest in many engineering
problems and industrial applications. Such models with static boundary
conditions have been extensively studied and the theory is well developed,
while papers dealing with heat models with dynamic boundary conditions are few
and still at an earlier stage.

In this paper, we have investigated the inverse problem of recovering unknown
initial temperatures of a heat equation with dynamic boundary conditions.
Using an optimal control approach, the existence, uniqueness and stability of
the minimizers are discussed for a general case of bounded domains
$\Omega\subset\mathbb{R}^{N}$ ($N\geq 1$), and the numerical results are
obtained by the conjugate gradient method for the 1-D and 2-D cases. The
presented method has the advantage of yielding a monotone scheme and then fast
numerical results. We expect that the proposed numerical approach is still
applicable for such a problem in the 3-D case of system (1.1), i.e. $N=3$.
However, the geometric complexity of the domain $\Omega$ and the presence of
the Laplace-Beltrami operator on the boundary will make the analysis more
difficult. The study of such a problem would be of much interest.

## References

  * [1] E. M. Ait Ben Hassi, S. E. Chorfi and L. Maniar,  An inverse problem of radiative potentials and initial temperatures in parabolic equations with dynamic boundary conditions,  _J. Inverse Ill-Posed Probl._ , 30 (2022), 363–378.
  * [2] E. M. Ait Ben Hassi, S. E. Chorfi, L. Maniar and O. Oukdach,  Lipschitz stability for an inverse source problem in anisotropic parabolic equations with dynamic boundary conditions,  _Evol Equ and Cont Theo_ , 10 (2021), 837–859.
  * [3] E. M. Ait Ben Hassi, S. E. Chorfi and L. Maniar,  Identification of source terms in heat equation with dynamic boundary conditions,  _Math. Meth. Appl. Sci._ , 45 (2022), 2364–2379.
  * [4] M. Alahyane, I. Boutaayamou, A. Chrifi, Y. Echarroudi and Y. Ouakrim,  Numerical study of inverse source problem for internal degenerate parabolic equation,  _Int. J. Comput. Methods_ , 18 (2021) 2050032.
  * [5] K. Ammari and M. Choulli,  Logarithmic stability in determining two coefficients in a dissipative wave equation. Extensions to clamped Euler–Bernoulli beam and heat equations,  _Differ. Equ._ , 259 (2015), 3344–3365.
  * [6] K. Ammari and M. Choulli,  Hölder stability in determining the potential and the damping coefficient in a wave equation,  _J. Evol. Equ._ , 19 (2019), 305–319.
  * [7] J. Cheng and M. Yamamoto,  One new strategy for a priori choice of regularizing parameters in Tikhonov’s regularization,  _Inverse Problems_ 16 (2000), L31–L38.
  * [8] S. E. Chorfi, G. El Guermai, L. Maniar and W. Zouhair,  Identification of source terms in wave equation with dynamic boundary conditions,  _Math. Meth. Appl. Sci._ , (2022), doi: 10.1002/mma.8556.
  * [9] H. W. Engl, M. Hanke and A. Neubauer,  _Regularization of Inverse Problems_ ,  Kluwer Academic Publishers, 2000.
  * [10] K.-J. Engel and R. Nagel,  _One-Parameter Semigroups for Linear Evolution Equations_ ,  Graduate Texts in Mathematics 194, Springer-Verlag, 1999.
  * [11] A. Favini, J. A. Goldstein, G. R. Goldstein and S. Romanelli,  The heat equation with generalized Wentzell boundary condition,  _J. Evol. Equ._ , 2 (2002), 1–19.
  * [12] G. R. Goldstein,  Derivation and physical interpretation of general boundary conditions,  _Adv. Diff. Equ._ , 11 (2006), 457–480.
  * [13] A. Hasanov,  An inverse source problem with single Dirichlet type measured output data for a linear parabolic equation,  _Applied Mathematics Letters_ , 24 (2011), 1269–1273.
  * [14] A. Hasanov,  Simultaneous determination of source terms in a linear parabolic problem from the final overdetermination: Weak solution approach,  _J. Math. Anal. Appl._ , 330 (2007), 766–779.
  * [15] A. Hasanov,  Simultaneously identifying the thermal conductivity and radiative coefficient in heat equation from Dirichlet and Neumann boundary measured outputs,  _J. Inverse Ill-Posed Probl._ , 29 (2021), 81–91.
  * [16] A. Hasanov Hasanoğlu and V. G. Romanov,  _Introduction to Inverse Problems for Differential Equations_ ,  Springer, second edition, 2021.
  * [17] M. J. Huntul,  Identifying an unknown heat source term in the third-order pseudo-parabolic equation from nonlocal integral observation,  _Int. Commun. Heat Mass Transf._ , 128 (2021), 105550.
  * [18] V. Isakov,  _Inverse Problems for Partial Differential Equations_ ,  Springer, Cham, 2017.
  * [19] A. Khoutaibi,  _Null controllability of linear and semilinear parabolic equations with dynamic boundary conditions_ ,  PhD thesis, Cadi Ayyad University, 2020.
  * [20] A. Khoutaibi and L. Maniar,  Null controllability for a heat equation with dynamic boundary conditions and drift terms,  _Evol. Equat. and Cont. Theo._ , 9 (2020), 535–559.
  * [21] J. Li, M. Yamamoto and J. Zou,  Conditional stability and numerical reconstruction of initial temperature,  _Comm. Pure. Appl. Anal._ , 8 (2009), 361–382.
  * [22] L. Maniar, M. Meyries and R. Schnaubelt,  Null controllability for parabolic equations with dynamic boundary conditions,  _Evol. Equat. and Cont. Theo._ , 6 (2017), 381–407.
  * [23] K. Masood and M. Yousuf,  Solution of the initial inverse problems in the heat equation using the finite difference method with positivity-preserving Padé schemes,  _Numer. Heat Transf. A: Appl._ , 57 (2010), 691–708.
  * [24] K. Masood and F. D. Zaman,  Investigation of the initial inverse problem in the heat equation,  _J. Heat Mass Transf._ , 126 (2004), 294–296.
  * [25] L. E. Payne,  _Improperly Posed Problems in Partial Differential Equations_ ,  SIAM, Philadelphia, 1975.
  * [26] S. S. Pereverzyev, R. Pinnau and N. Siedow,  Initial temperature reconstruction for nonlinear heat equation: application to a coupled radiative–conductive heat transfer problem,  _Inverse Probl. Sci. Eng._ 16 (2008), 55–67.
  * [27] W. Wang, M. Yamamoto and B. Han,  Two-dimensional parabolic inverse source problem with final overdetermination in reproducing kernel space.  _Chinese Ann. Math. Ser. B_ , 35 (2014), 469–482.
  * [28] L. Yang, Z-C Deng and Y. C. Hon,  Simultaneous identification of unknown initial temperature and heat source,  _Dyn. Syst. Appl._ , 25, (2016), 583–602.
  * [29] L. Yang, J-N Yu, G-W Luo and Z-C Deng,  Numerical identification of source terms for a two dimensional heat conduction problem in polar coordinate system,  _Appl. Math. Model._ , 37, (2013), 939–957.