Skip to content
Publicly Available Published by De Gruyter March 22, 2016

Error estimates for approximations of distributed order time fractional diffusion with nonsmooth data

  • Bangti Jin EMAIL logo , Raytcho Lazarov , Dongwoo Sheen and Zhi Zhou


In this work, we consider the numerical solution of a distributed order subdiffusion model, arising in the modeling of ultra-slow diffusion processes. We develop a space semidiscrete scheme based on the Galerkin finite element method, and establish error estimates optimal with respect to data regularity in L2(Ω) and H1(Ω) norms for both smooth and nonsmooth initial data. Further, we propose two fully discrete schemes, based on the Laplace transform and convolution quadrature generated by the backward Euler method, respectively, and provide optimal L2(Ω) error estimates, which exhibits exponential convergence and first-order convergence in time, respectively. Extensive numerical experiments are provided to verify the error estimates for both smooth and nonsmooth initial data.

1 Introduction

We consider the following initial boundary value problem for a distributed order time fractional diffusion equation for u(x, t):


where Ω is a bounded convex polygonal domain in ℝd (d = 1, 2, 3) with a boundary Ω, v is a given function on Ω, and T > 0 is a fixed value. Here,

Dt[μ]u denotes the distributed order fractional derivative of u in time t (with respect to μ) defined by



tαu, 0 < α < 1, is the Caputo derivative of order α in t, defined by [16, p. 91]


where Γ(·) denotes Euler’s Gamma function. In this paper we consider a weight function μC [0, 1] with 0 ≤ μ < 1 and μ (0) μ (1) > 0.

In the last few decades, fractional calculus has been extensively studied and successfully employed to model subdiffusion, in which the mean squared variance grows slower than that in a Gaussian process. The subdiffusion model involves a Caputo derivative

tα0u of order α0 ∈ (0, 1) in t


describes subdiffusion processes. The model (1.2) is recovered from (1.1) with a singular weight μ (α) = δ (αα0), the Dirac delta function at α0. Physically, subdiffusion can be characterized by a unique diffusion exponent known as the Hurst exponent [3]. In practice, the physical process may not possess a unique Hurst exponent, and the model (1.1) provides a flexible framework for describing a host of continuous and nonstationary signals [3, 4, 29]. It is often employed to describe ultraslow diffusion, where the mean squared variance grows only logarithmically with time.

In recent years, the theoretical study of problem (1.1) has attracted some attention. Kochubei [17] made some early contributions to the rigorous analysis of the model (1.1), by constructing fundamental solutions and establishing their positivity and subordination property. Mainardi et al. [22] studied the existence of a solution, asymptotic behavior, and positivity etc. Meerschaert and Scheffler [25] gave a stochastic model for ultraslow diffusion; see also [31]. Luchko [21] showed a weak maximum principle for the problem. Li et al [18] established a sharp asymptotic behavior of the solution for t → 0 and t → ∞, in the case of continuous density μ with μ (1) > 0; see also [1] for further discussions.

The solution to (1.1) is rarely available in closed form, which necessitates the development of efficient numerical schemes. Despite the extensive studies on (1.2), there are only very few studies [6, 15, 26, 9, 7] on (1.1). Diethelm and Ford [6] developed a numerical scheme for distributed order fractional ODEs. It approximates

Dt[μ]u(t) by quadrature. This technique was also employed to solve nonlinear distributed-order fractional ODEs in [15], but without any analysis. Recently, Morgado and Rebelo [26] developed an implicit finite difference method for (1.1) with a Lipschitz nonlinear source term. The scheme approximates
by quadrature together with the backward finite difference in time, and the second-order finite difference in space. Its stability, and an error estimate O (h2 + τ +(δα)2) (with h, τ and δα being the mesh size, time step size and step size for quadrature rule, respectively) were established by assuming that the solution u is C2 in time and C4 in space and μ (α) is sufficiently regular. See also [9] a high-order difference scheme with v ≡ 0, and also [7]. The regularity requirement is restrictive. The development of numerical schemes for (1.1) with nonsmooth data and their rigorous analysis have not been carried out, despite its importance, e.g., in inverse and optimal control problems [14].

In this work, we develop a Galerkin finite element method (FEM) for problem (1.1). It is based on the finite element space Xh of continuous piecewise linear functions over a family of shape regular quasi-uniform partitions {Th}0 < h < 1 of Ω into d-simplexes, where h is the mesh size. Then the space semidiscrete Galerkin FEM is given by: find uh (t) ∈ Xh such that


with uh(0) = vh, where (·, ·) denotes the L2(Ω)-inner product, a (u, w) = (∇u, ∇w) for

u,wH01(Ω), and vhXh is an approximation of v. Further, we develop two fully discrete schemes based on the Laplace transform and convolution quadrature generated by the backward Euler method.

Our main contributions are as follows. First, in Theorem 2.1, we establish sharp regularity estimates for problem (1.1). Second, in Theorems 3.1 and 3.2, we derive optimal error estimates for the semidiscrete scheme (1.3). Third, we develop a fully discrete scheme based on the Laplace transform, using a contour representation of the semidiscrete solution with a hyperbolic contour and trapezoidal quadrature. We show its exponential convergence for a fixed time t, cf. Theorem 4.1. Last, we develop a second fully discrete scheme based on convolution quadrature, generated by the backward Euler method, and in Theorem 5.3, establish a first-order convergence. All these error estimates are nearly optimal and expressed in terms of the data regularity directly.

The model (1.1) is closely related to parabolic equations with a positive type memory term, for which there are many studies on numerical schemes based on convolution quadrature [19, 20, 5] and Laplace transform [27, 28, 32]. These interesting works have inspired the current work. However, they do not cover (1.1), due to the general kernel function involved. Instead, we shall opt for a direct strategy by bounding the kernel function.

