A Galerkin finite element method for the modified distributed-order anomalous sub-diffusion equation

Journal Pre-proof A Galerkin finite element method for the modified distributed-order anomalous sub-diffusion equation Lang Li, Fawang Liu, Libo Feng, Ian Turner

PII: DOI: Reference:

S0377-0427(19)30594-1 https://doi.org/10.1016/j.cam.2019.112589 CAM 112589

To appear in:

Journal of Computational and Applied Mathematics

Received date : 10 January 2019 Revised date : 7 October 2019 Please cite this article as: L. Li, F. Liu, L. Feng et al., A Galerkin finite element method for the modified distributed-order anomalous sub-diffusion equation, Journal of Computational and Applied Mathematics (2019), doi: https://doi.org/10.1016/j.cam.2019.112589. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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.

© 2019 Elsevier B.V. All rights reserved.

Highlights (for review)

Journal Pre-proof Highlights

Jo

urn

al

Pr e-

p ro

of

1. A modified time distributed-order anomalous sub-diffusion equation is considered. 2. We transform the time distributed-order fractional diffusion equation into a modified multi-term anomalous sub-diffusion equation.. 3. A novel Galerkin finite element method is proposed. 4. The stability and convergence of the proposed method are discussed. 5. Numerical examples are given to verify the accuracy and stability of our scheme.

Manuscript Click here to view linked References

Journal Pre-proof

A Galerkin finite element method for the modified distributed-order anomalous sub-diffusion equation Lang Lia , Fawang Liub,c,∗, Libo Fengb , Ian Turnerb,d

of

a Department of Mathematics, South China Agricultural University, Guangzhou, 510642, China of Mathematical Sciences, Queensland University of Technology, GPO Box 2434, Brisbane, Qld. 4001, Australia c College of Mathematics and Computer Science, Fuzhou University, Fuzhou 350116, China d Australian Research Council Centre of Excellence for Mathematical and Statistical Frontiers (ACEMS), Queensland University of Technology, Brisbane, QLD 4001, Australia

p ro

b School

Abstract

Pr e-

In this paper, we consider a modified time distributed-order anomalous sub-diffusion equation for the description of processes becoming less anomalous over the course of time. First, we discretize the integral term using a mid-point quadrature rule and transform the time distributed-order fractional diffusion equation into a multi-term equation. Then applying the backward difference method in time and Galerkin finite element method in space, we present a fully discrete scheme to solve the time distributed-order fractional diffusion equation. Moreover, the stability and convergence of the proposed method are also discussed. Numerical examples are given to verify the accuracy and stability of our scheme. Keywords: The distributed-order anomalous sub-diffusion equation, Galerkin finite element method, Error analysis, Stability and Convergence

1. Introduction

urn

al

In the past two decades, much attention and effort have been given to fractional partial differential equations containing different fractional operators because of their ability to model anomalous phenomena in many scientific and engineering fields. Anomalous diffusion phenomena occuring in physics and other scientific applications can be modeled by fractional diffusion equations [1, 2, 3, 4, 5, 6]. Anomalous diffusion is characterized by the mean square displacement hr2 (t)i and a scaling parameter γ as well as a diffusion parameter Kγ as a nonlinear power-law in time, i.e. hr2 (t)i ∼ Kγ tγ . When 0 < γ < 1, we call it anomalous sub-diffusion. Experts have revealed sub-diffusion of proteins and lipids in a variety of cell membranes by single particle tracking experiments and photo-bleaching recovery experiments [7, 8, 9, 10]. It is noted that fractional kinetic equations played an important role in the description of anomalous sub-diffusion. In [11], an explicit finite difference method was proposed to solve the following fractional sub-diffusion equation ∂2 ∂ u(x, t) = K 0 Dt1−γ 2 u(x, t), t > 0, ∂t ∂x

(1)

Jo

where u(x, t) is the probability density function, the operator 0 Dt1−γ is the Riemann-Liouville derivative of order 1 − γ defined by 1−γ u(x, t) 0 Dt

1 ∂ = Γ(γ) ∂t

Z t

0

u(x, s) ds. (t − s)1−γ

There are also other papers that have considered this type of equation [12, 13, 14], etc. In [15, 16], a modified model (2) which included a second fractional time derivative acting on the diffusion term was proposed and studied for ∗ Corresponding

author Email address: [email protected] (Fawang Liu)

Preprint submitted to Journal of Computational and Applied Mathematics

October 7, 2019

Journal Pre-proof

describing some physical process that becomes less anomalous over the course of time. € Š ∂2 ∂ u(x, t), t > 0, u(x, t) = A 0 Dt1−α + B 0 Dt1−β ∂t ∂x2

(2)

∂u(x, t) = ∂t

Z 1

0



ω(γ)0 Dt1−γ

K

p ro

of

where 0 < α < β ≤ 1, A and B are constants. However, many physical processes lacking power-law scaling over the whole time domain cannot be characterized by a single or two scaling exponents. In [17, 18], it was noted that processes lacking scaling such as decelerating sub-diffusion (e.g. the Sinai diffusion [19]) and decelerating super-diffusion (truncated L´evy-flights [20]) can be described by the distributed-order fractional diffusion equation. Moreover, as it was shown in [21], the modified equation with the distributed-order derivatives on the right hand side can describe processes getting less anomalous over the course of time (decelerating super-diffusion and accelerating sub-diffusion). In this paper, we consider the following distributed-order anomalous sub-diffusion equation with initial and boundary conditions: ∂ 2 u(x, t) ∂x2

u(x, 0) = u0 (x), x ∈ Ω,

dγ + f (x, t), (x, t) ∈ Ω × (0, T ],

Pr e-

u(0, t) = 0, u(L, t) = 0, 0 ≤ t ≤ T,



(3) (4) (5)

where Ω = (0, L), and ω(γ) satisfies the following conditions

Z 1

0 ≤ ω(γ), ω(γ) 6≡ 0, γ ∈ [0, 1],

ω(γ)dγ = W > 0.

0

Jo

urn

al

The idea of distributed-order fractional differential equations, which can be viewed as a generalization of the single and the multi-term fractional differential equations, was first proposed by Caputo who applied them to infinite plates and dielectric spherical shells and to find the eigenfunctions of the torsional models of anelastic [22]. Recently, the research on this class of equations has become increasingly popular in many physical and biological applications, for example, distributed-order fractional Langevin-like equations [23], distributed-order membranes in the ear [24], viscoelastic oscillators [25], radial groundwater flow to or from a well [26], control and signal processing [27], dielectric induction and diffusion [28]. There are some related results about the analytical solutions of distributed-order fractional differential equations. For example, Mainardi et al. [29] obtained the Fourier-Laplace representation of the fundamental solutions to the time distributed diffusion equations. Gorenflo et al. [30] used Fourier and Laplace transforms to derive the representation of the fundamental solution for the Cauchy problem of a distributed-order time fractional diffusion wave equation. Uniqueness and existence of solutions to the generalized time fractional diffusion equation of distributed-order were obtained by using an approximate maximum principle [31]. However, for more general distributed-order fractional equations, it is not easy to obtain an analytical solution. Therefore, efficient numerical methods need to be developed to solve the distributed-order fractional differential equations. It appears that only a few papers have been devoted to the analysis of the numerical methods for solving partial differential equations with a distributed-order. Here we briefly review some of the recent advances on this topic. In [32], Ye et al. derived and analysed a compact difference scheme for a distributed-order time-fractional diffusionwave equation, where the order of the time fractional derivatives was distributed in the interval [1,2], and the Caputo time fractional derivatives were approximated by the L1 formula [33]. For the integral term in the distributed-order equation, it was approximated by the mid-point quadrature rule. In [34, 35], Gao et al. proposed some high order difference schemes for both one dimensional and two dimensional time distributed-order differential equations. They first used the composite Simpson quadrature formula or the composite trapezoid quadrature formula to discretize the distributed integral term and then approximated the fractional derivatives by the weighted and shifted Gr¨ unwaldLetnikov formula. For the space distributed-order fractional differential equations, Li et al. [36, 37] proposed a novel finite volume method for the Riesz space distributed-order advection-diffusion or diffusion equation. In [38], Zhang et al. presented a Crank-Nicolson alternating direction implicit Galerkin-Legendre spectral scheme to solve the two dimensional Riesz space distributed-order advection-diffusion equation. The Gauss quadrature formula was first used to discretize the integral term and then the time derivative was approximated by the Crank-Nicolson method, while the spatial derivatives were approximated by the ADI Galerkin-Legendre spectral scheme. Fan et al. 2

