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

Authors: Zhousheng Ruan; Zewen Wang; Wen Zhang
Journal: Mathematical Problems in Engineering
 (2015)

---

## 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.


---

## Body

## 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 u(x,T)∈H2(Ω)∩H01(Ω) 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α be a constant such that 0<α<1. We consider the following time-fractional diffusion equation:(1)∂αu(x,t)∂tα=(Lu)(x,t),(x,t)∈Ω×(0,T)with homogeneous boundary condition(2)ux,t∂Ω=0,t∈[0,T]and initial condition(3)ux,tt=0=φ(x),x∈Ω¯,where Ω is a bounded domain in Rd(d≥1) and L is symmetric uniformly elliptic operator given by(4)Lu=∑i=1d‍∂∂xi∑j=1d‍θi,j∂∂xjux-cxux;that is, there exists a constant v>0, such that v∑id‍ξi2≤∑i,j=1d‍θi,jξiξj,  x∈Ω¯,  and ξ∈Rd. The coefficients satisfy(5)θi,j∈C1(Ω¯),θi,j=θj,i,c(x)∈C(Ω¯),c(x)≥0,∀x∈Ω¯.Here, ∂αu(x,t)/∂tα is the Caputo fractional derivative which is defined by(6)∂αu(x,t)∂tαWW=1Γ1-α∫0t‍t-η-α∂u∂ηdη,0<α<1,where Γ(1-α) is the Gamma function.If the functionφ(x) 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 φ(x) by the additional data g(x) which is the measurement of the exact value u(x,T) and satisfies (7)u(·,T)-g(·)L2(Ω)≤δfor some known error level δ>0. 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 Section2, 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 u(·,t)∈H01(Ω) such that(8)∂αu∂tα,χ+au,χ=0,∀χ∈X,0<t<T,u(0)=φ(x),where X=H01(Ω), ·,· denotes the inner product in L2(Ω), and(9)a(u,χ)=∫Ω‍∑i,j=1d‍θi,j∂u(x)∂xi·χxj-c(x)u·χdx.Definition 1.
A functionu(x,t) is said to be a weak solution of the direct problem (1)–(3) if u∈C([0,T];L2(Ω))∩C((0,T]; H2(Ω)∩H01(Ω)) and the weak form (8) is satisfied.Lemma 2 (see [4]).
If0<α<1, φ∈L2(Ω), there exists a unique weak solution u∈C([0,T];L2(Ω))∩C((0,T];  H2(Ω)∩H01(Ω)) to problem (1)–(3), and the expression of the weak solution can be formulated by the following eigenfunction expansion:(10)u(x,t;φ)=∑k=1∞‍Eα,1(-λktα)(φ,χk)χk,where Eα,γ(z) is the double-parameterMittag-Leffler function and is defined by(11)Eα,γz∶=∑k=0∞‍zkΓαk+γ,z∈C,α>0,γ≥0;{λi}j=1∞  (0<λ1≤λ2≤⋯,i=1,…,∞) and {χi}j=1∞ are the Dirichlet eigenvalues and the orthonormal eigenfunctions of symmetric uniformly elliptic operator -L, respectively.The following two propositions will be used in the context.Proposition 3.
E
α
,
1
(
-
t
) is a completely monotonic decreasing function for t>0 and satisfies (12)1=Eα,1(0)>Eα,1(-λ1tα)>Eα,1(-λ2tα)>⋯>0.Proposition 4.
Let0<α<2, β∈R, and πα/2<μ<min⁡{π,πα}. Then there exists a constant C=C(α,β,μ)>0 such that(13)Eα,βz≤C1+z,μ≤argz≤π,z∈C.