The rest of the paper is organized as follows. In Section 2, we recall the solution theory of the model (1.1). In Section 3, we develop a space semidiscrete Galerkin scheme. Two fully discrete schemes are given in Sections 4 and 5. Finally, to verify the theory, we present in Section 6 some numerical experiments. Throughout, the notation c, with or without a subscript, denotes a generic constant, which may differ at different occurrences, but it is always independent of the mesh size h, the number N of quadrature points, and time step size τ . Further, we denote by ||·|| the L2(Ω) norm.

2 Solution Theory

In this part, we discuss the solution theory of problem (1.1). We denote by ˆ the Laplace transform. Next we denote by A the operator −Δ with a homogeneous Dirichlet boundary condition with a domain

D(A)=H01(Ω)H2(Ω). It is known that the operator A generates a bounded analytic semigroup of angle π/ 2, i.e., for any θ ∈ (π/2, π) [10, p. 321, Prop. C.4.2]



θ is a sector with the origin excluded, i.e., Σθ = {z ∈ ℂ : | arg (z) | < θ},
. Then by (1.1) and the relation



w(z)=01zα1μ(α)dα. Hence, u (t) is represented by


where H (z) = (zw(z) I + A )−1w(z), and Γθ,δ = {z ∈ ℂ : |z| = δ, |arg z| ≤ θ} ∪ {z ∈ ℂ : z = ρe±, ρδ}.

Now we give a few properties of w (z). The first is the sector-preserving property, since zα ∈ Σθ for any α ∈ (0, 1) and z ∈ Σθ.

Lemma 2.1

Let θ ∈ (π/2, π). Then zw(z) ∈ Σθ for all z ∈ Σθ.

The second result is an upper bound on w (z).

Lemma 2.2

Let μC [0, 1] with μ ≥ 0. Then there holds


The third result gives a lower bound on zw (z).

Lemma 2.3

Let θ ∈ (π/2, π) with μ(0) μ(1) > 0. Then there exists a constant c > 0 dependent only on θ and μ such that for any



Let z = reiϕ. Using μ (1) > 0 and μC[0, 1], we can find a small ϵ1 > 0 such that minαϵ[1 − ϵ1, 1]μ(α) ≥ δ1 > 0. Thus for r ≥ 1, there holds

01rαμ(α)dα1ϵ11rαμ(α)dαδ11ϵ11rαdαϵ1δ101rαdα. Similarly, we may find a small ϵ2 > 0 such that minαϵ[0, ϵ2]μ(α) ≥ δ2 > 0 and then for all r < 1, there holds
. Hence for ϕ ϵ(θπ, πθ), we get for c1 = min(ϵ1δ1, ϵ2δ2)


It suffices to consider ϕ ∈ [πθ, θ]. Since μ(0) > 0, there exists a small ϵ0 > 0 such that minα∈[0, ϵ0] cos(απ) μ(α) = δ0 > 0, and thus




Similarly, (2.3) holds for rr0 and ϕ ∈ [πθ, θ], thereby showing (2.3). The estimate (2.4) follows from

μC[0,1]01rαdα01rαμ(α)dα=|z|w(|z|) and the inequality |zw(z)| ≤ |z|w(|z|). □

Now we give the stability of problem (1.1).

Theorem 2.1

Let μC [0, 1] with μ ≥ 0 and μ(0) μ(1) > 0. Then the solution u to problem (1.1) satisfies for t ∈ (0, T] and ν = 0, 1:


where 1(t) = (ln(2T/t))−1, 2(t) = ln(max(t−1, 2) and cT > 0 may depend on d, Ω, μ, M, m and T.


The existence and uniqueness of a weak solution was shown in [18], and it suffices to show (2.5) and (2.6). First, by (2.1), we have ║H (z)║ = ║ (zw(z)I + A)−1║ |w(z)| ≤ M/|z| for all

zϵθ. Let t > 0, θ ∈ (π/2, π), δ > 0. We choose δ = t−1 and denote for short Γ = Γθ,δ. First we derive (2.5) for ν = 0 and m ≥ 0. By (2.2), we deduce


Next we prove (2.5) for ν = 1 and m ≥ 0, by taking δ = 2T/t in Γ. By applying A to (2.2) and differentiating with respect to t we arrive at


Since AH (z) = A (zw(z)I + A)−1w (z) = (I − zH (z))w (z), Lemma 2.2 yields ║AH(z)║ ≤ c|w(z)| ≤ c(|z| − 1)/(|z| ln |z|) for all

zΣ0. Hence (2.7) gives

AS(m)(t)cΓ|z|m|z|1|z|ln|z|e(z)t|dz|cδrm1r1ln rertcosθdr+cTtmδ1lnδθθe2Tcosψdψ=:I+II.

Since δ ≥ 2, we can bound the first term I by


The second term II is bounded by II = cTtm(δ − 1)/ln δcTtm− 11(t).

This shows (2.5). To prove (2.6) with ν = 0, we choose δ = t−1. Then,


Since zA−1H(z) = A−1 −(zw(z) I + A)−1 and ∫Γzm−1ezt dz = 0 for m ≥ 1, we have


By (2.1) and Lemma 2.3 we obtain ║zw(z) I + A )−1║ ≤ M|zw(z)|−1c(|z|−1)−1 ln |z|. This and the monotonicity of f(x) = − ln(x)/(1 − x ) on the positive real axis ℝ+ yield


For t−1 ≥ 2, ln(t−1)/(1 − t ) ≤ 2ln(t−1), and for t−1 < 2, ln(t−1)/(1 − t ) = ln(t)/(t − 1) ≤ 2ln(2). Thus ║S(m)(t)v║ ≤ ctm+12(t) ║Av║. Last, (2.6) with ν = 1 is equivalent to (2.5) with ν = 0 and v by Av. □

3 Semidiscrete Discretization by Galerkin FEM

