# Identification of heat conduction coefficient: application to nondestructive testing

Marc Bonnet, Huy Duong Bui, Hubert Maigre, Jacques Planchard

► To cite this version:

Marc Bonnet, Huy Duong Bui, Hubert Maigre, Jacques Planchard. Identification of heat conduction coefficient: application to nondestructive testing. IUTAM Symposium, 1992, Tokyo, Japan. pp.475-488, 10.1007/978-3-642-52439-4_47. hal-00122007

HAL Id: hal-00122007

https://hal.science/hal-00122007v1

Submitted on 17 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Identification of Heat Conduction Coefficient: Application to Nondestructive Testing.

Marc BONNET, Huy Duong BUI, Hubert MAIGRE

Laboratoire de Mécanique des Solides

(centre commun CNRS-Polytechnique-Ecole des Mines-Ecole des Ponts et Chaussées),
Ecole Polytechnique, 91128 Palaiseau Cedex, FRANCE.

Jacques PLANCHARD

Electricité de France, Direction des Etudes et Recherches,
1 av. du Général de Gaulle, 92141 Clamart cedex, FRANCE.

## 1 Introduction

Thermographic Non-Destructive Testing methods for the evaluation of materials and structures are of great interest. They are attractive because of rapid-scanning capabilities available today, Balageas et al [1]. The principle of photothermal methods are well known. The heat flux is applied on the sample by a laser pulse. Then, by means of a Infra-Red camera, one makes a full time record of the whole surface temperature field. There are very rich informations saved in the surface time history $\theta(x, t)$, which can be used for the reconstruction of unknown thermal conduction coefficient $\eta(x)$, defects, cracks etc... For example, the temperature response of an homogeneous body subjected to a pulse heating $Q$ is $\theta(t) = Q/e\sqrt{t}$, where $e = (\eta\rho c)^{1/2}$ is the effusivity. The slope $-1/2$ of the response in the log-log scale means the absence of defects. All perturbations of the straight line are indications of non-homogeneous conduction of heat. Qualitative and quantitative evaluations of defects or inhomogeneities from the knowledge of the surface temperature field belongs to the class of mathematically ill-posed problems.

The numerical solution of thermal inverse problems gives rise to a growing interest in the recent literature, see e.g. [2], [3], [7], [15]. Other different physical contexts share the same mathematical basic structure, notably electric conductivity inverse problems, [6], [8], [9], [10].

Generally, solutions to inverse problems are based on two steps. Firstly, objective functions may be derived for the optimized fitting between experimental and theoretical data. Secondly, one makes use of regularization methods for solving the ill-posed problems.

The paper by H.G. Natke [18] is devoted specifically to the second step, which is also dealt with e.g. in [13], [16], [17], [21], [22]. [23]. The present paper will concentrate on the first problem, which consists of deriving the so-called observation equation. We shall make a review of some mathematical methods for establishing the non-linear relation, which links the thermal coefficient $\eta$ to the surface temperature field, in transient as well as in stationary problems. In particular, we shall discuss problem of determining the optimal correction $\delta\eta$ of the model, in a sensitivity approach of the inverse problem.
## 2 The inverse heat conduction problem

Consider the heat conduction equations with inhomogeneous thermal conduction coefficient $k(x)$ in a three-dimensional solid $\Omega$ of external boundary $S$ ($S \neq \partial\Omega$ in case of a void defect).

$$ \frac{\partial \theta}{\partial t} - \operatorname{div}(k(\mathbf{x}) \nabla \theta) = 0 \quad \mathbf{x} \in \Omega \qquad (1) $$

$$ \theta(\mathbf{x}, t \le 0) = 0 \quad (\text{initial conditions}) \qquad (2) $$

with boundary conditions to be specified.

The direct problem is stated as follows: given the thermal coefficient $k(x)$ such that $k > 0$ and the boundary condition

$$ k(\mathbf{x}^0) \frac{\partial \theta}{\partial n} = f(\mathbf{x}^0, t) \quad (\mathbf{x}^0 \in \partial\Omega, t \ge 0) \qquad (3) $$

find the temperature field $\theta(x, t)$ inside the solid. We shall use the following notations:

* $\theta(x, t)$: the temperature field satisfying the heat equation (1) with unknown coefficient $k(x)$

* $u^n(x, t)$: the temperature field satisfying the heat equation (1) with a given coefficient $\eta(x)$.

We consider the inverse problem of determining the unknown heat conduction coefficient $k(x)$, from boundary data which consist of the heat flux (3) and the surface temperature $\theta(x, t)$. We suppose also that the heat conduction coefficient $k(x) > 0$ is known on the surface $\partial\Omega$ of the solid. Then the simplest formulation of the inverse problem is the minimization of a least-squares best-fit objective function $J(\eta) = J(u^n, \eta)$, where:

$$ J(u, \eta) = \frac{1}{2} \int_0^T \int_S \{u(\mathbf{x}, t) - \theta(\mathbf{x}, t)\}^2 dS_{\mathbf{x}} dt \qquad (4) $$

which is a functional of the difference between the surface temperature of the model $u^n$ and the experimental data $\theta|_S$

Another related inverse problem is the identification of a defect of known thermal characteristics but unknown location and shape.

First, two methods leading to integral equations on the coefficient perturbation $\delta\eta = k(x) - \eta(x)$ are discussed (sections 3 and 4). They allow discussion of uniqueness and stability with respect to the temperature data, as well as linearization methods in the case of small perturbation $\delta\eta(x)$.

Then we consider (section 5) the derivation of the gradient of data-fit functionals using the calculus of variations. This is very useful for many practical applications, since the actual numerical solution of the inverse problem involves a nonlinear optimization problem, which, after discretization, is best solved using conjugate gradient, BFGS or similar algorithms. All these algorithms make use of gradient information, and its computation using finite difference approaches is less efficient and less accurate than using an analytical derivation. Application to the shape identification problem is also discussed in section 6.
### 3 The Fréchet derivative method

Without loss of generality, and for simplicity reason, we consider in this section the particular geometry of a half space which is representative of most experiments in thermal inspections of materials. We shall then use the following additional notations:

$$ \bullet S : \mathbf{x}^0 | \mathbf{x}_3^0 = 0 $$

$$ \bullet \Omega : \mathbf{x} | \mathbf{x}_3 < 0 $$

* $V^η$: the fundamental solution of the heat equation with coefficient $η(\mathbf{x})$.

