# Identification of source terms in heat equation with dynamic boundary
conditions

E. M. Ait Ben Hassi  E. M. Ait Ben hassi, Cadi Ayyad University, Faculty of
Sciences Semlalia, LMDP, UMMISCO (IRD-UPMC), B.P. 2390, Marrakesh, Morocco
<EMAIL_ADDRESS>,  S. E. Chorfi  S. E. Chorfi, Cadi Ayyad University,
Faculty of Sciences Semlalia, LMDP, UMMISCO (IRD-UPMC), B.P. 2390, Marrakesh,
Morocco<EMAIL_ADDRESS>and  L. Maniar  L. Maniar, Cadi Ayyad University,
Faculty of Sciences Semlalia, LMDP, UMMISCO (IRD-UPMC), B.P. 2390, Marrakesh,
Morocco<EMAIL_ADDRESS>
###### Abstract.

We study an inverse parabolic problem of identifying two source terms in heat
equation with dynamic boundary conditions from a final time overdetermination
data. Using a weak solution approach by Hasanov, the associated cost
functional is analyzed, especially a gradient formula of the functional is
proved and given in terms of the solution of an adjoint problem. Next, the
existence and uniqueness of a quasi-solution are also investigated. Finally,
the numerical reconstruction of some heat sources in a 1-D equation is
presented to show the efficiency of the proposed algorithm.

###### Key words and phrases:

inverse source problem, parabolic problem, quasi-solution, adjoint problem,
Fréchet gradient, Lipschitz continuity.

###### Key words and phrases:

inverse source problem, parabolic problem, quasi-solution, adjoint problem,
Fréchet gradient, Lipschitz continuity.

###### 2020 Mathematics Subject Classification:

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

##  1\. Introduction

In this paper, we are interested in an inverse parabolic problem. It consists
of identifying two source terms in a heat equation with dynamic boundary
conditions, from a noisy measurement of the temperature at final time.

Let $T>0$ be a fixed final time and let $\Omega\subset\mathbb{R}^{N}$ be a
bounded domain ($N\geq 2$ is an integer) with boundary $\Gamma=\partial\Omega$
of class $C^{2}$. We denote $\Omega_{T}=(0,T)\times\Omega$ and
$\Gamma_{T}=(0,T)\times\Gamma$. We consider the following heat equation

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

for initial data $Y_{0}:=(y_{0},y_{0,\Gamma})\in L^{2}(\Omega)\times
L^{2}(\Gamma)$ and heat sources $F\in L^{2}(\Omega_{T})$ and $G\in
L^{2}(\Gamma_{T})$. The diffusion coefficients are positive constants
$d,\gamma>0$, and the spatial potentials are such that $a\in
L^{\infty}(\Omega)$ and $b\in L^{\infty}(\Gamma)$. $\Delta=\Delta_{x}$ denotes
the standard Laplace operator with respect to the space variable. By
$y_{\mid\Gamma},$ one denotes the trace of $y$, while the normal derivative is
denoted by $\partial_{\nu}y:=(\nabla y\cdot\nu)_{|\Gamma}$, where $\nu(x)$
stands for the unit outward normal vector to $\Gamma$ at $x$. The tangential
gradient $\nabla_{\Gamma}y$ is given by $\nabla_{\Gamma}y=\nabla
y-\left(\partial_{\nu}y\right)\nu$. Let $\mathrm{g}$ be the natural Riemannian
metric on $\Gamma$ inherited from $\mathbb{R}^{N}$. The Laplace-Beltrami
operator $\Delta_{\Gamma}$ (with respect to $\mathrm{g}$) is given locally by

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

where $\mathrm{g}=\left(\mathrm{g}_{ij}\right)$ is the metric tensor
corresponding to $\mathrm{g}$, $\mathrm{g}^{-1}=\left(\mathrm{g}^{ij}\right)$
its inverse and $|\mathrm{g}|=\det\left(\mathrm{g}_{ij}\right)$. In the
sequel, we mainly use the following 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.2)

where $\langle\cdot,\cdot\rangle_{\Gamma}$ is the Riemannian inner product of
tangential vectors on $\Gamma$.

Parabolic equations with dynamic boundary conditions have received a lot of
attention of many researchers in the last years [3, 9, 21, 22, 26, 27], since
they appear in several fields of applications including chemical engineering
such as chemical reactors, the chemistry of colloids and special flows in
hydrodynamics [30]. The scope of applications also includes heat transfer
problems where the diffusion takes place between a solid and a fluid in motion
[25]. We refer to the seminal paper [11], where the physical derivation and
interpretation of dynamic boundary conditions are discussed.

A variety of analytical and numerical techniques have been proposed for
solving inverse problems in partial differential equations (PDEs) with static
boundary conditions (Dirichlet, Neumann or Robin conditions); see, for
instance, [4, 8, 12, 18]. A well-known theoretical technique is the Carleman
estimate, which is a priori weighted inequality estimating the solution of a
PDE and its derivatives via the associated differential operator. Such an
estimate was introduced, for the first time, in the study of multidimensional
inverse problems by Bukhgeim and Klibanov [4], and became a powerful tool to
establish uniqueness and stability results. Recently in [2], the authors have
applied this method in the study of an inverse source problem for general
parabolic equations with dynamic boundary conditions from interior
measurements. They established a Lipschitz stability estimate for the source
terms. In [1], we have considered an inverse problem of radiative potentials
and initial temperatures for the same equation. We have proven a Lipschitz
stability estimate for the potentials and then obtained a logarithmic
stability result for the initial conditions by a logarithmic convexity method.
It is worthwhile to mention that Carleman estimates are also important in the
numerical study of inverse problems [23].

In the context of numerical studies, Hasanov [12] developed a weak solution
approach to study the simultaneous determination of source terms in heat
equation with static boundary conditions from the final overdetermination,
Dirichlet or Neumann types output measured data [14, 15, 20]. The underlying
method relies on reformulating a quasi-solution problem as a minimization
problem of Tikhonov functional, combined with an adjoint problem approach
introduced by DuChateau [5, 6]. This approach provides a monotone iteration
scheme to reconstruct unknown parameters in ill-posed problems, which implies
fast numerical results. It has also been successfully applied to source terms
in the cantilevered beam equation [13], Lotka Volterra system [10], and
recently for thermal conductivity and radiative potential in heat equation
from Dirichlet and Neumann boundary measured outputs [16]. We refer the
interested readers to the recent book by Hasanov and Romanov [17], which
presents a systematic study of mathematical and numerical methods used in
inverse problems within this framework.

In the present work, we extend the previous technique for the determination of
two source terms acting in both, in the domain and on the boundary, in heat
equation with dynamic boundary conditions. In this context, we mention [28]
for the identification of a time-dependent source term from the knowledge of a
space average using the backward Euler method. Also, in [19], the author deals
with a time-dependent source from an integral overdetermination condition in a
1-D heat equation by using the generalized Fourier method. These two papers
dealt with a basic dynamic boundary condition that only contains time and
normal derivatives on the boundary (known as Wentzell/Ventcel boundary
condition). To the best of our knowledge, the study of such problem by a weak
solution approach has not been addressed in the literature, and only equations
with classical boundary conditions (Dirichlet, Neumann and Robin) are
considered. Here, we deal with a multidimensional heat system of two equations
coupled through the boundary via the normal derivative, and which contains a
surface diffusion on the boundary. Thus, the present problem requires a
careful analysis and suitable combination between the two equations in order
to obtain the desired results. We also study numerically the reconstruction
for a heat source in the 1-D equation from a noisy terminal data.

The rest of the paper is organized as follows: in Section 2, we recall the
well-posedness and regularity results of the system (1.1). In Section 3, we
study the minimization problem associated to quasi-solutions. We also infer an
explicit gradient formula for the Tikhonov functional via the solution of an
appropriate adjoint problem. Then, we establish the Lipschitz continuity of
the gradient of the cost functional. In Section 4, we discuss the existence
and uniqueness of a quasi-solution. Finally, in Section 5, the gradient
formula is implemented via Landweber scheme for the numerical reconstruction
of an unknown source term in 1-D equation. Section 6 is devoted to some
conclusions.

