Finding the Time-dependent Term in 2D Heat Equation from Nonlocal Integral Conditions

The aim of this paper is to find the time-dependent term numerically in a two-dimensional heat equation using initial and Neumann boundary conditions and nonlocal integrals as over-determination conditions. This is a very interesting and challenging nonlinear inverse coefficient problem with important applications in various fields ranging from radioactive decay, melting or cooling processes, electronic chips, acoustics and geophysics to medicine. Unique solvability theorems of these inverse problems are supplied. However, since the problems are still ill-posed (a small modification in the input data can lead to bigger impact on the ultimate result in the output solution) the solution needs to be regularized. Therefore, in order to obtain a stable solution, a regularized objective function is minimized in order to retrieve the unknown coefficient. The two-dimensional inverse problem is discretized using the forward time central space (FTCS) finite-difference method (FDM), which is conditionally stable and recast as a nonlinear least-squares minimization of the Tikhonov regularization function. Numerically, this is effectively solved using the MATLAB subroutine lsqnonlin. Both exact and noisy data are inverted. Numerical results for a few benchmark test examples are presented, discussed and assessed with respect to the FTCS-FDM mesh size discretisation, the level of noise with which the input data is contaminated, and the choice of the regularization parameter is discussed based on the trial and error technique.


Introduction
The identification of coefficients in inverse heat conduction problems for the parabolic heat equation continues to receive significant attention in a variety of fields, such as heat transfer, oil recovery, groundwater flow, and finance. Some researchers investigated the case of simultaneous identification of coefficients in two-dimensional heat conduction problems, see Refs. [1][2][3][4][5][6][7][8] to mention only a few. Inverse identification problems with integral measurements arise naturally in various physics and engineering models, such as radioactive nuclear decay [9], reactive transport in fluid flows in porous media [10] and semiconductor devices [11].
The determination of physical properties such as thermal conductivity using measured temperature, nonlocal integral or heat flux values at wall sites is an important inverse problem. A common determination strategy is the indirect one where one can minimize the gap between a computed solution and the measured data (observations) via an iterative process [12]. The main difficulty in this kind of problem is that there are usually so few observations that one finds hard to evaluate the spatial derivative of temperature by simple numerical differentiation. Therefore, heavier and more time-consuming optimization techniques are needed to obtain reliable results.
The one-dimensional inverse problems involving the identification of the time-dependent thermal conductivity/diffusivity coefficients of the heat equation with nonlocal and integral over-determination conditions in different statements were studied in Refs. [13][14][15][16][17][18]. Additionally, Abbas et al. [19] employed a finite difference approach for approximate solution of one-dimensional wave equation while Abbas et al. [20] applied a cubic B-spline collocation scheme for the solution of a reaction-diffusion, with initial and Neumann boundary conditions. A two-time level implicit technique is proposed for the approximate solution of the nonclassical diffusion problem with nonlocal boundary condition in the study [21]. Nazir et al. [22], developed various numerical solution techniques of the advection-diffusion equation.
The estimation of thermal properties for the multi-dimensional is rather scarce in the literature [23,24]. The aim of this paper is to consider a two-dimensional coefficient identification problem to estimate the timedependent thermal conductivity component with initial and Neumann boundary conditions from non-local integral conditions. The inverse problems presented in this paper have already been showed to be locally uniquely solvable by Koval'chuk [25], but no numerical identification has been tried so far, therefore, the main goal of this work is to attempt numerical realization of this problem. Moreover, the novelty consists in the development of a convergent numerical optimization method for solving this nonlinear inverse coefficient problem for the heat equation. Numerically, the implementation is realised using the MATLAB subroutine lsqnonlin.
The rest of the paper is organized as follows. Section 2 describes the mathematical formulation of the inverse problems. The numerical forward time central space FDM discretization of the direct problem is described in Section 3. Section 4 introduces the regularized nonlinear minimization used for solving the inverse problems under investigation. In Section 5, numerical results and discussion are illustrated. Finally, conclusions are presented in Section 6.

