Moving Heat Source Reconstruction from Cauchy Boundary Data : The Cartesian Coordinates Case

We consider the problem of reconstruction of an unknown characteristic interval and block transient thermal source inside a domain. By exploring the definition of an Extended Dirichlet to Neumann map in the time space cylinder that has been introduced in Roberty and Rainha 2010a , we can treat the problemwithmethods similar to that used in the analysis of the stationary source reconstruction problem. Further, the finite difference θ-scheme applied to the transient heat conduction equation leads to a model based on a sequence of modified Helmholtz equation solutions. For each modified Helmholtz equation the characteristic interval and parallelepiped source function may be reconstructed uniquely from the Cauchy boundary data. Using representation formula we establish reciprocity functional mapping functions that are solutions of the modified Helmholtz equation to their integral in the unknown characteristic support. Numerical experiment for capture of an interval and an rectangular parallelepiped characteristic source inside a cubic box domain from boundary data are presented in threedimensional and onedimensional implementations. The problem of centroid determination is addressed and questions are discussed from an computational points of view.


Introduction
Inverse source transient heat problem has been studied by a huge number of authors.In the companion paper 1 , we have presented a brief review on the subject.By adoption of the reciprocity gap functional method to solve an sequence of stationary source problem, we have developed the methodology for transient heat characteristic source reconstruction presented in this work.The model is based on the modified Helmholtz Poisson equation that is obtained from the transient equation through the θ-scheme associated with a time finite difference discretization.Analysis of the related mathematical and computational work involved has been presented by the authors in national conferences, Roberty et al. 2-6 .Theoretical aspects of the problem can be found in the companion paper 1 .This paper will be structured as follows.Some definitions and mathematical results extracted from 1 that are important to the understanding of the inverse problems are presented in Section 2. There we introduce the concept of consistent Cauchy data, extended Dirichlet-to-Neumann map.The inverse transient heat source problem is introduced in Section 3. A basic lemma about the relative extended Dirichlet-to-Neumann map and reciprocity gap functional in the transient model extracted from 1 are discussed.These concepts present original aspects which have been introduced in the cited companion paper and are numerically investigated in the present work.We show that the transient problem can be studied with aid of results demonstrated for the modified Helmholtz Dirichlet source problem in Section 4.There the iterative source reconstruction scheme is presented.In Section 5, some problems related with the application of the Reciprocity Gap methodology for shape and centroid determination are resolved.We must point out that the nonlinear problem of shape parameter determination is central in the present work since they will change in time with the source support evolution.Finally, the numerical results for one-and threedimensional reconstruction of sources are presented and discussed in Section 6.We conclude by pointing out the advances introduced by the present work.

Direct Transient Heat Equation Problem in Cartesian Coordinates
By Ω ⊂ R d , d 1, 2, 3 we denote a bounded space domain with Cartesian coordinates boundary Γ ∂Ω.In this case the boundary is composed of two points if d 1, four intervals if d 2, or six rectangular faces if d 3.In the spatial surface Γ the normal ν is defined almost everywhere and the induced measure on the surface is denoted by dσ.In the timespace R d 1 , we consider the time interval I : 0, T , T > 0 to form the bounded cylinder whose basis is an interval, a rectangle, or a parallelepiped, Q : I × Ω, whose lateral timespace surface is Σ : I × Γ, where Γ depends on the space where Ω is embedded.A section in this space-time cylinder is Ω t : {t} × Ω, and the complete cylinder boundary is where Ω 0 and Ω T are, respectively, the cylinder bottom and top sections.At cylinder top and bottom there exist the corners The direct transient heat source initial boundary value problem consists in finding u t, x with t, x ∈ Q given a boundary input g t, x with t, x ∈ Σ, an initial input u 0 x with t, x ∈ Ω 0 , and a source distribution f t, x with t, x ∈ Q that verifies the problem and Dirichlet data compatibility condition, u 0 g, at the time-space cylinder corner Γ 0 .
For more information about the theoretical and abstract setting and the Hilbert space formulation of the problem please see the companion paper of this work 1 .There the important theorems, lemmas, and proposition to understand the problem may be found.Major trace problems presented in the space-time cylinder formulation for the problem disappear when the problem is treated in Cartesian coordinates and spatial and temporal variables separation is straightforward.For the Hilbert space framework we need to introduce, following Lions and Magenes 7 , anisotropic Sobolev spaces.A comprehensive presentation of these spaces in the contest of boundary integral operators related with the heat equation and the heat potential can be found in Costabel 8 .But since we are restricting the experiments presented in this work to problems in Cartesian coordinates, the main trace problems do not appear and Cauchy data at the boundary will always be consistent.
Remark 2.1 solution of the direct problem by Fourier sine series .When the external domain Ω is a box 0, 1 d ∈ R d , where d 1, 2, 3 is the physical domain, and the Dirichlet boundary condition is homogeneous, g 0, the transient heat problem P g,f 2.2 has an explicit Fourier sine solution: where 2.4