Now we discuss the space semidiscrete scheme (1.3) based on the Galerkin FEM. On the space Xh, we define the L2(Ω) projection Ph: L2(Ω) → Xh and the Ritz projection

Rh:H01(Ω)Xh, respectively, by


With the discrete Laplacian Δh : XhXh defined by − (Δhϕ, χ) = (∇ϕ, ∇χ) for all ϕ, χXh, and Ah = − Δh, (1.3) can be rewritten as


with uh(0) = vhXh. For the error analysis, we employ an operator trick [8]. To this end, we first represent uh by


Now we introduce the error function e(t) := u(t) − uh(t), which, by (2.2) and (3.2), is given by



φ^(z)=(zw(z)I+A)1v and
. The next lemma shows a bound on
. It follows from Lemma 2.1, similar to [2, Lemmas 3.3 and 3.4], and hence the proof is omitted.

Lemma 3.1

Let vL2(Ω), z ∈ Σθ with θ ϵ(π/2, π),

φ^(z)=(zw(z)I+A)1v and
. Then there holds


Now we can state an error estimate for nonsmooth data vL2(Ω).

Theorem 3.1

Let u and uh be the solutions of problem (1.1) and (3.1) with vL2(Ω) and vh = Phv, respectively. Then for t > 0 and 1(t) = ln(2T/t)−1, there holds


By (3.3), with δ = 2 T/t in Γθ,δ, and Lemmas 3.1 and 2.2,

e(t)chδert cos θr1r ln rdrv+chθθe2T cos ψδ1ln δdψv:=I+II.

Now the first term I can be bounded by

Ichδert cos θln rdrvchln δδert cos θdrvcTht11(t)v,

and the second term II by

IIcTht ln δθθe2T cos ψdψvcTht11(t)v.

Thus the bound on ║∇e(t)║ follows. The L2 estimate is similar. □

Next we turn to the case of smooth data AvL2(Ω).

Theorem 3.2

Let u and uh be the solutions of problem (1.1) and (3.1) with vD(A) and vh = Rhv, respectively. Then for t > 0,


Like before, we take θ ∈ (π/2, π) and δ = 1 /t in the contour Γθ,δ. Then the error eh(t) = u(t) − uh(t) can be represented by


Using w(z)(zw(z) I + A )−1 = z−1Iz−1(zw(z) I + A )−1A, we deduce



φ^(z)=(zw(z)I+A)1Av and

Then Lemma 3.1 and the identity AhRh = PhA give


Now it follows that

eh(t)ch2Av(t1ert cos θr1dr+θθecos ψdψ)ch2Av,

which gives the L2-estimate. The H1 estimate is analogous. □

Remark 3.1

The estimate for vL2(Ω) deteriorates like t−11(t) as t → 0+, which agrees with Theorem 2.1. The factor t−11(t) is different from that for subdiffusion [12] and multi-term counterpart [11]. For AvL2(Ω), the estimate is uniform in time.

4 Fully Discrete Scheme I: Laplace Transform

The first fully discrete scheme is based on Laplace transform, by applying a quadrature rule to the representation (3.2). We follow [24, 27, 28, 32] and deform the contour Γθ,δ to be a curve with a parametric representation


with λ > 0, ψ ∈ (0, π/2) and ξ ∈ ℝ. The optimal choices of λ and ψ will be given in Lemma 4.4. This deformation is valid since it does not transverse the poles of H(z)v = (zw(z) + Ah)−1w(z)v, cf., Lemmas 2.1 and 4.3. With z = x + iy, the contour (4.1) is the left branch of the hyperbola

((xλ)/(λ sin ψ))2(y/(λ cos ψ))2=1,

which intersects the real axis at x = λ(1 ± sin ψ ) and has asymptotes y = ±(λx) cot ψ. Now we can represent uh(t) by


with the integrand ĝ(ξ, t) being defined by


Now we describe the quadrature approximation. By setting zj = z(ξj) and

zj:=z(ξj) with ξj = jk and k being the step size, we have


and the truncated quadrature approximation


with φ̂j = (zjw(zj)I + Ah)−1w(zj)vh. To compute UN,h(t), we need to solve only N+1 elliptic problems, instead of 2N+1, by exploiting conjugacy:

zj=zj¯,w(zj)=w(zj)¯,ϕ^j=ϕ^¯j,j=1,...,N. Indeed, since
, with ζj = λcos(iξjψ), (4.4) is reduced to


Hence we solve the following complex–valued elliptic problems


Next we define a strip by 𝓢a,b = {p = ξ + iη : ξ ∈ ℝ, η ∈ (−b, a)}. The following lemma gives the quadrature error [23] [32, Theorem 2.1].

Lemma 4.1

Let g be an analytic function in a strip 𝓢a,b for some a, b > 0, and I and Ik, for k > 0, be defined by

, respectively. Furthermore, assume that g(z) → 0 uniformly as |z| → ∞ in 𝓢a,b, and that there exist M+,M > 0 such that


Then with E+ = M+/(e2/k − 1) and E = M/(e2/k − 1), there holds


The next lemma gives one crucial estimate on the map z(p) over the strip 𝓢a,b. Even though the hyperbolic contour (4.1) has been extensively used, the estimate on the map z(p) below seems to be new.

Lemma 4.2

Let p = ξ + iη with ξ, η ∈ ℝ. Then with a = π/2 − ψϵ and b = ψϵ, for small ϵ > 0, there holds


For p = ξ + iη with ξ, η ∈ ℝ, then the image z(p) in the parameterization (4.1) is given by z(p) = λ(1 − sin(ψ + η)cosh(ξ))+iλ cos(ψ + η) sinh(ξ), and its derivative z′(p) is given by z′(p) = λ cosh ξcos(ψ +η) − i sinh ξ sin(ψ + η). By writing z = x + iy, it can be expressed as the left branch of the hyperbola ((xλ)/(λ sin(ψ + η)))2−(y/(λcos(ψ +η )))2 = 1 . It intersects the real axis at x = λ(1− sin(ψ + η )) and has the asymptotes y = ± (xλ)cot(ψ + η). Next we show (4.5) and (4.6). First, for