##  2\. Wellposedness and regularity of the solution

We recall some results on the wellposedness and regularity of the solution to
system (1.1) needed in the sequel. The reader can refer to [26] for detailed
proofs.

We first introduce the following real spaces

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

$\mathbb{L}^{2}$ and $\mathbb{L}^{2}_{T}$ are Hilbert spaces equipped with the
scalar products given by

$\displaystyle\langle(y,y_{\Gamma}),(z,z_{\Gamma})\rangle_{\mathbb{L}^{2}}$
$\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}^{2}_{T}}$
$\displaystyle=\langle y,z\rangle_{L^{2}(\Omega_{T})}+\langle
y_{\Gamma},z_{\Gamma}\rangle_{L^{2}(\Gamma_{T})},$

respectively, where we denoted the Lebesgue measure on $\Omega$ by
$\mathrm{d}x$ and the surface measure on $\Gamma$ by $\mathrm{d}S$. We also
define the 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,$

equipped with the standard product norm.

For the regularity of the solution, we consider the following spaces

$\mathbb{E}_{1}(t_{0},t_{1}):=H^{1}\left(t_{0},t_{1};\mathbb{L}^{2}\right)\cap
L^{2}\left(t_{0},t_{1};\mathbb{H}^{2}\right)\text{ for }t_{1}>t_{0}\text{ in
}\mathbb{R},$
$\mathbb{E}_{2}(t_{0},t_{1}):=H^{1}\left(t_{0},t_{1};\mathbb{H}^{2}\right)\cap
H^{2}\left(t_{0},t_{1};\mathbb{L}^{2}\right)\text{ for }t_{1}>t_{0}\text{ in
}\mathbb{R}.$

In particular,

$\mathbb{E}_{1}:=\mathbb{E}_{1}(0,T)\;\text{ and
}\;\mathbb{E}_{2}:=\mathbb{E}_{2}(0,T).$

We rewrite (1.1) in the following abstract form

$\text{(ACP) }\;\begin{cases}\partial_{t}Y=\mathcal{A}Y+\mathcal{F},\quad
0<t\leq T,\\\ Y(0)=Y_{0}:=(y_{0},y_{0,\Gamma}),\end{cases}$

where $Y:=(y,y_{\Gamma})$, $\mathcal{F}=(F,G)$ and the linear operator
$\mathcal{A}\colon
D(\mathcal{A})\subset\mathbb{L}^{2}\longrightarrow\mathbb{L}^{2}$ is given 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}$ generates an analytic $C_{0}$-semigroup
$\left(\mathrm{e}^{t\mathcal{A}}\right)_{t\geq 0}$ on $\mathbb{L}^{2}$, see
[26] for more details.

In the sequel, we will adopt the following notions for the solutions to system
(1.1).

######  Definition 2.1.

Let $Y_{0}\in\mathbb{L}^{2}$ and $(F,G)\in\mathbb{L}^{2}_{T}$.

  1. (a)

A strong solution of (1.1) is a function $Y:=(y,y_{\Gamma})\in\mathbb{E}_{1}$
fulfilling (1.1) in $L^{2}(0,T;\mathbb{L}^{2})$.

  2. (b)

A mild solution of (1.1) is a function $Y:=(y,y_{\Gamma})\in
C([0,T];\mathbb{L}^{2})$ satisfying, for $t\in[0,T]$,

$Y(t,\cdot)=\mathrm{e}^{t\mathcal{A}}Y_{0}+\int_{0}^{t}\mathrm{e}^{(t-\tau)\mathcal{A}}[F(\tau,\cdot),G(\tau,\cdot)]\,\mathrm{d}\tau.$

  3. (c)

A distributional (weak) solution of (1.1) is a function $Y:=(y,y_{\Gamma})\in
L^{2}(0,T;\mathbb{L}^{2})$ such that for all
$(\varphi,\varphi_{\Gamma})\in\mathbb{E}_{1}$ with
$\varphi(T,\cdot)=\varphi_{\Gamma}(T,\cdot)=0$ we have

$\displaystyle\int_{\Omega_{T}}y(-\partial_{t}\varphi-d\Delta\varphi+a\varphi)\,\mathrm{d}x\,\mathrm{d}t+\int_{\Gamma_{T}}y_{\Gamma}(-\partial_{t}\varphi_{\Gamma}-\gamma\Delta_{\Gamma}\varphi_{\Gamma}+d\partial_{\nu}\varphi+b\varphi_{\Gamma})\,\mathrm{d}S\,\mathrm{d}t$
$\displaystyle=\int_{\Omega_{T}}F\varphi\,\mathrm{d}x\,\mathrm{d}t+\int_{\Gamma_{T}}G\varphi_{\Gamma}\,\mathrm{d}S\,\mathrm{d}t+\int_{\Omega}y_{0}\varphi(0,\cdot)\,\mathrm{d}x+\int_{\Gamma}y_{0,\Gamma}\varphi_{\Gamma}(0,\cdot)\,\mathrm{d}S.$

The following proposition shows the $\mathbb{L}^{2}$-regularity for the system
(1.1) and highlights the connections between different types of solutions. For
the proof, see [26, Propositions 2.4 and 2.5].

######  Proposition 2.2.

Let $F\in L^{2}(\Omega_{T})$ and $G\in L^{2}(\Gamma_{T})$.

  1. (i)

For all $Y_{0}:=(y_{0},y_{0,\Gamma})\in\mathbb{H}^{1}$, there exists a unique
strong solution of (1.1) such that $Y:=(y,y_{\Gamma})\in\mathbb{E}_{1}$.
Moreover, there exists a positive constant
$C=C(\|a\|_{\infty},\|b\|_{\infty})$ such that

$\|Y\|_{\mathbb{E}_{1}}\leq
C\left(\|F\|_{L^{2}(\Omega_{T})}+\|G\|_{L^{2}(\Gamma_{T})}+\|Y_{0}\|_{\mathbb{H}^{1}}\right).$
(2.1)

  2. (ii)

For all $Y_{0}:=(y_{0},y_{0,\Gamma})\in\mathbb{L}^{2}$, there exists a unique
mild solution of (1.1) $Y:=(y,y_{\Gamma})\in C([0,T];\mathbb{L}^{2})$ such
that for all $\tau\in(0,T)$,

$Y\in\mathbb{E}_{1}(\tau,T):=H^{1}\left(\tau,T;\mathbb{L}^{2}\right)\cap
L^{2}\left(\tau,T;\mathbb{H}^{2}\right).$

Moreover, if $\mathcal{F}=(F,G)\in H^{1}(0,T;\mathbb{L}^{2})$, then for all
$\tau\in(0,T)$, we have

$Y\in\mathbb{E}_{2}(\tau,T):=H^{1}\left(\tau,T;\mathbb{H}^{2}\right)\cap
H^{2}\left(\tau,T;\mathbb{L}^{2}\right),$

with initial data $Y(\tau)$.

  3. (iii)

A function $Y$ is a distributional solution of (1.1) if and only if it is a
mild solution.

##  3\. Fréchet differentiability and gradient formula of the cost functional

In this section, we consider the following inverse source problem.

Inverse Source Problem (ISP). A couple of source terms
$(F,G)\in\mathbb{L}^{2}_{T}$ in (1.1) is unknown and needs to be recovered
from the final temperature at $T$, namely,

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

which is not necessarily smooth due to the numerical noise.

Let $Y(t,\cdot,\mathcal{F})$ be the mild solution of (1.1) corresponding to
the source terms $\mathcal{F}=(F,G)\in\mathbb{L}^{2}_{T}$. We introduce the
input-output operator
$\Psi\colon\mathbb{L}^{2}_{T}\longrightarrow\mathbb{L}^{2}$ defined as follows

$(\Psi\mathcal{F})(\cdot)=Y_{T}(\cdot):=Y(T,\cdot,\mathcal{F})\qquad\text{ on
}\;\Omega\times\Gamma.$

Then the inverse source problem (ISP) can be reformulated as the following
operator equation

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

