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

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.

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 Ω⊂ℜ^𝑑,𝑑=1,2,3 we denote a bounded space domain with Cartesian coordinates boundary Γ=𝜕Ω. In this case the boundary is composed of two points if 𝑑=1, four intervals if 𝑑=2, or six rectangular faces if 𝑑=3. In the spatial surface Γ the normal 𝜈 is defined almost everywhere and the induced measure on the surface is denoted by 𝑑𝜎. In the time-space ℜ^𝑑+1,we consider the time interval 𝐼∶=(0,𝑇), 𝑇>0 to form the bounded cylinder whose basis is an interval, a rectangle, or a parallelepiped, 𝑄∶=𝐼×Ω, whose lateral time-space surface is Σ∶=𝐼×Γ, where Γ depends on the space where Ω is embedded. A section in this space-time cylinder is Ω_𝑡∶={𝑡}×Ω, and the complete cylinder boundary is 𝜕𝑄=Σ∪Ω₀∪Ω_𝑇, (2.1) where Ω₀ and Ω_𝑇 are, respectively, the cylinder bottom and top sections. At cylinder top and bottom there exist the corners Γ₀=Ω₀∩Σ⊂ℜ^𝑑−1 and Γ_𝑇=Ω_𝑇∩Σ⊂ℜ^𝑑−1, respectively.

The direct transient heat source initial boundary value problem consists in finding 𝑢(𝑡,𝑥) with (𝑡,𝑥)∈𝑄 given a boundary input 𝑔(𝑡,𝑥) with (𝑡,𝑥)∈Σ, an initial input 𝑢₀(𝑥) with (𝑡,𝑥)∈Ω₀, and a source distribution 𝑓(𝑡,𝑥) with (𝑡,𝑥)∈𝑄 that verifies the problem𝑃_{𝑢0,𝑔,𝑓}𝜕_𝑡𝑢−Δ𝑢=𝑓in𝑄,𝑢=𝑢₀inΩ₀,𝑢=𝑔onΣ (2.2) and Dirichlet data compatibility condition, 𝑢₀=𝑔, at the time-space cylinder corner Γ₀.

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)^𝑑∈𝑅^𝑑, where 𝑑=1,2,3 is the physical domain, and the Dirichlet boundary condition is homogeneous, 𝑔=0, the transient heat problem 𝑃_𝑔,𝑓 (2.2) has an explicit Fourier sine solution: 𝑢𝑡,𝑥₁,…,𝑥_𝑑=^𝑑𝑁_𝑖=1^𝑖_𝑛𝑖=1𝑐_{𝑛1⋯𝑛𝑑}(𝑡)^𝑑_𝑖=1sin𝑛_𝑖𝜋𝑥_𝑖, (2.3) where 𝑐_{𝑛1⋯𝑛𝑑}(𝑡)=exp−𝑡𝜋^2𝑑_𝑖=1𝑛²_𝑖¹₀⋯¹₀𝑢₀𝑥₁,…,𝑥_𝑑^𝑑_𝑖=1sin𝑛_𝑖𝜋𝑥_𝑖𝑑𝑥₁⋯𝑑𝑥_𝑑+^𝑡₀exp−(𝑡−𝜏)𝜋^2𝑑_𝑖=1𝑛²_𝑖¹₀⋯¹₀𝜒_𝜔(𝜏)𝑥₁,…,𝑥_𝑑^𝑑_𝑖=1sin𝑛_𝑖𝜋𝑥_𝑖𝑑𝑥₁⋯𝑑𝑥_𝑑𝑑𝜏. (2.4)

2.1. The Adjoint Transient Heat Problem

