Inverse Estimates for Nonhomogeneous Backward Heat Problems

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.

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 Section 2, 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 by 𝑓(𝑥).

The mathematical formulation of this problem is given by the following nonhomogeneous heat equation: 𝜕𝑢𝜕𝜕𝑡−𝑘²𝑢𝜕𝑥²()𝑢()𝑢()𝑢()(𝑥)=𝜑𝑥,𝑡,0<𝑥<𝐿,𝑡>0,0,𝑡=0,𝑡>0,𝐿,𝑡=0,𝑡>0,𝑥,0=𝑓,0≤𝑥≤𝐿, (1)(2) where 𝑢(𝑥,𝑡) (temperature), 𝜑(𝑥,𝑡) (source term), 𝑓(𝑥) (initial condition), 𝑥 (spatial variable), and 𝑡 (time variable) are dimensionless quantities and 𝑘 denotes the dispersion coefficient.

For the direct problem where the initial condition 𝑓(𝑥) is specified, the problem given by (1) and (2) is concerned with the determination of the temperature distribution 𝑢(𝑥,𝑡) in the interior region of the solid as a function of time and position.

2.2. Inverse Problem

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

3. Algorithm Analysis

The solution of the direct problem for a given initial condition 𝑓(𝑥) is explicitly obtained using separation of variables for 0<𝑥<𝐿,𝑡≥0. Consider the initial boundary problems (1) and (2); we look for solutions 𝑢 in the form 𝑢(𝑥,𝑡)=𝑇(𝑡)𝑋(𝑥). We look at the following eigenvalue problem: 𝑋^′′𝑋(0)(𝐿)+𝜆𝑋=0=𝑋=0. (4) The eigenvalues and the corresponding eigenfunctions are 𝜆=𝑛𝜋𝐿²,𝑋_𝑛(𝑥)=sin𝑛𝜋𝐿𝑥. (5) Now we set 𝑢(𝑥,𝑡)=^∞_𝑛=1𝑇_𝑛(𝑡)sin𝑛𝜋𝐿𝑥. (6) Formally computing 𝑢_𝑡 and 𝑢_𝑥𝑥 and substituting into (1), we get 𝜑(𝑥,𝑡)=𝑢_𝑡−𝑘𝑢_𝑥𝑥=^∞_𝑛=1𝑇^′_𝑛(𝑡)+𝑘𝑛𝜋𝐿²𝑇_𝑛(𝑡)sin𝑛𝜋𝑥𝐿. (7) Hence expanding 𝑓 and 𝜙 into the following Fourier series: 𝜑(𝑥,𝑡)=^∞_𝑛=1𝜑_𝑛(𝑡)sin𝑛𝜋𝑥𝐿,𝑓(𝑥)=^∞_𝑛=1𝑓_𝑛sin𝑛𝜋𝑥𝐿, (8) where 𝜑_𝑛(𝑡)=2𝐿^𝐿₀𝜑(𝑥,𝑡)sin𝑛𝜋𝑥𝐿𝑓𝑑𝑥,_𝑛=2𝐿^𝐿₀𝑓(𝑥)sin𝑛𝜋𝑥𝐿𝑑𝑥, (9) we get ^∞_𝑛=1𝑇^′_𝑛(𝑡)+𝑘𝑛𝜋𝐿²𝑇_𝑛(𝑡)sin𝑛𝜋𝑥𝐿=^∞_𝑛=1𝜑_𝑛(𝑡)sin𝑛𝜋𝑥𝐿. (10) The uniqueness of the Fourier expansion leads to the family of ODE's:𝑇^′_𝑛(𝑡)+𝑘𝑛𝜋𝐿²𝑇_𝑛(𝑡)=𝜑_𝑛(𝑡). (11) In addition, 𝑢(𝑥,0)=^∞_𝑛=1𝑇_𝑛(0)sin𝑛𝜋𝐿𝑥=^∞_𝑛=1𝑓_𝑛sin𝑛𝜋𝑥𝐿=𝑓(𝑥), (12) so that 𝑇_𝑛(0)=𝑓_𝑛,𝑛≥1. (13) Solving (11)–(13), we find 𝑇_𝑛(𝑡)=𝑓_𝑛𝑒^{−𝑘(𝑛𝜋/𝐿)2𝑡}+^𝑡₀𝑒^{−𝑘(𝑛𝜋/𝐿)2(𝑡−𝑠)}𝜑_𝑛(𝑠)𝑢()=𝑑𝑠,𝑥,𝑡^∞_𝑛=1𝑓_𝑛𝑒^{−𝑘(𝑛𝜋/𝐿)2𝑡}sin𝑛𝜋𝑥𝐿+^∞_𝑛=1sin𝑛𝜋𝑥𝐿^𝑡₀𝑒^{−𝑘(𝑛𝜋/𝐿)2(𝑡−𝑠)}𝜑_𝑛(𝑠)𝑢()−𝑑𝑠,𝑥,𝑡^∞_𝑛=1sin𝑛𝜋𝑥𝐿^𝑡₀𝑒^{−𝑘(𝑛𝜋/𝐿)2(𝑡−𝑠)}𝜑_𝑛(𝑠)=𝑑𝑠_𝐿0𝐾𝑥,𝑦,𝑡𝑓𝑦𝑑𝑦, (14) where 𝐾(𝑥,𝑦,𝑡)=(2/𝐿)∑^∞_𝑛=1𝑒^{−(𝑛𝜋)2𝑘𝑡/𝐿2}sin(𝑛𝜋𝑥/𝐿)sin(𝑛𝜋𝑦/𝐿. 𝐾(𝑥,𝑦,𝑡) 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 𝐾≈2𝑥,𝑦,𝑡𝐿²⁰⁰_𝑛=1𝑒^{−(𝑛𝜋)2𝑘𝑡/𝐿2}sin𝑛𝜋𝑥𝐿sin𝑛𝜋𝑦𝐿. (15) Thus initial inverse problem is reduced to solving integral equation of the first kind: 𝑤(𝑥)=^𝐿₀𝐾𝑥,𝑦,𝑇𝑓𝑦𝑑𝑦, (16) where 𝑤(𝑥)=𝑢(𝑥,𝑇)−^∞_𝑛=1𝑎_𝑛sin𝑛𝜋𝑥𝐿=𝑔(𝑥)−^∞_𝑛=1𝑎_𝑛sin𝑛𝜋𝑥𝐿,𝑎_𝑛=2𝐿^𝑇₀𝑒^{−𝑘(𝑛𝜋/𝐿)2(𝑇−𝑠)𝐿}₀𝜑(𝑥,𝑠)sin𝑛𝜋𝑥𝐿=𝑑𝑥𝑑𝑠_𝐿0𝑇02𝐿𝑒^{−𝑘(𝑛𝜋/𝐿)2(𝑇−𝑠)}𝜙(𝑥,𝑠)sin𝑛𝜋𝑥𝐿=𝑑𝑠𝑑𝑥_𝐿0𝑇0𝑔(𝑥,𝑠,𝑛)𝑑𝑠𝑑𝑥, (17) where 𝑔(𝑥,𝑠,𝑛)=(2/𝐿)𝑒^{−𝑘(𝑛𝜋/𝐿)2(𝑇−𝑠)}𝜙(𝑥,𝑠)sin(𝑛𝜋𝑥/𝐿). Since 𝑎_𝑛 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,𝐿] is equally subdivided into 𝑚 subintervals, {[𝑥_𝑖−1,𝑥_𝑖]}|^𝑚_𝑖=1, with equally spaced sampling points 𝑥_𝑖=𝑖ℎ where 𝑖=0,1,2,…,𝑚. Also, assume that the interval [0,𝑇] is equally subdivided into 𝑙 subintervals, {[𝑠_𝑗−1,𝑠_𝑗]}|^𝑙_𝑗=1, with equally spaced sampling points 𝑠_𝑗=𝑗𝑘 where 𝑗=0,1,2,…,𝑙. The composite Trapezoidal rule is 𝑎_𝑛≈14ℎ𝑘𝑔(0,0,𝑛)+𝑔(𝐿,0,𝑛)+𝑔(0,𝑇,𝑛)+𝑔(𝐿,𝑇,𝑛)+2^𝑚−1_𝑖=1𝑔𝑥_𝑖,0,𝑛+2^𝑚−1_𝑖=1𝑔𝑥_𝑖,𝑇,𝑛+2^𝑙−1_𝑗=1𝑔0,𝑠_𝑗,𝑛+2^𝑙−1_𝑗=1𝑔𝐿,𝑠_𝑗,𝑛+4^𝑙−1_𝑗=1^𝑚−1_𝑖=1𝑔𝑥_𝑖,𝑠_𝑗,𝑛. (18) Approximation of 𝑎_𝑛 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,𝐿] can be subdivided into 𝑁 equal intervals of width ℎ=Δ𝑦=𝐿/𝑁. Let 𝑦₀=0 and 𝑦_𝑗=𝑗Δ𝑦, and since the variable is either 𝑦 or 𝑥, let 𝑥₀=𝑦₀=0, 𝑥_𝑁=𝑦_𝑁=𝐿, and 𝑥_𝑖=𝑖Δ𝑦 (i.e., 𝑥_𝑖=𝑦_𝑖). Also denote 𝑓(𝑥_𝑖) as 𝑓_𝑖, 𝑤(𝑥_𝑖) as 𝑤_𝑖, and 𝐾(𝑥_𝑖,𝑦_𝑗,𝑇) as 𝐾_𝑖𝑗. Now if the trapezoid rule is used to approximate the given equation, then 𝑤(𝑥)=^𝐿₀𝐾𝑓𝑦1𝑥,𝑦,𝑇𝑑𝑦≈Δ𝑦2𝐾𝑥,𝑦₀,𝑇𝑓𝑦₀+𝐾𝑥,𝑦₁,𝑇𝑓𝑦₁+⋯+𝐾𝑥,𝑦_𝑁−1,𝑇𝑓𝑦_𝑁−1+12𝐾𝑥,𝑦_𝑁,𝑇𝑓𝑦_𝑁, (19) or, more tersely, 𝑤(𝑥)1≈Δ𝑦2𝐾𝑥,𝑦₀,𝑇𝑓₀+𝐾𝑥,𝑦₁,𝑇𝑓₁+⋯𝐾𝑥,𝑦_𝑁−1,𝑇𝑓_𝑁−1+12𝐾𝑥,𝑦_𝑁,𝑇𝑓_𝑁. (20) There are 𝑁+1 values of 𝑓_𝑖, as 𝑖=0,1,2,…,𝑁. Therefore the equation becomes a set of 𝑁+1 equations in 𝑓_𝑖: 𝑤_𝑖1=Δ𝑦2𝐾_𝑖0𝑓₀+𝐾_𝑖1𝑓₁+⋯𝐾_{𝑖(𝑁−1)}𝑓_𝑁−1+12𝐾_𝑖𝑁𝑓_𝑁,𝑖=0,1,2,…,𝑁, (21) which give the approximate solution to 𝑓(𝑥_𝑖) at 𝑥=𝑥_𝑖. This may also be written in matrix form: 𝐾𝐹=𝑊, (22) where 𝐾 is the matrix of coefficients,1𝐾=Δ𝑦⋅2𝐾𝑥₀,𝑦₀,𝑇𝐾𝑥₀,𝑦₁⋯1,𝑇2𝐾𝑥₀,𝑦_𝑁1,𝑇2𝐾𝑥₁,𝑦₀,𝑇𝐾𝑥₁,𝑦₁⋯1,𝑇2𝐾𝑥₁,𝑦_𝑁1,𝑇⋯⋯⋯⋯2𝐾𝑥_𝑁,𝑦₀,𝑇𝐾𝑥_𝑁,𝑦₁⋯1,𝑇2𝐾𝑥_𝑁,𝑦_𝑁,𝑇, (23)