The following lemma is needed in the sequel.

######  Lemma 3.1.

Let $Y_{0}\in\mathbb{L}^{2}$ and let $Y$ be the mild solution of (1.1)
corresponding to $\mathcal{F}=(F,G)$. Then the solution map
$\mathcal{F}\longmapsto Y$ is continuous from
$H^{1}\left(0,T;\mathbb{L}^{2}\right)$ to
$C\left([0,T];\mathbb{L}^{2}\right)\cap L^{2}\left(0,T;\mathbb{H}^{1}\right)$.

###### Proof.

Let $Y_{0}\in\mathbb{H}^{2}$ and $\delta\mathcal{F}$ be a small variation of
$\mathcal{F}$ such that $\mathcal{F}+\delta\mathcal{F}\in\mathcal{U}$.
Consider $\delta Y:=Y^{\delta}-Y$, where $Y^{\delta}$ is the mild solution of
(1.1) corresponding to $\mathcal{F}^{\delta}:=\mathcal{F}+\delta\mathcal{F}$.
Then $\delta Y\in C^{1}\left([0,T],\mathbb{L}^{2}\right)\cap
C\left(0,T;\mathbb{H}^{2}\right)$ and satisfies the following system

$\displaystyle\begin{aligned} &\partial_{t}(\delta y)-d\Delta(\delta
y)+a(x)(\delta y)=\delta F(t,x),&\qquad\text{in }\Omega_{T},\\\
&\partial_{t}(\delta y_{\Gamma})-\gamma\Delta_{\Gamma}(\delta
y_{\Gamma})+d\partial_{\nu}y+b(x)(\delta y_{\Gamma})=\delta
G(t,x),&\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}=(0,0),&\qquad\Omega\times\Gamma,\end{aligned}$ (3.2)

where

$\displaystyle\delta y(T,\cdot,\mathcal{F})$
$\displaystyle=y(T,\cdot,\mathcal{F}+\delta\mathcal{F})-y(T,\cdot,\mathcal{F}),$
$\displaystyle\delta y_{\Gamma}(T,\cdot,\mathcal{F})$
$\displaystyle=y_{\Gamma}(T,\cdot,\mathcal{F}+\delta\mathcal{F})-y_{\Gamma}(T,\cdot,\mathcal{F}).$

Multiplying $\eqref{0deq1to4}_{1}$ by $\delta y$, $\eqref{0deq1to4}_{2}$ by
$\delta y_{\Gamma}$ and using Green and surface divergence formulae, we obtain

