Reconstruction of initial heat distribution via Green function method

1.   Introduction

The thermal diffusion problem has long been studied. Applications can be found in various physical and engineering settings, in particular in hydrology ^[1], material sciences, heat transfer ^[2] and transport problems, etc. The inverse problems for determination of initial temperature, thermal conductivity or heat source from final over determination or other additional measurements have been widely considered by lots of researchers (see, e.g., ^{[3,4,5,6,7,8,9,10,11]}). In ^[11], the authors explored the Tikhonov regularization for simultaneous reconstruction of the initial temperature and heat radiative coefficient in a heat conductive system. They used the final observation and temperature in a small region as the additional data. In ^[9], an exact and analytical representation of the initial heat distribution was given by using only the measurements of temperature and heat flux at one point. In ^[12], the authors proved the uniqueness of the identification of unknown source locations in two-dimensional heat equations from scattered measurements and presented some numerical methods to identify the locations. Hon and Wei introduced the fundamental solution methods in solving inverse heat conduction problems ^[6,7]. In recent years a special kind of inverse heat conduction problem called sideways parabolic equations have also been widely studied ^{[13,14,15,16,17]}. Various regularization methods are presented to solve the severely ill-posed problem and optimal convergence results are obtained in different stopping rules.

Reconstruction of thermophysical properties such as thermal conductivity and heat capacity is another interesting and meaningful research area. Many theoretical and experimental methods have been developed in the literature, they include, among others, the steady-state method, the probe method ^[18], the periodic heating method ^[19], the least-squares method ^[20] and the pulse heating method. For a steady case, the layer potential technique is widely used for analysis of the solution and reconstruction of the thermal conductivity ^[21,22]. In ^[21], the authors construct the Generalized Polarization Tensors based on layer potential method in reconstruction of the shape of the homogeneous conductive body. Besides, layer potential techniques are widely used in wave prorogation problems, including Helmholtz problems, Maxwell problems, Elastic problems and so on (see ^{[23,24,25,26,27,28,29]} and there references there in).

In this paper, we shall consider the reconstruction of the initial heat distribution from final observation in the upper half plane. The Dirichlet boundary condition is given. We firstly use the layer potential method to represent the solution of the forward heat conduction problem in an integral form. Then by exploring the Green function, the Dirichlet boundary problem can be solved analytically without knowing the Neumann boundary condition. By using appropriate extension and Fourier transform, we derive a very simple relation between the initial heat distribution and the final observation. The severely ill-posed problem for reconstruction of initial temperature then is done by introducing a new iteration scheme. Convergence results are investigated under both the a priori and a posteriori stopping rules. The main contribution of this paper is to extend the layer potential techniques to the inverse heat diffusion problem and designing of a new iterative scheme in solving the related inverse heat distribution problem, theoretically. We shall consider the numerical implementation in a forthcoming work.

The organization of this paper is as follows. In the second section, we shall introduce the layer potential techniques together with Green's function in representing the integral solution to the heat diffusion problem, while proofs of the main theorems can be found in Appendix A. In Section 3, the attractive relation between initial heat distribution and the final observation shall be derived by making use of the properties of layer potentials. Base on the relation, a new iterative methods is designed in solving the inverse heat distribution problem. We then analyze some properties of the Fourier transform and derived the convergence results by using the aforementioned iterative methods, in Section 4 and 5, respectively. Some conclusions are made in Section 6.}

2.   Representation of solution

We shall first consider the forward heat conduction problem in the half plane

where $\mathbb{R}_+: = (0, +\infty)$ and $\mathbf{x}: = (x_1, x_2)$. The thermal diffusion coefficient $k$ is supposed to be a positive constant. The initial temperature and boundary temperature are given and denoted by $\mu_0(x_1, x_2)$ and $f(x_1, t)$, respectively. In this paper, we suppose that $\mu_0(\mathbf{x})\in L^2(\mathbb{R}\times \mathbb{R}_+)$ and $f(x, t)\in L^2(\mathbb{R}\times (0, T))$.