In order to apply the sensitivity method, one needs the calculation of the functional derivative $D_η J$ of the objective function, or the functional derivative $D_η u$ of the surface value $u^η(\mathbf{x}^0, t)$ with respect to $η$. The last derivative may be interpreted as the linear map $D_η u: \delta\eta \to \delta u$, or the Fréchet derivative defined by

$$ u^{\eta+\delta\eta}(\mathbf{x}^0, t) - u^{\eta}(\mathbf{x}^0, t) = D_{\eta}u^{\eta}(\mathbf{x}^0, t) \delta\eta + O(\delta\eta^2). \quad (5) $$

To determine explicitly the Fréchet derivative $D_η u$, let us consider the fundamental impulse heat flux solution $V^η$ with coefficient $η(\mathbf{x})$, associated to a given point $\mathbf{x}^0 \in S$, which solves:

$$ \frac{\partial V}{\partial t} - \operatorname{div}(\eta(\mathbf{x})\nabla V) = 0 \quad \mathbf{x} \in \Omega \quad (6) $$

$$ V(\mathbf{x}, t \le 0) = 0 \quad \text{(initial conditions)} \quad (7) $$

$$ k(\mathbf{x}^0) \frac{\partial V}{\partial n}(\mathbf{x}^0, t) = \delta(y^0 - x^0)\delta(t) \quad \mathbf{x}^0 \in S, t \ge 0 \quad (8) $$

$$ V \to 0 \text{ as } \|\mathbf{x}\| \to +\infty \quad (9) $$

This solution is identical to the solution of the instantaneous point heat source located at $\mathbf{x}^0$ in an infinite 3D solid, which can be found in textbooks (e.g. Carslaw and Jæger [5]).

The model solution $u^η(\mathbf{x}, t)$ satisfies the heat equation with the same boundary condition (3), with $k(\mathbf{x}^0) = \eta(\mathbf{x}^0)$:

$$ \frac{\partial u}{\partial t} - \operatorname{div}(\eta(\mathbf{x})\nabla u) = 0 \quad \mathbf{x} \in \Omega \quad (10) $$

$$ u(\mathbf{x}, t \le 0) = 0 \quad \text{(initial conditions)} \quad (11) $$

$$ k(\mathbf{x}^0) \frac{\partial u}{\partial n}(\mathbf{x}^0, t) = f(\mathbf{x}^0, t) \quad \mathbf{x}^0 \in S, t \ge 0 \quad (12) $$

The field $u^{η+δη}(\mathbf{x}, t)$ satisfies eqns. (10)-(12) with coefficient $η + δη$, and the same boundary conditions (12). Therefore, $R^η(\mathbf{x}, t; δη) = u^{η+δη}(\mathbf{x}, t) - u^η(\mathbf{x}, t)$ satisfies the equations

$$ \frac{\partial R}{\partial t} - \operatorname{div}(\eta(\mathbf{x})\nabla R) - \operatorname{div}(\delta\eta(\mathbf{x})\nabla u) = 0 \quad \mathbf{x} \in \Omega \quad (13) $$

$$ R(\mathbf{x}, t \le 0) = 0 \quad \text{(initial conditions)} \quad (14) $$

$$ k(\mathbf{x}^0) \frac{\partial R}{\partial n}(\mathbf{x}^0, t) = 0 \quad \mathbf{x}^0 \in S, t \ge 0 \quad (15) $$

$$ R \to 0 \text{ as } \|\mathbf{x}\| \to +\infty \quad (16) $$

Taking the time convolution of (13) with the fundamental solution $V^η$, and observing that $δη \equiv 0$ on $S$, we obtain the formula

$$ \delta u(\mathbf{x}^0, t; \eta) \equiv D_{\eta} u \delta \eta = -\frac{1}{k(\mathbf{x}^0)} \int_{\Omega} \delta \eta(\mathbf{x}) \nabla u^{\eta}(\mathbf{x}, t) * \nabla V^{\eta}(\mathbf{x}, t; \mathbf{x}^0) dV_{\mathbf{x}} \quad (17) $$
where (*) denotes the time convolution of scalar product. Eq. (17) gives the modification of the surface temperature at point x⁰, at time t, from perturbation δη(x). It confirms the "golden rule" according to which the determination of unknown δη(x) defined in a 3D subset of Ω requires data δu defined also on a 3D subspace, here S × [0, T].

Eq. (17) provides a linear relation of the form Aδη = δu, which can be used for determining δη which makes the best fit between data δu and the prediction Aδη. However the linear system Aδη - δu = 0 is a Fredholm integral equation of first kind. Hence it belongs to the class of ill-posed problems [22], and cannot be used directly in the inverse problem. Instead of solving exactly this equation, one can determine the steepest direction δη for minimizing the function J. Since

$$ \delta J = - \int_{0}^{T} \int_{S} \int_{\Omega} \frac{1}{k(y)} (u^{\eta}(y, t) - \theta(y, t)) \nabla u^{\eta}(x, t) * \nabla V^{\eta}(x, t; y) \delta \eta(x) dt dS_y dV_x \quad (18) $$

One can see from (18) that δJ = (β, δη) is a linear functional of δη. Therefore, for all perturbation of the same norm, the perturbation δη corresponding to the maximum variation of δJ is proportional to β. Consequently

$$ \delta \eta(\mathbf{x}) = \lambda \int_{0}^{T} \int_{S} \frac{1}{k(y)} (u^{\eta}(y, t) - \theta(y, t)) \nabla u^{\eta}(\mathbf{x}, t) * \nabla V^{\eta}(\mathbf{x}, t; y) \delta \eta(\mathbf{x}) dt dS_y \quad (19) $$

This result allows to calculate explicitly the steepest gradient of the objective function.

## 4 The adjoint fields method

For simplicity reason, we study the inverse problem for identifying the heat coefficient η(x) = 1 + h(x), with h(x⁰) = 0 on the surface S of the solid. The perturbation h(x) represents physically the damaged zone in the solid. Generally, damage by microcraking may change the heat coefficient, but neither the density ρ nor the specific heat c. Normalizing constants and variables, we assume that the unperturbed heat coefficient is equal to unity.

Let us consider a family (λ) of experiments on the actual solid, with the unknown coefficient k(x) = 1 + h(x), and the heat flux f_λ(x⁰, t). The response solution to the boundary condition is denoted by θ_λ(x, t), which satisfies the following equations