$\displaystyle\frac{1}{2}\partial_{t}\left(\int_{\Omega}|\delta
y|^{2}\,\mathrm{d}x\right)+d\int_{\Omega}|\nabla(\delta
y)|^{2}\,\mathrm{d}x-\frac{1}{2}\int_{\Gamma}d\partial_{\nu}(|\delta
y|^{2})\,\mathrm{d}S+\int_{\Omega}a|\delta y|^{2}\,\mathrm{d}x$
$\displaystyle\qquad=\int_{\Omega}\delta F\delta y\,\mathrm{d}x,$
$\displaystyle\frac{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+\frac{1}{2}\int_{\Gamma}d\partial_{\nu}(|\delta
y|^{2})\,\mathrm{d}S+\int_{\Gamma}b|\delta y_{\Gamma}|^{2}\,\mathrm{d}S$
$\displaystyle\qquad=\int_{\Gamma}\delta G\delta y_{\Gamma}\,\mathrm{d}S.$

Adding up the last two equalities, we arrive at

$\displaystyle\frac{1}{2}\partial_{t}\left(\int_{\Omega}|\delta
y|^{2}\,\mathrm{d}x+\int_{\Gamma}|\delta
y_{\Gamma}|^{2}\,\mathrm{d}S\right)+d\int_{\Omega}|\nabla(\delta
y)|^{2}\,\mathrm{d}x+\gamma\int_{\Gamma}|\nabla_{\Gamma}(\delta
y_{\Gamma})|^{2}\,\mathrm{d}S$ $\displaystyle\qquad+\int_{\Omega}a|\delta
y|^{2}\,\mathrm{d}x+\int_{\Gamma}b|\delta y_{\Gamma}|^{2}\,\mathrm{d}S$
$\displaystyle=\int_{\Omega}\delta F\delta y\,\mathrm{d}x+\int_{\Gamma}\delta
G\delta y_{\Gamma}\,\mathrm{d}S.$

Then, using the Cauchy-Schwarz inequality, we have

$\displaystyle\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\|\delta
Y(t)\|^{2}_{\mathbb{L}^{2}}+\min(d,\gamma)\|(\nabla
y(t),\nabla_{\Gamma}y_{\Gamma}(t))\|^{2}_{\mathbb{L}^{2}}$
$\displaystyle\qquad\leq\langle\delta
Y,\delta\mathcal{F}\rangle_{\mathbb{L}^{2}}+\max(\|a\|_{\infty},\|b\|_{\infty})\|\delta
Y(t)\|^{2}_{\mathbb{L}^{2}}$ $\displaystyle\qquad\leq\|\delta
Y(t)\|_{\mathbb{L}^{2}}\|\delta\mathcal{F}(t)\|_{\mathbb{L}^{2}}+\max(\|a\|_{\infty},\|b\|_{\infty})\|\delta
Y(t)\|^{2}_{\mathbb{L}^{2}}$ $\displaystyle\qquad\leq\frac{1}{2}\|\delta
Y(t)\|^{2}_{\mathbb{L}^{2}}+\frac{1}{2}\|\delta\mathcal{F}(t)\|^{2}_{\mathbb{L}^{2}}+\max(\|a\|_{\infty},\|b\|_{\infty})\|\delta
Y(t)\|^{2}_{\mathbb{L}^{2}},$ (3.3)

where $\mathcal{F}=(F,G)$. Consequently,

$\frac{\mathrm{d}}{\mathrm{d}t}\|\delta
Y(t)\|^{2}_{\mathbb{L}^{2}}\leq(1+2\max(\|a\|_{\infty},\|b\|_{\infty}))\|\delta
Y(t)\|^{2}_{\mathbb{L}^{2}}+\|\delta\mathcal{F}(t)\|^{2}_{\mathbb{L}^{2}}.$

By Gronwall inequality, we deduce

$\displaystyle\|\delta Y(t)\|^{2}_{\mathbb{L}^{2}}$
$\displaystyle\leq\mathrm{e}^{(1+2\max(\|a\|_{\infty},\|b\|_{\infty}))T}\left(\|\delta
Y(0)\|^{2}_{\mathbb{L}^{2}}+\|\delta\mathcal{F}\|^{2}_{\mathbb{L}^{2}_{T}}\right)$
$\displaystyle=\mathrm{e}^{(1+2\max(\|a\|_{\infty},\|b\|_{\infty}))T}\|\delta\mathcal{F}\|^{2}_{\mathbb{L}^{2}_{T}}$
(3.4)

for every $0\leq t\leq T$. On the other hand, from (3.3) and (3.4), we have

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\|\delta
Y(t)\|^{2}_{\mathbb{L}^{2}}+2\min(d,\gamma)\|(\nabla
y(t),\nabla_{\Gamma}y_{\Gamma}(t))\|^{2}_{\mathbb{L}^{2}}$
$\displaystyle\qquad\leq\left(1+2\max(\|a\|_{\infty},\|b\|_{\infty})\right)\|\delta
Y(t)\|^{2}_{\mathbb{L}^{2}}+\|\delta\mathcal{F}(t)\|^{2}_{\mathbb{L}^{2}}.$
$\displaystyle\qquad\leq\left(1+2\max(\|a\|_{\infty},\|b\|_{\infty})\right)\mathrm{e}^{(1+2\max(\|a\|_{\infty},\|b\|_{\infty}))T}\|\delta\mathcal{F}\|^{2}_{\mathbb{L}^{2}_{T}}+\|\delta\mathcal{F}(t)\|^{2}_{\mathbb{L}^{2}}.$

Integrating between $0$ and $T$ and using the fact that $\delta Y(0)=0$, we
deduce

$\displaystyle\|\delta
Y(T)\|^{2}_{\mathbb{L}^{2}}+2\min(d,\gamma)\int_{0}^{T}\|(\nabla
y(t),\nabla_{\Gamma}y_{\Gamma}(t))\|^{2}_{\mathbb{L}^{2}}\,\mathrm{d}t$
$\displaystyle\qquad\leq\left(1+2\max(\|a\|_{\infty},\|b\|_{\infty})\right)T\mathrm{e}^{(1+2\max(\|a\|_{\infty},\|b\|_{\infty}))T}\|\delta\mathcal{F}\|^{2}_{\mathbb{L}^{2}_{T}}+\|\delta\mathcal{F}\|^{2}_{\mathbb{L}^{2}_{T}}.$

Then we arrive at

$\displaystyle\int_{0}^{T}\|(\nabla
y(t),\nabla_{\Gamma}y_{\Gamma}(t))\|^{2}_{\mathbb{L}^{2}}\,\mathrm{d}t\leq
C_{T}\|\delta\mathcal{F}\|^{2}_{\mathbb{L}^{2}_{T}},$

with

$C_{T}:=\frac{(1+2\max(\|a\|_{\infty},\|b\|_{\infty}))T\mathrm{e}^{(1+2\max(\|a\|_{\infty},\|b\|_{\infty}))T}+1}{2\min(d,\gamma)}.$

The latter inequality and (3.4) imply that

$\displaystyle\sup_{t\in[0,T]}\|\delta
Y(t)\|^{2}_{\mathbb{L}^{2}}+\int_{0}^{T}\|\delta
Y(t)\|^{2}_{\mathbb{H}^{1}}\,\mathrm{d}t\leq
C\|\delta\mathcal{F}\|^{2}_{\mathbb{L}^{2}_{T}}$

for some generic constant $C=C(T,\|a\|_{\infty},\|b\|_{\infty},d,\gamma)>0$.
Hence,

$\displaystyle\|\delta Y\|^{2}_{C([0,T];\mathbb{L}^{2})}+\|\delta
Y\|^{2}_{L^{2}(0,T;\mathbb{H}^{1})}\leq
C\|\delta\mathcal{F}\|^{2}_{H^{1}(0,T;\mathbb{L}^{2})}.$

Since $\mathbb{H}^{2}$ is dense in $\mathbb{L}^{2}$, the same inequality holds
for any $Y_{0}\in\mathbb{L}^{2}$. This ends the proof. ∎

The following lemma highlights the compactness of the input-output operator
$\Psi\colon\mathbb{L}^{2}_{T}\longrightarrow\mathbb{L}^{2}$. By linearity, we
may assume here that $Y_{0}:=(y_{0},y_{0,\Gamma})=(0,0)$.

######  Lemma 3.2.

The input-output operator

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

is compact.

###### Proof.

Consider $(\mathcal{F}_{n})_{n\in\mathbb{N}}$ a bounded sequence in
$\mathbb{L}^{2}_{T}$. By estimate (2.1), the sequence of associated strong
solutions $(Y(\cdot,\cdot,\mathcal{F}_{n}))_{n\in\mathbb{N}}$ of (1.1) is
bounded in $\mathbb{E}_{1}$. Since the trace space of $\mathbb{E}_{1}$ at
$t=T$ equals $\mathbb{H}^{1}$ (see [26, Proposition 2.2]), the sequence
$(Y_{n})_{n\in\mathbb{N}}$ defined by $Y_{n}:=Y(T,\cdot,\mathcal{F}_{n})$ is
bounded in $\mathbb{H}^{1}$. By the compact embedding
$\mathbb{H}^{1}\hookrightarrow\mathbb{L}^{2}$ (see [26]), there exists a
subsequence of $(Y_{n})_{n\in\mathbb{N}}$ which converges in $\mathbb{L}^{2}$.
Thus, the input-output operator $\Psi$ is compact. ∎

It follows from Lemma 3.2 that the inverse problem (3.1) is ill-posed (in the
sense of Hadamard). A quasi-solution $\mathcal{F}_{*}$ of the ill-posed
problem (3.1) is defined as a solution of the following minimization problem

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

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$,
and

$\mathcal{U}:=\\{\mathcal{F}=(F,G)\in
H^{1}(0,T;\mathbb{L}^{2})\colon\|\mathcal{F}\|_{H^{1}(0,T;\mathbb{L}^{2})}\leq
R\\}$

is the set of admissible source terms. Clearly, $\mathcal{U}$ is a bounded,
closed and convex subset of $H^{1}(0,T;\mathbb{L}^{2})$, for any fixed $R>0$.

Due to the ill-posedness of (3.1), we usually use a Tikhonov regularization
approach and consider instead the following regularized functional

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

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

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{F}\in\mathcal{U}$ is given by

$\mathcal{J}^{\prime}(\mathcal{F})=\Phi,$ (3.7)

where $\Phi(t,\cdot,\mathcal{F})=(\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{F})-y_{T}^{\delta},&&\qquad\text{in
}\Omega,\\\
&\varphi_{\Gamma}\rvert_{t=T}=y_{\Gamma}(T,\cdot,\mathcal{F})-y_{T,\Gamma}^{\delta},&&\qquad\text{on
}\Gamma.\end{aligned}$ (3.8)

###### Proof.

We assume that $\mathcal{F},\mathcal{F}+\delta\mathcal{F}\in\mathcal{U}$. Let
us calculate the difference

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

We have

$\displaystyle\delta\mathcal{J}(\mathcal{F})$
$\displaystyle=\frac{1}{2}\|Y(T,\cdot,\mathcal{F}+\delta\mathcal{F})-Y_{T}^{\delta}\|_{\mathbb{L}^{2}}^{2}-\frac{1}{2}\|Y(T,\cdot,\mathcal{F})-Y_{T}^{\delta}\|_{\mathbb{L}^{2}}^{2},$
$\displaystyle=\frac{1}{2}\left(\|y(T,\cdot,\mathcal{F}+\delta\mathcal{F})-y_{T}^{\delta}\|_{L^{2}(\Omega)}^{2}+\|y_{\Gamma}(T,\cdot,\mathcal{F}+\delta\mathcal{F})-y_{T,\Gamma}^{\delta}\|_{L^{2}(\Gamma)}^{2}\right)$
$\displaystyle\qquad-\frac{1}{2}\left(\|y(T,\cdot,\mathcal{F})-y_{T}^{\delta}\|_{L^{2}(\Omega)}^{2}+\|y_{\Gamma}(T,\cdot,\mathcal{F})-y_{T,\Gamma}^{\delta}\|_{L^{2}(\Gamma)}^{2}\right)$
$\displaystyle=\frac{1}{2}\int_{\Omega}\left[(y(T,x,\mathcal{F}+\delta\mathcal{F})-y_{T}^{\delta}(x))^{2}-(y(T,x,\mathcal{F})-y_{T}^{\delta}(x))^{2}\right]\,\mathrm{d}x$
$\displaystyle\quad+\frac{1}{2}\int_{\Gamma}\left[(y_{\Gamma}(T,x,\mathcal{F}+\delta\mathcal{F})-y_{T,\Gamma}^{\delta}(x))^{2}-(y_{\Gamma}(T,x,\mathcal{F})-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{F})$
$\displaystyle=\int_{\Omega}(y(T,x,\mathcal{F})-y_{T}^{\delta}(x))\delta
y(T,x,\mathcal{F})\,\mathrm{d}x+\frac{1}{2}\int_{\Omega}[\delta
y(T,x,\mathcal{F})]^{2}\,\mathrm{d}x$ (3.9)
$\displaystyle+\int_{\Gamma}(y_{\Gamma}(T,x,\mathcal{F})-y_{T,\Gamma}^{\delta}(x))\delta
y_{\Gamma}(T,x,\mathcal{F})\,\mathrm{d}S+\frac{1}{2}\int_{\Gamma}[\delta
y_{\Gamma}(T,x,\mathcal{F})]^{2}\,\mathrm{d}S,$ (3.10)

where

$\displaystyle\delta y(T,\cdot,\mathcal{F})$
$\displaystyle=y(T,\cdot,\mathcal{F}+\delta\mathcal{F})-y(T,\cdot,\mathcal{F}),$
$\displaystyle\delta y_{\Gamma}(T,\cdot,\mathcal{F})$
$\displaystyle=y_{\Gamma}(T,\cdot,\mathcal{F}+\delta\mathcal{F})-y_{\Gamma}(T,\cdot,\mathcal{F}).$

By linearity of the systems, $\delta Y=(\delta y,\delta y_{\Gamma})$ is the
mild solution of the following system

$\displaystyle\begin{aligned} &\partial_{t}(\delta y)-d\Delta(\delta
y)+a(x)(\delta y)=\delta F(t,x),&\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})=\delta
G(t,x),&\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}=(0,0),&\qquad\Omega\times\Gamma.\end{aligned}$ (3.11)

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