## 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 Lemma2, formula (10) gives a uniquely weak solution u(x,t;φ) for any initial value φ(x)∈L2(Ω). Naturally, it defines a forward operator (14)M:φ(x)⟼u(x,T;φ).Clearly, the forward operator M is a linear map and has the following property.Lemma 5.
The operatorM is a well-defined bounded linear operator from L2(Ω) to L2(Ω). Moreover, it is injective and compact.Proof.
From Lemma2, the solution u(x,t;φ) can be represented by(15)u(x,t;φ)=∑k=1∞‍Eα,1(-λktα)(φ,χk)χk.From Proposition 3, we know that 0≤Eα,1(-t)<1 for t∈[σ,+∞); σ is a very small positive number. By the orthogonality of {χi}j=1∞ and Proposition 4, we obtain(16)Δu(·,t)L2(Ω)2=∑k=1∞‍λk2Eα,12-λktα2φ,χk2≤∑k=1∞‍C2σ2αφ,χk2=C2σ2αφL2(Ω)2,∀t∈[σ,+∞),which implies that uH2(Ω)≤CφL2(Ω). Then bySobolev embedding theorem, we conclude the compactness of the operator M.Results of Lemma5 show that the backward problem is ill-posed due to the compactness of operator M. Thus, regularization is necessary for recovering the initial value φ(x). To this end, we consider a Tikhonov functional as(17)J(φ)=12u(x,T;φ)-gL2(Ω)2+β2φL2(Ω)2,where φ(x)∈L2(Ω) and β is a regularization parameter balancing the fidelity term and the smoothness of the solution. Due to the L2-regularization term β/2φL2(Ω)2, the cost functional J(φ) is strongly convex. Subsequently, the unique existence of the minimizer can be obtained by standard arguments.Theorem 6.
There exists a unique minimizerφ⋆ to J(φ) for any given β>0.Now, we formulate the backward problem into the following minimization problem:(18)min⁡φ∈L2(Ω)J(φ).


### 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 Ω with a regular triangulation Th of simplicial elements; let {pi}i=0N1 be the set of the nodes, and define Vh to be the continuous piecewise linear finite element space defined over Th; that is,(19)Vh={v:v∈C0(Ω),vΔh∈P1(Δh),∀Δh∈Th}.Then any uh∈Vh can be repeated as uh=∑i=0N1‍uiψi, where ui is the value of uh(x) at point pi, and ψi is the pyramid function; that is,(20)ψi(pj)=1,i=j0,i≠j.Next, we need to consider the discretization of the bounded linear operatorM. We will adopt the discreteGalerkin method to solve the direct problem (1)–(3). The time interval [0,T] is partitioned into N2 equal subintervals by using nodal points 0=t0<t1<⋯<tN2-1<tN2=T, with tk=kτ, τ=T/N2. Then, the time-fractional derivative ∂αu(x,t)/∂tα at tk is estimated by(21)∂αux,t∂tαt=tkWW=1Γ(1-α)∫0tk‍tk-η-α∂u(x,η)∂ηdηWW=1Γ(1-α)∑l=1k‍∫tl-1tl‍tk-η-α∂u(x,η)∂ηdηWW≈τ-αΓ2-α∑l=1k‍ux,tl-ux,tl-1WWW×k+1-l1-α-k-l1-αWW=τ-αΓ(2-α)∑l=1k‍ωl(u(x,tk+1-l)-u(x,tk-l)),where ωl=l1-α-(l-1)1-α,  l=1,2,…,k,  k=1,2,…,N2. Denote by uhk∈Vh the approximation of u(·,tk) and(22)Duh,tα,k=τ-αΓ2-α∑l=1k‍ωluhk+1-l-uhk-l.Now we define the fully discrete finite element method by(23)Duh,tα,k,ψ+a(uhk,ψ)=0,u0=φ(x)for any ψ∈Vh∘, where a(uhk,ψ)=∫Ω‍(∑i,j=1d‍θi,j∂uhk/∂xi·ψxj-c(x)uhk·ψ)dx. The space Vh∘, in which all functions vanish on the boundary ∂Ω, is a subspace of Vh. Clearly, (23) is a linear system about uhk,  k=1,2,…,N2. Subsequently, there exists a discrete linear operator Mh such that(24)Mh:φ⟼uk,k=1,2,…,N2.Theorem 7.
Letu and uhk be the weak solution of (1)–(3) and the discreteGalerkin finite element solution of (23), respectively. Then there is a constant C>0 such that, for 0<α<1,(25)u(·,tk)-uhkL2(Ω)≤C(τ2-α+h),k=1,2,…,N2,where C is independent of h,α, and τ.The proof of Theorem7 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φ(x) can be written approximately in the finite element form of(26)φ≈φh=∑i=0N1‍φiψi(x),where φi∶=φ(xi). 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:(27)min⁡φh∈VhJφh=12∫Ω‍∑i=0N1‍φiuh,iN2-gh2dxW+β2∫Ω‍∑i=0N1‍φiψi(x)2dx,where uh,ik, k=1,2,…,N2, is the finite element solution of u(x,t) and satisfies(28)uh,i0=ψi,Duh,iα,k,ψj+a(uh,ik,ψj)=0,for any ψj∈Sh∘ and k=1,2,…,N2, where(29)Duh,iα,k=τ-αΓ(2-α)∑l=1k‍ωl(uh,ik+1-l-uh,ik-l),gh(x)=∑i=0N1‍giψi(x),a(uh,ik,ψj)=∫Ω‍∑i,j=1d‍θi,j∂uh,ik∂xi·ψjxj-c(x)uh,ik·ψjdx.Therefore, numerical solving of the backward problem is essential to determine the (N1+1)-dimensional real vector Φ=[φ0,…,φN1]T.From the necessary condition for minimizing the approximation functionJ(φh), that is,(30)∂J(φh)∂φi=0,i=0,1,…,N1,we obtain a linear algebraic system(31)AΦ=F,(32)A=∫Ω‍uh,iN2(x)uh,jN2(x)+βψiψjdx(N1+1)×(N1+1),Φ=φ0,φ1,…,φN1T,F=∫Ω‍gh(x)uh,iN2(x)dx(N1+1)×1.Let Φ∗=[φ0∗,φ1∗,…,φN1∗]T be the solution of (31) for a given regularization parameter β. Then, we obtain the approximation solution of φ as follows:(33)φh=∑i=0N1‍φi∗ψi(x).