We shall consider the analytic solution to the forward problem (2.1) by using the Green function method. Let $\Gamma(\mathbf{x}, t)$ be the fundamental solution to the heat equation, i.e.,

which satisfies

where $\delta_0$ is the Dirac function at the origin ^[31]. To derive the solution to (2.1), we firstly introduce the two well known quantities, single layer potential and double layer potential for heat equation. They are extensions of the classical layer potentials in electrostatic problem ^[32]. For a domain $D$, we denote by $\mathcal{S}_D^k$ (if $k = 1$ then we simply denote by $\mathcal{S}_D$) the single layer potential

for a density function $\varphi(\mathbf{x}, t)\in H^{-1/2, 1}(\partial D, \mathbb{R}_+)$ and denote by $\mathcal{D}_D^k$ (if $k = 1$ then denote by $\mathcal{D}_D$) the double layer potential

for a density funtion $\psi(\mathbf{x}, t)\in H^{1/2, 1}(\partial D, \mathbb{R}_+)$. We shall further define a volume potential by

for a function $\mu(\mathbf{x})\in L^2(D)$. Here $\mathbf{y} = (y_1, y_2)$. By using Green's formula, the solution to (2.1) can be written in the form of layer potentials. We state the result in the following theorem, for the sake of convenience to the reader, where proofs can be found in the appendix.

Theorem 2.1. The solution to the heat conduction problem

has the following form

We see from the solution form (2.6) that it contains both Dirichlet and Neumann boundary conditions. Thus, in order to get the solution of (2.1) by using the integral form (2.6), additional Neumann boundary condition is required. In this paper, we shall not use the Neumann boundary measurement thus we need to avoid the appearance of it in the solution form.

If for a bounded smooth domain $D$, we can introduce the jump formula for layer potentials. For the single layer potential we have on $\partial D$

where the $\pm$ signs mean the limits are taken from outside the domain $D$ and inside the domain $D$, respectively. The integral operator $\mathcal{K}_D^*$ is the adjoint operator of $\mathcal{K}_D$, which is defined by

where p.v. means the Cauchy principle value of the integral. For the jump formula of double layer potential we have

By taking the trace of the solution $u(\mathbf{x}, t)$ and using the jump formulas there holds

for $\mathbf{x}\in \partial D$, where $g: = k \partial u/ \partial \boldsymbol{\nu}$. The above integral equation can be used to solve the Neumann boundary condition $g$. The equation is a Volterra integral equation of the first kind with kernel $\Gamma_k(\mathbf{x}, t)$ and uniqueness of solution has been proved in ^[30].

However, for a unbounded domain $D$ the jump formula for the layer potentials may fail. Besides, it is also not easy and stable to solve the Volterra integral equation of the first kind. In what follows, we seek for the Green function to solve (2.1) without getting involved with the Neumann boundary condition. Denote by $G(\mathbf{x}, t)$ the Green function which satisfies

With the help of Green's function, solution to heat equation can be written by integral formulation.

Theorem 2.2. Let $\mathbf{x}\in D = \mathbb{R}\times \mathbb{R}_+$ then the solution $u(\mathbf{x}, t)$ to (2.1) has the form

The solution form by using Green function method gives us a direct way to compute the solution of (2.1). However, to compute the solution we need first compute the Green's function (2.10), which does not have an explicit form in most cases. However for $D = \mathbb{R}\times \mathbb{R}_+$ we can find the Green function by using the reflection method. We have the following theorem

Theorem 2.3. The solution to (2.1) can be represented by

The proofs of Theorem 2.2 and 2.3 can be found in Appendix. We have shown the solution to the forward problem (2.1) with nonhomogeneous boundary condition. The solution is an integral form associated with the Green function. We shall use this solution form to reconstruct the initial temperature distribution from final overdetermination.

3.   The inverse problem

The inverse problem we consider in this paper is the reconstruction of the initial temperature distribution $\mu_0(\mathbf{x})$ and generally the solution $u(\mathbf{x}, t)$ from the following problem

