# A Finite Element Method for the Multiterm Time-Space Riesz Fractional Advection-Diffusion Equations in Finite Domain

Authors: Jingjun Zhao; Jingyu Xiao; Yang Xu
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We present an effective finite element method (FEM) for the multiterm time-space Riesz fractional advection-diffusion equations (MT-TS-RFADEs). We obtain the weak formulation of MT-TS-RFADEs and prove the existence and uniqueness of weak solution by the Lax-Milgram theorem. For multiterm time discretization, we use the Diethelm fractional backward finite difference method based on quadrature. For spatial
discretization, we show the details of an FEM for such MT-TS-RFADEs. Then, stability and convergence of such numerical method are proved, and some numerical examples are given to match well with the main conclusions.


---

## Body

## 1. Introduction

Fractional differential equations are different from integer ones, in which the nature of the fractional derivative introduces the memory effect, thus increasing its modeling ability. Recently, many mathematical models with fractional derivatives have been successfully applied in biology, physics, chemistry, and biochemistry, hydrology, and finance [1–3]. The multiterm fractional differential equations have been widely studied in rheology, and, in many cases, the exact solutions are known [4, 5]. Summary of the fractional differential equations can be found in monographs [6–9]. As one of the main branch, fractional partial differential equations have attracted great attention. Therefore, the numerical treatment and supporting analysis of fractional order partial differential equations have become an important research topic that offers great potential.The FEM is one of the effective numerical methods for solving traditional partial differential equations. For fractional partial differential equations, FEM also can be a useful and effective numerical method. In recent years, some valuable papers are concerned with the FEM for fractional differential equations. Adolfsson et al. [10, 11] considered an efficient numerical method to integrate the constitutive response of fractional order, viscoelasticity based on the FEM. Roop and Ervin [12–15] investigated the theoretical framework for the Galerkin finite element approximation to some kinds of fractional partial differential equations. Li et al. [16] considered numerical approximation of fractional differential equations with subdiffusion and superdiffusion by using difference method and finite element method. Li and Xu [17, 18] proposed a time-space spectral method for time and time-space fractional partial differential equation based on a weak formulation, and a detailed error analysis was carried out. Jiang and Ma [19] considered a high-order FEM for time fractional partial differential equations and proved the optimal order error estimates. Ford et al. [20] studied an FEM for time fractional partial differential equations.Fractional advection-diffusion equations especially are important in describing and understanding the dispersion phenomena. Analytical solutions of such equations in finite domain have been obtained by Park in [21]. Also, the Riesz fractional advection-diffusion equations (RFADEs) with a symmetric fractional derivative (the Riesz fractional derivative) were derived from the kinetics of chaotic dynamics by [22] and summarized by [23]. Ciesielski and Leszczynski [24] presented a numerical solution for such equations based on the finite difference method.One often sees RFADEs defined in terms of the fractional Laplacian as follows:(1)ddtu(t,x)=-kα(-Δ)α/2u(t,x)-kβ(-Δ)β/2u(t,x),
for example, where u is a solute concentration, kα and kβ represent the dispersion coefficient and the average fluid velocity. Here, the fractional Laplacian operator -(-Δ)α/2 uses the Fourier transformation on an infinite domain, with a natural extension to include finite domains when the function u(t,x) is subject to the zero Dirichlet boundary conditions (see [9]). Due to Lemma  1 in [25], the fractional Laplacian operator -(-Δ)α/2 on an infinite domain x∈(-∞,∞) is equivalent to the Riesz fractional derivative operator RD|x|α. In particular, the Riesz fractional derivatives include both the left and the right Riemann-Liouville derivatives that allow the modeling of flow regime impacts from either side of the domain. Yang et al. [25] investigated the numerical treatment for the RFADE with the Riesz space fractional derivative as
(2)ddtu(t,x)=kαRD|x|αu(t,x)+kβRD|x|βu(t,x),
where RD|x|α is the Riesz space fractional operator defined in Section 2.To increase the modeling ability, some authors considered the equations with the fractional order in both time and spatial variables in RFADEs, which include more information and hence are more interesting. For the time-space fractional advection-dispersion equations, Shen et al. [26] presented the fundamental solution and numerical solution of the Riesz fractional advection-dispersion equation with initial and boundary conditions on a bounded domain, and derived the stability and convergence of their proposed numerical methods. Then, for fractional advection-diffusion equations, Shen et al. [27] presented an explicit difference approximation and an implicit difference approximation for the time-space Riesz-Caputo fractional advection-diffusion equations with initial and boundary conditions on a finite domain.All of the above papers only considered single-term fractional equations in time variable, where only one fractional differential operator appeared. In this paper, we consider the multiterm fractional differential equation, which includes more than one fractional derivative. For example, the so-called Bagley-Torvik equation [28] is
(3)ACD2y(t)+BCD3/2y(t)+Cy(t)=f(t),
where A,B, and C are certain constants and f is a given function. The Basset equation is, in [6],
(4)D1y(t)+bDαy(t)+cy(t)=f(t),y(0)=y0,
where 0<α<1, b, and c are positive real numbers. This equation describes the forces that occur when a spherical object sinks in an incompressible viscous fluid.Recently, some authors considered the applications of multiterm fractional differential equations [29] and the numerical methods for such equations [30–32]. At the same time, the multiterm fractional partial differential equations have been proposed in [33, 34]. The analytical solution and the numerical methods for multiterm time fractional wave-diffusion equations have been investigated in [35, 36]. This motivates us to consider the effective numerical solution for such multiterm fractional partial differential equations.In this paper, we consider MT-TS-RFADEs in finite domain with the zero Dirichlet boundary conditions. The analytical solution of such MT-TS-RFADEs has been investigated by Jiang et al. in [37]. Here, we present an FEM for a simplified MT-TS-RFADEs and obtain the optimal order error estimates both in semidiscrete and fully discrete cases and derive the stability of such FEM. As far as we are aware, there are few research papers in the published literature written on this topic.This paper is organized as follows. In Section2, the preliminaries of the fractional calculation are shown. Then, we give the weak formulation of MT-TS-RFADEs and prove the existence and uniqueness of this problems by the well-known Lax-Milgram theorem. In Section 3, we present the convergence rate of Diethelm’s fractional backward difference method (see [38, 39]) for time discretization. In Section 4, we propose a finite element method based on the weak formulations and carry out the error analysis. In Section 5, we prove the stability of such FEM for MT-TS-RFADEs. Finally, some numerical examples are considered in Section 6.


## 2. Existence and Uniqueness