The adjoint transient heat problem has a straightforward definition 𝑃^∗_{𝑣𝑇,𝑔,𝑓}−𝜕_𝑡𝑣−Δ𝑣=𝑓in𝑄,𝑣=𝑣_𝑇inΩ_𝑇,𝑣=𝑔onΣ (2.5) and Dirichlet data compatibility condition, 𝑣_𝑇=𝑔, at the time-space cylinder corner Γ_𝑇. The time reversal operator 𝜅_𝑇∶𝐻^𝑟,𝑠(𝑄)⟶𝐻^𝑟,𝑠(𝑄);𝑣(𝑡,𝑥)⟼𝜅[𝑣](𝑡,𝑥)=𝑣(𝑇−𝑡,𝑥) (2.6) can be used to change the changes of variables 𝑢^∗(𝑡,𝑥)=𝑣(𝑇−𝑡,𝑥) 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: 𝑣𝑡,𝑥₁,…,𝑥_𝑑=^𝑑𝑁_𝑖=1^𝑖_𝑛𝑖=1𝑐_{𝑛1…𝑛𝑑}(𝑇−𝑡)^𝑑_𝑖=1sin𝑛_𝑖𝜋𝑥_𝑖. (2.7)

Definition 2.3 (consistent Cauchy datum). By Consistent Cauchy datum associated with problem (2.2) one means the functions 𝑢₀,𝑔,𝑢_𝑇,𝑔_𝜈∈𝛾₀,𝛾,𝛾_𝑇,𝛾₁𝐻^𝑟,𝑠. (2.8) 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) when 𝑢∈𝐻^{2𝑟+2,𝑟+1}(𝑄) is the solution of problem (2.2) with initial data (𝑢₀,𝑔)=(𝑢|_Ω0,𝑢|_Σ).

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)

3. Inverse Transient Heat Equation Source Problem

The inverse source problem that we address consists in the recovery of the source 𝑓, knowing the extended Dirichlet-to-Neumann map Λ^𝑓_Ω,Σ. When 𝑟=0, the data are regular, Green's function exists, and 𝑓∈𝐿²(𝑄). 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^𝑓_{(𝑢0,𝑔),(𝑢𝑇,𝑔𝜈)} is to find 𝑓∈𝐿²(𝑄) such that IP^𝑓_{𝑢0,𝑢,𝑔𝑇,𝑔𝜈}𝑢_𝑇,𝑔^𝜈=Λ^𝑓_Ω,Σ𝑢₀,𝑔 (3.1) for all given data pair (𝑢₀,𝑔)×(𝑢_𝑇,𝑔^𝜈) 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 𝑃_{𝑢0,𝑔,𝑓} and 𝑃_𝑢0,𝑔,0, one with source 𝑓 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 𝑓∈𝐿²(𝑄) one means the application Λ^𝑓_Ω,Σ−Λ⁰_Ω,Σ∶𝐻¹Ω₀×𝐻^3/2,3/4(Σ)⟶𝐻¹Ω_𝑇×𝐻^1/2,1/4(Σ). (3.2)

Note that the consistence of data (𝑢₀,𝑔) is necessary to the existence of solution for the problems 𝑃_{𝑢0,𝑔,𝑓} and 𝑃_𝑢0,𝑔,0.