## 3.1. Regularized Optimization Functional

From results of Lemma2, formula (10) gives a uniquely weak solution u(x,t;φ) for any initial value φ(x)∈L2(Ω). Naturally, it defines a forward operator (14)M:φ(x)⟼u(x,T;φ).Clearly, the forward operator M is a linear map and has the following property.Lemma 5.
The operatorM is a well-defined bounded linear operator from L2(Ω) to L2(Ω). Moreover, it is injective and compact.Proof.
From Lemma2, the solution u(x,t;φ) can be represented by(15)u(x,t;φ)=∑k=1∞‍Eα,1(-λktα)(φ,χk)χk.From Proposition 3, we know that 0≤Eα,1(-t)<1 for t∈[σ,+∞); σ is a very small positive number. By the orthogonality of {χi}j=1∞ and Proposition 4, we obtain(16)Δu(·,t)L2(Ω)2=∑k=1∞‍λk2Eα,12-λktα2φ,χk2≤∑k=1∞‍C2σ2αφ,χk2=C2σ2αφL2(Ω)2,∀t∈[σ,+∞),which implies that uH2(Ω)≤CφL2(Ω). Then bySobolev embedding theorem, we conclude the compactness of the operator M.Results of Lemma5 show that the backward problem is ill-posed due to the compactness of operator M. Thus, regularization is necessary for recovering the initial value φ(x). To this end, we consider a Tikhonov functional as(17)J(φ)=12u(x,T;φ)-gL2(Ω)2+β2φL2(Ω)2,where φ(x)∈L2(Ω) and β is a regularization parameter balancing the fidelity term and the smoothness of the solution. Due to the L2-regularization term β/2φL2(Ω)2, the cost functional J(φ) is strongly convex. Subsequently, the unique existence of the minimizer can be obtained by standard arguments.Theorem 6.
There exists a unique minimizerφ⋆ to J(φ) for any given β>0.Now, we formulate the backward problem into the following minimization problem:(18)min⁡φ∈L2(Ω)J(φ).


