A modified quasi-boundary value method for an inverse source problem of the time-fractional diffusion equation

Accepted Manuscript A modified quasi-boundary value method for an inverse source problem of the time-fractional diffusion equation Ting Wei, Jungang Wang

PII: DOI: Reference:

S0168-9274(13)00172-4 10.1016/j.apnum.2013.12.002 APNUM 2775

To appear in:

Applied Numerical Mathematics

Received date: 13 June 2013 Revised date: 29 November 2013 Accepted date: 23 December 2013

Please cite this article in press as: T. Wei, J. Wang, A modified quasi-boundary value method for an inverse source problem of the time-fractional diffusion equation, Applied Numerical Mathematics (2013), http://dx.doi.org/10.1016/j.apnum.2013.12.002

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A modified quasi-boundary value method for an inverse source problem of the time-fractional diffusion equation Ting Wei ∗ , Jungang Wang School of Mathematics and statistics, Lanzhou University, Gansu 730000, P.R. China.

Abstract In this paper, we consider an inverse source problem for a time-fractional diffusion equation with variable coefficients in a general bounded domain. That is to determine a spacedependent source term in the time-fractional diffusion equation from a noisy final data. Based on a series expression of the solution, we can transform the original inverse problem into a first kind integral equation. The uniqueness and a conditional stability for the spacedependent source term can be obtained. Further, we propose a modified quasi-boundary value regularization method to deal with the inverse source problem and obtain two kinds of convergence rates by using an a priori and an a posteriori regularization parameter choice rule, respectively. Numerical examples in one-dimensional and two-dimensional cases are provided to show the effectiveness of the proposed method. Key words: Inverse source problem; Fractional diffusion equation; Quasi-boundary value method; Convergence analysis; A priori parameter choice; Morozov’s discrepancy principle.

1

Introduction

The fractional diffusion equations can be used to describe the anomalous diffusion phenomena instead of the classical diffusion procedure and have attracted wide attentions in recent years. The time fractional diffusion equation arises when replacing the standard time derivative with a time fractional derivative and can be used to describe superdiffusion and subdiffusion phenomena [2, 21, 34, 38]. The direct problems, i.e., initial value problem and initial boundary value problems for ∗. Corresponding author. Email address: [email protected] ( Ting Wei ).

Preprint submitted to Elsevier

December 27, 2013

the time fractional diffusion equation have been studied extensively in recent years, for instance, on maximum principle [17], on some uniqueness and existence results [16, 33], on numerical solutions by finite element methods [11, 13] and finite difference methods [14, 23, 35, 46], on analytic solutions [19, 20, 41]. However, for some practical problems, the part of boundary data, or initial data, or diffusion coefficients, or source term may not be given and we want to find them by additional measurement data which will yield to some fractional diffusion inverse problems. The early papers on inverse problems were provided by Murio in [22, 24, 25] for solving the sideways fractional heat equations by mollification methods. After that, some works on fractional inverse problems have been published. In [3], Cheng et al. considered an inverse problem for determining the order of fractional derivative and diffusion coefficient in a fractional diffusion equation and gave a uniqueness result. In [15], Liu et al. solved a backward problem for the time-fractional diffusion equation by a quasi-reversibility regularization method. Zheng et al. in [44, 45] solved the Cauchy problems for the time fractional diffusion equations on a strip domain by a Fourier truncation method and a convolution regularization method. Qian in [31] used a modified kernel method to deal with a sideways fractional equation inverse problem. In [4, 26, 33, 40, 42], some inverse source problems were investigated. Furthermore, the nonlinear fractional inverse problems have been considered recently in [12, 32]. To our knowledge, the research for inverse problems of fractional differential equations is a new topic and do not have too many results now. In this paper, we consider the following inverse source problem for a time fractional diffusion equation with variable coefficients in a general bounded domain. Let Ω be a bounded domain in Rd with sufficient smooth boundary ∂Ω. The inverse source problem we considered for the time-fractional diffusion problem is to determine an unknown source term f (x) from the following equations ⎧ α Dt u(x, t) = (Lu)(x, t) + f (x)q(t), ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ u(x, t) = 0, ⎪ ⎪ ⎪ u(x, 0) = 0, ⎪ ⎪ ⎪ ⎩ u(x, T ) = g(x),

x ∈ Ω, t ∈ (0, T ), 0 < α < 1, (1.1a) x ∈ ∂Ω, t ∈ (0, T ), (1.1b) x ∈ Ω, (1.1c) x ∈ Ω, (1.1d)

where Dαt is the Caputo fractional derivative of order α (0 < α ≤ 1) defined by ⎧  t ⎪ uτ (x, τ) 1 ⎪ ⎪ ⎪ dτ, 0 < α < 1, ⎨ α Γ(1 − α) 0 (t − τ)α Dt u(x, t) = ⎪ ⎪ ⎪ ⎪ ⎩ ut (x, t), α = 1; and −L is a symmetric uniformly elliptic operator defined on D(−L) = H 2 (Ω) 2

(1.2) 

H01 (Ω)

given by ⎞ ⎛ d d  ⎟⎟⎟ ∂ ∂ ⎜⎜⎜⎜ ⎜⎜⎝ ai j (x) Lu(x, t) = u(x, t)⎟⎟⎟⎠ + c(x)u(x, t), ∂xi ∂x j i=1

x ∈ Ω,

(1.3)

j=1

in which the coefficients satisfy ai j = a ji , 1 ≤ i, j ≤ d, ai j ∈ C 1 (Ω), d d   2 ξi ≤ ai j (x)ξi ξ j , x ∈ Ω, ξ ∈ Rd , for a constant ν > 0, ν i=1

(1.4) (1.5)

i, j=1

c(x) ≤ 0, x ∈ Ω, c(x) ∈ C(Ω).

(1.6)

We assume the time-dependent source term q(t) is given. The space-dependent source term f (x) is determined from a noisy final data gδ (x) which satisfies gδ (x) − g(x) ≤ δ,

(1.7)

where  ·  denotes the L2 (Ω) norm and δ > 0 is a noise level. If α = 1, the inverse source problem (1.1a)-(1.1d) is a classical ill-posed problem and has been studied in [7, 36]. However for the fractional inverse source problem, to our knowledge, there are very few works, for example, Sakamoto et al. in [33] used the data u(x0 , t)(x0 ∈ Ω) to determine q(t) when f (x) is given where the authors obtained a Lipschitz stability for q(t). Zhang et al. in [43] used the Fourier truncation method and Wang et al. in [39] used the Tikhonov regularization to solve an inverse source problem with q(t) ≡ 1 in (1.1a) for one-dimensional case with special coefficients. In this study, we fucus on a multi-dimensional problem with variable coefficients in a general bounded domain. The quasi-boundary value method, also called non-local boundary value method in [9], is a regularization technique by replacing the final condition or boundary condition by a new approximate condition. This method has been used to solve some inverse problems for parabolic equation [1, 5, 9, 10], hype-parabolic equations [37] and elliptic equations [6, 8]. The standard quasi-boundary value method to deal with the inverse source problem (1.1a)-(1.1d) is to modify the final condition (1.1d) to form an approximate nonlocal problem ⎧ α Dt v(x, t) = (Lv)(x, t) + f (x)q(t), ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ v(x, t) = 0, ⎪ ⎪ ⎪ v(x, 0) = 0, ⎪ ⎪ ⎪ ⎩ v(x, T ) = gδ (x) − μ f (x), 3

x ∈ Ω, t ∈ (0, T ), x ∈ ∂Ω, t ∈ (0, T ), x ∈ Ω, x ∈ Ω,

