# Inverse Estimates for Nonhomogeneous Backward Heat Problems

Authors: Tao Min; Weimin Fu; Qiang Huang
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

We investigate the inverse problem in the nonhomogeneous heat equation involving the recovery of the initial temperature from measurements of the final temperature. This problem is known as the backward heat problem and is severely ill-posed. We show that this problem can be converted into the first Fredholm integral equation, and an algorithm of inversion is given using Tikhonov's regularization method. The genetic algorithm for obtaining the regularization parameter is presented. We also present numerical computations that verify the accuracy of our approximation.


---

## Body

## 1. Introduction

Inverse heat conduction problems (IHCPs) appear in many important scientific and technological fields [1–6]. Hence analysis, design implementation, and testing of inverse algorithms are also great scientific and technological interest. In this paper we will examine some inversion techniques in order to estimate the initial temperature distribution of an inverse nonhomogeneous heat conduction problem; this problem is called the backward heat problem, the backward Cauchy problem, or the final value problem. As is known, the nonhomogeneous problem is severely ill-posed; that is, solutions do not always exist, and in the case of existence, these do not depend continuously on the given data. In fact, from small noise contaminated physical measurements, the corresponding solutions have large errors. It makes it difficult to do numerical calculations. Hence, a regularization is in order. Lattes and Lions, in [7], regularized the problem by adding a “corrector” to the main equation. Gajewski and Zacharias [8] considered a similar problem. Nakamura et al. [9] used transformation techniques to solve the initial inverse problem in heat conduction and Al-Khalidy [10] dealt with the problem numerically. For a comprehensive review of the literature and summary of various approaches in the field of inverse heat conduction problems, one can consult the books by Beck et al. [11] and by Hensel [12]. Although there are many papers on the linear homogeneous case of the backward problem, we only find a few papers on the nonhomogeneous case, such as in [13, 14].The plan of this paper is as follows. In Section2, we formulate a one-dimensional IHCP. In Section 3, we show that this problem can be converted into the first Fredholm integral equation, the numerical algorithm is derived, and an algorithm of inversion is also given using Tikhonov's regularization method. We use the genetic algorithm for the choice of the regularization parameter in Section 4. Section 5 presents numerical experiments to show the efficiency of the proposed method. Section 6 ends this paper with a brief discussion on some numerical aspects.


## 2. Mathematical Model

### 2.1. The Direct Problem

The direct (forward) problem consists of a transient heat conduction problem in a slab with adiabatic boundary condition and initially at a temperature denoted byf(x).The mathematical formulation of this problem is given by the following nonhomogeneous heat equation:(1)∂u∂t-k∂2u∂x2=φ(x,t),0<x<L,t>0,(2)u(0,t)=0,t>0,u(L,t)=0,t>0,u(x,0)=f(x),0≤x≤L,
where u(x,t) (temperature), φ(x,t) (source term), f(x) (initial condition), x (spatial variable), and t (time variable) are dimensionless quantities and k denotes the dispersion coefficient.For the direct problem where the initial conditionf(x) is specified, the problem given by (1) and (2) is concerned with the determination of the temperature distribution u(x,t) in the interior region of the solid as a function of time and position.


### 2.2. Inverse Problem

For the inverse problem, the initial conditionf(x) is regarded as being unknown. In addition, an overspecified condition is also considered available. To estimate the unknown coefficient f(x), the additional information,
(3)u(x,T)=g(x),
is given at time T, over a specified space interval 0≤x≤L. We note that the measured overspecified condition u(x,T)=g(x) should contain measurement errors. Therefore the inverse problem can be stated as follows: estimate the unknown function f(x) by utilizing the above-mentioned measured data.


## 2.1. The Direct Problem

The direct (forward) problem consists of a transient heat conduction problem in a slab with adiabatic boundary condition and initially at a temperature denoted byf(x).The mathematical formulation of this problem is given by the following nonhomogeneous heat equation:(1)∂u∂t-k∂2u∂x2=φ(x,t),0<x<L,t>0,(2)u(0,t)=0,t>0,u(L,t)=0,t>0,u(x,0)=f(x),0≤x≤L,
where u(x,t) (temperature), φ(x,t) (source term), f(x) (initial condition), x (spatial variable), and t (time variable) are dimensionless quantities and k denotes the dispersion coefficient.For the direct problem where the initial conditionf(x) is specified, the problem given by (1) and (2) is concerned with the determination of the temperature distribution u(x,t) in the interior region of the solid as a function of time and position.


## 2.2. Inverse Problem