$\displaystyle\int_{\Omega}(y(T,x,\mathcal{F})-y_{T}^{\delta}(x))\delta
y(T,x,\mathcal{F})\,\mathrm{d}x=\int_{\Omega}\varphi(T,x,\mathcal{F})\delta
y(T,x,\mathcal{F})\,\mathrm{d}x$
$\displaystyle=\bigintss_{\Omega}\left[\int_{0}^{T}\partial_{t}(\varphi(t,x,\mathcal{F})\delta
y(t,x,\mathcal{F}))\,\mathrm{d}t\right]\,\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$
$\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_{T}}\delta
F(t,x)\,\varphi(t,x)\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle=\int_{\Omega_{T}}-d[(\Delta\varphi)\delta y-\Delta(\delta
y)\varphi]\,\mathrm{d}x\,\mathrm{d}t+\int_{\Omega_{T}}\delta
F(t,x)\,\varphi(t,x)\,\mathrm{d}x\,\mathrm{d}t$
$\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_{T}}\delta
F(t,x)\,\varphi(t,x)\,\mathrm{d}x\,\mathrm{d}t.$ (3.12)

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

$\displaystyle\int_{\Gamma}(y_{\Gamma}(T,x,\mathcal{F})-y_{T,\Gamma}^{\delta}(x))\delta
y_{\Gamma}(T,x,\mathcal{F})\,\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.13)
$\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_{T}}\delta
G(t,x)\,\varphi_{\Gamma}(t,x)\,\mathrm{d}S\,\mathrm{d}t$
$\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_{T}}\delta
G(t,x)\,\varphi_{\Gamma}(t,x)\,\mathrm{d}S\,\mathrm{d}t$ (3.14)

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

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

For the second integrals in the right-hand side of (3.9) and (3.10), the
estimate (3.4) implies

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

This completes the proof. ∎

Next, we prove the Lipschitz continuity of the Fréchet gradient
$\mathcal{J}^{\prime}$, in particular, $\mathcal{J}\in C^{1}(\mathcal{U})$.

######  Lemma 3.4.

Let $\mathcal{F},\delta\mathcal{F}\in\mathcal{U}$. Then the Fréchet gradient
$\mathcal{J}^{\prime}$ is Lipschitz continuous,

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

where the Lipschitz constant $L>0$ depends on $T,\|a\|_{\infty}$ and
$\|b\|_{\infty}$ as follows

$L=\sqrt{2T\mathrm{e}^{(1+4\max(\|a\|_{\infty},\|b\|_{\infty}))T}}.$ (3.15)

###### Proof.

Let $\delta\Phi(t,\cdot,\mathcal{F}):=(\delta\varphi,\delta\varphi_{\Gamma})$
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.16)

Using Proposition 3.3, we have

$\displaystyle\|\mathcal{J}^{\prime}(\mathcal{F}+\delta\mathcal{F})-\mathcal{J}^{\prime}(\mathcal{F})\|^{2}_{\mathbb{L}^{2}_{T}}$
$\displaystyle=\int_{\Omega_{T}}|\delta\varphi|^{2}\,\mathrm{d}x\,\mathrm{d}t+\int_{\Gamma_{T}}|\delta\varphi_{\Gamma}|^{2}\,\mathrm{d}S\,\mathrm{d}t$
$\displaystyle\leq
2\left(\|\delta\varphi\|^{2}_{L^{2}(\Omega_{T})}+\|\delta\varphi_{\Gamma}\|^{2}_{L^{2}(\Gamma_{T})}\right)$
$\displaystyle=2\|\delta\Phi\|_{\mathbb{L}^{2}_{T}}^{2}.$ (3.17)

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

