finite element method for the 2D time-space Riesz fractional diffusion equation on irregular convex domains

Accepted Manuscript

Unstructured mesh finite difference/finite element method for the 2D time-space Riesz fractional diffusion equation on irregular convex domains Libo Feng, Fawang Liu, Ian Turner, Qianqian Yang, Pinghui Zhuang PII: DOI: Reference:

S0307-904X(18)30056-8 10.1016/j.apm.2018.01.044 APM 12160

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

16 February 2017 16 January 2018 25 January 2018

Please cite this article as: Libo Feng, Fawang Liu, Ian Turner, Qianqian Yang, Pinghui Zhuang, Unstructured mesh finite difference/finite element method for the 2D time-space Riesz fractional diffusion equation on irregular convex domains, Applied Mathematical Modelling (2018), doi: 10.1016/j.apm.2018.01.044

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

ACCEPTED MANUSCRIPT

Highlights • Establish and prove some new definitions and lemmas of fractional derivative space on convex domains; • An unstructured mesh finite element method is presented; • The stability and convergence of the method are discussed on different irregular convex domains; • Extend the theory and develop a computational model for the case of a multiply-connected domain;

AC

CE

PT

ED

M

AN US

CR IP T

• Apply to solve the coupled 2D fractional Bloch-Torrey equation on human brain-like domains.

1

ACCEPTED MANUSCRIPT

Unstructured mesh finite difference/finite element method for the 2D time-space Riesz fractional diffusion equation on irregular convex domains Libo Fenga , Fawang Liua,∗, Ian Turnera,b , Qianqian Yanga , Pinghui Zhuangc,d a School

of Mathematical Sciences, Queensland University of Technology, GPO Box 2434, Brisbane, QLD. 4001, Australia Research Council Centre of Excellence for Mathematical and Statistical Frontiers (ACEMS), Queensland University of Technology (QUT), Brisbane, Australia. c School of Mathematical Sciences, Xiamen University, Xiamen 361005, China d Fujian Provincial Key Laboratory of Mathematical Modeling and High-Performance Scientific Computation, Xiamen University, Xiamen 361005, China

CR IP T

b Australian

Abstract

M

AN US

Fractional differential equations are powerful tools to model the non-locality and spatial heterogeneity evident in many real-world problems. Although numerous numerical methods have been proposed, most of them are limited to regular domains and uniform meshes. For irregular convex domains, the treatment of the space fractional derivative becomes more challenging and the general methods are no longer feasible. In this work, we propose a novel numerical technique based on the Galerkin finite element method (FEM) with an unstructured mesh to deal with the space fractional derivative on arbitrarily shaped convex and non-convex domains, which is the most original and significant contribution of this paper. Moreover, we present a second order finite difference scheme for the temporal fractional derivative. In addition, the stability and convergence of the method are discussed and numerical examples on different irregular convex domains and non-convex domains illustrate the reliability of the method. We also extend the theory and develop a computational model for the case of a multiply-connected domain. Finally, to demonstrate the versatility and applicability of our method, we solve the coupled two-dimensional fractional Bloch-Torrey equation on a human brain-like domain and exhibit the effects of the time and space fractional indices on the behaviour of the transverse magnetization.

ED

Keywords: finite element method, unstructured mesh, Riesz fractional derivative, irregular domains, two-dimensional, time-space fractional diffusion equation

PT

1. Introduction

AC

CE

Over the past decades, fractional derivatives have been widely used in physics [1, 2], biology [3, 4], chemistry [5], hydrology [6, 7], finance [8] and the related theory has been expanding at a fast rate [9, 10, 11]. A considerable number of computational models emerged that are based on applying the finite element method (FEM) to fractional diffusion equations (FDE). This work dates back to Roop and Ervin [12, 13], who constructed appropriate fractional derivative spaces and presented a theoretical framework for the Galerkin finite element approximation to the fractional advection dispersion equation. Deng [14] developed the FEM for the numerical resolution of the space and time fractional Fokker-Planck equation. Zhang et al. [15] discussed the Galerkin finite element approximation of symmetric space-fractional partial differential equations. Li et al. [16] studied the Galerkin FEM for time-space fractional order nonlinear subdiffusion and superdiffusion equations. Liu et al. [17] investigated the finite element approximation for a modified anomalous subdiffusion equation. Zeng et al. [18] developed finite difference and finite element approaches for the time-fractional subdiffusion equation with Dirichlet boundary conditions. Jin et al. [19] studied the Galerkin FEM and lumped mass Galerkin FEM for the initial boundary value problem of a homogeneous time-fractional diffusion equation. Jin et al. also [20] considered the Galerkin FEM for the initial/boundary value problem involving multiple time-fractional derivatives on a bounded convex polyhedral domain. Bu et al. [21] considered the Galerkin FEM for two-dimensional Riesz space FDE. Liu et al. [22] presented a mixed FEM ∗ Corresponding

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

Preprint submitted to Applied Mathematical Modelling

February 15, 2018

ACCEPTED MANUSCRIPT

C γ 0 Dt u(x, y, t)

= K1

∂ 2α u(x, y, t) ∂|x|

+ K2

∂ 2β u(x, y, t) 2β

∂|y|

M

with the initial condition

2α

AN US

CR IP T

for a time-fractional fourth-order partial differential equation. Feng et al. [23] considered the FEM with a secondorder time scheme for a space-time FDE. Zhuang et al. [24] introduced the Galerkin FEM and error analysis for the fractional cable equation. Recently, Zhao et al. [25] established the nonconforming FEM for two-dimensional multi-term time fractional subdiffusion equations. Jin et al. [26] developed variational formulations of a PetrovGalerkin FEM type for one-dimensional fractional boundary value problems. Yang et al. [27] considered the FEM for nonlinear Riesz space fractional diffusion equations on irregular domains. Fan et al. [28] discussed the FEM for the two-dimensional time-space fractional wave equation on irregular domains. In fact, many problems from science and engineering involve mathematical models that must be computed on irregular domains and therefore seeking effective numerical methods to solve FDE on such domains is important. Although existing numerical methods for FDE are numerous, most of them are limited to regular domains and uniform meshes. Research involving unstructured meshes [28, 29] and irregular domains [27, 30, 31] is more sparse. For the classical diffusion equation with integer order derivatives, there is some theory and research involving unstructured meshes with the finite volume element method (see [32, 33, 34] and references therein). For the fractional case, recently, Karaa et al. [35] proposed a finite volume element method with unstructured mesh for approximating the anomalous subdiffusion equations with temporal fractional derivative. Le et al. [36] studied the finite element approximation for a time-fractional diffusion problem on a domain with a re-entrant corner. Fan et al. [28] discussed the unstructured mesh finite element method for the time-space fractional wave equation. They used the L2 scheme to approximate the temporal derivative with low accuracy and extended the properties of the fractional derivative space to the two-dimensional convex domain case without detailed proof. Here, we will consider a finite element method suitable for implementation with an unstructured mesh for the time-space fractional diffusion equation. In this paper, we consider the following two-dimensional time-space Riesz fractional diffusion equation (2D TSRFDE) on an irregular convex domain:

u(x, y, 0) = φ(x, y), and boundary condition

ED

u(x, y, t) = 0,

1 2

+ f (x, y, t),

(x, y, t) ∈ Ω × (0, T ],

(x, y) ∈ Ω,

(1)

(2)

(x, y, t) ∈ ∂Ω × (0, T ],

(3)

where 0 < γ < 1, < α, β ≤ 1, K1 > 0, K2 > 0, f (x, y, t) and φ(x, y) are two known smooth functions. The irregular convex domain Ω is defined as (see Figure 1): Ω = {(x, y)| c(y) < x < r(y), a1 < y < b1 } or Ω = {(x, y)| g(x) < y < m(x), c1 < x < d1 }, where a1 = min g(x), b1 = max m(x), c1 = min c(y), and

PT

(x,y)∈Ω

(x,y)∈Ω

(x,y)∈Ω

d1 = max r(y). Here the irregular convex domains means the convex domain in with a curved boundary, which is (x,y)∈Ω

y

y

b1 m(x)

AC

CE

compared with the regular convex domain such as a square or rectangle. In Eq.(1), the Caputo fractional derivative

c(y)

r(y)

g(x) a1

O

x

O

c1

d1

Figure 1: The boundaries of convex domain Ω .

3

x

ACCEPTED MANUSCRIPT

C γ 0 Dt u

is defined as [10] C γ 0 Dt u(x, y, t)

=

Z

1 Γ(1 − γ)

t

0

(t − s)−γ

∂u(x, y, s) ds. ∂s

Here, we define the Riesz fractional derivatives

∂ 2α u ∂|x|2α

and

∂ 2β u ∂|y|2β

as

CR IP T

As Ω is an irregular convex domain, the boundary value in the space fractional derivative is no longer a fixed constant, which is distinct from the common Riemann-Liouville fractional derivative definition [9] Z x ∂n 1 2α (x − s)n−2α−1 u(s, y, t) ds, a Dx u(x, y, t) = Γ(n − 2α) ∂xn a Z b (−1)n ∂ n 2α (s − x)n−2α−1 u(s, y, t) ds. x Db u(x, y, t) = Γ(n − 2α) ∂xn x   1 ˜ 2α ˜ 2α c(y) Dx u(x, y, t) + x Dr(y) u(x, y, t) , 2 cos(απ) ∂|x|   2β ∂ u(x, y, t) 1 ˜ y2β u(x, y, t) + y D ˜ 2β u(x, y, t) , = − D g(x) m(x) 2β 2 cos(βπ) ∂|y| 2α

=−

AN US

∂ 2α u(x, y, t)

˜ y2β u(x, y, t) ˜ 2α u(x, y, t), g(x) D ˜ x2α u(x, y, t), x D and the Riemann-Liouville fractional derivative with varying boundary c(y) D r(y) ˜ 2β u(x, y, t) (n − 1 < 2α, 2β < n) are given by and y D m(x)

˜ 2α c(y) Dx u(x, y, t) :=

1 ∂n Γ(n − 2α) ∂xn

Z

x

c(y)

(x − s)n−2α−1 u(s, y, t) ds,

Z r(y) (−1)n ∂ n (s − x)n−2α−1 u(s, y, t) ds, Γ(n − 2α) ∂xn x Z y 1 ∂n 2β ˜ D u(x, y, t) := (y − s)n−2β−1 u(x, s, t) ds, g(x) y Γ(n − 2β) ∂y n g(x) Z m(x) (−1)n ∂ n 2β ˜ (s − y)n−2β−1 u(x, s, t) ds. y Dm(x) u(x, y, t) := Γ(n − 2β) ∂y n y

ED

M

˜ 2α x Dr(y) u(x, y, t) :=

AC

CE

PT

The major contribution of this paper is as follows. Firstly, we establish and prove some new definitions and lemmas for convex domains, which extends the properties of the fractional derivative space from the one-dimensional case to the two-dimensional convex domain case. We believe this is a new contribution to the literature that was not discussed in [28, 37]. Secondly, we propose a novel technique utilizing FEM and unstructured triangular meshes to deal with the space fractional derivative on an irregular convex domain, which we believe is very flexible because our considered solution domain can be arbitrarily convex and we need less grid nodes to generate the meshes. For general convex domains, such as a human brain-like domain, the software Gmsh [38] can be used to partition and generate regular unstructured triangular meshes (see figure 2). Due to the irregular domain, the treatment of the space fractional derivative is not straightforward. Therefore, we reduce the calculation from the whole domain Ω to every single triangular element and deal with it approximately by the Gauss quadrature technique, which is elaborated in section 3. In [28], Fan et al. considered the time-space fractional wave equation, of which the order of time fractional derivative is 1 < γ < 2 and discretised the time fractional derivative using the L2 scheme, which requires discretisation of two time levels. Different to [28], we consider application of the FEM to time-space fractional diffusion equations, of which the order of time fractional derivative is 0 < γ < 1. The time fractional derivative is discretised using the L1 scheme, which only requires discretisation of the time fractional derivative at one temporal level and is more computational efficient. Furthermore, we present a second order numerical scheme for the temporal fractional derivative based on finite difference. Moreover, we extend the method to non-convex domains and solve a problem on a multiply-connected domain, respectively. Another important contribution of our work is the extension of the theory to allow the solution of the 2D TSRFDE on a multiply-connected domain and we illustrate the numerical results for a particular case. In order to demonstrate the applicability of our computational approach, we solve the two-dimensional coupled fractional Bloch-Torrey equation on a human brain-like domain numerically and analyse the impact of the space 4