pϵS¯a,0, i.e., η ∈ [0,a], z(p) lies in the sector Σπ−ψ. Since ϕ := η + ψ ∈ (ψ, π/2 − ϵ), and sin(π/2 − ϵ) ∼ 1 − ϵ2/2 ≤ 1 − ϵ2/3 for small ϵ, we have for all ξ ∈ ℝ


Hence (4.5) holds true. For

pS¯0,b, i.e., η ∈ [−b,0], z(p) lies in Σπ−(η+ψ) ⊂ Σπ−ϵ. Further, since ϕ := η + ψ ∈ (ϵ, ψ), for all ξ ϵ ℝ


Then the desired result (4.6) follows directly. □

The next result gives analyticity of and an estimate on ĝ(ξ, t) on 𝓢a,b.

Lemma 4.3

Let p = ξ + iη with ξ, η ϵ ℝ, a and b be defined as in Lemma 4.2, and ĝ(p, t) be defined by (4.3). Then ĝ(p, t) is analytic on the strip 𝓢a,b, and there holds


For p = ξ + iη with ξ, η ∈ ℝ, z(p) in (4.1) is given by


By Lemmas 4.2 and 2.1, z(p)w(z(p)) ∈ Σπ−ϵ, and thus the function


is analytic in 𝓢a,b. It sufficies to show (4.7). For

pϵS¯0,b, by (4.6), z(p) ∈ Σπ−ϵ. By Lemma 2.1, z(p)w(z(p)) ∈ Σπ−ϵ. By (2.1), for small ϵ > 0


For any p𝓢0,b, 𝓡(z(p)) = λ(1 − sin(ψ + η)cosh(ξ)). By Lemma 2.1,


This and (4.6) yield (4.7). The case

pϵS¯a,0 is direct: (4.6) and Lemma 2.1 imply z(p)w(z(p))ϵ Σθ. Then (4.7) follows from (4.5) and (2.1). □

Now we can give an error estimate on the approximation UN,h.

Lemma 4.4

Let uh(t) and UN,h(t) be defined in (4.2) and (4.4), respectively, and the contour be parametrically represented by (4.1). Then with the choice k = c0/N and λ = c1N/t, there holds


Let uhUN,h = (uhUh) + (UhUN,h) =: Eq + Et. Let a = π/2 − ψϵ and b = ψϵ. For p = ξ + ia, zw(z) lies in Σθ for some θ ∈ (π/2, π). Since cosh ξ ≥ 1 + ξ2/2 and 1 − sin(π/2 − ϵ) ≤ ϵ for small ϵ > 0, the choice λ = c1N/t and Lemma 4.3 yield


By Lemma 4.1, for k = c0/N we have


Next we bound the error due to the lower half. For the choice p = ξ − ib, λ = c1N/t and by the inequality cosh ξ ≥ 1 + ξ2/2, we deduce


Then for the choice k = c0/N, Lemma 4.1 yields the following estimate

Eqcϵ3/2N1/2e(2π(ψϵ)/c0c1(1ϵ))Nvh. Further, since cosh(ξ) ≥ cosh(c0)+sinh(c0)(ξc0) for ξc0, the error ║Et║ is bounded by


By disregarding terms ε, balancing the exponentials in

||Eq+||,||Eq|| and ║Et║, we get 2π(π/2−ψ)/c0 = 2πψ/c0−c1 = −c1(1−sin(ψ)cosh(c0)). Next we express c0 and c1 using ψ : c0 = cosh−1(2πψ/((4πψ − π2)sin ψ)) and c1 = (4πψ − π2)/ cosh−1(2πψ/((4πψ − π2)sin ψ)). Finally we minimize the ratio B(ψ) = c1−2πψ/c0 with respect to ψ, which achieves the minimum at ψ = 1.1721 and hence, c0 = 1.0818, c1 = 4.4920 and B(ψ) = − 2.32, which are identical to the values in [32]. Then collecting the balanced asymptotic bound and
, and choosing ϵ = 1/N yield


Since x−1 ln x ≤ 1/e for x ≥ 1, the L2-stability of Ph yields (4.8). □

Last, we give error estimates for the fully discrete scheme (4.4). It follows from Theorems 3.1 and 3.2, and Lemma 4.4.

Theorem 4.1

Let u(t) be the solution of problem (1.1), and UN,h(t) be the approximation (4.4), with the parameters chosen as in Lemma 4.4. Then with 1(t) = (ln 2T/t)−1, the following statements hold.

  1. If AvL2(Ω) and vh = Rhv, then

  2. If vL2(Ω) and vh = Phv, then


5 Fully Discrete Scheme II: Convolution Quadrature

Now we develop a second fully discrete scheme based on convolution quadrature. We divide the interval [0,T] into a uniform grid with a time step size τ = T/N, N ∈ ℕ, with tn = , n = 0, ..., N. Following [19, 5], we consider the convolution quadrature generated by the backward Euler method. The weights

{bj}j=0 are given by [19]
. Then the convolution quadrature approximation is given by: for n = 1, 2,..., N:
. Then we get a fully discrete scheme for (1.1): with
Uh0 = vh
, for n = 1, 2, ..., N


We denote the generating function

β˜ of

Now we analyze (5.1), following [20]. First we split the error into


By Theorems 3.1 and 3.2, it suffces to bound

uh(tn)Uhn. We write
, where yh(t) = uh(t) − vh and

First, we derive representations of yh and Yh.

Lemma 5.1



Then yh and

Yhncan be represented by