Journal Pre-proof

p ro

of

[39] proposed a numerical method to solve the two dimensional distributed-order space-fractional diffusion equation on an irregular domain. The distributed-order anomalous sub-diffusion equation (3)-(5) can be viewed as an extension to the anomalous sub-diffusion equation. To the best of our knowledge, no published paper takes into account the numerical method to this type of equation, which is important for their potential applications. The motivation of our work is to fill this gap. The main attention in this paper is focused on the numerical analysis for the distributed-order anomalous sub-diffusion equation, where the one dimensional problem in space and the order of fractional derivatives in time is distributed in a range [0,1]. First, we discretize the integral term in the time distributed-order anomalous sub-diffusion equation by the mid-point quadrature rule and then we transform it into multi-term fractional subdiffusion equation. Furthermore, we derive a Galerkin finite element scheme for space. Using the new energy method, we prove the stability and convergence of the solution. Finally, we give some numerical examples to show the effectiveness of our method and to explain the physics associated with the distributed-order sub-diffusion model. The rest of this paper is organized as follows. In Section 2, we derive a semi-discrete finite difference scheme and variation formulation. The stability and error analysis are shown by using a new energy method. In Section 3, we present the Galerkin finite element approximation for space and obtain the fully discrete Galerkin finite element scheme. The stability and convergence analysis are verified. In Section 4, we give some numerical examples.

Pr e-

2. Variation formulation and a semi-discrete finite difference scheme

First, we discrete the interval [0, 1], in which the order γ is changing, using the grid 0 = ξ0 < ξ1 < · · · < ξq = 1 (q ∈ N ), with the steps ∆ξs not necessarily equidistant. We obtain Z 1

0



ω(γ)0 Dt1−γ

∂ 2 u(x, t) ∂x2

K



q X

dγ ≈

quadrature rule to approximate the above integral. Let γs = ∂ 2 u(x, t) ∂x2



q X

al

0



ω(γ)0 Dt1−γ K

K

s=1

where γs ∈ (ξs−1 , ξs ], ∆ξs = ξs − ξs−1 , s = 1, 2, · · · , q. For simplicity, without loss of generality, we take ∆ξs =

Z 1



ω(γs )0 Dt1−γs

dγ =

1 q

= σ and ds =

ξs−1 +ξs 2



ds

∂ 2 u(x, t) ∂x2

=

2s−1 2q ,



1−γs 0 Dt

K

s=1



∆ξs ,

ω(γs ) q .

We can use the mid-point

s = 1, 2, · · · , q. Then we have

∂ 2 u(x, t) ∂x2



+ R1 ,

urn

where R1 = O(σ 2 ). Next we consider the following multi-term fractional sub-diffusion equation q

∂u(x, t) X = ds ∂t s=1





1−γs 0 Dt

∂ 2 u(x, t) K ∂x2



+ f (x, t) + R1 ,

(6)

Jo

with initial and boundary conditions



u(x, 0) = u0 (x), x ∈ Ω,

(7)

u(0, t) = 0, u(L, t) = 0, 0 ≤ t ≤ T.

(8)

In order to approximate 0 Dt1−γs K ∂

2

u(x,t) ∂x2



, we discretize the temporal domain [0, T ] by using a uniform grid.

Define tk = kτ , k = 0, 1, 2, . . . , n, where τ = T /n. Let 1 I y(t) = Γ(γ) γ

Z t

0

y(s) ds, t > 0, (t − s)1−γ

where 0 < γ < 1. Then we have the following lemmas.

3

Journal Pre-proof

Lemma 1. If y(t) ∈ C 1 [0, T ], then k−1

X τγ (1) (γ) b y(tk−j ) + rk,γ , Γ(γ + 1) j=0 j

I γ y(tk ) = (γ)

where bj

(1)

= (j + 1)γ − j γ (j = 0, 1, 2, ..., n) and |rk,γ | ≤ Ctγk τ, k = 1, 2, ..., n.

=

1 Γ(γ)

0

Z t k−1 j+1 X

j=0

k−1 P Rt j+1 y(s)−y(tj+1 ) tj (tk −s)1−γ ds. j=0 k−1 Z tj+1

1 X Γ(γ) j=0

tj

k−1

1 X Γ(γ) j=0

k−1

tj

k−1 Z tj+1

X y(s) − y(tj+1 ) 1 ds| ≤ max |y ′ (t)|τ 1−γ (tk − s) Γ(γ) 0≤t≤tk j=0

k−1

≤C

y(tj+1 ) (1) ds + rk,γ , (tk − s)1−γ

X X τγ y(tj+1 ) τγ (γ) (γ) bk−j−1 y(tj+1 ) = b y(tk−j ), ds = 1−γ (tk − s) Γ(γ + 1) j=0 Γ(γ + 1) j=0 j

k−1 Z tj+1

(1)

tj

tj

y(s) ds (tk − s)1−γ

After direct calculation, we have

and |rk,γ | =|

k−1 Z tj+1

y(s) 1 X ds = (tk − s)1−γ Γ(γ) j=0

Pr e-

(1)

where rk,γ =

1 Γ(γ)

Z t k

p ro

1 I y(tk ) = Γ(γ) γ

of

Proof. Note that

tj

1 ds (tk − s)1−γ

X 1 (γ) τ b τ γ ≤ Ctγk τ. Γ(γ + 1) j=0 k−j−1

Thus we complete the proof.

(1)

al

Since tk ≤ T and 0 ≤ γ ≤ 1, it follows that |rk,γ | ≤ Ctγk τ ≤ Cτ. (γ)

Lemma 2. ([40]) The coefficients bj (j = 0, 1, 2, ..., n) satisfy the following properties: = 1, bj

(γ)

(γ) bj

(γ) bj+1 ,

> 0, j = 0, 1, 2, . . . ;

urn

(γ)

1. b0

j = 0, 1, 2, . . . ; > 2. 3. there exists a positive constant C > 0 such that (γ)

τ ≤ Cbj τ γ , j = 1, 2, . . .

In the following, C is a positive constant which may change from line to line. Let ∆t y(tk ) = y(tk+1 ) − y(tk ), then we have the following lemma.

Jo

Lemma 3. For k = 1, 2, . . . , n,

∆t I γ y(tk−1 ) =

1 Γ(γ)

Z t k

0

1 y(s) ds − (tk − s)1−γ Γ(γ)

"

k−1

Z t k−1

0

y(s) ds (tk−1 − s)1−γ #

X τγ (2) (γ) (γ) = y(tk ) + (bj − bj−1 )y(tk−j ) + rk,γ , Γ(γ + 1) j=1

and "

1−γ y(tk ) 0 Dt

(2)

(γ)

(3)

k−1

#

X τ γ−1 (3) (γ) (γ) = y(tk ) + (bj − bj−1 )y(tk−j ) + rk,γ , Γ(γ + 1) j=1

(γ)

where |rk,γ | ≤ Cbk−1 τ γ+1 , |rk,γ | ≤ Cbk−1 τ γ . 4

Journal Pre-proof

Proof. Note that Z

Z

tk tk−1 1 y(s) y(s) 1 ds − ds 1−γ Γ(γ) 0 (tk − s) Γ(γ) 0 (tk−1 − s)1−γ Z τ  Z t k−1 y(s + τ ) − y(s) y(s) 1 = ds + ds Γ(γ) 0 (tk − s)1−γ (tk − s)1−γ 0

∆t I γ y(tk−1 ) =

of

= I1 + I2 . For I1 , we have 1 Γ(γ)

Z τ

0

Z

Assume that y(t) ∈ C 2 [0, T ], we obtain 1 Γ(γ)

Z τ

0



For I2 , applying Lemma 1, we have 1 I2 = Γ(γ)

Z t k−1

0

(γ)

y(s) − y(τ ) τ 1+γ bk−1 (γ) ds ≤ max |y ′ (t)| ≤ Cτ 1+γ bk−1 . (tk − s)1−γ Γ(γ + 1) 0≤t≤τ

Pr e-



Z