We suppose $f(x_1, t)\in L^2(\mathbb{R}\times(0, T))$ and $\mu_T(\mathbf{x})\in L^2(\mathbb{R}\times \mathbb{R}_+)$. For convenience to use the Fourier transform we also suppose those functions are in $L^1$ for any fixed time $t$. Define $\|\cdot\|$ as the $L^2$ norm with respect to $\mathbf{x}$ in $\mathbb{R}\times \mathbb{R}_+$. Denote by $\mathcal{G}_D$ the operator on $L^2(\mathbb{R}\times \mathbb{R}_+)$

where

By setting $t = T$ in the solution form in Theorem 2.3 we get

where

The inverse problem here is to solve the integral equation (3.3) to get the initial temperature distribution $\mu_0(\mathbf{x})$. We see that $\mathcal{G}_D$ is an integral operator with the kernel $G(\mathbf{x}, kT)$. The following properties on the operator $\mathcal{G}_D$ is straight forward

Lemma 3.1. The operator $\mathcal{G}_D$ in (3.3) is a self-adjoint and bounded operator on $L^2(\mathbb{R}\times \mathbb{R}_+)$.

Proof. Since the kernel of the operator $\mathcal{G}_D$ is symmetric we immediately get that $\mathcal{G}_D$ is self-adjoint on $L^2(\mathbb{R}\times \mathbb{R}_+)$. We shall show the bound of $\mathcal{G}_D$ on $L^2(\mathbb{R}\times \mathbb{R}_+)$. Define

Then there holds

where $\mathbf{x} = (x_1, x_2)$. It is easy to see that the right hand side above is a convolution of functions $\frac{1}{4k\pi T}e^{-\frac{| \mathbf{x}|^2}{4kT}}$ and ${\mu}(\mathbf{x})$. Thus by Young's inequality we have

which completes the proof.

We shall explore a method to solve (3.3). To proceed, we present some preliminary results

Lemma 3.2. There holds the following identities

and

Proof. For $n = 0$ the result is obvious. We assume $n\geq 1$. Define

then by integral by parts we have

and thus by recursion we obtain

Next by Taylor expansion of $\cos(m\pi y)$ we have

which completes the proof.

Lemma 3.3. There holds the following identity

where $i$ is the imaginary unit.

Proof. By using (3.6) and some elementary calculations we have

With (3.7) on hand, we can get the following result

Theorem 3.1. There holds the following relation

Proof. By using the definition of $\mathcal{G}_D$ in (3.2) we have

By using change of variables we have

thus by using (3.7) we finally get

which completes the proof.

The relation (3.8) inspires us that we can actually use the Fourier transform to solve (3.3). To do this, we shall actually extend the functions $\mu_0(\mathbf{x})$ and $h(\mathbf{x})$ to $\mathbb{R}^2$ by setting

Under this kind of extension we have

and

Define the Fourier transform operator $\mathcal{F}$ by

and the corresponding inverse Fourier transform $\mathcal{F}^{-1}$ by

where $\boldsymbol{\xi} = (\xi, \eta)$. Then by taking the Fourier transform on both sides of (3.10) and using (3.7) we obtain

Thus we have

and so

We mention that to get $\hat{\mu}_0(\boldsymbol{\xi})$ directly from (3.13) is quite unstable especially when the data $h(\mathbf{x})$ is not exactly given. We shall use the iteration schemes to solve this kind of problem when the final over determination $\mu_T(\mathbf{x})$ and the boundary measurement $f(x, t)$ are not exactly given. For the sake of simplicity we denote by $h^{\epsilon}(\mathbf{x})$ the perturbation of the function $h(\mathbf{x})$ in (3.3) and satisfies

where $\epsilon$ is a small number which can be treated as the noise level of the measurement. Similar to ^[13,14], we introduce the following iteration scheme to solve (3.13)

where

and ${{\chi }_{\varsigma }}$ is the characteristic function of $\varsigma = \{\boldsymbol{\xi}|\ k|\boldsymbol{\xi}|^2\le \rho \}$, that is