respectively, with the contour Γτ = {z ∈ Γθ,δ : |𝔉(z)| ≤ π/τ}.


By definition, yh satisfies

Dt[μ]yh+Ahyh=Ahvh, with yh(0) = 0. The Laplace transform gives zw(z)ŷh(z) + Ah ŷh(z)= −(z)−1Ahυh. Hence, ŷ(z) = K(z)υh, with K(z) = −z−1(zw(z)I + Ah)−1Ah, and thus follows the representation of yh(t). Likewise,


Multiplying both sides by ξn and summing in n from 1 to ∞ yield


Using the condition

Yh0=0, we have


Thus, simple calculation shows Ỹh(ξ) = (ξ/τ)K((1 − ξ))υh and it is analytic at ξ = 0. Then Cauchy theorem implies that for small ϱ


Now, by changing variable ξ = e−zτ, we obtain


where the contour Γ0 = {z = −ln(ϱ)/τ + iy : |y| ≤ π/τ} is oriented counterclockwise. The desired representation follows by deforming the contour Γ0 to Γτ = {z ∈ Γθ,δ : |𝔍 (z)| ≤ π/τ}.

By Lemma 5.1, we have

yh(tn)Yhn=I+II, where


Since |e| is uniformly bounded on Γτ, we have

K(z)ezτK(χ(z))|ezτ|K(z)K(χ(z))+|1ezτ|K(z)                                    cK(z)K(χ(z))+cτ|z|K(z)                                    cK(z)K(χ(z))+cτ,

where the last line, using (2.1), follows from the inequality


Hence, it remains to bound the term ║K(z) − K(χ(z))║. First we recall a bound on χ(z) = τ−1(1 − e−zτ) [13, Lemma 3.1].

Lemma 5.2

Let χ(z) = τ−1(1 − e). Then for all z ∈ Γτ, there hold


and χ(z) lies in a sector Σθ′ for some θ′ ∈ (π/2, π).

Next we give one estimate on the approximation χ(z)ω(χ(z)) to (z).

Lemma 5.3

For z ∈ Γτ, there holds


By the intermediate value theorem, for z ∈ Γτ, we have


where zη = ηχ(z) + (1 − η)z with η ∈ [0, 1]. Next we claim |zη|−1c|z|−1 for η ∈ [0, 1]. We split Γτ into

Γτ=Γτ+ΓτcΓτ, with
being the rays in the upper and lower half plane, and
is the circular arc. For
, Taylor expansion gives
. Since || ≤ 1 for
, |zη|−1c|z|−1 for
. It remains to consider
. For z = rei(πθ) with ∈ (δ, π/ sinθ) we have
, and since sin θπ, 𝔍(χ(z)) ≥ 0. Then Lemma 5.2 yields
. This shows the desired claim. Hence, by Lemma 5.2, for z ∈ Γτ there holds

Next we give an error estimate on the approximation K(χ(z)) to K(z).

Lemma 5.4

For z ∈ Γτ, there holdsK(z) − K(χ(z))║ ≤ .


Let B(z) = zK(z). Simple computation shows

B(z)B(χ(z))=zw(z)((zw(z)I+Ah)1(χ(z)w(χ(z))I+Ah)1)          +(zw(z)χ(z)w(χ(z)))(χ(z)w(χ(z))I+Ah)1:=I+II.

First, by Lemmas 2.3 and 5.2, there holds |χ(z)w(χ(z))| ≥ c|χ(z)|w(|χ(z)|) ≥ c|z|ω(|z|). Further, by Lemma 2.1, (2.1) and Lemma 2.3, we have


Likewise, in view of Lemmas 5.2 and 2.1 and (2.1), we have


Now, by the identity


Lemma 2.3, (5.4) and (5.5), the first term I can be bounded by


Likewise, by Lemma 5.3 and (5.5)


Hence, ║B(z) − B(χ(z))║ ≤ cτ |z|. Last, by Lemma 5.2 and B (z) ≤ c,


Now we can state an error estimate for nonsmooth data vL2(Ω).

Theorem 5.1

Let uh and

Uhnbe the solutions of problems (3.1) and (5.1) with vL2(Ω), vh = Phv and f ≡ 0, respectively. Then there holds


It suffices to bound the terms I and II in (5.3). By choosing

δ=tn1 in Γδ, θ and (2.1), we bound I by


Using Lemma 5.4, we deduce



yh(tn)Yhncτtn1vh, and the estimate follows from
and the L2 stability of Ph. □

Next we turn to smooth data AvL2(Ω).

Lemma 5.5

Let Ks(z) = −z−1(zw(z)I + Ah)−1. Then for any z ∈ Γτ


Let Bs(z) = − (zw(z)I + Ah)−1. Then by the identity Bs(z) − Bs(χ(z)) = (χ(z)w(χ(z)) − zw(z)) (zw(z)I + Ah)−1 (χ(z)w(χ(z))I + Ah)−1, Lemma 5.3, and (5.4) and (5.5), we deduce Bs(z) − Bs(χ(z)) ≤ |w (z)| −1. Appealing to (5.4) again yields Bs(z) ≤ c|zw(z)|−1, and thus


Then the desired result follows from Lemma 2.3.

Theorem 5.2

Let uh and

Uhnbe the solutions of problems (3.1) and (5.1) with AvL2(Ω) and v h = Rhv, respectively. Then with 𝓁2(t) = ln max(t−1, 2)),


Let Ks(z) = − z−1 (zw(z) I + Ah)−1. Then


By Lemma 5.5, ║Ks(z) − eKs (χ(z))║ ≤ ln |z| / (|z| − 1) for z ∈ Γτ. With

δ=tn1, the monotonicity of f(x) = (1 − x)−1 ln x yields


By (2.1), for z ∈ Γθ, δ, ║Ks(z)≤ c|z|−1 |zw(z)| −1. Now Lemma 2.3 yields


