Probing for inclusions in heat conductive bodies

This work deals with an inverse boundary value problem arising from the equation of heat conduction. Mathematical theory and algorithm is described in dimensions 1–3 for probing the discontinuous part of the conductivity from local temperature and heat flow measurements at the boundary. The approach is based on the use of complex spherical waves, and no knowledge is needed about the initial temperature distribution. In dimension two we show how conformal transformations can be used for probing deeper than is possible with discs. Results from numerical experiments in the one-dimensional case are reported, suggesting that the method is capable of recovering locations of discontinuities approximately from noisy data.


Introduction.
1.1. Inverse heat conductivity problem. Let Ω be a bounded open set in R d with Lipschitz boundary ∂Ω, and consider the following boundary value problem   where γ(x) ∈ L ∞ (Ω) such that γ(x) > c for a constant c > 0. Let v f be the unique solution to the above equation. The time-dependent Dirichlet-to-Neumann (DN) map Λ is then defined by where ν is the outer unit normal to ∂Ω.
Physically, we consider a heat-conducting body modelled by the set Ω and the strictly positive heat conductivity distribution γ inside the body. The function v 0 is the initial temperature distribution in Ω over which we do not have control. We perform boundary measurements by applying the temperature f (x, t) at the boundary ∂Ω during the time 0 < t < T and measuring the resulting heat flux γ(x) ∂v f ∂ν | ∂Ω through the boundary. The DN map Λ defined in (2) is an ideal model of all possible infinite-precision measurements of the above type.
We study the ill-posed inverse problem of detecting conductivity inclusions inside Ω from the local knowledge of Λ. Namely, we assume that only a part Γ ⊂ ∂Ω is available for measurements and take as data the restrictions γ(x) ∂v f ∂ν | Γ of heat fluxes corresponding to functions f that are supported in Γ. An inclusion is defined as a subdomain Ω 1 ⊂ Ω such that γ(x) is perturbed on Ω 1 from some known γ 0 (x). Our aim is to find the location of ∂Ω 1 from the local measurements. The typical situation is drawn in the figure 1. The above inverse boundary value problem is related to nondestructive testing where one looks for anomalous materials inside a known material. One such example is monitoring a blast furnace used in ironmaking: the corroded thickness of the accreted refractory wall based on temperature and heat flux measurement on the accessible part of the furnace wall [11].
We remark that although Λ depends on both f and v 0 , our method uses no information of the initial data v 0 to detect ∂Ω 1 .
Moreover, we shall assume that γ 1 (x) − γ 0 (x) has a constant sign on Ω 1 . Our main purpose is to study discontinuous perturbations, however, we allow γ(x) to be continuous. Hence, we impose the following assumption.

PROBING FOR INCLUSIONS IN HEAT CONDUCTIVE BODIES
We extend γ 0 (x) smoothly outside Ω so that γ 0 (x) = 1 for large |x|. First, we consider the 1-dimensional problem. Namely, for the heat equation on (a, b), we try to detect dist(a, ∂Ω 1 ) from the measurement at a.
Take x 0 ∈ R arbitrarily. Then there exists a real fucntion ϕ λ (x) ∈ C ∞ (R) depending on a large parameter λ > 0 having the following properties.
(1) It satisfies (3) Take 0 < T 1 ≤ T arbitrarily, and let where y(x) is defined by For higher dimensions, we need the following notations. By dist(x 0 , A), we mean the distance between a point x 0 ∈ R d and a set A we take an open subset Γ on the boundary ∂Ω, and try to recover the location of ∂Ω 1 from a measurement on Γ. Take x 0 ∈ R d arbitrarily. The roles of half-lines (−∞, x 0 ) and (x 0 , ∞) are played by B(x 0 , R) and B(x 0 , R) c , respectively.
and fix constants T 1 and µ such that Then there exists ϕ λ (x) ∈ C ∞ (R d \ {x * }) depending on a large parameter λ > 0 having the following properties.
(1) It satisfies (2) For any compact sets dσ(x) being the induced measure on ∂Ω. Then we have In Theorems 1.1 and 1.2, ϕ λ (x) is constructed by using γ 0 (x) only. We put Let us remark that Then by Theorem 1.2, if ±(γ 1 − γ 0 ) > 0 on Ω 1 and R < R 1 , then ±I(λ) tends to 0 exponentially as λ → ∞, and if R > R 1 , then ±I(λ) tends to ∞ exponentially as λ → ∞. This means, we can see whether B(x 0 , R) touches Ω 1 or not using only the knowledge of γ 0 (x). Our results also hold for 4 ≤ d ≤ 6. However, to fix the idea, we shall consider the case 1 ≤ d ≤ 3.

