Error estimates for approximations of distributed order time fractional diffusion with nonsmooth 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,

where

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

describes subdiffusion processes. The model (1.2) is recovered from (1.1) with a singular weight μ (α) = δ (α − α₀), the Dirac delta function at α₀. 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)

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

with u_h(0) = v_h, where (·, ·) denotes the L²(Ω)-inner product, a (u, w) = (∇u, ∇w) for

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 L²(Ω) 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

where

Σθ′ = Σθ\{0}

∂tαu^=zαu^−zα−1u(0)

with

where H (z) = (zw(z) I + A )⁻¹w(z), and Γ_θ,δ = {z ∈ ℂ : |z| = δ, |arg z| ≤ θ} ∪ {z ∈ ℂ : z = ρe^±iθ, ρ ≥ δ}.

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

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

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

z∈∑′θ

(2.3)|zw(z)|≥c(|z|−1)/ln  |z|,

(2.4)|z|w(|z|)≥|zw(z)|≥c|z|w(|z|).

Proof

Let z = re^iϕ. Using μ (1) > 0 and μ ∈ C[0, 1], we can find a small ϵ₁ > 0 such that min_{αϵ[1 − ϵ1, 1]}μ(α) ≥ δ₁ > 0. Thus for r ≥ 1, there holds

∫01rαμ(α)dα≥∫1−ϵ11rαμ(α)dα ≥ δ1∫1−ϵ11rα dα  ≥ ϵ1δ1∫01rαdα. Similarly, we may find a small ϵ₂ > 0 such that min_{αϵ[0, ϵ2]}μ(α) ≥ δ₂ > 0 and then for all r < 1, there holds

∫01rαμ(α)dα≥ϵ2δ2∫01rαdα. Hence for ϕ ϵ(θ −π, π − θ), we get for c₁ = min(ϵ₁δ₁, ϵ₂δ₂)

|zω(z)|≥ℜ(zω(z))≥cos(π−θ)∫01rαμ(α)dα≥c1cos(π−θ)∫01rαdα.

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

ℜ(zw(z))≥∫0ϵ0rαcos(αφ)μ(α)dα−∫ϵ01rα|cos(αφ)|μ(α)dα≥δ0∫0ϵ0rαdα−‖μ‖C[0,1]∫ϵ01rαdα≥−(In r)−1[δ0−rϵ0(δ0+‖μ‖C[0,1])].

Consequently,

|zw(z)|≥ℜ(zw(z))≥δ02∫01rαdα  ∀r ≤ r0=:(δ0/(2(δ0+‖μ‖C[0,1])))1/ϵ0.

Similarly, (2.3) holds for r ≥ r₀ 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).

Proof

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)⁻¹║ |w(z)| ≤ M/|z| for all

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

‖S(m)(t)‖ =‖12πi∫ΓzmeztH(z)dz‖≤c∫Γ|z|meℜ(z)t‖H(z)‖|dz|≤c(∫t−1∞rm−1ert  cos θdr+∫−θθecos  ψt−mdψ)≤ ct−m.

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

(2.7)AS(m)(t)=12πi∫ΓzmeztAH(z)dz.