Lemma 3.2. Let 𝑢_𝑗,𝑗=1,2,3,… be different solutions of problem (2.2) with the same source 𝑓∈𝐿²(𝑄) and different initial time and Dirichlet data (𝑢_0𝑗,𝑔_𝑗),𝑗=1,2,3,…, respectively. Then,(i)the relative extended Dirichlet-to-Newman operator Λ^𝑓_Ω,Σ−Λ⁰_Ω,Σ is an operator whose functional value depends only on the source function 𝑓∈𝐿²(𝑄) but is independent of the initial time and Dirichlet data (𝑢₀,𝑔); (ii)for all solutions of consistent data problems 𝑃_{𝑓,𝑢0𝑗,𝑔𝑗}, 𝑗=1,2,3,…, with the same source, the source satisfies the systems of integral equations _𝑄𝑓()𝐺(),()𝜏,𝜁𝑇,𝑥,𝜏,𝜁𝜕𝐺𝑡,𝑥,𝜏,𝜁𝜕𝜈_(𝑡,𝑥)𝑑𝜁𝑑𝜏=Λ^𝑓_Ω,Σ−Λ⁰_Ω,Σ𝑢_0𝑗,𝑔_𝑗=Λ^𝑓_Ω,Σ[0,0], (3.3) which depend only on the relative extended Dirichlet-to-Neumann map. Here 𝐺(𝑡,𝑥,𝜏𝜁) is the causal Dirichlet Green function for the transient heat problem;(iii)for all test functions 𝑣 in 𝐻^2,1_{−𝜕𝑡−Δ}(𝑄)=𝑣∈𝐻^2,1(𝑄)∣−𝜕_𝑡𝑣−Δ𝑣=0 (3.4) the source 𝑓(𝑡,𝑥) satisfies the transient heat reciprocity gap equation _𝑄𝑓𝑣𝑑𝑥𝑑𝑡=−_Ω𝑇Λ^𝑓_Ω,•[0,0]𝛾_𝑇[𝑣]𝑑𝑥−_ΣΛ^𝑓_•,Σ[0,0]𝛾[𝑣]𝑑𝜎_(𝑡,𝑥). (3.5)

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 𝑢_𝑇,𝜕_{𝜈(𝑡,𝑥)}𝑢=Λ^𝑓_Ω,Σ𝑢₀,𝑔=Λ^𝑓_Ω,Σ[0,0] (3.6) corresponding to some specific initial Dirichlet data (𝑢₀,𝑔), 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 𝑄∗ containing 𝑄 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 Ω_𝑇 is null, that is, 𝛾_𝑇[𝑣]=0. For these functions we have (iii*)for all test functions 𝑣 in 𝐻^2,1_{−𝜕𝑡−Δ}(𝑄)=𝑣∈𝐻^2,1(𝑄)∣−𝜕_𝑡𝑣−Δ𝑣=0∧𝛾_𝑇[𝑣]=0 (3.7) the source 𝑓(𝑡,𝑥) satisfies the transient heat lateral reciprocity gap equation _𝑄𝑓𝑣𝑑𝑥𝑑𝑡=−_ΣΛ^𝑓_•,Σ[0,0]𝛾[𝑣]𝑑𝜎_(𝑡,𝑥). (3.8)

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 𝑓(𝑡,𝑥)=𝜒_𝜔(𝑡)(𝑥)in𝑄, (4.1) where 𝜔(𝑡), 𝑡∈[0,𝑇], 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,𝑇] into 𝑁 subintervals of length 𝜏>0. Let {𝑡₀,𝑡₁,𝑡₂,…,𝑡_𝑛,𝑡_𝑛+1,…𝑡_𝑁} be the knots of this partition, with 𝑡₀=0 and 𝑡_𝑁=𝑇. For 𝑡_𝑛<𝑡<𝑡_𝑛+1,𝑛=0,1,𝑁−1 we use the 𝜃-scheme approach, [1] for the discretization of (2.2).