Now the identity AhRh = PhA completes the proof. □

The next theorem gives error estimates for the scheme (5.1), which follow from Theorems 3.1, 3.2, 5.1 and 5.2.

Theorem 5.3

Let u and

Uhnbe the solutions of problems (1.1) and (5.1) with
respectively. Then the following statements hold.
  1. If AvL2(Ω) and vh = Rhv, then with 𝓁1(t) = ln(2 T/ t)−1

  2. If vL2(Ω) and vh = Phv, then with ℓ2(t) = ln(max(t−1, 2))


6 Numerical Experiments and Discussions

Now we present numerical results to verify the convergence theory. Throughout, let the domain Ω = (0, 1) and consider

  1. v(x)=sin(2πx)H2(Ω)H01(Ω)

  2. v = 𝓧(0, 1/2)H1/2 − ε (Ω) with ε ∈ (0, 1/2), and 𝓧S the characteristic function of a set S;

  3. v(x) = x−1/4H1/4 − ε(Ω) with ∈ (0, 1/4).

We measure the temporal error by the L2 errors ║u(tn) − UN, h(tn)║ / ║v║ or

u(tn)Uhn/υ, and the spatial error by ║ u(t) − uh(t)║ /║v║ and ║∇ (u(t) − uh(t))║/║v║. Throughout we divide the domain Ω into M equally spaced subintervals with a mesh size h = 1 / M. Since the exact solution u(t) is not unavailable, we compute a reference solution using a finer mesh.

6.1 Numerical results for the semidiscrete scheme

First we examine the convergence behavior of (1.3). Here we fix N = 10 in the Laplace transform approach such that the error due to time discretization is negligible. The numerical results are given in Table 1. In the table, rate denotes the empirical rates when the mesh size h halves, and the numbers in the bracket denote the theoretical rates. For all cases, the L2 and H1 norms of the error exhibit O(h2) and O(h) convergence, respectively, agreeing with Theorems 3.1 and 3. 2. The scheme is robust since the convergence rates hold for both smooth and nonsmooth data. The error increases as t → 0, due to the weak solution singularity around t → 0, cf. Theorem 2.1.

6.2 Numerical results for the fully discrete scheme I

Next we illustrate the convergence of the scheme (4.4). To make the spatial discretization error negligible, we fix h at h = 10−5. The numerical results are summarized in Tables 2 and 3 for μ(α) = (α − 1/2)2 and μ(α) = 𝓧[1/2, 1](α), respectively. The notation r denotes the exponent in the estimate

||u(N,h)nu(tn)||CerN. The results indicate an exponential convergence with respect to N, with a rate e− 2.15N and e− 2.14N for μ(α) = (α − 1/2)2 and μ(α) = 𝓧[1/2, 1](α), respectively, which agree with Theorem 4.1. Note that even though μ (α) = 𝓧[1/2, 1](α) does not satisfy the condition μ(0) μ(1) > 0, the empirical rates still agree well with the theoretical one, which calls for further study. The convergence rate is independent of t, and thus the scheme is robust. Interestingly, the smoothness of the initial data does not affect much the time discretization errors, even for small time, cf. Table 4.

Table 1

Numerical results for the semidiscrete scheme (1.3) with μ(α) = (α − 1/2)2.

Table 2

The L2 errors for (a)-(c) with h = 10−5 and μ(α) = (α − 1/2)2, by the Laplace transform method.

caset \ N35791113r

One salient feature of the fully discrete scheme I is that it allows computing the solution at large time. This allows one to examine the asymptotic behavior of the solution as the time t → ∞; see Table 5 and Fig. 1. In particular, one clearly observes the logarithmic decay of the solution [18, Theorem 2.1]; see also Fig. 1. This numerically verifies the ultraslow decay asymptotics for distributed order diffusion process.

Table 3

The L2 errors for (a)-(c) with h = 10 − 5 and μ(α) = χ[1/2, 1](α), by the Laplace transform method.

caset \ N35791113r
Table 4

The L2 errors for (b) and (c) with h = 10−5, μ(α) = (α−1/2)2 and N = 5 at time t = 10k, k = 4,...,9, by the Laplace transform method.

Table 5

The L2 norm of the solution for (a) and (c) with h = 10−5, μ (α) = (α − 1/2)2 and N = 10 at time t = 10k, k = 6, 8, ···, 18, by the Laplace transform method.

(a)3.33e-42.70e-42.26e-41.95e-41.71e-41.52e-41.37e-41 / k
(c)1.06e-38.54e-47.17e-46.17e-45.41e-44.82e-44.34e-41 / k
Figure 1 The L2 norm of the solution for (a) and (c) at t = 10k, k = 6, 8,....., 18, by the Laplace transform method.
Figure 1

The L2 norm of the solution for (a) and (c) at t = 10k, k = 6, 8,....., 18, by the Laplace transform method.

6.3 Numerical results for the fully discrete scheme II

Last we verify the convergence of the fully discrete scheme II, i.e., (5.1). The results are shown in Tables 6 and 7 for μ(α) = (α − 1/2)2 and μ (α) = 𝓧[1/2, 1](α), respectively. An O(τ) convergence is always observed, cf. Theorem 5.3.

To examine more closely its convergence behavior, we consider t = 10k, k = 4,...,9, and at each time t, divide the interval [0, t] into N = 10 subintervals. The scheme works well for the smooth initial data, however, it works poorly for the singular initial data, cf. Table 8. This behavior is predicted by Theorem 5.1: the error is dominated by the factor τ/t for L2 initial data. In Fig. 2, we plot the ratio

Uh1u(τ)/τ versus ln τ for smooth initial data. Theorem 5.2 predicts
. The log factor 2(t) in Theorem 5.2 is confirmed by Fig. 2, and thus it is sharp.