The parameter $\rho$ here is a selective positive number which can be chosen in a flexible way depending on the stopping rule that is used for terminating the iteration scheme. The parameter $n$ is a selective nature number which is used to control the iteration speed. Generally speaking, the lager $n$ is the less iteration steps are required to get an approximate solution. The iterative steps actually decrease exponentially with the parameter $n$ increase normally. However, the error estimation for the solution behaviors badly as $n$ increases as well.

4.   Fourier transform of $h$

The iteration scheme (3.15) requires us to get the Fourier transform of the perturbation $h^{\epsilon}(\mathbf{x})$. In fact, the perturbation comes from the noise of the final observation $\mu_T$ and the boundary condition $f$ in (3.4). In this section we consider the Fourier transform of $h$. To extend $h$ to $\mathbb{R}^2$ we only need to extend the final observation $\mu_T(\mathbf{x})$ by

We observe that (3.4) can be written as follows

Taking the Fourier transform on both sides of (4.1) and using (3.7) we get

where by $\hat{f}(\xi, s)$ we mean the one dimensional Fourier transform of $f(x_1, s)$ on $x_1$, i.e.,

and so the inverse one dimensional Fourier transform is

5.   Error estimates

In this section, we consider the convergence results of iteration scheme (3.15) by using both the a priori and the a posteriori stopping rules. Before presenting the error estimate of the initial temperature distribution $\mu_0(\mathbf{x})$, we suppose the a priori bound on $\mu_0(\mathbf{x})$

where

We see that if $p = 0$ then the norm $\|\cdot\|_p$ turns to $L^2$ norm.

5.1.   The a priori stopping rule

By using the prior information (5.1) we can design the following stopping rule

where $\lfloor a\rfloor$ denotes the largest integer that is not greater than $a$. We note that in the stopping rule (5.2), if $n$ increases then the iteration steps $j$ are greatly reduced.

Theorem 5.1. Let $\mu_0(\mathbf{x})$ be the initial temperature distribution in (3.1).Suppose ${{\mu }_T^{\epsilon }}(\mathbf{x})$ is the measured temperature at $t = T$ and $f^{\epsilon}(x_1, t)$ is themeasured boundary condition. Suppose the corresponding function $h^{\epsilon}(\mathbf{x})$ by (3.4) satisfies(3.14). Let $\mu_{j}^{\epsilon }(\mathbf{x})$ be the iterates defined by (3.15) withinitial iterate $\mu_{0}^{\epsilon }(\mathbf{x}) = 0$.If (5.1) is satisfied for $M > 0$, $p > 0$ and we select $j\sim\left\lfloor {{(M/\epsilon)}^{1/n}} \right\rfloor$ as the prior stopping rule and

then we have the following error estimate

where $C$ is a constant depending on $j$ and $n$.

Proof. Set $\alpha = \sqrt[n]{v(\boldsymbol{\xi})}$ then $0 < \alpha \le 1$. By (3.15) we obtain

Take ${{p}_{j}}(\alpha) = \sum\limits_{i = 0}^{j-1}{{{(1-\alpha)}^{i}}}$, ${{r}_{j}}(\alpha) = 1-\alpha {{p}_{j}}(\alpha) = {{(1-\alpha)}^{j}}$, we have the following elementary results (cf.^[33]):

where

From (3.13) we have

then there holds the estimates

Next we derive the estimates for ${{I}_{1}}$ and ${{I}_{2}}$, respectively.

By using the definition of $\varsigma$ and the choice of $\rho$ we thus have

With similar strategy and (5.2) we estimate $I_2$

where (5.5) and the choice of $\rho$ are taken into account and the selection of constant $C$ is only dependent on $j$ and $n$. We come to the conclusion by combining the estimates of $I_1$ and $I_2$ and using Parseval equality.

5.2.   The a posteriori stopping rule

The prior stopping rule needs the a priori bound of $\mu_0$ to stop the iteration process. However if such bound can not be obtained accurately, the iteration may not stop at a "good time" and thus the result may not be satisfactory. Thus an a posteriori stopping rule is required. We introduce the widely-used "discrepancy principle" due to Morozov ^[34]. We set the discrepancy principle in the following form