Detection algorithm.
Suppose for the sake of simplicity that γ 0 ≡ 1. The 1-dimensional case is rather easy. We have only to use (5) to compute a 1 . The detection algorithm for d = 2, 3 is as follows. 1. Take an open set Γ ⊂ ∂Ω, and x * outside Ω, but close to Γ. 2. Draw a straight line l with end point x * , take x 0 ∈ l, and the ball B = B(x 0 , R) such that (B ∩ ∂Ω) ⊂ Γ. 3. For large λ > 0, compute I(λ). The function will give an approximation of the probing data. 4. If I(λ) → 0 as λ → ∞ we infer that B does not intersect the inclusion. 5. If I(λ) → ∞ (or −∞) as λ → ∞ we infer that B intersects the inclusion and 1.4. Transformation of exponentially growing solutions. The main idea of the proof (for d = 2, 3) consists in using the probing data of the form e λt ϕ λ (x), where ϕ λ (x) is exponentially growing as λ → ∞ in the ball B(x 0 , R) and exponentially decaying outside B(x 0 , R). Such a method was found and used in [4] for the case of the elliptic problem: ∇ · (γ(x)∇u) = 0 by passing through the hyperbolic space. The essential feature of this idea is to use conformal transformation which maps a sphere to a plane. In this paper, we shall work entirely in the Euclidean space, and use inversion with respect to a sphere (this is also the case in [4]), which maps the ball B(x 0 , R) to the half-space {y · y * < 0}. We are then led to consider the equation of the form where ζ ∈ C d . Although in our case q λ (x) grows up linearly in λ, suitable change of variables enables us to reduce the construction of solutions to now standard method of Sylvester-Uhlmann [9]. Another important difference is that in [4], the point x 0 , the center of the ball B(x 0 , R), is taken outside the covex hull of Ω, while in our method it should satisfy (C-1) and that x * ∈ ∂B \ Ω. This makes it possible to probe regions deeper than that of [4].
It is interesting to find other conformal transformations mapping planes (lines) to some surfaces (curves) which are useful for the probing problem. Putting trivial transformations such as translation and orthogonal transformation aside, for d = 3, the spherical inversion is essentially the unique conformal transformation for the Euclidean Laplacian. It can be seen by embedding the problem into the hyperbolic space, as was done in [4]. However, in 2-dimensions, the complex function theory provides us with lots of conformal transformations. We shall give in this paper one of such examples and use it for the inclusion detection problem in §5.
1.5. Literature review. There are lots of works on inverse problem of heat conductivity. In a recent article of [6], the following heat equation is considered for ν being the unit normal on ∂D. Assuming that D and ρ are unknown, they recover the location and shape of D from the knowledge of u(x, t) and ∂u(x, t)/∂ν on ∂Ω.
The proof is based on the exponentially growing solution of the form exp(τ x · (ω + iω ⊥ )) with ω, ω ⊥ ∈ S 2 , ω ⊥ ω ⊥ . This is an extreme case of the problem treated in this paper by taking the heat conducting coefficient γ(x) = 0 on Ω 1 . In [2] the interior boundary curve of an arbitrary-shaped annulus is reconstructed from overdetermined Cauchy-data on the exterior boundary curve, by assuming γ = 1 and using the Newton method to a boundary integral equation approach. In [1] the shape of an inaccessible portion of the boundary is reconstructed by linearisation. In [3] the position of the boundary is reconstructed as a function of time, using the sideways heat equation. The proposed method of the present paper differs from these; it is a direct reconstruction method, and it finds inclusions in non-constant background conductivities. For other related results, see [11,12,8] and the survey article [10].
The paper is organized as follows. In §2 and §3, we shall constuct the probing data ϕ λ (x). Theorems 1.1 and 1.2 are proved in §4. The use of conformal mapping in 2-dimensions is explained in §5. In §6, we give some numerical results for the one dimensional case.
where y(x) is defined by (6) and 3. Multi-dimensional trial function. The construction of ϕ λ (x) is more complicated in multi-dimensions than in the one-dimensional case. We shall prove the following theorem.
and let ζ ∈ C d be a complex vector such that where µ > 0 is a constant satisfying the condition (8). Let y = y(x) be the inversion defined by Then for large λ > 0, there exists a solution ϕ λ (x) to the equation (9) in the region 0 < |x − x * | < 2 √ λR having the following form: where for any for a constant C δ > 0, which also depends on R and x 0 but is independent of µ and large λ. Moreover