The Adjoint Transient Heat Problem
The adjoint transient heat problem has a straightforward definition and Dirichlet data compatibility condition, v T g, at the time-space cylinder corner Γ T .The time reversal operator when u ∈ H 2r 2,r 1 Q is the solution of problem 2.2 with initial data u 0 , g u| Ω 0 , u| Σ .
Note that this operator can be viewed as a combination of the standard Dirichlet-to-Neumann map in the spatial boundary with the input-to-output map in the time boundary, that is, in initial and final interval times, found in control theory.For more information, please see 1 .

Inverse Transient Heat Equation Source Problem
The inverse source problem that we address consists in the recovery of the source f, knowing the extended Dirichlet-to-Neumann map Λ f Ω,Σ .When r 0, the data are regular, Green's function exists, and f ∈ L 2 Q .Let us investigate this situation.And then, we will show that the unique information available in this inverse problem is given only by one measurement, say, the bottom and top Dirichlet data and lateral cylinder Cauchy boundary data.The inverse problem IP for all given data pair u 0 , g × u T , g ν corresponding to different solutions for the direct problem.
For a specific source and an appropriate dimension, the synthetic final time and Neumann data 3.1 to be used in the reconstruction inverse problem can be calculated as the normal derivative and final time value of solution 2.3 .Definition 3.1 relative extended Dirichlet-to-Neumann map .Consider two problems P u 0 ,g,f and P u 0 ,g,0 , one with source f and the other with zero source, but both with the same consistent initial time and Dirichlet data.By the relative extended Dirichlet-to-Neumann map Note that the consistence of data u 0 , g is necessary to the existence of solution for the problems P u 0 ,g,f and P u 0 ,g,0 .

an operator whose functional value depends only on the source function
ii for all solutions of consistent data problems P f,u 0 j ,g j , j 1, 2, 3, . .., with the same source, the source satisfies the systems of integral equations

which depend only on the relative extended Dirichlet-to-Neumann map. Here G t, x, τζ is the causal Dirichlet Green function for the transient heat problem;
iii for all test functions v in the source f t, x satisfies the transient heat reciprocity gap equation Proof.See 1 .
Remark 3.3.Note that in this case the unique information available for source reconstruction is given by only one measurement, that is, the final Neumann boundary measurement corresponding to some specific initial Dirichlet data u 0 , g , which may be assumed as zero without loss of generality.
Remark 3.4.Note that functions in the test space 3.4 are solutions to problem 2.7 with domain Q * containing Q and arbitrary boundary conditions.So there are plenty of functions, regular and singular.An important subclass of these test functions are those for which the trace on Ω T is null, that is, γ T v 0. For these functions we have iii * for all test functions v in the source f t, x satisfies the transient heat lateral reciprocity gap equation With 3.8 we can pose the problem of reconstructing the source using only lateral Cauchy data.

The θ-Scheme and the Modified Helmholtz Model for the Transient Heat Problem
We present now an algorithm for moving transient source reconstruction in the heat equation based on this result.Let the source be given by where ω t , t ∈ 0, T , is a representation of the star-shaped source boundary.For onedimensional problems it is a set with two points.For two-or three-dimensional problems it is a moving Lipschitz parametric curve or surface in which the parameter has been omitted, but, in order to make the implementation simpler, we are considering that the source is a characteristic rectangular parallelepiped or a rectangle.Consider a partition of the time interval 0, T into N subintervals of length τ > 0. Let {t 0 , t 1 , t 2 , . . ., t n , t n 1 , . . .t N } be the knots of this partition, with t 0 0 and t N T .For t n < t < t n 1 , n 0, 1, N − 1 we use the θ-scheme approach, 1 for the discretization of 2.2 .
By denoting u n 1 x with x ∈ Ω, the approximate solution at the time step t n 1 , the transient system 2.2 is approximated by the following sequence of stationary problems: for n 0, 1, 2, . . ., N.Here λ 1/τθ and Note that Δu n χ ω t n ∂u n /∂t and that the initial Poisson problem determining the u 0 and