ACCEPTED MANUSCRIPT

fractional index on the diffusion behavior. In clinical settings, diffusion-weighted imaging (DWI) is increasingly utilised to study heterogeneous water diffusion in the human brain, which is thought to be anomalous. In [39], the researchers used the stretched exponential model to investigate the diffusion behavior in complex biological tissues, such as human brain gray matter, glioblastoma tissues and erythrocyte ghosts. In [40], Magin et al. showed that the model was a fundamental extension of the classical Bloch-Torrey equation through application of the operators of fractional calculus. The time and space fractional Bloch-Torrey equation has the following form [40, 41]: ã Å 2β ∂ 2β Mxy (r, t) ∂ 2β Mxy (r, t) α 2(β−1) ∂ Mxy (r, t) + λ(t)Mxy (r, t), (4) ω α−1C + + D M (r, t) = Dµ xy 0 t 2β 2β 2β ∂|x| ∂|y| ∂|z|

AN US

CR IP T

where λ(t) = −iγ(r · G(t)) , G(t) is the magnetic field gradient, γ and D are the gyromagnetic ratio and the √ diffusion coefficient respectively, ω and µ are time and space constants, r = (x, y, z), i = −1 and Mxy (r, t) = Mx (r, t) + iMy (r, t). Therefore, models using fractional-order calculus have emerged as promising tools to analyse diffusion images of the human brain to provide new insights into the investigations of tissue structures and the microenvironment. For the two dimensional case, we equate real and imaginary components to express equation (4) as a coupled system for the components Mx and My , i.e.,  ã Å 2β  ∂ 2β Mx (x,y,t) α 2(β−1) ∂ Mx (x,y,t)  + + λ0 (t)My (x, y, t), D M (x, y, t) = Dµ  ω α−1C x 2β 2β t 0 ∂|x| ∂|y| Å ã (5) 2β  ∂ 2β My (x,y,t) 2(β−1) ∂ My (x,y,t) α  + − λ (t)M (x, y, t).  ω α−1C 0 x 2β 2β 0 Dt My (x, y, t) = Dµ ∂|x| ∂|y|

ED

M

In [41], Yu et al. proposed an implicit numerical method for the two dimensional time-space fractional Bloch-Torrey

Figure 2: The sectional view of human brain [42] and unstructured mesh partition.

AC

CE

PT

equation on a finite rectangular domain. However, this method may have low spatial accuracy (O(τ α + h2 + ρh), 0 ≤ ρ ≤ 1) for problems where the solution domain is irregular [43], for example a human brain-like domain (see Figure 2(a)). The outline of this paper is as follows. In section 2, some new definitions and properties of the fractional derivative space and fractional Sobolev space on irregular convex domains are introduced. In section 3, we derive the fully discrete variational formulation of problem (1)-(3) and describe how the finite element method implemented using an unstructured mesh can be used to solve the 2D TSRFDE on an arbitrarily convex domain and present a second order temporal numerical scheme. In section 4, we discuss the stability and convergence of the method. In section 5, we discuss the FEM for the 2D TSRFDE on a non-convex domain. In section 6, some numerical examples on irregular domains, including convex domains and a multiply-connected domain, are presented. We also solve a coupled fractional Bloch-Torrey equation to verify the effectiveness of the method. Finally, some conclusions of the work are drawn. 2. Preliminary knowledge At first, we introduce some definitions and lemmas on the fractional derivative space, which were first established by Roop and Ervin [12, 13] in the one-dimensional case. In the two-dimensional rectangular domains, these results are also applicable [37]. Here, we extend them to convex domains. Referring to Fig.1, we define Z d1 Z m(x) Z b1 Z r(y) (u, v)Ω := u(x, y)v(x, y)dxdy = u(x, y)v(x, y)dydx a1

c(y)

c1

5

g(x)

ACCEPTED MANUSCRIPT

1/2

and ||u||L2 (Ω) := (u, u)Ω . Definition 2.1. (Left fractional derivative space) For µ > 0, denote the semi-norm and the norm respectively as 1/2 1/2   ˜µ 2 ˜ µ u||2 , ||u||J˜µ (Ω) := ||u||2L2 (Ω) + |u|2J˜µ (Ω) , |u|J˜µ (Ω) := ||c(y) D x L2 (Ω) + ||g(x) Dy u||L2 (Ω) L

L

L

µ and define J˜Lµ (Ω) (J˜L,0 (Ω)) as the closure of C ∞ (Ω) (C0∞ (Ω)) with respect to || · ||J˜µ (Ω) . L

Definition 2.2. (Right fractional derivative space) For µ > 0, denote the semi-norm and the norm respectively as:

R

and define

µ J˜R (Ω)

R

µ (J˜R,0 (Ω))

CR IP T

 1/2  1/2 2 2 ˜ µ u||2 ˜ µ u||2 µ |u|J˜µ (Ω) := ||x D + || + |u| D , ||u|| := ||u|| , µ ˜ y ˜ L2 (Ω) L2 (Ω) L2 (Ω) J (Ω) r(y) m(x) J (Ω) R

as the closure of C ∞ (Ω) (C0∞ (Ω)) with respect to || · ||J˜µ (Ω) . R

Definition 2.3. (Symmetric fractional derivative space) For µ > 0, µ 6= n − the norm respectively as:

1 2,

n ∈ N denote the semi-norm and

 1/2  1/2 ˜ xµ u, x D ˜ µ u)Ω | + |(g(x) D ˜ yµ u,y D ˜ µ u)Ω | , ||u||J˜µ (Ω) := ||u||2L2 (Ω) + |u|2J˜µ (Ω) |u|J˜µ (Ω) := |(c(y) D , r(y) m(x) S

AN US

S

S

µ (Ω)) as the closure of C ∞ (Ω) (C0∞ (Ω)) with respect to || · ||J˜µ (Ω) . and define J˜Sµ (Ω) (J˜S,0 S

Definition 2.4. (Fractional Sobolev space) For µ > 0, denote the semi-norm and the norm respectively as:  1/2 |u|H µ (Ω) := || |ξ|µ F(ˆ u)(ξ)||L2 (R2 ) , ||u||H µ (Ω) := ||u||2L2 (Ω) + |u|2H µ (Ω) ,

M

where F(ˆ u)(ξ) is the Fourier transform of u ˆ, which is the zero extension of u outside Ω and define H µ (Ω) (H0µ (Ω)) ∞ ∞ as the closure of C (Ω) (C0 (Ω)) with respect to || · ||H µ (Ω) .

AC

CE

PT

ED

Define the following fractional derivative and integral operators: Z x ∂n 1 µ (x − s)n−µ−1 u(s, y, t) ds, −∞ Dx u(x, y, t) = Γ(n − µ) ∂xn −∞ Z +∞ (−1)n ∂ n µ (s − x)n−µ−1 u(s, y, t) ds, x D+∞ u(x, y, t) = Γ(n − µ) ∂xn x Z x 1 −µ (x − s)µ−1 u(s, y, t) ds, −∞ Ix u(x, y, t) = Γ(µ) −∞ Z (−1)n +∞ −µ (s − x)µ−1 u(s, y, t) ds, x I+∞ u(x, y, t) = Γ(µ) x Z x 1 −µ ˜ I u(x, y, t) := (x − s)µ−1 u(s, y, t) ds, c(y) x Γ(µ) c(y) Z (−1)n r(y) ˜−µ (s − x)µ−1 u(s, y, t) ds. x Ir(y) u(x, y, t) := Γ(µ) x

The definitions of the operators in the y direction are similar. Remark 2.1. If supp(u) ⊂ Ω, then ˜−µ . xI

µ −∞ Dx u

=

µ ˜µ c(y) Dx u, x D+∞ u

˜ µ u, = xD r(y)

−µ −∞ Ix u

=

˜−µ c(y) Ix u

−µ and x I+∞ u=

r(y)

Lemma 2.1. Let µ > 0, define the operators: (i) −∞ Ix−µ : L2 (Ω) → L2 (Ω), (ii) −∞ Dxµ : J˜Lµ (Ω) → L2 (Ω), (iii) −µ µ ˜µ x I+∞ : L2 (Ω) → L2 (Ω), (iv) x D+∞ : JR (Ω) → L2 (Ω), then all the operators are bounded operators.

6

ACCEPTED MANUSCRIPT

Proof. Using Young’s theorem [44] ||v ∗ w||L2 (Ω) ≤ ||v||L1 (Ω) ||w||L2 (Ω) and noting that denotes convolution, we have

−µ −∞ Ix u

=

xµ−1 Γ(µ)

∗ u, where ∗

1 1 ||xµ−1 ∗ u||L2 (Ω) ≤ ||xµ−1 ||L1 (Ω) ||u||L2 (Ω) Γ(µ) Γ(µ) Z b1 Z r(y) Z b1 Z d1 1 1 |x|µ−1 dxdy||u||L2 (Ω) |x|µ−1 dxdy||u||L2 (Ω) ≤ = Γ(µ) a1 c(y) Γ(µ) a1 c1

||−∞ Ix−µ u||L2 (Ω) =

(b1 − a1 )(|d1 |µ + |c1 |µ ) ||u||L2 (Ω) ≤ C||u||L2 (Ω) . Γ(µ + 1)

CR IP T

≤

By the definition of J˜Lµ (Ω), we have

1

||−∞ Dxµ u||L2 (Ω) ≤ (||u||2L2 (Ω) + ||−∞ Dxµ u||2L2 (Ω) + ||−∞ Dyµ u||2L2 (Ω) ) 2

1 ˜ µ u||2 ˜µ 2 2 = (||u||2L2 (Ω) + ||c(y) D x L2 (Ω) + ||g(x) Dy u||L2 (Ω) ) = ||u||J˜µ (Ω) . L

The proofs of (iii) and (iv) are similar.

AN US

µ µ Lemma 2.2. For u ∈ J˜L,0 (Ω) ∩ J˜R,0 (Ω) and 0 < s < µ, we have

˜ xµ u||L (Ω) , ˜ xs u||L (Ω) ≤ C2 ||c(y) D ||u||L2 (Ω) ≤ C1 ||c(y) D 2 2 s µ ˜ ˜ ||u||L (Ω) ≤ C3 ||g(x) D u||L (Ω) ≤ C4 ||g(x) D u||L (Ω) , y

2

y

2

2

where C1 , C2 , C3 and C4 are some positive constants independent of u. Proof. Combining Lemma 2.1, we have

M

˜ s u||L (Ω) , ||u||L2 (Ω) = ||−∞ Ix−s −∞ Dxs u||L2 (Ω) ≤ C1 ||−∞ Dxs u||L2 (Ω) = C1 ||c(y) D x 2

˜ xs u||L (Ω) = ||−∞ Dxs u||L (Ω) = ||−∞ Ix−(µ−s) −∞ Dxµ u||L (Ω) ||c(y) D 2 2 2 ˜ µ u||L (Ω) . ≤ C2 ||−∞ Dµ u||L (Ω) = C2 ||c(y) D x

2

2

ED

x

The second inequality can be proved similarly.

PT

µ Lemma 2.3. If µ > 0, then J˜Lµ (Ω), J˜R (Ω) and H µ (Ω) are equivalent with equivalent norms and semi-norms; if µ µ µ 1 ˜ µ > 0, µ 6= n − 2 , n ∈ N, then JL,0 (Ω), J˜R,0 (Ω), J˜S,0 (Ω) and H0µ (Ω) are equivalent with equivalent norms and semi-norms.

Proof. The proof is similar to the one-dimensional case in [12], therefore we omit it here.

CE

2µ 2µ Lemma 2.4. If µ ∈ (1/2, 1), u, v ∈ J˜L,0 (Ω) ∩ J˜R,0 (Ω), then

AC





˜ 2µ c(y) Dx u(x, y), v(x, y)



=



˜µ ˜µ c(y) Dx u(x, y), x Dr(y) v(x, y)



, Ω  Ω ˜ 2µ ˜ µ u(x, y), c(y) D ˜ xµ v(x, y) . = xD x Dr(y) u(x, y), v(x, y) r(y) Ω

Ω

Proof. Combining the formula [12] (a Ix−µ w, v)L2 (a,b) = (w, x Ib−µ v)L2 (a,b) , we have     ˜ 2µ = Dx2 c(y) I˜x−(2µ−2) u(x, y), v(x, y) c(y) Dx u(x, y), v(x, y) Ω Ω     −(2µ−2) −(2µ−2) ˜ ˜ = Dxc(y) Ix u(x, y), −Dx v(x, y) = c(y) Ix Dx u(x, y), −Dx v(x, y) Ω    Ω −(µ−1) ˜ xµ u(x, y), x D ˜ µ v(x, y) . = c(y) I˜x−(µ−1) Dx u(x, y), −x I˜r(y) Dx v(x, y) = c(y) D r(y) Ω

The proof of the second identity is similar.

7

Ω

ACCEPTED MANUSCRIPT

3. Variational formulation 3.1. The fully discrete variational formulation For convenience, in the subsequent sections, we suppose that C, C1 , C2 , ... are some positive constants, which may take distinct values according to different contexts discussed throughout this paper. T be the time step and tn = nτ , n = 0, 1, 2, ..., N . Using the finite difference method we have Let τ = N 1 = Γ(1 − γ)

Z

tn

0

n

(tn − s)−γ

∂u(x, y, s) τ −γ X n ds = bk u(x, y, tk ) + Rtn , ∂s Γ(2 − γ) k=0

CR IP T

C γ 0 Dt u(x, y, tn )

where bnn = 1, bn0 = (n − 1)1−γ − n1−γ < 0, bnk = (n − k + 1)1−γ − 2(n − k)1−γ + (n − k − 1)1−γ < 0, k = 1, 2, ..., n − 1. Denote n

∇γt u(x, y, tn ) =

τ −γ X n bk u(x, y, tk ), Γ(2 − γ)

n = 1, 2, ..., N,

k=0

then [45]

(7)

AN US

γ γ 2−γ ||Rtn ||0 = ||C . 0 Dt u(x, y, tn ) − ∇t u(x, y, tn )||0 ≤ Cτ

(6)

Denote V = H0α (Ω) ∩ H0β (Ω). We divide the domain Ω into a number of regular triangular regions. Let Th be this triangulation and h be the maximum diameter of the triangular elements. We define the finite element subspace as: n o Vh := vh |vh ∈ C(Ω) ∩ V, vh |K is linear for all K ∈ Th and vh |∂Ω = 0 .

M

Assume that unh ∈ Vh is the approximation of u(x, y, t) at t = tn . Then, by Lemma 2.4, we obtain the fully discrete scheme associated with the variational form of Eq.(1) is: find unh ∈ Vh for each t = tn (n = 1, 2, ..., N ) such that ß (∇γt unh , vh )Ω + A(unh , vh )Ω = (f n , vh )Ω , ∀vh ∈ Vh , (8) u0h = u0h , K1 2 cos(απ) ,

Ky =

K2 2 cos(βπ)

and

 o   n ˜ xα v) ˜ α v + xD ˜ α u, c(y) D ˜ xα u, x D := Kx c(y) D r(y) r(y) n Ω   Ωo β β β β ˜ ˜ ˜ ˜ +Ky g(x) Dy u, y Dm(x) v + y Dm(x) u, g(x) Dy v .