Mathematical Formulation of the Inverse Problems
In the domain D ¼ Â 0; T ð Þ, where is the rectangle 0; l ð Þ Â 0; h ð Þ, we consider the inverse problem of finding the time-dependent conductivity a t ð Þ > 0 in the two-dimensional heat equation where u ¼ u x; y; t ð Þ is a unknown temperature, subject to the initial condition the Neumann boundary conditions u x 0; y; t ð Þ¼ l 1 y; t ð Þ; u x l; y; t ð Þ ¼ l 2 y; t ð Þ; y 2 0; h ; t 2 ½0; T ½ ; u y x; 0; t ð Þ¼v 1 x; t ð Þ; u y x; h; t ð Þ¼v 2 x; t ð Þ; x 2 0; l ; t 2 ½0; T ½ ; and the nonlocal integral condition where ', l 1 , l 2 , v 1 , v 2 , j are the given functions, while the functions u and a are unknown.
The uniqueness of solution and continuous dependence on the input data for the solution pair a; u ð Þ to this problem and other related inverse problem have been stated in the next two subsections.

Theorem 1 (Uniqueness of solution of IP1)
Assume that the following conditions are satisfied: where Then, the IP1 Eqs. (1) to (5) has a unique solution a t ð Þ; u x;

Inverse Problem 2 (IP2)
A related inverse problem (termed IP2) which requires the identification of the time-dependent coefficient a t ð Þ > 0 and the temperature u x; y; t ð Þ, which satisfy the two-dimensional heat equation in Eq. (1), the conditions Eqs. (2) to (4) and nonlocal integral over-specification condition The IP2 was previously investigated in Koval'chuk [25], where its unique solvability has been proved and given as follows.

Discretization of the Forward Problem
Consider the forward (direct) initial boundary value problem given by Eqs.
(1) to (4), when a t ð Þ, ' x; y ð Þ; We apply the FTCS-FDM to solve the Eq. (1) which is conditionally stable, [26]. So we obtain u nþ1 i;j À u n i;j Dt ¼ a n u n iÀ1;j À 2u n i;j þ u n iþ1;j ðDxÞ 2 þ u n i;jÀ1 À 2u n i;j þ u n i;jþ1 for i ¼ 1; M 1 À 1, j ¼ 1; M 2 À 1 and n ¼ 0; N . In order to obtain explicit expression for u nþ1 i;j , the Eq. (7) is rearranged as u nþ1 i;j ¼ u n i;j þ a n Dt ðDxÞ 2 u n iÀ1;j À 2u n i;j þ u n iþ1;j þ a n Dt ðDyÞ 2 u n i;jÀ1 À 2u n i;j þ u n i;jþ1 : The initial condition Eq. (2) gives the Neumann boundary conditions Eqs. (3) and (4) give where u n À1;j , u n M 1 þ1;j , u n i;À1 and u n i;M 2 þ1 are fictitious values situated outside the domain. These values can be obtained as follows: The stability condition of the explicit FDM Eq. (8) is given as [26] aDt whereã is a maximum value of a t ð Þ.
Finally, the trapezoidal rule is applied for all integrals conditions in Eqs. (5) and (6) as where and

Numerical Solution of the Inverse Problems
The numerical solution of the inverse problems Eqs. (1) to (5) or, Eqs. (1) to (4) and (6) is obtanied by minimizing the nonlinear Tikhonov regularization function for IP1, and for IP2, where u solves the direct problem Eqs. (1) to (4) for given a t ð Þ, and > 0 is regularization parameter to be prescribed. The unregularized case, i.e., ¼ 0, yields the ordinary nonlinear least-squares method which is usually producing unstable solutions when noisy data are inverted. The minimization of F 1 , or F 2 , is performed using the MATLAB subroutine lsqnonlin [27]. This subroutine attempts to find the minimum of a sum of squares by starting from a given initial guess. Simple bounds on the variable are allowed and the explicit calculation (analytical or numerical) of the gradient is not required to be supplied by the user.