We consider the MT-TS-RFADEs with multiterm time fractional derivatives and the Riesz space fractional derivatives in the following form:(5)P(0Dt)u(t,x)=kβDR|x|βu(t,x)+kγDR|x|γu(t,x)+f(t,x),
where 0<β<1, 1<γ<2, x∈[0,X], and t∈[0,T] are respectively the space and time variables and kβ, kγ are positive constants. We consider this problem with the zero Dirichlet boundary value conditions and the initial value condition defined as follows:
(6)u(t,0)=u(t,X)=0,t∈[0,T],u(0,x)=ϕ(x),x∈(0,X).
For nonzero boundary value conditions, we need to transform the problem into one with zero boundary value conditions before using the method in this paper.Here, we consider the multiterm time fractional differential operator which has the subdiffusion process (see [16]). It is different from (2), which only has one integer order differential operator in time.Note that the analytical solutions for MT-TS-RFADEs have been studied in [37], in which this problem is well defined. For this problem, some new techniques have been used, such as a spectral representation of the fractional Laplacian operator and the equivalent relationship between the fractional Laplacian operator and the Riesz fractional derivative.For convenience, we introduce the following definitions and properties. The space derivativesDR|x|βu(t,x) and DR|x|γu(t,x) are the Riesz space fractional derivatives of order β and γ, respectively. The definitions of them can be found in [40].LetΓ(·) denote the gamma function. For any positive integer n and real number θ (n-1<θ<n), there are different definitions of fractional derivatives with order θ in [8]. During this paper, we consider the left, (right) Caputo derivative and left (right) Riemann-Liouville derivative defined as follows:(i)
the left Caputo derivative:(7)0CDtθv(t)=1Γ(n-θ)∫0t‍1(t-τ)θ-n+1(dndτnv(τ))dτ,(ii)
the right Caputo derivative:(8)tCDTθv(t)=(-1)nΓ(n-θ)∫tT‍1(τ-t)θ-n+1(dndτnv(τ))dτ,(iii)
the left Riemann-Liouville derivative:(9)0RDtθv(t)=1Γ(n-θ)dndtn∫0t‍v(τ)(t-τ)θ-n+1dτ,(iv)
the right Riemann-Liouville derivative:(10)tRDTθv(t)=(-1)nΓ(n-θ)dndtn∫tT‍v(τ)(τ-t)θ-n+1dτ.The Riesz fractional operatorsDR|x|β and DR|x|γ in (5) can be defined by the left and the right Riemann-Liouville fractional derivatives.Definition 1 (see [40]).
The Riesz fractional derivatives of orderθ for n∈N,n-1<θ<n, on a finite interval 0≤x≤X is defined as
(11)DR|x|θv(t,x)=-Cθ(0RDxθv(t,x)+xRDXθv(t,x)),x∈[0,X],
where Cθ=1/(2cos(πθ/2)).The multiterm fractional operatorP(0Dt) in (5) is
(12)P(0Dt)u(t,x)=(0CDtα+∑i=1s‍ai0CDtαi)u(t,x),
where 0<αs<αs-1<⋯<α1<α<1 or 1<αs<αs-1<⋯<α1<α<2 or 0<αs<⋯<αs0<1<αs0-1<⋯<α1<α<2, and ai∈R, i=1,…,s,s∈N. Here, 0CDtα and 0CDtαi(i=1,…,s) denote the left Caputo fractional derivatives with respect to the time variable t of order 0<α<1 and 0<αs<⋯<α1<1. There are three cases for multiterm time-space fractional derivatives.Case 1.
If0<αs<⋯<α1<α<1, (5) is a generalized MT-TS-RFADE with multiterm time fractional diffusion terms with initial conditions given as (6). And especially, if αi=0, i=1,…,s and α=1, then (5) becomes a space fractional advection-diffusion equation with the Riesz space fractional derivatives, which was discussed by Yang et al. [25].Case 2.
If1<αs<⋯<α1<α<2, (5) is a generalized MT-TS-RFADE with multiterm time fractional wave terms. In this case, the initial conditions are given as follows:
(13)u(0,x)=ϕ(x),ut(0,x)=ψ(x),x∈(0,X).Case 3.
If0<αs<⋯<αs0<1<αs0-1<⋯<α1<α<2, (5) becomes a generalized MT-TS-RFADE, which we refer to as a multiterm time-space fractional mixed wave-diffusion equation. In this case, the initial conditions are also given by (13).In this paper, we just consider Case1 of (5) with 0<αs<⋯<α1<α<1. Another two cases with 1<αs<⋯<α1<α<2 and 0<αs<⋯<αs0<1<αs0-1<⋯<α1<α<2 will be studied in our following work.In order to establish the weak formulation of the problem (5), we need some preparatory work. We use definitions of functional spaces and derive some properties related to these spaces. Let C∞(0,T) denote the space of infinitely differentiable functions on (0,T), and let C0∞(0,T) denote the space of infinitely differentiable functions with compact support in (0,T). Let L2(𝒬) be the space of measurable functions whose square is the Lebesgue integrable in 𝒬, which may denote a domain 𝒬=I or Ω or I×Ω, where I=[0,T] denotes the time domain and Ω=[0,X] denotes the space domain. The inner product and norm of L2(𝒬) are defined by
(14)(u,v)L2(𝒬)=∫𝒬‍uvd𝒬,∥u∥L2(𝒬)=(u,u)L2(𝒬)1/2,∀u,v∈L2(𝒬).For any realσ>0, we define the spaces lH0σ(𝒬) and rH0σ(𝒬) to be the closure of C0∞(𝒬) with respect to the norms ∥v∥lH0σ(𝒬) and ∥v∥rH0σ(𝒬), respectively, where
(15)∥v∥lH0σ(𝒬)≔(∥v∥L2(𝒬)2+|v|lH0σ(𝒬)2)1/2,|v|lH0σ(𝒬)2:=∥0RDtσv∥L2(𝒬)2,∥v∥rH0σ(𝒬)≔(∥v∥L2(𝒬)2+|v|rH0σ(𝒬)2)1/2,|v|rH0σ(𝒬)2:=∥tRDTσv∥L2(𝒬)2.
In the usual Sobolev space H0σ(𝒬), we also have the definition
(16)∥v∥H0σ(𝒬)≔(∥v∥L2(𝒬)2+|v|H0σ(𝒬)2)1/2,|v|H0σ(𝒬)2:=((0RDtσv,tRDTσv)L2(𝒬)cos(πσ)).
From [18], for σ>0,σ≠n-1/2, the spaces lH0σ(𝒬), rH0σ(𝒬), and H0σ(𝒬) are equal, and their seminorms are all equivalent to |·|H0σ(𝒬).We now give some results for fractional operators on these spaces.Lemma 2 (see [8, 18]).
(
1
) For real 0<θ<1,0<δ<1, if v(0)=0,x∈(0,X), then
(17)0RDxθ+δv(x)=(0RDxθ)(0RDxδ)v(x)=(0RDxδ)(0RDxθ)v(x),∀v∈Hθ+δ(0,X).
(
2
)  Let 0<θ<1. Then, one has
(18)(0RDxθw,v)L2(0,X)=(w,xRDXθv)L2(0,X),∀w∈Hθ(0,X),v∈C0∞(0,X).Lemma 3 (see [18]).
Let0<θ<2,θ≠1. Then, for any w,v∈H0θ/2(0,X), then
(19)(0RDxθw,v)L2(0,X)=(0RDxθ/2w,tRDXθ/2v)L2(0,X).Following from the definitions and lemmas above, we define the space(20)Bα/2,γ/2(I×Ω):=Hα/2(I,L2(Ω))∩Hα1/2(I,L2(Ω))∩⋯∩Hαs/2(I,L2(Ω))∩L2(I,H0β/2(Ω))∩L2(I,H0γ/2(Ω)),
when 0<αs<⋯<α1<α<1,0<β<1 and 1<γ<2.Then, one obtains(21)Bα/2,γ/2(I×Ω)=Hα/2(I,L2(Ω))∩L2(I,H0γ/2(Ω)).
Here, Bα/2,γ/2(I×Ω) is a Banach space with respect to the following norm:
(22)∥v∥Bα/2,γ/2(I×Ω):=(∥v∥Hα/2(I,L2(Ω))2+∥v∥L2(I,H0γ/2(Ω))2)1/2,∀v∈Bα/2,γ/2(I×Ω),
where
(23)Hα/2(I,L2(Ω)):={v;∥v(t,·)∥L2(Ω)∈Hα/2(I)},
endowed with the norm
(24)∥v∥Hα/2(I,L2(Ω)):=∥∥v(t,·)∥L2(Ω)∥Hα/2(I).For obtaining a suitable weak solution for problem (5) with the Caputo time fractional derivation, we consider the connection between the Caputo and the Riemann-Liouville fractional definition. Based on the definitions of the Caputo and the Riemann-Liouville fractional differential operators, we have an immediate consequence, for any real order θ>0, in [38],
(25)0CDtθf(t)=1Γ(-θ)∫0t‍(t-s)-θ-1(f(s)-Tn-1[f;0](s))ds=0RDtθf(t)-0RDtθTn-1[f;0](t),
where Tn-1[f;0](t) denotes the Taylor polynomial for f of order n-1, centered at 0,
(26)Tn-1[f;0](t)=∑k=0n‍tkk!dkdtkf(0),
where f(t)∈Cn-1[0,T].Based on (25), we make
(27)f-(t,x)=f(t,x)+u0(x)t-αΓ(1-α)+∑i=1s‍aiu0(x)t-αiΓ(1-αi).
Therefore, that led to the following weak formulation of (5). Let Bα/2,γ/2(I×Ω)′ be the dual space of Bα/2,γ/2(I×Ω). For f-∈Bα/2,γ/2(I×Ω)′, we find u(t,x)∈Bα/2,γ/2(I×Ω) such that
(28)𝒜(u,v)=ℱ(v),v∈Bα/2,γ/2(I×Ω),
where the bilinear form 𝒜(·,·) is, based on Lemmas 2 and 3,
(29)𝒜(u,v):=(0RDtα/2u,tRDTα/2v)L2(𝒬)+∑i=1s‍ai(0RDtαi/2u,tRDTαi/2v)L2(𝒬)+kβCβ[(0RDxβ/2u,xRDXβ/2v)L2(𝒬)+kβCβ+(xRDXβ/2u,0RDxβ/2v)L2(𝒬)]+kγCγ[(0RDxγ/2u,xRDXγ/2v)L2(𝒬)+kγCγ+(xRDXγ/2u,0RDxγ/2v)L2(𝒬)],
and the functional ℱ(·) is given by
(30)ℱ(v):=〈f-,v〉L2(𝒬)′.Lemma 4 (see [17]).
For realθ>0, v∈C0∞(R), then
(31)(-∞RDtθv(t),tRD∞θv(t))L2(R)=cos(πθ)∥-∞RDtθv(t)∥L2(R)2,(-∞RDtθv(t),tRD∞θv(t))L2(R)=cos(πθ)∥tRD∞θv(t)∥L2(R)2.Based on Lemma4, we can prove the following existence and uniqueness theorem. During this paper, we use the expression A≲B(A≳B) to mean that there exists a positive real number c such that A≤cB (A≥cB). At the same time, we denote A≅B to mean that A≲B≲A, which means there exist positive real numbers c1,c2 such that A≤c1B and B≤c2A (i.e., (1/c1)A≤B≤c2A).Theorem 5.
Assume that0<α1<⋯<αs<α<1,0<β<1, 1<γ<2 and f∈Bα/2,γ/2(I×Ω)′. Then, system (28) has a unique solution in Bα/2,γ/2(I×Ω). Furthermore,
(32)∥u∥Bα/2,γ/2(I×Ω)≲∥f-∥Bα/2,γ/2(I×Ω)′.Proof.
The existence and uniqueness of the solution is guaranteed by the well-known Lax-Milgram theorem. First, from the equivalence oflH0α(I×Ω), rH0α(I×Ω), and H0α(I×Ω), for all u,v∈Bα/2,γ/2(I×Ω), it follows that
(33)|𝒜(u,v)|=|+(xRDXγ/2u,0RDxγ/2v)L2(𝒬)](0RDtα/2u,tRDTα/2v)L2(𝒬)+∑i=1s‍ai(0RDtαi/2u,tRDTαi/2v)L2(𝒬)+kβCβ[(0RDxβ/2u,xRDXβ/2v)L2(𝒬)+kβCβ+(xRDXβ/2u,0RDxβ/2v)L2(𝒬)]+kγCγ[(0RDxγ/2u,xRDXγ/2v)L2(𝒬)+kγCγ+(xRDXγ/2u,0RDxγ/2v)L2(𝒬)]|≲∥0RDtα/2u∥L2(𝒬)∥0RDtα/2v∥L2(𝒬)+∑i=1s‍ai∥0RDtαi/2u∥L2(𝒬)∥0RDtαi/2v∥L2(𝒬)+kβCβ[∥0RDxβ/2u∥L2(𝒬)∥xRDXβ/2v∥L2(𝒬)+kβCβ+∥xRDXβ/2u∥L2(𝒬)∥0RDxβ/2v∥L2(𝒬)]+kγCγ[∥0RDxγ/2u∥L2(𝒬)∥xRDXγ/2v∥L2(𝒬)+kγCγ+∥xRDXγ/2u∥L2(𝒬)∥0RDxγ/2v∥L2(𝒬)]≲∥u∥Hα/2(I,L2(Ω))∥v∥Hα/2(I,L2(Ω))+∑i=1s‍ai∥u∥Hαi/2(I,L2(Ω))∥v∥Hαi/2(I,L2(Ω))+kβ∥u∥L2(I,Hβ/2(Ω))∥v∥L2(I,Hβ/2(Ω))+kγ∥u∥L2(I,Hγ/2(Ω))∥v∥L2(I,Hγ/2(Ω))≲∥u∥Bα/2,γ/2(I×Ω)∥v∥Bα/2,γ/2(I×Ω),(34)|ℱ(v)|≲∥f-∥Bα/2,γ/2(I×Ω)′∥v∥Bα/2,γ/2(I×Ω).
This implies the continuity of the bilinear form𝒜(·,·) and the right-hand side function ℱ(v).
We next prove the coercivity of the bilinear operator𝒜(·,·). Note that
(35)(0RDtα/2φ,tRDTα/2φ)L2(𝒬)=(-∞RDtα/2φ~,tRD∞α/2φ~)L2(R)=cos(πα2)·∥-∞RDtα/2φ~∥L2(R)=cos(πα2)·∥0RDtα/2φ∥L2(𝒬),
for all φ∈C0∞(𝒬), where φ~ is the extension of φ by zero outside of (0,T). Thus, we find that (0RDtα/2v,tRDTα/2v)L2(𝒬) is nonnegative for v∈Hα/2(0,T),0<α<1 since cos(πα/2) is nonnegative for 0<α<1. That is the same for fractional operators with αi,i=1,…,s.
From the above analysis, we have(36)𝒜(v,v)=(0RDtα/2v,tRDTα/2v)L2(𝒬)+∑i=1s‍ai(0RDtαi/2v,tRDTαi/2v)L2(𝒬)+2kβCβ(0RDxβ/2v,xRDXβ/2v)L2(𝒬)+2kγCγ(0RDxγ/2v,xRDXγ/2v)L2(𝒬)≅cos(πα2)(0RDtα/2v,0RDtα/2v)L2(𝒬)+∑i=1s‍aicos(παi2)(0RDtαi/2v,0RDtαi/2v)L2(𝒬)+kβ(0RDxβ/2v,0RDxβ/2v)L2(𝒬)+kγ(0RDxγ/2v,0RDxγ/2v)L2(𝒬)≳∥v∥Bα/2,γ/2(I×Ω)2.
By using the well-known Lax-Milgram theorem, there exists a unique solutionu∈Bα/2,γ/2(I×Ω) such that (28) holds.
To prove the stability estimate (32), by using (34) and (36), we take v=u in (28) and obtain
(37)∥u∥Bα/2,γ/2(I×Ω)2≲𝒜(u,u)=ℱ(u)≲∥f-∥Bα/2,γ/2(I×Ω)′∥u∥Bα/2,γ/2(I×Ω),
which implies that ∥u∥Bα/2,γ/2(I×Ω)≲∥f-∥Bα/2,γ/2(I×Ω)′.