PT

A(u, v)Ω

ED

where u0h ∈ Vh is a reasonable approximation of u0 and Kx =

Ω

Ω

AC

CE

3.2. The implementation of FEM with an unstructured mesh We consider the computation process for piecewise linear polynomials on the triangular element ep , p = 1, 2, ..., Ne , where Ne is the total number of triangles. Then, within element ep , the field function up (x, y) can be written as up (x, y) =

3 X

uj0 ϕj0 (x, y),

j0 =1

where the triangle vertices are numbered in a counter-clockwise order as 1, 2, 3 and the basis function ϕj0 (x, y) is defined by 1 ϕj0 (x, y) = (aj0 x + bj0 y + cj0 ), ϕj0 (x, y) = 0, 2∆ep (x,y)∈e / p (x,y)∈ep a1 = y2 − y3 , a2 = y3 − y1 , a3 = y1 − y2 ,

b1 = x3 − x2 , b2 = x1 − x3 , b3 = x2 − x1 ,

c1 = x2 y3 − x3 y2 , c2 = x3 y1 − x1 y3 , c3 = x1 y2 − x2 y1 , 8

ACCEPTED MANUSCRIPT

where ∆ep is the area of triangle element p. It is well-known that ϕj0 (xi0 , yi0 ) = δi0 j0 ,

i0 , j0 = 1, 2, 3,

where δ is the Kronecker function. With these local field functions and basis functions, we can obtain the formulation of u(x, y) for the whole triangulation: u(x, y) =

Np X

ui li (x, y),

CR IP T

i=1

where li (x, y) is the new basis function whose support domain is Ωei (see figure 3) and Np is the total number of vertices on the convex domain Ω. Now, we rewrite unh in the form of unh

=

Np X

uni li (x, y),

i=1

(9)

Np X i=1

AN US

where uni are the coefficients that are to be solved for. Substituting Eq.(9) into Eq.(8) with vh = lj (x, y), j = 1, 2, . . . , Np gives Np n−1 h i X X uni (li , lj )Ω + ω0 A(li , lj )Ω = − bnk uki (lk , lj )Ω − bn0 (u0 , lj )Ω + ω0 (f n , lj )Ω ,

(10)

i=1 k=1

where ω0 = τ γ Γ(2 − γ). Eq.(10) can be expressed in matrix form as n−1 X

M

(M + ω0 B)U n = −M

k=1

bnk U k − G0 + ω0 F1n ,

(11)

ED

P e where M is the mass matrix with elements Mij = (lj , li )Ω = N p=1 (lj , li )ep , B is the stiffness matrix with elements n n n n T 0 Bij = A(lj , li )Ω and U = [u1 , u2 , ..., uNp ] . The vectors G and F1n are given by G0 = bn0 [(u0 , l1 )Ω , (u0 , l2 )Ω , ..., (u0 , lNp )Ω ]T =

Ne X

bn0 [(u0 , l1 )ep , (u0 , l2 )ep , ..., (u0 , lNp )ep ]T ,

PT

p=1

F1n = [(f n , l1 )Ω , (f n , l2 )Ω , ..., (f n , lNp )Ω ]T =

Ne X

[(f n , l1 )ep , (f n , l2 )ep , ..., (f n , lNp )ep ]T .

p=1

AC

CE

Due to the non-local property of the fractional derivative, matrix B is the most difficult part to calculate. For matrix B, the (i, j) entry is given by n    o ˜ xα lj (x, y), x D ˜ α li (x, y) + x D ˜ α lj (x, y), c(y) D ˜ xα li (x, y) Bij = A(lj , li )Ω = Kx c(y) D r(y) r(y) Ω n Ω   o β β β β ˜ ˜ ˜ ˜ + Ky g(x) Dy lj (x, y), y Dm(x) li (x, y) + y Dm(x) lj (x, y), g(x) Dy li (x, y) . (12) Ω

Ω

˜ xα l(x, y), x D ˜ α l(x, y), g(x) D ˜ yβ l(x, y), Here, we first discuss the support domains of the four fractional derivatives c(y) D r(y) β L R D U ˜ l(x, y), and denote them as Ωe , Ωe , Ωe , Ωe , respectively. We adopt the notation that L, R, D, U correspond yD m(x)

i

i

i

i

to the ’left’, ’right’, ’down’ and ’up’ directions for the four fractional derivatives. They do not correspond to the actual location in the domain Ω. From figure 3, we see node i is surrounded by the elements e1 , e2 , e3 , e4 , e5 and the boundary ∂Ωei is constituted by ∂Ωlei and ∂Ωrei , where ∂Ωlei is made up of line segments B1 B2 , B2 B3 , B3 B4 and ∂Ωrei is connected by the line segments B4 B5 , B5 B1 , respectively. Then ∂Ωlei , y = yB4 (x ≥ xB4 ), y = yB1 (x ≥ xB1 ) R D U and ∂Ω form the support domain ΩL ei . Similarly, we can obtain the support domains Ωei , Ωei and Ωei (see figure 3). In view of the similarity of the four terms in the right hand side of Eq.(12), we illustrate the computation of 9

ACCEPTED MANUSCRIPT

ΩD ei

B1 B2 e1 e2 Ωeei5 e B3 e3 4 B5

∂Ω

B1 B2 e1 e2 Ωeei5 e B3 e3 4 B5

∂Ω

B4

B4

ΩLei

CR IP T

ΩR ei

∂Ω

ΩUei

∂Ω

U D R Figure 3: The illustration of support domains ΩL ei , Ωei , Ωei , Ωei .



˜α ˜α c(y) Dx lj , x Dr(y) li



Ω

=

Ne  X

AN US