𝐹 is the matrix of solutions, 𝐹=𝑓𝑦₀,…,𝑓𝑦_𝑖,…,𝑓𝑦_𝑁^𝑇, (24) and 𝑊 is the matrix of the nonhomogeneous part: 𝑊=𝑤𝑥₀,…,𝑤𝑥_𝑖,…,𝑤𝑥_𝑁^𝑇,𝑤𝑥_𝑖=𝑢𝑥_𝑖,𝑇−^∞_𝑛=1𝑎_𝑛sin𝑛𝜋𝑥_𝑖𝐿=𝑔𝑥_𝑖−^∞_𝑛=1𝑎_𝑛sin𝑛𝜋𝑥_𝑖𝐿,𝑖=0,1,2,…,𝑁. (25)

The problem (22) is ill-posed in the sense that the inverse operator 𝐾⁻¹ of 𝐾 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 𝐾 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) where, instead of the exact data 𝑊, noisy data 𝑊_𝛿 is available with 𝑊−𝑊_𝛿≤𝛿. (27) In order to find a solution in stable manner, Tikhonov proposed to solve the following: 𝐹_𝛼,𝛿=argmin_{𝐹∈𝑅𝑁+1}𝐽_𝛼(𝐹)=𝐾𝐹−𝑊_𝛿²₂+𝛼‖𝐹‖²₂, (28) where 𝛼 is the regularization parameter. The computation of the approximate solution 𝐹_𝛼,𝛿 consists in solving the Euler equation corresponding to the functional 𝐽_𝛼. This equation has the following form: 𝐾^∗𝐾+𝛼𝐼𝐹_𝛼,𝛿=𝐾^∗𝑊_𝛿, (29) where 𝐾^∗ is the adjoint operator of 𝐾 and 𝐼 is the identity operator. The regularization solution 𝐹_𝛼,𝛿 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 of  𝑊, 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 𝐾𝐹_𝛼,𝛿−𝑊_𝛿2=𝛿. (30)