## 3. Time Discretization

In this section, we consider the Diethelm fractional backward difference method based on quadrature, which was independently introduced by [39], for ordinary fractional differential equations. Here, we consider this method for the time discretization of (5) and derive the convergence rate for the time-discretization of MT-TS-RFADEs.Lemma 6 (see [25]).
Forn-1<θ<n(n=1,2,…) and a function u(x) defined on an unbounded domain (-∞<x<∞), the following equality holds:
(38)RD|x|θu(x)=-12cos(πθ/2)(-∞RDxθu+xRD∞θu)=-(-Δ)θ/2u(x),
where Δ is the Laplacian.Definition 7 (see [26]).
Suppose the Laplacian has a complete set of orthonormal eigenfunctionsϕn corresponding to eigenvalues λn2 on a bounded region 𝒬, that is,-Δϕn=λn2ϕn on a bounded region 𝒬; B(ϕ)=0 on ∂𝒬 is one of the standard three homogeneous boundary conditions. Let
(39)F={f=∑n=1∞‍cnϕn,cn=〈f,ϕn〉,∑n=1∞‍|cn|2|λn|θ≤∞},
then for any f∈F, (-Δ)θ/2 is defined by
(40)(-Δ)θ/2f=∑n=1∞‍cn(λn2)θ/2ϕn.Here, for anyt∈I, we make a continuous prolongation in space for function u(t,x) as follows:
(41)u~(t,x)={u(t,x),x∈[0,X],0,x∈(-∞,0)∪(X,∞).
Let u(t) denote the one-variable function of u~(t,x), and let f(t) and u(0) denote the one-variable functions as f(t,·) and u~(0,·), respectively. By the definition (12) of P(0Dt) and Lemma 6, (5) is equivalent to
(42)0CDtαu(t)+∑i=1s‍ai0CDtαiu(t)+kβ(-Δ)β/2u(t)+kγ(-Δ)γ/2u(t)=f(t).LetA=-Δ, D(A)=H01(Ω)∩H2(Ω). Then, the system (5) and the initial condition (6) can be rewritten in the abstract form, for 0≤t≤T,
(43)0CDtαu(t)+∑i=1s‍ai0CDtαiu(t)+kβAβ/2u(t)+kγAγ/2u(t)=f(t),(44)u(0)=ϕ0.
Based on (25), (43) can be rewritten as
(45)0RDtα[u-ϕ0](t)+∑i=1s‍ai0RDtαi[u-ϕ0](t)+kβAβ/2u(t)+kγAγ/2u(t)=f(t),
where 0<t<T,0<αs<⋯<α1<α<1.Let0=t0<t1<⋯<tN=T be a partition of [0,T]. Without loss of generality, for 0<α<1 and fix tn,n=1,2,…,N, we have
(46)0RDtα[u-ϕ0](tj)=tj-αΓ(-α)∫01‍g(θ)θ-1-αdθ,
where g(θ)=u(tj-tjθ)-ϕ0. Here, the integral is a Hadamard finite-part integral, in [38, 39].Now, for everyn, we replace the integral by a first-degree compound quadrature formula with equispaced nodes 0,1/n,2/n,…,1 and obtain
(47)∫01‍g(θ)θ-1-αdθ=∑k=0n‍αkn(α)g(kn)+Rn(α)(g),0<α<1,
where the weights αkn(α) are
(48)α(1-α)n-ααkn(α)={-1,k=0,2k1-α-(k-1)1-α-(k+1)1-α,fork=1,2,…,n-1,(α-1)k-α-(k-1)1-α+k1-α,fork=n,
and the remainder term Rn(α)(g) satisfies
(49)∥Rn(α)(g)∥≤γαnα-2sup0≤t≤T∥g′′(t)∥,
where γα>0 is a constant.Thus, we have(50)0RDtα[u-ϕ0](tn)=tn-αΓ(-α)(∑k=0n‍αkn(α)(u(tn-tk)-ϕ0)+Rn(α)(g))=Δt-α∑k=0n‍ωkn(α)(u(tn-tk)-ϕ0)+tn-αΓ(-α)Rn(α)(g),
where ωkn(α)=n-ααkn(α)/Γ(-α).Lett=tn, we can write (45) as
(51)Δt-α∑k=0n‍ωkn(α)(u(tn-tk)-ϕ0)+∑i=1s‍aiΔt-αi∑k=0n‍ωkn(αi)(u(tn-tk)-ϕ0)+kβAβ/2u(tn)+kγAγ/2u(tn)=f(tn)-tn-αΓ(-α)Rn(α)(g)-∑i=1s‍aitn-αiΓ(-αi)Rn(αi)(g),
for n=1,2,3,…, where ωkn(αi)=n-αiαkn(αi)/Γ(-αi).DenoteUn as the approximation of u(tn). We can define the following time stepping method:
(52)Δt-α∑k=0n‍ωkn(α)(Un-k-U0)+∑i=1s‍aiΔt-αi∑k=0n‍ωkn(αi)(Un-k-U0)+kβAβ/2Un+kγAγ/2Un=fn,
where fn=f(tn).Lemma 8 (see [38]).
For0<α<1, let the sequence {dj}j=1,2,… be given by d1=1 and dj=1+α(1-α)j-α∑k=1j-1αkjdj-k,j=2,3,…. Then, 1≤dj≤(sin(πα))/(πα(1-α))jα,j=1,2,3,….Leten=Un-u(tn) denote the error in tn. Let ∥·∥ be any given norm. Then, we have the following error estimate.Theorem 9.
Letu(tn) and Un be the solutions of (43) and (52), respectively, and Δt is the time step size. Then, one has
(53)maxn=0,1,…,N∥Un-u(tn)∥=O(Δt2-α+∑i=0s‍aiΔt2-αi).Proof.
Subtracting (52) from (51), we obtain the error equation
(54)Δt-α∑k=0n‍ωkn(α)(en-k-e0)+∑i=1s‍aiΔt-αi∑k=0n‍ωkn(αi)(en-k-e0)+kβAβ/2en+kγAγ/2en=-tn-αΓ(-α)Rn(α)(g)-∑i=1s‍aitn-αiΓ(-αi)Rn(αi)(g).
Note that e0=ϕ0-U0=0 and denote
(55)en=(Δt-αω0n(α)+∑i=1s‍aiΔt-αiω0n(αi)+kβAβ/2+kγAγ/2∑i=1s‍)-1×(Δt-α∑k=1n‍ωkn(α)en-k×+∑i=1s‍aiΔt-αi∑k=1n‍ωkn(αi)en-k×-tn-αΓ(-α)Rn(α)(g)×-∑i=1s‍aitn-αiΓ(-αi)Rn(αi)(g))=(-∑i=1s‍aiΓ(-α)Γ(-αi)α0n(α)+∑i=1s‍aiΔtα-αiΓ(-α)Γ(-αi)α0n(αi)+kβAβ/2tnαΓ(-α)+kγAγ/2tnαΓ(-α)+∑i=1s‍aiΔtα-αiΓ(-α)Γ(-αi))-1×(-∑i=1s‍aiΓ(-α)Γ(-αi)∑k=1n‍αkn(α)en-k+∑i=1s‍aiΔtα-αiΓ(-α)Γ(-αi)××∑k=1n‍αkn(αi)en-k-Rn(α)(g)×-∑i=1s‍aiΓ(-α)Γ(-αi)tnα-αiRn(αi)(g)).
Thus, we have
(56)∥en∥≤∥+kγAγ/2tnαΓ(-α)+∑i=1s‍aiΔtα-αiΓ(-α)Γ(-αi))-1(α0n(α)+∑i=1s‍aiΔtα-αiΓ(-α)Γ(-αi)α0n(αi)+kβAβ/2tnαΓ(-α)+kγAγ/2tnαΓ(-α)+∑i=1s‍aiΔtα-αiΓ(-α)Γ(-αi))-1∥×(∑i=1s‍aiΓ(-α)Γ(-αi)∑k=1n‍αkn(α)∥en-k∥×+∑i=1s‍aiΔtα-αiΓ(-α)Γ(-αi)∑k=1j‍αkn(αi)∥en-k∥×+∥Rn(α)(g)∥+∑i=1s‍aiΓ(-α)Γ(-αi)tnα-αi∥Rn(αi)(g)∥).
From Definition7, note that A is a positive definite elliptic operator with all of eigenvalues λ>0. Since α0n(αi)<0,i=1,…,s,α0n(α)<0 and Γ(-α)<0 when 0<α<1, we have
(57)∥+kβAβ/2tnαΓ(-α)+kγAγ/2tnαΓ(-α)∑i=1s‍aiΔtα-αiΓ(-α)Γ(-αi))-1(α0n(α)+∑i=1s‍aiΔtα-αiΓ(-α)Γ(-αi)α0n(αi)+kβAβ/2tnαΓ(-α)+kγAγ/2tnαΓ(-α)∑i=1s‍aiΔtα-αiΓ(-α)Γ(-αi))-1∥=supλ>0∥(+λβ/2tnαΓ(-α)+λγ/2tnαΓ(-α)∑i=1s‍aiΔtα-αiΓ(-α)Γ(-αi))-1(α0n(α)+∑i=1s‍aiΔtα-αiΓ(-α)Γ(-αi)α0n(αi)=supλ>0+λβ/2tnαΓ(-α)+λγ/2tnαΓ(-α)∑i=1s‍aiΔtα-αiΓ(-α)Γ(-αi))-1∥<(-α0n(α)-∑i=1s‍aiΔtα-αiΓ(-α)Γ(-αi)α0n(αi))-1.
Hence,
(58)∥en∥≤(-α0n(α))-1∑k=1n‍αkn(α)∥en-k∥+∑i=1s‍(-α0n(αi))-1∑k=1n‍αkn(αi)∥en-k∥+(-α0n(α))-1∥Rn(α)(g)∥+∑i=1s‍(-α0n(αi))-1∥Rn(αi)(g)∥≤α(1-α)n-α∑k=1n‍αkn(α)∥en-k∥+∑i=1s‍αi(1-αi)n-αi∑k=1n‍αkn(αi)∥en-k∥+α(1-α)γαN-2sup0≤t≤T∥u′′∥+∑i=1s‍αi(1-αi)γαiN-2sup0≤t≤T∥u′′∥.
Note that e0=ϕ0-U0=0. Denote d1=1 and
(59)dn=1+α(1-α)n-α∑k=1n-1‍αkn(α)dn-k,n=2,3,…,N,dni=1+αi(1-αi)n-αi∑k=1j-1‍αkn(αi)dn-k,n=2,3,…,N,i=1,2,…,s.
Then, by induction and Lemma8, we have
(60)∥en∥≲α(1-α)N-2sup0≤t≤T∥u′′(t)∥·dn+∑i=1s‍aiαi(1-αi)N-2sup0≤t≤T∥u′′(t)∥·dni≲N-2sup0≤t≤T∥u′′(t)∥sin(πα)πnα+∑i=1s‍aiN-2sup0≤t≤T∥u′′(t)∥sin(παi)πnαi=O(Δt2-α+∑i=1s‍aiΔt2-αi).