$$ \frac{\partial}{\partial t} \theta_{\lambda}(\mathbf{x}, t) - \mathrm{div}(k(\mathbf{x}) \nabla \theta_{\lambda}(\mathbf{x}, t)) = 0 \quad \mathbf{x} \in \Omega \quad (20) $$

$$ \theta_{\lambda}(\mathbf{x}, t)(\mathbf{x}, t \le 0) = 0 \quad (\text{initial condition}) \quad (21) $$

$$ \frac{\partial}{\partial n} \theta_{\lambda}(\mathbf{x}^0, t)(\mathbf{x}, t) = f_{\lambda}(\mathbf{x}^0, t) \quad \mathbf{x}^0 \in S, t \ge 0 \quad (22) $$

where h(x⁰) = 0 (i.e. k(x⁰) = 1) on S has been taken into account. The inverse problem under consideration is to determine h(x) from surface measurements of f_λ(x⁰, t) and θ_λ(x⁰, t), during the time interval [0, T]. Here, we have much more data than necessary, so that the inverse problem is overdetermined.

Let us introduce two auxiliary fields φ_λ and ψ_μ. The first field, for each λ, is solution of the unperturbed problem with the unit heat coefficient and the same surface data (22).

$$ \frac{\partial}{\partial t} \phi_{\lambda}(\mathbf{x}, t) - \mathrm{div}(\nabla \phi_{\lambda}(\mathbf{x}, t)) = 0 \quad \mathbf{x} \in \Omega \quad (23) $$

$$ \phi_{\lambda}(\mathbf{x}, t)(\mathbf{x}, t \le 0) = 0 \quad (\text{initial condition}) \quad (24) $$

$$ \frac{\partial}{\partial n} \phi_{\lambda}(\mathbf{x}^0, t)(\mathbf{x}, t) = f_{\lambda}(\mathbf{x}^0, t) \quad \mathbf{x}^0 \in S, t \ge 0 \quad (25) $$
The second field depending on another family of parameter (μ) is solution of the adjoint problem

$$ -\frac{\partial}{\partial t}\psi_{\mu}(\mathbf{x}, t) - \operatorname{div}(\nabla\psi_{\mu}(\mathbf{x}, t)) = 0 \quad \mathbf{x} \in \Omega \qquad (26) $$

$$ \psi_{\mu}(\mathbf{x}, t \ge T) = 0 \quad (\text{final condition}) \qquad (27) $$

$$ \frac{\partial}{\partial n}\psi_{\mu}(\mathbf{x}^0, t) = g_{\mu}(\mathbf{x}^0, t) \quad \mathbf{x}^0 \in S, t \le T \qquad (28) $$

The unperturbed direct problems (23)-(25) and the unperturbed adjoint problems (26)-(28) are well-posed problems, which can be solved by the finite element method. We assume that the auxiliary fields are known for all $\lambda$ and $\mu$. Multiplying (20) by $\psi_\mu$ and (26) by $\theta_\lambda$, then combining the results, yields

$$ \int_0^T \int_{\Omega} \operatorname{div}\{\theta_{\lambda} \nabla \psi_{\mu} - \psi_{\mu} \nabla \theta_{\lambda}\} dV_{\mathbf{x}} dt - \int_0^T \int_{\Omega} \psi_{\mu} \operatorname{div}\{h \nabla \theta_{\lambda}\} dV_{\mathbf{x}} dt = 0 \qquad (29) $$

Upon suitable application of the divergence formula and taking into account the boundary conditions (22), (28), eq. (29) can be arranged in the following form

$$ \int_{\Omega} h(\mathbf{x}) K^{NL}(\mathbf{x}; \lambda, \mu, T) dV_{\mathbf{x}} = B(\lambda, \mu, T) \qquad (30) $$

where

$$ K^{NL}(\mathbf{x}; \lambda, \mu, T) = \int_{0}^{T} \nabla \psi_{\mu} \cdot \nabla \theta_{\lambda} dt \qquad (31) $$

$$ B(\lambda, \mu, T) = \int_0^T \int_S (\theta_\lambda g_\mu - \psi_\mu f_\lambda) dt dS_x \qquad (32) $$

The integrand in the right-hand side of (32) is a known surface quantity. But a look at (31) shows that the kernel $K^{NL}$ depends on the unknown field $\theta_\lambda$, hence eq. (31) is non-linear with respect to $h(\mathbf{x})$ (the superscript NL means non-linear kernel).

An alternative expression of (32) may be derived as follows:

$$ \begin{aligned} \int_0^T \int_S \psi_\mu f_\lambda dS_\mathbf{x} dt &= \int_0^T \int_S \psi_\mu \frac{\partial}{\partial n} \phi_\lambda dS_\mathbf{x} dt \\ &= \int_0^T \int_\Omega \operatorname{div} \psi_\mu \nabla \phi_\lambda dV_\mathbf{x} dt \\ &= \int_0^T \int_\Omega \left\{ \phi_\lambda \frac{\partial}{\partial t} \psi_\mu + \nabla \phi_\lambda \cdot \nabla \psi_\mu \right\} dV_\mathbf{x} dt \\ &= \int_0^T \int_S \phi_\lambda \frac{\partial}{\partial n} \psi_\mu dS_\mathbf{x} dt + \int_0^T \int_\Omega \frac{\partial}{\partial t} \psi_\mu \phi_\lambda dV_\mathbf{x} dt \\ &= \int_0^T \int_S \phi_\lambda \frac{\partial}{\partial n} \psi_\mu dS_\mathbf{x} dt \end{aligned} \qquad (33) $$

$$ B(\lambda, \mu, T) = \int_{0}^{T} \int_{S} (\theta_{\lambda} - \phi_{\lambda}) g_{\mu} dt dS_{\mathbf{x}} \qquad (34) $$

In particular, application of (34) to the special case $\psi_\mu(\mathbf{x}, t) = V(\mathbf{x}, T - t; x_\mu)$ (with $x_\mu \in S$, $t \in [0, T]$) gives:

$$ B(\lambda, \mu, T) = \theta_{\lambda}(\mathbf{x}, t) - \phi_{\lambda}(\mathbf{x}, t) \qquad (35) $$

Then, (34) and (37) provide a generalization of the Frechet derivative (17) for perturbations $h(\mathbf{x})$ of arbitrary amplitude:

$$ \theta_{\lambda}(\mathbf{x}, t) - \phi_{\lambda}(\mathbf{x}, t) = \int_{\Omega} h(\mathbf{x}) \int_{0}^{T} \nabla \psi_{\mu} \cdot \nabla \theta_{\lambda} dt dV_{\mathbf{x}} \qquad (36) $$
The adjoint method is also interesting when the perturbation $h(x)$ is small enough to allow the linearization process $\theta = \phi + O(h)$, or $\theta \sim \phi$. Substituting $\theta_\lambda$ by $\phi_\lambda$ in the kernel (31) yields the linear equation for $h(x)$:

$$ \int_{\Omega} h(x) K(x; \lambda, \mu, T) dV_x = B(\lambda, \mu, T) \quad (37) $$

where $B$ is still given by (32) and $K$ is a known kernel:

$$ K(x; \lambda, \mu, T) = \int_{0}^{T} \nabla \psi_{\mu} \cdot \nabla \phi_{\lambda} dt \quad (38) $$

Equations (37), (38) is nothing but the generalization to transient heat problem of the method suggested by Calderon (see [6], [10]) for stationary problem. Furthermore, eq. (36) then becomes identical to (17) (apart from notation differences).

# 5 Variational formulation of the inverse conductivity problem

Let us reformulate the model problem (10)-(12) in the weak form. Set the following definitions and notations

$$ \mathcal{V} = \{v \in H^1(\Omega), v(x, 0) = 0 (x \in \Omega), v(x^0, t) = 0 (x^0 \in S, t \ge 0)\} \quad (39) $$

$$ u' = \frac{\partial u}{\partial t} \quad (40) $$

$$ \langle w, v \rangle_{\Omega} = \int_{\Omega} w(x)v(x) dV_x \quad (41) $$

$$ (f, v)_S = \int_S f(x)v(x) dS_x \quad (42) $$

$$ a_{\Omega}^{\eta}(u, v) = \int_{\Omega} \eta(x) \nabla u \cdot \nabla v \, dV_x \quad (43) $$

Then the model problem with coefficient $\eta$ is the solution of the variational problem ('state equation')

$$ \int_{0}^{T} \langle u', v \rangle_{\Omega} dt + \int_{0}^{T} a_{\Omega}^{\eta}(u, v) dt = \int_{0}^{T} \langle f, v \rangle_{S} dt \quad \forall v \in \mathcal{V} \quad (44) $$

The optimization problem (min$_\eta$ $J(\eta)$) with respect to $\eta$ belongs to the optimal control theory for partial derivative equations, cf. Lions [14]. It consists of minimizing $J(u; \eta)$ with the constraints on $u$ which satisfies the variational equation (44). The classical approaches to this minimization problem, with constraints, make use of the Lagrange multiplier $\psi(x, t)$ and the Lagrangian $\mathcal{L}$