On how to solve the equation ‖𝐾𝐹_𝛼,𝛿−𝑊_𝛿‖₂=𝛿, 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), 𝜆(𝛼)=𝐾𝐹_𝛼,𝛿−𝑊_𝛿2−𝛿, (31) 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 𝑝_𝑖=𝜆_𝑖∑^𝑛_𝑖=1𝜆_𝑖,𝑖=1,2,…,𝑛, (32) where  𝜆_𝑖 is the fitness function, defined in (31), and 𝑛 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 𝑝_𝑐 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 𝑝_𝑚 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 determining 𝑓(𝑥) in the inverse problems (1)–(3). We use the 𝐿_∞ error norm and the relative error to measure the difference between the numerical and analytical solutions. The 𝐿_∞ error norm is defined by 𝐿_∞=max_{0≤𝑗≤𝑁}𝑓𝑥_𝑗−𝑓𝑥_𝑗, (33) and the relative error (RE) is defined by RE=∑^𝑁_𝑗=0𝑓(𝑥_𝑗)−𝑓(𝑥_𝑗)²∑^𝑁_𝑗=0𝑓(𝑥_𝑗)², (34) where 𝑥_𝑗 are test points and 𝑁 is the total number of uniformly distributed points on [0,1]. 𝑓(𝑥) is the exact solution and 𝑓(𝑥) is the numerical solution.