τ τ y(s) y(τ ) y(s) − y(τ ) 1 1 ds = ds + ds 1−γ 1−γ (tk − s) Γ(γ) 0 (tk − s) Γ(γ) 0 (tk − s)1−γ Z τ 1 τγ y(s) − y(τ ) (γ) ds. bk−1 y(τ ) + = Γ(γ + 1) Γ(γ) 0 (tk − s)1−γ

p ro

I1 =

k−1

X y(s + τ ) − y(s) τγ (4) (γ) ds = [y(tk+1−j ) − y(tk−j )] + rk,γ , b (tk − s)1−γ Γ(γ + 1) j=1 j−1

where (4)

τ tγ max |y ′ (t + τ ) − y ′ (t)| Γ(γ + 1) k−1 0≤t≤tk−1 τ = tγ max |y ′′ (ξ)τ | Γ(γ + 1) k−1 0≤ξ≤tk−1

|rk,γ | ≤

≤ Cτ 2 .

al

Hence, we have

"

#

k−1

X τγ (γ) (γ) (2) bj−1 (y(tk+1−j ) − y(tk−j )) + rk,γ bk−1 y(τ ) + ∆t I y(tk−1 ) = Γ(γ + 1) j=1

γ

"

#

k−1

urn

X τγ (γ) (γ) (2) = (bj − bj−1 )y(tk−j ) + rk,γ , y(tk ) + Γ(γ + 1) j=1

where

(2)

(γ)

rk,γ ≤ C(bk−1 τ γ+1 + τ 2 ).

(γ)

Jo

Using Lemma 2, we have τ 2 ≤ Cbk−1 τ γ+1 . Then we have Since 0 Dt1−γ y(tk ) =

∂ γ ∂t I y(tk ),

(2)

(γ)

rk,γ ≤ Cbk−1 τ γ+1 .

we obtain 1−γ y(tk ) 0 Dt

=

1 (5) ∆t I γ y(tk−1 ) + rk,γ , τ

(5)

where |rk,γ | ≤ Cτ. Applying Lemma 2, we have 1−γ y(tk ) 0 Dt (3)

(γ)

"

k−1

#

X τ γ−1 (γ) (γ) (3) y(tk ) + (bj − bj−1 )y(tk−j ) + rk,γ , = Γ(γ + 1) j=1

(γ)

where |rk,γ | ≤ Cbk−1 τ γ + Cτ ≤ Cbk−1 τ γ . 5

Journal Pre-proof

Without loss of generality, we assume f (x, t) = 0 in the construction of semi-discrete scheme and its numerical analysis. As a result of Lemma 3, from (6), we have "

q

#

k−1

q

2 X u(x, tk ) − u(x, tk−1 ) X Kτ γs −1 ∂ 2 u(x, tk ) X (γs ) (γs ) ∂ u(x, tk−j ) + = ) − b (b + ds ds Rsk + R1 , j−1 j τ Γ(γs + 1) ∂x2 ∂x2 s=1 s=1 j=1

(9)

(γ )

uk = uk−1 +

q X

"

Ks

#

k−1

2 k−j ∂ 2 uk X (γs ) (γs ) ∂ u , + ) − b (b j−1 j ∂x2 ∂x2 j=1

p ro

s=1

of

s where |Rsk | ≤ Cbk−1 τ γs . Let uk (x) be a numerical approximation to u(x, tk ). We obtain the following fractional implicit difference approximation

γs

(10)

Kτ where k = 0, 1, · · · , n − 1, Ks = ds Γ(γ . The variation formulation of equation (10) subjects to the boundary s +1) conditions

uk (0) = uk (L) = 0 reads as follows: find uk ∈ H01 (Ω) such that k

Š

€

u ,ϕ = u

k−1

Š

,ϕ −

"



Pr e-

€

q X

Ks

s=1

∂uk ∂ϕ , ∂x ∂x

+

k−1 X€

(γ ) bj s

j=1

−

(γs ) bj−1

Š



∂uk−j ∂ϕ , ∂x ∂x

#

,

holds for ∀ ϕ ∈ H01 (Ω). In the following, we will discuss the stability and error estimates for the semi-discrete problem. 2.1. Stability analysis for the semi-discrete scheme Now we state one of our main results in Theorem 1.

al

Theorem 1. The semi-discrete difference scheme is unconditionally stable in the sense for all τ , the following inequalities hold Ek ≤ Ek−1 ≤ · · · ≤ ||u0 ||2 , k = 1, 2, · · · , n,

Proof. For k = 1,

q P s=1

Ks

k−1 P (γs ) k−j bj || ∂u∂x j=0

||2 .

urn

where Ek = ||uk ||2 +

(u1 , ϕ) = (u0 , ϕ) −

q X



Ks

s=1

∂u1 ∂ϕ , ∂x ∂x



.

Jo

Taking ϕ = u1 , we have

||u1 ||2 = (u0 , u1 ) −

q X s=1

Ks ||

∂u1 2 || . ∂x

Using Cauchy Schwarz inequality (u, w) ≤ 12 (||u||2 + ||w||2 ), we obtain q

||u1 ||2 ≤

X 1 ∂u1 2 (||u0 ||2 + ||u1 ||2 ) − || . Ks || 2 ∂x s=1

Hence ||u1 ||2 + 2

q X s=1

Ks ||

6

∂u1 2 || ≤ ||u0 ||2 . ∂x

Journal Pre-proof

Furthermore, q X

||u1 ||2 +

s=1

Ks ||

∂u1 2 || ≤ ||u0 ||2 . ∂x

i.e., E1 ≤ ||u0 ||2 . For k > 1, taking ϕ = uk , the following formula is obtained. q X

k

,u ) −

s=1

"

Ks

k−1

∂uk 2 X (γs ) (γs ) || + ) (bj − bj−1 || ∂x j=1

Using Cauchy Schwarz inequality, we have q

||uk ||2 ≤

q

X 1 ∂uk 2 X Ks (||uk−1 ||2 + ||uk ||2 ) − || + Ks || 2 ∂x 2 s=1 s=1

(γs )

(γ )

s > bj where we have used the fact bj−1

k−1 X j=1



q

Pr e-

q X

q X

1 k−1 2 Ks ∂uk 2 Ks ||u || − || || + 2 2 ∂x 2 s=1 s=1

Furthermore, we have q X

Ks

j=0

||

Ek = ||uk ||2 +

Thus

.



) ||

∂uk 2 ∂uk−j 2 || + || || , ∂x ∂x

j=1

j=1

(γ )

(γs )

(γ )

s − bj (bj−1

(γs )

s − bj (bj−1

)||

)||

∂uk−j 2 || ∂x

∂uk−j 2 || . ∂x

k−1

k−j X X ∂uk−j 2 (γs ) ∂u || ≤ ||uk−1 ||2 + ||2 || bj−1 Ks ∂x ∂x s=1 j=1

urn

Let

k−1 X

k−1 X

q

(γs )

bj

al

s=1

k−1 X

(γs )

(γ )

#

for j = 1, 2, · · · . After simplification, we have

q

||uk ||2 +

∂uk−j ∂uk , ∂x ∂x