˜ α li )Ω as an example. By applying Gauss quadrature, we have ˜ xα lj , x D (c(y) D r(y) ˜α ˜α c(y) Dx lj , x Dr(y) li

p=1

≈

Ne X

X

ωq

p=1 (xˆq ,yˆq )∈GK



=

ep

Ne Z X p=1



˜α c(y) Dx lj

(xˆq ,yˆq )

ep

˜α ˜α c(y) Dx lj x Dr(y) li

˜α x Dr(y) li

(xˆq ,yˆq )

dxdy

,

M

where GK stands for the set of all Gauss points in element ep and ωq are the weights associated with the Gauss point R ˜α P (xˆq , yˆq ) (see figure 4). When point P (xˆq , yˆq ) is out of the support domains ΩL ei and Ωei , we have c(y) Dx lj (x, y) = 0 ˜ xα lj ˜ α li (x, y) = 0. To evaluate c(y) D , suppose that segment y = yˆq , c(yˆq ) ≤ x ≤ xˆq intersects nq and x D r(y)

(xˆq ,yˆq )

n

PT

ED

points with the triangular element of Ωej , and these points are numbered as x0q < x1q < x2q <, ..., xq q , then by the definition of the fractional derivative, we have Z x Å ã 1 ∂ α α −α ˜ ˜ D l = D l (x, y ˆ ) = (x − ξ) l (ξ, y ˆ )dξ q j q c(y) x j c(yˆq ) x j Γ(1 − α) ∂x c(yˆq ) (xˆq ,yˆq ) x=xˆq x=xˆq k Z Å ã n i x X i 1 ∂ . (13) = (x − ξ)−α lj (ξ, yˆq )dξ Γ(1 − α) ∂x xk−1 x=xˆq i k=1

AC

CE

Then, there are three different cases that need to is only located in ΩL ej , and we have       lj (x, y0 ) =     

be discussed (see figure 4). In case I, the Gauss point P0 (x0 , y0 ) 0, ϕj2 (x, y0 ), ϕj1 (x, y0 ), ϕj5 (x, y0 ), 0,

x00 x10 x20 x30

≤ x ≤ x10 , ≤ x ≤ x20 , ≤ x ≤ x30 , ≤ x ≤ x40 , x40 ≤ x,

where ϕjp (x, y) is the basis function of node j on element p and c(y0 ) = x00 . Case II, the Gauss point P1 (x1 , y1 ) is R L α ˜ x lj only located in the Ωej \Ωej , then we have c(y) D = 0. Case III, the Gauss point P2 (x2 , y2 ) is located in (x1 ,y1 )

R ΩL ej ∩ Ωej = Ωej , then we have

   

0, ϕj2 (x, y2 ), lj (x, y2 ) = ϕ (x, y2 ),    j3 ϕj4 (x, y2 ), 10

x02 x12 x22 x32

≤ x ≤ x12 , ≤ x ≤ x22 , ≤ x ≤ x32 , ≤ x ≤ x42 .

ACCEPTED MANUSCRIPT

Case I

x00

Case II

x10 x20 x30 x40

e1 e P0 (x0 , y0 ) 5 j e e4

Case III

P2 (x2 , y2 )

e1 e 5 j e e4

e2

e2

3

x02

3

x31

P1 (x1 , y1 ) x0 x1 x21 1 1

ΩLej

ΩLej

ΩR ej

x52

ΩLej

ΩR ej

CR IP T

ΩR ej

e1 e

j 5 x12 e2 e e 4 x22 3 x32 x42

Figure 4: The points of intersection by y = yˆq with the triangle element of Ωej and ∂Ω.

II,

In case III,

(xˆq ,yˆq )

˜ α lj , we also need to consider three cases. In case I, x D r(y)    

(x0 ,y0 )

= 0. In case

0, x ≤ x01 , ϕj3 (x, y1 ), x01 ≤ x ≤ x11 , lj (x, y1 ) = ϕj4 (x, y1 ), x11 ≤ x ≤ x21 ,    0, x21 ≤ x ≤ x31 .

AN US

˜ α lj Similarly, to evaluate x D r(y)

M

 ϕj2 (x, y2 ),    ϕj3 (x, y2 ), lj (x, y2 ) = ϕ (x, y2 ),    j4 0,

x12 x22 x32 x42

≤ x ≤ x22 , ≤ x ≤ x32 , ≤ x ≤ x42 , ≤ x ≤ x52 .

ED

As lj (ξ, yˆq ) is a linear function on [xk−1 , xki ], k = 1, 2, ..., ni , Eq.(13) can be evaluated using integration by parts. i Finally, we summarise the whole computation process in the following algorithm (see Algorithm 1). Algorithm 1 Compiling fractional derivative using FEM on an unstructured mesh

6: 7:

(xˆq ,yˆq )

(xˆq ,yˆq )

end for Form stiffness matrix B; end for Calculate (lj , li )ep on each triangle element ep to form the mass matrix M ; Calculate (u0 , lk )ep and (f n , lk )ep , k = 1, 2, ..., Np and obtain G0 , F n ; Solve the linear system (11) and obtain U n .

AC

8:

PT

2: 3: 4: 5:

Partition the convex domain Ω with unstructured triangular elements ep and save the element information (node number, coordinates, and element number ); for p = 1, 2, · · · , Ne do Find the Gauss points (xˆq , yˆq ) and weights ωi on each triangle element ep ; for j = 1, 2, · · · , Np do Find the support domain Ωej ; ˜ xα lj ˜ α li , xD ; Find the points of intersection by y = yˆq with Ωej and calculate c(y) D r(y) ( xˆq ,yˆq ) (xˆq ,yˆq ) ˜ β li ˜ yβ lj Find the points of intersection by x = xˆq with Ωej and calculate g(x) D , yD ; m(x)

CE

1:

9: 10: 11: 12: 13:

3.3. Second order temporal numerical scheme In this part, we will give a second order temporal numerical scheme. Firstly, when 0 < γ < 1, for a function W (t) for which the Caputo fractional derivative of order γ exists, we have the following relationship between its

11

ACCEPTED MANUSCRIPT

Caputo fractional derivative and Riemann-Liouville fractional derivative [10] C γ 0 Dt W (t)

W (0) . Γ(1 − γ)tγ

= 0 Dtγ [W (t) − W (0)] = 0 Dtγ W (t) −

(14)

Then, we consider the shifted Gr¨ unwald formula [46] to approximate the Riemann-Liouville fractional derivative. Define

(γ)

where gk = (−1)k lemma [47, 48]


γ k

∞ X

k=0

(γ)

gk W (t − (k − l)τ ),

CR IP T

Aγτ,l W (t) = τ −γ

for k ≥ 0. To obtain the second order temporal numerical scheme, we need the following

Lemma 3.1. Let W (t) ∈ L1 (R), −∞ Dtγ+2 W and its Fourier transform belong to L1 (R), and define the weighted and shifted Gr¨ unwald difference operator by RL

∇γt W (t) =

2p − γ γ γ − 2q γ Aτ,p W (t) + A W (t), 2(p − q) 2(p − q) τ,q

RL

∇γt W (t) = −∞ Dtγ W (t) + O(τ 2 ),

for t ∈ R, where p and q are integers and p 6= q. By choosing (p, q) = (0, −1) in (15), we obtain RL

∇γt W (tn ) =

AN US

then we have

(15)

n n−1 γ X (γ) 2 + γ X (γ) g W (t ) − gk W (tn−k−1 ) n−k k 2τ γ 2τ γ k=0

k=0

M

n n 1 X (γ) 1 X (γ) = γ ωk W (tn−k ) = γ ωn−k W (tk ), τ τ k=0

=

2+γ (γ) 2 g0 ,

(γ)

ωk

=

2+γ (γ) 2 gk

(γ)

− γ2 gk−1 , k ≥ 1. Now, we transform (1) into the following form

ED

(γ)

where ω0

k=0

γ 0 Dt u(x, y, t)

= K1

∂ 2α u(x, y, t) ∂|x|

2α

+ K2

∂ 2β u(x, y, t) 2β

∂|y|

+ f˜(x, y, t),

(16)

CE

PT

u(x,y,0) where f˜(x, y, t) = f (x, y, t) + Γ(1−γ)t γ . Then, we obtain the fully discrete scheme associated with the variational n form of Eq.(16) is: to find uh ∈ Vh for each t = tn (n = 1, 2, ..., N ) such that ß RL γ n ( ∇t uh , vh )Ω + A(unh , vh )Ω = (f˜n , vh )Ω , ∀vh ∈ Vh , (17) u0h = u0h ,

Substituting Eq.(9) into Eq.(17) with vh = lj (x, y), j = 1, 2, . . . , Np gives

AC

Np X i=1

uni

Np n−1 h i X X (γ) (γ) γ ω0 (li , lj )Ω + τ A(li , lj )Ω = − ωn−k uki (lk , lj )Ω + τ γ (f˜n , lj )Ω .

(18)

i=1 k=0

Writing (18) in matrix form, we have (γ)

(ω0 M + τ γ B)U n = −M

n−1 X

(γ)

ωn−k U k + τ γ F2n ,

(19)

k=0

where F2n = [(f˜n , l1 )Ω , (f˜n , l2 )Ω , ..., (f˜n , lNp )Ω ]T . Remark 3.1. When the solution u(x, y, t) is not smooth enough, the modified weighted shifted Gr¨ unwald-Letnikov formula with appropriate correction terms can be used [49]. 12

ACCEPTED MANUSCRIPT

4. Stability and convergence of the fully discrete scheme Before giving the proof, we introduce some new definitions and lemmas. Here, we simplify the notations (·, ·)Ω , || · ||L2 (Ω) and || · ||H s (Ω) as (·, ·), || · ||0 , || · ||s , respectively. Let σ = max{α, β}, we define the semi-norm | · |(α,β) and norm ||| · |||(α,β) as follows  1 ˜ α u||2 + K2 ||g(x) D ˜ β u||2 2 , |u|(α,β) := K1 ||c(y) D x 0 y 0

Then we have the following lemma.

 12  |||u|||(α,β) := ||u||20 + |u|2(α,β) .

CR IP T

Lemma 4.1. Assume that u ∈ H0α (Ω) ∩ H0β (Ω) and Ω is a convex domain, then C1 |||u|||(α,β) ≤ |u|(α,β) ≤ |||u|||(α,β) ≤ C2 |u|H σ (Ω) ,

where positive constants C1 < 1 and C2 are independent of u, i.e., the semi-norm |u|(α,β) and ||u||(α,β) are equivalent. Proof. We have that |u|(α,β) ≤ |||u|||(α,β) . By Lemma 2.2, ˜ xs u||0 ≤ C|u|(α,β) , ||u||0 ≤ C||c(y) D

AN US

then |||u|||2(α,β) ≤ (C 2 + 1)|u|2(α,β) , therefore,

C1 |||u|||(α,β) ≤ |u|(α,β) ,

C1 = √

1 . C2 + 1

Since |u|(α,β) ≤ C3 |u|JLσ (Ω) , by Lemma 2.3, we obtain

p p C 2 + 1|u|(α,β) ≤ C3 C 2 + 1|u|JLσ (Ω) ≤ C2 |u|H σ (Ω) .

M

|||u|||(α,β) ≤

By the above lemma, it is straightforward to obtain, ∀ u ∈ H0α (Ω) ∩ H0β (Ω),

ED

A(u, v) ≤ C|||u|||(α,β) |||v|||(α,β) ,

A(u, u) ≥ C|||u|||2(α,β) .

Here we choose the interpolation operator Ih to satisfy the approximation properties of the subspace of H s+1 (Ω) [50], then Ih : H s+1 (Ω) → Vh satisfies

PT

||u − Ih u||l ≤ Chµ−l ||u||µ ,

∀u ∈ H µ (Ω), 0 ≤ l < µ ≤ s + 1.

CE

We define the projection operator Ph : V → Vh satisfying u ∈ V, ∀v ∈ Vh .

A(Ph u, v) = A(u, v),

Furthermore, we can derive the approximation property of Ph .

AC

Lemma 4.2. If u ∈ H 2 (Ω) ∩ V , σ = max(α, β), then |u − Ph u|(α,β) ≤ Ch2−σ ||u||2 ,

in which the constant C is independent of u and h. Proof. Since |u − Ph u|2(α,β) = A(u − Ph u, u − Ph u) = A(u − Ph u, u − Ih u) and A(u − Ph u, u − Ih u) ≤ C|||u − Ph u|||(α,β) |||u − Ih u|||(α,β) ≤ C1 |u − Ph u|(α,β) |u − Ih u|(α,β) , therefore, via the approximation properties (20) and Lemma 4.1, we have |u − Ph u|(α,β) ≤ C1 |u − Ih u|(α,β) ≤ C2 |u − Ih u|σ ≤ C3 h2−σ ||u||2 .

13

(20)

ACCEPTED MANUSCRIPT

As A(u, v) is continuous and coercive, we can derive the existence and uniqueness of scheme (8). Now, we discuss stability and the convergence of this scheme. Theorem 4.3. The fully discrete variational scheme (8) is unconditionally stable. Proof. Assume that zhn (n = 1, 2, ..., N ) is another solution of the fully scheme (8), and let Ehn = unh − zhn , then (∇γt Ehn , vh ) + A(Ehn , vh ) = 0, i.e., n−1 X

bnk Ehk , vh ).

k=0

CR IP T