$\displaystyle-\frac{1}{2}\partial_{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}\partial_{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 used integration by parts and the surface divergence formula (1.2).
Adding up the last two equalities, we obtain

$\displaystyle\frac{1}{2}\partial_{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.$

Then

$\displaystyle\frac{1}{2}\partial_{t}\left(\int_{\Omega}|\delta\varphi|^{2}\,\mathrm{d}x+\int_{\Gamma}|\delta\varphi_{\Gamma}|^{2}\,\mathrm{d}S\right)+\max(\|a\|_{\infty},\|b\|_{\infty})\left(\int_{\Omega}|\delta\varphi|^{2}\,\mathrm{d}x+\int_{\Gamma}|\delta\varphi_{\Gamma}|^{2}\,\mathrm{d}S\right)\geq
0.$

This inequality implies that the function $H$ defined by

$H(t)=\mathrm{e}^{2\max(\|a\|_{\infty},\|b\|_{\infty})t}\left(\int_{\Omega}|\delta\varphi|^{2}\,\mathrm{d}x+\int_{\Gamma}|\delta\varphi_{\Gamma}|^{2}\,\mathrm{d}S\right)$

is nondecreasing on $[0,T]$. Using $\eqref{2aeq1to4}_{4}$, it follows, for all
$t\in[0,T]$, that

$\displaystyle\|\delta\Phi\|^{2}_{\mathbb{L}^{2}_{T}}=\int_{\Omega_{T}}|\delta\varphi(t,x,\mathcal{F})|^{2}\,\mathrm{d}x\,\mathrm{d}t+\int_{\Gamma_{T}}|\delta\varphi_{\Gamma}(t,x,\mathcal{F})|^{2}\,\mathrm{d}S\,\mathrm{d}t$
$\displaystyle\leq
T\mathrm{e}^{2\max(\|a\|_{\infty},\|b\|_{\infty})T}\left(\int_{\Omega}|\delta\varphi(T,x,\mathcal{F})|^{2}\,\mathrm{d}x+\int_{\Gamma}|\delta\varphi_{\Gamma}(T,x,\mathcal{F})|^{2}\,\mathrm{d}S\right)$
$\displaystyle\leq
T\mathrm{e}^{2\max(\|a\|_{\infty},\|b\|_{\infty})T}\left(\int_{\Omega}|\delta
y(T,x,\mathcal{F})|^{2}\,\mathrm{d}x+\int_{\Gamma}|\delta
y_{\Gamma}(T,x,\mathcal{F})|^{2}\,\mathrm{d}S\right).$

Using this last inequality and inequality (3.4), we obtain

$\displaystyle 2\|\delta\Phi\|^{2}_{\mathbb{L}^{2}_{T}}$ $\displaystyle\leq
2T\mathrm{e}^{2\max(\|a\|_{\infty},\|b\|_{\infty})T}\mathrm{e}^{(1+2\max(\|a\|_{\infty},\|b\|_{\infty}))T}\|\delta\mathcal{F}\|^{2}_{\mathbb{L}^{2}_{T}}$
$\displaystyle\leq
2T\mathrm{e}^{(1+4\max(\|a\|_{\infty},\|b\|_{\infty}))T}\|\delta\mathcal{F}\|^{2}_{\mathbb{L}^{2}_{T}}.$

With the estimate (3.17) this implies that

$\|\mathcal{J}^{\prime}(\mathcal{F}+\delta\mathcal{F})-\mathcal{J}^{\prime}(\mathcal{F})\|^{2}_{\mathbb{L}^{2}_{T}}\leq
2T\mathrm{e}^{(1+4\max(\|a\|_{\infty},\|b\|_{\infty}))T}\|\delta\mathcal{F}\|^{2}_{\mathbb{L}^{2}_{T}}.$

This yields the desired result. ∎

Next, we consider the Landweber method given by the following iteration

$\mathcal{F}_{k+1}=\mathcal{F}_{k}-\alpha_{k}\mathcal{J}^{\prime}(\mathcal{F}_{k}),\quad
k=0,1,2,\dots,$ (3.18)

where $\mathcal{F}_{0}\in\mathcal{U}$ is a given initial iteration and
$\alpha_{k}$ is a relaxation parameter defined by the minimization problem

$h_{k}(\alpha):=\inf_{\alpha\geq 0}h_{k}(\alpha),\quad
h_{k}(\alpha):=\mathcal{J}\left(\mathcal{F}_{k}-\alpha\mathcal{J}^{\prime}(\mathcal{F}_{k})\right),\quad
k=0,1,2,\dots$

We refer to [7] for a detailed exposition of this method.

######  Remark 3.5.

At this level, some remarks should be made:

  * $\bullet$

The Lipschitz continuity of the gradient $\mathcal{J}^{\prime}$ implies that
the sequence $\left(\mathcal{J}(\mathcal{F}_{k})\right)$ is decreasing, where
$\left(\mathcal{F}_{k}\right)$ is defined by (3.18). This fact yields fast
numerical results; see [17, Lemma 3.4.4] for more details.

  * $\bullet$

In some situations, the choice of the iteration parameter $\alpha_{k}>0$ is
difficult. However, the Lipschitz continuity of $\mathcal{J}^{\prime}$ allows
us to estimate this parameter by the Lipschitz constant $L$ via the estimate

$0<\lambda_{0}\leq\alpha_{k}\leq\frac{2}{L+2\lambda_{1}}$

for arbitrary parameters $\lambda_{0},\lambda_{1}>0$; see [12, Section 3.4.3].

  * $\bullet$

If $\alpha_{k}=\alpha=\mathrm{const}>0$ for all $k$, the optimal value of
$\alpha$ (which corresponds to $\lambda_{0}=\frac{1}{L}$ and
$\lambda_{1}=\frac{L}{2}$) is $\alpha_{*}=\frac{1}{L}$. Hence, if $L$ is
large, the step parameter $\alpha_{*}$ will be small. This fact illustrates
the importance of a sharp Lipschitz constant $L$.

The next lemma follows the same ideas in Corollary 4.1 and Theorem 4.1 in [12]
(one can also refer to [29]). We denote by $\mathcal{U}_{*}$ the set of all
quasi-solutions of (3.5)-(3.6).

######  Lemma 3.6.

Let $(\mathcal{F}_{k})\subset\mathcal{U}$ be the sequence defined by (3.18).
If the iteration parameter $\alpha_{k}=\alpha_{*}$ for all $k$. Then the
following assertions hold.

  1. (i)

The sequence $(\mathcal{J}(\mathcal{F}_{k}))$ is monotone decreasing and
convergent. Moreover,

$\lim_{k\to\infty}\|\mathcal{J}^{\prime}(\mathcal{F}_{k})\|_{\mathbb{L}^{2}_{T}}=0.$

Moreover,

$\|\mathcal{F}_{k+1}-\mathcal{F}_{k}\|_{\mathbb{L}^{2}_{T}}^{2}\leq\frac{2}{L}[\mathcal{J}(\mathcal{F}_{k})-\mathcal{J}(\mathcal{F}_{k+1})],\quad
k=1,2,3,\ldots,$

where $L$ is the Lipschitz constant in (3.15);

  2. (ii)

for any given initial iteration $\mathcal{F}_{0}\in\mathcal{U}$, the sequence
$(\mathcal{F}_{k})$ converges weakly in $\mathbb{L}^{2}$ to a quasi-solution
$\mathcal{F}_{*}\in\mathcal{U}_{*}$ of (ISP). The rate of convergence of
$(\mathcal{J}(\mathcal{F}_{k}))$ can be estimated as follows

$0\leq\mathcal{J}(\mathcal{F}_{k})-\mathcal{J}_{*}\leq
2L\frac{\beta^{2}}{k},\quad k=1,2,3,\ldots,$

where
$\mathcal{J}_{*}:=\lim\limits_{k\to\infty}\mathcal{J}\left(\mathcal{F}_{k}\right)$
and
$\beta:=\sup\\{\left\|\mathcal{F}_{k}-\mathcal{F}_{*}\right\|_{\mathbb{L}^{2}_{T}}\colon\mathcal{F}_{k}\in\mathcal{U},\mathcal{F}_{*}\in\mathcal{U}_{*}\\}.$

##  4\. Existence and Uniqueness of the solution to (ISP)

In the following, we use some tools from the calculus of variations to study
the existence and uniqueness of the solution to (ISP). First, we prove the
following lemma.

######  Lemma 4.1.

For the cost functional $\mathcal{J}\in C^{1}(\mathcal{U})$, the following
formula holds

$\langle\mathcal{J}^{\prime}(\mathcal{F}+\delta\mathcal{F})-\mathcal{J}^{\prime}(\mathcal{F}),\delta\mathcal{F}\rangle_{\mathbb{L}^{2}_{T}}=2\|\delta
Y(T,\cdot,\mathcal{F})\|^{2}_{\mathbb{L}^{2}},\quad\forall\mathcal{F},\delta\mathcal{F}\in\mathcal{U},$
(4.1)

where $\delta Y(T,\cdot,\mathcal{F})$ is the solution of (3.11).

###### Proof.

Let $(\delta\varphi,\delta\varphi_{\Gamma})$ be the solution of (3.16). By the
gradient formula (3.7), we have

$\displaystyle\langle\mathcal{J}^{\prime}(\mathcal{F}+\delta\mathcal{F})-\mathcal{J}^{\prime}(\mathcal{F}),\delta\mathcal{F}\rangle_{\mathbb{L}^{2}_{T}}=\int_{\Omega_{T}}\delta
F\delta\varphi\,\mathrm{d}x\,\mathrm{d}t+\int_{\Gamma_{T}}\delta
G\delta\varphi_{\Gamma}\,\mathrm{d}S\,\mathrm{d}t$
$\displaystyle=\int_{\Omega_{T}}\partial_{t}(\delta
y)\delta\varphi\,\mathrm{d}x\,\mathrm{d}t+\int_{\Gamma_{T}}\partial_{t}(\delta
y_{\Gamma})\delta\varphi_{\Gamma}\,\mathrm{d}S\,\mathrm{d}t+\int_{\Omega_{T}}[-d\Delta(\delta
y)+a(x)(\delta y)]\delta\varphi\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle\qquad+\int_{\Gamma_{T}}[-\gamma\Delta(\delta
y_{\Gamma})+d\partial_{\nu}(\delta y)+b(x)(\delta
y_{\Gamma})]\delta\varphi_{\Gamma}\,\mathrm{d}S\,\mathrm{d}t$
$\displaystyle=\|\delta y(T,\cdot)\|^{2}_{L^{2}(\Omega)}+\|\delta
y_{\Gamma}(T,\cdot)\|^{2}_{L^{2}(\Gamma)}-\int_{\Omega_{T}}\partial_{t}(\delta\varphi)\delta
y\,\mathrm{d}x\,\mathrm{d}t-\int_{\Gamma_{T}}\partial_{t}(\delta\varphi_{\Gamma})\delta
y_{\Gamma}\,\mathrm{d}S\,\mathrm{d}t$
$\displaystyle+\int_{\Omega_{T}}[-d\Delta(\delta\varphi)+a(x)(\delta\varphi)]\delta
y\,\mathrm{d}x\,\mathrm{d}t+\int_{\Gamma_{T}}[-\gamma\Delta(\delta\varphi_{\Gamma})+d\partial_{\nu}(\delta\varphi)+b(x)(\delta\varphi_{\Gamma})]\delta
y_{\Gamma}\,\mathrm{d}S\,\mathrm{d}t$ $\displaystyle=\|\delta
y(T,\cdot)\|^{2}_{L^{2}(\Omega)}+\|\delta
y_{\Gamma}(T,\cdot)\|^{2}_{L^{2}(\Gamma)},$

where we used integration by parts with respect to $t$, the Green formula in
$\Omega$, and the surface divergence formula on $\Gamma$, together with the
system (3.11). ∎

The monotonicity of
$\mathcal{J}^{\prime}\colon\mathcal{U}\rightarrow\mathbb{L}^{2}_{T}$ in Lemma
4.1 implies the convexity of the functional $\mathcal{J}$. As a direct
consequence of Lemma 3.1 and [32, Theorem 25.C], we have the following
existence result.

######  Corollary 4.2.

The cost functional $\mathcal{J}$ is continuous and convex on $\mathcal{U}$.
There exists then a minimizer $\hat{\mathcal{F}}\in\mathcal{U}$ such that

$\mathcal{J}\left(\hat{\mathcal{F}}\right)=\min_{\mathcal{F}\in\mathcal{U}}\mathcal{J}(\mathcal{F}).$

Since the strict convexity of $\mathcal{J}$ is characterized by the strict
monotonicity of $\mathcal{J}^{\prime}$, the equality (4.1) yields a sufficient
condition for uniqueness.

######  Lemma 4.3.

If the positivity condition

$\int_{\Omega}(\delta
y(T,x;\mathcal{F}))^{2}\,\mathrm{d}x+\int_{\Gamma}(\delta
y_{\Gamma}(T,x;\mathcal{F}))^{2}\,\mathrm{d}S>0,\qquad\forall\mathcal{F}\in\mathcal{U},$

holds, then the problem (ISP) admits at most one solution.

The above lemma is a consequence of the uniqueness theorem for strictly convex
functionals defined on convex sets (see, e.g.,[32, Corollary 25.15]).

##  5\. Numerical simulation for 1-D internal heat source

In this section, we analyze the numerical reconstruction for an unknown source
term that depends only on the space variable. More precisely, we consider the
reconstruction of $f(x)$ in the following 1-D heat equation with dynamic
boundary conditions

$\displaystyle\begin{aligned}
&y_{t}(t,x)-y_{xx}(t,x)=f(x)r(t,x),&&\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)=y_{0}(x),\quad(y(0,0),y(0,\ell))=(a,b),&&\qquad
x\in(0,\ell),\end{aligned}$ (5.1)

where $T>0$ is a final time, $\ell>0$ and $\left(y_{0},a,b\right)\in
L^{2}(0,\ell)\times\mathbb{R}^{2}$ is an initial condition. The function $r\in
C^{1}\left([0,T];C\left([0,\ell]\right)\right)$ is assumed to be known.

The 1-D heat equation with dynamic boundary conditions has attracted special
attention as a model of heat conduction problems for a metal bar of length
$\ell$. We mention [24], where the authors have studied the approximate
controllability problem from the boundary.

We shall apply the Landweber iteration scheme discussed in a previous section
to system (5.1). Let $Y(t,x,f):=\left(y(t,x),y(t,0),y(t,\ell)\right)$ be the
solution of (5.1). The input-output operator $\Psi\colon
L^{2}(0,\ell)\longrightarrow L^{2}(0,\ell)\times\mathbb{R}^{2}$ is given by

$(\Psi f)(x):=Y(T,x,f)=\left(y(T,x),y(T,0),y(T,\ell)\right),\qquad
x\in(0,\ell),$

and the Tikhonov functional by

$\displaystyle J_{\varepsilon}(f)$
$\displaystyle=\frac{1}{2}\left\|Y(T,\cdot,f)-Y_{T}^{\delta}\right\|_{L^{2}(0,\ell)\times\mathbb{R}^{2}}^{2}+\frac{\varepsilon}{2}\|f\|_{L^{2}(0,\ell)}^{2},\qquad
f\in L^{2}(0,\ell),$
$\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\|f\|_{L^{2}(0,\ell)}^{2}\right),$

