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

Authors: Nilson C. Roberty; Denis M. de Sousa; Marcelo L. S. Rainha
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

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 problem with methods 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 one-dimensional implementations. The problem of centroid determination is addressed and questions are discussed from an computational points of view.


---

## Body

## 1. 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 three-dimensional reconstruction of sources are presented and discussed in Section 6. We conclude by pointing out the advances introduced by the present work.


## 2. Direct Transient Heat Equation Problem in Cartesian Coordinates

ByΩ⊂ℜ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 time-space ℜ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 time-space 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 (2.1)∂Q=Σ¯∪Ω0∪ΩT,
where Ω0 and ΩT are, respectively, the cylinder bottom and top sections. At cylinder top and bottom there exist the corners Γ0=Ω0¯∩Σ¯⊂ℜd-1 and ΓT=ΩT¯∩Σ¯⊂ℜd-1, respectively.The direct transient heat source initial boundary value problem consists in findingu(t,x) with (t,x)∈Q given a boundary input g(t,x) with (t,x)∈Σ¯, an initial input u0(x) with (t,x)∈Ω0, and a source distribution f(t,x) with (t,x)∈Q that verifies the problem(2.2)(Pu0,g,f){∂tu-Δu=finQ,u=u0inΩ0,u=gon⁡Σ
and Dirichlet data compatibility condition, u0=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∈Rd,where d=1,2,3 is the physical domain, and the Dirichlet boundary condition is homogeneous, g=0, the transient heat problem Pg,f (2.2) has an explicit Fourier sine solution:
(2.3)u(t,x1,…,xd)=∑i=1d∑ni=1Nicn1⋯nd(t)∏i=1dsin⁡(niπxi),
where
(2.4)cn1⋯nd(t)=exp⁡(-tπ2∑i=1dni2)∫01⋯∫01u0(x1,…,xd)∏i=1dsin⁡(niπxi)dx1⋯dxd+∫0texp⁡(-(t-τ)π2∑i=1dni2)∫01⋯∫01χω(τ)(x1,…,xd)∏i=1dsin⁡(niπxi)dx1⋯dxddτ.


### 2.1. The Adjoint Transient Heat Problem

The adjoint transient heat problem has a straightforward definition(2.5)(PvT,g,f*){-∂tv-Δv=finQ,v=vTinΩT,v=gon⁡Σ
and Dirichlet data compatibility condition, vT=g, at the time-space cylinder corner ΓT. The time reversal operator (2.6)κT:Hr,s(Q)⟶Hr,s(Q);v(t,x)⟼κ[v](t,x)=v(T-t,x)
can be used to change the changes of variables u*(t,x)=v(T-t,x) and convert the adjoint problem into an equivalent direct problem.Remark 2.2.
Note that reciprocally solutions to the direct problem (2.3) can be converted into solutions to the adjoint problem:
(2.7)v(t,x1,…,xd)=∑i=1d∑ni=1Nicn1…nd(T-t)∏i=1dsin⁡(niπxi).Definition 2.3 (consistent Cauchy datum).
By Consistent Cauchy datum associated with problem (2.2) one means the functions
(2.8)(u0,g,uT,gν)∈(γ0,γ,γT,γ1)[Hr,s].
Now, since for Cartesian exact problems Cauchy data are consistent, they are in the range of the trace operators, and the nonhomogeneous problem will be always well posed.Definition 2.4 (extended Dirichlet-to-Neumann map).
One calls the extended Dirichlet-to-Neumann map for the problem (2.2) the function defined by
(2.9)ΛΩ,Σf[(u0,g)]=(u|ΩT,∂νu|∂Ω),
when u∈H2r+2,r+1(Q) is the solution of problem (2.2) with initial data (u0,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].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(2.10)ΛΩ,Σf[u0,g]=(γT,γ1)S[u0,g,f]=(γT∘S,γ1∘S)[u0,g,f]=(uT,gν).


## 2.1. The Adjoint Transient Heat Problem

The adjoint transient heat problem has a straightforward definition(2.5)(PvT,g,f*){-∂tv-Δv=finQ,v=vTinΩT,v=gon⁡Σ
and Dirichlet data compatibility condition, vT=g, at the time-space cylinder corner ΓT. The time reversal operator (2.6)κT:Hr,s(Q)⟶Hr,s(Q);v(t,x)⟼κ[v](t,x)=v(T-t,x)
can be used to change the changes of variables u*(t,x)=v(T-t,x) and convert the adjoint problem into an equivalent direct problem.Remark 2.2.
Note that reciprocally solutions to the direct problem (2.3) can be converted into solutions to the adjoint problem:
(2.7)v(t,x1,…,xd)=∑i=1d∑ni=1Nicn1…nd(T-t)∏i=1dsin⁡(niπxi).Definition 2.3 (consistent Cauchy datum).
By Consistent Cauchy datum associated with problem (2.2) one means the functions
(2.8)(u0,g,uT,gν)∈(γ0,γ,γT,γ1)[Hr,s].
Now, since for Cartesian exact problems Cauchy data are consistent, they are in the range of the trace operators, and the nonhomogeneous problem will be always well posed.Definition 2.4 (extended Dirichlet-to-Neumann map).
One calls the extended Dirichlet-to-Neumann map for the problem (2.2) the function defined by
(2.9)ΛΩ,Σf[(u0,g)]=(u|ΩT,∂νu|∂Ω),
when u∈H2r+2,r+1(Q) is the solution of problem (2.2) with initial data (u0,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].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(2.10)ΛΩ,Σf[u0,g]=(γT,γ1)S[u0,g,f]=(γT∘S,γ1∘S)[u0,g,f]=(uT,gν).


## 3. Inverse Transient Heat Equation Source Problem

The inverse source problem that we address consists in the recovery of the sourcef, knowing the extended Dirichlet-to-Neumann map ΛΩ,Σf. When r=0, the data are regular, Green's function exists, and f∈L2(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(u0,g),(uT,gν)f is to find f∈L2(Q) such that (3.1)(IP(u0,g),(uT,gν)f){(uT,gν)=ΛΩ,Σf(u0,g)}
for all given data pair (u0,g)×(uT,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 problemsPu0,g,f and Pu0,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 for f∈L2(Q) one means the application
(3.2)ΛΩ,Σf-ΛΩ,Σ0:H1(Ω0)×H3/2,3/4(Σ)⟶H1(ΩT)×H1/2,1/4(Σ).Note that the consistence of data(u0,g) is necessary to the existence of solution for the problems Pu0,g,f and Pu0,g,0.Lemma 3.2.
Letuj,j=1,2,3,… be different solutions of problem (2.2) with the same source f∈L2(Q) and different initial time and Dirichlet data (u0j,gj),j=1,2,3,…,  respectively. Then,(i)
the relative extended Dirichlet-to-Newman operatorΛΩ,Σf-ΛΩ,Σ0 is an operator whose functional value depends only on the source function f∈L2(Q) but is independent of the initial time and Dirichlet data (u0,g);(ii)
for all solutions of consistent data problemsPf,u0j,gj, j=1,2,3,…, with the same source, the source satisfies the systems of integral equations
(3.3)∫Qf(τ,ζ)(G(T,x,τ,ζ),∂G(t,x,τ,ζ)∂ν(t,x))dζdτ=(ΛΩ,Σf-ΛΩ,Σ0)[u0j,gj]=ΛΩ,Σf[0,0],
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 functionsv in
(3.4)H-∂t-Δ2,1(Q)={v∈H2,1(Q)∣-∂tv-Δv=0}
the source f(t,x) satisfies the transient heat reciprocity gap equation
(3.5)∫Qfvdxdt=-∫ΩTΛΩ,•f[0,0]γT[v]dx-∫ΣΛ•,Σf[0,0]γ[v]dσ(t,x).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(3.6)(uT,∂ν(t,x)u)=ΛΩ,Σf[u0,g]=ΛΩ,Σf[0,0]
corresponding to some specific initial Dirichlet data (u0,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 functionsv in
(3.7)H-∂t-Δ2,1(Q)={v∈H2,1(Q)∣-∂tv-Δv=0}∧{γT[v]=0}
the source f(t,x) satisfies the transient heat lateral reciprocity gap equation
(3.8)∫Qfvdxdt=-∫ΣΛ•,Σf[0,0]γ[v]dσ(t,x).With (3.8) we can pose the problem of reconstructing the source using only lateral Cauchy data.


## 4. 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(4.1)f(t,x)=χω(t)(x)inQ,
where ω(t), t∈[0,T], is a representation of the star-shaped source boundary. For one-dimensional 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 {t0,t1,t2,…,tn,tn+1,…tN} be the knots of this partition, with t0=0 and tN=T. For tn<t<tn+1,n=0,1,N-1 we use the θ-scheme approach, [1] for the discretization of (2.2).By denotingun+1(x) with x∈Ω, the approximate solution at the time step tn+1, the transient system (2.2) is approximated by the following sequence of stationary problems: (4.2)(Hχωn+1){-Δun+1+λun+1=fn+χω(tn+1)inΩ,un+1=0on⁡Γ,withgν(tn+1):=∂νun+1on⁡Γ,
for n=0,1,2,…,N. Here λ=1/τθ and (4.3)fn=un+τ(1-θ)Δun+τθχω(tn)(x)τθ.
Note that Δun+χω(tn)=∂un/∂t and that the initial Poisson problem determining the u0 and χω(0) is (4.4)(Hg,χω0){-Δu0=χω(0)(x)inΩ,u0=0on⁡Γ,withgν(0):=∂νu0on⁡Γ.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 χω(tn)(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.


### 4.1. Iterative Source Reconstruction Scheme

The source at timetn may be further calculated as (4.5)fn=λθun-1-θθfn-1,forn=1,2,…
with f-1=0 and f0=λu0. 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, fn and χω(tn+1). Let yn+1,n=-1,0,1,…, be a solution of(4.6)(Hfnn+1){-Δyn+1+λyn+1=fninΩ,yn+1=0on⁡Γ,withgyν(tn+1):=∂νyn+1on⁡Γ,
and let wn+1,n=-1,0,1,…, be solution of (4.7)(Hχωn+1){-Δwn+1+λwn+1=χω(tn+1)(x)inΩ,wn+1=0on⁡Γ,withgwν(tn+1):=∂νwn+1on⁡Γ.Then, by the superposition principle, the solution of (4.2) is un+1=wn+1+yn+1 and the Neumann data will be the sum of the decomposed parts gν(tn+1)=gwν(tn+1)+gyν(tn+1). The Y-problems Equations (4.6) form a discrete sequence of problems with continuous source fn that may be solved before the time increment at tn begins. Its normal derivatives may be calculated (4.8)gyν(tn+1):=∂νyn+1on⁡Γ,
and subtracted from the synthetic transient Neumann data at knot tn+1(4.9)gwν(tn+1)=gν(tn+1)-gyν(tn+1),
to produce the data for the modified Helmholtz equations (4.7) that will be used in the reconstruction of the source χω(tn+1) at time tn+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.


### 4.2. 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. Letv be the space of Helmholtz functions in (4.10)Hλ2(Ω)={v∈H2(Ω):-Δv+λv=0}.
The reciprocity gap functional, [9] for the Cauchy data in the sequence of Helmholtz problem equations (4.2) is (4.11)Rfn+χω(tn+1)λ(v)=∫Ωfnvdx+∫Ωχω(tn+1)vdxforv∈Hλ2(Ω).
It is a direct consequence of Green's theorem that (4.12)Rfn+χω(tn+1)λ(v)=-∫Γvgν(tn+1)dσforv∈Hλ2(Ω).
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 fn giving (4.13)Rχω(tn+1)λ(v)=∫Ωχω(tn+1)vdx=-∫Γvgw,ν(tn+1)dσforv∈Hλ2(Ω).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:
(4.14)Q=(tn+1,tn+1+τ)×Ω
for a small time increment τ and consider that a field v∈H-∂t-Δ2,1(Q)∧γtn+1+τ[v]=0. The equation in the definition of H-∂t-Δ2,1(Q) is then integrated in this time interval:
(4.15)-∂tv-Δv=0⟹v(tn+1,⋅)-Δ∫tn+1tn+1+τv(t,⋅)dt=0.
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 (4.16)∫tn+1tn+1+τv(t,⋅)dt=θτv(tn+1,⋅).
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 H-∂t-Δ2,1(Q)∧γtn+1+τ[v]=0. The same averaging process can be done with the transient reciprocity gap equation (3.5):
(4.17)∫Ω∫tn+1tn+1+τfvdtdx=-∫Γ∫tn+1tn+1+τΛ•,Σf[0,0]γ[v]dtdσx.
With
(4.18)∫tn+1tn+1+τf(t,⋅)v(t,⋅)dt≈θ1τf(tn+1,⋅)v(tn+1,⋅),∫tn+1tn+1+τΛ•,Σf[0,0]γ[v]dt≈θ2τΛ•,Σf[0,0](tn+1,⋅)v(tn+1,⋅)
we reproduce an equation that approximates (4.13): (4.19)∫Ωf(tn+1,x)v(tn+1,x)dx≈-θ2θ1∫ΓΛ•,Σf[0,0](tn+1,x)v(tn+1,x)dσx.
Note that the data preprocessing for the determination of Λ•,Σf[0,0](tn+1,.) involves the time homogenization of initial values and is similar to that used for the determination of gw,ν(tn+1,·).


## 4.1. Iterative Source Reconstruction Scheme

The source at timetn may be further calculated as (4.5)fn=λθun-1-θθfn-1,forn=1,2,…
with f-1=0 and f0=λu0. 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, fn and χω(tn+1). Let yn+1,n=-1,0,1,…, be a solution of(4.6)(Hfnn+1){-Δyn+1+λyn+1=fninΩ,yn+1=0on⁡Γ,withgyν(tn+1):=∂νyn+1on⁡Γ,
and let wn+1,n=-1,0,1,…, be solution of (4.7)(Hχωn+1){-Δwn+1+λwn+1=χω(tn+1)(x)inΩ,wn+1=0on⁡Γ,withgwν(tn+1):=∂νwn+1on⁡Γ.Then, by the superposition principle, the solution of (4.2) is un+1=wn+1+yn+1 and the Neumann data will be the sum of the decomposed parts gν(tn+1)=gwν(tn+1)+gyν(tn+1). The Y-problems Equations (4.6) form a discrete sequence of problems with continuous source fn that may be solved before the time increment at tn begins. Its normal derivatives may be calculated (4.8)gyν(tn+1):=∂νyn+1on⁡Γ,
and subtracted from the synthetic transient Neumann data at knot tn+1(4.9)gwν(tn+1)=gν(tn+1)-gyν(tn+1),
to produce the data for the modified Helmholtz equations (4.7) that will be used in the reconstruction of the source χω(tn+1) at time tn+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.


## 4.2. 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. Letv be the space of Helmholtz functions in (4.10)Hλ2(Ω)={v∈H2(Ω):-Δv+λv=0}.
The reciprocity gap functional, [9] for the Cauchy data in the sequence of Helmholtz problem equations (4.2) is (4.11)Rfn+χω(tn+1)λ(v)=∫Ωfnvdx+∫Ωχω(tn+1)vdxforv∈Hλ2(Ω).
It is a direct consequence of Green's theorem that (4.12)Rfn+χω(tn+1)λ(v)=-∫Γvgν(tn+1)dσforv∈Hλ2(Ω).
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 fn giving (4.13)Rχω(tn+1)λ(v)=∫Ωχω(tn+1)vdx=-∫Γvgw,ν(tn+1)dσforv∈Hλ2(Ω).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:
(4.14)Q=(tn+1,tn+1+τ)×Ω
for a small time increment τ and consider that a field v∈H-∂t-Δ2,1(Q)∧γtn+1+τ[v]=0. The equation in the definition of H-∂t-Δ2,1(Q) is then integrated in this time interval:
(4.15)-∂tv-Δv=0⟹v(tn+1,⋅)-Δ∫tn+1tn+1+τv(t,⋅)dt=0.
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 (4.16)∫tn+1tn+1+τv(t,⋅)dt=θτv(tn+1,⋅).
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 H-∂t-Δ2,1(Q)∧γtn+1+τ[v]=0. The same averaging process can be done with the transient reciprocity gap equation (3.5):
(4.17)∫Ω∫tn+1tn+1+τfvdtdx=-∫Γ∫tn+1tn+1+τΛ•,Σf[0,0]γ[v]dtdσx.
With
(4.18)∫tn+1tn+1+τf(t,⋅)v(t,⋅)dt≈θ1τf(tn+1,⋅)v(tn+1,⋅),∫tn+1tn+1+τΛ•,Σf[0,0]γ[v]dt≈θ2τΛ•,Σf[0,0](tn+1,⋅)v(tn+1,⋅)
we reproduce an equation that approximates (4.13): (4.19)∫Ωf(tn+1,x)v(tn+1,x)dx≈-θ2θ1∫ΓΛ•,Σf[0,0](tn+1,x)v(tn+1,x)dσx.
Note that the data preprocessing for the determination of Λ•,Σf[0,0](tn+1,.) involves the time homogenization of initial values and is similar to that used for the determination of gw,ν(tn+1,·).


## 5. Determining Centroid and Shape

The necessary functions for definition of source centroid, constant function and position function{1,xi,i=1,…,d}∈H-∂t-Δ2,1(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(5.1)xi¯=∫Qfxidxdt∫Qf1dxdt,i=1,…,d.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.2)xi¯=∫ΩTΛΩ,•f[0,0]xidx+∫ΣΛ•,Σf[0,0]xidσ(t,x)∫ΩTΛΩ,•f[0,0]dx+∫ΣΛ•,Σf[0,0]dσ(t,x).A problem appears when we work with the modified Helmholtz approximation (4.10), since functions {1,xi,i=1,…,d}∉Hλ2(Ω)={v∈H2(Ω):-Δv+λv=0}. In this case, if we first define κ=(λ)=(1/θτ), rooted parameter in the auxiliary modified Helmholtz parameter, we have a special set of Hλ(Ω) functions, that is, (5.3){exp⁡(κ(∑i=1dlixi))}for(l1,…ld)∈Sd-1.
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 lj∈Sd-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:
(5.4){sinh⁡(κ(∑i=1dli(xi-x̅i)))κ;cosh⁡(κ(∑i=1dli(xi-x̅i)))}for(l1,…ld)∈Sd-1
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,(x1,…,xd)} 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.Definition 5.3 (meta centroid).
Letω⊂Ω∈Rd.By meta centroid x̅=(x̅1,…,x̅d) of this subdomain one means
(5.5)x̅i=∫ωxi(sinh⁡(κ(xi-x̅i)))/(κ(xi-x̅i))dx∫ω(sinh⁡(κ(xi-x̅i)))/(κ(xi-x̅i))dx,fori=1,…,d.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 functionsinh 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.


### 5.1. 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 coordinatesx̅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
(5.6)x̅iκ=argmin⁡{|∫Γsinh⁡(κ(xi-xic))κgνdσ(x)|2∣xc∈Ω}fori=1,…,d.


### 5.2. Determining Shape Parameters

Once we have reconstructed the meta centroid, we may proceed with the shape parameter determination with the same modified Helmholtz data:(5.7)ω̅=argmin⁡{|∫ωcosh⁡(κ(∑i=1dli(xi-x̅iκ)))dx-IΓ(cosh⁡,κ,l,x̅iκ,gν)|2:ω⊂Ω;l∈Sd-1},
where the set of directions l:=(l1,…ld)∈Sd-1 is used to generate linearly independent functionals of the trial shape and (5.8)IΓ(cosh⁡,κ,l,x̅iκ,gν):=∫Γcosh⁡(κ(∑i=1dli(xi-x̅iκ)))gν(x)dσ(x)
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 thecosh 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
(5.9)|∫ωcosh⁡(κ(xi-x̅i(tn+1)))dx-∫Γcosh⁡(κ(xi-x̅i(tn+1)))gw,ν(tn+1)dσ(x)|2
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.


## 5.1. 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 coordinatesx̅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
(5.6)x̅iκ=argmin⁡{|∫Γsinh⁡(κ(xi-xic))κgνdσ(x)|2∣xc∈Ω}fori=1,…,d.


## 5.2. Determining Shape Parameters

Once we have reconstructed the meta centroid, we may proceed with the shape parameter determination with the same modified Helmholtz data:(5.7)ω̅=argmin⁡{|∫ωcosh⁡(κ(∑i=1dli(xi-x̅iκ)))dx-IΓ(cosh⁡,κ,l,x̅iκ,gν)|2:ω⊂Ω;l∈Sd-1},
where the set of directions l:=(l1,…ld)∈Sd-1 is used to generate linearly independent functionals of the trial shape and (5.8)IΓ(cosh⁡,κ,l,x̅iκ,gν):=∫Γcosh⁡(κ(∑i=1dli(xi-x̅iκ)))gν(x)dσ(x)
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 thecosh 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
(5.9)|∫ωcosh⁡(κ(xi-x̅i(tn+1)))dx-∫Γcosh⁡(κ(xi-x̅i(tn+1)))gw,ν(tn+1)dσ(x)|2
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.


## 6. 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.


### 6.1. 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 xc=.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)=gwν,(n+1)+gyν,(n+1) 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) (6.1)Rχωλ(sinh⁡(κ(xi-x̅i(tn+1)))κ)=∫Γsinh⁡(κ(xi-x̅i(tn+1)))κgw,ν(tn+1)dσ(x)=0,
for i=1,…,d in the metacentroid calculations and as
(6.2)Rχωκ(cosh⁡(κ(xi-x̅i(tn+1))))(tn+1)=∫ωcosh⁡(κ(x-x̅i(tn+1)))dx=∫Γcosh⁡(κ(xi-x̅i(tn+1)))gw,ν(tn+1)dσ(x),fori=1,…,d
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 xc of the centroid is obtained by minimizing the following functional (6.3)|Rχωλ(sinh⁡(κ(xi-x̅i(tn+1)))κ)|2=|∫Γsinh⁡(κ(xi-x̅i(tn+1)))κgw,ν(x,tn+1)dσ(x)|2,
for i=1,…,d.Some results have been selected to show the methodology.Figure1 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.Heat flux on the transient one-dimensional model.Transient source reconstruction for times0,2, and 4.Transient interval source centroid and size: exact blue captured as red.


### 6.2. Three-Dimensional Numerical Simulations

The model case studied here is a source inside the domainΩ=(0,1)d∈ℝd with a parallelepiped block (voxel) shape. It is supported with a harmonic centroid evolution following the parametric curve (6.4)(.5(1+.6sin⁡(2πt)),.5(1+.6cos⁡(2πt)),.25+.25t)
and deforming equally in all directions by the following time rule: (6.5)hx=hy=hz=h=.15(1+.25|cos⁡(πt2)|),
where block edge is 2h. The number of harmonics in the Fourier sine series is 20, the Δx for spacial collocation is .01, and θ is chosen as .8. The evolution is calculated for various 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 Figure4, 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.Figure 4
3D metacentroid movement with time increment of.05.In Figure5, 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.3D centroid coordinates and block sizes with time increment of.05.3D centroid coordinates and block sizes with time increment of.02.In Figures7 and 8, we show the absolute error in the centroid and in the size calculation, respectively.Figure 7
3D absolute centroid calculation error with time increment of.05 sec.Figure 8
3D size error with time increment of.05 sec.Finally, the reciprocity gap functional atcosh(κ(xi-xci)) and sinh(κ(xi-xci)) for i=1,…,d is shown in Figure 9.Reciprocity gap functional atsinh and cosh with time increment of .05 sec.


## 6.1. 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 xc=.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)=gwν,(n+1)+gyν,(n+1) 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) (6.1)Rχωλ(sinh⁡(κ(xi-x̅i(tn+1)))κ)=∫Γsinh⁡(κ(xi-x̅i(tn+1)))κgw,ν(tn+1)dσ(x)=0,
for i=1,…,d in the metacentroid calculations and as
(6.2)Rχωκ(cosh⁡(κ(xi-x̅i(tn+1))))(tn+1)=∫ωcosh⁡(κ(x-x̅i(tn+1)))dx=∫Γcosh⁡(κ(xi-x̅i(tn+1)))gw,ν(tn+1)dσ(x),fori=1,…,d
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 xc of the centroid is obtained by minimizing the following functional (6.3)|Rχωλ(sinh⁡(κ(xi-x̅i(tn+1)))κ)|2=|∫Γsinh⁡(κ(xi-x̅i(tn+1)))κgw,ν(x,tn+1)dσ(x)|2,
for i=1,…,d.Some results have been selected to show the methodology.Figure1 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.Heat flux on the transient one-dimensional model.Transient source reconstruction for times0,2, and 4.Transient interval source centroid and size: exact blue captured as red.


## 6.2. Three-Dimensional Numerical Simulations

The model case studied here is a source inside the domainΩ=(0,1)d∈ℝd with a parallelepiped block (voxel) shape. It is supported with a harmonic centroid evolution following the parametric curve (6.4)(.5(1+.6sin⁡(2πt)),.5(1+.6cos⁡(2πt)),.25+.25t)
and deforming equally in all directions by the following time rule: (6.5)hx=hy=hz=h=.15(1+.25|cos⁡(πt2)|),
where block edge is 2h. The number of harmonics in the Fourier sine series is 20, the Δx for spacial collocation is .01, and θ is chosen as .8. The evolution is calculated for various 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 Figure4, 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.Figure 4
3D metacentroid movement with time increment of.05.In Figure5, 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.3D centroid coordinates and block sizes with time increment of.05.3D centroid coordinates and block sizes with time increment of.02.In Figures7 and 8, we show the absolute error in the centroid and in the size calculation, respectively.Figure 7
3D absolute centroid calculation error with time increment of.05 sec.Figure 8
3D size error with time increment of.05 sec.Finally, the reciprocity gap functional atcosh(κ(xi-xci)) and sinh(κ(xi-xci)) for i=1,…,d is shown in Figure 9.Reciprocity gap functional atsinh and cosh with time increment of .05 sec.


## 7. 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.



---

*Source: 130834-2011-10-24.xml*