Schrödinger equation and spherical inversion.
To prove Theorem 3.1, we need two changes of (in)dependent variables. For a solution ϕ λ to (9), we put Next we pass to the inversion (16). Letting y ′ = y − (y · y * )y * , and using The inverse map : y → x is given by Direct computation yields the following lemma.
For a solution v(x) to (20), we put w(y) = |y − y * | 2−d v(x(y)), where x(y) is given by (22). Then by the above lemma, w(y) satisfies (−∆ y + q λ (y)) w(y) = 0, 3.2. Proof of Theorem 3.1. Taking ζ from (15), we look for a solution w(y) of (23) of the form w(y) = (1 + φ(y))e −ζ·y . Since This is the equation treated in [9], with the difference that the potential q λ grows up linearly in λ. The remedy is to make the following change of variables and parameter : Letting Φ(Y ) = φ(y), we then have

PROBING FOR INCLUSIONS IN HEAT CONDUCTIVE BODIES 9
The potential Q λ (Y ) is now bounded in λ, however, it has a singularity at Y = 0. We take We prove that if |η| = µ √ λ is sufficiently large, the equation admits a unique solution Ψ ∈ L 2,σ , where the weighted L 2 -space L 2,σ is defined by with obvious norm.
The following lemma can be proved easily, and the proof is omitted.
where the constant C does not depend on µ or on λ ≥ 1.
The above lemma implies which yields the following corollary.
Thus, w(t) = e −tAD (v 0 − u), where A D denotes the self-adjoint operator A with the homogeneous Dirichlet boundary condition. Therefore, for any t > 0 and m ≥ 0, w(t) ∈ D((A D ) m ). By the regularity theorem for elliptic operators, we then have Let φ be the solution of Aφ = 0 on Ω, φ = u 0 on ∂Ω. Then φ H 1 (Ω) ≤ C u 0 H 1/2 (∂Ω) and u − φ satisfies where C is independent of λ > 0. Hence we obtain for 0 < t ≤ T and m ≥ 0

and for a bounded open set
Note that a(O) = a B (O) for d = 2, 3 by (7) and (21). We put Let ϕ λ (x) be as in Theorem 2.1 or Theorem 3.1. Thanks to (13) and (17), we have We then put (2) We have lim Taking the logarithm, and letting κ → ∞, we get (38). Let O be a bounded open set such that O ⊂ Ω 1 . Then in view of (37), we have which implies by (38) Since O ⊂⊂ Ω 1 can be taken arbitrarily, we get (39).
We put By Corollary 4.2 and (39), we have Before entering into the probing problem, we need to study geometric properties of the map : x → y 1 (x). The next lemma follows from a direct computation using R = |x 0 − x * | and (21).

