A Directly Numerical Algorithm for a Backward Time-Fractional Diffusion Equation Based on the Finite Element Method

Abstract

We study a backward problem for a time-fractional diffusion equation, which is formulated into a regularized optimization problem. After solving a sequence of well-posed direct problems by the finite element method, a directly numerical algorithm is proposed for solving the regularized optimization problem. In order to obtain a reasonable regularization solution, we utilize the discrepancy principle with decreasing geometric sequence to choose regularization parameters. One- and two-dimensional examples are given to verify the efficiency and stability of the proposed method.

Zhousheng Ruan 1, 2, 3 and Zewen Wang 2 and Wen Zhang 1, 2

Academic Editor:Valery G. Yakhno

1, Fundamental Science on Radioactive Geology and Exploration Technology Laboratory, East China Institute of Technology, Nanchang, Jiangxi 330013, China
2, School of Science, East China Institute of Technology, Nanchang, Jiangxi 330013, China
3, School of Mathematics and Statistics, Wuhan University, Wuhan, Hubei 430072, China

Received 8 April 2014; Accepted 22 July 2014

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Nowadays, there is increasing attention on fractional diffusion equations which can be used to describe anomalous diffusion phenomena instead of classical diffusion process. These new fractional-order models are more efficient than the integer-order models, because the fractional-order derivatives and integrals enable the description of the memory and hereditary properties of different substance [1]. By an argument similar to the derivation of the classical diffusion equation from Brownian motion, one can derive a fractional diffusion equation from continuous-time random walk. For example, in paper [2] the authors illustrated a fractional diffusion with respect to a non-Markovian diffusion process, while the authors discussed continuous-time random walks on fractals in paper [3].

We notice that mathematical and numerical analysis of the direct problems of the time-fractional diffusion equations has aroused wide concern in recent years; see [4-10] and references therein. At the same time, the inverse problems for the time-fractional diffusion equations have attracted more and more attention, not only for theoretical analysis but also for popular applications. The authors concluded that there exists a unique weak solution for the backward time-fractional diffusion equation problem under the overdetermined condition [figure omitted; refer to PDF] in paper [4]. The authors of papers [11-13] considered the backward problem of the time-fractional diffusion equation and proposed, respectively, a quasi-reversibility method, an optimization method, and a data regularization method for reconstructing the initial value. Inverse source problems for time-fractional diffusion equations were studied by using the method of the eigenfunction expansion [14], the integral equation method [15], and the separation of variables method [16], respectively, for recovering the space-dependent or time-dependent source term. In [17], the authors recovered the temperature function from one measured temperature at one interior point of a one-dimensional semi-infinite fractional diffusion equation based on Dirichlet kernel mollification techniques. The authors studied an inverse problem of identifying a spatially varying potential term in a one-dimensional time-fractional diffusion equation from the flux measurements taken at a single fixed time corresponding to a given set of input sources in [18]. Recently, for determining the space-dependent source in a parabolic equation, the authors [19] proposed a regularized optimization method together with the linear model function method [19, 20] for choosing regularization parameters. Inspired by this noniterative optimization method, we develop it to solve the backward problem for a time-fractional diffusion equation in this paper.