## 4. Space Discretization

In this section, we consider the space discretization of (5) with homogeneous boundary condition. Using the FEM, we obtain the numerical approximation solution in a finite domain. Then, we prove the convergence rate of this method. Let Ω=[0,X] be an interval in one-dimensional space. All of the results in this section can be generalized into the cases of high dimension.Based on the time discretization in Section3 and (25), we need to find u(t,·)∈H0γ/2(Ω) such that
(61)(0RDtαu(t,x),v)L2(Ω)+∑i=1s‍ai(0RDtαiu(t,x),v)L2(Ω)+kβCβ[(0RDxβ/2u(t,x),xRDXβ/2v)L2(Ω)+kβCβ+(xRDXβ/2u(t,x),0RDxβ/2v)L2(Ω)]+kγCγ[(0RDxγ/2u(t,x),xRDXγ/2v)L2(Ω)+kγCγ+(xRDXγ/2u(t,x),0RDxγ/2v)L2(Ω)]=(f-(t,x),v)L2(Ω),∀v∈H0γ/2(Ω).Leth denote the maximal length of intervals in Ω, and let r be any nonnegative integer. We denote the norm in Hr(Ω) by ∥·∥Hr(Ω). Let Sh⊂H0γ/2 be a family of finite element space with the accuracy of order r≥2;thatis,Sh consists of continuous functions on the closure Ω- of Ω which are polynomials of degree at most r-1 in each interval and which vanish outside Ωh such that for small h, v∈Hη(Ω)∩H0γ/2(Ω),
(62)infχ∈Sh(∥(v-χ)∥L2(Ω)+h∥∇(v-χ)∥L2(Ω))≲hη∥v∥Hη(Ω),1≤η≤r.The semidiscrete problem of (5) is to find the approximate solution uh(t)=uh(t,·)∈Sh for each t such that
(63)(0RDtαuh(t),χ)L2(Ω)+∑i=1s‍ai(0RDtαiuh(t),χ)L2(Ω)+kβCβ[(0RDxβ/2uh(t),xRDXβ/2χ)L2(Ω)+kβCβ+(xRDXβ/2uh(t),0RDxβ/2χ)L2(Ω)]+kγCγ[(0RDxγ/2uh(t),xRDXγ/2χ)L2(Ω)+kγCγ+(xRDXγ/2uh(t),0RDxγ/2χ)L2(Ω)]=(f-(t),χ)L2(Ω),∀χ∈Sh.
Let Rh:Hγ/2(Ω)→Sh be the elliptic projection defined by
(64)kβCβ(0RDxβ/2Rhu,xRDXβ/2χ)L2(Ω)+kβCβ(xRDXβ/2Rhu,0RDxβ/2χ)L2(Ω)+kγCγ(0RDxγ/2Rhu,xRDXγ/2χ)L2(Ω)+kγCγ(xRDXγ/2Rhu,0RDxγ/2χ)L2(Ω)=kβCβ(0RDxβ/2u,xRDXβ/2χ)L2(Ω)+kβCβ(xRDXβ/2u,0RDxβ/2χ)L2(Ω)+kγCγ(0RDxγ/2u,xRDXγ/2χ)L2(Ω)+kγCγ(xRDXγ/2u,0RDxγ/2χ)L2(Ω),∀χ∈Sh.Lemma 10.
Assume that (62) holds. Then, for Rh defined by (64) and any v∈Hη(Ω)∩H0θ/2(Ω), one has
(65)∥0RDxθ/2(Rhv-v)∥L2(Ω)≲hη-θ/2∥v∥Hη(Ω).Proof.
LetIh be the projection operator from Hη(Ω)∩H0θ/2(Ω) to Sh. From the definition of L2-norm, we obtain
(66)(0RDxθ/2(Rhv-v),xRDXθ/2(Rhv-v))L2(Ω)+(xRDXθ/2(Rhv-v),0RDxθ/2(Rhv-v))L2(Ω)=(0RDxθ/2(Rhv-v),xRDXθ/2(Rhv-Ihv))L2(Ω)+(xRDXθ/2(Rhv-v),0RDxθ/2(Rhv-Ihv))L2(Ω)+(0RDxθ/2(Rhv-v),xRDXθ/2(Ihv-v))L2(Ω)+(xRDXθ/2(Rhv-v),0RDxθ/2(Ihv-v))L2(Ω).
Let χ=Rhv-Ihv∈Sh. From (64), we obtain
(67)(0RDxθ/2(Rhv-v),xRDXθ/2(Rhv-v))L2(Ω)+(xRDXθ/2(Rhv-v),0RDxθ/2(Rhv-v))L2(Ω)≤∥0RDxθ/2(Rhv-v)∥L2(Ω)∥xRDXθ/2(Ihv-v)∥L2(Ω)+∥xRDXθ/2(Rhv-v)∥L2(Ω)∥0RDxθ/2(Ihv-v)∥L2(Ω)≤(∥0RDxθ/2(Rhv-v)∥L2(Ω)+∥xRDXθ/2(Rhv-v)∥L2(Ω))×(∥xRDXθ/2(Ihv-v)∥L2(Ω)+∥0RDxθ/2(Ihv-v)∥L2(Ω)).
Note that(68)(0RDxθ/2(Rhv-v),xRDXθ/2(Rhv-v))L2(Ω)+(xRDXθ/2(Rhv-v),0RDxθ/2(Rhv-v))L2(Ω)≳(∥0RDxθ/2(Rhv-v)∥L2(Ω)+∥xRDXθ/2(Rhv-v)∥L2(Ω))2.
Thus, combining (67) and (68), we obtain
(69)∥0RDxθ/2(Rhv-v)∥L2(Ω)+∥xRDXθ/2(Rhv-v)∥L2(Ω)≲∥0RDxθ/2(Ihv-v)∥L2(Ω)+∥xRDXθ/2(Ihv-v)∥L2(Ω).
By (62), we see that
(70)∥v-Ihv∥L2(Ω)+h∥∇(v-Ihv)∥L2(Ω)≲hη|v|Hη(Ω).
By interpolation properties, we obtain
(71)∥0RDxθ/2(v-Ihv)∥L2(Ω)≤|(v-Ihv)|Hθ/2(Ω)≲hη-θ/2|v|Hη/2(Ω).
Similarly, we have
(72)∥xRDXθ/2(v-Ihv)∥L2(Ω)≤|(v-Ihv)|Hθ/2(Ω)≲hη-θ/2|v|Hη/2(Ω).
Combining (71) and (72), we can obtain (65).Theorem 11.
For0<αs<⋯<α1<α<1, let uh and u be the solutions of (63) and (61), respectively. Then, it holds
(73)∥u(t,x)-uh(t,x)∥Hγ/2(Ω)=O(Δxη-γ/2).Proof.
We write(74)uh-u=ϑ+ρ,
where ϑ=uh-Rhu, ρ=Rhu-u. The second term is easily bounded by Lemma 10 and has the obvious estimate
(75)∥ρ(x)∥Hγ/2(Ω)2≲h2η-γ∥u∥Hη(Ω)2.
In order to estimateϑ, for all χ∈Sh, we get
(76)(0RDtαϑ(t),χ)L2(Ω)+∑i=1s‍ai(0RDtαiϑ(t),χ)L2(Ω)+kβCβ[(0RDxβ/2ϑ(t),xRDXβ/2χ)L2(Ω)+kβCβ+(xRDXβ/2ϑ(t),0RDxβ/2χ)L2(Ω)]+kγCγ[(0RDxγ/2ϑ(t),xRDXγ/2χ)L2(Ω)+kγCγ+(xRDXγ/2ϑ(t),0RDxγ/2χ)L2(Ω)]=(0RDtα(uh(t)-Rhu(t)),χ)L2(Ω)+∑i=1s‍ai(0RDtαi(uh(t)-Rhu(t)),χ)L2(Ω)+kβCβ[(0RDxβ/2(uh(t)-Rhu(t)),xRDXβ/2χ)L2(Ω)+kβCβ+(xRDXβ/2(uh(t)-Rhu(t)),0RDxβ/2χ)L2(Ω)]+kγCγ[(0RDxγ/2(uh(t)-Rhu(t)),xRDXγ/2χ)L2(Ω)+kγCγ+(xRDXγ/2(uh(t)-Rhu(t)),0RDxγ/2χ)L2(Ω)]=(0RDtαuh(t),χ)L2(Ω)+∑i=1s‍ai(0RDtαiuh(t),χ)L2(Ω)+kβCβ(0RDxβ/2uh(t),xRDXβ/2χ)L2(Ω)+kβCβ(xRDXβ/2uh(t),0RDxβ/2χ)L2(Ω)+kγCγ(0RDxγ/2uh(t),xRDXγ/2χ)L2(Ω)+kγCγ(xRDXγ/2uh(t),0RDxγ/2χ)L2(Ω)-(+∑i=1s‍(0RDtαRhu(t),χ)L2(Ω)-+∑i=1s‍ai(0RDtαiRhu(t),χ)L2(Ω)-+kβCβ(0RDxβ/2Rhu(t),xRDXβ/2χ)L2(Ω)-+kβCβ(xRDXβ/2Rhu(t),0RDxβ/2χ)L2(Ω)-+kγCγ(0RDxγ/2Rhu(t),xRDXγ/2χ)L2(Ω)-+kγCγ(xRDXγ/2Rhu(t),0RDxγ/2χ)L2(Ω)+∑i=1s‍)=(f-(t),χ)L2(Ω)-(+∑i=1s‍(0RDtαRhu(t),χ)L2(Ω)-+∑i=1s‍ai(0RDtαiRhu(t),χ)L2(Ω)-+kβCβ(0RDxβ/2Rhu(t),xRDXβ/2χ)L2(Ω)-+kβCβ(xRDXβ/2Rhu(t),0RDxβ/2χ)L2(Ω)-+kγCγ(0RDxγ/2Rhu(t),xRDXγ/2χ)L2(Ω)-+kγCγ(xRDXγ/2Rhu(t),0RDxγ/2χ)L2(Ω)+∑i=1s‍)=(P(0Dt)u(t),χ)L2(Ω)-(P(0Dt)Rhu(t),χ)L2(Ω)=((I-Rh)P(0Dt)u(t),χ)L2(Ω)=-(P(0Dt)ρ(t),χ)L2(Ω).
Choosingχ=ϑ(t) and integrating on both sides with respect to t on [0,T], we obtain
(77)(P(0Dt)ϑ(t),ϑ(t))L2(Ω)+kβCβ(0RDxβ/2ϑ(t),xRDXβ/2ϑ(t))L2(Ω)+kβCβ(xRDXβ/2ϑ(t),0RDxβ/2ϑ(t))L2(Ω)+kγCγ(0RDxγ/2ϑ(t),xRDXγ/2ϑ(t))L2(Ω)+kγCγ(xRDXγ/2ϑ(t),0RDxγ/2ϑ(t))L2(Ω)=-(P(0Dt)ρ(t),ϑ(t))L2(Ω).
By Lemmas2–4, for any small ϵ>0, we have
(78)Aα(0RDtα/2ϑ(t),0RDtα/2ϑ(t))L2(Ω)+Aαi∑i=1s‍ai(0RDtαi/2ϑ(t),0RDtαi/2ϑ(t))L2(Ω)+kβCβ[(0RDxβ/2ϑ(t),xRDXβ/2ϑ(t))L2(Ω)+kβCβ+(xRDXβ/2ϑ(t),0RDxβ/2ϑ(t))L2(Ω)]+kγCγ[(xRDXγ/2ϑ(t),0RDxγ/2ϑ(t))L2(Ω)(0RDxγ/2ϑ(t),xRDXγ/2ϑ(t))L2(Ω)+kγCγ+(xRDXγ/2ϑ(t),0RDxγ/2ϑ(t))L2(Ω)]=Aα(0RDtα/2ρ(t),0RDtα/2ϑ(t))L2(Ω)+∑i=1s‍aiAαi(0RDtαi/2ρ(t),0RDtαi/2ϑ(t))L2(Ω)≤Cϵ(+∑i=1s‍aiAαi∥0RDtαi/2ρ(t)∥L2(Ω)2Aα∥0RDtα/2ρ(t)∥L2(Ω)2≤Cϵ+∑i=1s‍aiAαi∥0RDtαi/2ρ(t)∥L2(Ω)2)+ϵ(∑i=1s‍aiAαi∥0RDtαi/2ϑ(t)∥L2(Ω)2Aα∥0RDtα/2ϑ(t)∥L2(Ω)2+ϵ+∑i=1s‍aiAαi∥0RDtαi/2ϑ(t)∥L2(Ω)2),
where Cϵ is a constant respect to ϵ, Aα=cos(πα/2) and Aαi=cos((παi)/2),i=1,…,s.
For sufficiently smallϵ>0, by (75), we obtain
(79)(P(0Dt)ϑ(t),ϑ(t))L2(Ω)≤Cϵ(+∑i=1s‍ai∥0RDtαi/2ρ(t)∥L2(Ω)2∥0RDtα/2ρ(t)∥L2(Ω)2≤Cϵ+∑i=1s‍ai∥0RDtαi/2ρ(t)∥L2(Ω)2)≤Cϵ(hη-γ/2)2(+∑i=1s‍ai∥0RDtαi/2u(t)∥Hη(Ω)2∥0RDtα/2u(t)∥Hη(Ω)2≤Cϵ(hη-γ/2)2+∑i=1s‍ai∥0RDtαi/2u(t)∥Hη(Ω)2).
Combining (75) with (79), we obtain that
(80)∥u-uh∥Hγ/2(Ω)=∥ρ+ϑ∥Hγ/2(Ω)=O(hη-γ/2).For the reason that the time and space fractional derivatives, we introduce the complete form of this FEM. In view of space discretization, we first pose the finite-dimensional problem to finduh(t,·)∈Sh such that (63) holds.In terms of the basis{ψi}i=1M-1⊆Sh, write uh(t,x)=∑j=1M-1‍uj(t)ψj(x), and insert it into (63). After the time discretization on tn, in Section 3, we get
(81)∑j=1M-1‍Δt-α∑k=1n‍ωkn(α)Ujn-k(ψj,ψi)L2(Ω)+∑i=1s‍ai∑j=1M-1‍Δt-αi∑k=1n‍ωkn(αi)Ujn-k(ψj,ψi)L2(Ω)+kβCβ∑j=1M-1‍Ujn[(0RDxβ/2ψj,xRDXβ/2ψi)L2(Ω)+kβCβ∑j=1M-1‍Ujn+(xRDXβ/2ψj,0RDxβ/2ψi)L2(Ω)]+kγCγ∑j=1M-1‍Ujn[(0RDxγ/2ψj,xRDXγ/2ψi)L2(Ω)+kγCγ∑j=1M-1‍Ujn+(xRDXγ/2ψj,0RDxγ/2ψi)L2(Ω)]=(f-,ψi)L2(Ω),∀ψi∈Sh,i=1,…,M-1.To obtain the value of{Ujn}j=1M-1, Let Uj=(Uj1,Uj2,…,Ujn)T. From (81), we obtain a vector equation
(82)Δt-αω0n(α)MUn+∑i=1s‍aiΔt-αiω0n(αi)MUn+kβCβ(S1(β/2)+S2(β/2))Un+kγCγ(S1(γ/2)+S2(γ/2))Un=MFn-Δt-α∑k=1n‍ωkn(α)MUn-k-∑i=1s‍Δt-αi∑k=1n‍ωkn(αi)MUn-k,
where M={(ψj,ψi)L2(Ω)}j,i=1M-1 is the mass matrix; S1(θ) and S2(θ) are stiffness matrices with θ=β/2,γ/2, as follows:
(83)S1(θ)={(0RDxθψj,xRDXθψi)L2(Ω)}j,i=1M-1,S2(θ)={(xRDXθψj,0RDxθψi)L2(Ω)}j,i=1M-1,
and FN=(f-1,…,f-M-1)T is a vector valued function. Then, from (82), we can obtain the solution Un.