Figure 2 The L2 errors for (a) at t1 = τ = 10−k, k = 5,..., 12, by convolution quadrature.
Figure 2

The L2 errors for (a) at t1 = τ = 10k, k = 5,..., 12, by convolution quadrature.

Table 6

The L2 errors for (a)-(c) with h = 10−4 and μ(α) = (α − 1/2)2, by convolution quadrature.

caset \ N10204080160320rate
(a)11.82e-58.78e-64.31e-62.12e-61.01e-64.74e-71.05 (1.00)
10−28.64e-43.91e-41.88e-49.20e-54.55e-52.26e-51.05 (1.00)
10−32.17e-21.10e-25.51e-32.76e-31.38e-36.92e-40.99 (1.00)
(b)14.81e-52.32e-51.14e-55.60e-62.67e-61.26e-61.05 (1.00)
10−28.11e-33.87e-31.88e-39.29e-44.61e-42.30e-41.03 (1.00)
10−31.48e-27.46e-33.74e-31.88e-39.39e-44.70e-41.00 (1.00)
(c)15.81e-52.81e-51.38e-56.76e-63.23e-61.52e-61.05 (1.00)
10−21.01e-24.80e-32.34e-31.15e-35.72e-42.85e-41.03 (1.00)
10−37.35e-33.66e-31.82e-39.11e-44.55e-42.27e-41.00 (1.00)
Table 7

The L2 errors for (a)-(c) with h = 10−4 and μ(α) = χ[1/2, 1](α), by convolution quadrature.

caset \ N10204080160320rate
(a)12.20e-41.06e-45.20e-52.58e-51.28e-56.40e-61.02 (1.00)
10−21.76e-28.81e-34.40e-32.20e-31.10e-35.49e-41.00 (1.00)
10−33.92e-31.98e-39.95e-44.99e-42.50e-41.25e-40.99 (1.00)
(b)16.52e-43.11e-41.52e-47.53e-53.74e-51.87e-51.03 (1.00)
10−21.25e-26.26e-33.13e-31.56e-37.82e-43.91e-41.00 (1.00)
10−35.76e-32.88e-31.44e-37.18e-43.59e-41.79e-41.00 (1.00)
(c)17.92e-43.78e-41.85e-49.14e-54.54e-52.27e-51.03 (1.00)
10−27.40e-33.71e-31.86e-39.28e-34.64e-42.32e-41.00 (1.00)
10−36.10e-33.06e-31.53e-37.65e-43.83e-41.91e-41.00 (1.00)
Table 8

The L2 errors with h = 10−5 and N = 10, at t = 10k, k = 4,..., 9, by convolution quadrature.

case \ t10−410−510−610−710−810−9rate
(a)2.42e-31.03e-47.87e-67.59e-77.58e-87.44e-91.01 (1.00)
(c)7.44e-35.67e-34.30e-33.27e-32.49e-31.88e-30.12 (0.12)


The research of B. Jin was partially supported by UK EPSRC Grant EP/M025160/1, that of R. Lazarov and Z. Zhou by NSF Grant DMS-1016525, and that of D. Sheen by NRF-2014R1A2A1A11052429.


[1] E. Bazhlekova, Completely monotone functions and some classes of fractional evolution equations. Integral Transforms Spec. Funct. 26, No 9 (2015), 737–752; DOI:10.1080/10652469.2015.1039224.10.1080/10652469.2015.1039224Search in Google Scholar

[2] E. Bazhlekova, B. Jin, R. Lazarov, Z. Zhou, An analysis of the Rayleigh-Stokes problem for the generalized second grade fluid. Numer. Math. 131, No 1 (2015), 1–31; DOI:10.1007/s00211-014-0685-2.10.1007/s00211-014-0685-2Search in Google Scholar

[3] A. V. Chechkin, R. Gorenflo, I. M. Sokolov, Retarding subdiffusion and accelerating superdiffusion governed by distributed-order fractional diffusion equations. Phys. Rev. E66 (2002), 046129; DOI:10.1103/PhysRevE.66.04612910.1103/PhysRevE.66.046129Search in Google Scholar

[4] A. V. Chechkin, R. Gorenflo, I. M. Sokolov, V. Y. Gonchar, Distributed order time fractional diffusion equation. Fract. Calc. Appl. Anal. 6, No 3 (2003), 259–279.Search in Google Scholar

[5] E. Cuesta, C. Lubich, C. Palencia, Convolution quadrature time discretization of fractional diffusion-wave equations. Math. Comp. 75, No 254 (2006), 673–696; DOI: 10.1090/S0025-5718-06-01788-1.10.1090/S0025-5718-06-01788-1Search in Google Scholar

[6] K. Diethelm, N. J. Ford, Numerical analysis for distributed-order differential equations. J. Comput. Appl. Math. 225, No 1 (2009), 96–104; DOI: 10.1016/ in Google Scholar

[7] N. J. Ford, M. L. Morgado, M. Rebelo, An implicit finite difference approximation for the solution of the diffusion equation with distributed order in time. Electron. Trans. Numer. Anal. 44 (2015), 289–305.Search in Google Scholar

[8] H. Fujita, T. Suzuki, Evolution problems. In: Handbook of Numerical Analysis , Vol. II, North-Holland, Amsterdam (1991), 789–928.10.1016/S1570-8659(05)80043-2Search in Google Scholar

[9] G. Gao, H. Sun, Z. Sun, Some high-order difference schemes for the distributed-order differential equations. J. Comput. Phys. 298, No 1 (2015), 337–359; DOI: 10.1016/ in Google Scholar

[10] M. Hasse, The Functional Calculus for Sectorial Operators. Birkhäuser, Basel (2006).10.1007/3-7643-7698-8Search in Google Scholar