$$ \mathcal{L}(u; \eta, \psi) = J(u; \eta) + \int_{0}^{T} \{\langle u', \psi \rangle_{\Omega} + a_{\Omega}^{\eta}(u, v) - \langle f, \psi \rangle_{S}\} dt \quad (45) $$

It is clear that $\mathcal{L} \equiv J$ when (44) is satisfied. The optimal solution $u$ minimizing $J$ with the constraints (44) is the stationary point of $\mathcal{L}$. One observes that the stationarity condition $\delta\mathcal{L} = 0$ under fixed $u$ and $\eta$, and arbitrary $\delta\psi$ yields the variational equation ((44)) with $v \equiv \delta\psi$.

Consider now the variation of $\mathcal{L}$ due to $\delta u$ and $\delta\eta$:

$$ \delta \mathcal{L} = \frac{\partial \mathcal{L}}{\partial u} \delta u + \frac{\partial \mathcal{L}}{\partial \eta} \delta \eta \quad (46) $$

$$ \frac{\partial \mathcal{L}}{\partial u} \delta u = \int_{0}^{T} \{ \langle u - \theta, \delta u \rangle_{S} + \langle \delta u', \psi \rangle_{\Omega} + a_{\Omega}^{\eta}(\delta u, \psi) \} dt \quad (47) $$

$$ \frac{\partial \mathcal{L}}{\partial \eta} \delta \eta = \int_{0}^{T} a_{\Omega}^{\eta}(u, \psi) dt \quad (48) $$
We now restrict the choice of $\psi$ in such a way that $\delta\mathcal{L} = 0$ for $\delta\eta \equiv 0$, that is, we put:

$$ \frac{\partial \mathcal{L}}{\partial u} \delta u = 0 \quad \forall \delta u \in \mathcal{V} \qquad (49) $$

By integrating the time integral by parts in (47), we obtain

$$ \int_0^T \{ (u - \theta, \delta u)_S + a_{\Omega}^{\eta}(\delta u, \psi) - \langle \delta u, \psi' \rangle_{\Omega} \} dt + \langle \delta u, \psi' \rangle_{\Omega} \Big|_0^T = 0 \quad \forall \delta u \in \mathcal{V} \qquad (50) $$

Since $\delta u(x, 0) = 0$, one arrives at the adjoint backward heat equation for the lagrangian multiplier $\psi(x, t)$

$$ -\frac{\partial \psi}{\partial t} - \text{div}(\eta(\mathbf{x})\nabla\psi) = 0 \quad \mathbf{x} \in \Omega \qquad (51) $$

$$ \psi(x, t \ge T) = 0 \quad (\text{final condition}) \qquad (52) $$

$$ k(x^0) \frac{\partial \psi}{\partial n}(x^0, t) = (u^\eta - \theta)(x^0, t) \quad x^0 \in S, t \le T \qquad (53) $$

where $u^\eta$ is the solution of the state equation (44). This system corresponds to the back diffusion equation, with the final condition (52) and boundary condition (53), which is a well-posed problem. Its variational formulation reads:

$$ \int_0^T -\langle \psi', w \rangle_{\Omega} dt + \int_0^T a_{\Omega}^{\eta}(\psi, w) dt = \int_0^T \langle u^{\eta} - \theta, w \rangle_S dt \quad \forall w \in \mathcal{V}' \qquad (54) $$

where $\mathcal{V}'$ is a test function space similar to $\mathcal{V}$ but with a final condition $w(x, T) = 0, \forall w \in \mathcal{V}'$.

Finally, this particular choice $\psi^\eta(x, t)$ for the Lagrange multiplier $\psi(x, t)$ yields the following formula

$$ \delta J(\eta) = \delta \mathcal{L}(u^\eta; \psi^\eta, \eta) = \int_0^T a_{\Omega}^{\delta\eta}(u^\eta, \psi^\eta) dt \qquad (55) $$

Comparison with formula (18) shows that $\psi$ is related to the mismatch $(u^\eta - \theta)$ of surface data by the fundamental solution $\mathcal{V}$. Although formulae (55) and (18) are different in their presentation, they represent essentially the same result.

## Some remarks about the result (55)

Equation (55) expresses analytically the gradient of $J(\eta)$ with respect to $\eta(x)$. It is therefore a very valuable tool for usual nonlinear optimization strategies, since:

1. It is an exact expression (provided the state and adjoint equations are solved exactly).

2. The whole gradient of $J$ is computed using only one adjoint equation per objective function considered, whereas a direct differentiation (either analytical or numerical) would require to set up as many auxiliary boundary/initial value problems as the number of design variables describing the conductivity field $k(x)$.

3. The adjoint variational equation (54) involves the same bilinear forms $\langle \cdot, \cdot \rangle_{\Omega}$ and $a_{\Omega}^{\eta}(\cdot, \cdot)$ than the state equation. In a finite element discretization approach, it means that the corresponding matrices are assembled only once. The numerical setting up of the adjoint equation needs only the building of a new right-hand side.
## 6 Variational formulation of the defect shape inverse problem

Let us consider the case where the physical nature of the defect is known (e.g. cavity, inclusion of a known material...). For simplicity of presentation, the present discussion is restricted to volumic defects, although other types may be considered (e.g. interface defects). Accordingly we put $\Omega = \Omega_1 \cup \Omega_2$, where $\Omega_1$ and $\Omega_2$ (the defect) have known thermal conductivities $k_1(x)$ and $k_2(x)$ respectively. We consider a variant of the variational approach of section 5 above, in which the unknown is the defect boundary $\Gamma = \partial\Omega_2$ instead of $k(x)$. Accordingly, the least-squares functional (4) is considered as a function of $\Gamma: J \equiv J(\Gamma)$.

As in the previous section, the minimization of $J(\Gamma)$ is constrained: $u$ must solve

$$ \frac{\partial u}{\partial t} - \operatorname{div}(\eta_i(x)\nabla u) = 0 \quad x \in \Omega_i, i = 1, 2 \qquad (56) $$

$$ u(x, t \le 0) = 0 \quad (\text{initial condition}) \qquad (57) $$

$$ k(x^0) \frac{\partial u}{\partial n}(x^0, t) = f(x^0, t) \quad x^0 \in S, t \le T \qquad (58) $$

$$ [u](x, t) = [\eta \frac{\partial u}{\partial n}](x, t) = 0 \quad x \in \Gamma \text{ (continuity accross } \Gamma) \qquad (59) $$

where $[\cdot](\mathbf{x}) = (\cdot)_2(\mathbf{x}) - (\cdot)_1(\mathbf{x})$ ($\mathbf{x} \in \Gamma$) denotes the jump accross $\Gamma$. The Lagrangian $\mathcal{L}$ for the problem under consideration accordingly reads:

$$ \mathcal{L}(u; \psi, \Gamma) = J(u; \Gamma) + \int_{0}^{T} \{ (u', \psi)_{\Omega} + a_{\Omega}^{n}(u, v) - (f, \psi)_{S} \} dt \quad \psi \in \mathcal{V} \qquad (60) $$

where $\mathcal{V}$ is defined by (40) (in particular, each $\psi \in \mathcal{V}$ is continuous accross $\Gamma$).

The derivative of $\mathcal{L}(u; \psi, \Gamma)$ with respect to $\Gamma$, necessary for the minimization of $J(\Gamma)$ using standard algorithms, is provided by the shape differentiation approach [19]. Let $\Gamma$ denote the current location of the unknown boundary during the minimization process, and assume a further evolution of the surface described by means of a time-like parameter $\tau$ and a normal 'velocity' field $v_n$:

$$ \Gamma(\tau) = \Gamma + v_n n \tau \qquad (61) $$

while the external boundary $S$ remains fixed ($v_n(x) = 0 \quad \forall x \in S$). Then the derivative $\frac{d\mathcal{L}}{d\tau}$ of a functional $\mathcal{L}$ is a linear form of the field $v_n$. In what follows, all derivatives with respect to $\tau$ are implicitly taken for $\tau = 0$.

Various formulas are given in the literature (see e.g. [19]) for the derivative of integrals with respect to variable volumes $\Omega$ or surfaces $\Gamma$, among which:

$$ \frac{d}{d\tau} \int_{\Omega} a(\mathbf{x}, \tau) dV_{\mathbf{x}} = \int_{\Omega} \frac{\partial}{\partial \tau} a(\mathbf{x}, \tau) dV_{\mathbf{x}} + \int_{\partial \Omega} a(\mathbf{x}, \tau) v_n(\mathbf{x}) dS_{\mathbf{x}} \qquad (62) $$

$$ \frac{d}{d\tau} \int_{\Gamma} a(\mathbf{x}, \tau) dS_{\mathbf{x}} = \int_{\Gamma} \left\{ \frac{\partial}{\partial \tau} a(\mathbf{x}, \tau) + \left( \frac{\partial}{\partial n} a(\mathbf{x}, \tau) - 2K(\mathbf{x}) a(\mathbf{x}, \tau) \right) v_n(\mathbf{x}) \right\} dS_{\mathbf{x}} \qquad (63) $$

where $K(\mathbf{x})$ denotes the mean curvature at $\mathbf{x} \in \Gamma$. Equation (63) holds only for a closed smooth surface, while in equation (62) $v_n$ refers to the unit normal $n$ directed towards the exterior of $\Omega$. Application of formulas (62)-(63) above gives:
$$ \frac{d}{d\tau} \mathcal{L} = \frac{\partial \mathcal{L}}{\partial u} \frac{\partial u}{\partial \tau} + \frac{\partial \mathcal{L}}{\partial \Gamma} v_n \quad (64) $$

$$ \frac{\partial \mathcal{L}}{\partial u} \frac{\partial u}{\partial \tau} = \int_{0}^{T} \left\{ \langle u - \theta, \frac{\partial u}{\partial \tau} \rangle s + \left\langle \frac{\partial u'}{\partial \tau}, \psi \right\rangle_{\Omega} + a_{\Omega_1}^{\eta_1} \left( \frac{\partial u}{\partial \tau}, \psi \right) + a_{\Omega_2}^{\eta_2} \left( \frac{\partial u}{\partial \tau}, \psi \right) \right\} dt \quad (65) $$

$$ \frac{\partial \mathcal{L}}{\partial \Gamma} v_n = \int_0^T \langle [u'\psi + \eta \nabla u \cdot \nabla \psi], v_n \rangle_{\Gamma} dt \quad (66) $$

where the condition $v_n \equiv 0$ on $S$ has been taken into account.

Now the choice of $\psi$ is restricted in such a way that $\frac{d}{d\tau} \mathcal{L} = 0$ for $v_n \equiv 0$, that is, we put:

$$ \frac{\partial \mathcal{L}}{\partial u} \frac{\partial u}{\partial \tau} = 0 \quad \forall \frac{\partial u}{\partial \tau} \in \mathcal{V} \quad (67) $$

By integrating the time integral by parts in (65), (67) and (65) give:

$$ \left( \forall \frac{\partial u}{\partial \tau} \in \mathcal{V} \right) \int_0^T \left\{ \langle u - \theta, \frac{\partial u}{\partial \tau} \rangle s + a_{\Omega_1}^{\eta_1} \left( \frac{\partial u}{\partial \tau}, \psi \right) + a_{\Omega_2}^{\eta_2} \left( \frac{\partial u}{\partial \tau}, \psi \right) - \left\langle \frac{\partial u}{\partial \tau}, \psi' \right\rangle_{\Omega} \right\} dt + \left\langle \frac{\partial u}{\partial \tau}, \psi' \right\rangle_{\Omega} \Big|_0^T = 0 \quad (68) $$

The statement of the (backward) adjoint problem readily follows:

$$ -\frac{\partial \psi}{\partial t} - \operatorname{div}(\eta_i(\mathbf{x})\nabla \psi) = 0 \quad \mathbf{x} \in \Omega_i, i = 1, 2 \quad (69) $$

$$ \psi(\mathbf{x}, t \ge T) = 0 \quad \text{(final condition)} \quad (70) $$

$$ k(\mathbf{x}^0) \frac{\partial \psi}{\partial n}(\mathbf{x}^0, t) = (u_{\Gamma} - \theta)(\mathbf{x}^0, t) \quad \mathbf{x}^0 \in S, t \le T \quad (71) $$

$$ \langle [\psi](\mathbf{x}, t) = [\eta \frac{\partial \psi}{\partial n}] (\mathbf{x}, t) = 0 \quad \mathbf{x} \in \Gamma \text{ (continuity accross } \Gamma \quad (72) $$

Finally, the jump conditions (59) on $\Gamma$ imply that

$$ [u'_{\Gamma} \psi + \eta \nabla u_{\Gamma} \cdot \nabla \psi] = \eta_1 \frac{\partial u_{\Gamma}}{\partial n} \left[ \frac{\partial \psi}{\partial n} \right] \quad (73) $$

and, as a result, the derivative $\frac{d}{d\tau} J(\Gamma)$ is given by:

$$ \frac{d}{d\tau} J(\Gamma) = \frac{d}{d\tau} \mathcal{L}(u_{\Gamma}; \Gamma, \psi_{\Gamma}) = \int_{0}^{T} \langle \eta_1 \frac{\partial u_{\Gamma}}{\partial n} \left[ \frac{\partial \psi}{\partial n} \right], v_n \rangle_{\Gamma} dt \quad (74) $$

## Some remarks about the result (74)

1. Expression (74) gives at once the whole gradient of $J(\Gamma)$, and (linearly) depends upon the design variables which describe the current surface $\Gamma$ through the normal 'velocity' $v_n$.

2. Contrarily to the result (55) of the previous section, the derivative of $J(\Gamma)$ is expressed using only boundary integrals, as is always the case in shape differentiation approach [19].

3. It is worth noticing that, in the case of piecewise constant material properties (i.e. $\eta_i(\mathbf{x}) \equiv \eta_i$, $i = 1, 2$), the temperature field $u^{\eta}(\mathbf{x}, t)$ and the adjoint field $\psi^{\eta}(\mathbf{x}, t)$ may be conveniently solved using boundary integral equations (BIE). This, combined to the 'boundary only' character of (74), allows a 'boundary only' treatment of the shape identification problem. The BIE formulation is well-known (see [4] among many references) and will not be repeated here.
7 Use of the constitutive equation error

The so-called ‘constitutive equation error’ *E* is a special type of objective function which has been
considered e.g. in [11] for electric conductivity inverse problems or in [20] for elastic FEM model
updating. For thermal problems, let us assume that the boundary condition (3) holds and that
the temperature field of the solid with defect δη is measured on the entire domain Ω: uη+δη = U
on Ω. The following functional J(u, q; η) is introduced:

$$
\begin{align}
J(u, q; \eta) &= \frac{1}{2} \int_0^T \left\{ \int_{\Omega} \frac{1}{\eta} (q - \eta \nabla u) \cdot (q - \eta \nabla u) \, dV_x + \frac{\gamma}{2} \int_{\Omega} \eta \nabla (u - U) \cdot \nabla (u - U) \, dV_x \right\} dt \nonumber \\
&= \frac{1}{2} \left( \frac{1}{\eta} (q - \eta \nabla u), q - \eta \nabla u \right)_{\Omega} + \frac{\gamma}{2} a_{\Omega}^{\eta}(u - U, u - U) \tag{75}
\end{align}
$$

where *u*, *q* denote 'admissible' temperature and heat flux fields respectively: *u* ∈ V and *q* ∈ B,
with

$$
B = \{ q \mid \operatorname{div}(\eta q) - u' = 0, q \cdot n(\mathbf{x}) = f(\mathbf{x}), \mathbf{x} \in S \} \quad (76)
$$

and $\gamma$ is an adjustable weighing constant, which expresses the expected degree of accuracy of the
measured field $U$. Let $u_\Omega$, $q_\Omega$ and $J(\eta)$ be defined as:

$$
\begin{align}
J(\eta) &= \min_{u \in V, q \in B} J(u, q; \eta) \notag \\
(u_{\Omega}, q_{\Omega}) &= \operatorname{Arg min} J(u, q; \eta) \tag{77}
\end{align}
$$

Then the constitutive equation error *E*(*η*) is defined as:

$$
E(\eta) = \frac{1}{2} \int_{0}^{T} \left\{ \left\langle \frac{1}{\eta}(q_{\Omega} - \eta \nabla u_{\Omega}), q_{\Omega} - \eta \nabla u_{\Omega} \right\rangle_{\Omega} \right\} dt \quad (78)
$$

Let us examine how *E*(*η*), together with the gradient of *J* with respect to *η* are computed in
practice. Due to the constraint *q* ∈ *B*, the following lagrangian *L*, with multiplier field *w*(*x*, *t*), is
introduced:

$$
\mathcal{L}(u, q, w; \eta) = J(u, q; \eta) + \int_{0}^{T} \{ \langle u', w \rangle_{\Omega} + b_{\Omega}(q, w) - \langle f, w \rangle_{S} \} \, dt \quad (79)
$$

$$
b_{\Omega}(q, w) = \int_{\Omega} q \cdot \nabla w \, dV_{x} \qquad (80)
$$

Its variation is given by:

$$
\delta \mathcal{L} = \frac{\partial \mathcal{L}}{\partial q} \delta q + \frac{\partial \mathcal{L}}{\partial u} \delta u + \frac{\partial \mathcal{L}}{\partial \eta} \delta \eta \quad (81)
$$

$$
\frac{\partial \mathcal{L}}{\partial q} \delta q = \int_{0}^{T} \left\{ b_{\Omega}(\delta q, w) + \left\langle \frac{1}{\eta}(\delta q, q - \eta \nabla u) \right\rangle_{\Omega} \right\} dt \quad (82)
$$

$$
\frac{\partial \mathcal{L}}{\partial u} \delta u = \int_{0}^{T} \left\{ \gamma a_{\Omega}^{\eta} (u - U, u - U) + \langle \delta u', w \rangle_{\Omega} - b_{\Omega} (q - \eta \nabla u, \delta u) \right\} dt \quad (83)
$$

$$
\frac{\partial \mathcal{L}}{\partial \eta} \delta \eta = \int_{0}^{T} \left\langle \delta \eta (\nabla u - \frac{q}{\eta}), \nabla u - \frac{q}{\eta} \right\rangle_{\Omega} a_{\Omega}^{\delta \eta} \left( \nabla u - \frac{q}{\eta}, \nabla u + \frac{q}{\eta} \right) dt \quad (84)
$$

First, $J(\eta) = L(u_\Omega, q_\Omega, w_\Omega; \eta)$, where $(u_\Omega, q_\Omega, w_\Omega)$ solve

$$
\forall \delta q \in B, \forall T > 0 \quad \frac{\partial \mathcal{L}}{\partial q} \delta q = 0 \qquad (85)
$$

$$
\forall \delta u \in V, \forall T > 0 \quad \frac{\partial \mathcal{L}}{\partial u} \delta u = 0 \qquad (86)
$$
and the constraint (76). This leads to the following equations:

$$q = \eta \nabla(u - w) \qquad (87)$$

and

$$\forall \delta u \in \mathcal{V} \quad \gamma a_{\Omega}^{\eta}(u - U, \delta u) - \langle \delta u', w \rangle_{\Omega} + b_{\Omega}(\eta \nabla w, \delta u)$$

$$\forall v \in \mathcal{V} \quad \langle u', v \rangle_{\Omega} + a_{\Omega}^{\eta}(u, v) - \langle f, w \rangle_{S} = a_{\Omega}^{\eta}(w, v)$$

The last two equations above can be rewritten as follows, putting $u = u^{\eta} + \Delta u$, $U = u^{\eta} + \Delta U$ with $u^{\eta}$ the solution of the unperturbed direct problem:

$$\forall \delta u \in \mathcal{V} \quad \gamma a_{\Omega}^{\eta}(\Delta u - \Delta U, \delta u) + \langle \delta u, w' \rangle_{\Omega} + b_{\Omega}(\eta \nabla w, \delta u)$$

$$\forall v \in \mathcal{V} \quad \langle \Delta u', v \rangle_{\Omega} + a_{\Omega}^{\eta}(\Delta u, v) = a_{\Omega}^{\eta}(w, v)$$

or, using 'stiffness' and 'mass' operators $K, M$ (e.g. within the finite element method framework):

$$Kw - Mw' = -\gamma K(\Delta u - \Delta U) \qquad (88)$$

$$Kw = K\Delta u + M\Delta u' \qquad (89)$$

Then $\Delta u$ is expressed in terms of $w$ in terms of $w$ using (89):

$$\Delta u = \Delta U - \frac{1}{\gamma} (w - K^{-1}Mw') \qquad (90)$$

and the result is inserted in (88), giving:

$$Kw + \frac{1}{\gamma} (Kw - MK^{-1}Mw'') = K\Delta U + M\Delta U' \qquad (91)$$

Now let us recall that $U = u^{\eta} + \Delta U = u^{\eta+\delta\eta}$ ($u^{\eta}$ is the temperature field solution of the perturbed direct problem), so that the right-hand side of (91) becomes:

$$K\Delta U + M\Delta U' = -\Delta KU \qquad (92)$$

where $\Delta K$ is the perturbation of $K$ induced by the conductivity perturbation $\delta\eta$.

Let us now consider the case $\gamma \gg 1$ in eqn. (91): the measured field $U$ is considered very accurate, and accordingly given a large weight in the functional (75). Then, from eqns. (91) and (92), one has for the solution $w_{\Omega}$ of the coupled equations (88), (89):

$$K w_{\Omega} = -\Delta K U + O\left(\frac{1}{\gamma}\right) \qquad (93)$$

which means that $K w_{\Omega}$, which is computed without actual knowledge of $\Delta K$, takes nonzero values only at points (or on elements, in a FEM approach) where $\delta\eta$ is nonzero. This allows, at least in the rather idealized situation considered here where $U$ is known with great confidence over the entire $\Omega$, the geometrical localization of the defect. Similarly, one can see from (78), (87) and (93) that the value of the constitutive equation error is given by:

$$\begin{aligned} E(\eta) &= \frac{1}{2} a_{\Omega}^{\eta}(w, w) \\ &= -\frac{1}{2} a_{\Omega}^{\delta\eta}(w, U) + O\left(\frac{1}{\gamma}\right) \end{aligned} \qquad (94)$$
Expression (94) is additive with respect to $\Omega$ and can then be split into a sum over a partition of $\Omega$ in subdomains (e.g. finite elements), which can be used to indicate which are the subdomains with nonzero $\delta\eta$.

Then, the subsequent minimization of $J(\eta)$ can be given a reduced size by inspecting the distribution of the density of $E(\eta)$ over $\Omega$ and deciding in advance where $\eta$ has to be corrected. As usual, it can be desirable to use an exact expression of the gradient of $J(\eta)$ with respect to $\eta$:

$$ 
\begin{aligned}
\delta J(\eta) &= \frac{\partial J}{\partial \eta} \delta\eta \\
&= \delta L(u_{\Omega}, q_{\Omega}, w_{\Omega}; \eta) \\
&= \frac{\partial L}{\partial \eta} \delta\eta \\
\frac{\partial L}{\partial \eta} \delta\eta &= \int_{0}^{T} a_{\Omega}^{\delta\eta}(w_{\Omega}, 2u_{\Omega} - w_{\Omega}) dt.
\end{aligned} 
\quad (95) $$

## Comments about the above analysis

1. In this approach, the objective function to be minimized is $J(\eta)$, but the constitutive equation error $E(\eta)$, which contributes to $J$, is used in order to restrict the geometrical area over which a nonzero correction $\delta\eta$ is sought.

2. The error localization property of $E$ have been initially studied and applied in [20], [12] for elastic FEM model updating, allowing substantial reduction of the size of the inversion problem. To our best knowledge, it has not yet been applied to thermal inverse problems.

3. In more realistic situations where $U$ is only known over a subset of $\Omega$ and with small but not infinitesimal $1/\gamma$, the localization property is expected to hold, though obviously in an approximate manner.

4. A similar kind of error functional has also been introduced by for transient thermal inverse problems, but with no attempt to geometrically localize the conductivity error.

5. Using (95), any conventional optimization algorithm using gradients can be used. As in other cases discussed in this paper, the computation of the variation $\delta J$ of the error functional $J(\eta)$ uses two temperature fields $u_{\Omega}$, $w_{\Omega}$. However, their computation is somewhat more complicated due to their coupling through eqns. (88)-(89).

## References

[1] Balageas D.L., Deom A.A., Boscher D.M. (1987) - Characterization of Non Destructive Testing of Carbon-Epoxy Composites by a pulsed photothermal method, Materials Evaluation (45)4, pp.461-465.

[2] Banks H.T., Kojima F. (1989) - Boundary Shape Identification Problems in Two-dimensional Domains Related to Thermal Testing of Materials. Quart. Appl. Math. 47(2), pp 273-293.

[3] Bonnet M., Bui H.D., Planchard J. (1989) - Problème inverse pour l'équation de la chaleur: applications au contrôle non destructif thermique; Report Electricité de France no. HI-70-6391, pp.1-13, Clamart, France.

[4] Brebbia C.A., Telles J.C.F., Wrobel L.C. (1984) - Boundary Element Techniques. Theory and Application in Engineering. Springer - Verlag, 1984.
[5] Carslaw H.S., Jaeger J.C. (1973) - Conduction of heat in solids; Oxford University Press, Oxford, UK.

[6] Calderon A.P. (1980) - On an inverse boundary value problem; Seminar on Numerical Analysis and its applications to Continuum Physics, Soc. Brasilian de Matematica, Rio de Janeiro, pp.65-73.

[7] Connolly T.J., Wall D.J.N. (1988) - On an Inverse Problem, With Boundary Measurements, for the Steady State Diffusion Equation. Inverse Problems, 4 pp 995-1012.

[8] Friedman A., Vogelius M. (1989) - Determining cracks by boundary measurements; Indiana Univ. Math. J.; 38, pp.527-556.

[9] Kohn R., Vogelius M; (1984) - Determining conductivity by boundary measurements; Comm. Pure Appl. Math., 37, pp.289-298.

[10] Isaacson D., Isaacson E.L. (1989) - Comment on Calderon's paper "On an Inverse Boundary Value Problem". Math. Comput. 52, pp.553-559.

[11] Kohn R., McKenney A. - Numerical implementation of a variational method for electric impedance tomography. Inverse Problems, 6 pp 389-414, 1990.

[12] Ladevèze P., Reynier M., Nedjar D. - Parametric correction of finite element models using modal tests. *IUTAM Symposium on Inverse Problems in Engineering Mechanics (Tokyo, 11-15 may 1992)*, H.D. Bui & M. Tanaka, eds., Springer-Verlag.

[13] Lavrentiev M.M. (1967) - Some improperly posed problems of mathematical physics, Springer-Verlag.

[14] Lions J.L. (1968) - Contrôle optimal de systèmes gouvernés par des équations aux dérivées partielles, Dunod, Paris.

[15] Lund J., Vogel C.R. (1990) - A fully-Galerkin approach for the numerical solution of an inverse problem in a parabolic partial differential equation. Inverse Problems, 6 pp. 205-217.

[16] Marchuk G.I. (1982) - Methods of numerical mathematics, chapter 7 of *Numerical methods for some inverse problems*, pp 312-351, Springer Verlag.

[17] Menke W. - Geophysical data analysis : discrete inverse theory. Academic Press, 1984.

[18] H.G. Natke (1992) - On Regularization Methods within System Identification, In *IUTAM Symposium on Inverse Problems in Engineering Mechanics*, H.D. Bui & M. Tanaka, eds., Springer-Verlag.

[19] Petryk H., Mroz Z. (1986) - Time derivatives of integrals and functionals defined on varying volume and surface domains. Arch. Mech. 38(5-6), pp.697-724.

[20] Reynier M. - Sur le contrôle de modélisations éléments finis: recalage à partir d'essais dynamiques. PhD thesis, Ecole Normale Supérieure de Cachan, France, 1990.

[21] Tarantola A. (1987) - Inverse problem theory. Elsevier.

[22] Tikhonov A.N., Arsenin V.Y. (1977) - Solutions to ill-posed problems, Winston Wiley, New York, 1977.

[23] Vogel C.R. (1987) - An overview of numerical methods for nonlinear ill-posed problems, in *Inverse and ill-posed problems*, H.W. Engl and C.W. Groetsch, eds., Academic Press.