## 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 Ω with a regular triangulation Th of simplicial elements; let {pi}i=0N1 be the set of the nodes, and define Vh to be the continuous piecewise linear finite element space defined over Th; that is,(19)Vh={v:v∈C0(Ω),vΔh∈P1(Δh),∀Δh∈Th}.Then any uh∈Vh can be repeated as uh=∑i=0N1‍uiψi, where ui is the value of uh(x) at point pi, and ψi is the pyramid function; that is,(20)ψi(pj)=1,i=j0,i≠j.Next, we need to consider the discretization of the bounded linear operatorM. We will adopt the discreteGalerkin method to solve the direct problem (1)–(3). The time interval [0,T] is partitioned into N2 equal subintervals by using nodal points 0=t0<t1<⋯<tN2-1<tN2=T, with tk=kτ, τ=T/N2. Then, the time-fractional derivative ∂αu(x,t)/∂tα at tk is estimated by(21)∂αux,t∂tαt=tkWW=1Γ(1-α)∫0tk‍tk-η-α∂u(x,η)∂ηdηWW=1Γ(1-α)∑l=1k‍∫tl-1tl‍tk-η-α∂u(x,η)∂ηdηWW≈τ-αΓ2-α∑l=1k‍ux,tl-ux,tl-1WWW×k+1-l1-α-k-l1-αWW=τ-αΓ(2-α)∑l=1k‍ωl(u(x,tk+1-l)-u(x,tk-l)),where ωl=l1-α-(l-1)1-α,  l=1,2,…,k,  k=1,2,…,N2. Denote by uhk∈Vh the approximation of u(·,tk) and(22)Duh,tα,k=τ-αΓ2-α∑l=1k‍ωluhk+1-l-uhk-l.Now we define the fully discrete finite element method by(23)Duh,tα,k,ψ+a(uhk,ψ)=0,u0=φ(x)for any ψ∈Vh∘, where a(uhk,ψ)=∫Ω‍(∑i,j=1d‍θi,j∂uhk/∂xi·ψxj-c(x)uhk·ψ)dx. The space Vh∘, in which all functions vanish on the boundary ∂Ω, is a subspace of Vh. Clearly, (23) is a linear system about uhk,  k=1,2,…,N2. Subsequently, there exists a discrete linear operator Mh such that(24)Mh:φ⟼uk,k=1,2,…,N2.Theorem 7.
Letu and uhk be the weak solution of (1)–(3) and the discreteGalerkin finite element solution of (23), respectively. Then there is a constant C>0 such that, for 0<α<1,(25)u(·,tk)-uhkL2(Ω)≤C(τ2-α+h),k=1,2,…,N2,where C is independent of h,α, and τ.The proof of Theorem7 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φ(x) can be written approximately in the finite element form of(26)φ≈φh=∑i=0N1‍φiψi(x),where φi∶=φ(xi). 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:(27)min⁡φh∈VhJφh=12∫Ω‍∑i=0N1‍φiuh,iN2-gh2dxW+β2∫Ω‍∑i=0N1‍φiψi(x)2dx,where uh,ik, k=1,2,…,N2, is the finite element solution of u(x,t) and satisfies(28)uh,i0=ψi,Duh,iα,k,ψj+a(uh,ik,ψj)=0,for any ψj∈Sh∘ and k=1,2,…,N2, where(29)Duh,iα,k=τ-αΓ(2-α)∑l=1k‍ωl(uh,ik+1-l-uh,ik-l),gh(x)=∑i=0N1‍giψi(x),a(uh,ik,ψj)=∫Ω‍∑i,j=1d‍θi,j∂uh,ik∂xi·ψjxj-c(x)uh,ik·ψjdx.Therefore, numerical solving of the backward problem is essential to determine the (N1+1)-dimensional real vector Φ=[φ0,…,φN1]T.From the necessary condition for minimizing the approximation functionJ(φh), that is,(30)∂J(φh)∂φi=0,i=0,1,…,N1,we obtain a linear algebraic system(31)AΦ=F,(32)A=∫Ω‍uh,iN2(x)uh,jN2(x)+βψiψjdx(N1+1)×(N1+1),Φ=φ0,φ1,…,φN1T,F=∫Ω‍gh(x)uh,iN2(x)dx(N1+1)×1.Let Φ∗=[φ0∗,φ1∗,…,φN1∗]T be the solution of (31) for a given regularization parameter β. Then, we obtain the approximation solution of φ as follows:(33)φh=∑i=0N1‍φi∗ψi(x).