Since AH (z) = A (zw(z)I + A)⁻¹w (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∫δ∞rm−1r−1ln rert cos  θdr+cTt−mδ−1lnδ∫−θθe2Tcosψdψ=:I+II. 

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

I ≤  c∫δ∞ rm (ln r)−1ert cos θ dr   ≤ cℓ1 (t)∫δ∞  rm ert cos θdr ≤ cTt−m−1  ℓ1 (t).

The second term II is bounded by II = c_Tt^−m(δ − 1)/ln δ ≤ c_Tt^{−m− 1}ℓ₁(t).

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

S(m)(t)υ=12πi∫ΓzmeztH(z)υdz=12πi∫Γzm−1eztzA−1H(z)Aυdz.

Since zA⁻¹H(z) = A⁻¹ −(zw(z) I + A)⁻¹ and ∫_Γz^m−1e^zt dz = 0 for m ≥ 1, we have

S(m)(t)v=12πi∫Γzm−1eztvdz−12πi∫Γzm−1ezt(zw(z)I+A)−1dzAv=−12πi∫Γzm−1ezt(zw(z)I+A)−1dzAv.

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

‖s(m)(t)v‖≤c(∫Γ|z|m−1eℜ(z)t‖zw(z)I+A−1‖|dz|)‖Av‖≤c(∫t−1∞ertcosθrm−1lnrr−1dr+t−mlnt−1t−1−1∫−θθecosψdψ)‖Av‖≤ct−m+1ln(t−1)(1−t)−1‖Av‖.

For t⁻¹ ≥ 2, ln(t⁻¹)/(1 − t ) ≤ 2ln(t⁻¹), and for t⁻¹ < 2, ln(t⁻¹)/(1 − t ) = ln(t)/(t − 1) ≤ 2ln(2). Thus ║S^(m)(t)v║ ≤ ct^−m+1 ℓ₂(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 X_h, we define the L²(Ω) projection P_h: L²(Ω) → X_h and the Ritz projection

With the discrete Laplacian Δ_h : X_h → X_h defined by − (Δ_hϕ, χ) = (∇ϕ, ∇χ) for all ϕ, χ ∈ X_h, and A_h = − Δ_h, (1.3) can be rewritten as

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

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

with

φ^h(z)=(zw(z)I+Ah)−1Phv

φ^h-φ^

Lemma 3.1

Let v ∈ L²(Ω), z ∈ Σ_θ with θ ϵ(π/2, π),

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

φ^h(z)=(zw(z)I+Ah)−1Phv. Then there holds

‖φ^(z)−φ^h(z)‖+h‖∇(φ^(z)−φ^h(z))‖≤ch2‖v‖.

Now we can state an error estimate for nonsmooth data v ∈ L²(Ω).

Next we turn to the case of smooth data Av ∈ L²(Ω).

Proof

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

eh(t)=12πi∫Γθ,δeztw(z)((zw(z)I+A)−1−(zw(z)I+Ah)−1Rh)v dz.

Using w(z)(zw(z) I + A )⁻¹ = z⁻¹I − z⁻¹(zw(z) I + A )⁻¹A, we deduce

eh(t)=12πi(∫Γθ,δeztz−1(φ^(z)−φ^(z))dz+∫Γθ,δeztz−1(v−Rhv)dz),

where

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

φ^h(z)=(zw(z)I+Ah)−1AhRhv.

Then Lemma 3.1 and the identity A_hR_h = P_hA give

‖φ^(z)−φ^(z)‖+h‖∇(φ^(z)−φ^h(z))‖≤ch2‖Av‖.

Now it follows that

‖eh(t)‖≤ch2‖Av‖(∫t−1∞ert cos θr−1dr+∫−θθecos ψdψ)≤ch2‖Av‖,

which gives the L²-estimate. The H¹ estimate is analogous. □

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) + A_h)⁻¹w(z)v, cf., Lemmas 2.1 and 4.3. With z = x + iy, the contour (4.1) is the left branch of the hyperbola

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

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

Now we describe the quadrature approximation. By setting z_j = z(ξ_j) and

and the truncated quadrature approximation

with φ̂_j = (z_jw(z_j)I + A_h)⁻¹w(z_j)v_h. To compute U_N,h(t), we need to solve only N+1 elliptic problems, instead of 2N+1, by exploiting conjugacy:

z′j=z′(ξj)=iλcos(iξj−ψ)

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 I_k, for k > 0, be defined by

I=∫−∞∞g(x)dxand

Ik=k∑j=−∞∞g(jk), respectively. Furthermore, assume that g(z) → 0 uniformly as |z| → ∞ in 𝓢_a,b, and that there exist M₊,M₋ > 0 such that

limr→a− ∫−∞∞|g(x+ir)|dx   ≤ M+      and     lims→b−∫−∞∞|g(x+is)|dx   ≤ M−.

Then with E⁺ = M₊/(e^2aπ/k − 1) and E⁻ = M₋/(e^2bπ/k − 1), there holds

|I−Ik|≤E++E−.

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.

Proof

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(ψ + η)))²−(y/(λcos(ψ +η )))² = 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 ≤ 1 − ϵ²/3 for small ϵ, we have for all ξ ∈ ℝ

|z′(p)z(p)|2= |cos(φ)cosh(ξ)−isin(φ)sin(ξ)(1−cosh(ξ)sin(φ))+isinh(ξ)cos(φ)|2=cosh2(ξ)−sin2(φ)(cosh(ξ)−sin(φ))2=cosh(ξ)+sin(φ)cosh(ξ)−sin(φ)≤1+sin(φ)1−sin(φ)≤21−(1−ϵ2/3)=6ϵ−2

Hence (4.5) holds true. For

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

|z′(p)|2|z(p)|−2≤(1+sin(φ))/(1−sin(φ))≤(1+sin(ψ))/(1−sin(ψ)).

Then the desired result (4.6) follows directly. □

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

Proof

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

z(p)=λ(1−sin(ψ+η)cosh(ξ))+iλcos(ψ+η)sinh(ξ).

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