For the inverse problem, the initial conditionf(x) is regarded as being unknown. In addition, an overspecified condition is also considered available. To estimate the unknown coefficient f(x), the additional information,
(3)u(x,T)=g(x),
is given at time T, over a specified space interval 0≤x≤L. We note that the measured overspecified condition u(x,T)=g(x) should contain measurement errors. Therefore the inverse problem can be stated as follows: estimate the unknown function f(x) by utilizing the above-mentioned measured data.


## 3. Algorithm Analysis

The solution of the direct problem for a given initial conditionf(x) is explicitly obtained using separation of variables for 0<x<L,t≥0. Consider the initial boundary problems (1) and (2); we look for solutions u in the form u(x,t)=T(t)X(x). We look at the following eigenvalue problem:
(4)X′′+λX=0X(0)=X(L)=0.
The eigenvalues and the corresponding eigenfunctions are
(5)λ=(nπL)2,Xn(x)=sinnπLx.
Now we set
(6)u(x,t)=∑n=1∞Tn(t)sinnπLx.
Formally computing ut and uxx and substituting into (1), we get
(7)φ(x,t)=ut-kuxx=∑n=1∞[Tn′(t)+k(nπL)2Tn(t)]sin(nπxL).
Hence expanding f and ϕ into the following Fourier series:
(8)φ(x,t)=∑n=1∞φn(t)sin(nπxL),f(x)=∑n=1∞fnsin(nπxL),
where
(9)φn(t)=2L∫0Lφ(x,t)sin(nπxL)dx,fn=2L∫0Lf(x)sin(nπxL)dx,
we get
(10)∑n=1∞[Tn′(t)+k(nπL)2Tn(t)]sin(nπxL)=∑n=1∞φn(t)sin(nπxL).
The uniqueness of the Fourier expansion leads to the family of ODE's:(11)Tn′(t)+k(nπL)2Tn(t)=φn(t).
In addition,
(12)u(x,0)=∑n=1∞Tn(0)sinnπLx=∑n=1∞fnsin(nπxL)=f(x),
so that
(13)Tn(0)=fn,n≥1.
Solving (11)–(13), we find
(14)Tn(t)=fne-k(nπ/L)2t+∫0te-k(nπ/L)2(t-s)φn(s)ds,u(x,t)=∑n=1∞fne-k(nπ/L)2tsin(nπxL)+∑n=1∞sin(nπxL)∫0te-k(nπ/L)2(t-s)φn(s)ds,u(x,t)-∑n=1∞sin(nπxL)∫0te-k(nπ/L)2(t-s)φn(s)ds=∫0LK(x,y,t)f(y)dy,
where K(x,y,t)=(2/L)∑n=1∞e-(nπ)2kt/L2sin⁡(nπx/L)sin(nπy/L. K(x,y,t) is an infinite series; numerically we cannot handle infinite sums. Limit the sum to a finite number of expansion terms 200 which guarantees the convergence of the series. So
(15)K(x,y,t)≈2L∑n=1200e-(nπ)2kt/L2sin⁡(nπxL)sin⁡(nπyL).
Thus initial inverse problem is reduced to solving integral equation of the first kind:
(16)w(x)=∫0LK(x,y,T)f(y)dy,
where
(17)w(x)=u(x,T)-∑n=1∞ansin(nπxL)=g(x)-∑n=1∞ansin(nπxL),an=2L∫0Te-k(nπ/L)2(T-s)[∫0Lφ(x,s)sin(nπxL)dx]ds=∫0L∫0T2Le-k(nπ/L)2(T-s)ϕ(x,s)sin(nπxL)dsdx=∫0L∫0Tg(x,s,n)dsdx,
where g(x,s,n)=(2/L)e-k(nπ/L)2(T-s)ϕ(x,s)sin(nπx/L). Since an is a double integral, it is difficult to determine the exact value. The authors try to seek its approximation by complex trapezoid formula.Given that the interval[0,L] is equally subdivided into m subintervals, {[xi-1,xi]}|i=1m, with equally spaced sampling points xi=ih where i=0,1,2,…,m. Also, assume that the interval [0,T] is equally subdivided into l subintervals, {[sj-1,sj]}|j=1l, with equally spaced sampling points sj=jk where j=0,1,2,…,l. The composite Trapezoidal rule is
(18)an≈14hk(∑i=1m-1g(xi,T,n)g(0,0,n)+g(L,0,n)+g(0,T,n)+g(L,T,n)mmm+2∑i=1m-1g(xi,0,n)+2∑i=1m-1g(xi,T,n)mmm+2∑j=1l-1g(0,sj,n)+2∑j=1l-1g(L,sj,n)mmm+4∑j=1l-1(∑i=1m-1g(xi,sj,n))).
Approximation of an can be obtained by the above formula.The first step in the numerical treatment used in this research consists in discretization of (16) by the quadrature formula. The interval [0,L] can be subdivided into N equal intervals of width h=Δy=L/N. Let y0=0 and yj=jΔy, and since the variable is either y or x, let x0=y0=0, xN=yN=L, and xi=iΔy (i.e., xi=yi). Also denote f(xi) as fi, w(xi) as wi, and K(xi,yj,T) as Kij. Now if the trapezoid rule is used to approximate the given equation, then
(19)w(x)=∫0LK(x,y,T)f(y)dy≈Δy[12K(x,y0,T)f(y0)+K(x,y1,T)f(y1)mmm+⋯+K(x,yN-1,T)f(yN-1)mmm+12K(x,yN,T)f(yN)],
or, more tersely,
(20)w(x)≈Δy[12K(x,y0,T)f0+K(x,y1,T)f1mim+⋯K(x,yN-1,T)fN-1+12K(x,yN,T)fN].
There are N+1 values of fi, as i=0,1,2,…,N. Therefore the equation becomes a set of N+1 equations in fi:
(21)wi=Δy[12Ki0f0+Ki1f1+⋯Ki(N-1)fN-1+12KiNfN],mmmmmmmmmmmmmmmmmmmmi=0,1,2,…,N,
which give the approximate solution to f(xi) at x=xi. This may also be written in matrix form:
(22)KF=W,
where K is the matrix of coefficients,(23)K=Δy·[12K(x0,y0,T)K(x0,y1,T)⋯12K(x0,yN,T)12K(x1,y0,T)K(x1,y1,T)⋯12K(x1,yN,T)⋯⋯⋯⋯12K(xN,y0,T)K(xN,y1,T)⋯12K(xN,yN,T)],F is the matrix of solutions,
(24)F=[f(y0),…,f(yi),…,f(yN)]T,
and W is the matrix of the nonhomogeneous part:
(25)W=[w(x0),…,w(xi),…,w(xN)]T,w(xi)=u(xi,T)-∑n=1∞ansin(nπxiL)=g(xi)-∑n=1∞ansin(nπxiL),i=0,1,2,…,N.The problem (22) is ill-posed in the sense that the inverse operator K-1 of K exists but is not continuous. Hence, although the problem (22) has a unique solution, solving it directly will not give a right solution. Indeed, the linear operator K is so badly conditioned that any numerical attempt to directly solve (22) may fail.In this note we describe the Tikhonov regularization for finding a stable approximate solution to a linear ill-posed problem represented in the form of an operator equation:(26)KF=W,
where, instead of the exact data W, noisy data Wδ is available with
(27)∥W-Wδ∥≤δ.
In order to find a solution in stable manner, Tikhonov proposed to solve the following:
(28)Fα,δ=argmin⁡F∈RN+1Jα(F)=(∥KF-Wδ∥22+α∥F∥22),
where α is the regularization parameter. The computation of the approximate solution Fα,δ consists in solving the Euler equation corresponding to the functional Jα. This equation has the following form:
(29)(K*K+αI)Fα,δ=K*Wδ,
where K* is the adjoint operator of K and I is the identity operator. The regularization solution Fα,δ obviously depends on the parameter α, the value of which directly affects the degree of approximation and the stability of the solution. From the viewpoint of the approximation degree, the smaller the better for the value of α; however, from the viewpoint of the solution stability, the larger the better. So the key point for solving the problem is getting the value of α.In order to obtain an effectively approximate solution to the original ill-posed problem, it is particularly important to choose suitable parameter. If we do not have the information ofW, we have to use an “a posteriori” parameter choice, where one uses results of the actual computation like the residual for computing a regularization parameter. The most widely used method is Morozov’s discrepancy principle, where one takes α such that
(30)∥KFα,δ-Wδ∥2=δ.On how to solve the equation∥KFα,δ-Wδ∥2=δ, some researchers have developed some theories and methods to determine α [13, 15–19]. In the present paper, the genetic algorithm for obtaining the regularization parameter is presented.


## 4. The Genetic Algorithm

Genetic algorithms belong to the general category of stochastic global optimization methods [20–22]. They have their philosophical basis in a process found in nature related to the evolution of the different biological species. Based on Darwin's theory of evolution where the fittest individuals survive, while the less fit are eliminated, a simple genetic algorithm is basically composed of three operators: (i) selection, (ii) crossover, and (iii) mutation. The computational implementation of a genetic algorithm starts with random generation of a population of individuals where each individual is actually a chain composed of representations of binary numbers, sequentially, random values for design variables. Each individual in the population corresponds to a value of the objective function (fitness function),
(31)λ(α)=|∥KFα,δ-Wδ∥2-δ|,
which is associated with a function of adaptability that is minimized.With the fitness function values calculated for each individual in the population, the probability that this individual will be selected for crossover is determined. This probability is given by(32)pi=λi∑i=1nλi,i=1,2,…,n,
where  λi is the fitness function, defined in (31), and n is the population size. The fittest individuals have higher probability of being selected for crossover. After selecting a pair of individuals, crossover is done by switching between the parts of the binary strings of these individuals. The population size remains constant through all the generations. This exchange-digit binary string corresponds in nature to the exchange of genetic information between parents, which represents the process by which it promotes the improvement of the species. The crossover for each pair of selected individuals is also subject to a probability pc to be performed. Mutation is a safeguard against losses of useful genetic material during reproduction and crossover. The process of mutation is simply to choose a few members from the population based on a probability of mutation pm and switch a 0 to a 1 and vice versa at a randomly selected site on the selected member. The process consisting of selection, crossover, and mutation is then repeated by the total number of generations, and at the end the most fitted individual is selected, that is, the one that leads to lower value of the objective function (31). The specific implementation process is shown in Figure 1.Figure 1
Flow chart for the calculation of regularization parameter based on GA.


## 5. Numerical Examples

In this section, we are going to demonstrate some numerical results for determiningf(x) in the inverse problems (1)–(3). We use the L∞ error norm and the relative error to measure the difference between the numerical and analytical solutions. The L∞ error norm is defined by
(33)L∞=max⁡0≤j≤N|f(xj)-f~(xj)|,
and the relative error (RE) is defined by
(34)RE=∑j=0N|f(xj)-f~(xj)|2∑j=0N|f(xj)|2,
where xj are test points and N is the total number of uniformly distributed points on [0,1]. f(x) is the exact solution and f~(x) is the numerical solution.In our computations, we takeN=40. The noisy data {wδ(xj)}|j=0N were assumed to contain some random errors. Assume the observed data has the following noised form:
(35)Wδ=W+δ×randn(size(W)),
where W=(w(x0),w(x1),…,w(x40))T, “randn(·)” is a normally distributed random variable with zero mean and unit standard deviation, and δ dictates the level of noise. “randn(size(W)” returns an array of random entries that has the same size as W.Example 1.
In this example let us consider the following inverse problem:(36)∂u∂t=∂2u∂x2+φ(x,t),0<x<1,t>0,u(0,t)=0,t>0,u(1,t)=0,t>0,u(x,0)=f(x),0≤x≤1.
The overspecified condition is as follows:
(37)u(x,T)=g(x)=e-Tsin(πx).
The analytical solution of this example is u(x,t)=e-tsin(πx) and φ(x,t)=(π2-1)e-tsin(πx), f(x)=sin(πx).The regularization parameterα is chosen using the genetic algorithm. In the genetic algorithm implemented, a tournament selection scheme by two individuals was used. A uniform crossover operator is applied with a probability of 0.95 and a uniform mutation operator where each individual has a probability 0.05 of being mutated. Population size is 20, number of generations is 100. The search domain is [0,100] for α. The objective function λ(α), the L∞ error norm, and relative error RE are presented in Tables 1, 2, 3, and 4 corresponding to T=0.1,0.5,1, and 1.2, respectively. Especially when T=1, the corresponding errors between the analytical values and the estimated results by function f(x) are listed in Table 5  (δ=0) and Table 6  (δ=0.1).Table 1
T
=
0.1, the objective function λ(α), the L∞, and the relative errors for f(x), with δ=0,0.01,0.05,0.1.
δ
α
λ
(
α
)
L
∞
RE
0
6.7341
×
10
-
10
8.6679
×
10
-
9
5.8643
×
10
-
6
5.8822
×
10
-
6
0.01
1.9422
×
10
-
3
0
1.6741
×
10
-
3
1.5696
×
10
-
3
0.05
9.7098
×
10
-
3
0
1.9019
×
10
-
3
1.7996
×
10
-
3
0.1
1.9419
×
10
-
2
0
1.9352
×
10
-
3
1.8332
×
10
-
3Table 2
T
=
0.5, the objective function λ(α), the L∞, and the relative errors for f(x), with δ=0,0.01,0.05,0.1.
δ
α
λ
(
α
)
L
∞
RE
0
3.1837
×
10
-
13
1.7318
×
10
-
7
7.7189
×
10
-
8
4.9543
×
10
-
8
0.01
3.7471
×
10
-
5
0
3.2502
×
10
-
6
3.2502
×
10
-
6
0.05
1.8735
×
10
-
4
0
3.2502
×
10
-
6
3.2502
×
10
-
6
0.1
3.7471
×
10
-
4
0
3.2502
×
10
-
6
3.2502
×
10
-
6Table 3
T
=
1, the objective function λ(α), the L∞, and the relative errors for f(x), with δ=0,0.01,0.05,0.1.
δ
α
λ
(
α
)
L
∞
RE
0
8.9302
×
10
-
18
3.7471
×
10
-
6
8.1049
×
10
-
6
8.0927
×
10
-
6
0.01
2.6955
×
10
-
7
0
2.5360
×
10
-
4
2.5360
×
10
-
4
0.05
1.3477
×
10
-
6
0
2.5464
×
10
-
4
2.5464
×
10
-
4
0.1
2.6955
×
10
-
6
0
2.5477
×
10
-
4
2.5477
×
10
-
4Table 4
T
=
1.2, the objective function λ(α), the L∞, and the relative errors for f(x), with δ=0,0.01,0.05,0.1.
δ
α
λ
(
α
)
L
∞
RE
0
6.1296
×
10
-
16
2.9210
×
10
-
6
7.4720
×
10
-
7
7.4718
×
10
-
7
0.01
3.7452
×
10
-
8
0
4.8134
×
10
-
4
4.8134
×
10
-
4
0.05
1.8726
×
10
-
7
0
4.8590
×
10
-
4
4.8590
×
10
-
4
0.1
3.7453
×
10
-
7
0
4.8647
×
10
-
4
4.8647
×
10
-
4Table 5
T
=
1, the analytical and numerical results for the f(xj) in xj=0.1j, where δ=0.
j
Exactf(xj)
Numericalf~(xj)
|
f
(
x
j
)
-
f
~
(
x
j
)
|
1
0.309016994374947
0.309019506872923
2.51249797578179
×
10
-
6
2
0.587785252292473
0.587790008583465
4.75629099216945
×
10
-
6
3
0.809016994374947
0.809023584278051
6.58990310320462
×
10
-
6
4
0.951056516295154
0.951064106758891
7.59046373755012
×
10
-
6
5
1
1.00000810910255
8.10910254811503
×
10
-
6
6
0.951056516295154
0.951064262563094
7.74626794020072
×
10
-
6
7
0.809016994374947
0.809023549839842
6.5554648942534
×
10
-
6
8
0.587785252292473
0.587789998909465
4.74661699212398
×
10
-
6
9
0.309016994374948
0.309019477607716
2.48323276808016
×
10
-
6
10
1.22464679914735
×
10
-
16
6.12328388409843
×
10
-
17
6.1231841073751
×
10
-
17Table 6
T
=
1, the analytical and numerical results for the f(xj) in xj=0.1j, where δ=0.1.
j
Exactf(xj)
Numericalf~(xj)
|
f
(
x
j
)
-
f
~
(
x
j
)
|
1
0.309016994374947
0.308938263749532
7.8730625415413
×
10
-
5
2
0.587785252292473
0.587635497743805
1.49754548667946
×
10
-
4
3
0.809016994374947
0.808810874921651
2.06119453296294
×
10
-
4
4
0.951056516295154
0.950814208345435
2.42307949718912
×
10
-
4
5
1
0.999745222344235
2.54777655765204
×
10
-
4
6
0.951056516295154
0.950814208345429
2.42307949724685
×
10
-
4
7
0.809016994374947
0.808810874921642
2.06119453305731
×
10
-
4
8
0.587785252292473
0.587635497743796
1.49754548676939
×
10
-
4
9
0.309016994374948
0.308938263749526
7.87306254214082
×
10
-
5
10
1.22464679914735
×
10
-
16
6.12167393253357
×
10
-
17
6.12479405893996
×
10
-
17From this numerical example, it can be seen that the numerical results are quite satisfactory. As shown in Tables1–6, when the noise level of input data δ=0 with T=0.1,0.5,1, and 1.2, the numerical solutions obtained using our method are of high degree of accuracy. Even with the noise level of input data up to δ=0.1, the numerical solutions are still in good agreement with the exact solutions.


## 6. Conclusion

This paper deals with the effective algorithms for solving the backward heat problem and the following results are obtained.(1)
The present study successfully applies the numerical method involving Tikhonov's regularization method in conjunction with the first Fredholm integral equation to the inverse nonhomogeneous heat conduction problems.(2)
From the illustrated example, it can be seen that the proposed numerical method is efficient and accurate to estimate the initial conditionf(x).



---

*Source: 529618-2014-03-02.xml*