(1.8a) (1.8b) (1.8c) (1.8d)

where μ plays a role of regularization parameter. It can be proved by the similar method in Section 4 that the best convergence rate for the source term f (x) is O(δ1/2 ) under an a priori choice of regularization parameter and an a priori bound assumption to the exact f (x) which is not best and can be improved. In this study, we propose a modified version of quasi-boundary value method to solve the inverse source problem (1.1a)-(1.1d), i.e. replacing the final condition (1.1d) with a perturbed condition containing the value (L f )(x) as follows ⎧ α (1.9a) Dt v(x, t) = (Lv)(x, t) + f (x)q(t), x ∈ Ω, t ∈ (0, T ), ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ x ∈ ∂Ω, t ∈ (0, T ), (1.9b) ⎨ v(x, t) = 0, ⎪ ⎪ ⎪ v(x, 0) = 0, x ∈ Ω, (1.9c) ⎪ ⎪ ⎪ ⎩ δ x ∈ Ω. (1.9d) v(x, T ) = g (x) + μ(L f )(x), For this modification, we can obtain a better convergence rate O(δ2/3 ) under an a priori choice of regularization parameter and the same a priori bound assumption. Meanwhile by using Morozov’s discrepancy principle for choosing a regularization parameter, we can also obtain a convergence rate which is not available in the standard quasi-boundary value method (1.8). Another motivation to propose this modified quasi-boundary value method is that we can solve the regularized problem (1.9) by finite difference methods or finite element methods without computing the eigenvalues and eigenfunctions of the operate L which are usually required in the Tikhonov regularization and spectral truncation method. The paper is organized as follows. In Section 2, we give some preliminary results. The uniqueness, ill-posedness and a conditional stability for the inverse source problem (1.1a)-(1.1d) are provided in Section 3. In Section 4, we propose a modified quasi-boundary value method. The convergence estimates under an a priori assumption for the exact solution and two regularization parameter choice rules are obtained. Numerical implementation methods for (1.9) are given in Section 5. Finally two numerical examples in one-dimensional and two-dimensional cases respectively are tested in Section 6 for showing that the proposed methods are effective and stable. A brief conclusion is given in Section 7.

2

Preliminaries

Throughout this paper, we use the following definition and lemmas. Definition 2.1 [28] The Mittag-Leffler function is Eα,β (z) =

∞  k=0

zk , Γ(αk + β)

where α > 0 and β ∈ R are arbitrary constants. 4

z ∈ C,

Lemma 2.1 [28] Let λ > 0, then we have d Eα,1 (−λtα ) = −λtα−1 Eα,α (−λtα ), dt

t > 0,

0 < α < 1.

(2.1)

Lemma 2.2 [30] For 0 < α < 1, η > 0, we have 0 < Eα,1 (−η) < 1. Moreover, Eα,1 (−η) is completely monotonic, that is (−1)n

dn Eα,1 (−η) ≥ 0, dηn

∀n ∈ N ∪ {0}.

Lemma 2.3 For 0 < α < 1, η > 0, we have 0 ≤ Eα,α (−η) ≤ Eα,α (−η) is a monotonic decreasing function with η > 0.

1 . Γ(α)

Moreover,

d Eα,1 (−η), then from Lemma 2.2, the Proof: It is easy to prove that Eα,α (−η) = −α dη proof is clear.

Lemma 2.4 [28] For α > 0 and β ∈ R, we have Eα,β (z) = zEα,α+β (z) +

1 , Γ(β)

z ∈ C.

Lemma 2.5 [18] Assume q(t) ∈ C[0, T ], for 0 < α < 1 and λ > 0, we have

Dαt



t

q(τ)(t − τ)α−1 Eα,α (−λ(t − τ)α )dτ 0  t q(τ)(t − τ)α−1 Eα,α (−λ(t − τ)α )dτ, = q(t) − λ

(2.2) 0 ≤ t ≤ T.

0

We provide the following lemmas which will be used in the proofs of convergence estimates. Lemma 2.6 For any λn satisfying λn ≥ λ1 > 0, there exist positive constants C depending on α, T, λ1 such that C 1 ≤ Eα,α+1 (−λn T α ) ≤ α . α T λn T λn (−λ T α )

(2.3)

α,1 n Proof: By Lemma 2.4, we know Eα,α+1 (−λn T α ) = . From Lemma 2.2, λn T α we have 1 − Eα,1 (−λ1 T α ) 1 ≤ Eα,α+1 (−λn T α ) ≤ . λn T α λn T α Denote C = 1 − Eα,1 (−λ1 T α ), then the proof is completed.

1−E

5

Lemma 2.7 For constants p > 0, μ > 0, β > 0, s ≥ λ1 > 0, we have ⎧ p p ⎪ ⎪ μs2− 2 ⎨C1 μ 4 , 0 < p < 4, ≤⎪ F(s) = 2 ⎩C2 μ, μs + β ⎪ p ≥ 4,

(2.4)

where C1 = C1 (p, β) > 0, C2 = C2 (p, β, λ1 ) > 0 are independent of s. Proof: If 0 < p < 4, then we have lim s→0 F(s) = lim s→∞ F(s) = 0, thus F(s) ≤ sup F(s) ≤ F(s0 ), s∈(0,+∞)