where $k_*$ is the first iteration step which satisfies the left inequality of (5.6).

Theorem 5.2. Let $\mu_0(\mathbf{x})$ be the initial temperature distribution in (3.1).Suppose ${{\mu }_T^{\epsilon }}(\mathbf{x})$ is the measured temperature at $t = T$ and $f^{\epsilon}(x, t)$ is themeasured boundary condition. Suppose the corresponding function $h^{\epsilon}(\mathbf{x})$ by (3.4) satisfies(3.14). Let $\mu_{j}^{\epsilon }(\mathbf{x})$ be the iterates defined by (3.15) withinitial iterate $\mu_{0}^{\epsilon }(\mathbf{x}) = 0$.If (5.1) is satisfied for $M > 0$, $p > 0$ and we select the discrepancy principle (5.6) as the stopping rule and

then we have the following error estimate

where $C$ is a constant depending on $j$ and $n$.

Proof. It sufficient to prove ${{j}_{*}} \sim {{(M/\epsilon)}^{1/n}}$ by Theorem 5.1. We have

furthermore

thus we obtain

Set $\tau > 3$ and according to (5.6) we have $(\tau -3)\epsilon < {{n}^{n}}{{j}^{-n}}M$ and thus $j\sim\left\lfloor {{(\frac{M}{\epsilon })}^{1/n}} \right\rfloor$.

We remark that although the parameter $\rho$ appears in (3.15) requires the a priori bound on $\mu_0$, we do not care much on the choice. In fact the choice of $\rho$ is just for technique setup. We can actually set a very large number $M$ such that it is lager than the sharp upper bound on $\mu_0$. By using the discrepancy principle (5.6) as the stopping rule, the convergence rate and reconstruction performance are still ensured.

6.   Conclusions

We have introduced the layer potential technique to get the analytic solution of the heat diffusion problem in the half plane. The Green function has been constructed to avoid the usage of the Neumann boundary condition. We have combined the Fourier transform and the Green function method to get the relation between the initial temperature distribution and the given final observation. We designed an iteration scheme to reconstruct the initial temperature distribution. Although the whole analysis is on two dimensional space heat conduction problem, it can be similar used to analyze higher dimensional problem. The method can also be used to reconstruct the initial heat distribution on any kind of domain and numerical realization of two and three dimensional inverse heat conduction problem will be a forthcoming work.

Acknowledgements

The work of X. Fang was supported by NSF grant China No. 72001077, Humanities and Social Sciences Foundation of the Ministry of Education No. 20YJC910005, NSF grant of Hunan No. 2021JJ30192. The work of Y. Deng was supported by NSF grant of China No. 11971487 and NSF grant of Hunan No. 2020JJ2038. The work of Z. Zhang was supported by Hunan Provincial Natural Science Foundation of China Grant No. 2021JJ30297 and Research and Innovation team of Hunan Institute of Science and Technology (Grant No. 2019-TD-15).

A.   Proofs of Theorem 2.1 to 2.3

In this section, we show the main proofs to some Theorems for the sake of convenience to the reader.

Proof of Theorem 2.1 By using the Green's theorem, we have

where we have changed the order of integration and used the integral by parts. Thus there holds

which completes the proof.

Proof of Theorem 2.2 We first show the relation between $G(\mathbf{x}, t)$ and $\Gamma(\mathbf{x}, t)$. If fact, let $v(\mathbf{x}, t)$ be the solution to

then it is easy to see that

By using the Green's formula there holds for $u(\mathbf{x}, t)$ and $v(\mathbf{x}, t)$

By subtracting the above equation from (2.6) and using the fact that

we immediately get the conclusion.

Proof of Theorem 2.3 Let $\mathbf{x} = (x_1, x_2)$, $\mathbf{y} = (y_1, y_2)$. By using the reflection method, it is easy to find the Green function of (2.10) for $D = \mathbb{R}\times \mathbb{R}_+$

To use the result in Theorem 2.2, we compute

By substituting the Green function and the above formula into the solution form in Theorem 2.2 we come to the conclusion.

Conflict of interest

The authors declare there is no conflict of interest.