4.4
The sequence of modified Helmholtz source inverse problem equations 4.2 starting with stationary problem equations 4.4 may be used to model a scheme for the reconstruction of star-shaped sources χ ω t n x , for a time knot sequence, showing its movement and deformation in the external domain Ω.For this, we only need to know the transient Neumann data with zero Dirichlet datum on the external boundary Γ ∂Ω.Since we do not have experimental data, we will solve the direct problem with a different method, adding noise, and do an experimental data synthesis.

Iterative Source Reconstruction Scheme
The source at time t n may be further calculated as with f −1 0 and f 0 λu 0 .Since the discretized direct problem equations 4.2 and 4.4 are linear, the problem may be decomposed into two subproblems separating the known part of the source from the part to be reconstructed, that is, f n and χ ω t n 1 .Let y n 1 , n −1, 0, 1, . .., be a solution of and let w n 1 , n −1, 0, 1, . .., be solution of

4.7
Then, by the superposition principle, the solution of 4.2 is u n 1 w n 1 y n 1 and the Neumann data will be the sum of the decomposed parts g ν t n 1 g ν w t n 1 g ν y t n 1 .The Y-problems Equations 4.6 form a discrete sequence of problems with continuous source f n that may be solved before the time increment at t n begins.Its normal derivatives may be calculated to produce the data for the modified Helmholtz equations 4.7 that will be used in the reconstruction of the source χ ω t n 1 at time t n 1 .Note that by using the reciprocity gap functional the characteristic star-shaped source may be reconstructed without solving the direct problem equations 4.7 .By using the second Green's formula, this inverse problem is modeled with a nonlinear Fredholm integral equation of first kind.

Reciprocity Gap Functional for the Helmholtz Problem
The reciprocity gap functional for the Helmholtz problem depends only on boundary values of the solution, and its properties are derived from elementary properties of Green's theorem.
Let v be the space of Helmholtz functions in The reciprocity gap functional, 9 for the Cauchy data in the sequence of Helmholtz problem equations 4.2 is It is a direct consequence of Green's theorem that The combination of these equations for test functions in H λ Ω will form the nonlinear system of equations for the source reconstruction inverse problem, 9 .We may improve the implementation of the reciprocity gap functional 4.11 by subtracting the Cauchy data of the auxiliary modified Helmholtz problem equations 4.8 .This eliminates the source term f n giving We can compare this reciprocity gap equation 4.13 with the transient heat reciprocity gap equation 3.5 .For this, let us consider the following time space cylinder: for a small time increment τ and consider that a field v ∈ H 2,1 Since the interval is sufficiently small and the field is zero in its upper extremity, by the mean value theorem, we find 0 ≤ θ ≤ 1 such that Let us define λ κ 2 1/θτ.By noting the definition of the space of Helmholtz functions H 2 λ Ω given by 4.10 , we can see that it is a θ weight average of functions of the space The same averaging process can be done with the transient reciprocity gap equation 3.5 : 4.17 With we reproduce an equation that approximates 4.13 :

Mathematical Problems in Engineering
Note that the data preprocessing for the determination of Λ f •,Σ 0, 0 t n 1 , .involves the time homogenization of initial values and is similar to that used for the determination of g w,ν t n 1 , • .

Determining Centroid and Shape
The necessary functions for definition of source centroid, constant function and position function {1, x i , i 1, . . ., d} ∈ H 2,1 −∂ t −Δ Q and we can use 3.5 to define: Definition 5.1 average centroid .By the average source centroid in the time-space cylinder one means , i 1, . . ., d.

5.1
Based on this definition, we can enunciate the following trivial lemma consequence of the definition.
Lemma 5.2 relation between the extended Dirichlet-to-Neumann map and the average centroid .