## 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 0<r<1 and β0>0, we consider a geometric sequence of regularization parameters(34)βk=β0rk,k∈N.Then, we employ the discrepancy principle to choose a regularization parameter βk∗ after k∗ steps with(35)∫Ω‍uhk∗x,T-gh2dx≤δ2<∫Ω‍uhkx,T-gh2dx,0≤k<k∗,where uhk(x,t) is the finite element solution with respect to φh and βk.


## 5. Numerical Examples

In all one-dimensional examples,Ω=[0,1], we divide Ω into 100 equal subintervals which means that there are 100 elements and 101 nodes, N2=100, (Lu)(x,t)=d/dx(θ(x)d/dxu)-c(x)u. In all two-dimensional examples, Ω=[0,1]×[0,1], we divide Ω into 1024 equal triangle element, N2=80, (Lu)(x,t)=∇(θ(x,y)∇u)-c(x,y)u. In the computational process, the measurement vector gh is obtained actually at the points of the mesh grid and added by randomly distributed perturbations with relative noise level δ^; that is, gh=uh(x,T)∗(1+δ^∗(2∗rand(size(uh(x,T)))-1)). We take β0=1 and r=0.1 in all numerical examples. The relative error of the inverse solutions is defined by (36)RelError=φexact-φinversionL2(Ω)φexactL2(Ω).Example 1.
We takeα=0.9, θ(x)=1, and c(x)=0. Let the exact initial value for problem (1)–(3) be sin(πx). 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 Examples1–4.
Examples
δ
^
β
RelError
Example1
0.01
1
e
-
6
3.05
e
-
3
0.05
1
e
-
5
1.61
e
-
2
Example2
0.01
1
e
-
5
2.64
e
-
2
0.05
1
e
-
4
6.87
e
-
2
Example3
0.01
1
e
-
7
2.29
e
-
2
0.05
1
e
-
6
6.02
e
-
2
Example4
0.01
1
e
-
3
3.7
e
-
2
0.05
1
e
-
2
5.01
e
-
2Figure 1
Comparison between exact solution and inverse solution for Example1.
(a)
δ
^
=
0.01
(b)
δ
^
=
0.05Example 2.
Let the exact initial value for problem (1)–(3) be(37)φ(x)=x,x∈[0,12](1-x),x∈[12,1],θ(x)=x, c(x)=-0.001x, and α=0.6. Numerical results with the relative noise levels 1% and 5% are shown in Figure 2 and listed in Table 1.Figure 2
Comparison between exact solution and inverse solution for Example2.
(a)
δ
^
=
0.01
(b)
δ
^
=
0.05Example 3.
Let the exact initial value for problem (1)–(3) be 10x3y2(1-x)(1-y). And we take α=0.8, θ(x,y)=1, and c(x,y)=0. 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δ^=0.01 for Example 3.
(a)
Inverse solution
(b)
Error surfaceFigure 4
Results withδ^=0.05 for Example 3.
(a)
Inverse solution
(b)
Error surfaceExample 4.
In this example, the exact initial value for problem (1)–(3) is taken as φ(x,y)=sin(2πx)sin(2πy). And let α=0.8, θ(x,y)=0.01(x+y), and c(x,y)=-0.001. 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δ^=0.01 for Example 4.
(a)
Inverse solution
(b)
Error surfaceFigure 6
Results withδ^=0.05 for Example 4.
(a)
Inverse solution
(b)
Error surface



---

*Source: 414727-2015-03-25.xml*