(Ehn , vh ) + τ γ Γ(2 − γ)A(Ehn , vh ) = −(

Using the Cauchy-Schwarz inequality and taking vh = Ehn and noting the positivity of A(·, ·), we have ||Ehn ||0 ≤ −

n−1 X k=0

bnk ||Ehk ||0 .

Pn−1

Utilising mathematical induction and noticing − k=0 bnk = 1, it is readily concluded that ||Ehn ||0 ≤ ||Eh0 ||0 . Therefore, the fully discrete scheme (8) is unconditional stable.

AN US

Theorem 4.4. Suppose that u(tn ), unh are the exact solution and numerical solution of problem (1)-(3) at t = tn γ ∞ s+1 respectively, and u, utt , C (Ω)). Then for n = 1, 2, ..., N , the following error estimate holds 0 Dt u ∈ L (0, T ; H n
o γ 2 |||unh − u(tn )|||2(α,β) ≤ C τ 4−2γ + ||u0h − u(t0 )||2σ + h2s+2−2σ ||u(t0 )||2s+1 + ||u(tn )||2s+1 + max ||C . 0 Dt u(t)||s+1 ∀v ∈ Vh .

(21)

M

Proof. As u(x, y, t) is the exact solution of problem (1)-(3) , then     C γ n 0 Dt u(tn ), vh + A u(tn ), vh = (f , vh ),

0≤t≤T

ED

Denote en = u(tn ) − unh , and subtract (8) from (21), we have   γ γ (∇γt en , vh ) + A(en , vh ) = − C 0 Dt u(tn ) − ∇t u(tn ), vh .

Define en = ρn + θn , ρn = u(tn ) − Ph u(tn ), θn = Ph u(tn ) − unh and let vh = ∇γt θn , we obtain (22)

PT

(∇γt θn , ∇γt θn ) + A(θn , ∇γt θn ) = −(∇γt ρn , ∇γt θn ) − (Rtn , ∇γt θn ). 1 1 ||∇γt ρn ||20 + ||∇γt θn ||20 , 2 2 1 1 |(Rtn , ∇γt θn )| ≤ ||Rtn ||20 + ||∇γt θn ||20 . 2 2

|(∇γt ρn , ∇γt θn )| ≤

CE

We have

Substituting them into (22) leads to (23)

Å X ã ï ò n n−1 n−1 X X 1 A θn , bnk θk = A(θn , θn ) + bnk A(θk , θk ) − bnk A(θk − θn , θk − θn ) . 2

(24)

AC

Å X ã n τ γ Γ(2 − γ) τ γ Γ(2 − γ) n 2 A θn , bnk θk ≤ ||∇γt ρn ||20 + ||Rt ||0 . 2 2

Note that

k=0

k=0

k=0

k=0

As A(θk − θn , θk − θn ) ≥ 0, bnk < 0, k = 0, 1, ..., n − 1, substituting (24) into (23), gives A(θn , θn ) ≤ −

n−1 X k=0

bnk A(θk , θk ) + τ γ Γ(2 − γ)||∇γt ρn ||20 + τ γ Γ(2 − γ)||Rtn ||20 . 14

(25)

ACCEPTED MANUSCRIPT

From (7) we have γ n γ n 2−γ ||Rtn ||20 ≤ Cτ 4−2γ max ||utt ||20 , ||C max ||ρtt ||0 , 0 Dt ρ − ∇t ρ ||0 ≤ Cτ 0≤t≤T

0≤t≤T

then γ n 2 4−2γ ||∇γt ρn ||20 ≤ ||C max ||ρtt ||20 . 0 Dt ρ ||0 + Cτ 0≤t≤T

CR IP T

Using Lemma 2.2 and Lemma 4.2, we obtain ™ ß γ n 2 4−2γ . ||∇γt ρn ||20 ≤ C h2s+2−2σ ||C D ρ || + τ t 0 s+1

By Lemma 4.1, Lemma 4.2 and A(u, u) ≤ C||u||2(α,β) , we have n o n o A(θ0 , θ0 ) ≤ C||θ0 ||2(α,β) ≤ C ||u(t0 ) − u0h ||2(α,β) + ||ρ0 ||2(α,β) ≤ C ||u(t0 ) − u0h ||2σ + h2s+2−2σ ||u(t0 )||2s+1 . Then, from (25), we obtain

n o A(θn , θn ) ≤ Γ(2 − γ)(−bn0 )C ||u(t0 ) − u0h ||2σ + h2s+2−2σ ||u(t0 )||2s+1

AN US

ß ™ n−1 X γ 2 + Γ(2 − γ)τ γ C τ 4−2γ + h2s+2−2σ ||C bnk A(θk , θk ). 0 Dt u(tn )||s+1 −

(26)

k=1

M

Utilising mathematical induction, it is readily concluded that n o A(θn , θn ) ≤ Γ(1 − γ)C1 ||u(t0 ) − u0h ||2σ + h2s+2−2σ ||u(t0 )||2s+1 ß ™ γ 2 + Γ(1 − γ)nγ τ γ C2 τ 4−2γ + h2s+2−2σ max ||C D u(t)|| t 0 s+1 , 0≤t≤T

ED

Since A(θn , θn ) ≥ C||θn ||2(α,β) , the above mathematical induction gives ß h i™ n 2 4−2γ 0 2 2s+2−2σ 2 C γ 2 ||θ ||(α,β) ≤ C τ + ||u(t0 ) − uh ||σ + h ||u(t0 )||s+1 + max ||0 Dt u(t)||s+1 ,

and

0≤t≤T

(27)

PT

  ||en ||2(α,β) ≤ C ||ρn ||2(α,β) + ||θn ||2(α,β) .

CE

Therefore, using Lemma 4.2, we obtain n
o γ 2 ||unh − u(tn )||2(α,β) ≤ C τ 4−2γ + ||u(t0 ) − u0h ||2σ + h2s+2−2σ ||u(t0 )||2s+1 + ||u(tn )||2s+1 + max ||C . 0 Dt u(t)||s+1 0≤t≤T

AC

Remark 4.1. It can deduced from Theorem 4.4 that when a triangular linear basis function (s = 1) is used, the error satisfies |||unh − u(tn )|||(α,β) ≤ C(τ 2−γ + h2−σ ).

Remark 4.2. It is well-known that the study of the regularities of the elliptic space-fractional PDEs is still under development. Wang and Yang [51] analyzed the regularity of the solution of fractional-order derivatives involving a variable coefficient in H¨ older spaces, and established the well-posedness of a Petrov-Galerkin formulation. Jin et al. [52] developed variational formulations for boundary value problems involving either Riemann-Liouville or Caputo fractional derivatives and established the Sobolev regularity of the variational solutions. Recently, Ervin et al. [53] investigated the regularity solution of the steady-state fractional diffusion equation on a bounded domain. As the theoretical analysis of regularities for the elliptic space-fractional PDEs on convex domains with both left and right sided fractional derivatives is challenging, it is still an open research problem and needs further exploration. 15

ACCEPTED MANUSCRIPT

5. FEM for the 2D TSRFDE on non-convex domains In this section, we will discuss the FEM for the two-dimensional time-space Riesz fractional diffusion equation on non-convex domains. When the considered solution domain is nonconvex, the definition of the fractional derivative should change according to the shape of theSsolution domain. For example, we consider the multiply-connected domain Ω0 = [a, b] × [c, d]\[a∗ , b∗ ] × [c∗ , d∗ ] = 8i=1 Ωi shown in Figure 5. In this case, the fractional derivative must be modified to take the central hole into consideration.

Ω1

Ω2

Ω3

d∗ Ω4

Ω5

c

∗

Ω6

Ω7

Ω8

c a∗

b∗

b

AN US

a

CR IP T

d

Figure 5: The illustration of a multiply-connected domain.

M

Definition 5.1. The definition of the fractional derivative in the x direction for the domain Ω0 illustrated in figure 5 is defined as: ß 2α x ∈ Ω0 \Ω5 , a Dx u(x, y, t), 2α a Dx u(x, y, t) = 2α x ∈ Ω5 . b∗ Dx u(x, y, t), ß 2α x Db u(x, y, t), x ∈ Ω0 \Ω4 , 2α x Db u(x, y, t) = 2α x ∈ Ω4 . D x a∗ u(x, y, t),

ED

The definition of the fractional derivative in the y direction is defined in a similar manner.

CE

PT

Following Algorithm 1, we can solve the problem on a multiply-connected domain Ω0 . We now provide the theoretical analysis when the solution is nonconvex. Here we still take the multi-connected domain Ω0 shown in S figure 5 as an example. The domain Ω0 can be divided into a finite union of convex domains Ωi : Ω0 = 8i=1 Ωi . Next, we define the inner product and its induced norm as: Ã 8 8 X X 1/2 (u, v)Ωi , ||u||L2 (Ω0 ) := (u, u)Ω0 = (u, v)Ω0 := ||u||2L2 (Ωi ) , i=1

i=1

AC

where (u, v)Ωi and ||u||L2 (Ωi ) are the general inner product and norm defined on convex domain Ωi . Then we define the semi-norm and the norm of the left fractional derivative space on Ω0 as: 1/2  1/2  ˜ xµ u||2 ˜µ 2 |u|J˜µ (Ω0 ) := ||a D , ||u||J˜µ (Ω0 ) := ||u||2L2 (Ω0 ) + |u|2J˜µ (Ω0 ) . L2 (Ω0 ) + ||c Dy u||L2 (Ω0 ) L

L

L

It is straightforward to derive

à |u|J˜µ (Ω0 ) = L

8 X i=1

à |u|2J˜µ (Ω ) , ||u||J˜µ (Ω0 ) = L

i

L

8 X i=1

||u||2J˜µ (Ω ) , L

i

where |u|J˜µ (Ωi ) and ||u||J˜µ (Ωi ) are defined in definition 2.1. Similarly, we can define other norm operators for the L L fractional derivative space on domain Ω0 . In the following, we will prove Lemmas 2.1 to 2.4 also hold on domain Ω0 . 16

ACCEPTED MANUSCRIPT

Lemma 5.1. Let µ > 0, define the operators: (i) −∞ Ix−µ : L2 (Ω0 ) → L2 (Ω0 ), (ii) −∞ Dxµ : J˜Lµ (Ω0 ) → L2 (Ω0 ), (iii) −µ µ ˜µ x I+∞ : L2 (Ω0 ) → L2 (Ω0 ), (iv) x D+∞ : JR (Ω0 ) → L2 (Ω0 ), then all the operators are bounded operators. Proof. By the definition we have

à ||−∞ Ix−µ u||L2 (Ω0 ) =

8 X i=1

||−∞ Ix−µ u||2L2 (Ωi ) .

−µ −∞ Ix

i=1

i=1

i=1

where C0 = max{Ci }, i = 1, 2, . . . , 8. Then we obtain

||−∞ Ix−µ u||L2 (Ω0 ) ≤ C||u||L2 (Ω0 ) .

By the definition of J˜Lµ (Ω0 ), we have

CR IP T

According to Lemma 2.1, the operator : L2 (Ωi ) → L2 (Ωi ) is bounded. Then we have à à à 8 8 8 X X X p ||−∞ Ix−µ u||2L2 (Ωi ) ≤ Ci ||u||2L2 (Ωi ) ≤ C0 ||u||2L2 (Ωi ) = C0 ||u||L2 (Ω0 ) ,

1

AN US

||−∞ Dxµ u||L2 (Ω0 ) ≤ (||u||2L2 (Ω0 ) + ||−∞ Dxµ u||2L2 (Ω0 ) + ||−∞ Dyµ u||2L2 (Ω0 ) ) 2

1

˜µ 2 ˜ xµ u||2 2 = (||u||2L2 (Ω0 ) + ||c(y) D L2 (Ω0 ) + ||g(x) Dy u||L2 (Ω0 ) ) = ||u||J˜µ (Ω0 ) .

The proofs of (iii) and (iv) are similar.

L

µ µ Lemma 5.2. For u ∈ J˜L,0 (Ω0 ) ∩ J˜R,0 (Ω0 ) and 0 < s < µ, we have

M

˜ xµ u||L (Ω ) , ˜ xs u||L (Ω ) ≤ C2 ||c(y) D ||u||L2 (Ω0 ) ≤ C1 ||c(y) D 2 0 2 0 ˜ yµ u||L (Ω ) , ˜ ys u||L (Ω ) ≤ C4 ||g(x) D ||u||L2 (Ω0 ) ≤ C3 ||g(x) D 2 0 2 0 where C1 , C2 , C3 and C4 are some positive constants independent of u.

˜ xs u||L (Ω ) = ||c(y) D 2 0

i=1

8 X

Ã

˜ xs u||2 ||c(y) D L2 (Ωi ) ≤

i=1

i=1

8 X i=1

˜ xµ u||2 ˜µ C2i ||c(y) D L2 (Ωi ) ≤ C2 ||c(y) Dx u||L2 (Ω0 ) ,

√ where C1 = max{ C1i } and C2 = max{ C2i }, i = 1, 2, . . . , 8. The second inequality can be proved similarly.

CE

√

Ã

PT

i=1

ED

Proof. Combining Lemma 2.2, we have à à à 8 8 8 X X X ˜s ˜ xs u||2 ˜ xs u||2 ||u||L2 (Ω0 ) = ≤ C ||u||2L2 (Ωi ) ≤ C1i ||c(y) D ||c(y) D 1 L2 (Ωi ) L2 (Ωi ) = C1 ||c(y) Dx u||L2 (Ω0 ) ,

AC

2µ 2µ Lemma 5.3. If µ ∈ (1/2, 1), u, v ∈ J˜L,0 (Ω0 ) ∩ J˜R,0 (Ω0 ), then     ˜ 2µ ˜ µ u(x, y), x D ˜ µ v(x, y) = c(y) D , c(y) Dx u(x, y), v(x, y) x r(y) Ω Ω   0   0 ˜ 2µ ˜ µ u(x, y), c(y) D ˜ xµ v(x, y) = xD . x Dr(y) u(x, y), v(x, y) r(y) Ω0

Proof. Combining Lemma 2.4, we have 

 ˜ 2µ c(y) Dx u(x, y), v(x, y)

à =

8  X

Ω0

=

Ω0

à 8  X i=1



˜µ ˜µ c(y) Dx u(x, y), x Dr(y) v(x, y)

i=1



˜ 2µ c(y) Dx u(x, y), v(x, y)

Ωi

The proof of the second identity is similar. 17

=



Ωi

˜µ ˜µ c(y) Dx u(x, y), x Dr(y) v(x, y)



Ω0

.

ACCEPTED MANUSCRIPT

Then we can obtain the fully discrete scheme associated with the variational form of Eq.(1) on the multiconnected domain Ω0 : find unh ∈ Vh for each t = tn (n = 1, 2, ..., N ) such that ß (∇γt unh , vh )Ω0 + A(unh , vh )Ω0 = (f n , vh )Ω0 , ∀vh ∈ Vh , (28) u0h = u0h , Now we will prove the stability of scheme (28). Theorem 5.4. The fully discrete variational scheme (28) is unconditionally stable.

(∇γt Ehn , vh )Ω0 + A(Ehn , vh )Ω0 = 0, i.e., (Ehn , vh )Ω0 + τ γ Γ(2 − γ)A(Ehn , vh )Ω0 = −( Then we have

i=1

(Ehn , vh )Ωi + τ γ Γ(2 − γ)

8 X

A(Ehn , vh )Ωi = −

n−1 X

bnk Ehk , vh )Ω0 .

k=0

8 n−1 X X ( bnk Ehk , vh )Ωi .

AN US

8 X

CR IP T

Proof. Assume that zhn (n = 1, 2, ..., N ) is another solution of the fully scheme (28), and let Ehn = unh − zhn , then

i=1

i=1 k=0

According to Theorem 4.3, we have ||Ehn ||L2 (Ωi ) ≤ ||Eh0 ||L2 (Ωi ) . Then we obtain Ã

8 X i=1

Ã

||Ehn ||2L2 (Ωi ) ≤

8 X i=1

||Eh0 ||2L2 (Ωi ) ,

M

i.e., ||Ehn ||L2 (Ω0 ) ≤ ||Eh0 ||L2 (Ω0 ) , which means the fully discrete scheme (28) is unconditional stable.

PT

6. Numerical examples

ED

Following a similar process of proof, we can obtain that Theorem 4.4 also hold for multi-connected domain Ω0 , which means that we can extend the analysis of stability and convergence from convex domain Ωi to multi-connected domain Ω0 . The extension of our theory to multi-connected domain Ω0 is studied numerically in Example 4. We remark that the results from this initial investigation appear very promising and we believe the theory can be generalised to more complicated domain in a straightforward manner.

AC

CE

In this section, we present some numerical examples to verify the effectiveness of our theoretical analysis presented in section 4. We adopt the linear polynomials on triangles, where h is the maximum length of the triangles and Ne is the number of triangles in Th . By Theorem 4.4, it is expected that ||u(tn ) − unh ||0 ∼ O(h2 ), |||u(tn ) − unh |||(α,β) ∼ O(h2−σ ), σ = max(α, β) in the spatial direction and ||u(tn ) − unh ||0 ∼ O(τ 2−γ ) in the temporal direction. Here, we use the following formulation to calculate the convergence order: ( log(||E(h )|| /||E(h )|| ) 1 0 2 0 , in space, log(h1 /h2 ) Order = log(||E(τ1 )||0 /||E(τ2 )||0 ) , in time. log(τ1 /τ2 ) Example 1

Firstly, we consider the following 2D TSRFDE on an elliptical domain  2β ∂ 2α u(x,y,t) u(x,y,t) γ  + K2 ∂ ∂|y| + f (x, y, t), (x, y, t) ∈ Ω × (0, T ],  C 2β 0 Dt u(x, y, t) = K1 ∂|x|2α 1 2 2 2 u(x, y, 0) = (4x + y − 1) , (x, y) ∈ Ω, 10   u(x, y, t) = 0, (x, y, t) ∈ ∂Ω × (0, T ],

18

ACCEPTED MANUSCRIPT

where Ω = {(x, y)| 4x2 + y 2 < 1}, K1 = 1, K2 = 1, T = 1, f (x, y, t) =

t2−γ t2 + 1 n (4x2 + y 2 − 1)2 + 16[f1 (x, a0 , 2α) + g1 (x, b0 , 2α)] 5Γ(3 − γ) 20 cos(απ)

CR IP T

o + 8(y 2 − 1)[f2 (x, a0 , 2α) + g2 (x, b0 , 2α)] + (y 4 − 2y 2 + 1)[f3 (x, a0 , 2α) + g3 (x, b0 , 2α)] t2 + 1 n [f1 (y, c0 , 2β) + g1 (y, d0 , 2β)] + 8(x2 − 1)[f2 (y, c0 , 2β) + g2 (y, d0 , 2β)] + 20 cos(βπ) o + (16x4 − 8x2 + 1)[f3 (y, c0 , 2β) + g3 (y, d0 , 2β)] , p p 1p 1p 1 − y 2 , b0 = 1 − y 2 , c0 = − 1 − 4x2 , d0 = 1 − 4x2 , a0 = − 2 2 f1 (x, a, α) = a Dxα (x4 ), f2 (x, a, α) = a Dxα (x2 ), f3 (x, a, α) = a Dxα (1), g1 (x, b, α) = x Dbα (x4 ), g2 (x, b, α) = x Dbα (x2 ), g3 (x, b, α) = x Dbα (1). t2 +1 2 10 (4x

+ y 2 − 1)2 . 1

1

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.2

0.2

0

0

-0.2

-0.2

-0.4

-0.4

-0.6

-0.6

-0.8

-0.8

-1 -0.5

0

0.5

x

-1 -0.5

AN US

1

0.8

0.4

0.2

y

y

y

The exact solution is u(x, y, t) =

0

-0.2

-0.4

-0.6

-0.8

0

0.5

-1 -0.5

x

0

0.5

x

M

Figure 6: The unstructured triangular meshes used in the calculation for h = 0.12558, 0.08391 and 0.04531.

PT

ED

Figure 6 shows the elliptical domain partitioned by different unstructured triangular meshes. The corresponding numerical results are given in Table 1, which illustrates the L(α,β) error, L2 error and corresponding convergence order of h. Table 2 displays the L2 error and the convergence order of τ . From these two tables we can see that the expected convergence orders O(h2−σ ), O(h2 ), and O(τ 2−γ ) are attained. Table 3 shows the L2 error and the convergence order of τ = h for the second order temporal numerical scheme. We can see that the numerical results are in excellent agreement with the exact solution and we attain the second order, which demonstrates the effectiveness of the numerical method. Table 1: The L(α,β) error, L2 error and convergence order of h for different α, β at t = 1 with γ = 0.7, τ = 1/1000.

CE

γ = 0.7

AC

α = 0.75 β = 0.95

Example 2 domain Ω

α = 0.8 β = 0.8

Ne 70 468 1142 4324 70 468 1142 4324

h 3.0312E-01 1.2558E-01 8.3913E-02 4.5308E-02 3.0312E-01 1.2558E-01 8.3913E-02 4.5308E-02

L(α,β) error 8.9507E-02 3.5831E-02 2.1668E-02 1.0303E-02 9.0145E-02 3.1969E-02 1.8311E-02 8.1488E-03

Order – 1.04 1.25 1.21 – 1.18 1.38 1.31

L2 error 9.1296E-03 1.7483E-03 6.9743E-04 1.8504E-04 8.6868E-03 1.4612E-03 5.5341E-04 1.3987E-04

Order – 1.88 2.28 2.15 – 2.02 2.41 2.23

Next, we consider the following 2D TSRFDE without a source term on a general irregular convex  2β ∂ 2α u(x,y,t) u(x,y,t) γ  + ∂ ∂|y| , (x, y, t) ∈ Ω × (0, T ],  C 2β 0 Dt u(x, y, t) = ∂|x|2α π  u(x, y, 0) = 10 cos( 2 xy), (x, y) ∈ Ω,  u(x, y, t) = 0, (x, y, t) ∈ ∂Ω × (0, T ]. 19

ACCEPTED MANUSCRIPT

Table 2: The L2 error and convergence order of τ for γ = 0.7 at t = 1 with α = β = 0.8 and h2 ≈ τ 2−γ .

Ne 276 1142 1738

τ

h 1.8428E-01 8.3913E-02 6.9134E-02

1 14 1 46 1 61

L2 error 2.2122E-03 5.5444E-04 3.7308E-04

Order – 1.16 1.40

Table 3: The L2 error and convergence order of τ = h for the second order numerical scheme with γ = 0.7, α = β = 0.8 at t = 1.

h 1.8428E-01 8.3913E-02 6.9134E-02 4.5308E-02

L2 error 2.2182E-03 5.7249E-04 3.9132E-04 1.6187E-04

Order – 1.72 1.96 2.09

CR IP T

Ne 276 1142 1738 4324

2

1.6

1.4

y

1

0.8

0.6

0.4

0.2

1.8

1.6

1.6

1.4

1.4

1.2

1.2

1

ED

y

1.2

1.8

y

1.8

0

2

M

2

AN US

Figure 7 shows the different triangular unstructured meshes used to partition the domain Ω. Figure 8 shows a diffusion behaviour of u(x, y, t) at different time t = 0.2, 0.4, 0.8 that decays with increasing time. As we cannot determine the exact solution for this problem, we use an approximate solution unh derived by choosing a very fine mesh (h = 0.0557). To observe the convergence behaviour, we choose a set of points in the domain and derive the L2 error for different h. The corresponding numerical results are given in Table 4. Again, we see that the numerical results exhibits a convergence order that attains the expected value of O(h2 ), which means that the numerical method is effective on general irregular convex domains. As other arbitrarily shaped convex domains can be partitioned similarly, we can conclude that the method is applicable to other arbitrarily shaped convex domains as well.

0.8

0.6

0.6

0.4

0.4

0.2

0

0.2

0.4

0.6

0.8

1

1.4

1.6

1.8

2

0

0

0.2

0.2

0.4

0.6

0.8

1

x

1.2

1.4

1.6

1.8

2

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x

PT

x

1.2

1

0.8

CE

Figure 7: The unstructured triangular meshes used in the calculation for h = 0.2825, 0.1454 and 0.0557.

AC

Example 3 Then, we consider the following two dimensional coupled fractional Bloch-Torrey diffusion equation with a time-varying magnetic field gradient on human brain-like domain shown in Figure 2, in which we choose 41 sample points on the boundary of the human brain (Figure 2) and connect the adjacent points by a line to form an approximate boundary of the human brain (Figure 2(b)).  ã Å 2β  ∂ 2β Mx (x,y,t) γ−1C γ 2(β−1) ∂ Mx (x,y,t)  ω D M (x, y, t) = Dµ + + λ1 (t)My (x, y, t), (x, y, t) ∈ Ω × (0, T ],  x t 0  ∂|x|2β ∂|y|2β   Å ã  2β ∂ 2β My (x,y,t) γ 2(β−1) ∂ My (x,y,t) ω γ−1C + − λ1 (t)Mx (x, y, t), (x, y, t) ∈ Ω × (0, T ], (29) 2β 0 Dt My (x, y, t) = Dµ ∂|x| ∂|y|2β     Mx (x, y, 0) = 0, My (x, y, 0) = 100, x, y ∈ Ω,    Mx (x, y, t)|∂Ω = 0, My (x, y, t)|∂Ω = 0, (x, y, t) ∈ ∂Ω × (0, T ].

Here, we choose ω = 2, D = 1 × 10−3 , µ = p 15, λ1 (t) = t, T = 50, τ = 1/20 to observe the behaviour of the transverse magnetization |Mxy (x, y, t)| = Mx (x, y, t)2 + My (x, y, t)2 . The related numerical scheme for the system (29) is given in the appendix. Figure 9 shows the different triangular unstructured meshes used to partition the domain Ω. Figure 10 displays the solution behaviour for Mx (x, y, t), My (x, y, t) at the randomly chosen point 20

ACCEPTED MANUSCRIPT

t = 0.2

t = 0.4

t = 0.8

2.5

2.5

2

2

1.5

1.5

1 0.5 0

u(x, y, t)

2 1.5

u(x, y, t)

u(x, y, t)

2.5

1 0.5 0

1 0.5 0

-0.5

-0.5

-0.5

-1

-1

-1

-1.5 0

-1.5 0

2

0.5

x

2

0.5

2

0.5

2

0.5

1

1.5

x

y

0

1.5

1

1

1.5

x

y

0

1.5

1

1

1.5

-1.5 0

2

0.5

1.5

1

0.5

2

y

0

CR IP T

Figure 8: The diffusion profiles of u(x, y, t) at different t with γ = 0.9, α = β = 0.80, h = 0.1454, τ = 1/1000 at t = 1.0.

Table 4: The L2 error and convergence order of h for α = β = 0.8, γ = 0.9, τ = 1/1000 at t = 1.0.

Ne 244 958 2476

h 2.8250E-01 1.4536E-01 8.9184E-02

L2 error 5.0708E-03 1.6899E-03 5.4918E-04

Order – 1.65 2.30

ED

M

AN US

(x∗ , y ∗ ) = (0.5702, 0.8548) with different values of γ. We observe that the effects of γ on the solution behaviour is significant and the smaller γ is, the faster it decays from (0, 100) to (0, 0). Figure 11 highlights a new finding of the effect of β on the solution behavior, which we believe is an original contribution to the literature. Although not obvious, the effects of β on the solution behaviour is similar to that of γ, which can be reflected clearly in Figure 12(b). Figure 12 exhibits the normalized decay of the transverse magnetization versus t at the point (x∗ , y ∗ ) = (0.5702, 0.8548) for different γ and β. We conclude that the smaller γ is, the sharper the decay rate is at the beginning time, but more time consuming in the process of decaying to 0 at latter times, which conforms with Figure10. The effects of β is similar, which are aligned with Figure 11. Finally, we choose a very fine grid mesh (h = 0.03278) to observe the discrete L2 error of Mxy for different h. The corresponding numerical results are given in Table 5 exhibiting the convergence order of O(h2 ), which shows the effectiveness of the method on the coupled fractional Bloch-Torrey equation.

1.2

1

β = 0.60

100

β = 1.00 100

0.8

80

40

My (x∗ , y∗ , t)

PT

0.6

40

My (x∗ , y∗ , t)

60

y

60

80

20

0.4

0

-20 -40

0.2

-60

CE

-80

-80

-60

-80 0

-40

-20

0

20

Mx (x∗ , y∗ , t)

40

60

0 -20 -40

-60

-100 -100

20

80

100

0

0.1

0.2

0.3

0.4

0.5

x

0.6

0.7

0.8

0.9

1

-100 -100

-80

-60

-40

-20

0

20

40

60

80

100

Mx (x∗ , y∗ , t)

AC

Figure 9: The unstructured triangular meshes used in the calculation for h = 0.11624, 0.06943 and 0.03278.

Example 4 To further demonstrate the flexibility and effectiveness of our method, we consider the following 2D TSRFDE without a source term on a multiply-connected domain (see Figure 13)  2β ∂ 2α u(x,y,t) u(x,y,t) γ  , (x, y, t) ∈ Ω × (0, T ], + ∂ ∂|y|  C 2β 0 Dt u(x, y, t) = ∂|x|2α π u(x, y, 0) = 10 cos( 2 xy), (x, y) ∈ Ω,   u(x, y, t) = 0, (x, y, t) ∈ ∂Ω × (0, T ].

Figure 14 illustrates the behaviour of u(x, y, t) at different times t = 0.2, 0.4, 0.8 on a multiply-connected domain, which decays with increasing time. We also choose a very fine grid mesh (h = 0.0693) to observe the L2 error of u(x, y, t) for different h. The corresponding numerical results are given in Table 6 exhibiting the convergence order of O(h2 ), which is in agreement with the analysis presented in Section 5. 21

ACCEPTED MANUSCRIPT

γ = 0.80 1

1

γ=0.99 γ=0.90 γ=0.80

0.9 0.8

0.8

80 60

0.7

0.6 0.5 0.4

0.6 0.5 0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0

5

10

15

20

25

30

35

40

45

50

40

My (x∗ , y∗ , t)

|Mxy /M0 |

0.7

|Mxy /M0 |

100

β=1.0 β=0.8 β=0.6

0.9

20 0 -20 -40 -60 -80

0

5

10

15

20

t

25

30

35

40

45

-100 -100

50

-80

-60

-40

t

-20

0

20

40

60

80

100

Mx (x∗ , y∗ , t)

β = 0.80

β = 1.00

100

100

80

80

80

60

60

60

40

40

40

20 0 -20 -40

20 0 -20 -40

-60

-80

-60

-40

-20

0

20

40

60

80

100

Mx (x , y , t) ∗

-80

-80

-60

-40

-20

0

20

Mx (x , y , t)

∗

0 -20

-60

-80 -100 -100

20

-40

-60

-80 -100 -100

My (x∗ , y∗ , t)

100

My (x∗ , y∗ , t)

My (x∗ , y∗ , t)

β = 0.60

∗

∗

CR IP T

Figure 10: Plots of Mx (x, y, t), My (x, y, t) at point (x∗ , y ∗ ) = (0.5702, 0.8548) for different γ with β = 1.0, h = 0.06943.

40

60

80

100

-100 -100

-80

-60

-40

-20

0

20

40

60

80

100

Mx (x∗ , y∗ , t)

AN US

Figure 11: Plots of Mx (x, y, t), My (x, y, t) at point (x∗ , y ∗ ) = (0.5702, 0.8548) for different β with γ = 0.99, h = 0.06943.

7. Conclusions

Acknowledgment

PT

ED

M

In this paper, we considered the Galerkin FEM to a class of two-dimensional time-space Riesz fractional diffusion equation on irregular convex domains. We partitioned the irregular convex domain into a sum of unstructured triangular meshes. Then utilising FEM, we obtained the variation formulation of the problem and the associated discrete scheme with the accuracy of O(τ 2−γ + h2 ). Furthermore, we reduced the computation of inner products from the whole domain Ω to a single triangular element and evaluated it approximately by the Gauss quadrature technique. Moreover, we derived a second order temporal numerical scheme for the problem. Finally, numerical examples on irregular convex domains were conducted, which verified the effectiveness and reliability of the method. Furthermore, with our numerical method, we are able to exhibit the effects of the time and space fractional indices for the coupled two-dimensional fractional Bloch-Torrey equation. We concluded that the numerical method can be extended to other arbitrarily shaped convex domains and even some non-convex domains. In future work, we shall investigate the FEM to other fractional problems on irregular convex domains, such as the two-dimensional FDE with variable coefficients.

AC

CE

Authors Turner and Liu wish to acknowledge that this research was partially supported by the Australian Research Council (ARC) via the Discovery Project (DP180103858), authors Liu wishes to acknowledge that this research was partially supported by Natural Science Foundation of China (Grant No.11772046), and author Zhuang wishes to acknowledge that this research was partially supported by Natural Science Foundation of China (Grant No. 11771364). All the authors wish to thank the editor and anonymous referees for their many constructive comments and suggestions that resulted in an improved final version of the paper. Appendix

In this section, we outline the numerical scheme for solving the coupled fractional Bloch-Torrey equation (29). Assume that Xhn ∈ Vh and Yhn ∈ Vh are the approximations of Mx (x, y, t) and My (x, y, t) at t = tn , respectively. Then the fully discrete scheme associated with the variational form of (29) is: find Xhn ∈ Vh and Yhn ∈ Vh for each t = tn (n = 1, 2, ..., N ) such that ß (∇γt Xhn , vh )Ω + A(Xhn , vh )Ω = λ2 (tn )(Y n , vh )Ω , ∀vh ∈ Vh , (30) (∇γt Yhn , vh )Ω + A(Yhn , vh )Ω = −λ2 (tn )(X n , vh )Ω , ∀vh ∈ Vh , 22

ACCEPTED MANUSCRIPT

t = 0.4

0.7

0.7

0.6

0.6

u(x, y, t)

u(x, y, t)

t = 0.2

0.5 0.4 0.3

0.5 0.4 0.3

0.2

0.2

0.1

0.1

0 1

0 1 0.5

0.5

1

y

1

0.5

0 -0.5

-0.5 -1

-1

0.5

0

0

y

x

0

-0.5

-0.5 -1

-1

x

CR IP T

Figure 12: Normalized decay of the transverse magnetization versus t at point (x∗ , y ∗ ) = (0.5702, 0.8548) for different γ (with fixed β = 1.0) and β (with fixed γ = 0.99).

Table 5: The L2 error and convergence order of Mxy for different h with β = 0.8, γ = 0.95, τ = 1/20 at t = 50.

h 1.3485E-01 6.9427E-01 5.1856E-02

L2 error 4.0010E-04 9.9154E-05 5.5050E-05

Order – 2.10 2.02

AN US

Ne 221 882 1777

PNp n PNp n (tn ) where λ2 (tn ) = λω1γ−1 . Substituting Xhn = i=1 Xi li (x, y) and Yhn = i=1 Yi li (x, y) into (30) with vh = lj (x, y) leads to  N i Np Np n−1 Pp n h P P P n k   Xi (li , lj )Ω + ω0 A(li , lj )Ω − ω0 λ2 (tn ) Yin (li , lj )Ω = − bk Xi (lk , lj )Ω − bn0 (X 0 , lj )Ω ,  i=1 i=1 i=1 k=1 (31) h i Np Np Np n−1 P P P P n k    Yin (li , lj )Ω + ω0 A(li , lj )Ω + ω0 λ2 (tn ) Xin (li , lj )Ω = − bk Yi (lk , lj )Ω − bn0 (Y 0 , lj )Ω . i=1

i=1 k=1

M

i=1

ED

where ω0 = τ γ Γ(2 − γ). Equation (31) can be written in matrix form as  n−1 P n k   bk X − G0X ,  (M + ω0 B)X n − ω0 λ2 (tn )M Y n = −M k=1

n−1 P n k   bk Y − G0Y .  ω0 λ2 (tn )M X n + (M + ω0 B)Y n = −M

(32)

k=1

AC

CE

PT

where M is the mass matrix with elements Mij = (lj , li )Ω , B is the stiffness matrix with elements Bij = n T A(lj , li )Ω , X n = [X1n , X2n , ..., XN ] , Y n = [Y1n , Y2n , ..., YNnp ]T . The vectors G0X and G0Y are given by G0X = p bn0 [(X 0 , l1 )Ω , (X 0 , l2 )Ω , ..., (X 0 , lNp )Ω ]T and G0Y = bn0 [(Y 0 , l1 )Ω , (Y 0 , l2 )Ω , ..., (Y 0 , lNp )Ω ]T respectively. Finally, equation (32) can be recast into the form Ü ê n−1 P n k Å ãÅ n ã bk X − G0X −M M + ω0 B −ω0 λ2 (tn )M X k=1 = . n−1 P n k ω0 λ2 (tn )M M + ω0 B Yn −M bk Y − G0Y k=1

This linear system is then solved using general iterative methods. References

[1] I.M. Sokolov, J. Klafter, A. Blumen, Fractional kinetics, Phys. Today 55 (2002) 48-54. [2] G.M. Zaslavsky, Chaos, fractional kinetics, and anomalous transport, Phys. Rep. 371 (2002) 461-580. [3] R.L. Magin, Fractional calculus in bioengineering, Redding: Begell House, 2006.

23

ACCEPTED MANUSCRIPT

1

y

0.5

0

-0.5

-1 -1

-0.5

0

0.5

1

CR IP T

x

Figure 13: The illustration of a multiply-connected domain with triangular partition.

γ = 0.90

γ = 0.80

100

80

80

60

60

60

40

40

40

20 0 -20

My (x∗ , y∗ , t)

100

80

My (x∗ , y∗ , t)

20 0 -20

20 0 -20

-40

-40

-40

-60

-60

-60

-80

-80

-100 -80

-60

-40

-20

0

20

40

60

80

100

-100 -100

Mx (x∗ , y∗ , t)

AN US

My (x∗ , y∗ , t)

γ = 0.99 100

-80

-80

-60

-40

-20

0

20

Mx (x∗ , y∗ , t)

40

60

80

100

-100 -100

-80

-60

-40

-20

0

20

40

60

80

100

Mx (x∗ , y∗ , t)

Figure 14: The diffusion profiles of u(x, y, t) on a multiply-connected domain at different t with γ = 0.9, α = β = 0.80, τ = 1/1000 at t = 1.0.

M

[4] L. Liu, L. Zheng, F. Liu, X. Zhang, Anomalous convection diffusion and wave coupling transport of cells on comb frame with fractional Cattaneo-Christov flux, Commun. Nonlinear Sci. Numer. Simulat. 38 (2016) 45-58.

ED

[5] S.B. Yuste, L. Acedo, K. Lindenberg, Reaction front in an A + B → C reaction-subdiffusion process, Phys. Rev. E 69 (2004) 036126. [6] D.A. Benson, S.W. Wheatcraft, M.M. Meerschaert, Application of a fractional advection-dispersion equation, Water Resouc. Res. (2000) 1403-1412.

PT

[7] F. Liu, V. Anh, I. Turner, Numerical solution of the space fractional Fokker-Planck equation, J. Comput. Appl. Math. 166 (2004) 209-219.

CE

[8] E. Scalas, R. Gorenflo and F. Mainardi, Fractional calculus and continuous-time finance, Phys. A 284 (2000) 376-384. [9] I. Podlubny, Fractional differential equations, Academic Press, San Diego (1999).

AC

[10] A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and applications of fractional differential equations, Elsevier (North-Holland) Science Publishers, Amsterdam, 2006. [11] K. Diethelm, The analysis of fractional differential equations: an application-oriented exposition using differential opereators of Caputo type, Springer-Verlag, Berlin 2010. [12] J.P. Roop, Variational solution of the fractional advection-dispersion equation, PhD thesis, Clemson University, 2004. [13] V.J. Ervin, J.P. Roop, Variational formulation for the stationary fractional advection dispersion equation, Numer. Methods Partial Differential Equation 22 (2006) 558-576. [14] W. Deng, Finite element method for the space and time fractional Fokker-Planck equation, SIAM J. Numer. Anal. 47 (2008) 204-226. 24

ACCEPTED MANUSCRIPT

Table 6: The L2 error and convergence order of h for α = β = 0.8, γ = 0.9, τ = 1/1000 at t = 1.0 on the multiply-connected domain.

Ne 418 714 1710

h 1.7470E-01 1.3974E-01 8.7473E-02

L2 error 3.4886E-03 2.0598E-03 7.1365E-04

Order – 2.36 2.26

CR IP T

[15] H. Zhang, F. Liu, V. Anh, Galerkin finite element approximation of symmetric space-fractional partial differential equations, Appl. Math. Comput. 217 (2010) 2534-2545. [16] C. Li, Z. Zhao, Y.Q. Chen, Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion, Comput. Math. Appl. 62 (2011) 855-875. [17] Q. Liu, F. Liu, I. Turner, V. Anh, Finite element approximation for a modified anomalous subdiffusion equation, Appl. Math. Model. 35 (2011) 4103-4116.

AN US

[18] F. Zeng , C. Li, F. Liu, I. Turner, The use of finite difference/element approaches for solving the timefractional subdiffusion equation, SIAM J. Sci. Comput. 35 (2013) A2976-A3000. [19] B. Jin, R. Lazarov, Z. Zhou, Error estimates for a semidiscrete finite element method for fractional order parabolic equations, SIAM Journal on Numerical Analysis 51 (2013) 445-466. [20] B. Jin, R. Lazarov, Y. Liu, Z. Zhou, The Galerkin finite element method for a multi-term time-fractional diffusion equation, Journal of Computational Physics 281 (2015) 825-843. [21] W. Bu, Y. Tang, J. Yang, Galerkin finite element method for two-dimensional Riesz space fractional diffusion equations, J. Comput. Phys. 276 (2014) 26-38.

M

[22] Y. Liu, Z. Fang, H. Li, S. He, A mixed finite element method for a time-fractional fourth-order partial differential equation, Appl. Math. Comput. 243 (2014) 703-717.

ED

[23] L.B. Feng, P. Zhuang, F. Liu, I. Turner, Y.T. Gu, Finite element method for space-time fractional diffusion equation, Numer. Algorithms 72 (2016) 749-767. [24] 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-66.

PT

[25] Y.M. Zhao, Y.D. Zhang, F. Liu, I. Turner, D.Y. Shi, Analytical solution and nonconforming finite element approximation for the 2D multi-term fractional subdiffusion equation, Appl. Math. Model. 40 (2016) 88108825.

CE

[26] B. Jin, R. Lazarov, Z. Zhou, A Petrov-Galerkin finite element method for fractional convection-diffusion equations, SIAM J. Numer. Anal. 54 (2016) 481-503.

AC

[27] Z. Yang, Z. Yuan, Y. Nie, J. Wang, X. Zhu, F. Liu, Finite element method for nonlinear Riesz space fractional diffusion equations on irregular domains, J. Comput. Phys. 330 (2017) 863-883. [28] W. Fan, F. Liu, X. Jiang, I. Turner, A novel unstructured mesh finite element method for solving the time-space fractional wave equation on a two-dimensional irregular convex domain, Fractional Calculus and Applied Analysis 20 (2017) 352-383. [29] L.L. Qiu, W.H. Deng, J.S. Hesthaven, Nodal discontinuous Galerkin methods for fractional diffusion equations on 2D domain with triangular meshes, J. Comput. Phys. 298 (2015) 678-694. [30] Q. Yang, I. Turner, T. Moroney, F. Liu, A finite volume scheme with preconditioned Lanczos method for two-dimensional space-fractional reaction-diffusion equations, Appl. Math. Model. 38 (2014) 3755-3762.

25

ACCEPTED MANUSCRIPT

[31] F. Liu, P. Zhuang, I. Turner, V. Anh, K. Burrage, A semi-alternating direction method for a 2-D fractional FitzHugh-Nagumo monodomain model on an approximate irregular domain, J. Comput. Phys. 293 (2015) 252-263. [32] Z. Cai, On the finite volume element method, Numer. Math. 58 (1991) 713-735. [33] E. S¨ uli, Convergence of finite volume schemes for Poisson equation on nonuniform meshes, SIAM J. Numer. Anal. 28 (1991) 1419-1430.

CR IP T

[34] P. Chatzipantelidis, R.D. Lazarov, V. Thom´ee, Parabolic finite volume element equations in nonconvex polygonal domains, Numer. Methods Partial Differential Equation 25 (2009) 507–525. [35] S. Karaa, K. Mustapha, A.K. Pani, Finite volume element method for two-dimensional fractional subdiffusion problems, IMA Journal of Numerical Analysis 37 (2016) 945-964. [36] K.N. Le, W. McLean, B. Lamichhane, Finite element approximation of a time-fractional diffusion problem for a domain with a re-entrant corner, The ANZIAM Journal 59 (2017) 61-82 [37] W.P. Bu, Y.F. Tang, Y.C. Wu, J.Y. Yang, Finite difference/finite element method for two-dimensional space and time fractional Bloch-Torrey equations, J. Comput. Phys. 293 (2015) 264-279.

AN US

[38] C. Geuzaine, J.F. Remacle, Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities, International Journal for Numerical Methods in Engineering 79 (2009) 1309-1331. [39] K.M. Bennett, J.S. Hyde and K.M. Schmainda, Water diffusion heterogeneity index in the human brain is insensitive to the orientation of applied magnetic field gradients, Magn. Reson. Med. 56 (2006) 235-239. [40] R.L. Magin, O. Abdullah, D. Baleanu, and X.J. Zhou, Anomalous diffusion expressed through fractional order differential operators in the Bloch-Torrey equation, Journal of Magnetic Resonanc 190 (2008) 255-270.

M

[41] Q. Yu, F. Liu, I. Turner and K. Burrage, Stability and convergence of an implicit numerical method for the space and time fractional Bloch-Torrey equation, Phil. Trans. R. Soc. A 371 (2013) 20120150.

ED

[42] S. Qin, F. Liu, I. Turner, Q. Yu, Q. Yang, V. Vegh, Characterization of anomalous relaxation using the time-fractional Bloch equation and multiple echo T2∗ -weighted magnetic resonance imaging at 7 T, Magn. Reson. Med. 77 (2017) 1485-1494.

PT

[43] S. Qin, F. Liu, I. Turner, Q. Yang, Q. Yu, Modelling anomalous diffusion using fractional BlochTorrey equations on approximate irregular domains, Computers and Mathematics with Applications (2017) http://dx.doi.org/10.1016/j.camwa.2017.08.032 [44] R.A. Adams, Sobolev spaces, Academic Press, New York, 1975.

CE

[45] Z. Sun, X. Wu, A fully discrete difference scheme for a diffusion-wave system, Appl. Numer. Math. 56 (2006) 193-209.

AC

[46] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for fractional advection-dispersion flow equations, J. Comput. Appl. Math. 172 (2004) 65-77. [47] Z. Wang, S. Vong, Compact difference schemes for the modified anomalous fractional sub-diffusion equation and the fractional diffusion-wave equation, J. Comput. Phys. 277 (2014) 1-15. [48] W. Tian, H. Zhou, W. Deng, A class of second order difference approximations for solving space fractional diffusion equations, Math. Comput. 84 (2015) 1703-1727. [49] F. Zeng, Z. Zhang, G.E. Karniadakis, Second-order numerical methods for multi-term fractional differential equations: smooth and non-smooth solutions, Computer Methods in Applied Mechanics and Engineering 327 (2017) 478-502 [50] S.C. Brenner, L.R. Scott, The mathematical theory of finite element methods, Springer-Verlag, New York, 1994. 26

ACCEPTED MANUSCRIPT

[51] H. Wang, D. Yang, Wellposedness of variable-coefficient conservative fractional elliptic differential equations, SIAM J. Numer. Anal. 51 (2013) 1088-1107. [52] B. Jin, R. Lazarov, J. Pasciak, W. Rundell, Variational formulation of problems involving fractional order differential operators, Math. Comp. 84 (2015), 2665-2700.

AC

CE

PT

ED

M

AN US

CR IP T

[53] V. J. Ervin, N. Heuer, J. P. Roop, Regularity of the solution to 1-D fractional order diffusion equations, Math. Comp. (2017) https://doi.org/10.1090/mcom/3295.

27