2-dimensional conformal map.
In this section, we identify x = (x 1 , x 2 ) ∈ R 2 with a complex number z = x 1 +ix 2 ∈ C. Consider a univalent holomorphic function F (z) defined in a neighborhood of Ω, and the conformal map : z → w = y 1 + iy 2 = F (z). As in §3, we look for a solution of (9) in the form ϕ λ (x) = (1 + φ(y))e −ζ·F (x) , where φ satisfies (24). Note that ζ · F (z) is identified with µλ √ 2 F (z). Assume that 0 ∈ Ω and let us choose F (z) = z −n − ρ −n with ρ > 0 and n = 1, 2, . . .. The case n = 1 corresponds to the inversion we have studied in §3, where R corresponds to 2/ρ. Consider the curve Re F (z) = 0 that separates the set of x where ϕ λ (x) has an exponential growth from the set of x where ϕ λ (x) has an exponential decay as λ → ∞. Writing z = re iα , we have that Re F (z) = 0 ⇐⇒ r = ρ(cos(nα)) 1/n . It appears that if Ω is convex, then n = 3 is a better choice than n = 1, 2 since the curve Re F (z) = 0 goes more deeply inside Ω (cf. figure 3). 6. Numerical tests in the one-dimensional case. We take γ 0 = 1 which implies The equation (45) is calculated with a finite λ, so we have the approximative reconstruction equation which is the more accurate the larger the parameter λ is. In practice larger values here implies very large values of f λ on the boundary, so numerical errors will be a problem. However for initial testing purposes we proceed by using λ = 1.0, 1.5, 2.0, . . . , 20.0. Also we choose a = 0, b = 1, x 0 = 0 and T = 1. Four test cases of the following heat conductivity are considered: where the values of a 1 , b 1 and γ 1 are given in table 6. These test cases are pictured in figure 4. The solutions v in (1) are calculated with the finite element method (FEM) so that there are N x = 200 points of x ∈ [0, 1] and N t = 100 points of t ∈ [0, 1]. We consider three different initial temperature distributions: v(x, 0) = 1, v(x, 0) = 2 sin(πx), v(x, 0) = 4| sin(4πx)|.
In figure 5 we display the solution function v(x, t), v(x, 0) = 1, λ = 2 calculated by FEM for the test case 1, smaller numbers of N x and N t are only used for printing purposes.  Table 1. Parameters a 1 , b 1 and γ 1 for the test cases To simulate measurement noise we define the vector where L is the vector containing the values of Λf λ (T 1 , 0) = γ(x) ∂v f λ ∂ν x=0 at the FEM mesh points and r is a normally distributed random vector (white noise) of the same size and the parameter c is adjusted separately for each λ so that the relative noise is approximately two percent: The logarithm of the indicator function (44) is calculated for each t ∈ [0, 1] and λ = 1.0, 1.5, 2.0, . . . , 20.0. This information can be presented in several ways. In figure 6 one can see how the logarithm of the indicator function, noisy and nonnoisy, behaves as a function of T 1 , for λ = 1 and λ = 10, in test case 1. We conclude that the value λ = 1 is too small to get accurate indicator function values, and that λ = 10 provides good values throughout the time interval [T /2, T ]. In figure 7 we have a similar graph for λ = 10, but for each initial temperature distribution. We conclude that the initial distribution v(x, 0) does not change the indicator function values in the interval [T /2, T ] and thus it will not affect the reconstructions. In figure 8 one can see how the logarithm of the indicator function, noisy and nonnoisy, behaves as a function of λ, for several choices of T 1 , in test case 1. Notice that in these graphs the analytical value based on (46) is also shown.
Based on the observation that the values of the indicator function are good between the time interval [T /2, T ] we use the following indicator function value: where I(T 1 , λ) is defined by (44) and N ⋆ t is the number of FEM mesh points t ∈ [0, 1] between the time values T /2 and T . The value I ⋆ (λ) is then inserted into (46), from which the reconstruction a 1 is calculated with λ = 1.0, 1.5, 2.0, . . . , 20.0. In figure 9 the reconstructed a 1 , noisy and non-noisy, is shown as a function of λ, for each test case.