Numerical Results and Discussion
This section presents two benchmark test examples in order to test the accuracy and stability of the FDM numerical procedure introduced in Section 3. The root mean square errors (RMSE) is used to evaluate the accuracy of the numerical results as follows: Let us take l ¼ h ¼ T ¼ 1, for simplicity. The lower and upper bounds for the time-dependent conductivity coefficient a t ð Þ > 0 is taken as 10 À3 and 10 3 , respectively.
The inverse problems are solved subject to both exact and noisy input data which is numerically simulated as follows: j e 1 t n ð Þ ¼ j t n ð Þ þ e 1 n ; n ¼ 0; N; j e 2 t n ð Þ ¼j t j À Á þ e 2 n ; n ¼ 0; N ; where e 1 n , e 2 n are random variables generated from a Gaussian normal distribution with mean zero and standard deviations r 1 and r 2 , respectively, given by where p represents the percentage of noise. In the case of noisy data Eqs. (18) and (19), we replace j t n ð Þ by j e 1 t n ð Þ for n ¼ 0; N in Eq. (15) , andj t j À Á byj e 2 t n ð Þ for n ¼ 0; N in Eq. (16).

Example 1 (for IP1)
Consider the IP1 Eqs. (1) to (5) with unknown time-dependent coefficient a t ð Þ, and solve it with the input data: We remark that the conditions A1 ð Þ-A5 ð Þ of Theorem 1 are holds and thus, the uniqueness of the solution is ensured. It can be easily verified that the analytical solution u x; y; t ð Þ; a t ð Þ ð Þis given by Let us solve the direct problem Eqs. (1)-(4) with the data Eq. (22), when a t ð Þ is known and given by Eq. (25), using the FDM with the mesh sizes Dx ¼ Dy ¼ 0:05 and Dt ¼ 0:0125 ensure that the stability condition Eq. (12) is always satisfied. The analytical Eq. (24) and numerical solutions for u x; y; t ð Þ is plotted in Fig. 1(a). The exact Eq. (23) and numerical solutions for j t ð Þ is shown in Fig. 1(b). From this figure, one can observe that an excellent agreement is obtained. (1) to (5) with the input Eqs. (22) and (23) using the subroutine lsqnonlin minimization of the objective functional F 1 in Eq. (15) with the initial guess for the vector a ¼ ða t n ð ÞÞ n¼1;N given by The mesh size is taken as the direct problem above and we start the examination for reconstructing the time-dependent coefficient a t ð Þ, when p ¼ 0 in Eq. (20). The objective function F 1 is illustrated, as a function of the number of iterations, in Fig. 2(a), where a convergence, which is monotonically decreasing, is obtained in about 31 iterations to reach a very low prescribed tolerance of O 10 À29 À Á . Fig. 2(b) illustrates the corresponding numerical results for the time-dependent term a t ð Þ; obtaining RMSE(a)=1.9E-04. From this figure, we found an excellent agreement, between the analytical and the numerical solutions. Now, we study the stability of the approximate solution with respect to various levels of p 2 0:01%; 0:1% f g noise in Eq. (20) included in the input data j t ð Þ. The decreasing monotonic convergence of the objective function F 1 , without and with regularization is shown in Fig. 3. The reconstruction of the estimated a t ð Þ is shown in Figs. 4 and 5. From Figs. 4(a) and 5(a) one can notice that unstable (highly oscillation) and inaccurate solutions for a t ð Þ are obtained with RMSE a ð Þ 2 0:0100; 0:0272 f g , if no regularization, i.e. ¼ 0, is employed. This is expected since the IP1 is ill-posed. Therefore, regularization is required in order to retrieve the loss of stability. From all regularization parameters that were selected, we conclude that ¼ 10 À4 , for p ¼ 0:01% and ¼ 10 À3 , for p ¼ 0:1% noise gives a stable and reasonable accurate numerical solution for the time-depndent term a t ð Þ; obtaining RMSE a ð Þ 2 f2:2E-04; 0:0014g, respectively. From Figs. 2(b), 4, 5 and Tab. 1 it is observed that when p decreasing from 0:1% to 0:01% and then to zero, the numerically achieved results become more accurate and stable.