where s0 ∈ (0, +∞) such that F (s0 ) = 0. It is easy to prove that s0 = ( β(4−p) ) 2 > 0, pμ thus we have 4−p 4−p ) 4 ( μ1 ) 4 μ ( β(4−p) p p F(s) ≤ =: C1 (p, β)μ 4 . β(4−p) +β p 1

If p ≥ 4, then for s ≥ λ1 > 0, we have F(s) ≤

μ 1 μ 1 ≤ =: C2 (p, β, λ1 )μ. p p −2 + β s2 β λ 2 −2 1

μs2

Lemma 2.8 For constants p > 0, μ > 0, β > 0, s ≥ λ1 > 0, we have ⎧ 2+p p ⎪ ⎪ μs1− 2 ⎨C3 μ 4 , 0 < p < 2, ≤⎪ G(s) = 2 ⎩C4 μ, μs + β ⎪ p ≥ 2,

(2.5)

where C3 = C3 (p, β) > 0, C4 = C4 (p, β, λ1 ) > 0 are independent of s. Proof: If 0 < p < 2, then we have lim s→0 G(s) = lim s→∞ G(s) = 0, thus we know G(s) ≤ sup G(s) ≤ G(s0 ), s∈(0,+∞) β(2−p) ) 2 > 0, where s0 ∈ (0, +∞) such that G (s0 ) = 0. It is easy to prove that s0 = ( (2+p)μ thus we have 1

G(s) ≤ G(s0 ) =

) ( β(2−p) 2+p

2−p 4

β(2−p) (2+p)

( μ1 )

2−p 4

+β

μ

=: C3 (p, β)μ

2+p 4

If p ≥ 2, then we have G(s) ≤

1 μ 1 μ ≤ =: C4 (p, β, λ1 )μ. p p + β s 2 −1 β λ 2 −1 1

μs2

6

.

3

Uniqueness, ill-posedness and a conditional stability for the inverse source problem (1.1a)-(1.1d)

Denote the eigenvalues of the operator −L as λn and the corresponding eigenfunc tions as ϕn (x) ∈ H 2 (Ω) H01 (Ω), then we have Lϕn = −λn ϕn .

(3.1)

Since −L is a symmetric uniformly elliptic operator, counting according to the multiplicities, we can assume [33] 0 < λ1 ≤ λ2 ≤ · · · ≤ λn ≤ · · · ,

lim λn = +∞,

n→∞

(3.2)

and {ϕn (x)} is an orthonormal basis in space L2 (Ω). Define γ

D((−L) ) = {ψ ∈ L (Ω); 2

∞ 

2 λ2γ n |(ψ, ϕn )| < ∞},

(3.3)

n=1

where (·, ·) is the inner product in L2 (Ω), then D((−L)γ ) is a Hilbert space with the norm ∞  2 12 ψD((−L)γ ) = ( λ2γ (3.4) n |(ψ, ϕn )| ) . n=1

From Theorem 2.2 in [33], we know there exists a unique weak solution u ∈  L2 ((0, T ); H 2 (Ω) H01 (Ω)) for the direct problem (1.1a)-(1.1c) if we know a source function f (x)q(t) ∈ L∞ (0, T ; L2 (Ω)). However the inverse source problem is ill-posed, see the following statements. We will provide the uniqueness and a conditional stability in Theorems 3.1 and 3.2, respectively. In this paper, we focus mainly on the regularization method for solving (1.1a)-(1.1d). By the separation of variables and Lemma 2.5, we know the formal solution for (1.1a)-(1.1c) can be expressed by u(x, t) =

∞ 



t

q(τ)(t − τ)α−1 Eα,α (−λn (t − τ)α )dτϕn (x),

(3.5)

q(τ)(T − τ)α−1 Eα,α (−λn (T − τ)α )dτϕn (x).

(3.6)

fn 0

n=1

where fn = ( f, ϕn ). Let t = T , we have g(x) =

∞  n=1

 fn

T

0

7

Denote gn = (g, ϕn ) and Qn (t) =

t 0

q(τ)(t − τ)α−1 Eα,α (−λn (t − τ)α )dτ, then we have

gn = fn Qn (T ).

(3.7)

To find f (x), we just need to solve the following first kind integral equation  (K f )(x) = k(x, ξ) f (ξ)dξ = g(x), x ∈ Ω, (3.8) Ω

where the kernel is k(x, ξ) =

∞ 

Qn (T )ϕn (x)ϕn (ξ).

n=1

Theorem 3.1 If q(t) ∈ C[0, T ] satisfying q(t) ≥ q0 > 0 for all t ∈ [0, T ], then the solution u(x, t), f (x) of problem (1.1a)-(1.1d) is unique. Proof: By Lemma 2.3, we know Eα,α (−λn (t − τ)α ) ≥ 0, for τ ≤ t. By Lemmas 2.1 and 2.4, we have  T Qn (T ) ≥ q0 (T − τ)α−1 Eα,α (−λn (T − τ)α )dτ = q0 T α Eα,α+1 (−λn T α ). (3.9) 0

Further, by Lemma 2.6, we have Qn (T ) > 0, thus from (3.7), we know if g(x) = 0, we have f (x) = 0. Further, by (3.5), we know u(x, t) = 0. This ends the proof. If we take gk (x) =

ϕk (x) √ λk

in (1.1a)-(1.1d), then we know gk  =

The corresponding source terms are fk (x) = By Lemma 2.6, we have  Qk (T ) ≤ qC[0,T ]

T

ϕk (x) √ , Qk (T ) λk

√1 λk

→ 0, as k → ∞.

and we have  fk  =

1√ . Qk (T ) λk

(T − τ)α−1 Eα,α (−λk (T − τ)α )dτ

0

qC[0,T ] , = qC[0,T ] T Eα,α+1 (−λk T ) ≤ λk α

(3.10)

α

√

λk thus  fk  = Q (T1) √λ ≥ qC[0,T → ∞, as k → ∞. Therefore we know the inverse ] k k source problem (1.1a)-(1.1d) is ill-posed.

Remark 3.1 By (3.8), we know the singular values for the operate K are Qn (T ). q ] . From (3.9) and Lemma 2.6, we know Qn (T ) ≥ From (3.10), we have Qn (T ) ≤ λC[0,T n q0 C . λn

2

By λn = O(n d ), we know that the inverse source problem (1.1) is a moderately ill-posed problem. We give a conditional stability in the following theorem. 8

Theorem 3.2 If q(t) ∈ C[0, T ] satisfying q(t) ≥ q0 > 0, t ∈ [0, T ]. Let f (x) ∈ p D((−L) 2 ) satisfy an a priori bound condition  f D((−L) 2p ) ≤ E,

p > 0,

(3.11)

then we have p

2

 f  ≤ C5 E p+2 g p+2 , where C5 = (q0C)

p − p+2

p > 0,

is a constant depending on α, T, p, λ1 , and q0 .

Proof: From (3.7) and the H¨older inequality, we have f = 2

∞ 

fn2

=

∞ 

n=1 ∞ 

n=1

g2n Q2n (T )

4 p+2

2p gn p+2 g n Q2n (T ) n=1 ∞ ∞   p 2 g2n p+2 ( ≤( g2n ) p+2 . ) p+2 n=1 Qn (T ) n=1

=

(3.12)

Applying Lemma 2.6 and (3.9), we obtain ∞  n=1

g2n Qnp+2 (T )

≤ =

∞ 

fn2 p Qn (T ) n=1 ∞  fn2 λnp (q0C)−p n=1

(3.13) = f

2

p

D((−L) 2 )

−p

(q0C) .

Combining (3.12)-(3.13), we get 2p

4

 f 2 ≤ (q0C)− p+2  f  p+2

p D((−L) 2 )

2p

g p+2 .

The proof is completed.

4

A modified quasi-boundary value method and convergence rates

In this section, we propose a modified quasi-boundary value method to solve problem (1.1a)-(1.1d) and give two convergence estimates under an a priori regularization parameter choice rule and an a posteriori regularization parameter choice rule, respectively. 9

Let {uδμ (x, t) , fμδ (x)} be the solution of the following regularized problem ⎧ α δ ⎪ Dt uμ (x, t) = (Luδμ )(x, t) + fμδ (x)q(t), x ∈ Ω, t ∈ (0, T ), ⎪ ⎪ ⎪ ⎪ ⎪ δ ⎪ x ∈ ∂Ω, t ∈ (0, T ), ⎪ ⎨ uμ (x, t) = 0, ⎪ δ ⎪ ⎪ uμ (x, 0) = 0, x ∈ Ω, ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ uδ (x, T ) = gδ (x) + μ(L f δ )(x), x ∈ Ω, μ μ

(4.1a) (4.1b) (4.1c) (4.1d)

where μ > 0 is a regularization parameter. In following subsections, we will prove that fμδ (x) is a regularized solution for the source term f (x) in (1.1a). By the separation of variables and Lemma 2.5, we know uδμ (x, t) has the following form ∞  uδμ (x, t) = ( fμδ )n Qn (t)ϕn (x), (4.2) where

( fμδ )n

=

n=1

( fμδ (x), ϕn (x)).

From (4.1d), we get

( fμδ )n Qn (T ) + μ( fμδ )n λn = gδn ,

where gδn = (gδ , ϕn ). Thus, ( fμδ )n = uδμ (x, t) =

gδn . Qn (T )+μλn

∞  n=1

and fμδ (x) =

gδn Qn (t)ϕn (x), Qn (T ) + μλn

∞  n=1

Denote fμ (x) =

Substituting ( fμδ )n into (4.2), we get

∞  n=1

(4.3)

gδn ϕn (x). Qn (T ) + μλn

(4.4)

gn ϕn (x). Qn (T ) + μλn

(4.5)

In the following, we give two convergence rate estimates for  fμδ (x)− f (x) by using an a priori and an a posteriori choice rule for the regularization parameter.

4.1

Convergence rate estimate under an a priori regularization parameter choice rule

Theorem 4.1 Let q(t) ∈ C[0, T ] satisfying q(t) ≥ q0 > 0, t ∈ [0, T ]. Suppose the a priori condition (3.11) and the noise assumption (1.7) hold, then, 4 (1) If 0 < p < 4 and choose μ = ( Eδ ) p+2 , we have a convergence rate estimate ⎛ ⎞ ⎜ ⎟⎟ 2 p 1 ⎜  fμδ (x) − f (x) ≤ ⎜⎜⎜⎝  + C1 ⎟⎟⎟⎠ E p+2 δ p+2 ; 2 q0 C 10

(2) If p ≥ 4 and choose μ = ( Eδ ) 3 , we have a convergence rate estimate 2

 fμδ (x)

⎛ ⎞ ⎜⎜⎜ 1 ⎟⎟ 1 2 − f (x) ≤ ⎜⎜⎝  + C2 ⎟⎟⎟⎠ E 3 δ 3 , 2 q0 C

where C1 , C2 are positive constants depending on p, α, T , λ1 and q0 . Proof: By the triangle inequality, we know  fμδ (x) − f (x) ≤  fμδ (x) − fμ (x) +  fμ (x) − f (x).

(4.6)

We firstly give an estimate for the first term. From (4.4), (4.5) and (1.7), we have ∞ 

 fμδ (x) − fμ (x)2 =

(

n=1

gδn − gn )2 ≤ δ2 (sup A(n))2 , Qn (T ) + μλn n

where A(n) =

(4.7)

1 . Qn (T ) + μλn

By Lemma 2.6 and (3.9), we get 1 1 . ≤  q0C/λn + μλn 2 q0Cμ

A(n) ≤

Substituting it into (4.7), we obtain  fμδ (x) − fμ (x) ≤

1 δ  √ . 2 q0 C μ

(4.8)

Now we estimate the second term in (4.6). By (4.5), we can deduce that fμ (x) − f (x) = =

∞  n=1

gn gn − )ϕn (x) Qn (T ) + μλn Qn (T )