s − bj (bj−1

X K 1 ∂uk 2 X Ks 1 k 2 s (γs ) ||u || ≤ ||uk−1 ||2 − (1 + bk−1 )|| || + 2 2 2 ∂x 2 s=1 s=1

≤



of

||u || = (u

k−1

p ro

k 2

≤ ||uk−1 ||2 +

q X

Ks

s=1

k−1 X

(γs )

bj

j=0

q X s=1

||

Ks

k−2 X j=0

(γs )

bj

||

∂uk−1−j 2 || . ∂x

∂uk−j 2 || . ∂x

Ek ≤ Ek−1 .

Jo

By induction, the proof is complete.

2.2. Error estimates for the semi-discrete scheme Before we give the error analysis for the solution of the semi-discrete difference scheme, we introduce the following lemma. Lemma 4. (Poincare-Friedrich’s inequality) If y(x) ∈ H01 (0, L), then ||y(x)||2 ≤

L2 dy 2 || || . 2 dx

Let u(x, t) and {uk (x)}nk=0 be the exact solution and time-discrete solution respectively. Denote η k = u(x, tk ) − u , we have the following theorem. k

7

Journal Pre-proof

Theorem 2. There exists a positive constant C˜ > 0 such that k 2

Yk = ||η || +

q X

Ks

s=1

k−1 X

(γs )

bj

j=0

||

∂η k−j 2 ˜ + σ 2 )2 , || ≤ C(τ ∂x

consequently,

of

È

˜ + σ 2 ), k = 1, 2, · · · , n. C(τ

||η k || ≤ Proof. From (9)-(10), we obtain "

Ks

s=1

∂η k ∂ϕ , ∂x ∂x



+

k−1 X



(γ ) (bj s

j=1

−

∂η k−j ∂ϕ , ∂x ∂x

#

p ro

(η k , ϕ) = (η k−1 , ϕ) −

q X

(γs ) ) bj−1

+ (τ

q X

ds Rsk , ϕ) + (τ R1 , ϕ).

s=1

Taking ϕ = η k , and then using Cauchy-Schwarz inequality, we have ||η || = (η ≤

k−1

k

,η ) −

q X

"

k−1

∂η k 2 X (γs ) (γs ) || || + ) (bj − bj−1 ∂x j=1

Ks

s=1

q X

q X

1 ∂η k 2 Ks (||η k ||2 + ||η k−1 ||2 ) − || + Ks || 2 ∂x 2 s=1 s=1

+ |(τ

q X s=1



∂η k−j ∂η k , ∂x ∂x

Pr e-

k 2

k−1 X j=1

(γ )

#

+ (τ



(γs )

s − bj (bj−1

) ||

q X

ds Rsk , η k ) + (τ R1 , η k )

s=1

∂η k 2 ∂η k−j 2 || + || || ∂x ∂x



ds Rsk , η k )| + |(τ R1 , η k )|.

Let

k 2

Yk = ||η || +

Ks

s=1

k−1 X

(γs )

bj

j=0

al

then we have

q X

q X s=1

(γ )

s Ks bk−1 ||

∂η k−j 2 || , ∂x

q

X ∂η k 2 ds Rsk , η k )| + 2|(τ R1 , η k )|. || + 2|(τ ∂x s=1

urn

Yk ≤ Yk−1 −

||

Using Lemma 4, we can deduce that

≤ Yk−1 − ≤ Yk−1 +

q X s=1 q X s=1

≤ Y0 + C

q

(γ )

q

(γ )

q

q

(γ )

s X Ks b s X L2 d2 X Ks b s 2Ks bk−1 k−1 k−1 s k 2 k 2 2 k 2 ||η || + ||η || + ||η k ||2 + τ (R ) + s 2 2 (γs ) L2 L L K b s s=1 s=1 s=1

q X

L2 d2s 2 k 2 τ (Rs ) (γs ) s=1 Ks bk−1

≤ Yk−1 + C ≤ ···

(γ )

s X 2Ks bk−1 ds Rsk , η k )| + 2|(τ R1 , η k )| ||η k ||2 + 2|(τ 2 L s=1

Jo

Yk ≤ Yk−1 −

q X

+

1 (γs ) Pq Ks bk−1 s=1 L2

s=1

j=1 s=1

τ 2 |R1 |2 1

(γ )

s ds Γ(γs + 1)τ 2 bk−1 τ γs + C

q k X X

k−1

(γ )

Pq

s τ γs + C ds Γ(γs + 1)τ 2 bj−1

(γ )

s Ks bk−1 s=1 L2

k X

τ 2 σ4

1

(γs ) Pq Ks bj−1 j=1 s=1 L2

8

τ 2 σ4 .

1 (γs ) Pq Ks bk−1 s=1 L2

τ 2 |R1 |2

Journal Pre-proof

Since Y0 = 0, it follows

≤C (γs )

q X s=1

q L2 X

K

ds Γ(γs + 1)τ 2

k X

(γ )

s τ γs + C bj−1

j=1

k X j=1

ds Γ(γs + 1)τ 2 (kτ )γs + Ck

s=1

L2 K

q X

≤C

L2 K

q X

(γ )

s Ks bj−1 s=1 L2

1 (γs ) Pq Ks bk−1 s=1 L2

. Thus

Yk ≤ C

1 Pq

τ 2 σ4

τ 2 σ4 ,

1

p ro

(γ )

s ≥ bj owing to the fact bj−1

L2 K

of

Yk ≤ C

ds Γ(γs + 1)τ 2 T γs + CT

s=1

(γs ) Pq ds Kτ γs bk−1 s=1 Γ(γs +1)L2

ds Γ(γs + 1)τ 2 T γs + CT Pq

1

ds Kτ s=1 Γ(γs +1)L2

s=1

τ σ4

τ σ4

(γ )

Pr e-

T L2 4 L2 ˜ 2 + σ 4 ) ≤ C(τ ˜ + σ 2 )2 , σ ≤ C(τ ≤ C C1 τ 2 + C K KC2 s , where we have used the facts τ ≤ Cτ γs bk−1

q X s=1

ds Γ(γs + 1)T γs →

and q X

Consequently,

ω(γ)Γ(γ + 1)T γ dγ = C1 ,

0

Z 1

ds /Γ(γs + 1) →

ω(γ)/Γ(γ + 1)dγ = C2 .

0

al

s=1

Z 1

È

||η k || ≤

urn

Thus we complete the proof.

˜ + σ 2 ), k = 1, 2, · · · , n. C(τ

3. Galerkin finite element approximation in the space and error estimates for full discretization

Jo

In this section, we first introduce some notations and lemmas that are needed in the following. Denote by (·, ·) the inner product defined on the space L2 (Ω) with the norm k · k0 . Define Pl (Ω) as the space of polynomials defined on Ω with degree at most l . Let Vh be a uniform partition of Ω, which is given by 0 = x0 < x1 < x2 < · · · < xm = L,

where m is a positive integer. Let h = L/m and Ωi = [xi−1 , xi ] for i = 1, 2, · · · , m. Denote Sh to be a family of piecewise polynomials on mesh xi , which can be expressed by Sh = {v : v|Ωi ∈ Pl (Ωi ), v ∈ C(Ω)}.

When l = 1, the nodal basis functions φ0 , φ1 , · · · , φm of Sh can be expressed in the form 8

x−xi−1 , h xi+1 −x , h : <

φi (x) =

0,

9

x ∈ [xi−1 , xi ], x ∈ [xi , xi+1 ], elsewhere,

Journal Pre-proof

§ x −x 1 h ,

φ0 (x) =

x ∈ [x0 , x1 ], elsewhere,

0,



x−xm−1 , h

φm (x) =

x ∈ [xm−1 , xm ], elsewhere.

0,

Lemma 5. [41] For i = 1, 2, · · · , m − 1, we have

of

where i = 1, 2, · · · , m − 1 and

8

p ro

1, |j − i| = 1, h< 4, j = i, (φi (x), φj (x)) = 6: 0, elsewhere and ∂φi (x) ∂φj (x) , ∂x ∂x

8

−1, |j − i| = 1, 1< 2, j = i, h: 0, elsewhere.

‹

=

Pr e-



Introduce the Ritz projection Πh : H l+1 (Ω) ∩ H01 (Ω) → Sh by the orthogonal relation (∂x (u − Πh u), ∂x χ) = 0, ∀ χ ∈ Sh , u ∈ H l+1 (Ω) ∩ H01 (Ω).

(11)

We have the well known that approximated property [42]:

||u(x, tk ) − Πh u(x, tk )|| + h||(u(x, tk ) − Πh u(x, tk ))′ || ≤ Chl+1 , k = 0, 1, 2, · · · , n. Let

1 ∂u(x, tk ) (1) = (Πh u(x, tk ) − Πh u(x, tk−1 )) + Rk , ∂t τ

al

(12)

where

(1)

∂u(x, tk ) 1 − (Πh u(x, tk ) − Πh u(x, tk−1 )). ∂t τ

urn

Rk = (1)

We also can rewrite Rk as the following form (1) Rk

Jo

Thus

∂u(x, tk ) = (I − Πh ) + Πh ∂t



‹

∂u(x, tk ) u(x, tk ) − u(x, tk−1 ) − . ∂t τ

(1)

||Rk || ≤ C(τ + hl+1 ).

In the same way, we also have

1−γs u(x, tk ) 0 Dt

"

#

k−1

X τ γs −1 (s) (γs ) (γ ) Πh u(x, tk ) + )Πh u(x, tk−j ) + Rk , (bj s − bj−1 = Γ(γs + 1) j=1

(13)

where (

(s) Rk

= (I −

Πh )0 Dt1−γs u(x, tk )

+ Πh

"

1−γs u(x, tk ) 0 Dt

k−1

#)

X τ γs −1 (γs ) (γ ) − )u(x, tk−j ) u(x, tk ) + (bj s − bj−1 Γ(γs + 1) j=1

10

,

Journal Pre-proof

and (s)

(γ )

s ||Rk || ≤ C(bk−1 τ γs + hl+1 ).

Applying Lemma 3, and using (11)-(13), the weak form of (6)-(8) at tk is given by q

k−1 X

(γs )

(bj

j=1

(γ )

s )(∂x Πh u(x, tk−j ), ∂x ϕ)] + (fk , ϕ) + (Rk , ϕ), ∀ϕ ∈ H01 (Ω), − bj−1

p ro

+

of

X Kτ γs −1 1 [(∂x Πh u(x, tk ), ∂x ϕ) ( (Πh u(x, tk ) − Πh u(x, tk−1 )), ϕ) = − ds τ Γ(γs + 1) s=1

where q X

||Rk || ≤ C

!

(γs ) γs ds bk−1 τ

+ σ 2 + τ + hl+1

s=1

Since Sh ⊂ H01 (Ω), we can choose ϕ = χ ∈ Sh , and obtain q

, k = 1, 2, · · · , n.

+

k−1 X

(γs )

(bj

j=1

Pr e-

X 1 Kτ γs −1 [(∂x Πh u(x, tk ), ∂x χ) ( (Πh u(x, tk ) − Πh u(x, tk−1 )), χ) = − ds τ Γ(γs + 1) s=1

(γ )

s )(∂x Πh u(x, tk−j ), ∂x χ)] + (fk , χ) + (Rk , χ), ∀χ ∈ Sh , k = 1, 2, · · · , n. − bj−1

Thus, we have

(Πh u(x, tk ), χ) = (Πh u(x, tk−1 ), χ) − k−1 X

(γs )

(bj

Ks [(∂x Πh u(x, tk ), ∂x χ)

s=1

(γ )

s )(∂x Πh u(x, tk−j ), ∂x χ)] + (τ fk , χ) + (τ Rk , χ), ∀χ ∈ Sh , k = 1, 2, · · · , n. − bj−1

al

+

q X

j=1

(14)

Denote by U k ∈ Sh , the approximation of u(·, tk ), then the fully discrete finite element method is given by q X

"

urn

k

(U , χ) = (U In particular, let U k =

k−1

Pm j=0

, χ) −

k

Ks (∂x U , ∂x χ) +

#

k−1 X

s=1

(γ ) (bj s

j=1

−

(γs ) )(∂x U k−j , ∂x χ) bj−1

+ (τ fk , χ).

ukj φj (x) and χ = φi (x)(i = 1, 2, · · · , m − 1). Using Lemma 5, we have q

h k−1 1X h k k−1 Ks [uki−1 − 2uki + uki+1 ] (ui−1 + 4uki + uki+1 ) = (ui−1 + 4uik−1 + ui+1 )+ 6 6 h s=1

Denote

"

Jo

q

1X + Ks h s=1 0

A=

h 6

B B B B B B B 

k−1 X

(γ ) (bj s

j=1

4 1 1 4 0 1 .. .. . . 0 0 0 0

−

#

(γs ) k−j )(ui−1 bj−1

0 1 4 .. .

··· ··· ··· .. .

0 0 0 .. .

0 0 0 .. .

0 0

··· ···

4 1

1 4

−

2uik−j

+

k−j ui+1 )

1

0

C C C C C, C C A

B B B B B B B 

B=

11

1 h

−2 1 0 .. . 0 0

+ τ (f k , φi (x)).

1 0 −2 1 1 −2 .. .. . . 0 0 0 0

··· ··· ··· .. . ··· ···

0 0 0 .. .

0 0 0 .. .

−2 1 1 −2

1 C C C C C, C C A

(15)

Journal Pre-proof

2

3

(f k , φ1 (x)) (f k , φ2 (x)) .. .

6 6

Fk = 6 4

2

7 7 7, 5

3

uk1 uk2 .. .

6 6

Wk = 6 4

7 7 7. 5

ukm−1

(f k , φm−1 (x))

q X

(A −

s=1

k−1 X

(γ ) (bj s

−

(γs ) )BW k−j bj−1

=

!

q X

W k = AW k−1 +

Ks B

k−1 X j=1

it follows that

k−1 X

s=1

s=1

(γs )

(bj

j=1

(γ )

s )BW k−j + τ F k . − bj−1

Ks

(γ )

(γ )

s s (bk−j − bk−j−1 )BW j ,

k−1 X j=1

(γ )

(γ )

s s (bk−j − bk−j−1 )BW j + τ F k .

Pr e-

A−

Ks

s=1

j=1

q X

q X

p ro

Since

Ks B)W k = AW k−1 +

of

We have

(16)

Remark 1. When quadratic basis functions are chosen, the implementation process of the FEM is similar. For example, we take 8 <

φi (x) =

:

(2 xih−x − 1)( xih−x − 1), x ∈ [xi−1 , xi ], x−xi i (2 x−x h − 1)( h − 1), x ∈ [xi , xi+1 ], 0, elsewhere,

where i = 1, 2, · · · , m − 1. Then we have

0

1 15 4 15 1 15 −1 30

0 0 .. .

0

0

0

1 15 8 15 1 15 1 15 −1 30

−1 30 1 15 4 15 8 15 1 15

.. . 0

.. . 0

0 0 0 1 15 1 15 4 15

.. . 0

··· ··· ··· ··· ··· ··· .. .

0 0 0

0 0 0 0 0 0 .. .

0 0 0 0 0 0 .. .

−1 30

0

1 15

.. . 0

···

1 15

8 15

4(x−xi ) (1 h

− 0,

φi+ 12 (x) =

1

0

C C C C C C C, C C C C A

al

A=

8 15 B 1 B 15 B 0 B B 0 B hB 0 B B 0 B B .  ..

urn

0



B=

1 h

16 3 B −8 B 3 B 0 B B 0 B B 0 B B 0 B B .  ..

−8 3 14 3 −8 3 1 3

0 0 .. .

0

0

0 −8 3 16 3 −8 3 −8 3 1 3

.. . 0

x−xi h ),

0 0 −8 3 14 3 16 3 −8 3

0 0 0

0 0 0

−8 3 −8 3 14 3

1 3

.. . 0

.. . 0

x ∈ [xi , xi+1 ], elsewhere,

0 −8 3

.. . 0

··· ··· ··· ··· ··· ··· .. .

0 0 0 0 0 0 .. .

0 0 0 0 0 0 .. .

···

−8 3

16 3

1 C C C C C C C. C C C C A

Theorem 3. Let u(x, tk ) be the exact solution and be sufficiently smooth. Denote U k the solution of fully discrete element method.

Jo

||u(·, tk ) − U k || ≤ C(τ + σ 2 + hl+1 ), k = 1, 2, · · · , n,

where C is a positive constant independent of τ and h. Proof. Let εk = U k − Πh u(x, tk ). From (14)-(15), we have (εk , χ) = (εk−1 , χ) −

q X

"

Ks

s=1

∂εk ∂χ , ∂x ∂x



+

k−1 X

(γ ) (bj s

j=1



−

(γs ) ) bj−1

∂εk−j ∂χ , ∂x ∂x

#

+ (τ Rk , χ).

Taking χ = εk , we have ||εk ||2 = (εk−1 , εk ) −

q X s=1

"

Ks

k−1

∂εk 2 X (γs ) (γs ) || + ) || (bj − bj−1 ∂x j=1 12



∂εk−j ∂εk , ∂x ∂x

#

+ (τ Rk , εk ).

Journal Pre-proof

By using the Cauchy-Schwarz inequality, then we have q

q



k−1

X X ∂εk 2 1 ∂εk−j 2 ∂εk 2 1 X (γ ) (γs ) || + || + || || − bj s ) || ||ε || ≤ (||εk−1 ||2 + ||εk ||2 ) − (bj−1 Ks || Ks 2 ∂x 2 s=1 ∂x ∂x s=1 j=1



k 2

+ (τ Rk , |εk |).

Zk = ||εk ||2 +

q X

Ks

s=1

k−1 X

(γs )

bj

j=1

||

of

Let ∂εk−j 2 || , ∂x

Zk ≤ Zk−1 −

q X s=1

p ro

then applying Lemma 4, we have k

2

k

Note that

s=1

(γ )

s ds bk−1 τ 1+γs + τ σ 2 + τ 2 + τ hl+1 ) ≤ C

q X

(γ )

s ds bk−1 τ γs (τ + σ 2 + hl+1 ).

s=1

Pr e-

q X

τ ||Rk || ≤ C(

(γ )

q X

2b s Ks k−1 ||εk ||2 + (τ Rk , |εk |). || + (τ Rk , |ε |) ≤ Zk−1 − ∂x L2 s=1

(γs ) ∂ε Ks bk−1 ||

Therefore, by Cauchy-Schwarz inequality, we have (γ )

q X

Ks

s=1

≤ Zk−1 + C ≤ Zk−1 + C

q X s=1 q X

(γ )

≤

s=1

Hence we have

(γ )

s=1

ds Γ(γs +

(γs ) γs τ (τ 1)bj−1

2

+σ +h

l+1 2

j=0 s=1

q X

(γ )

s ds Γ(γs + 1)bk−1 τ γs (τ + σ 2 + hl+1 )2 .

urn

Zk ≤ C

q

s d2s L2 bk−1 τ 2γs (τ + σ 2 + hl+1 )2 2Ks

Since Z0 = 0, it follows that q k−1 X X

(γ )

q

s s X X d2 L2 b s 2bk−1 2bk−1 s k−1 2γs k 2 k 2 K ||ε || + ||ε || + C τ (τ + σ 2 + hl+1 )2 s 2 L2 L 2K s s=1 s=1

al

Zk ≤ Zk−1 −

) ≤

q X

ds Γ(γs + 1)(kτ )γs (τ + σ 2 + hl+1 )2

s=1

ds Γ(γs + 1)T γs (τ + σ 2 + hl+1 )2 ≤ C(τ + σ 2 + hl+1 )2 .

Jo

||U k − Πh u(x, tk )|| ≤ C(τ + σ 2 + hl+1 ).

Therefore, we have

||u(x, tk ) − U k || ≤ ||u(x, tk ) − Πh u(x, tk )|| + ||U k − Πh u(x, tk )|| ≤ C(τ + σ 2 + hl+1 ), k = 0, 1, 2, · · · , n. In the following, we want to increase the order of convergence in time and obtain a new finite element method. Now we first state the following lemma, which will be used later. Lemma 6. [41] If the function v(x, t) is sufficiently smooth, then v(xi , η) =

(tj+1 − η)v(xi , tj ) + (η − tj )v(xi , tj+1 ) + O(τ 2 ). τ 13

Journal Pre-proof

Next, integrating both sides of (6) with respect to t, we have u(xi , tk+1 ) − u(xi , 0) =

q X

Z t k+1

ds 0

s=1

K ∂ 2 u(xi , η) dη + (tk+1 − η)γs −1 Γ(γs ) ∂x2

Z t k+1

(f (xi , η) + R1 )dη.

(17)

0

u(xi , tk+1 ) − u(xi , 0) ds

s=1

k Z tj+1 X tj

j=0

∂ 2 u(x ,t )

K (tj+1 − η) ∂x2i j + (η − tj ) Γ(γs ) τ (tk+1 − η)1−γs

Thus q X

(1)

(1) ∂

2

[Cj

s=1

j=0

where Cj

k X

Ks

= (j + 1)γs −

Z t k+1

dη +

(f (xi , η) + R1 )dη + O(τ 2 ).

0

2 u(xi , tk−j ) (2) ∂ u(xi , tk−j+1 ) + C ]+ j ∂x2 ∂x2

Z t k+1

(f (xi , η) + R1 )dη + O(τ 2 ),

0

1 1 (2) [(j + 1)γs +1 − j γs +1 ], Cj = [(j + 1)γs +1 − j γs +1 ] − j γs . γs + 1 γs + 1

Pr e-

u(xi , tk+1 ) − u(xi , 0) =

∂ 2 u(xi ,tj+1 ) ∂x2

p ro

=

q X

of

Applying Lemma 6, we have

The weak form of (6)-(8) subject to the boundary condition

uk+1 (0) = 0, uk+1 (L) = 0, reads: find uk+1 ∈ H01 (Ω) such that (u(xi , tk+1 ), ϕ) − (u(xi , 0), ϕ) =−

q X s=1

Ks

k • X



(1) Cj

j=0

∂u(xi , tk−j ) ∂ϕ , ∂x ∂x

‹



+

(2) Cj

∂u(xi , tk−j+1 ) ∂ϕ , ∂x ∂x

Z t k+1

‹˜

+



(f (xi , η) + R1 )dη, ϕ 0

0

, χ) − (U , χ) = −

for ∀χ ∈ Sh .

q X

Ks

k  X



(1) Cj

urn

(U

k+1

al

for ∀ ϕ ∈ H01 (Ω). Therefore, let U k+1 be the approximation of u(x, tk+1 ), we define the Galerkin finite element method for (6)-(8) by

s=1

j=0

∂U k−j ∂χ , ∂x ∂x





+

(2) Cj

∂U k−j+1 ∂χ , ∂x ∂x



Z t k+1

+



f (x, η)dη, χ

(18)

0

4. Numerical Results

Jo

In this section, we give two numerical examples to verify the efficiency of the proposed numerical method for the fractional diffusion equation. Example 1 Consider the following distributed-order anomalous sub-diffusion equation: 8 < :

∂u(x,t) ∂t

R1

2

= 0 ω(γ)0 Dt1−γ (K ∂∂xu2 )dγ + f (x, t), u(x, 0) = 0, 0 ≤ x ≤ 1, u(0, t) = u(1, t) = 0, 0 ≤ t ≤ 21 , 2

0 < x < 1, 0 < t ≤ 12 ,

2 where K = 1, ω(γ) = Γ(2 + γ) and f (x, t) = 2tx2 (1 − x)2 − 4 tln−t t (6x − 6x + 1). The exact solution is given by 2 2 2 u(x, t) = x (1 − x) t . We first choose the piecewise linear polynomials (l = 1) as the basis function in the finite element approximation. Table 1 displays the L2 error and convergence order of σ for τ = 1/1000 and h = σ. Table 2 shows the L2 error and convergence order of h for τ = 1/1000 and σ = 1/40. In Table 3, we obtain the L2 error and convergence order of τ for h = 1/100 and σ = 1/40. Clearly, from Tables 1-3, it can be seen that the numerical solutions fit well with

14

Journal Pre-proof

of

the exact solutions, the second-order convergence in σ and space, and first-order convergence in time are observed, which is in good agreement with the theoretical analysis. Next, the quadratic polynomials (l = 2) are adopted in the finite element approximation and the numerical results are displayed in Table 4. From Table 4, it can be observed that the accuracy in space is O(h3 ), which is in agreement with the theoretical analysis. Furthermore, we applied the high order time scheme (18) to the equation and the related numerical results are illustrated in Table 5. We can see that the convergence order of time has been improved from first order to second order. Finally, we present a comparison of the numerical solution and the exact solution in Fig 1 for σ = 1/40, h = 1/50, τ = 1/100 at t = 1/2. Table 1: The L2 error and convergence order of σ for l = 1, τ = 1/1000 and h = σ. L2 error 2.9579E-04 7.5349E-05 1.8379E-05 4.4427E-06

Order – 1.97 2.04 2.05

p ro

σ 1/8 1/16 1/32 1/64

Pr e-

Table 2: The L2 error and convergence order of h for l = 1, τ = 1/1000 and σ = 1/40. h 1/8 1/16 1/32 1/64

L2 error 2.9820E-04 7.5873E-05 1.8424E-05 4.4551E-06

Order – 1.97 2.04 2.05

Table 3: The L2 error and convergence order of τ for l = 1, σ = 1/40, and h = 1/100. L2 error 3.6261E-04 1.6740E-04 7.7406E-05 3.5872E-05

Example 2 source term:

urn

al

τ 1/10 1/20 1/40 1/80

Order – 1.12 1.11 1.11

Next, we consider the following distributed-order anomalous sub-diffusion equation without a 8 <

R1

2

= 0 ω(γ)0 Dt1−γ (K ∂∂xu2 )dγ, u(x, 0) = δ(x − 21 ), 0 ≤ x ≤ 1, u(0, t) = u(1, t) = 0, 0 ≤ t ≤ 1,

Jo

:

∂u(x,t) ∂t

0 < x < 1, 0 < t ≤ 1,

(19)

1 1 1 where K = 1 and δ(x) is the Kronecker delta function. In the simulation, we choose h = 100 , σ = 100 , and τ = 100 . Here, we perform the numerical tests for three different ω(γ) cases. Firstly, we consider the case ω(γ) = δ(γ − γ0 ). Then the equation (19) reduces to the single term anomalous sub-diffusion equation (1). Now we apply the finite element method to the equation to observe the diffusion behaviour of the solution. Figure 2 shows the profiles of the solution for different γ0 at t = 1. We can see that the solution decays with increasing γ0 . Particularly, when γ0 = 1, equation (1) is a normal diffusion equation, which means the solution becomes less anomalous with increasing γ0 . Figure 3 displays the profiles of the solution at different time with γ0 = 0.705. We can observe that the solution decays with increasing t and becomes less anomalous over the course of time. Next, we consider the case ω(γ) = k1 δ(γ − γ1 ) + (1 − k1 )δ(γ − γ2 ), 0 < k1 < 1. For this case, the equation (19) can be simplified as the two term anomalous sub-diffusion equation (2). In the following, we study the diffusion behaviour of the solution in this case. Figure 4 shows the profiles of the solution for different pairs of (γ1 , γ2 ) with k1 = 0.3 at t = 1. Figure

15

Journal Pre-proof

Table 4: The L2 error and convergence order of h for l = 2, σ = 1/40 and τ = h3 . L2 error 1.0021E-03 1.5018E-04 1.9582E-05 2.5233E-06

Order – 2.74 2.94 2.96

of

h 1/2 1/4 1/8 1/16

Table 5: The L2 error and convergence order of τ for the high-order time scheme with l = 1, σ = 1/40, and h = τ . L2 error 7.9210E-05 1.9765E-05 4.7818E-06 1.0898E-06

Order – 2.00 2.05 2.13

p ro

τ 1/16 1/32 1/64 1/128

5. Conclusions

urn

al

Pr e-

5 displays the profiles of the solution at different time with (γ1 , γ2 ) = (0.505, 0.995). Similarly, we can observe that the solution decays with increasing t and becomes less anomalous over the course of time. We note that this case shows a more complex behavior, a non-Gausssian behavior at small times and approximately to Gausssian one at large times. This is the case of accelerating sub-diffusion, which means that the process shows anomalous sub-diffusion with a diffusion exponent 0.505 at small times and sub-diffusion with a diffusion exponent 0.995 at large times. As for the physical meaning of the modified distributed-order anomalous sub-diffusion equation, it can be viewed as a composition of two independent continuous time random walk (CTRW) processes. Analogously, if the process governed by the modified distributed-order anomalous sub-diffusion equation with N delta functions, it can be viewed as a sum of N independent CTRW processes, which correspond to the pooling of CTRWs [43]. Now we consider some special cases of weight functions ω(γ) include delta function, linear function, quadratic function and exponential function. In Figure 6, we present a comparison of the numerical solution profiles for different ω(γ) 1−γ 2 γ−1 . cases: ω1 (γ) = δ(γ − 0.705), ω2 (γ) = 0.3δ(γ − 0.305) + 0.7δ(γ − 0.705), ω3 (γ) = 1−γ 2 , ω4 (γ) = 2 , ω5 (γ) = 2 We can conclude that the numerical solution profiles are different for choosing distinct ω(γ). It can be seen that the numerical solution decays faster with the exponential weight function ω5 (γ) = 2γ−1 than the other weight functions. The figure may give us some hints to select right weight functions in the modified distributed-order anomalous subdiffusion equation to model many physical processes. Finally, we study the numerical solution profiles for the case ω(γ)= 1−γ 2 at different time. Again, the solution decays with increasing t.

Jo

In this paper, we have dealt with the numerical solutions to the distributed-order anomalous sub-diffusion equation. By using the mid-point quadrature rule, we discretize the integral term in the distributed-order anomalous sub-diffusion equation and transform it into a multi-term fractional diffusion equation. Then we present the variation formulation and construct a semi-discrete finite difference scheme. The stability and convergence of the semidiscrete finite difference scheme are proved by using new energy method. Furthermore, we derive a Galerkin finite element scheme for space and obtain the fully discrete scheme. The convergence orders in discrete L2 norm are O(τ 2 + σ 2 + hl+1 ). Numerical results show the efficiency and numerical accuracy of the proposed scheme. Our results allow us to suggest that the modified distributed-order anomalous sub-diffusion equation could be a useful tool to describe the situations demonstrating anomalous sub-diffusion at small times and normal diffusion at large times. Acknowledgments Author Li wishes to acknowledge that this research was supported by the National Natural Science Foundation of China (Grant No.11401223), the Science and Technology Program of Guangzhou (No.201707010031) and the

16

Journal Pre-proof

0.016 Exact Numerical

0.014 0.012 0.01

of

0.008 0.006 0.004

0 0

0.1

0.2

0.3

0.4

0.5

x

p ro

0.002

0.6

0.7

0.8

0.9

1

Figure 1: The comparison of the numerical solution and exact solution for σ =

0.05

1 40 ,

h=

1 50 ,

τ=

1 100

at t = 21 .

Pr e-

γ 0=0.105

0.045

0.04

0.035

u(x, T )

0.03

0.025

0.02

0.015

0.01

γ 0=0.305 γ 0=0.505 γ 0=0.705 γ 0=0.905 γ 0=0.995

0.005

al

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x

urn

Figure 2: The numerical solution profiles of u(x, t) for different γ0 with ω(γ) = δ(γ − γ0 ) at T = 1.

0.09

t=0.1 t=0.2 t=0.5 t=0.8 t=1.0

0.08

0.07

u(x, t)

Jo

0.06

0.05

0.04

0.03

0.02

0.01

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x

Figure 3: The numerical solution profiles of u(x, t) at different time with γ0 = 0.705 and ω(γ) = δ(γ − γ0 ).

17

Journal Pre-proof

0.045 γ 1 =0.105,γ 2 =0.305 γ 1 =0.305,γ 2 =0.505

0.04

γ 1 =0.505,γ 2 =0.705

of

γ 1 =0.705,γ 2 =0.905 γ 1 =0.905,γ 2 =0.995

0.035

0.025

0.02

0.015

0.01

0.005

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Pr e-

x

p ro

u(x, T )

0.03

al

Figure 4: The numerical solution profiles of u(x, t) for different (γ1 ,γ2 ) with ω(γ) = k1 δ(γ − γ1 ) + (1 − k1 )δ(γ − γ2 ), k1 = 0.3 at T = 1.

0.09

t=0.1 t=0.2 t=0.5 t=0.8 t=1.0

urn

0.08 0.07 0.06 0.05 0.04 0.03

Jo

0.02 0.01 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figure 5: The numerical solution profiles of u(x, t) at different time with (γ1 , γ2 ) = (0.505, 0.995) and ω(γ) = k1 δ(γ − γ1 ) + (1 − k1 )δ(γ − γ2 ), k1 = 0.3.

18

Journal Pre-proof

0.04 ω(γ)=ω 1 (γ) ω(γ)=ω 2 (γ) ω(γ)=ω 3 (γ)

of

0.035

ω(γ)=ω 4 (γ) ω(γ)=ω 5 (γ)

0.03

0.02

0.015

0.01

0.005

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Pr e-

x

p ro

u(x, T )

0.025

al

Figure 6: The numerical solution profiles of u(x, t) for different ω(γ) with at T = 1.

0.012

urn

t=0.1 t=0.2 t=0.5 t=0.8 t=1.0

0.01

u(x, t)

0.008

0.006

Jo

0.004

0.002

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x

Figure 7: The numerical solution profiles of u(x, t) at different time for the case ω(γ)= 1−γ 2 .

19

Journal Pre-proof

China Scholarship Council. Author Liu wishes to acknowledge that this research was partially supported by the Australian Research Council (ARC) via the Discovery Projects (DP180103858 and DP190101889) and National Natural Science Foundation of China (Grant No.11772046). Author Feng wishes to acknowledge that this research was partially supported by the National Natural Science Foundation of China (Grant No. 11801060). Author Turner wishes to acknowledge that this research was supported by the ARC Discovery Project (DP180103858).

of

References

[1] A. Kilbas, H.M. Srivastava, J. Trujillo, Theory and Applications of Fractional Differential Equations, Elsevier Science, 2006.

p ro

[2] I. Podlubny, Fractional Differential Equations, first ed., in: Mathematics in Science and Engineering, vol. 198, Academic Press, New York, 1998. [3] K. Diethelm, The Analysis of Fractional Differential Equations, Springer, Berlin, 2010. [4] F. Liu, P. Zhuang, Q. Liu, Numerical Methods of Fractional Partial Differential Equations and Applications, Science Press, Beijing, 2015.

Pr e-

[5] H. Zhang, X. Jiang, C. Wang, W. Fan, Galerkin-Legendre spectral schemes for nonlinear space fractional Schr¨odinger equation, Numer. Algorithms 78 (2018) 337-356. [6] H. Zhang, X. Jiang, X. Yang, A time-space spectral method for the time-space fractional Fokker-Planck equation and its inverse problem, Appl. Math. Comput. 320 (2018) 302-318. [7] R. Ghosh, W. Webb, Automated detection and tracking of individual and clustered cell surface low density lipoprotein receptor molecules, Biophys. J. 66 (1994) 1301-1318. [8] T. Feder, I. Brust-Mascher, J. Slattery, B. Baird, W. Webb, Constrained diffusion or immobile fraction on cell surfaces: a new interpretation, Biophys. J. 70 (1996) 2767-2773.

al

[9] E. Sheets, G. Lee, R. Simson, K. Jacobson, Transient confinement of a glycosylphosphatidylinositol-anchored protein in the plasma membrane, Biochemistry 36 (1997) 12449-12458. [10] E. Brown, E. Wu, W. Zipfel, W. Webb, Measurement of molecular diffusion in solution by multiphoton fluorescence photobleaching recovery, Biophys. J. 77 (1999) 2837-2849

urn

[11] S.B. Yuste, L. Acedo, An explicit finite difference method and a new vonneumann-type stability analysis for fractional diffusion equations, SIAM J. Numer. Anal. 42 (2005) 1862-1874. [12] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: A fractional dynamic approach, Phys. Rep. 339 (2000) 1-77. [13] P. Zhuang, F. Liu, I. Turner, V. Anh, Galerkin finite element method and error analysis for the fractional cable equation, Numer. Algorithms 72 (2016) 447-466.

Jo

[14] S.B. Yuste, K. Lindenberg, Subdiffusion limited A+A reactions, Phys. Rev. Lett. 87 (2001) 118301. [15] T.A.M. Langlands, Solution of a modified fractional diffusion equation, Physica A, 367 (2006) 136-144. [16] F. Liu, C. Yang, K. Burrage, Numerical method and analytical technique of the modified anomalous subdiffusion equation with a nonlinear source term, J. Comput. Appl. Math. 231 (2009) 160-176. [17] A.V. Chechkin, R. Gorenflo, I.M. Sokolov, Retarding subdiffusion and accelerating superdiffusion governed by distributed-order fractional diffusion equations, Phys. Rev. E 66 (2002) 046129. [18] I. Sokolov, A. Chechkin, J. Klafter, Distributed-order fractional kinetics, Acta. Phys. Pol. B 35 (2004) 1323-1341. [19] Ya. G. Sinai, Limit behaviour of one-dimensional random walks in random environments, Theor. Prob. Appl. 27 (1982) 247-258. 20

Journal Pre-proof

[20] R.N. Mantegna, H.E. Stanley, Stochastic process with ultraslow convergence to a Gaussian: The truncated L´evy flight, Phys. Rev. Lett. 73 (1994) 2946. [21] A.V. Chechkin, V.Yu. Gonchar, R. Goreflo, N. Korabel, I.M. Sokolov, Generalized fractional diffusion equations for accelerating subdiffusion and truncated L´evy flights, Phys. Rev. E 78 (2008) 021111.

of

[22] M. Caputo, Elasticit´ ae Dissipazione, Zanichelli, Bologna, 1969. [23] C.H. Eab, S.C. Lim, Fractional Langevin equations of distributed order, Phys. Rev. E 83 (2011) 031136.

p ro

[24] M. Naghibolhosseini, Estimation of outer-middle ear transmission using DPOAEs and fractional-order modeling of human middle ear, Ph.D thesis, City University of New York, NewYork, 2015. [25] T. Atanackovic, M. Budincevic, S. Pilipovic, On a fractional distributed-order oscillator. J. Phys. A Math. Gen. 38 (2005) 6703. [26] Ninghu Su, Paul N. Nelson, Sarah Connor, The distributed-order fractional diffusion-wave equation of groundwater flow: Theory and application to pumping and slug tests, J.Hydrol. 529 (2015) 1262-1273. [27] Z. Jiao, Y. Chen, I. Podlubny, Distributed-order dynamic systems, Springer, Berlin, 2012.

Pr e-

[28] M. Caputo, Distributed order differential equations modelling dielectric induction and diffusion, Fract. Calc. Appl. Anal. 4 (2001) 421-442. [29] F. Mainardi, G. Pagnini, A. Mura, R. Gorenflo, Time-fractional diffusion of distributed order, J. Vib.Control 14 (2008) 1267-1290. [30] R. Gorenflo, Y. Luchko, M. Stojanove´c, Fundamental solution of a distributed order time-fractional diffusionwave equation as probability density, Fract. Calc. Appl. Anal. 16 (2013) 297-316. [31] Y. Luchko, Boundary value problems for the generalized time-fractional diffusion equation of distributed order, Fract. Calc. Appl. Anal. 12 (2009) 409-422.

al

[32] H. Ye, F. Liu, V. Anh, Compact difference scheme for distributed-order time-fractional diffusion-wave equation on bounded domains, J. Comput. Phys. 298 (2015) 652-660. [33] Z. Sun, X. Wu, A fully discrete difference scheme for a diffusion-wave system, Appl. Numer. Math. 56 (2006) 193-209.

urn

[34] G. Gao, Z. Sun, Two alternating direction implicit difference schemes with the extrapolation method for the two-dimensional distributed-order differential equations, Comput. Math. Appl. 69 (2015) 926-948. [35] G. Gao, H. Sun, Z. Sun, Some high-order difference schemes for the distributed-order differential equations, J. Comput. Phys. 298 (2015) 337-359. [36] J. Li, F. Liu, L. Feng, I. Turner, A novel finite volume method for the Riesz space distributed-order diffusion equation, Comput. Math. Appl. 74 (2017) 772-783.

Jo

[37] J. Li, F. Liu, L. Feng, I. Turner, A novel finite volume method for the Riesz space distributed-order advectiondiffusion equation, Appl. Math. Model. 46 (2017) 536-553. [38] H. Zhang, F. Liu, X. Jiang, F. Zeng, I. Turner, A Crank-Nicolson ADI Galerkin-Legendre spectral method for the two-dimensional Riesz space distributed-order advection-diffusion equation, Comput. Math. Appl. 76 (2018) 2460-2476. [39] W. Fan, F. Liu, A numerical method for solving the two-dimensional distributed order space-fractional diffusion equation on an irregular convex domain, Appl. Math. Lett. 77 (2018) 114-121. [40] P. Zhuang, F. Liu, V. Anh, I. Turner, New solution and analytical techniques of the implicit numerical method for the anomalous subdiffusion equation, SIAM J. Num. Anal. 46 (2008) 1079-1095.

21

Journal Pre-proof

[41] L. Feng, P. Zhuang, F. Liu, I.Turner, Y. Gu, Finite element method for space-time fractional diffusion equation, Numer. Algorithms 72 (2016) 749-767. [42] V. Thom´ee, Galerkin Finite Element Mehtods for Parabolic Problems, Springer, Berlin, 2006.

Jo

urn

al

Pr e-

p ro

of

[43] T. Sandev, A.V. Chechkin, N. Korabel, H. Kantz, I.M. Sokolov, R. Metzler, Distributed-order diffusion equations and multifractality: Models and solutions, Phys. Rev. E 92 (2015) 042117.

22