5.3
These functions are in the domain of the modified Helmholtz equation and form a base for the space H 2 λ Ω .We may construct a dense set by choosing some discrete set of directions l j ∈ S d−1 appropriately.An appropriate modification of this set will be obtained by substituting these exponentials with hyperbolic functions sinh and cosh.These functions are, respectively, skew symmetric and symmetric with respect to origin of the coordinates system.If we know the star-shaped source centroid, it is best to choose the origin in the centroid and set the following basis: to have a more balanced system of test functions to use in the reciprocity gap functional equation 4.13 .As already mentioned, contrary to the classical Novikov's star-shaped source reconstruction with boundary data problem for the Laplace operator κ 0 , in which the centroid and the source volume may be obtained as zero and first-order moments of the Neumann data at the boundary, the necessary functions for centroid calculations {1, x 1 , . . ., x d } are in the space H 2 λ Ω .Fortunately, in this generic case of κ / 0 we may introduce a concept that we are naming metacentroid, κ-centroid, or λ-space centroid.It may also be estimated from Neumann data in the boundary, and in the case in which the star-shaped source is a Cartesian domain interval, rectangle, or parallelepiped rectangular voxel, that this κ-centroid is equivalent to the κ 0 centroid, that is, the harmonic centroid in Novikov's problem, in the sense that if the source domain is star-shaped with respect to one centroid, it also is star-shaped with respect to the other metacentroid.

5.5
Lemma 5.4.Suppose that the star-shaped source characteristic support border curve is symmetric with respect to the ordinates and the abscissa axis passing through the centroid.Then the metacentroid coincides with the harmonic centroid.
Proof.In fact, since the function sinh is skew symmetric, expression 5.5 will calculate zero in the coordinates system for with the harmonic source centroid is the origin.
Since in the transient problem the source is moving inside the box, which means that its centroid position may vary with time, the capacity of centroid position determination is fundamental for the solution of the source reconstruction problem.

Determining the Metacentroid
Since the Neumann data are frequently noisy, the least square nonlinear method may be used to formulate an unconstrained minimizing problem for the determination of coordinates x i of the centroid.If necessary, classical regularization methods, such as the method of Tikhonov, may be adapted for the stabilization and improvement of the algorithm.Without any regularization other than truncation, the problem of centroid determination in the modified Helmholtz equation with boundary Dirichlet data zero and g ν Neumann data on the boundary is

Determining Shape Parameters
Once we have reconstructed the meta centroid, we may proceed with the shape parameter determination with the same modified Helmholtz data: where the set of directions l : l 1 , . . .l d ∈ S d−1 is used to generate linearly independent functionals of the trial shape and may be computed by using only the just calculated metacentroid coordinates and the already known Cauchy data on the boundary.
Remark 5.5.Note also that we have shown the cosh dependency of the integral equation 5.8 to stress the fact that we may prefer or it may be more convenient to use other kind of functions such as exponentials and modified Bessel and develop other numerical schemes in the shape determination.For the Cartesian source interval, rectangle, or parallelepiped inside the unitary box, these integrals may be evaluated with a symbolic solver such as the Mathematica or the Maple.The formal solution is given by a huge combination of exponentials forming the hyperbolic functions.This function is then implemented as a least square nonlinear minimization problem to be numerically solved.The main question with this minimizing problem is the behavior of the hyperbolic function for high values of κ.If the time step is decreased in order to improve the direct problem solution, the parameter κ in the modified Helmholtz equation model for the inverse problem will increase and the least square nonlinear method will be inadequate for the centroid search.Our experiments show that problems start for κ 6.So we restrict the time increment to τ .05, in the present work, in order to avoid these problematic higher values κ.
The same treatment was adopted for 6.2 , used for the determination of thickness.The associated least square nonlinear functional to be minimized is for i 1, . . ., d.Note that, at these calculation stages, the centroid has been reconstructed for the present time increment with the minimizing problem given by the functional equatoin 5.6 .Now only parameters associated with thickness are to be reconstructed.Obviously, a good reconstruction of centroid coordinates is fundamental in the thickness determination with 5.9 .We have observed that, for the same Helmholtz parameter κ, this second reconstruction runs worse than the first one.