## 5. Stability of the Numerical Method

In this section, we analyze the stability of the FEM for MT-TS-RFADEs (5). Now, we do some preparation before proof. Based on the definition of coefficients ωkn(α) in Section 3, we can obtain the following lemma easily.Lemma 12.
For0<α<1, the coefficients ωkn(α), (k=1,…,n) satisfy the following properties: (i)
ω0n(α)>0 and ωkn(α)<0 for k=1,2,…,n,(ii)
Γ(2-α)∑k=1n‍ωkn(α)=(1-α)n-α+1.Now, we report the stability theorem of this FEM for MT-TS-RFADEs as follows.Theorem 13.
The FEM defined in (81) is unconditionally stable.Proof.
LetUn denote the approximation to uh at t=tn and χ(·)=Un(·) and the right-hand side f-=0. From (63), we have
(84)Δt-α{1Γ(2-α)(Un,Un)L2(Ω)Δt-α+∑k=1n-1‍ωkn(α)(Un-k,Un)L2(Ω)Δt-α+ωnn(α)(U0,Un)L2(Ω)1Γ(2-α)}+∑i=1s‍aiΔt-αi{1Γ(2-αi)(Un,Un)L2(Ω)+∑i=1s‍aiΔt-αi+∑k=1n-1‍ωkn(αi)(Un-k,Un)L2(Ω)+∑i=1s‍aiΔt-αi+ωnn(αi)(U0,Un)L2(Ω)1Γ(2-αi)}+kβCβ[(xRDXβ/2Un,0RDxβ/2Un)L2(Ω)(0RDxβ/2Un,xRDXβ/2Un)L2(Ω)+kβCβ+(xRDXβ/2Un,0RDxβ/2Un)L2(Ω)]+kγCγ[(xRDXγ/2Un,0RDxγ/2Un)L2(Ω)(0RDxγ/2Un,xRDXγ/2Un)L2(Ω)+kγCγ+(xRDXγ/2Un,0RDxγ/2Un)L2(Ω)]=0.
Using the Cauchy-Schwarz inequality, ±(Un-k,Un)≤(1/2)(∥Un-k∥L2(Ω)2+∥Un∥L2(Ω)2) for k=0,1,2,…,n, by Definition 1 and Lemma 12, we obtain
(85)(+∑i=1s‍aiΔt-αi2Γ(2-αi)Δt-α2Γ(2-α)(1+(1-α)n-α)+∑i=1s‍aiΔt-αi2Γ(2-αi)(1+(1-αi)n-αi))∥Un∥L2(Ω)2+kβ∥0RDxβ/2Un∥L2(Ω)2+kγ∥0RDxγ/2Un∥L2(Ω)2≤Δt-α2[-∑k=1n-1‍ωkn(α)∥Un-k∥L2(Ω)2-ωnn(α)∥U0∥L2(Ω)2]+∑i=1s‍aiΔt-αi2[-∑k=1n-1‍ωkn(αi)∥Un-k∥L2(Ω)2-ωnn(αi)∥U0∥L2(Ω)2].
We prove the stability of (63) by induction. From the beginning, we have
(86)(+∑i=1s‍aiΔt-αi2Γ(2-αi)Δt-α2Γ(2-α)(1+(1-α))+∑i=1s‍aiΔt-αi2Γ(2-αi)(1+(1-αi)))∥U1∥L2(Ω)2≤(+∑i=1s‍aiΔt-αi2Γ(2-αi)Δt-α2Γ(2-α)(1-(1-α))≤+∑i=1s‍aiΔt-αi2Γ(2-αi)(1-(1-αi)))∥U0∥L2(Ω)2.
The induction basis ∥U1∥L2(Ω)≤∥U0∥L2(Ω) is presupposed. For the induction step, we have ∥Un∥L2(Ω)≤∥Un-1∥L2(Ω)≤⋯≤∥U0∥L2(Ω). Then, using this result, by Lemma 12, we obtain
(87)(+∑i=1s‍aiΔt-α(1+(1-α)(n+1)-α)2Γ(2-α)+∑i=1s‍aiΔt-αi(1+(1-αi)(n+1)-αi)2Γ(2-αi))∥Un+1∥L2(Ω)2≤(+∑i=1s‍Δt-α(1-(1-α)(n+1)-α)2Γ(2-α)≤+∑i=1s‍aiΔt-αi(1-(1-αi)(n+1)-αi)2Γ(2-αi))∥U0∥L2(Ω)2.
Here, 0<1-α<1 and 0<1-αi<1 for i=1,…,s. After squaring at both sides of the above inequality, we obtain ∥Un+1∥L2(Ω)≤∥U0∥L2(Ω).