By denoting 𝑢^𝑛+1(𝑥) with 𝑥∈Ω, the approximate solution at the time step 𝑡_𝑛+1, the transient system (2.2) is approximated by the following sequence of stationary problems: 𝐻_𝜒^𝑛+1_𝜔−Δ𝑢^𝑛+1+𝜆𝑢^𝑛+1=𝑓_𝑛+𝜒_{𝜔(𝑡𝑛+1)}𝑢inΩ,_𝑛+1=0onΓ,with𝑔^𝜈𝑡_𝑛+1∶=𝜕_𝜈𝑢^𝑛+1onΓ, (4.2) for 𝑛=0,1,2,…,𝑁. Here 𝜆=1/𝜏𝜃 and 𝑓_𝑛=𝑢^𝑛+𝜏(1−𝜃)Δ𝑢^𝑛+𝜏𝜃𝜒_{𝜔(𝑡𝑛)}(𝑥)𝜏𝜃. (4.3) Note that Δ𝑢^𝑛+𝜒_{𝜔(𝑡𝑛)}=𝜕𝑢^𝑛/𝜕𝑡 and that the initial Poisson problem determining the 𝑢₀ and 𝜒_𝜔(0) is 𝐻⁰_{𝑔,𝜒𝜔}−Δ𝑢⁰=𝜒_𝜔(0)(𝑥)𝑢inΩ,₀=0onΓ,with𝑔^𝜈(0)∶=𝜕_𝜈𝑢⁰onΓ. (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 𝜒_{𝜔(𝑡𝑛)}(𝑥), 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 time 𝑡_𝑛 may be further calculated as 𝑓_𝑛=𝜆𝜃𝑢^𝑛−1−𝜃𝜃𝑓_𝑛−1,for𝑛=1,2,… (4.5) with 𝑓₋₁=0 and 𝑓₀=𝜆𝑢⁰. 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, 𝑓_𝑛 and 𝜒_{𝜔(𝑡𝑛+1)}. Let 𝑦^𝑛+1,𝑛=−1,0,1,…, be a solution of𝐻_𝑓^𝑛+1_𝑛−Δ𝑦^𝑛+1+𝜆𝑦^𝑛+1=𝑓_𝑛𝑦inΩ,_𝑛+1=0onΓ,with𝑔^𝜈_𝑦𝑡_𝑛+1∶=𝜕_𝜈𝑦^𝑛+1onΓ, (4.6) and let 𝑤^𝑛+1,𝑛=−1,0,1,…, be solution of 𝐻_𝜒^𝑛+1_𝜔−Δ𝑤^𝑛+1+𝜆𝑤^𝑛+1=𝜒_{𝜔(𝑡𝑛+1)}(𝑥)𝑤inΩ,_𝑛+1=0onΓ,with𝑔^𝜈_𝑤𝑡_𝑛+1∶=𝜕_𝜈𝑤^𝑛+1onΓ. (4.7)

Then, by the superposition principle, the solution of (4.2) is 𝑢^𝑛+1=𝑤^𝑛+1+𝑦^𝑛+1 and the Neumann data will be the sum of the decomposed parts 𝑔^𝜈(𝑡_𝑛+1)=𝑔^𝜈_𝑤(𝑡_𝑛+1)+𝑔^𝜈_𝑦(𝑡_𝑛+1). The Y-problems Equations (4.6) form a discrete sequence of problems with continuous source 𝑓_𝑛 that may be solved before the time increment at 𝑡_𝑛 begins. Its normal derivatives may be calculated 𝑔^𝜈_𝑦𝑡_𝑛+1∶=𝜕_𝜈𝑦^𝑛+1onΓ, (4.8) and subtracted from the synthetic transient Neumann data at knot 𝑡_𝑛+1𝑔^𝜈_𝑤𝑡_𝑛+1=𝑔^𝜈𝑡_𝑛+1−𝑔^𝜈_𝑦𝑡_𝑛+1, (4.9) to produce the data for the modified Helmholtz equations (4.7) that will be used in the reconstruction of the source 𝜒_{𝜔(𝑡𝑛+1)} at time 𝑡^𝑛+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. Let 𝑣 be the space of Helmholtz functions in 𝐻²_𝜆(Ω)=𝑣∈𝐻²(Ω)∶−Δ𝑣+𝜆𝑣=0. (4.10) The reciprocity gap functional, [9] for the Cauchy data in the sequence of Helmholtz problem equations (4.2) is 𝑅^𝜆_{𝑓𝑛+𝜒)𝜔(𝑡𝑛+1}(𝑣)=_Ω𝑓_𝑛𝑣𝑑𝑥+_Ω𝜒_{𝜔(𝑡𝑛+1)}𝑣𝑑𝑥for𝑣∈𝐻²_𝜆(Ω). (4.11) It is a direct consequence of Green's theorem that 𝑅^𝜆_{𝑓𝑛+𝜒)𝜔(𝑡𝑛+1}(𝑣)=−_Γ𝑣𝑔^𝜈𝑡_𝑛+1𝑑𝜎for𝑣∈𝐻²_𝜆(Ω). (4.12) The combination of these equations for test functions in 𝐻_𝜆(Ω) 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 𝑓_𝑛 giving 𝑅^𝜆_{𝜒)𝜔(𝑡𝑛+1}(𝑣)=_Ω𝜒_{𝜔(𝑡𝑛+1)}𝑣𝑑𝑥=−_Γ𝑣𝑔^𝑤,𝜈𝑡_𝑛+1𝑑𝜎for𝑣∈𝐻²_𝜆(Ω). (4.13)

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: 𝑄=𝑡_𝑛+1,𝑡_𝑛+1+𝜏×Ω (4.14) for a small time increment 𝜏 and consider that a field 𝑣∈𝐻^2,1_{−𝜕𝑡−Δ}(𝑄)∧𝛾_{𝑡𝑛+1+𝜏}[𝑣]=0. The equation in the definition of 𝐻^2,1_{−𝜕𝑡−Δ}(𝑄) is then integrated in this time interval: −𝜕_𝑡𝑣−Δ𝑣=0⟹𝑣𝑡_𝑛+1,⋅−Δ^𝑡𝑛+1_𝑡^+𝜏_𝑛+1𝑣(𝑡,⋅)𝑑𝑡=0. (4.15) 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 ^𝑡𝑛+1_𝑡^+𝜏_𝑛+1𝑣(𝑡,⋅)𝑑𝑡=𝜃𝜏𝑣𝑡_𝑛+1,⋅. (4.16) Let us define 𝜆=𝜅²=1/𝜃𝜏. By noting the definition of the space of Helmholtz functions 𝐻²_𝜆(Ω) given by (4.10), we can see that it is a 𝜃 weight average of functions of the space 𝐻^2,1_{−𝜕𝑡−Δ}(𝑄)∧𝛾_{𝑡𝑛+1+𝜏}[𝑣]=0. The same averaging process can be done with the transient reciprocity gap equation (3.5): _Ω^𝑡𝑛+1_𝑡^+𝜏_𝑛+1𝑓𝑣𝑑𝑡𝑑𝑥=−_Γ^𝑡𝑛+1_𝑡^+𝜏_𝑛+1Λ^𝑓_•,Σ[0,0]𝛾[𝑣]𝑑𝑡𝑑𝜎_𝑥. (4.17) With ^𝑡𝑛+1_𝑡^+𝜏_𝑛+1𝑓(𝑡,⋅)𝑣(𝑡,⋅)𝑑𝑡≈𝜃₁𝜏𝑓𝑡_𝑛+1,⋅𝑣𝑡_𝑛+1,⋅_{𝑡𝑛+1𝑡+𝜏𝑛+1},Λ^𝑓_•,Σ[0,0]𝛾[𝑣]𝑑𝑡≈𝜃₂𝜏Λ^𝑓_•,Σ[0,0]𝑡_𝑛+1,⋅𝑣𝑡_𝑛+1,⋅ (4.18) we reproduce an equation that approximates (4.13): _Ω𝑓𝑡_𝑛+1,𝑥𝑣𝑡_𝑛+1𝜃,𝑥𝑑𝑥≈−²𝜃_1ΓΛ^𝑓_•,Σ[0,0]𝑡_𝑛+1,𝑥𝑣𝑡_𝑛+1,𝑥𝑑𝜎_𝑥. (4.19) Note that the data preprocessing for the determination of Λ^𝑓_•,Σ[0,0](𝑡_𝑛+1,.) involves the time homogenization of initial values and is similar to that used for the determination of 𝑔^𝑤,𝜈(𝑡_𝑛+1,⋅).

5. Determining Centroid and Shape

The necessary functions for definition of source centroid, constant function and position function {1,𝑥_𝑖,𝑖=1,…,𝑑}∈𝐻^2,1_{−𝜕𝑡−Δ}(𝑄) 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 𝑥_𝑖=∫_𝑄𝑓𝑥_𝑖𝑑𝑥𝑑𝑡∫_𝑄𝑓1𝑑𝑥𝑑𝑡,𝑖=1,…,𝑑. (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). 𝑥_𝑖=∫_Ω𝑇Λ^𝑓_Ω,•[0,0]𝑥_𝑖𝑑𝑥+∫_ΣΛ^𝑓_•,Σ[0,0]𝑥_𝑖𝑑𝜎_(𝑡,𝑥)∫_Ω𝑇Λ^𝑓_Ω,•[0,0]𝑑𝑥+∫_ΣΛ^𝑓_•,Σ[0,0]𝑑𝜎_(𝑡,𝑥). (5.2)

A problem appears when we work with the modified Helmholtz approximation (4.10), since functions {1,𝑥_𝑖,𝑖=1,…,𝑑}∉𝐻²_𝜆(Ω)={𝑣∈𝐻²(Ω)∶−Δ𝑣+𝜆𝑣=0}. In this case, if we first define 𝜅=(𝜆)=(1/𝜃𝜏), rooted parameter in the auxiliary modified Helmholtz parameter, we have a special set of 𝐻_𝜆(Ω) functions, that is, exp𝜅^𝑑_𝑖=1𝑙_𝑖𝑥_𝑖for𝑙₁,…𝑙_𝑑∈𝑆^𝑑−1. (5.3) These functions are in the domain of the modified Helmholtz equation and form a base for the space 𝐻²_𝜆(Ω). We may construct a dense set by choosing some discrete set of directions 𝑙^𝑗∈𝑆^𝑑−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: sinh𝜅∑^𝑑_𝑖=1𝑙_𝑖𝑥_𝑖−𝑥_𝑖𝜅;cosh𝜅^𝑑_𝑖=1𝑙_𝑖𝑥_𝑖−𝑥_𝑖for𝑙₁,…𝑙_𝑑∈𝑆^𝑑−1 (5.4) 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,(𝑥₁,…,𝑥_𝑑)} are in the space 𝐻²_𝜆(Ω). 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 𝜔⊂Ω∈𝑅^𝑑. By meta centroid 𝑥=(𝑥₁,…,𝑥_𝑑) of this subdomain one means 𝑥_𝑖=∫_𝜔𝑥_𝑖sinh𝜅𝑥_𝑖−𝑥_𝑖/𝜅𝑥_𝑖−𝑥_𝑖𝑑𝑥∫_𝜔sinh𝜅𝑥_𝑖−𝑥_𝑖/𝜅𝑥_𝑖−𝑥_𝑖𝑑𝑥,for𝑖=1,…,𝑑. (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.

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 coordinates 𝑥_𝑖 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 𝑔^𝜈 Neumann data on the boundary is 𝑥^𝜅_𝑖=argmin_Γsinh𝜅𝑥_𝑖−𝑥^𝑐_𝑖𝜅𝑔^𝜈𝑑𝜎(𝑥)²∣𝑥^𝑐∈Ωfor𝑖=1,…,𝑑. (5.6)

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:𝜔=argmin_𝜔cosh𝜅^𝑑_𝑖=1𝑙_𝑖𝑥_𝑖−𝑥^𝜅_𝑖𝑑𝑥−𝐼_Γcosh,𝜅,𝑙,𝑥^𝜅_𝑖,𝑔^𝜈2∶𝜔⊂Ω;𝑙∈𝑆^𝑑−1, (5.7) where the set of directions 𝑙∶=(𝑙₁,…𝑙_𝑑)∈𝑆^𝑑−1 is used to generate linearly independent functionals of the trial shape and 𝐼_Γcosh,𝜅,𝑙,𝑥^𝜅_𝑖,𝑔^𝜈∶=_Γcosh𝜅^𝑑_𝑖=1𝑙_𝑖𝑥_𝑖−𝑥^𝜅_𝑖𝑔^𝜈(𝑥)𝑑𝜎(𝑥) (5.8) 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 cosh𝜅_𝜔𝑥_𝑖−𝑥_𝑖𝑡_𝑛+1𝑑𝑥−_Γcosh𝜅𝑥_𝑖−𝑥_𝑖𝑡_𝑛+1𝑔^𝑤,𝜈𝑡_𝑛+1𝑑𝜎(𝑥)² (5.9) for 𝑖=1,…,𝑑. 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 𝑥_𝑐=.4+.2sin(2𝜋𝑡/𝑇) and size ℎ=.1+.05|sin(2𝜋𝑡/𝑇)|. The total time 𝑇=2𝜋 is incremented with 𝜏=.05, and 𝜃=.7 is chosen for the 𝜃-scheme. The synthetic net flux data 𝑔^{𝜈,(𝑛+1)}=𝑔_𝑤^{𝜈,(𝑛+1)}+𝑔_𝑦^{𝜈,(𝑛+1)} at the two points is 10% noisy, and is substituted in the reciprocity gap functional equatoin (4.13) with the test functions in 𝐻_𝜆(Ω) chosen as in (5.6) 𝑅^𝜆_𝜒𝜔sinh𝜅𝑥_𝑖−𝑥_𝑖𝑡_𝑛+1𝜅=_Γsinh𝜅𝑥_𝑖−𝑥_𝑖𝑡_𝑛+1𝜅𝑔^𝑤,𝜈𝑡_𝑛+1𝑑𝜎(𝑥)=0, (6.1) for 𝑖=1,…,𝑑 in the metacentroid calculations and as 𝑅^𝜅_𝜒𝜔cosh𝜅𝑥_𝑖−𝑥_𝑖𝑡_𝑛+1𝑡_𝑛+1=_𝜔cosh𝜅𝑥−𝑥_𝑖𝑡_𝑛+1𝑑𝑥=_Γcosh𝜅𝑥_𝑖−𝑥_𝑖𝑡_𝑛+1𝑔^𝑤,𝜈𝑡_𝑛+1𝑑𝜎(𝑥),for𝑖=1,…,𝑑 (6.2) for the interval size calculation. Here, as a one-dimensional model, 𝑑=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 𝑥^𝑐 of the centroid is obtained by minimizing the following functional 𝑅^𝜆_𝜒𝜔sinh𝜅𝑥_𝑖−𝑥_𝑖𝑡_𝑛+1𝜅²=_Γsinh𝜅𝑥_𝑖−𝑥_𝑖𝑡_𝑛+1𝜅𝑔^𝑤,𝜈𝑥,𝑡_𝑛+1𝑑𝜎(𝑥)², (6.3) for 𝑖=1,…,𝑑.

Some results have been selected to show the methodology.

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 𝑡=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.

Figure 1 

Heat flux on the transient one-dimensional model.

Figure 2 

Transient source reconstruction for times 0,2, and 4.

Figure 3 

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)^𝑑∈ℝ^𝑑 with a parallelepiped block (voxel) shape. It is supported with a harmonic centroid evolution following the parametric curve (.5(1+.6sin(2𝜋𝑡)),.5(1+.6cos(2𝜋𝑡)),.25+.25𝑡) (6.4) and deforming equally in all directions by the following time rule: ℎ_𝑥=ℎ_𝑦=ℎ_𝑧=ℎ=.151+.25cos𝜋𝑡2, (6.5) where block edge is 2ℎ. The number of harmonics in the Fourier sine series is 20, the Δ𝑥 for spacial collocation is .01, and 𝜃 is chosen as .8. The evolution is calculated for various time steps 𝜏=[.1,.01] in the interval 𝑡∈[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.

Figure 4 

3D metacentroid movement with time increment of .05.

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.

Figure 5 

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

Figure 6 

3D centroid coordinates and block sizes with time increment of .02.

In Figures 7 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 at cosh(𝜅(𝑥_𝑖−𝑥𝑐_𝑖)) and sinh(𝜅(𝑥_𝑖−𝑥𝑐_𝑖)) for 𝑖=1,…,𝑑 is shown in Figure 9.

Figure 9 

Reciprocity gap functional at sinh 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.

Acknowledgments

N. C. Roberty work is partially supported by the Brazilian agencies CNPq and Coppetec Foundation. D. M. de Sousa and M. L. S. Rainha are supported by CNPq. This work is part of the project CAPES/FCT-305/11.