n=1

gn −μλn ϕn (x). Qn (T ) Qn (T ) + μλn

∞ 

(

(4.9)

Applying the a priori bound condition (3.11), we obtain  fμ (x) − f (x)2 =

∞  n=1 2

(

gn 2 p μλn 1 ) λn ( )2 p Qn (T ) Qn (T ) + μλn λn

≤ E (sup B(n)) , 2

n

11

(4.10)

where

1− p

μλn 2 B(n) = . Qn (T ) + μλn By Lemma 2.6 and 2.7, we have 2− p

μλn 2 B(n) ≤ 2 μλ + q0C ⎧ n p ⎪ ⎪ ⎨C1 (p, q0C)μ 4 , ≤⎪ ⎪ ⎩C2 (p, q0C, λ1 )μ,

0 < p < 4, p ≥ 4.

Substituting the above inequality into (4.10) and using (4.8), we obtain ⎧ p ⎪ ⎪ δ ⎨C1 Eμ 4 , 0 < p < 4, δ  fμ (x) − f (x) ≤  +⎪ ⎩C2 Eμ, p ≥ 4. 2 q0Cμ ⎪ Choose the regularization parameter μ by ⎧ 4 ⎪ ⎪ ⎨( Eδ ) p+2 , 0 < p < 4, μ=⎪ ⎪ ⎩( δ ) 23 , p ≥ 4, E

(4.11)

(4.12)

(4.13)

then, we have  ⎧ p ⎪ 2 ⎪ 1 ⎪ ⎪ √ + C1 E p+2 δ p+2 , 0 < p < 4, ⎪ ⎪ 2 q C ⎨ 0   uδμ (x, 0) − u(x, 0) ≤ ⎪ ⎪ ⎪ 1 2 ⎪ 1 ⎪ √ + C p ≥ 4. ⎪ 2 E 3 δ3 , ⎩ 2

(4.14)

q0 C

The proof is completed.

4.2

Convergence rate estimate under an a posteriori regularization parameter choice rule

In this subsection, we use an a posterior regularization parameter choice, i.e., Morozov’s discrepancy principle to choose the regularization parameter μ in (4.4). Based on the conditional stability estimate in Theorem 3.2, we can obtain a convergence rate for the regularized solution (4.4). Morozov’s discrepancy principle for our case is to find μ such that K fμδ (x) − gδ (x) = τδ,

(4.15)

where τ > 1 is a constant. According to the following lemma, we know there exists a unique solution for (4.15) if gδ  > τδ > 0. 12

Lemma 4.1 Set ρ(μ) = K fμδ (x) − gδ (x). If gδ  > δ > 0, then the following results hold: (a) ρ(μ) is a continuous function; (b) limμ→0 ρ(μ) = 0; (c) limμ→+∞ ρ(μ) = gδ (x); (d) ρ(μ) is a strictly increasing function over (0, ∞). Proof: The proofs are straightforward results by the expression of ρ(μ) = (

∞  n=1

(

μλn 1 )2 (gδn )2 ) 2 . Qn (T ) + μλn

Theorem 4.2 Let q(t) ∈ C[0, T ] satisfying q(t) ≥ q0 > 0, t ∈ [0, T ]. Suppose the a priori condition (3.11) and the noise assumption (1.7) hold, and there exists τ > 1 such that gδ  > τδ > 0. The regularization parameter μ > 0 is chosen by Morozov’s discrepancy principle (4.15). Then (1) If 0 < p < 2, we have a convergence rate estimate ⎛ ⎞ ⎜ p 2 ⎟ 1 q C ⎜ ⎟ 2 p 1 3 ) 2+p ⎟⎟⎟⎠ E p+2 δ p+2 ;  fμδ (x) − f (x) ≤ ⎜⎜⎜⎝C5 (τ + 1) p+2 +  ( 2 q0 C τ − 1 (2)If p ≥ 2, we have a convergence rate estimate  fμδ (x)

⎛ ⎞ ⎜⎜⎜ 1 q C 1 1⎟ ⎟ 1 1 1 4 − f (x) ≤ ⎜⎜⎝C5 (τ + 1) 2 +  ( ) 2 ⎟⎟⎟⎠ E 2 δ 2 , 2 q0C τ − 1

where C3 , C4 , C5 are positive constants depending on p, α, T , λ1 and q0 . Proof: Similar to (4.6), we have  fμδ (x) − f (x) ≤  fμδ (x) − fμ (x) +  fμ (x) − f (x).

(4.16)

We firstly give an estimate for the second term. From (4.9), we know ∞ 

−μλn ϕn (x) Q (T ) + μλn n n=1 ∞ ∞   −μλn −μλn (gn − gδn ) ϕn (x) + gδn ϕn (x). = Qn (T ) + μλn Qn (T ) + μλn n=1 n=1

K( fμ (x) − f (x)) =

gn

(4.17)

Using (1.7) and (4.15), we get K( fμ (x) − f (x)) ≤ δ + τδ = (τ + 1)δ. 13

(4.18)

Applying the a priori bound condition for f (x), we obtain  fμ (x) − f (x)

p D((−L) 2 )

∞  p −μλn gn 1 =( ( λn2 )2 ) 2 Qn (T ) Qn (T ) + μλn n=1 ∞  gn 1 )2 λnp ) 2 ≤ E. ≤( ( Q (T ) n n=1

(4.19)

By Theorem 3.2, we deduce that p

2

p

 fμ (x) − f (x) ≤ C5 (τ + 1) p+2 E p+2 δ p+2 ,

∀p > 0.

(4.20)

Now we give the bound for the first term. Similar to (4.8), we have  fμδ (x) − fμ (x) ≤

δ 1  √ . 2 q0 C μ

(4.21)

From (4.15), there holds τδ = 

∞  n=1

≤

∞  n=1

μλn gδ ϕn (x) Qn (T ) + μλn n ∞  μλn μλn (gδn − gn )ϕn (x) +  gn ϕn (x) Qn (T ) + μλn Qn (T ) + μλn n=1

(4.22)

≤ δ + J. Using again the a priori bound condition for f (x), we obtain J=

∞  p μλn Qn (T ) 1 gn λn2 ϕn (x) p Qn (T ) + μλn λ 2 Qn (T ) n=1

(4.23)

n

≤ E sup C(n).

(4.24)

n

From Lemma 2.6, we have C(n) =

μλn Qn (T ) 1 p Qn (T ) + μλn λ 2 n

≤

μλn λq1n q0 C λn

1− p

q1 μλn 2 . p ≤ q0C + μλ2n + μλn λn2 1

(4.25)

where we denote q1 = qC[0,T ] . By Lemma 2.8, (4.25) becomes ⎧ 2+p ⎪ ⎪ ⎨q1C3 (p, q0C)μ 4 , 0 < p < 2, C(n) ≤ ⎪ ⎪ ⎩q1C4 (p, q0C, λ1 )μ, p ≥ 2.

14

(4.26)

Combining (4.22)-(4.26), we obtain ⎧ 2+p ⎪ ⎪ ⎨q1C3 Eμ 4 , 0 < p < 2, (τ − 1)δ ≤ ⎪ ⎪ ⎩q1C4 Eμ, p ≥ 2.

(4.27)

This yields ⎧ qC 4 4 1 3 2+p E 2+p ⎪( τ−1 ) ( δ ) , 0 < p < 2, 1 ⎪ ⎨ ≤⎪ ⎪ q1C4 E , μ ⎩ p ≥ 2. τ−1 δ

(4.28)

Substituting (4.28) into (4.21), we get p

2

p

 f δ (x) − f (x) ≤ C5 (τ + 1) p+2 E p+2 δ p+2 ⎧μ p 2 2 1 C 3 2+p ⎪ ⎪ √1 ( qτ−1 ) E 2+p δ 2+p , 0 < p < 2, ⎪ ⎪ ⎨ 2 q0C +⎪ q1 C4 1 1 1 ⎪ 1 ⎪ p ≥ 2. ⎪ ⎩ √ ( τ−1 ) 2 E 2 δ 2 , 2

(4.29)

q0 C

This ends the proof. Remark 4.1 From Theorems 4.1 and 4.2, we can see that the best convergence 2 rate for the a priori parameter choice method is O(δ 3 ) and for the a posteriori 1 parameter choice method is O(δ 2 ). Both of them have the saturation phenomenon. However, because the a priori bound E is difficult to obtain in the a priori parameter choice method, we prefer to use the a posteriori choice rule to get the regularization parameter.

Remark 4.2 The proposed method can be extended to a more general case, that is, the final value condition (4.1d) can be replaced by uδμ (x, T ) = gδ (x)+μ(L2m−1 fμδ )(x), where m = 1, 2, 3, · · · . Then the best convergence rate could be improved. For an a priori parameter choice rule, when p ≥ 4m, we can get the best convergence rate 2m as O(δ 2m+1 ). For the Morozov discrepancy principle, when p ≥ 4m − 2, the best 2m−1 convergence rate is O(δ 2m ).

5

Numerical implementations

Since the analytic solution of problem (1.1a) is difficult to obtain, we construct the final data g(x) by solving the following forward problem ⎧ α ⎪ Dt u(x, t) = (Lu)(x, t) + f (x)q(t), x ∈ Ω, t ∈ (0, T ), 0 < α < 1 ⎪ ⎪ ⎪ ⎪ ⎨ u(x, t) = 0, x ∈ ∂Ω, t ∈ (0, T ), ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ u(x, 0) = 0, x ∈ Ω. 15

(5.1)

with the given data f (x), q(t). Numerical schemes for the forward problem (5.1) and the regularized problem (4.1a)-(4.1d) are sketched as follows. 5.1 One-dimensional case Without loss of generality, we assume Ω = (0, 1). Take the grid sizes for time and space variable in the finite difference algorithm are Δt = NT and Δx = M1 , respectively. The grid points on the time interval [0, T ] are labeled tn = nΔt, n = 0, 1, · · · , N, and the grid points in the space interval [0, 1] are xi = iΔx, i = 0, 1, 2, · · · , M. The approximate values of function u at the grid points are denoted uni ≈ u(xi , tn ). The time-fractional derivative is approximated by (Δt)−α  j+1 j ≈ w j (un− − un− i i ), Γ(2 − α) j=1 n

Dαt u(xi , tn )

(5.2)

for i = 1, · · · , M − 1, n = 1, · · · , N where w j = j1−α − ( j − 1)1−α . This scheme was used in [27, 46]. We use the scheme Lu(xi , tn ) ≈

1 (a 1 un − (ai+ 12 + ai− 12 )uni + ai− 12 uni−1 ) + c(xi )uni (Δx)2 i+ 2 i+1

for i = 1, · · · , M −1, n = 1, · · · , N to approach the value of Lu at point (xi , tn ) where ai+ 12 = a(xi+ 12 ) with xi+ 12 = (xi + xi+1 )/2. By the initial condition and boundary condition in (5.1), it is easy to get a numerical solution for forward problem (5.1) from the finite difference scheme (Δt)−α  j+1 j w j (un− − un− i i ) Γ(2 − α) j=1 n

=

(5.3)

1 (a 1 un − (ai+ 12 + ai− 12 )uni + ai− 12 uni−1 ) + c(xi )uni + f (xi )q(tn ). (Δx)2 i+ 2 i+1

Denote U n = (un1 , un2 , · · · , unM−1 )T , F = ( f (x1 ), f (x2 ), · · · , f (x M−1 ))T , then the scheme (5.3) leads to the following iteration scheme AU 1 = q(t1 )F, AU n = σ(β1 U n−1 + β2 U n−2 + · · · + βn−1 U 1 ) + q(tn )F, 16

(5.4) (5.5)

−α

(Δt) where σ = Γ(2−α) , β j = w j − w j+1 and A is a tridiagonal matrix given by Aii = σ + 1 1 r(ai+ 2 + ai− 2 ) − c(xi ) for i = 1, 2 · · · , M − 1 and Ai,i−1 = −rai− 12 for i = 2, 3 · · · , M − 1 1 and Ai,i+1 = −rai+ 12 for i = 1, 2, · · · , M − 2, where r = (Δx) 2.

We take g = U N as the ”exact” final data. For the regularized problem (4.1a)-(4.1d), we discretize equation (4.1a) by the finite difference scheme mentioned above, then get AV 1 = q(t1 )Fμδ , AV = σβ1 V n

n−1

(5.6) + σβ2 V

n−2

+ · · · + σβn−1 V + 1

q(tn )Fμδ ,

n = 2, 3, · · · , N(5.7)

where V n = (vn1 , vn2 , · · · , vnM−1 )T with vni is the approximate value of uδμ (xi , tn ) and δ )T with fiδ is the approximate value of fμδ (xi ). Furthermore, Fμδ = ( f1δ , f2δ , · · · , f M−1 we can deduce the following relationship V N = BFμδ . The condition (4.1d) is approximated by   1 N δ δ δ δ δ 1 1 1 1 vi = g (xi ) + μ (a f − (ai+ 2 + ai− 2 ) fi + ai− 2 fi−1 ) + c(xi ) fi (Δx)2 i+ 2 i+1

(5.8)

(5.9)

for i = 1, · · · , M − 1. The matrix form is V N = Gδ − μHFμδ ,

(5.10)

where Gδ = (gδ (x1 ), · · · , gδ (x M−1 ))T and H is a tridiagonal matrix given by H ii = r(ai+ 12 +ai− 12 )−c(xi ) for i = 1, 2 · · · , M −1 and H i,i−1 = −rai− 12 for i = 2, 3 · · · , M −1 and H i,i+1 = −rai+ 12 for i = 1, 2, · · · , M − 2. . Combining (5.8) and (5.10), we know the unknown source vector satisfies (B + μH)Fμδ = Gδ .

(5.11)

5.2 Two-dimensional case For the forward problem, we compute firstly some eigenfunctions of operator −L and then compute the Fourier coefficients of f (x), by using the formula (3.6) with a suitable truncation to get the function g(x) as the ”exact” final data. In our computations, we use the matlab command pdeeig to generate eigenfunctions. To solve the regularized problem (4.1a)-(4.1d) in two dimensional case, we use the finite difference scheme (5.2) to approach the time fractional derivative and then 17

employ a finite element method to discretize the resulted elliptic problem at each time step. By the nonlocal boundary condition (4.1d), we can finally deduce a linear system of equations for the space-dependent source term on mesh nodes. The detail is given in the following. Denote vn (x) = uδμ (x, tn ). By the finite difference scheme (5.2), we can deduce the following Dirichlet problems for elliptic equations ⎧ 1 1 δ ⎪ ⎪ ⎨ σv (x) − Lv (x) = fμ (x)q(t1 ), ⎪ ⎪ ⎩ v1 (x) = 0, x ∈ ∂Ω,

x ∈ Ω,

(5.12)

and ⎧ n−1  ⎪ ⎪ ⎪ n n ⎪ β j vn− j (x) + fμδ (x)q(tn ), x ∈ Ω, ⎪ ⎨ σv (x) − Lv (x) = σ ⎪ j=1 ⎪ ⎪ ⎪ ⎪ ⎩ vn (x) = 0, x ∈ ∂Ω,

(5.13)

for n = 2, · · · , N. Let x j , j = 1, 2, · · · , m be the mesh nodes located in Ω and ψ j be the corresponding finite element basis functions. Denote V n = (uδμ (x1 , tn ), · · · , uδμ (xm , tn ))T , Fμδ = ( fμδ (x1 ), · · · , fμδ (xm ))T , then by the standard finite element procedure, we can deduce the following linear equations for (5.12) and (5.13) as (σM + S + D)V 1 = q(t1 )MFμδ , and (σM + S + D)V n = σ 2

n−1 

β j MV n− j + q(tn )MFμδ ,

(5.14)

(5.15)

j=1

where M = ((ψi , ψ j ))m×m , S = ( k,l=1 (akl ∂ xk ψi , ∂ xl ψ j ))m×m , D = ((−c(x)ψi , ψ j ))m×m are mass matrix and stiff matrixes in which (·, ·) is the L2 inner product. By a recursion process, we can find a relationship between Fμδ and V N with V N = BFμδ ,

(5.16)

where B is a matrix. The nonlocal boundary condition (4.1d) in finite element scheme becomes MV N + μ(S + D)Fμδ = MGδ , (5.17) where Gδ = (gδ (x1 ), · · · , gδ (xm ))T . Substituting (5.16) into (5.17), we can get (MB + μ(S + D))Fμδ = MGδ . 18

(5.18)

Solving the above linear system of equations, we finally obtain the regularized so lution fμδ (x) = mj=1 fμδ (x j )ψ j (x).

6

Numerical experiments

In this section, we present some numerical results for two examples in one-dimensional and two-dimensional cases to show the effectiveness and stability of our proposed method. The noisy data is generated by adding a random perturbation, i.e., gδ (xi ) = g(xi ) + εg(xi ) · (2rand(i) − 1). The corresponding noise level is calculated by δ = gδ (x) − g(x). To show the accuracy of numerical solution, we compute the approximate L2 error denoted by e( f, ε) =  f (x) − fμδ (x), (6.1) 2 and the approximate relative error in L norm denoted by er ( f, ε) =  f (x) − fμδ (x)/ f (x).

(6.2)

To verify the convergence rate, we use the following definition Convergence Order = log2

e( f, 2ε) . e( f, ε)

(6.3)

In our numerical experiments, we always fix T = 1. For computing Mittag-Leffler function, we need a better algorithm, e.g. see [29]. In general, the a priori bound E is difficult to obtain, thus we only give the numerical results by using the a posteriori parameter choice rule which is independent of E. The regularization parameter is chosen by (4.15) with τ = 1.1. Here we use the matlab command fzero to find μ. Example 1. Let d = 1, Ω = (0, 1), a(x) = x2 + 1, c(x) = −(x + 1). Take a source function q(t) = e−t and f (x) = [x(1 − x)]α sin 5πx. In order to avoid ‘inverse crime’, we use a finer grid to compute the forward problem, i.e. take M = 100, N = 200 and choose M = 50, N = 100 for solving the regularized inverse problem. The numerical results by using the a posteriori parameter choice rule for various noise levels ε = 0.005, 0.01, 0.05 in the case of α = 0.2, 0.8 are showed in Figure 1 19

0.8

0.4 exact ε=0.005 ε=0.01 ε=0.05

0.2

δ

f(x) and its approximation fμ(x)

0.4 0.2 0 −0.2 −0.4

0.1

0

−0.1

−0.2

−0.6 −0.8

exact ε=0.005 ε=0.01 ε=0.05

0.3

δ

f(x) and its approximation fμ(x)

0.6

0

0.2

0.4

0.6

0.8

−0.3

1

0

0.2

0.4

x

0.6

0.8

1

x

(a) α = 0.2

(b) α = 0.8

Figure 1. The exact and regularized source terms given by the a posteriori parameter choice rule for Example 1.

in which we use μ = 1.08e − 8, 2.06e − 8, 9.31e − 8 for α = 0.2 and μ = 1.4e − 9, 2.69e − 8, 1.25e − 7 for α = 0.8. We can see that the numerical results are in good agreement with the exact shape. In Table 1, we show the numerical errors of Example 1 for different α with a fixed ε = 0.01. It can be seen that the numerical accuracy is very stable to the fractional order α. Table 1 Numerical results of Example 1 for different α with ε = 0.01. α

0.05

0.1

0.3

0.7

0.9

0.95

e( f, 0.01)

0.0198

0.0180

0.0132

0.0079

0.0060

0.0055

er ( f, 0.01)

0.0308

0.0306

0.0323

0.0377

0.0393

0.0396

The numerical errors and convergence orders for Example 1 (α = 0.6) with different δ are shown in Table 2, from which we can see that the numerical error is decreasing as the level of noise becomes smaller and the convergence order is near about 0.5. This is consistent with our convergence estimate. Table 2 Numerical results of Example 1 for different ε with α = 0.6. ε

0.0005

0.001

0.002

0.004

0.008

0.016

0.032

0.064

e( f, ε)

0.0024

0.0031

0.0043

0.0060

0.0082

0.0110

0.0149

0.0207

er ( f, ε)

0.0096

0.0128

0.0176

0.0243

0.0332

0.0448

0.0606

0.0841

0.4124

0.4651

0.4645

0.4477

0.4333

0.4369

0.4731

Order

Example 2. Consider the case of d = 2 and denote the coordinate as (x, y). Let Ω 20

be a bounded domain with a smooth boundary given by parametrization ∂Ω = {(x, y) = r(θ)(cos θ, sin θ), θ ∈ [0, 2π]},

(6.4)

where r(θ) = 1 + 0.1 cos 5θ, see Figure 2 for the configuration of Ω. We take a11 (x, y) = x2 +3, a12 (x, y) = a21 (x, y) = x+y+1, a22 (x, y) = y2 +3, c(x, y) = −(x + y)2 . The exact time-dependent source term is q(t) = e−t + tα + 1 and the spacedependent source term is α f (x, y) = exp( √ ), 2 x +y2 ( r(θ) )2 − 1 where θ is the polar angle of (x, y). For (x, y) = (0, 0), we define f (0, 0) = exp(−α). In solving the forward problem, we use the matlab command pdeeig and spend 40 minutes to compute the first 338 eigenvalues in interval [0, 6000] and the corresponding eigenfunctions. To obtain the regularized solution fμδ (x, y), along the time direction, we take N = 20. The mesh in domain Ω is generated by PDE toolbox in matlab and consists of 2205 nodes and 4256 triangles. Domain Ω and its mesh 1.5

1

y

0.5

0

−0.5

−1

−1.5 −1.5

−1

−0.5

0 x

0.5

1

1.5

(a) Domain Figure 2. Domain and mesh

In Figure 3, we show the exact source term and the regularized approximations given by the a posteriori parameter choice rule with ε = 0.01 in the case of α = 0.2, 0.8 in which we use μ = 2.5304e−5 for α = 0.2, μ = 4.6147e−6 for α = 0.8. It can be observed that the proposed method gives accurate numerical reconstructions and is also effective for two dimensional example. In Table 3, we show the numerical results of Example 2 for various α with ε = 0.005 from which we can see the numerical results depend very little on α. The numerical results and convergence orders of Example 2 (α = 0.6) for different δ are shown in Table 4. We can also see that the numerical errors become smaller if decreasing amounts of random noise. The convergence order is about 0.2 which 21

0.5 0.4

0.6

0.3 f(x,y)

0.4 0.2

0.2 0.1

0 2

0 2 1

1

1.5

−1

−1

x

(b) Exact f (x, y) for α = 0.8

1

0.5

0.8

0.4

0.6

0.3

μ

fδ (x,y)

0.4 0.2

0.2 0.1

0 2

0 2 1

1.5

1

0.02

0.0 0.0 0 4 1 .0 3

.0

1

03 0.

102 0.0 0.0

02

0.03 0.04

2

y

4

0.03

02

0.01

−0

0.02 0.0 2

0.02

2

.0

0.04

0.0

−−0 0 0.1 0..106.02 −0−.0 8 −0.0 0.0 60 2

0

0−.000.02 20.04

−0.

0.01 02

(e) Numerical error for α = 0.2

−0.5

0.03 0.02 0.010 −0.01

0.05

−1

0.021 00.0 0.03 2

−1

.0

0 0.01 0.02 0.03 0.04 0.03 0.04 0.01 .05 00.0 2

0.03

−0

2

−0.8

.03 02 0

0.

0.

0.01

1 .0 −−0 0.02

−0.4

1

1 0.0 02

4

0.0

−0.6

0.

1

−0.2

0 1

0.0 0

0.0

0.5

0

0.0 0.0 0 3 00.0.1032 −0 .01

0

5

0 x

0.0

−0.5

.0

−1

−0

0 0.02 61 .00. −0− 2 .0 −0

−1

.1 −082 .0 −0.0 −002 0.

04 0. 4

8

−0.8

06

−0.6

0.

4 0.0 − 0.04 02 0.0 −0.0 0.02 0. 4 −− 0.00.0 −204.08 6 −0.0−40.0 0 −0.0 2 0 6−0.08 .1 1 0.02 0.02 −− −0 −.0068. 00 .0.0 0 0 .0 204.0 −−0 2

−0.4

0.2

0.0

.10.0 −0−

0.06

0

0.4

0.03 0. 1 3 02 0.0 00. .002 0.01 0 0.02 5 4 0 0.0.0 −0 0. .01 05

0.03 0.01 01 −0.0

1

0

0.06

x

0.0

2

.1

04

0.

0.06

.01

0.6

−0

.0 −0 8 −−00.0 .1 .0024 00.0.02 −0−.00 2 4.0 4 0 6 0.06

0.06

02 −0.0 −0.04 −0.008.1 .1 − −00.06 0−

−0.2

0.04

0

0.02

6 .024 −0−.0−00.0

0.02

0.8

02 −0. −0.03

0.0 02 −0 − .0 −04 −.006 .08 .02 0240 −0− .1 0.06 00..0 −0.04 −0.0 0 2 0.02 0.04 0 .06 −0

−0.02 0.−004.08 −0 6.04 −00−.0 .002 −0

−1

(d) Numerical fμδ (x, y) for α = 0.8

−0

1

−0.5

−2

y

x

0.03

−1

(c) Numerical fμδ (x, y) for α = 0.2

0.2

0.5 0

−1

−0.5

−2

y

0.6

1

0

0.5 0

−1

0.8

1.5

1

0

0.

1

0.01

μ

−0.5

−2

y

x

(a) Exact f (x, y) for α = 0.2

fδ (x,y)

0.5 0

−1

−0.5

−2

y

y

1

0

0.5 0

−1

0.4

1.5

1

0

0.0

f(x,y)

1 0.8

0 x

0.5

1

(f) Numerical error for α = 0.8

Figure 3. The exact and regularized source terms given by the a posteriori parameter choice rule for Example 2. Table 3 Numerical results for Example 2 with various α and a fixed ε = 0.005. α

0.05

0.1

0.3

0.7

0.9

0.95

e( f, 0.005)

0.2070

0.1340

0.0603

0.0403

0.0325

0.0307

er ( f, 0.005)

0.1374

0.0997

0.0646

0.0784

0.0829

0.0835

22

is somewhat lower than the one in Example 1. One possible reason is that Example 2 is not smooth enough as Example 1. Table 4 Numerical results for Example 2 with α = 0.6. ε

0.0005

0.001

0.002

0.004

0.008

0.016

0.032

0.064

e( f, ε)

0.0326

0.0354

0.0390

0.0430

0.0473

0.0525

0.0597

0.0715

er ( f, ε)

0.0552

0.0598

0.0659

0.0727

0.0800

0.0887

0.1009

0.1208

0.1161

0.1401

0.1414

0.1382

0.1489

0.1859

0.2607

Order

7

Conclusion

In this paper, we investigate an inverse source problem for a time-fractional diffusion equation with variable coefficients defined in a general domain. The conditional stability is given. We propose a modified quasi-boundary value method for obtaining a regularized solution. Based on an a priori assumption for the exact solution, the error estimates are obtained under an a priori regularization parameter choice rule and Morozov’s discrepancy principle, respectively. Numerical examples show that our proposed method is effective and stable.

Acknowledgments The authors would like to thank the editor and the anonymous referees for their valuable comments and helpful suggestions that improve the quality of our paper. This work is supported by the NSF of China (11371181, 11171136) and the Fundamental Research Funds for the Central Universities (lzujbky-2013-k02).

References [1] K.A. Ames and L.E. Payne. Asymptotic behavior for two regularizations of the Cauchy problem for the backward heat equation. Math. Models Methods Appl. Sci., 8(1):187, 1998. [2] B. Berkowitz, H. Scher, and S.E. Silliman. Anomalous transport in laboratory-scale, heterogeneous porous media. WATER RESOUR. RES. , 36(1):149–158, JAN 2000. 23

[3] J. Cheng, J. Nakagawa, M. Yamamoto, and T. Yamazaki. Uniqueness in an inverse problem for a one-dimensional fractional diffusion equation. Inverse Problems, 25(11):115002, 16, 2009. [4] G. Chi, G. Li, and X. Jia. Numerical inversions of a source term in the fade with a dirichlet boundary condition using final observations. Comput. Math. Appl., 62(4):1619–1626, 2011. [5] M. Denche and K. Bessila. A modified quasi-boundary value method for ill-posed problems. J. Math. Anal. Appl., 301(2):419–426, 2005. [6] X.L. Feng, L. Eld´en, and C.L. Fu. A quasi-boundary-value method for the cauchy problem for elliptic equations with nonhomogeneous neumann data. J. Inverse Ill-Posed Probl., 18(6):617–645, 2010. [7] F.Z. Geng and Y.Z. Lin. Application of the variational iteration method to inverse heat source problems. Comput. Math. Appl., 58(11-12):2098–2102, 2009. [8] D.N. H`ao, N.V. Duc, and D. Lesnic. A non-local boundary value problem method for the Cauchy problem for elliptic equations. Inverse Problems, 25:055002, 2009. [9] D.N. H`ao, N.V. Duc, and D. Lesnic. Regularization of parabolic equations backward in time by a non-local boundary value problem method. IMA J. Appl. Math., 75(2):291–315, 2010. [10] D.N. H`ao, N.V. Duc, and H. Sahli. A non-local boundary value problem method for parabolic equations backward in time. J. Math. Anal. Appl., 345(2):805–815, 2008. [11] Y. J. Jiang and J. T. Ma. High-order finite element methods for time-fractional partial differential equations. J. Comput. Appl. Math., 235(11):3285–3290, APR 1 2011. [12] B.T. Jin and W. Rundell. An inverse problem for a one-dimensional timefractional diffusion problem. Inverse Problems, 28(7), JUL 2012. [13] B.T. Jin,R. Lazarov and Z. Zhou. Error estimates for a semidiscrete finite element method for fractional order parabolic equations. SIAM J. Numer. Anal., 51(1): 445-466. 2013. [14] Y.M. Lin and C.J. Xu. Finite difference/spectral approximations for the timefractional diffusion equation. J. Comput. Phys., 225(2):1533–1552, 2007. [15] J.J. Liu and M. Yamamoto. A backward problem for the time-fractional diffusion equation. Appl. Anal., 89(11):1769–1788, 2010. [16] Y. Luchko. Some uniqueness and existence results for the initial-boundaryvalue problems for the generalized time-fractional diffusion equation. Comput. Math. Appl., 59(5):1766–1772, 2010. [17] Y. Luchko. Maximum principle and its application for the time-fractional diffusion equations. Fract. Calc. Appl. Anal. , 14(1):110–124, 2011. [18] A.A. Kilbas, H.M. Srivastava, and J.J. Trujillo. Theory and applications of fractional differential equations, volume 204 of North-Holland Mathematics Studies. Elsevier Science B.V., Amsterdam, 2006. [19] F. Mainardi, G. Pagnini, and R. Gorenflo. Some aspects of fractional diffusion equations of single and distributed order. Appl. Math. Comput., 187:295–305, 24

2007. [20] R. Metzler and J. Klafter. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Phys. Rep., 339:1–77, 2000. [21] R. Metzler and J. Klafter. Subdiffusive transport close to thermal equilibrium: From the Langevin equation to fractional diffusion. Phys. Rev. E, 61(6, Part a):6308–6311, JUN 2000. [22] D. A. Murio. Stable numerical solution of a fractional-diffusion inverse heat conduction problem. Comput. Math. Appl., 53:1492–1501, 2007. [23] D. A. Murio. Implicit finite difference approximation for time fractional diffusion equations. Comput. Math. Appl., 56:1138–1145, 2008. [24] D. A. Murio. Time fractional IHCP with Caputo fractional derivatives. Comput. Math. Appl., 56:2371–2381, 2008. [25] D. A. Murio. Stable numerical evaluation of Gr¨unwald-Letnikov fractional derivatives applied to a fractional ihcp. Inverse Probl. Sci. Eng., 17:229–243, 2009. [26] D. A. Murio and C. E. Mej´ıa. Source terms identification for time fractional diffusion equation. Rev. Colombiana Mat., 42(1):25–46, 2008. [27] D.A. Murio. Implicit finite difference approximation for time fractional diffusion equations. Comput. Math. Appl., 56(4):1138–1145, 2008. [28] I. Podlubny. Fractional differential equations, volume 198 of Mathematics in Science and Engineering. Academic Press Inc., San Diego, CA, 1999. [29] I. Podlubny and M. Kacenak. Mittag-leffler function. The MATLAB routine, http://www. mathworks. com/matlabcentral/fileexchange, 2006. [30] H. Pollard. The completely monotonic character of the mittag-leffler function Eα (−x). Bull. Amer. Math. Soc., 54(12):1115–1116, 1948. [31] Z. Qian. Optimal modified method for a fractional-diffusion inverse heat conduction problem. Inverse Probl. Sci. Eng., 18(4):521–533, 2010. [32] W. Rundell, X. Xu, and L. H. Zuo. The determination of an unknown boundary condition in a fractional diffusion equation. Appl. Anal., 1–16, 2012. [33] K. Sakamoto and M. Yamamoto. Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems. J. Math. Anal. Appl., 382(1):426–447, OCT 1 2011. [34] E. Scalas, R. Gorenflo, and F. Mainardi. Fractional calculus and continuoustime finance. Phys. A, 284:376–384, 2000. [35] R. Scherer, S. L. Kalla, L. Boyadjiev, and B. Al-Saqabi. Numerical treatment of fractional heat equations. Appl. Numer. Math., 58:1212–1223, 2008. [36] A. Shidfar, A. Babaei, and A. Molabahrami. Solving the inverse problem of identifying an unknown source term in a parabolic equation. Comput. Math. Appl., 60(5):1209–1213, 2010. [37] R.E. Showalter. Cauchy problem for hyper-parabolic partial differential equations. North-Holland Mathematics Studies, 110:421–425, 1985. [38] I. M. Sokolov and J. Klafter. From diffusion to anomalous diffusion: A century after Einstein’s Brownian motion. Chaos, 15:1–7, 2005. [39] J.G. Wang, Y.B. Zhou, and T. Wei. Two regularization methods to identify a space-dependent source for the time-fractional diffusion equation. Appl. 25

Numer. Math., 68:39–57, 2013. [40] H. Wei, W. Chen, H. G. Sun, and X. C. Li. A coupled method for inverse source problem of spatial fractional anomalous diffusion equations. Inverse Probl. Sci. Eng., 18(7):945–956, 2010. [41] W. Wyss. The fractional diffusion equation. J. Math. Phys., 27:2782–2785, 1986. [42] Y. Zhang and X. Xu. Inverse source problem for a fractional diffusion equation. Inverse Problems, 27(3):035010, 12, 2011. [43] Z.Q. Zhang and T. Wei. Identifying an unknown source in time-fractional diffusion equation by a truncation method. Appl. Math. Comput., 219(11):5972– 5983, 2013. [44] G. H. Zheng and T. Wei. Spectral regularization method for a Cauchy problem of the time fractional advection-dispersion equation. J. Comput. Appl. Math., 233(10):2631–2640, 2010. [45] G. H. Zheng and T. Wei. A new regularization method for a Cauchy problem of the time fractional diffusion equation. Adv. Comput. Math., 36(2):377–398, 2012. [46] P. Zhuang and F. Liu. Implicit difference approximation for the time fractional diffusion equation. J. Appl. Math. Comput., 22(3):87–99, 2006.

26