where
$Y_{T}^{\delta}:=\left(y_{T}^{\delta},y_{T}^{0,\delta},y_{T}^{\ell,\delta}\right)\in
L^{2}(0,\ell)\times\mathbb{R}^{2}$. The adjoint system corresponding to (3.8)
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.2)

Similarly to calculations in Section 3, the gradient of $J_{\varepsilon}$ is
given by

$J_{\varepsilon}^{\prime}(f)(x)=\int_{0}^{T}\varphi(t,x,f)r(t,x)\,\mathrm{d}t+\varepsilon
f(x),\quad f\in L^{2}(0,\ell),\;x\in(0,\ell).$ (5.3)

The gradient formula (5.3) for the Tikhonov functional $J_{\varepsilon}$
allows us to implement classical versions of the Conjugate Gradient Algorithm
with different conjugation coefficients. Here, we shall restrict ourselves to
the following basic algorithm.

1 Set $k=0$ and choose an initial source $f_{0}$;

2 Solve the direct problem (5.1) to obtain $Y(t,x,f_{k})$;

3 Knowing the computed $Y(T,x,f_{k})$ and the measured $Y_{T}^{\delta}$, solve
the adjoint problem (5.2) to obtain $\varphi(t,x,f_{k})$;

4 Compute the descent direction $p_{k}=J_{\varepsilon}^{\prime}(f_{k})$ using
(5.3);

5 Solve the direct problem (5.1) with source $p_{k}$ to get the solution $\Psi
p_{k}$;

6 Compute the relaxation parameter
$\displaystyle\alpha_{k}=\frac{\|p_{k}\|_{L^{2}(0,\ell)}^{2}}{\|\Psi
p_{k}\|_{L^{2}(0,\ell)\times\mathbb{R}^{2}}^{2}}$ (see [17, Lemma 3.4.1]);

7 Find the next iteration $f_{k+1}=f_{k}-\alpha_{k}p_{k}$;

8 Stop the iteration process if the stopping criterion
$J_{\varepsilon}(f_{k+1})<e_{J}$ holds. Otherwise, set $k:=k+1$ and go to Step
2;

Algorithm 1 Landweber iteration scheme

Since we deal with dynamic boundary conditions that contain time derivative,
the numerical solutions for the direct problem (5.1) and the adjoint problem
(5.2) will be obtained by using the method of lines. This method is
implemented in our case with help of the Mathematica system [31] via the
function NDSolve for solving ordinary differential systems. The noisy terminal
data will be generated from the exact output data as follows

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

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

In the subsequent numerical tests, we will choose, for simplicity, the
following values for the known parameters

$T=1,\quad\ell=1,\quad r=1,\quad y_{0}=0,\quad a=b=0.$

Example 1
We consider the following basic source term $f(x)=x(1-x),\;x\in(0,1)$. First,
we compute the numerical solution $Y(t,\cdot,f)$ to generate the data
$Y(1,\cdot,f)$. Next, we plot the corresponding solution.

Figure 1. Solution $Y(t,x,f)$ of (5.1). Figure 2. Exact and recovered $f(x)$
by using Algorithm 1, for $p\in\\{1\%,3\%,5\%\\}$, respectively.

The initial iteration is chosen as $f_{0}=0$ with regularization parameter and
stopping parameter $\varepsilon=e_{J}=10^{-6}$. The algorithm stops at
iterations $k\in\\{2,2,2\\}$, for $p\in\\{1\%,3\%,5\%\\}$, respectively.

######  Remark 5.1.

The optimal regularization parameter $\varepsilon^{\mathrm{opt}}$ can be
defined from the conditions $\varepsilon<1$ and
$\frac{\delta^{2}}{\varepsilon}<1$ depending on the noise level $\delta$. The
optimal stopping parameter $e_{J}^{\mathrm{opt}}$ can be defined by analyzing
the behavior of the convergence error and the accuracy error defined
respectively by

$\displaystyle e(k,f_{k})$ $\displaystyle:=\|\Psi
f_{k}-Y_{T}\|_{L^{2}(0,1)\times\mathbb{R}^{2}}^{2}$ (5.4) $\displaystyle
E(k,f_{k})$ $\displaystyle:=\|f-f_{k}\|_{L^{2}(0,1)},$

