Scaling, Self-Similarity, and Systems of Fractional OrderView this Special Issue
Research Article | Open Access
A Discretized Tikhonov Regularization Method for a Fractional Backward Heat Conduction Problem
We propose a numerical reconstruction method for solving a time-fractional backward heat conduction problem. Based on the idea of reproducing kernel approximation, we reconstruct the unknown initial heat distribution from a finite set of scattered measurements of transient temperature at a fixed final time. The standard Tikhonov regularization technique using the norm of reproducing the kernel Hilbert space as the penalty term is adopted to provide a stable solution when the measurement data contains noise. Numerical results indicate that the proposed method is efficient.
Let , , be a bounded domain with sufficiently smooth boundary . Consider the following initial boundary value problem for time-fractional diffusion equation (TFDE): Here denotes the Caputo fractional derivative with respect to and is defined by is the Gamma function and is a symmetric uniformly elliptic operator and is a fixed final time, and are constants, and is the outward unit normal vector of the domain . In what follows, let and be given by where , for and , , . Moreover, we assume that there exists a positive constant such that
Recently, people are shifting their partial focus to fractional-order differential equations (FDEs) with the realization that the use of fractional-order derivatives and integrals leads to formulas of certain physical processes (for instance, some anomalous diffusion processes) which is more economical and useful than the classical approach in terms of Fick’s laws of diffusion. In this paper, we consider the fractional-order partial differential equation (FPDE) in (1), which is obtained from the standard diffusion equation by replacing the first-order time derivative with a fractional derivative of order , with . Different models using this kind of FDEs have been proposed [1–3], and there has been significant interest in developing numerical schemes for their solution. Physically, the time-fractional partial differential equations describe the continuous time random walk problems (the non-Markovian process). The physical interpretation of the fractional derivative is that it represents a degree of memory in the diffusing material . Actually, the convolution integral in the definition of the Caputo fractional-order derivatives for at time requires all the knowledge of classical derivative for , which reflects the “memory effect” of fractional derivatives. The utilization of the memory effect of fractional derivatives comes with a high cost regarding numerical solvability. Any algorithm using a discretization of a noninteger derivative has to take into account its nonlocal structure which means in general a high storage requirement and great overall complexity of the algorithm. Similar convolution model can be used to depict the tumor-immune system . Fractional dimensional model is also used to characterize the binary images of DNA  and oscillators [7, 8]. In , a general approach is proposed to approximate ideal filters based on fractional calculus from the point of view of systems of fractional order. Recently, numerous attempts to solve TFDE can be found in the literature [4, 10–14]. We can refer to [15–18] for more specified knowledge about fractional calculus.
The backward problem governed by time-fractional partial differential equation in (1) is to recover the heat distribution at any earlier time from the measurement , written as , which is noise-contaminated data for the exact temperature : for some known noise level . Such inverse problems have been considered by several authors. Based on the eigenfunction expansions, Sakamoto and Yamamoto  establish the unique existence of the weak solution and the asymptotic behavior as the time goes to for the forward problem and prove the stability and uniqueness in the backward problem in time. For the one-dimensional case, Liu and Yamamoto  propose a regularizing scheme by the quasi-reversibility to restore the stability of the backward problem. In , a regularization by projection is applied to the same problem as in  and the corresponding convergence rates are obtained under a priori and a posteriori parameter choice strategies, respectively. Here, we pay our attention to the situation of stable reconstruction of the initial heat distribution from some scattered noisy data of . More specifically, the data are collected only at a finite set of points . We then reconstruct the initial temperature distribution from the scattered noisy data at . For solving the backward diffusion problem, we employ a discretized Tikhonov regularization by the Ritz approach coupled with the reproducing kernel Hilbert space (RKHS), which is proposed in .
The Tikhonov regularization method has been widely studied and applied to all varieties of ill-posed problems [23, 24]. The discretized Tikhonov regularization method and its relative theories are also explored in detail . We adopt the Tikhonov regularization method by a reproducing kernel Hilbert space into the backward problem (1). As we know, the theory and practice of reproducing kernel are a fast growing research area. The numerical methods by RKHS have been also rapidly developed in recent years [25, 26]. These developments are due to the increasing interest in the use of reproducing kernel for the solution of mathematical and engineering problems, for instance, machine learning , signal processing , stochastic processes , wavelet transforms , and so forth. For the details of RKHS, we are able to refer to . However, to the authors’ knowledge, there are few applications of RKHS to inverse problems. We provide the partial list of the recent works. Takeuchi and Yamamoto  prove the convergence of the discretized Tikhonov regularization method by RKHS. Hon and Takeuchi  apply this method into a backward heat conduction problem for parabolic-type partial differential equations. In reproducing the kernel Hilbert space settings, an inverse source identification problem for parabolic equations is considered in . In , Saitoh discusses comprehensively the corresponding applications of RKHS in inverse problems.
The remainder of this paper is composed of five sections. In Section 2, we discuss Green’s function for problem (1) and use it to construct the reproducing kernel. In Section 3, we state the reconstruction method of the fractional backward diffusion problem by using the reproducing kernel. In Section 4, some numerical examples are given to illustrate the effectiveness of our method. A summary is made in the Section 5. Finally, we list some existing knowledge about the reproducing kernel Hilbert space in the Appendix.
2. Green’s Function and the Reproducing Kernel
In this section, we explore Green’s function of system (1) and use it to construct the reproducing kernel. Firstly, let function be Green’s function of the system (1); that is, satisfies the following equations in distribution’s sense: It is easy to know that problem (6) is equivalent to the following problem: Let us employ the Laplace transform to solve the system (6). The Laplace transform of a function on is defined by The Laplace transform of the Caputo fractional derivative is given by  where . The Caputo fractional derivative appears more suitable to be treated by the Laplace transforming in that it requires the knowledge of initial values of the function and of its integer derivatives of order . By the Laplace transforming about the time variable , the system (6) becomes where denotes the Laplace transform of Green’s function . It should be noted here that the system in (1) is only defined on , not on . When we make the Laplace transform, some necessary preprocess, for example, the function continuation technique, needs to be done on the solution of (1) to satisfy the condition of the Laplace transform. Because we do not use the value of for , the condition can be satisfied easily.
Applying this technique of eigenfunction expansions to problems (10), we have that where is the th orthonormal eigenfunction and is the corresponding eigenvalue to the Sturm-Liouville problem subject to the boundary conditions Taking the Laplace inverse of (11), we have that where is the Mittag-Leffler function defined by For the details of the Mittag-Leffler function, one can refer to . Subsequently, one can easily verify that the unique solution of system (1) with initial value can be written as where is Green’s function defined by (14). With the aid of (16), for each , introducing the operator by we may formulate the inverse problem as an integral equation of the first kind
The symmetry of Green’s function indicates that the operator is self-adjoint. We proceed by giving a brief account on the close connections of the ill posedness of operator equation (18) with the singular value system of operator . According to the symmetry of operator , we only need to discuss its eigensystem. It is easy to know that the eigensystem of is given by where is the eigenvalue and is the corresponding eigenfunction. We can see the decay of the eigenvalues with the increase of from the following asymptotic behavior of the Mittag-Leffler function which can be found in [20, 21] or by the results in .
Lemma 1 (see ). Let . Then there exist constants depending only on and such that These estimates are uniform for all .
Next, we consider the construction of reproducing kernel using Green’s function . Define for some . The symmetry about the space variable , of Green’s function indicates that is symmetric. Now that is a symmetric positive definite kernel, a unique RKHS in which the given kernel acts as the reproducing kernel can be constructed (see  for details). Henceforth we denote by the RHKS generated by the kernel . Actually, according to [26, Chapter 10], the inner product and norm on are defined by respectively. The space is actually given by The second inequalities in Lemma 1 show that, as the function of , , where denotes the collection of slowly increasing functions  defined by Hence, according to the theoretical results of [26, Chapter 10], we assert that the space is consistent with some Sobolev space for and the norm on is equivalent to the norm .
3. Formulation of the Inverse Problem and the Reconstruction Method
In order to find the initial temperature distribution , we would like to determine the solution of minimization problem However, the minimization element of (25) generally is a poor approximation of the desired initial function due to the error in and the ill posedness of operator equation (18). The Tikhonov regularization replaces the minimization problem (25) by the solution of a penalized least-squares problem with regularization parameter . We can see from [23, Proposition 3.11] that the convergence of any regularization method can be arbitrarily slow in general. Actually, convergence rates can be given only on some subset of , that is, under a priori assumptions on the exact data. Here, we assume the exact solution belongs to the set of source conditions where is the a priori bound and is a constant. As in , we know that there exists a constant , named the “qualification” of the regularization method, such that . For the Tikhonov regularization method, the qualification . However, according to (21), it is easy to know that as and . Therefore, here we only need consider the case of for the Tikhonov method (26).
In order to solve the minimization problem (26), some discretization scheme needs to be given. A natural way to obtain such discretization is to generate a finite dimensional approximation to the minimal element of the Tikhonov functional . For this, we define a subspace , where . This approximation by discretization is equivalent to finding the minimal norm least-squares solution of the equation where is the projection operator. Moreover, we produce a finite dimensional approximation to by minimizing the Tikhonov functional (26) over the finite dimensional space . Denote by the minimizer of for noise input data . It is well known that satisfy [23, 24] Provided is an expanding sequence, the convergence of the Tikhonov regularized solutions is proved . Takeuchi and Yamamoto  show that the discretized Tikhonov regularized solutions converge to the exact solution without the monotonicity of under an a priori choice strategy for and . However, one can see that, from the existing results, in both cases, the regularized solutions converging to the exact solution depends on whether , . Moreover, the convergence of requires the knowledge of the fill distance , which is defined by  The fill distance can be interpreted in various geometrical ways. For example, we can consider it as the radius of the largest ball which is completely contained in and which does not contain a data site. In this sense describes the largest data-site-free hole in .
Here, we utilize the same proof as in  to give the following lemma.
Lemma 2. Consider as .
It is obvious that the kernel function is sufficiently smooth. Therefore, according to the error estimate in the Appendix, we can find a positive constant such that the estimate holds for all . Meanwhile, it is easy to know where is a constant. Moreover, the property of RKHS, , leads to Combining (32)–(34), we get where is a constant. Consequently, substituting the estimate (35) into (31) and letting , we complete the proof.
According to the classical results on the Tikhonov regularization for linear ill-posed problem (see, e.g., [23, 24]) and in view of (1), it holds that As in , we can estimate the noise-free term as follows: where is a constant. In view of the best possible error bound being , the term has to be chosen such that From the above discussions, we have the following theorem.
4. Numerical Tests
In this section, we present numerical results to illustrate the feasibility of the reconstruction method as described in the previous section.
In practical situation, we only can get the scattered noisy data of , that is, . As a result, instead of solving (25) we intend to deal with the following problem: where Since , the minimizer can be written as From the definition of RKHS, it follows that In addition, we know that for and . Now, it is easy to see that the coefficient vector satisfies the following linear system: where is an matrix, is the transpose of , and is an matrix defined by If we truncate Green’s function to the former terms, matrix can be decomposed to the following product: where denotes the conjugate transpose of and is given by Moreover, we turn to search for the minimizer of the following functional: After obtaining the vector , we substitute it into (43) and then get the regularized approximation .
In our tests, the measurement points , which are randomly generated by using the Matlab function rand, are scattered in the domain . Now we generate the final measurement data at with noise by where are the measurement points and generates a standard dimensional random vector. To evaluate the proposed method, we compute the relative error of the reconstructed solutions denoted by : where denotes the norm and denotes the norm. Before we proceed, it is natural that we have to determine , , measurement points , and observation time to define , . According to the convergence theorem, smaller yields better numerical solution, which implies that we need to choose as large as possible. However, the ill posedness of the backward diffusion problem results in the ill condition of the matrix , which causes us not to implement the numerical computation of the inverse of when the regularizing term is also ill conditioned. We use the following one-dimensional example to depict the change trend of the condition number of , cond, with respect to .
Example 1. Consider the following Dirichlet boundary value problem:
The forward problem has a unique solution
where the coefficient
To clarify the numerical influence of some relative parameters but , we fix parameter firstly.
In this test, we fix firstly. In Figure 1, for the cases of , we plot cond versus the number of running from to , respectively. The displayed results in Figure 1 show that the condition number cond increases exponentially as increases. Nevertheless, we can remove such influence of the ill condition of on numerical computation through modifying its small singular value as a fixed small constant . And in doing so, the numerical precision does not change significantly. Therefore, we can implement the proposed method without worrying so much about the size limitation of . For the choice of , we compare the computational results for some different ’s and ’s using exact final data for Example 1. For this, a preset value for the regularization parameter needs to be provided. Here, we simply choose as in , where for matrix . In the following computation, the Matlab code developed by Hansen [36, 37] is used to obtain the approximation solution for solving the discrete system (46). In addition, note that the Mittag-Leffler function is numerically realized by implementing the Matlab toolbox by Podlubny . In Figure 2, we plot the relative errors for different , , versus the number of running from to for and . The computational results show that, using smaller , we have less relative error. In addition, we also see from Figure 2 that, with the increase of the number of , the relative error becomes smaller firstly and then it remains steady after arriving to a certain extent.
We also need to consider the effect of the number of truncation term of Green’s function and the number of measurement points on the numerical precision. With , we display the numerical results for several and as the noise level in Table 2. It can be seen that when and become sufficiently large, the relative errors almost remain at the level .
Next, in the absence of the a priori information, we only evaluate the proposed algorithm in (46) for noisy data by using L-curve parameter choice method instead of that in Theorem 3. In Example 1, we fix , , and . Table 1 reports the relative errors of for different noise levels and final measurement times . These numerical results for all noisy cases are satisfactory. In general, it can been seen from Table 1 that, at the smaller and , the numerical effects are better. In addition, when the measurement time , the exact solution and the numerical solution with the relative noise level are displayed in Figure 3. It can be observed that the method works even for the case of with noise level as well.
Finally, we hope to use Example 1 to show that the proposed algorithm is robust for order . For varying , we plot the relative error versus in Figure 4. The displayed results show that the numerical method is robust when is varying.
Example 2. Consider the following Neumann boundary value problem: The unique solution to (56) is given by In this example, we show the numerical results for under the setting , , and . Figure 5 shows that the proposed method is capable of giving satisfactory results for the case of the Neumann boundary condition.
Example 3. We consider a two-dimensional fractional diffusion problem with the Dirichlet boundary value in :
The initial distribution is to be recovered by using the exact solution We firstly display the choice of the set in Figure 6. According to the analysis for 1-dimensional case about the number of and , we only deal with the case of and . The measurement points are scattered in , which are generated using the Matlab command rand. We report the relative error of for different final times and noise levels in Table 3 for . The numerical comparison between exact solution and regularized solution is shown in Figure 7. The numerical results show that the proposed method is acceptable for the 2-dimensional example. We also consider the influence of varying on the numerical stability. The relative error versus is plotted in Figure 8, from which we can see that the proposed method is robust about parameter .
In this paper, in a reproducing kernel Hilbert space setting, we propose a numerical reconstruction method, namely, the discretized Tikhonov regularization method, to recover the initial temperature distribution of the backward fractional diffusion problem. The implementation of the proposed method is simple and easy. Numerical tests show that the method is efficient.
Reproducing the Kernel Hilbert Spaces and Positive Definite Kernels
Most of the material in this appendix can be found in the excellent monograph . For the readers’ convenience we would like to repeat the theoretical results of RKHS. We are interested in linear vector spaces consisting of functions defined on a connected domain of .
Definition A.1. Let be a Hilbert space consisting of functions . is called a reproducing kernel Hilbert space and a kernel is called a reproducing kernel for if (i) for all ,(ii) for all and all ,
where is the inner product of .
The reproducing kernel of a RKHS is uniquely determined. According to , is a RKHS if and only if the point evaluation functionals are continuous, that is, for all . Also  discloses the connection between RKHS and positive definite kernels. Here, we call symmetric if for all .
Definition A.2. A continuous symmetric is called positive definite on if, for all , all sets of pairwise distinct centers ; the quadratic form
If is a symmetric positive definite kernel, then a unique RKHS in which the given kernel acts as the reproducing kernel can be constructed. Now, it follows from the definition of RKHS that(i) for all ,(ii) for all in the form of with .
For a symmetric positive definite kernel, introduce integral operator by
By [26, Proposition 10.28], maps continuously into the RKHS and is the adjoint of the embedding operator of the RKHS into . For such an operator, Mercer’s theorem  shows that can be represented as where are the nonnegative eigenvalues and are the eigenfunctions of . This allows us to derive the final characterization for RKHS .
Theorem A.3. Suppose is a symmetric positive definite kernel on a compact set . Then the RKHS is given by and the inner product has the representation
For a finite set of points and , consider the finite sum as an approximation of , which is actually the interpolant function of . We also can consider in the following way: define a subspace . We define the projection operator by