[11] B. Jin, R. Lazarov, Y. Liu, Z. Zhou, The Galerkin finite element method for a multi-term time-fractional diffusion equation. J. Comput. Phys. 281 (2015), 825–843; DOI: 10.1016/ in Google Scholar

[12] B. Jin, R. Lazarov, Z. Zhou, Error estimates for a semidiscrete finite element method for fractional order parabolic equations. SIAM J. Numer. Anal. 51, No 1 (2013), 445–466; DOI:10.1137/12087398410.1137/120873984Search in Google Scholar

[13] B. Jin, R. Lazarov, Z. Zhou, An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data. IMA J. Numer. Anal. 36, No 1 (2016), 197–221; DOI: 10.1093/imanum/dru063.10.1093/imanum/dru063Search in Google Scholar

[14] B. Jin, W. Rundell, A tutorial on inverse problems for anomalous diffusion processes. Inverse Problems31, No 3 (2015), 035003; DOI: 10.1088/0266-5611/31/3/035003.10.1088/0266-5611/31/3/035003Search in Google Scholar

[15] J. T. Katsikadelis. Numerical solution of distributed order fractional differential equations. J. Comput. Phys. 259 (2014), 11–22; DOI: 10.1016/ in Google Scholar

[16] A. Kilbas, H. Srivastava, J. Trujillo, Theory and Applications of Fractional Differential Equations . Elsevier, Amsterdam (2006).Search in Google Scholar

[17] A. N. Kochubei, Distributed order calculus and equations of ultra-slow diffusion. J. Math. Anal. Appl. 340, No 1 (2008), 252–281; DOI: 10.1016/j.jmaa.2007. in Google Scholar

[18] Z. Li, Y. Luchko, M. Yamamoto, Asymptotic estimates of solutions to initial-boundary-value problems for distributed order time-fractional diffusion equations. Fract. Calc. Appl. Anal. 17, No 4 (2014), 1114–1136; DOI: 10.2478/s13540-014-0217-x; in Google Scholar

[19] C. Lubich, Convolution quadrature and discretized operational calculus, I. Numer. Math. 52, No 2 (1988), 129–145; DOI: 10.1007/BF01398686.10.1007/BF01398686Search in Google Scholar

[20] C. Lubich, I. H. Sloan, V. Thomée, Nonsmooth data error estimates for approximations of an evolution equation with a positive-type memory term. Math. Comp. 65, No 213 (1996), 1–17; DOI: 10.1090/S0025-5718-96-00677-1.10.1090/S0025-5718-96-00677-1Search in Google Scholar

[21] Y. Luchko, Boundary value problems for the generalized time fractional diffusion equation of distributed order. Fract. Calc. Appl. Anal. ,12, No 4 (2009), 409–422; at fcaa.∼fcaaSearch in Google Scholar

[22] F. Mainardi, A. Mura, G. Pagnini, R. Gorenflo, Time-fractional diffusion of distributed order. J. Vibr. Control14, No 9–10 (2008), 1267– 1290; DOI: 10.1177/1077546307087452.10.1177/1077546307087452Search in Google Scholar

[23] E. Martensen, Zur numerischen Auswertung uneigenlicher Integrale. Z. Angew. Math. Mech.48 (1968), T83–T85.Search in Google Scholar

[24] W. McLean, V. Thomée, Numerical solution via Laplace transforms of a fractional order evolution equation. J. Integral Equations Appl. 22, No 1 (2010), 57–94; DOI: 10.1216/JIE-2010-22-1-57.10.1216/JIE-2010-22-1-57Search in Google Scholar

[25] M. M. Meerschaert, H.-P. Scheffer, Stochastic model for ultraslow diffusion. Stochastic Process. Appl. 116, No 9 (2006), 1215–1235; DOI: 10.1016/ in Google Scholar

[26] M. L. Morgado, M. Rebelo, Numerical approximation of distributed order reaction diffusion equations. J. Comput. Appl. Math. 275 (2015), 216–227; DOI: 10.1016/ in Google Scholar

[27] D. Sheen, I. Sloan, V. Thomée, A parallel method for time-discretization of parabolic problems based on contour integral representation and quadrature. Math. Comp. 69, No 229 (2000), 177–195; DOI: 10.1090/S0025-5718-99-01098-4.10.1090/S0025-5718-99-01098-4Search in Google Scholar

[28] D. Sheen, I. Sloan, V. Thomée, A parallel method for time discretization of parabolic equations based on Laplace transformation and quadrature. IMA J. Numer. Anal. 23, No 2 (2003), 269–299; DOI:10.1093/imanum/ in Google Scholar

[29] I. M. Sokolov, A. V. Chechkin, J. Klafter, Distributed-order fractional kinetics. Acta Phys. Polon. B35, No 4 (2004), 1323–1341.Search in Google Scholar

[30] V. Thomée, Galerkin Finite Element Methods for Parabolic Problems . Springer-Verlag, Berlin (2006).Search in Google Scholar

[31] S. Umarov, Continuous time random walk models associated with distributed order diffusion equations. Fract. Calc. Appl. Anal. 18, No 3 (2015), 821–837; DOI: 10.1515/fca-2015-0049; in Google Scholar

[32] J. A. C. Weideman, L. N. Trefethen, Parabolic and hyperbolic contours for computing the Bromwich integral. Math. Comp. 76, No 259 (2007), 1341–1356; DOI: 10.1090/S0025-5718-07-01945-X.10.1090/S0025-5718-07-01945-XSearch in Google Scholar

  1. Please cite to this paper as published in:

    Fract. Calc. Appl. Anal., Vol. 19, No 1 (2016), pp. 69–93, DOI: 10.1515/fca-2016-0005

Received: 2015-6-5
Published Online: 2016-3-22
Published in Print: 2016-3-1

© 2016 Diogenes Co., Sofia

Downloaded on 2.12.2023 from
Scroll to top button