## 6. Numerical Tests

Based on the above analysis, we present three numerical examples for MT-TS-RFADEs to demonstrate the efficiency of our theoretical analysis. The main purpose is to check the convergence behavior of numerical solutions with respect to time step sizeΔt and space step size Δx, which have been shown in Theorems 9 and 11.Example 1.
Consider MT-TS-RFADE, int∈[0,0.5],x∈[0,1],
(88)0CDtαu(t,x)+0CDtα1u(t,x)-RD|x|βu(t,x)-RD|x|γu(t,x)=f(t,x),u(0,x)=x2(1-x)2,x∈(0,1),u(t,0)=u(t,1)=0,t∈[0,0.5].
The exact solution of (88) is u(t,x)=(2t-1)2x2(1-x)2. From the definition of the Riemann-Liouville differential operator, it holds
(89)aRDtα(t-a)β=Γ(β+1)Γ(β+1-α)(t-a)β-α,
where a is a positive constant. From (89), we can choose right-hand side function f(t,x) to satisfy (88).
Choosingα=0.9 and α1=0.2 in time fractional operators and β=0.4γ=1.6 in the space Riesz fractional operators, we can obtain the numerical approximation to the exact solution of (88) on finite domain [0,0.5]×[0,1], with space step size Δx=0.05 and time step size Δt=0.01. In Figure 1, one can see that the numerical solution matches well with the exact solution.Figure 1
Numerical solution and exact solution for Example1.Example 2.
Consider MT-TS-RFADE (88) with conditions as follows:
(90)u(0,x)=0,x∈[0,1],u(t,0)=u(t,1)=0,t∈[0,0.5].
Let the exact solution is u(t,x)=t0.9x2(1-x)2. We choose α1=0.3,α=0.4,β=0.4,γ=1.7 and obtain the numerical solution and exact solution when t=0.05,0.1,0.25. The results have been shown in Figure 2, where the exact solution is noted by lines and numerical solution is noted by squares. Here, space step size is Δx=0.1; time step size is Δt=0.01.Figure 2
Numerical solution and exact solution for Example2.Example 3.
Consider MT-TS-RFADE (88) with the zero Dirichlet boundary conditions, for t∈[0,1], x∈[0,1]. We require that the exact solution is u(t,x)=sin(2πt)x2(1-x)2.For this example, in the first test, we obtain the numerical solution and exact solution whent=0.05, 0.1, 0.25, 0.4 in Figure 3, where we choose time step size Δt=0.01 and space step size Δx=0.01 with α1=0.3, α=0.5, and β=0.9, γ=1.9.Figure 3
Numerical solution and exact solution for Example3.In the second test, we check the convergence rates of numerical solutions with respect to the fractional ordersα1,α,β,  and γ. We fix α1=0.2,α=0.8,β=0.8,  and γ=1.8 and choose Δx=0.001 which is small enough such that the space discretization errors are negligible as compared with the time errors. Choosing step size Δt=1/2i(i=1,…,5), we present Table 1 with the convergence rate which is equal to 1.2, as Theorem 9 predicted. Table 2 shows the spatial approximate convergence rate, by fixing Δt=0.001 and choosing Δx=1/2i(i=1,…,5). From Theorem 11, the convergence rate should be equal to or less than 1.1 (i.e., η-γ/2 for η=2 and γ=1.8). In Table 2, the numerical results match well with such conclusion. Here, we also report both the L2-norm and H1-norm of errors in Figure 4.Table 1
Convergence rate in time for Example3.
Δ
x
Δ
t
H
1-norm
L
2-norm
cvge. rate
0.001
1/2
7.9601
×
1
0
-
3
1.0925
×
1
0
-
3
0.001
1/4
3.4906
×
1
0
-
3
4.7906
×
1
0
-
4
1.1893
0.001
1/8
1.5238
×
1
0
-
3
2.0914
×
1
0
-
4
1.1957
0.001
1/16
6.6273
×
1
0
-
4
9.0956
×
1
0
-
5
1.2012
0.001
1/32
2.8667
×
1
0
-
4
3.9343
×
1
0
-
5
1.2090Table 2
Convergence rate in space for Example3.
Δ
t
Δ
x
H
1-norm
L
2-norm
cvge. rate
0.001
1/2
3.2591
×
1
0
-
2
4.4731
×
1
0
-
3
0.001
1/4
1.6288
×
1
0
-
2
2.2356
×
1
0
-
3
1.0006
0.001
1/8
8.1364
×
1
0
-
3
1.1167
×
1
0
-
3
1.0014
0.001
1/16
4.0602
×
1
0
-
3
5.5723
×
1
0
-
4
1.0029
0.001
1/32
2.0220
×
1
0
-
3
2.7751
×
1
0
-
4
1.0057H
1-norm and L2-norm of errors for Example 3; here, Δx=0.001 (a) and Δt=0.001 (b).
(a)
(b)



---

*Source: 868035-2013-04-15.xml*