g^(p,t)=(2πi)−1ez(p)t(z(p)w(z(p))I+Ah)−1w(z(p))z′(p)vh

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

‖(zI+Ah)−1‖≤c/|ℑ(z)|≤c/|zsin(π−ϵ)|≤c/(|z|ϵ)   ∀z ∈Σ′π−ϵ.

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

‖g^(p,t)‖≤ceℜ(z(p))t|z′(p)w(z(p))|‖(z(p)w(z(p))+A)−1‖‖vh‖≤cϵ−1eλ(1−sin(ψ+η)cosh(ξ))t|z′(p)||z(p)|−1‖vh‖.

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 U_N,h.

Proof

Let u_h − U_N,h = (u_h − U_h) + (U_h − U_N,h) =: E_q + E_t. Let a = π/2 − ψ − ϵ and b = ψ − ϵ. For p = ξ + ia, zw(z) lies in Σ_θ for some θ ∈ (π/2, π). Since cosh ξ ≥ 1 + ξ²/2 and 1 − sin(π/2 − ϵ) ≤ ϵ for small ϵ > 0, the choice λ = c₁N/t and Lemma 4.3 yield

‖∫−∞∞|g^(ξ+ia)|dξ‖≤cϵ−1∫−∞∞ec1N(1−sin(π/2−ϵ)cosh(ξ))dξ‖vh‖≤cϵ−1ec1Nε∫0∞e−c1Nsin(π/2−ϵ)ξ2/2dξ‖vh‖≤cϵ−1N−1/2ec1Nϵ‖vh‖.

By Lemma 4.1, for k = c₀/N we have

‖Eq+‖≤cϵ−1N−12e−(2π(π2−ψ−ϵ)/c0−ϵc1)N‖vh‖.

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

‖∫−∞∞|g^(ξ−ib)|dξ‖≤cϵ−1∫0∞ec1N(1−sin(ϵ)cosh(ϵ))dξ‖vh‖≤cϵ−1ec1N(1−sin(ϵ))∫0∞e−c1Nsin(ϵ)ξ2/2dξ‖vh‖≤cϵ−3/2N−1/2ec1N(1−ϵ)‖vh‖.

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

‖Eq−‖≤cϵ−3/2N−1/2e−(2π(ψ−ϵ)/c0−c1(1−ϵ))N‖vh‖. Further, since cosh(ξ) ≥ cosh(c₀)+sinh(c₀)(ξ − c₀) for ξ ≥ c₀, the error ║E_t║ is bounded by

‖Et‖≤cϵ−1∫c0∞ec1N(1−sin(ψ)cosh(ξ))dξ‖vh‖≤cϵ−1ec1N(1−sin(ψ)cosh(c0))∫c0∞e−c1Nsin(ψ)sinh(c0)(ξ−c0)dξ‖vh‖=c/(c1sin(ψ)sinh(c0)ϵ)N−1ec1N[1−sin(ψ)cosh(c0)]‖vh‖.

By disregarding terms ε, balancing the exponentials in

||Eq+||,||Eq−|| and ║E_t║, we get 2π(π/2−ψ)/c₀ = 2πψ/c₀−c₁ = −c₁(1−sin(ψ)cosh(c₀)). Next we express c₀ and c₁ using ψ : c₀ = cosh⁻¹(2πψ/((4πψ − π²)sin ψ)) and c₁ = (4πψ − π²)/ cosh⁻¹(2πψ/((4πψ − π²)sin ψ)). Finally we minimize the ratio B(ψ) = c₁⁻²πψ/c₀ with respect to ψ, which achieves the minimum at ψ = 1.1721 and hence, c₀ = 1.0818, c₁ = 4.4920 and B(ψ) = − 2.32, which are identical to the values in [32]. Then collecting the balanced asymptotic bound and

||Eq+||,||Eq−||, and choosing ϵ = 1/N yield

‖uh(t)−UN,h(t)‖≤cϵ−1(N−12+ϵ−12N−12+N−1)e−[2.32−(2πc0+c1)ϵ]N‖vh‖≤ce(−2.32+ln NN)N‖vh‖.

Since x⁻¹ ln x ≤ 1/e for x ≥ 1, the L²-stability of P_h 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.

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 t_n = nτ, n = 0, ..., N. Following [19, 5], we consider the convolution quadrature generated by the backward Euler method. The weights

ω˜(ξ)=∑j=0∞bjξj=(1−ξτ)w(1−ξτ)

Qn(φ)=∑j=0nbn−jφ(tj)

Uh0 = vh

We denote the generating function

{βj}j=0∞

β˜(ξ)=∑j=0∞βjξj

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