In our computations, we take 𝑁=40. The noisy data {𝑤^𝛿(𝑥_𝑗)}|^𝑁_𝑗=0 were assumed to contain some random errors. Assume the observed data has the following noised form: 𝑊^𝛿=𝑊+𝛿×randn(size(𝑊)), (35) where 𝑊=(𝑤(𝑥₀),𝑤(𝑥₁),…,𝑤(𝑥₄₀))^𝑇, “randn(⋅)” is a normally distributed random variable with zero mean and unit standard deviation, and 𝛿 dictates the level of noise. “randn(size(𝑊)” returns an array of random entries that has the same size as 𝑊.

Example 1. In this example let us consider the following inverse problem: 𝜕𝑢=𝜕𝜕𝑡²𝑢𝜕𝑥²()𝑢()𝑢()𝑢()(𝑥)+𝜑𝑥,𝑡,0<𝑥<1,𝑡>0,0,𝑡=0,𝑡>0,1,𝑡=0,𝑡>0,𝑥,0=𝑓,0≤𝑥≤1. (36) The overspecified condition is as follows: 𝑢(𝑥,𝑇)=𝑔(𝑥)=𝑒^−𝑇sin(𝜋𝑥). (37) The analytical solution of this example is 𝑢(𝑥,𝑡)=𝑒^−𝑡sin(𝜋𝑥) and 𝜑(𝑥,𝑡)=(𝜋²−1)𝑒^−𝑡sin(𝜋𝑥), 𝑓(𝑥)=sin(𝜋𝑥).

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 𝐿_∞ error norm, and relative error RE are presented in Tables 1, 2, 3, and 4 corresponding to 𝑇=0.1,0.5,1, and 1.2, respectively. Especially when 𝑇=1, the corresponding errors between the analytical values and the estimated results by function 𝑓(𝑥) are listed in Table 5  (𝛿=0) and Table 6  (𝛿=0.1).

Table 1 

𝑇=0.1, the objective function 𝜆(𝛼), the 𝐿_∞, and the relative errors for 𝑓(𝑥), with 𝛿=0,0.01,0.05,0.1.

Table 2 

𝑇=0.5, the objective function 𝜆(𝛼), the 𝐿_∞, and the relative errors for 𝑓(𝑥), with 𝛿=0,0.01,0.05,0.1.

Table 3 

𝑇=1, the objective function 𝜆(𝛼), the 𝐿_∞, and the relative errors for 𝑓(𝑥), with 𝛿=0,0.01,0.05,0.1.

Table 4 

𝑇=1.2, the objective function 𝜆(𝛼), the 𝐿_∞, and the relative errors for 𝑓(𝑥), with 𝛿=0,0.01,0.05,0.1.

Table 5 

𝑇=1, the analytical and numerical results for the 𝑓(𝑥_𝑗) in 𝑥_𝑗=0.1𝑗, where 𝛿=0.

Table 6 

𝑇=1, the analytical and numerical results for the 𝑓(𝑥_𝑗) in 𝑥_𝑗=0.1𝑗, where 𝛿=0.1.

From this numerical example, it can be seen that the numerical results are quite satisfactory. As shown in Tables 1–6, when the noise level of input data 𝛿=0 with 𝑇=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 condition 𝑓(𝑥).

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The work of the authors is supported by the Special Funds of the National Natural Science Foundation of China (nos. 51190093 and 51179151). The authors would like to thank the referees for constructive suggestions and comments.