Example 2 (for IP2)
In this example, we consider the IP2 given by Eqs. (1) to (4) and (6) with unknown time-dependent coefficient a t ð Þ, and solve it with the input data: RMSE(a) 0.0272 0.0049 0.0021 0.0014 0.0090 The data j t ð Þ given by Eq. (5) is replaced by the measured dataj t ð Þ given by Eq. (6) as We observe that the conditions A1 ð Þ, A2 ð Þ, B1 ð Þ-B3 ð Þ of Theorems 1 and 2 are fulfilled and thus, the uniqueness condition of the solution is guaranteed. In fact, the exact solution u x; y; t ð Þ; a t ð Þ ð Þof the IP2 is given by The initial guess for the vector a is taken as 0.05, namely, We take a mesh size with M 1 ¼ M 2 ¼ 10 and N ¼ 60, which together with the upper bound 10 3 for the time-dependent term a satisfying the stability condition Eq. (12).
As in Example 1, first consider the case where there is no noise (p ¼ 0) in the input dataj t ð Þ in Eq. (21). The objective function F 2 Eq. (16), with ¼ 0 is plotted in Fig. 6(a), where a convergence, which is monotonically decreasing, is obtained in about 31 iterations to achieve a low tolerance of O 10 À29 À Á . The analytical Eq. (30) and numerical solutions for a t ð Þ is illustrated in Fig. 6(b), where the reconstructed timewise term is in very excellent agreement with their corresponding exact solutions, obtaining with RMSE a ð Þ ¼2.9E-3.
Next, we add p 2 0:01%; 0:1% f gnoise in the input dataj t ð Þ, as in Eq. (21). The corresponding exact Eq. (30) and numerical solutions for a t ð Þ are illustrated in Figs. 7 and 9 for various regularizations. When ¼ 0, we achieve unstable and inaccurate numerical solutions, with RMSE a ð Þ ¼ 0:0281 for p ¼ 0:01%, see Fig. 7(a), and with RMSE a ð Þ ¼ 0:1590 for p ¼ 0:1%, see Fig. 9(a). We apply the Tikhonov regularization method in order to overcome this instability. We deduce that ¼ 10 À5 , for p ¼ 0:01%, see Fig. 7(b), and ¼ 10 À4 , for p ¼ 0:1%, see Fig. 9(b), provides a accurate and stable approximate solutions for a t ð Þ, with RMSE a ð Þ ¼0.0038 and RMSE a ð Þ ¼0.0116, respectively. Also, from Figs. 7, 9 and Tab. 2 it is observed that when p decreasing from 0:1% to 0:01% and then to zero, the numerically achieved results become more accurate and stable. The related absolute errors between the approximate and analytical Eq. (29) solutions for u x; t ð Þ with p 2 0:01%; 0:1% f g noise, without and with regularization parameters, are illustrated in Figs. 8 and 10. From these figures it can be noticed that the temperatures u x; t ð Þ component is stable and accurate by adding a penalty term > 0 as in Eq. (16) to stabilize the solution. The same conclusions as those obtained for Example 1 can be carried out about the stable reconstructions for the time-dependent coefficient a t ð Þ by observing these figures.

Conclusions
A couple of inverse problems which require finding a time-dependent thermal conductivity coefficient a t ð Þ along with the temperature u x; y; t ð Þ in the two-dimensional heat equation Eq. (1) with initial and Neumann boundary conditions from the nonlocal integrals as over-determination conditions Eqs. (5) and (6) have been investigated numerically. A direct solver based on the forward time central space FDM has been developed. The inverse solver based on a nonlinear least-squares minimization has been solved computationally using the MATLAB subroutine lsqnonlin. The Tikhonov regularization has been used in order to obtain stable solutions since the inverse problem is ill-posed and very sensitive to noise. The main difficulty in regularization when we solve the ill-posed problem is how to choose an appropriate regularization parameter which compromises between accuracy and stability. However, one can use techniques such as the L-curve method [28] or Morozov's discrepancy principle [29] to find such a  Figure 10: The absolute error between the analytical Eq. (29) and numerical solutions for u x; y; t ð Þwith: (a) ¼ 0, (b) ¼ 10 À5 , (c) ¼ 10 À4 and (d) ¼ 10 À3 , with p ¼ 0:1% noise, for Example 2 parameter, but in our work we have used trial and error. The numerical results have been presented and discussed for the two inverse problems, showing that accurate and stable approximate solutions have been achieved.