Let [figure omitted; refer to PDF] be a constant such that [figure omitted; refer to PDF] . We consider the following time-fractional diffusion equation: [figure omitted; refer to PDF] with homogeneous boundary condition [figure omitted; refer to PDF] and initial condition [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is a bounded domain in [figure omitted; refer to PDF] and [figure omitted; refer to PDF] is symmetric uniformly elliptic operator given by [figure omitted; refer to PDF] that is, there exists a constant [figure omitted; refer to PDF] , such that [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . The coefficients satisfy [figure omitted; refer to PDF] Here, [figure omitted; refer to PDF] is the Caputo fractional derivative which is defined by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the Gamma function.

If the function [figure omitted; refer to PDF] and the coefficients in (1) are all known, problem (1)-(3) is the so-called direct problem that can be solved stably by the finite element method, the finite difference, the spectrum method, and so forth. Here, we focus on the backward problem; that is, we try to determine the initial value [figure omitted; refer to PDF] by the additional data [figure omitted; refer to PDF] which is the measurement of the exact value [figure omitted; refer to PDF] and satisfies [figure omitted; refer to PDF] for some known error level [figure omitted; refer to PDF] . As we all know, the backward problem is ill-posed, which means that the solution does not depend continuously on the given data and any small perturbation in the given data may cause large change to the solution. For overcoming the ill-posedness we will adopt Tikhonov regularization in our treatment.

The rest of the paper is organized as follows. In Section 2, we reformulate the direct problem in a weak and variational sense. Then we formulate the inverse problem into a regularized optimization problem in Section 3. In Section 4, we give implementations of the regularized optimization method. Finally, numerical results are given to illustrate the efficiency and stability of the proposed method.

2. Weak Form and Weak Solution

The weak form of problem (1)-(3) is finding [figure omitted; refer to PDF] such that [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , [figure omitted; refer to PDF] denotes the inner product in [figure omitted; refer to PDF] , and [figure omitted; refer to PDF]

Definition 1.

A function [figure omitted; refer to PDF] is said to be a weak solution of the direct problem (1)-(3) if [figure omitted; refer to PDF] ; [figure omitted; refer to PDF] and the weak form (8) is satisfied.

Lemma 2 (see [4]).

If [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , there exists a unique weak solution [figure omitted; refer to PDF] ; [figure omitted; refer to PDF] to problem (1)-(3), and the expression of the weak solution can be formulated by the following eigenfunction expansion: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the double-parameter Mittag-Leffler function and is defined by [figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the Dirichlet eigenvalues and the orthonormal eigenfunctions of symmetric uniformly elliptic operator [figure omitted; refer to PDF] , respectively.

The following two propositions will be used in the context.

Proposition 3.

[figure omitted; refer to PDF] is a completely monotonic decreasing function for [figure omitted; refer to PDF] and satisfies [figure omitted; refer to PDF]

Proposition 4.

Let [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . Then there exists a constant [figure omitted; refer to PDF] such that [figure omitted; refer to PDF]

3. The Regularized Optimization Problem

In this section, we will propose a regularized optimization method together with its implementations for solving the considered backward problem.

3.1. Regularized Optimization Functional

From results of Lemma 2, formula (10) gives a uniquely weak solution [figure omitted; refer to PDF] for any initial value [figure omitted; refer to PDF] . Naturally, it defines a forward operator [figure omitted; refer to PDF] Clearly, the forward operator [figure omitted; refer to PDF] is a linear map and has the following property.

Lemma 5.

The operator [figure omitted; refer to PDF] is a well-defined bounded linear operator from [figure omitted; refer to PDF] to [figure omitted; refer to PDF] . Moreover, it is injective and compact.

Proof.

From Lemma 2, the solution [figure omitted; refer to PDF] can be represented by [figure omitted; refer to PDF] From Proposition 3, we know that [figure omitted; refer to PDF] for [figure omitted; refer to PDF] ; [figure omitted; refer to PDF] is a very small positive number. By the orthogonality of [figure omitted; refer to PDF] and Proposition 4, we obtain [figure omitted; refer to PDF] which implies that [figure omitted; refer to PDF] . Then by Sobolev embedding theorem, we conclude the compactness of the operator [figure omitted; refer to PDF] .

Results of Lemma 5 show that the backward problem is ill-posed due to the compactness of operator [figure omitted; refer to PDF] . Thus, regularization is necessary for recovering the initial value [figure omitted; refer to PDF] . To this end, we consider a Tikhonov functional as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] is a regularization parameter balancing the fidelity term and the smoothness of the solution. Due to the [figure omitted; refer to PDF] -regularization term [figure omitted; refer to PDF] , the cost functional [figure omitted; refer to PDF] is strongly convex. Subsequently, the unique existence of the minimizer can be obtained by standard arguments.

Theorem 6.

There exists a unique minimizer [figure omitted; refer to PDF] to [figure omitted; refer to PDF] for any given [figure omitted; refer to PDF] .

Now, we formulate the backward problem into the following minimization problem: [figure omitted; refer to PDF]

3.2. Finite Element Method Approximation

Obviously, problem (18) is a function space minimization problem. Here, we use the finite element method to approximate it. Similar to that done in [19, 21], we first triangulate the domain [figure omitted; refer to PDF] with a regular triangulation [figure omitted; refer to PDF] of simplicial elements; let [figure omitted; refer to PDF] be the set of the nodes, and define [figure omitted; refer to PDF] to be the continuous piecewise linear finite element space defined over [figure omitted; refer to PDF] ; that is, [figure omitted; refer to PDF] Then any [figure omitted; refer to PDF] can be repeated as [figure omitted; refer to PDF] , where [figure omitted; refer to PDF] is the value of [figure omitted; refer to PDF] at point [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] is the pyramid function; that is, [figure omitted; refer to PDF]

Next, we need to consider the discretization of the bounded linear operator [figure omitted; refer to PDF] . We will adopt the discrete Galerkin method to solve the direct problem (1)-(3). The time interval [figure omitted; refer to PDF] is partitioned into [figure omitted; refer to PDF] equal subintervals by using nodal points [figure omitted; refer to PDF] , with [figure omitted; refer to PDF] , [figure omitted; refer to PDF] . Then, the time-fractional derivative [figure omitted; refer to PDF] at [figure omitted; refer to PDF] is estimated by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] . Denote by [figure omitted; refer to PDF] the approximation of [figure omitted; refer to PDF] and [figure omitted; refer to PDF]

Now we define the fully discrete finite element method by [figure omitted; refer to PDF] for any [figure omitted; refer to PDF] , where [figure omitted; refer to PDF] . The space [figure omitted; refer to PDF] , in which all functions vanish on the boundary [figure omitted; refer to PDF] , is a subspace of [figure omitted; refer to PDF] . Clearly, (23) is a linear system about [figure omitted; refer to PDF] , [figure omitted; refer to PDF] . Subsequently, there exists a discrete linear operator [figure omitted; refer to PDF] such that [figure omitted; refer to PDF]

Theorem 7.

Let [figure omitted; refer to PDF] and [figure omitted; refer to PDF] be the weak solution of (1)-(3) and the discrete Galerkin finite element solution of (23), respectively. Then there is a constant [figure omitted; refer to PDF] such that, for [figure omitted; refer to PDF] , [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is independent of [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] .

The proof of Theorem 7 follows the same lines as the proof of Theorem 2.1 in [22]. So, we omit it.

3.3. Implementations of the Regularized Optimization Method

Applying the interpolation of finite element, the initial value function [figure omitted; refer to PDF] can be written approximately in the finite element form of [figure omitted; refer to PDF] where [figure omitted; refer to PDF] . Due to the linearity of the homogeneous governing equation and the homogeneous boundary condition, we easily see that problem (1)-(3) satisfies the principle of superposition. Here, we also use this principle of superposition to formulate the continuous problem (18) into the following discrete problem: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , is the finite element solution of [figure omitted; refer to PDF] and satisfies [figure omitted; refer to PDF] for any [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , where [figure omitted; refer to PDF] Therefore, numerical solving of the backward problem is essential to determine the [figure omitted; refer to PDF] -dimensional real vector [figure omitted; refer to PDF] .

From the necessary condition for minimizing the approximation function [figure omitted; refer to PDF] , that is, [figure omitted; refer to PDF] we obtain a linear algebraic system [figure omitted; refer to PDF] Let [figure omitted; refer to PDF] be the solution of (31) for a given regularization parameter [figure omitted; refer to PDF] . Then, we obtain the approximation solution of [figure omitted; refer to PDF] as follows: [figure omitted; refer to PDF]

4. Method for Choosing Regularization Parameters

As we all know, the backward problem for determining the initial value is an ill-posed problem; that is, the round-off errors and the measurement noises may be highly amplified due to the choice of an unreasonable regularization parameter, therefore making the regularization solution completely useless [19, 20]. Because of the important role of regularization parameters, a good strategy for selecting regularization parameters should be taken in the computational process. For a fixed [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , we consider a geometric sequence of regularization parameters [figure omitted; refer to PDF] Then, we employ the discrepancy principle to choose a regularization parameter [figure omitted; refer to PDF] after [figure omitted; refer to PDF] steps with [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the finite element solution with respect to [figure omitted; refer to PDF] and [figure omitted; refer to PDF] .

5. Numerical Examples

In all one-dimensional examples, [figure omitted; refer to PDF] , we divide [figure omitted; refer to PDF] into 100 equal subintervals which means that there are 100 elements and 101 nodes, [figure omitted; refer to PDF] , [figure omitted; refer to PDF] . In all two-dimensional examples, [figure omitted; refer to PDF] , we divide [figure omitted; refer to PDF] into 1024 equal triangle element, [figure omitted; refer to PDF] , [figure omitted; refer to PDF] . In the computational process, the measurement vector [figure omitted; refer to PDF] is obtained actually at the points of the mesh grid and added by randomly distributed perturbations with relative noise level [figure omitted; refer to PDF] ; that is, [figure omitted; refer to PDF] . We take [figure omitted; refer to PDF] and [figure omitted; refer to PDF] in all numerical examples. The relative error of the inverse solutions is defined by [figure omitted; refer to PDF]

Example 1.

We take [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . Let the exact initial value for problem (1)-(3) be [figure omitted; refer to PDF] . Numerical results for relative noise levels 1% and 5% are shown and listed in Figure 1 and Table 1.

Table 1: Some numerical results for Examples 1-4.

Examples	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	RelError
Example 1	0.01	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
0.05	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
Example 2	0.01	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
0.05	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
Example 3	0.01	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
0.05	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
Example 4	0.01	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
0.05	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]

Figure 1: Comparison between exact solution and inverse solution for Example 1.

(a) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

Example 2.

Let the exact initial value for problem (1)-(3) be [figure omitted; refer to PDF] [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . Numerical results with the relative noise levels [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are shown in Figure 2 and listed in Table 1.

Figure 2: Comparison between exact solution and inverse solution for Example 2.

(a) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

Example 3.

Let the exact initial value for problem (1)-(3) be [figure omitted; refer to PDF] . And we take [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . Numerical results are listed in Table 1 and shown in Figures 3 and 4 with relative noise levels 1% and 5%, respectively.

Figure 3: Results with [figure omitted; refer to PDF] for Example 3.

(a) Inverse solution

[figure omitted; refer to PDF]

(b) Error surface

[figure omitted; refer to PDF]

Figure 4: Results with [figure omitted; refer to PDF] for Example 3.

(a) Inverse solution

[figure omitted; refer to PDF]

(b) Error surface

[figure omitted; refer to PDF]

Example 4.

In this example, the exact initial value for problem (1)-(3) is taken as [figure omitted; refer to PDF] . And let [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . Numerical results are listed in Table 1 and shown in Figures 5 and 6 with relative noise levels 1% and 5%, respectively.

Figure 5: Results with [figure omitted; refer to PDF] for Example 4.

(a) Inverse solution

[figure omitted; refer to PDF]

(b) Error surface

[figure omitted; refer to PDF]

Figure 6: Results with [figure omitted; refer to PDF] for Example 4.

(a) Inverse solution

[figure omitted; refer to PDF]

(b) Error surface

[figure omitted; refer to PDF]

Acknowledgments

This work is supported by National Natural Science Foundation of China (11161002), Young Scientists Training Project of Jiangxi Province (no. 20122BCB23024), Natural Science Foundation of Jiangxi Province of China (no. 20142BAB201008), Ground Project of Science and Technology of Jiangxi Universities (no. KJLD14051), and National High-Tech R&D Program of China (2012AA061504).

Conflict of Interests

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

References

[1] I. Podlubny Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Applications , vol. 198, of Mathematics in Science and Engineering, Academic Press, San Diego, Calif, USA, 1999.

[2] R. Metzler, J. Klafter, "Boundary value problems for fractional diffusion equations," Physica A: Statistical Mechanics and its Applications , vol. 278, no. 1-2, pp. 107-125, 2000.

[3] H. E. Roman, P. A. Alemany, "Continuous-time random walks and the fractional diffusion equation," Journal of Physics A Mathematical and General , vol. 27, no. 10, pp. 3407-3410, 1994.

[4] K. Sakamoto, M. Yamamoto, "Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems," Journal of Mathematical Analysis and Applications , vol. 382, no. 1, pp. 426-447, 2011.

[5] B. Jin, R. Lazarov, Z. Zhou, "Error estimates for a semidiscrete finite element method for fractional order parabolic equations," SIAM Journal on Numerical Analysis , vol. 51, no. 1, pp. 445-466, 2013.

[6] S. D. Eidelman, A. N. Kochubei, "Cauchy problem for fractional diffusion equations," Journal of Differential Equations , vol. 199, no. 2, pp. 211-255, 2004.

[7] Y. Lin, C. Xu, "Finite difference/spectral approximations for the time-fractional diffusion equation," Journal of Computational Physics , vol. 225, no. 2, pp. 1533-1552, 2007.

[8] F. Mainardi, "The fundamental solutions for the fractional diffusion-wave equation," Applied Mathematics Letters , vol. 9, no. 6, pp. 23-28, 1996.

[9] X. J. Li, C. J. Xu, "A space-time spectral method for the time fractional diffusion equation," SIAM Journal on Numerical Analysis , vol. 47, no. 3, pp. 2108-2131, 2009.

[10] O. P. Agrawal, "Solution for a fractional diffusion-wave equation defined in a bounded domain," Nonlinear Dynamics , vol. 29, no. 1-4, pp. 145-155, 2002.

[11] J. J. Liu, M. Yamamoto, "A backward problem for the time-fractional diffusion equation," Applicable Analysis , vol. 89, no. 11, pp. 1769-1788, 2010.

[12] X. T. Xiong, J. X. Wang, M. Li, "An optimal method for fractional heat conduction problem backward in time," Applicable Analysis , vol. 91, no. 4, pp. 823-840, 2012.

[13] L. Y. Wang, J. J. Liu, "Data regularization for a backward time-fractional diffusion problem," Computers & Mathematics with Applications , vol. 64, no. 11, pp. 3613-3626, 2012.

[14] Y. Zhang, X. Xu, "Inverse source problem for a fractional diffusion equation," Inverse Problems , vol. 27, pp. 1-12, 2011.

[15] T. Wei, Z. Q. Zhang, "Reconstruction of a time-dependent source term in a time-fractional diffusion equation," Engineering Analysis with Boundary Elements , vol. 37, no. 1, pp. 23-31, 2013.

[16] J. G. Wang, Y. B. Zhou, T. Wei, "Two regularization methods to identify a space-dependent source for the time-fractional diffusion equation," Applied Numerical Mathematics , vol. 68, pp. 39-57, 2013.

[17] Z.-L. Deng, X.-M. Yang, X.L. Feng, "A mollification regularization method for a fractional-diffusion inverse heat conduction problem," Mathematical Problems in Engineering , vol. 2013, 2013.

[18] B. T. Jin, W. Rundell, "An inverse problem for a one-dimensional time-fractional diffusion problem," Inverse Problems , vol. 28, no. 7, pp. 1-19, 2012.

[19] Z. Wang, D. Xu, "On the linear model function method for choosing Tikhonov regularization parameters in linear ill-posed problems," Chinese Journal of Engineering Mathematics , vol. 30, no. 3, pp. 451-466, 2013.

[20] Z. W. Wang, J. J. Liu, "New model function methods for determining regularization parameters in linear inverse problems," Applied Numerical Mathematics , vol. 59, no. 10, pp. 2489-2506, 2009.

[21] Q. Chen, J. J. Liu, "Solving an inverse parabolic problem by optimization from final measurement data," Journal of Computational and Applied Mathematics , vol. 193, no. 1, pp. 183-203, 2006.

[22] Y. J. Jiang, J. T. Ma, "High-order finite element methods for time-fractional partial differential equations," Journal of Computational and Applied Mathematics , vol. 235, no. 11, pp. 3285-3290, 2011.