Numerical Simulations
We have investigated a series of numerical experiments in one, two, and three dimensions and selected some results to show the proposed methodology for the transient star-shaped source reconstruction problem.

One-Dimensional Numerical Simulations for the Transient Heat Source Problem
In the one-dimensional box 0, 1 , slab with unitary thickness, we consider an interval source evolving with time, with centroid x c .4 .2sin 2πt/T and size h .1 .05|sin 2πt/T |.The total time T 2π is incremented with τ .05,and θ .7 is chosen for the θ-scheme.The synthetic net flux data g ν, n 1 g ν, n 1 w g ν, n 1 y at the two points is 10% noisy, and is substituted in the reciprocity gap functional equatoin 4.13 with the test functions in H λ Ω chosen as in 5.6  for i 1, . . ., d in the metacentroid calculations and as for the interval size calculation.Here, as a one-dimensional model, d 1, but we left the generic expression with d arbitrary since in multidimensional case the parallelepiped voxel case is calculated in a similar way with the interval case.The computational implementation of the solution was done by the least square nonlinear method, and the root x c of the centroid is obtained by minimizing the following functional Figure 1 presents the synthetic heat flux at the one-dimensional boundary used for transient reconstruction of the moving interval source presented in three special times t 0, 2, 4, in Figure 2.For the same model, Figure 3 shows the exact and the reconstructed evolution of the centroid position and interval thickness.

Three-Dimensional Numerical Simulations
The model case studied here is a source inside the domain Ω 0, 1 d ∈ R d with a parallelepiped block voxel shape.It is supported with a harmonic centroid evolution following the parametric curve .5    time steps τ .1,.01 in the interval t ∈ 0, 2 .For these values of time increment, the modified Helmholtz equation parameter varies as κ 3, . . ., 12.As the κ value approaches to 6, the reconstruction starts to become worse, so we may here observe that, for this special set of heat equation coefficients, which means a thermal inertia equal to one, the minimum time increment for the present methodology without any kind of regularization procedure is approximately τ 0.05.
Note that we face here the classical problem associated with ill-posedness.If we try to improve the transient problem calculations by adopting a small time increment, the associate inverse source reconstruction problem starts to present error propagation and stability problems.In Figure 4, we show the moving centroid, exact blue captured as red in its spiral trajectory inside the unitary box.For time increments smaller than τ 0.05, the results become worse.
In Figure 5, we show the centroid coordinates and block size evolution; again, the exact blue is reconstructed as red.The same parameters are presented in Figure 6, with time increment τ 0.02, to put in evidence the bad behavior of this time increment.This obviously means that, if for some special reason we need to use higher values for the parameter κ in the auxiliary Helmholtz equation, a special regularization methodology must be developed.
In Figures 7 and 8, we show the absolute error in the centroid and in the size calculation, respectively.
Finally, the reciprocity gap functional at cosh κ x i − xc i and sinh κ x i − xc i for i 1, . . ., d is shown in Figure 9.

Conclusions
We have presented a methodology for star-shaped source reconstruction in the transient heat problem by using one set of Cauchy data history.With the adoption of an anisotropic Sobolev Hilbert mathematical framework, we can treat the problem with a methodology analogous to that used to study stationary elliptic problems.Therefore, by introducing a finite difference time θ-scheme, we developed an algorithm based on a modified Helmholtz system, for which we have already studied computationally the inverse source reconstruction problem.An original methodology for centroid and shape capture is introduced.Numerical experiments in Cartesian geometry involving an interval and a rectangular parallelepiped are investigated to stress-associated difficulties.

Remark 2 . 5
extended Dirichlet-to-Neumann map is composition of trace and solution .The extended Dirichlet-to-Neumann map is a composition of the final time trace and the lateral boundary normal trace with the solution operator

Figure 1 :
Figure 1: Heat flux on the transient one-dimensional model.

Figure 3 :
Figure 3: Transient interval source centroid and size: exact blue captured as red.

Figure 5 :
Figure 5: 3D centroid coordinates and block sizes with time increment of .05.

3 for i 1 ,Figure 6 :
Figure 6: 3D centroid coordinates and block sizes with time increment of .02.

Figure 9 :
Figure 9: Reciprocity gap functional at sinh and cosh with time increment of .05sec.
.8 and subtracted from the synthetic transient Neumann data at knot t n 1