in terms of the iteration number $k$. We refer to [17] for more details.

Example 2
We take the exact source term as $f(x)=\sin(\pi x),\;x\in(0,1)$.

Figure 3. Exact and recovered $f(x)$ by using Algorithm 1, for
$p\in\\{1\%,3\%,5\%\\}$, respectively.

The initial iteration is chosen as $f_{0}=0$ with regularization parameter and
stopping parameter $\varepsilon=e_{J}=10^{-8}$. The algorithm stops at
iterations $k\in\\{3,4,4\\}$, for $p\in\\{1\%,3\%,5\%\\}$, respectively.

Example 3
We take the exact source term as
$f(x)=\exp\left(-8\left(x-\frac{1}{2}\right)^{2}\right),\;x\in(0,1)$.

Figure 4. Exact and recovered $f(x)$ by using Algorithm 1 for
$p\in\\{1\%,3\%,5\%\\}$, respectively.

The initial iteration is chosen as $f_{0}=0$ with regularization parameter and
stopping parameter $\varepsilon=e_{J}=10^{-8}$. The algorithm stops at
iterations $k\in\\{5,5,35\\}$ for $p\in\\{1\%,3\%,5\%\\}$, respectively.

Table 1. Errors depending on the iteration number $k$ for noise free data
($p=0\%$).

$\begin{array}[]{cccccc}\hline\cr k&1&2&3&4&5\\\ \hline\cr\\\\[-8.53581pt]
e\left(k,f_{k}\right)&3.421\times 10^{-1}&1.101\times 10^{-2}&4.517\times
10^{-4}&2.767\times 10^{-4}&1.413\times 10^{-4}\\\ \hline\cr\\\\[-8.53581pt]
E\left(k,f_{k}\right)&2.839\times 10^{-1}&3.471\times 10^{-2}&3.469\times
10^{-2}&9.412\times 10^{-3}&9.409\times 10^{-3}\\\ \hline\cr\end{array}$

######  Remark 5.2.

The above numerical experiments show that the Landweber scheme yields stable,
accurate and fast results for the reconstruction of unknown source terms in
heat equation with dynamic boundary conditions. We clearly see that the
recovery of the source $f(x)$ becomes more accurate as the noise level $p$
decreases. This approach can be adapted for the simultaneous recovery of
internal and boundary source terms that depend on both time and space, but it
might require a large number of iterations. A performance analysis for this
case will be treated in a forthcoming paper.

##  6\. Conclusions

In this paper, we have considered an inverse problem for determining internal
and boundary source terms from final time data in heat equation with dynamic
boundary conditions. Adapting the weak solution approach, a minimization
problem for the Tikhonov functional is analyzed, and a gradient formula of the
functional is established via the solution of an appropriate adjoint system.
Then the Lipschitz continuity of the Fréchet gradient is proved. Using
calculus of variations techniques, the existence and the uniqueness of a
quasi-solution are proved. In particular, a sufficient condition for the
uniqueness is presented. Finally, some numerical tests are provided for
recovering an internal heat source in the one-dimensional case.

## 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._ , (2021), doi: 10.1515/jiip-2020-0067.
  * [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. Equat. and Cont. Theo._ , (2020), doi: 10.3934/eect.2020094.
  * [3] I. Boutaayamou, S. E. Chorfi, L. Maniar and O. Oukdach,  The cost of approximate controllability of heat equation with general dynamical boundary conditions,  _Portugal Math._ , 78 (2021), 65–99.
  * [4] A. L. Bukhgeim and M. V. Klibanov,  Global Uniqueness of a class of multidimensional inverse problems,  _Sov. Math. Dokl._ , 24 (1981), 244–247.
  * [5] P. DuChateau,  An introduction to inverse problems in partial differential equations for for physicists, scientists and engineers,  in _Parameter Identification and Inverse Problems in Hydrology, Geology and Ecology_ (Karlsruhe, 1995), (1996), 3–50.
  * [6] P. DuChateau, R. Thelwell and G. Butters,  Analysis of an adjoint problem approach to the identification of an unknown diffusion coefficient,  _Inverse Problems_ , 20 (2004), 601–625.
  * [7] H. W. Engl, M. Hanke and A. Neubauer,  _Regularization of Inverse Problems_ ,  Kluwer Academic Publishers, 2000.
  * [8] A. Erdem,  A simultaneous approach to inverse source problem by Green’s function,  _Math. Meth. Appl. Sci._ , 38 (2015), 1393–1404.
  * [9] 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.
  * [10] S. Gnanavel, N. Barani Balan and K. Balachandran,  Identification of source terms in the Lotka–Volterra system,  _J. Inverse Ill-Posed Probl._ , 20 (2012), 287–312.
  * [11] G. R. Goldstein,  Derivation and physical interpretation of general boundary conditions,  _Adv. Diff. Equ._ , 11 (2006), 457–480.
  * [12] 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.
  * [13] A. Hasanov,  Identification of an unknown source term in a vibrating cantilevered beam from final overdetermination,  _Inverse Problems_ , 25 (2009), 115015 (19pp).
  * [14] 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.
  * [15] A. Hasanov, M. Otelbaev and B. Akpayev,  Inverse heat conduction problems with boundary and final time measured output data,  _Inverse Probl. Sci. Eng._ , 19 (2011), 985–1006.
  * [16] 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.
  * [17] A. Hasanov and V. G. Romanov,  _Introduction to Inverse Problems for Differential Equations_ ,  Springer, New York, 2017.
  * [18] V. Isakov,  _Inverse Source Problems_ ,  Mathematical Surveys and Monographs, vol. 34, Amer. Math. Soc., Providence, RI, 1990.
  * [19] M. I. Ismailov,  Inverse source problem for heat equation with nonlocal Wentzell boundary condition,  _Results Math_ , 73 (2018), 68.
  * [20] M. Kaya, A. Erdem,  Simultaneous reconstruction of the source term and the surface heat transfer coefficient,  _Math. Meth. Appl. Sci._ , 38 (2015], 517–526.
  * [21] 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.
  * [22] A. Khoutaibi, L. Maniar, D. Mugnolo and A. Rhandi,  Parabolic equations with dynamic boundary conditions and drift terms,  To appear in _Mathematische Nachrichten_.
  * [23] M. V. Klibanov,  Carleman estimates for global uniqueness, stability and numerical methods for coefficient inverse problems,  _J. Inverse Ill-Posed Probl._ , 21 (2013), 477–560.
  * [24] M. Kumpf and G. Nickel,  Dynamic boundary conditions and boundary control for the one-dimensional heat equation,  _J. Dynam. Control Systems_ , 10 (2004), 213–225.
  * [25] R. E. Langer,  A problem in diffusion or in the flow of heat for a solid in contact with a fluid,  _Tohoku Math. J._ , 35 (1932), 260–275.
  * [26] L. Maniar, M. Meyries and R. Schnaubelt,  Null controllability for parabolic equations with dynamic boundary conditions of reactive-diffusive type,  _Evol. Equat. and Cont. Theo._ , 6 (2017), 381–407.
  * [27] A. Miranville and S. Zelik,  Exponential attractors for the Cahn-Hilliard equation with dynamic boundary conditions,  _Math. Meth. Appl. Sci._ , 28 (2005), 709–735.
  * [28] M. Slodička,  A parabolic inverse source problem with a dynamical boundary condition,  _Appl. Math. Comput._ , 256 (2015), 529–539.
  * [29] F. P. Vasil’ev,  _Methods for Solving Extremal Problems_ ,  Nauka, Moscow, 1981.
  * [30] R. D. Vold and M. J. Vold,  _Colloid and Interface Chemistry_ ,  Addison-Wesley, Reading, Mass, 1983.
  * [31] S. Wolfram,  _The Mathematica Book_ ,  Wolfram Media, Cambridge, UK, 2005.
  * [32] E. Zeidler,  _Nonlinear Functional Analysis and Its Applications, II/B Nonlinear Monotone Operators_ ,  Springer, New York, 1990.