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

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.

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: 𝑑𝑢()𝑑𝑡𝑡,𝑥=−𝑘_𝛼(−Δ)^𝛼/2𝑢(𝑡,𝑥)−𝑘_𝛽(−Δ)^𝛽/2𝑢(𝑡,𝑥), (1) for example, where 𝑢 is a solute concentration, 𝑘_𝛼 and 𝑘_𝛽 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 𝑢(𝑡,𝑥) 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 𝑥∈(−∞,∞) is equivalent to the Riesz fractional derivative operator ^𝑅𝐷^𝛼_|𝑥|. 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) where ^𝑅𝐷^𝛼_|𝑥| 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/2𝑦(𝑡)+𝐶𝑦(𝑡)=𝑓(𝑡), (3) where 𝐴,𝐵, and 𝐶 are certain constants and 𝑓 is a given function. The Basset equation is, in [6], 𝐷¹𝑦(𝑡)+𝑏𝐷^𝛼𝑦(𝑡)+𝑐𝑦(𝑡)=𝑓(𝑡),𝑦(0)=𝑦₀, (4) where 0<𝛼<1, 𝑏, and 𝑐 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 Section 2, 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) where 0<𝛽<1, 1<𝛾<2, 𝑥∈[0,𝑋], and 𝑡∈[0,𝑇] are respectively the space and time variables and 𝑘_𝛽, 𝑘_𝛾 are positive constants. We consider this problem with the zero Dirichlet boundary value conditions and the initial value condition defined as follows: 𝑢()()[],𝑢()(𝑥)().𝑡,0=𝑢𝑡,𝑋=0,𝑡∈0,𝑇0,𝑥=𝜙,𝑥∈0,𝑋 (6) 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 derivatives ^𝑅𝐷^𝛽_|𝑥|𝑢(𝑡,𝑥) and ^𝑅𝐷^𝛾_|𝑥|𝑢(𝑡,𝑥) 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 𝑛 and real number 𝜃 (𝑛−1<𝜃<𝑛), 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: ^𝐶₀𝐷^𝜃_𝑡𝑣(𝑡)=1Γ(𝑛−𝜃)^𝑡₀1(𝑡−𝜏)^{𝜃−𝑛+1}𝑑^𝑛𝑑𝜏^𝑛𝑣(𝜏)𝑑𝜏, (7)(ii)the right Caputo derivative: ^𝐶_𝑡𝐷^𝜃_𝑇𝑣(𝑡)=(−1)^𝑛Γ(𝑛−𝜃)^𝑇_𝑡1(𝜏−𝑡)^{𝜃−𝑛+1}𝑑^𝑛𝑑𝜏^𝑛𝑣(𝜏)𝑑𝜏, (8)(iii) the left Riemann-Liouville derivative: ^𝑅₀𝐷^𝜃_𝑡𝑣(𝑡)=1Γ()𝑑𝑛−𝜃𝑑𝑡^𝑛𝑛𝑡₀𝑣(𝜏)(𝑡−𝜏)^{𝜃−𝑛+1}𝑑𝜏, (9)(iv) the right Riemann-Liouville derivative: ^𝑅_𝑡𝐷^𝜃_𝑇𝑣(𝑡)=(−1)^𝑛Γ()𝑑𝑛−𝜃𝑑𝑡^𝑛𝑛𝑇_𝑡𝑣(𝜏)(𝜏−𝑡)^{𝜃−𝑛+1}𝑑𝜏. (10)

The Riesz fractional operators ^𝑅𝐷^𝛽_|𝑥| and ^𝑅𝐷^𝛾_|𝑥| 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 𝑛∈𝐍,𝑛−1<𝜃<𝑛, on a finite interval 0≤𝑥≤𝑋 is defined as ^𝑅𝐷^𝜃_|𝑥|𝑣(𝑡,𝑥)=−𝐶_𝜃^𝑅₀𝐷^𝜃_𝑥𝑣(𝑡,𝑥)+^𝑅_𝑥𝐷^𝜃_𝑋𝑣(𝑡,𝑥),𝑥∈[0,𝑋], (11) where 𝐶_𝜃=1/(2cos(𝜋𝜃/2)).

The multiterm fractional operator 𝑃(₀𝐷_𝑡) in (5) is 𝑃₀𝐷_𝑡𝑢(𝑡,𝑥)=^𝐶₀𝐷^𝛼_𝑡+^𝑠_𝑖=1𝑎_𝑖^𝐶₀𝐷^𝛼𝑖_𝑡𝑢(𝑡,𝑥), (12) where 0<𝛼_𝑠<𝛼_𝑠−1<⋯<𝛼₁<𝛼<1 or 1<𝛼_𝑠<𝛼_𝑠−1<⋯<𝛼₁<𝛼<2 or 0<𝛼_𝑠<⋯<𝛼_𝑠0<1<𝛼_𝑠0−1<⋯<𝛼₁<𝛼<2, and 𝑎_𝑖∈𝐑, 𝑖=1,…,𝑠,𝑠∈𝐍. Here, ^𝐶₀𝐷^𝛼_𝑡 and ^𝐶₀𝐷^𝛼𝑖_𝑡(𝑖=1,…,𝑠) denote the left Caputo fractional derivatives with respect to the time variable 𝑡 of order 0<𝛼<1 and 0<𝛼_𝑠<⋯<𝛼₁<1. There are three cases for multiterm time-space fractional derivatives.

Case 1. If 0<𝛼_𝑠<⋯<𝛼₁<𝛼<1, (5) is a generalized MT-TS-RFADE with multiterm time fractional diffusion terms with initial conditions given as (6). And especially, if 𝛼_𝑖=0, 𝑖=1,…,𝑠 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. If 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: 𝑢(0,𝑥)=𝜙(𝑥),𝑢_𝑡(0,𝑥)=𝜓(𝑥),𝑥∈(0,𝑋). (13)

Case 3. If 0<𝛼_𝑠<⋯<𝛼_𝑠0<1<𝛼_𝑠0−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 Case 1 of (5) with 0<𝛼_𝑠<⋯<𝛼₁<𝛼<1. Another two cases with 1<𝛼_𝑠<⋯<𝛼₁<𝛼<2 and 0<𝛼_𝑠<⋯<𝛼_𝑠0<1<𝛼_𝑠0−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 𝐶^∞(0,𝑇) denote the space of infinitely differentiable functions on (0,𝑇), and let 𝐶^∞₀(0,𝑇) denote the space of infinitely differentiable functions with compact support in (0,𝑇). Let 𝐿₂(𝒬) be the space of measurable functions whose square is the Lebesgue integrable in 𝒬, which may denote a domain 𝒬=𝐼 or Ω or 𝐼×Ω, where 𝐼=[0,𝑇] denotes the time domain and Ω=[0,𝑋] denotes the space domain. The inner product and norm of 𝐿₂(𝒬) are defined by (𝑢,𝑣)_𝐿2(𝒬)=_𝒬‖𝑢‖𝑢𝑣𝑑𝒬,_𝐿2(𝒬)=(𝑢,𝑢)_𝐿^1/2_2(𝒬),∀𝑢,𝑣∈𝐿₂(𝒬). (14)

For any real 𝜎>0, we define the spaces ^𝑙𝐻^𝜎₀(𝒬) and ^𝑟𝐻^𝜎₀(𝒬) to be the closure of 𝐶^∞₀(𝒬) with respect to the norms ‖𝑣‖_𝑙𝐻^𝜎_0(𝒬) and ‖𝑣‖_𝑟𝐻^𝜎_0(𝒬), respectively, where ‖𝑣‖_𝑙𝐻^𝜎_0(𝒬)∶=‖𝑣‖²_𝐿2(𝒬)+|𝑣|²_{𝑙𝐻𝜎0(𝒬)}^1/2,|𝑣|²_{𝑙𝐻𝜎0(𝒬)}∶=^𝑅₀𝐷^𝜎_𝑡𝑣²_𝐿2(𝒬),‖𝑣‖_𝑟𝐻^𝜎_0(𝒬)∶=‖𝑣‖²_𝐿2(𝒬)+|𝑣|²_𝑟𝐻^𝜎_0(𝒬)^1/2,|𝑣|²_𝑟𝐻^𝜎_0(𝒬)∶=^𝑅_𝑡𝐷^𝜎_𝑇𝑣²_𝐿2(𝒬). (15) In the usual Sobolev space 𝐻^𝜎₀(𝒬), we also have the definition ‖𝑣‖_𝐻^𝜎_0(𝒬)∶=‖𝑣‖²_𝐿2(𝒬)+|𝑣|²_𝐻^𝜎_0(𝒬)^1/2,|𝑣|²_𝐻^𝜎_0(𝒬)∶=^𝑅0𝐷^𝜎_𝑡𝑣,^𝑅_𝑡𝐷^𝜎_𝑇𝑣_𝐿2(𝒬)().cos𝜋𝜎 (16) From [18], for 𝜎>0,𝜎≠𝑛−1/2, the spaces ^𝑙𝐻^𝜎₀(𝒬), ^𝑟𝐻^𝜎₀(𝒬), and 𝐻^𝜎₀(𝒬) are equal, and their seminorms are all equivalent to |⋅|_𝐻^𝜎_0(𝒬).

We now give some results for fractional operators on these spaces.

Lemma 2 (see [8, 18]). (1) For real 0<𝜃<1,0<𝛿<1, if 𝑣(0)=0,𝑥∈(0,𝑋), then ^𝑅₀𝐷_𝑥^𝜃+𝛿𝑣(𝑥)=^𝑅₀𝐷^𝜃_𝑥^𝑅₀𝐷^𝛿_𝑥𝑣(𝑥)=^𝑅₀𝐷^𝛿_𝑥^𝑅₀𝐷^𝜃_𝑥𝑣(𝑥),∀𝑣∈𝐻^𝜃+𝛿(0,𝑋). (17)
(2) Let 0<𝜃<1. Then, one has ^𝑅₀𝐷^𝜃_𝑥𝑤,𝑣_{𝐿2(0,𝑋)}=𝑤,^𝑅_𝑥𝐷^𝜃_𝑋𝑣_{𝐿2(0,𝑋)},∀𝑤∈𝐻^𝜃(0,𝑋),𝑣∈𝐶^∞₀(0,𝑋). (18)

Lemma 3 (see [18]). Let 0<𝜃<2,𝜃≠1. Then, for any 𝑤,𝑣∈𝐻₀^𝜃/2(0,𝑋), then ^𝑅₀𝐷^𝜃_𝑥𝑤,𝑣_{𝐿2(0,𝑋)}=^𝑅₀𝐷_𝑥^𝜃/2𝑤,^𝑅_𝑡𝐷_𝑋^𝜃/2𝑣_{𝐿2(0,𝑋)}. (19)

Following from the definitions and lemmas above, we define the space 𝐵^{𝛼/2,𝛾/2}(𝐼×Ω)∶=𝐻^𝛼/2𝐼,𝐿₂(Ω)∩𝐻^𝛼1/2𝐼,𝐿₂(Ω)∩⋯∩𝐻^𝛼𝑠/2𝐼,𝐿₂(Ω)∩𝐿₂𝐼,𝐻₀^𝛽/2(Ω)∩𝐿₂𝐼,𝐻₀^𝛾/2(Ω), (20) when 0<𝛼_𝑠<⋯<𝛼₁<𝛼<1,0<𝛽<1 and 1<𝛾<2.

Then, one obtains 𝐵^{𝛼/2,𝛾/2}(𝐼×Ω)=𝐻^𝛼/2𝐼,𝐿₂(Ω)∩𝐿₂𝐼,𝐻₀^𝛾/2(Ω). (21) Here, 𝐵^{𝛼/2,𝛾/2}(𝐼×Ω) is a Banach space with respect to the following norm: ‖𝑣‖_𝐵^{𝛼/2,𝛾/2}_(𝐼×Ω)∶=‖𝑣‖²_{𝐻𝛼/2𝐼,𝐿2(Ω)}+‖𝑣‖²_{𝐿2𝐼,𝐻0𝛾/2(Ω)}^1/2,∀𝑣∈𝐵^{𝛼/2,𝛾/2}(𝐼×Ω), (22) where 𝐻^𝛼/2𝐼,𝐿₂(Ω)∶=𝑣;‖𝑣(𝑡,⋅)‖_𝐿2(Ω)∈𝐻^𝛼/2(𝐼), (23) endowed with the norm ‖𝑣‖_𝐻^𝛼/2_{(𝐼,𝐿2(Ω))}∶=‖𝑣(𝑡,⋅)‖_{𝐿2(Ω)𝐻𝛼/2(𝐼)}. (24)

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], ^𝐶₀𝐷^𝜃_𝑡𝑓(𝑡)=1Γ(−𝜃)^𝑡₀(𝑡−𝑠)^−𝜃−1𝑓(𝑠)−𝑇_𝑛−1𝑠)=𝑓;0(𝑑𝑠_𝑅0𝐷^𝜃_𝑡𝑓(𝑡)−^𝑅₀𝐷^𝜃_𝑡𝑇_𝑛−1𝑓;0(𝑡), (25) where 𝑇_𝑛−1[𝑓;0](𝑡) denotes the Taylor polynomial for 𝑓 of order 𝑛−1, centered at 0, 𝑇_𝑛−1𝑓;0(𝑡)=^𝑛_𝑘=0𝑡^𝑘𝑑𝑘!𝑑𝑡^𝑘𝑘𝑓(0), (26) where 𝑓(𝑡)∈𝐶^𝑛−1[0,𝑇].

Based on (25), we make 𝑓()()+𝑢𝑡,𝑥=𝑓𝑡,𝑥⁰(𝑥)𝑡^−𝛼Γ()+1−𝛼^𝑠_𝑖=1𝑎_𝑖𝑢₀(𝑥)𝑡^−𝛼𝑖Γ1−𝛼_𝑖. (27) Therefore, that led to the following weak formulation of (5). Let 𝐵^{𝛼/2,𝛾/2}(𝐼×Ω)^′ be the dual space of 𝐵^{𝛼/2,𝛾/2}(𝐼×Ω). For 𝑓∈𝐵^{𝛼/2,𝛾/2}(𝐼×Ω)^′, we find 𝑢(𝑡,𝑥)∈𝐵^{𝛼/2,𝛾/2}(𝐼×Ω) such that 𝒜(𝑢,𝑣)=ℱ(𝑣),𝑣∈𝐵^{𝛼/2,𝛾/2}(𝐼×Ω), (28) where the bilinear form 𝒜(⋅,⋅) is, based on Lemmas 2 and 3, 𝒜(𝑢,𝑣)∶=^𝑅₀𝐷_𝑡^𝛼/2𝑢,^𝑅_𝑡𝐷_𝑇^𝛼/2𝑣_𝐿2(𝒬)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝑢,^𝑅_𝑡𝐷^𝛼𝑖_𝑇^/2𝑣_𝐿2(𝒬)+𝑘_𝛽𝐶_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝑢,^𝑅_𝑥𝐷_𝑋^𝛽/2𝑣_𝐿2(𝒬)+^𝑅_𝑥𝐷_𝑋^𝛽/2𝑢,^𝑅₀𝐷_𝑥^𝛽/2𝑣_𝐿2(𝒬)+𝑘_𝛾𝐶_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝑢,^𝑅_𝑥𝐷_𝑋^𝛾/2𝑣_𝐿2(𝒬)+^𝑅_𝑥𝐷_𝑋^𝛾/2𝑢,^𝑅₀𝐷_𝑥^𝛾/2𝑣_𝐿2(𝒬), (29) and the functional ℱ(⋅) is given by ℱ(𝑣)∶=𝑓,𝑣_{𝐿2(𝒬)′}. (30)

Lemma 4 (see [17]). For real 𝜃>0, 𝑣∈𝐶^∞₀(𝐑), then ^𝑅_−∞𝐷^𝜃_𝑡𝑣(𝑡),^𝑅_𝑡𝐷^𝜃_∞𝑣(𝑡)_𝐿2(𝐑)=cos(𝜋𝜃)^𝑅_−∞𝐷^𝜃_𝑡𝑣(𝑡)²_𝐿2(𝐑),_𝑅−∞𝐷^𝜃_𝑡𝑣(𝑡),^𝑅_𝑡𝐷^𝜃_∞𝑣(𝑡)_𝐿2(𝐑)=cos(𝜋𝜃)^𝑅_𝑡𝐷^𝜃_∞𝑣(𝑡)²_𝐿2(𝐑). (31)

Based on Lemma 4, we can prove the following existence and uniqueness theorem. During this paper, we use the expression 𝐴≲𝐵(𝐴≳𝐵) to mean that there exists a positive real number 𝑐 such that 𝐴≤𝑐𝐵 (𝐴≥𝑐𝐵). At the same time, we denote 𝐴≅𝐵 to mean that 𝐴≲𝐵≲𝐴, which means there exist positive real numbers 𝑐₁,𝑐₂ such that 𝐴≤𝑐₁𝐵 and 𝐵≤𝑐₂𝐴 (i.e., (1/𝑐₁)𝐴≤𝐵≤𝑐₂𝐴).

Theorem 5. Assume that 0<𝛼₁<⋯<𝛼_𝑠<𝛼<1,0<𝛽<1, 1<𝛾<2 and 𝑓∈𝐵^{𝛼/2,𝛾/2}(𝐼×Ω)^′. Then, system (28) has a unique solution in 𝐵^{𝛼/2,𝛾/2}(𝐼×Ω). Furthermore, ‖𝑢‖_𝐵^{𝛼/2,𝛾/2}_(𝐼×Ω)≲𝑓_{𝐵𝛼/2,𝛾/2(𝐼×Ω)′}. (32)

Proof. The existence and uniqueness of the solution is guaranteed by the well-known Lax-Milgram theorem. First, from the equivalence of ^𝑙𝐻^𝛼₀(𝐼×Ω), ^𝑟𝐻^𝛼₀(𝐼×Ω), and 𝐻^𝛼₀(𝐼×Ω), for all 𝑢,𝑣∈𝐵^{𝛼/2,𝛾/2}(𝐼×Ω), it follows that |𝒜(𝑢,𝑣)|=^𝑅₀𝐷_𝑡^𝛼/2𝑢,^𝑅_𝑡𝐷_𝑇^𝛼/2𝑣_𝐿2(𝒬)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝑢,^𝑅_𝑡𝐷^𝛼𝑖_𝑇^/2𝑣_𝐿2(𝒬)+𝑘_𝛽𝐶_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝑢,^𝑅_𝑥𝐷_𝑋^𝛽/2𝑣_𝐿2(𝒬)+^𝑅_𝑥𝐷_𝑋^𝛽/2𝑢,^𝑅₀𝐷_𝑥^𝛽/2𝑣_𝐿2(𝒬)+𝑘_𝛾𝐶_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝑢,^𝑅_𝑥𝐷_𝑋^𝛾/2𝑣_𝐿2(𝒬)+^𝑅_𝑥𝐷_𝑋^𝛾/2𝑢,^𝑅₀𝐷_𝑥^𝛾/2𝑣_𝐿2(𝒬)≲^𝑅₀𝐷_𝑡^𝛼/2𝑢_𝐿2(𝒬)^𝑅₀𝐷_𝑡^𝛼/2𝑣_𝐿2(𝒬)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝑢_𝐿2(𝒬)^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝑣_𝐿2(𝒬)+𝑘_𝛽𝐶_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝑢_𝐿2(𝒬)^𝑅_𝑥𝐷_𝑋^𝛽/2𝑣_𝐿2(𝒬)+^𝑅_𝑥𝐷_𝑋^𝛽/2𝑢_𝐿2(𝒬)^𝑅₀𝐷_𝑥^𝛽/2𝑣_𝐿2(𝒬)+𝑘_𝛾𝐶_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝑢_𝐿2(𝒬)^𝑅_𝑥𝐷_𝑋^𝛾/2𝑣_𝐿2(𝒬)+^𝑅_𝑥𝐷_𝑋^𝛾/2𝑢_𝐿2(𝒬)^𝑅₀𝐷_𝑥^𝛾/2𝑣_𝐿2(𝒬)≲‖𝑢‖_𝐻^𝛼/2_{(𝐼,𝐿2(Ω))}‖𝑣‖_𝐻^𝛼/2_{(𝐼,𝐿2(Ω))}+^𝑠_𝑖=1𝑎_𝑖‖𝑢‖_𝐻^𝛼_𝑖^/2_{(𝐼,𝐿2(Ω))}‖𝑣‖_𝐻^𝛼_𝑖^/2_{(𝐼,𝐿2(Ω))}+𝑘_𝛽‖𝑢‖_𝐿²_(𝐼,𝐻^𝛽/2_(Ω))‖𝑣‖_𝐿²_(𝐼,𝐻^𝛽/2_(Ω))+𝑘_𝛾‖𝑢‖_𝐿²_(𝐼,𝐻^𝛾/2_(Ω))‖𝑣‖_𝐿²_(𝐼,𝐻^𝛾/2_(Ω))≲‖𝑢‖_𝐵^{𝛼/2,𝛾/2}_(𝐼×Ω)‖𝑣‖_𝐵^{𝛼/2,𝛾/2}_(𝐼×Ω),|ℱ(𝑣)|≲𝑓_{𝐵𝛼/2,𝛾/2(𝐼×Ω)′}‖𝑣‖_𝐵^{𝛼/2,𝛾/2}_(𝐼×Ω). (33)(34)
This implies the continuity of the bilinear form 𝒜(⋅,⋅) and the right-hand side function ℱ(𝑣).
We next prove the coercivity of the bilinear operator 𝒜(⋅,⋅). Note that ^𝑅₀𝐷_𝑡^𝛼/2𝜑,^𝑅_𝑡𝐷_𝑇^𝛼/2𝜑_𝐿2(𝒬)=^𝑅_−∞𝐷_𝑡^∼𝛼/2𝜑,^𝑅_𝑡𝐷_∞^∼𝛼/2𝜑_𝐿2(𝐑)=cos𝜋𝛼2⋅^𝑅_−∞𝐷_𝑡^∼𝛼/2𝜑_𝐿2(𝐑)=cos𝜋𝛼2⋅^𝑅₀𝐷_𝑡^𝛼/2𝜑_𝐿2(𝒬), (35) for all 𝜑∈𝐶^∞₀(𝒬), where 𝜑^∼ is the extension of 𝜑 by zero outside of (0,𝑇). Thus, we find that (^𝑅₀𝐷_𝑡^𝛼/2𝑣,^𝑅_𝑡𝐷_𝑇^𝛼/2𝑣)_𝐿2(𝒬) is nonnegative for 𝑣∈𝐻^𝛼/2(0,𝑇),0<𝛼<1 since cos(𝜋𝛼/2) is nonnegative for 0<𝛼<1. That is the same for fractional operators with 𝛼_𝑖,𝑖=1,…,𝑠.
From the above analysis, we have 𝒜(𝑣,𝑣)=^𝑅₀𝐷_𝑡^𝛼/2𝑣,^𝑅_𝑡𝐷_𝑇^𝛼/2𝑣_𝐿2(𝒬)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝑣,^𝑅_𝑡𝐷^𝛼𝑖_𝑇^/2𝑣_𝐿2(𝒬)+2𝑘_𝛽𝐶_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝑣,^𝑅_𝑥𝐷_𝑋^𝛽/2𝑣_𝐿2(𝒬)+2𝑘_𝛾𝐶_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝑣,^𝑅_𝑥𝐷_𝑋^𝛾/2𝑣_𝐿2(𝒬)≅cos𝜋𝛼2^𝑅0𝐷_𝑡^𝛼/2𝑣,^𝑅₀𝐷_𝑡^𝛼/2𝑣_𝐿2(𝒬)+^𝑠_𝑖=1𝑎_𝑖cos𝜋𝛼_𝑖2^𝑅0𝐷^𝛼𝑖_𝑡^/2𝑣,^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝑣_𝐿2(𝒬)+𝑘_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝑣,^𝑅₀𝐷_𝑥^𝛽/2𝑣_𝐿2(𝒬)+𝑘_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝑣,^𝑅₀𝐷_𝑥^𝛾/2𝑣_𝐿2(𝒬)≳‖𝑣‖²_{𝐵𝛼/2,𝛾/2(𝐼×Ω)}. (36)
By using the well-known Lax-Milgram theorem, there exists a unique solution 𝑢∈𝐵^{𝛼/2,𝛾/2}(𝐼×Ω) such that (28) holds.
To prove the stability estimate (32), by using (34) and (36), we take 𝑣=𝑢 in (28) and obtain ‖𝑢‖²_{𝐵𝛼/2,𝛾/2(𝐼×Ω)}()(𝑢)≲≲𝒜𝑢,𝑢=ℱ𝑓_{𝐵𝛼/2,𝛾/2(𝐼×Ω)′}‖𝑢‖_𝐵^{𝛼/2,𝛾/2}_(𝐼×Ω), (37) which implies that ‖𝑢‖_𝐵^{𝛼/2,𝛾/2}_(𝐼×Ω)≲‖𝑓‖_{𝐵𝛼/2,𝛾/2(𝐼×Ω)′}.

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]). For 𝑛−1<𝜃<𝑛(𝑛=1,2,…) and a function 𝑢(𝑥) defined on an unbounded domain (−∞<𝑥<∞), the following equality holds: ^𝑅𝐷^𝜃_|𝑥|𝑢(𝑥)1=−2cos(𝜋𝜃/2)^𝑅−∞𝐷^𝜃_𝑥𝑢+^𝑅_𝑥𝐷^𝜃_∞𝑢=−(−Δ)^𝜃/2𝑢(𝑥), (38) where Δ is the Laplacian.

Definition 7 (see [26]). Suppose the Laplacian has a complete set of orthonormal eigenfunctions 𝜙_𝑛 corresponding to eigenvalues 𝜆²_𝑛 on a bounded region 𝒬, that is,−Δ𝜙_𝑛=𝜆²_𝑛𝜙_𝑛 on a bounded region 𝒬; 𝐵(𝜙)=0 on 𝜕𝒬 is one of the standard three homogeneous boundary conditions. Let 𝐅=𝑓=^∞_𝑛=1𝑐_𝑛𝜙_𝑛,𝑐_𝑛=𝑓,𝜙_𝑛,^∞_𝑛=1𝑐_𝑛²𝜆_𝑛^𝜃≤∞, (39) then for any 𝑓∈𝐅, (−Δ)^𝜃/2 is defined by (−Δ)^𝜃/2𝑓=^∞_𝑛=1𝑐_𝑛𝜆²_𝑛^𝜃/2𝜙_𝑛. (40)

Here, for any 𝑡∈𝐼, we make a continuous prolongation in space for function 𝑢(𝑡,𝑥) as follows: 𝑢^∼()=𝑢()[],()∪().𝑡,𝑥𝑡,𝑥,𝑥∈0,𝑋0,𝑥∈−∞,0𝑋,∞ (41) Let 𝑢(𝑡) denote the one-variable function of 𝑢^∼(𝑡,𝑥), and let 𝑓(𝑡) and 𝑢(0) denote the one-variable functions as 𝑓(𝑡,⋅) and 𝑢^∼(0,⋅), respectively. By the definition (12) of 𝑃(₀𝐷_𝑡) and Lemma 6, (5) is equivalent to ^𝐶₀𝐷^𝛼_𝑡𝑢(𝑡)+^𝑠_𝑖=1𝑎_𝑖^𝐶₀𝐷^𝛼𝑖_𝑡𝑢(𝑡)+𝑘_𝛽(−Δ)^𝛽/2𝑢(𝑡)+𝑘_𝛾(−Δ)^𝛾/2𝑢(𝑡)=𝑓(𝑡). (42)

Let 𝐴=−Δ, 𝐷(𝐴)=𝐻¹₀(Ω)∩𝐻²(Ω). Then, the system (5) and the initial condition (6) can be rewritten in the abstract form, for 0≤𝑡≤𝑇, ^𝐶₀𝐷^𝛼_𝑡𝑢(𝑡)+^𝑠_𝑖=1𝑎_𝑖^𝐶₀𝐷^𝛼𝑖_𝑡𝑢(𝑡)+𝑘_𝛽𝐴^𝛽/2𝑢(𝑡)+𝑘_𝛾𝐴^𝛾/2𝑢(𝑡)(𝑡),𝑢(0)=𝑓=𝜙₀. (43)(44) Based on (25), (43) can be rewritten as ^𝑅₀𝐷^𝛼_𝑡𝑢−𝜙₀(𝑡)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡𝑢−𝜙₀(𝑡)+𝑘_𝛽𝐴^𝛽/2𝑢(𝑡)+𝑘_𝛾𝐴^𝛾/2𝑢(𝑡)=𝑓(𝑡), (45) where 0<𝑡<𝑇,0<𝛼_𝑠<⋯<𝛼₁<𝛼<1.

Let 0=𝑡₀<𝑡₁<⋯<𝑡_𝑁=𝑇 be a partition of [0,𝑇]. Without loss of generality, for 0<𝛼<1 and fix 𝑡_𝑛,𝑛=1,2,…,𝑁, we have ^𝑅₀𝐷^𝛼_𝑡𝑢−𝜙₀𝑡_𝑗=𝑡_𝑗^−𝛼Γ(−𝛼)¹₀𝑔(𝜃)𝜃^−1−𝛼𝑑𝜃, (46) where 𝑔(𝜃)=𝑢(𝑡_𝑗−𝑡_𝑗𝜃)−𝜙₀. Here, the integral is a Hadamard finite-part integral, in [38, 39].

Now, for every 𝑛, we replace the integral by a first-degree compound quadrature formula with equispaced nodes 0,1/𝑛,2/𝑛,…,1 and obtain ¹₀𝑔(𝜃)𝜃^−1−𝛼𝑑𝜃=^𝑛_𝑘=0𝛼^(𝛼)_𝑘𝑛𝑔𝑘𝑛+𝑅_𝑛^(𝛼)𝑔,0<𝛼<1, (47) where the weights 𝛼^(𝛼)_𝑘𝑛 are 𝛼(1−𝛼)𝑛^−𝛼𝛼^(𝛼)_𝑘𝑛=−1,𝑘=0,2𝑘^1−𝛼−(𝑘−1)^1−𝛼−(𝑘+1)^1−𝛼,()𝑘for𝑘=1,2,…,𝑛−1,𝛼−1^−𝛼−(𝑘−1)^1−𝛼+𝑘^1−𝛼,for𝑘=𝑛, (48) and the remainder term 𝑅_𝑛^(𝛼)(𝑔) satisfies 𝑅_𝑛^(𝛼)𝑔≤𝛾_𝛼𝑛^𝛼−2sup_{0≤𝑡≤𝑇}𝑔^′′(𝑡), (49) where 𝛾_𝛼>0 is a constant.

Thus, we have ^𝑅₀𝐷^𝛼_𝑡𝑢−𝜙₀𝑡_𝑛=𝑡_𝑛^−𝛼Γ(−𝛼)^𝑛_𝑘=0𝛼^(𝛼)_𝑘𝑛𝑢𝑡_𝑛−𝑡_𝑘−𝜙₀+𝑅_𝑛^(𝛼)𝑔=Δ𝑡^𝑛−𝛼_𝑘=0𝜔^(𝛼)_𝑘𝑛𝑢𝑡_𝑛−𝑡_𝑘−𝜙₀+𝑡_𝑛^−𝛼Γ()𝑅−𝛼^𝑛(𝛼)𝑔, (50) where 𝜔^(𝛼)_𝑘𝑛=𝑛^−𝛼𝛼^(𝛼)_𝑘𝑛/Γ(−𝛼).

Let 𝑡=𝑡_𝑛, we can write (45) as Δ𝑡^𝑛−𝛼_𝑘=0𝜔^(𝛼)_𝑘𝑛𝑢𝑡_𝑛−𝑡_𝑘−𝜙₀+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{−𝛼𝑖𝑛}_𝑘=0𝜔^(𝛼𝑖)_𝑘𝑛𝑢𝑡_𝑛−𝑡_𝑘−𝜙₀+𝑘_𝛽𝐴^𝛽/2𝑢𝑡_𝑛+𝑘_𝛾𝐴^𝛾/2𝑢𝑡_𝑛=𝑓𝑡_𝑛−𝑡_𝑛^−𝛼Γ()𝑅−𝛼^𝑛(𝛼)𝑔−^𝑠_𝑖=1𝑎_𝑖𝑡^−𝛼𝑖_𝑛Γ−𝛼_𝑖𝑅^(𝛼𝑖)_𝑛𝑔, (51) for 𝑛=1,2,3,…, where 𝜔^(𝛼𝑖)_𝑘𝑛=𝑛^−𝛼𝑖𝛼^(𝛼𝑖)_𝑘𝑛/Γ(−𝛼_𝑖).

Denote 𝑈^𝑛 as the approximation of 𝑢(𝑡_𝑛). We can define the following time stepping method: Δ𝑡^𝑛−𝛼_𝑘=0𝜔^(𝛼)_𝑘𝑛𝑈^𝑛−𝑘−𝑈⁰+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{−𝛼𝑖𝑛}_𝑘=0𝜔^(𝛼𝑖)_𝑘𝑛𝑈^𝑛−𝑘−𝑈⁰+𝑘_𝛽𝐴^𝛽/2𝑈^𝑛+𝑘_𝛾𝐴^𝛾/2𝑈^𝑛=𝑓_𝑛, (52) where 𝑓_𝑛=𝑓(𝑡_𝑛).

Lemma 8 (see [38]). For 0<𝛼<1, let the sequence {𝑑_𝑗}𝑗=1,2,… be given by 𝑑₁=1 and 𝑑_𝑗=1+𝛼(1−𝛼)𝑗^−𝛼∑^𝑗−1_𝑘=1𝛼_𝑘𝑗𝑑_𝑗−𝑘,𝑗=2,3,…. Then, 1≤𝑑_𝑗≤(sin(𝜋𝛼))/(𝜋𝛼(1−𝛼))𝑗^𝛼,𝑗=1,2,3,….

Let 𝑒^𝑛=𝑈^𝑛−𝑢(𝑡_𝑛) denote the error in 𝑡_𝑛. Let ‖⋅‖ be any given norm. Then, we have the following error estimate.

Theorem 9. Let 𝑢(𝑡_𝑛) and 𝑈^𝑛 be the solutions of (43) and (52), respectively, and Δ𝑡 is the time step size. Then, one has max_{𝑛=0,1,…,𝑁}𝑈^𝑛−𝑢𝑡_𝑛=𝑂Δ𝑡^2−𝛼+^𝑠_𝑖=0𝑎_𝑖Δ𝑡^2−𝛼𝑖. (53)

Proof. Subtracting (52) from (51), we obtain the error equation Δ𝑡^𝑛−𝛼_𝑘=0𝜔^(𝛼)_𝑘𝑛𝑒^𝑛−𝑘−𝑒⁰+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{−𝛼𝑖𝑛}_𝑘=0𝜔^(𝛼𝑖)_𝑘𝑛𝑒^𝑛−𝑘−𝑒⁰+𝑘_𝛽𝐴^𝛽/2𝑒^𝑛+𝑘_𝛾𝐴^𝛾/2𝑒^𝑛𝑡=−^𝑛−𝛼Γ()𝑅−𝛼^𝑛(𝛼)𝑔−^𝑠_𝑖=1𝑎_𝑖𝑡^−𝛼𝑖_𝑛Γ−𝛼_𝑖𝑅^(𝛼𝑖)_𝑛𝑔. (54) Note that 𝑒⁰=𝜙₀−𝑈⁰=0 and denote 𝑒^𝑛=Δ𝑡^−𝛼𝜔^(𝛼)_0𝑛+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖𝜔^(𝛼𝑖)_0𝑛+𝑘_𝛽𝐴^𝛽/2+𝑘_𝛾𝐴^𝛾/2−1×Δ𝑡^𝑛−𝛼_𝑘=1𝜔^(𝛼)_𝑘𝑛𝑒^𝑛−𝑘+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{−𝛼𝑖𝑛}_𝑘=1𝜔^(𝛼𝑖)_𝑘𝑛𝑒^𝑛−𝑘−𝑡_𝑛^−𝛼Γ()𝑅−𝛼^𝑛(𝛼)𝑔−^𝑠_𝑖=1𝑎_𝑖𝑡^−𝛼𝑖_𝑛Γ−𝛼_𝑖𝑅^(𝛼𝑖)_𝑛𝑔=𝛼^(𝛼)_0𝑛+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{𝛼−𝛼𝑖}Γ(−𝛼)Γ−𝛼_𝑖𝛼^(𝛼𝑖)_0𝑛+𝑘_𝛽𝐴^𝛽/2𝑡^𝛼_𝑛Γ(−𝛼)+𝑘_𝛾𝐴^𝛾/2𝑡^𝛼_𝑛Γ(−𝛼)⁻¹×^𝑛_𝑘=1𝛼^(𝛼)_𝑘𝑛𝑒^𝑛−𝑘+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{𝛼−𝛼𝑖}Γ(−𝛼)Γ−𝛼_𝑖×^𝑛_𝑘=1𝛼^(𝛼𝑖)_𝑘𝑛𝑒^𝑛−𝑘−𝑅_𝑛^(𝛼)𝑔−^𝑠_𝑖=1𝑎_𝑖Γ(−𝛼)Γ−𝛼_𝑖𝑡^{𝛼−𝛼𝑖}_𝑛𝑅^(𝛼𝑖)_𝑛𝑔. (55) Thus, we have 𝑒^𝑛≤𝛼^(𝛼)_0𝑛+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{𝛼−𝛼𝑖}Γ(−𝛼)Γ−𝛼_𝑖𝛼^(𝛼𝑖)_0𝑛+𝑘_𝛽𝐴^𝛽/2𝑡^𝛼_𝑛Γ(−𝛼)+𝑘_𝛾𝐴^𝛾/2𝑡^𝛼_𝑛Γ(−𝛼)⁻¹×^𝑛_𝑘=1𝛼^(𝛼)_𝑘𝑛𝑒^𝑛−𝑘+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{𝛼−𝛼𝑖}Γ(−𝛼)Γ−𝛼_𝑖^𝑗_𝑘=1𝛼^(𝛼𝑖)_𝑘𝑛𝑒^𝑛−𝑘+𝑅_𝑛^(𝛼)𝑔+^𝑠_𝑖=1𝑎_𝑖Γ(−𝛼)Γ−𝛼_𝑖𝑡^{𝛼−𝛼𝑖}_𝑛𝑅^(𝛼𝑖)_𝑛𝑔. (56)
From Definition 7, note that 𝐴 is a positive definite elliptic operator with all of eigenvalues 𝜆>0. Since 𝛼^(𝛼𝑖)_0𝑛<0,𝑖=1,…,𝑠,𝛼^(𝛼)_0𝑛<0 and Γ(−𝛼)<0 when 0<𝛼<1, we have 𝛼^(𝛼)_0𝑛+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{𝛼−𝛼𝑖}Γ(−𝛼)Γ−𝛼_𝑖𝛼^(𝛼𝑖)_0𝑛+𝑘_𝛽𝐴^𝛽/2𝑡^𝛼_𝑛Γ(−𝛼)+𝑘_𝛾𝐴^𝛾/2𝑡^𝛼_𝑛Γ(−𝛼)⁻¹=sup_𝜆>0𝛼^(𝛼)_0𝑛+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{𝛼−𝛼𝑖}Γ(−𝛼)Γ−𝛼_𝑖𝛼^(𝛼𝑖)_0𝑛+𝜆^𝛽/2𝑡^𝛼_𝑛Γ(−𝛼)+𝜆^𝛾/2𝑡^𝛼_𝑛Γ(−𝛼)⁻¹<−𝛼^(𝛼)_0𝑛−^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{𝛼−𝛼𝑖}Γ(−𝛼)Γ(−𝛼_𝑖)𝛼^(𝛼𝑖)_0𝑛⁻¹. (57) Hence, 𝑒^𝑛≤−𝛼^(𝛼)_0𝑛^𝑛−1_𝑘=1𝛼^(𝛼)_𝑘𝑛𝑒^𝑛−𝑘+^𝑠_𝑖=1−𝛼^(𝛼𝑖)_0𝑛^𝑛−1_𝑘=1𝛼^(𝛼𝑖)_𝑘𝑛𝑒^𝑛−𝑘+−𝛼^(𝛼)_0𝑛⁻¹𝑅_𝑛^(𝛼)𝑔+^𝑠_𝑖=1−𝛼^(𝛼𝑖)_0𝑛⁻¹𝑅^(𝛼𝑖)_𝑛𝑔≤𝛼(1−𝛼)𝑛^𝑛−𝛼_𝑘=1𝛼^(𝛼)_𝑘𝑛𝑒^𝑛−𝑘+^𝑠_𝑖=1𝛼_𝑖1−𝛼_𝑖𝑛^{−𝛼𝑖𝑛}_𝑘=1𝛼^(𝛼𝑖)_𝑘𝑛𝑒^𝑛−𝑘+𝛼(1−𝛼)𝛾_𝛼𝑁⁻²sup_{0≤𝑡≤𝑇}𝑢^′′+^𝑠_𝑖=1𝛼_𝑖1−𝛼_𝑖𝛾_𝛼𝑖𝑁⁻²sup_{0≤𝑡≤𝑇}𝑢^′′. (58) Note that 𝑒⁰=𝜙₀−𝑈⁰=0. Denote 𝑑₁=1 and 𝑑_𝑛=1+𝛼(1−𝛼)𝑛^{−𝛼𝑛−1}_𝑘=1𝛼^(𝛼)_𝑘𝑛𝑑_𝑛−𝑘,𝑑𝑛=2,3,…,𝑁,_𝑖𝑛=1+𝛼_𝑖1−𝛼_𝑖𝑛^{−𝛼𝑖𝑗−1}_𝑘=1𝛼^(𝛼𝑖)_𝑘𝑛𝑑_𝑛−𝑘,𝑛=2,3,…,𝑁,𝑖=1,2,…,𝑠. (59)
Then, by induction and Lemma 8, we have 𝑒^𝑛≲𝛼(1−𝛼)𝑁⁻²sup_{0≤𝑡≤𝑇}𝑢^′′(𝑡)⋅𝑑_𝑛+^𝑠_𝑖=1𝑎_𝑖𝛼_𝑖1−𝛼_𝑖𝑁⁻²sup_{0≤𝑡≤𝑇}𝑢^′′(𝑡)⋅𝑑^𝑖_𝑛≲𝑁⁻²sup_{0≤𝑡≤𝑇}𝑢^′′(𝑡())sin𝜋𝛼𝜋𝑛^𝛼+^𝑠_𝑖=1𝑎_𝑖𝑁⁻²sup_{0≤𝑡≤𝑇}𝑢^′′(𝑡)sin𝜋𝛼_𝑖𝜋𝑛^𝛼𝑖=𝑂Δ𝑡^2−𝛼+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^2−𝛼𝑖. (60)

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,𝑋] 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 Section 3 and (25), we need to find 𝑢(𝑡,⋅)∈𝐻₀^𝛾/2(Ω) such that ^𝑅₀𝐷^𝛼_𝑡𝑢(𝑡,𝑥),𝑣_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡𝑢(𝑡,𝑥),𝑣_𝐿2(Ω)+𝑘_𝛽𝐶_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝑢(𝑡,𝑥),^𝑅_𝑥𝐷_𝑋^𝛽/2𝑣_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝛽/2𝑢(𝑡,𝑥),^𝑅₀𝐷_𝑥^𝛽/2𝑣_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝑢(𝑡,𝑥),^𝑅_𝑥𝐷_𝑋^𝛾/2𝑣_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝛾/2𝑢(𝑡,𝑥),^𝑅₀𝐷_𝑥^𝛾/2𝑣_𝐿2(Ω)=𝑓(𝑡,𝑥),𝑣_𝐿2(Ω),∀𝑣∈𝐻₀^𝛾/2(Ω). (61)

Let ℎ denote the maximal length of intervals in Ω, and let 𝑟 be any nonnegative integer. We denote the norm in 𝐻^𝑟(Ω) by ‖⋅‖_𝐻^𝑟_(Ω). Let 𝑆_ℎ⊂𝐻₀^𝛾/2 be a family of finite element space with the accuracy of order 𝑟≥2;thatis,𝑆_ℎ consists of continuous functions on the closure Ω of Ω which are polynomials of degree at most 𝑟−1 in each interval and which vanish outside Ω_ℎ such that for small ℎ, 𝑣∈𝐻^𝜂(Ω)∩𝐻₀^𝛾/2(Ω), inf_{𝜒∈𝑆ℎ}𝑣−𝜒_𝐿2(Ω)+ℎ∇𝑣−𝜒_𝐿2(Ω)≲ℎ^𝜂‖𝑣‖_𝐻^𝜂_(Ω),1≤𝜂≤𝑟. (62)

The semidiscrete problem of (5) is to find the approximate solution 𝑢_ℎ(𝑡)=𝑢_ℎ(𝑡,⋅)∈𝑆_ℎ for each 𝑡 such that ^𝑅₀𝐷^𝛼_𝑡𝑢_ℎ(𝑡),𝜒_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡𝑢_ℎ(𝑡),𝜒_𝐿2(Ω)+𝑘_𝛽𝐶_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝑢_ℎ(𝑡),^𝑅_𝑥𝐷_𝑋^𝛽/2𝜒_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝛽/2𝑢_ℎ(𝑡),^𝑅₀𝐷_𝑥^𝛽/2𝜒_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝑢_ℎ(𝑡),^𝑅_𝑥𝐷_𝑋^𝛾/2𝜒_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝛾/2𝑢_ℎ(𝑡),^𝑅₀𝐷_𝑥^𝛾/2𝜒_𝐿2(Ω)=𝑓(𝑡),𝜒_𝐿2(Ω),∀𝜒∈𝑆_ℎ. (63) Let 𝑅_ℎ∶𝐻^𝛾/2(Ω)→𝑆_ℎ be the elliptic projection defined by 𝑘_𝛽𝐶_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝑅_ℎ𝑢,^𝑅_𝑥𝐷_𝑋^𝛽/2𝜒_𝐿2(Ω)+𝑘_𝛽𝐶_𝛽^𝑅_𝑥𝐷_𝑋^𝛽/2𝑅_ℎ𝑢,^𝑅₀𝐷_𝑥^𝛽/2𝜒_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝑅_ℎ𝑢,^𝑅_𝑥𝐷_𝑋^𝛾/2𝜒_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅_𝑥𝐷_𝑋^𝛾/2𝑅_ℎ𝑢,^𝑅₀𝐷_𝑥^𝛾/2𝜒_𝐿2(Ω)=𝑘_𝛽𝐶_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝑢,^𝑅_𝑥𝐷_𝑋^𝛽/2𝜒_𝐿2(Ω)+𝑘_𝛽𝐶_𝛽^𝑅_𝑥𝐷_𝑋^𝛽/2𝑢,^𝑅₀𝐷_𝑥^𝛽/2𝜒_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝑢,^𝑅_𝑥𝐷_𝑋^𝛾/2𝜒_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅_𝑥𝐷_𝑋^𝛾/2𝑢,^𝑅₀𝐷_𝑥^𝛾/2𝜒_𝐿2(Ω),∀𝜒∈𝑆_ℎ. (64)

Lemma 10. Assume that (62) holds. Then, for 𝑅_ℎ defined by (64) and any 𝑣∈𝐻^𝜂(Ω)∩𝐻₀^𝜃/2(Ω), one has ^𝑅₀𝐷_𝑥^𝜃/2(𝑅_ℎ𝑣−𝑣)_𝐿2(Ω)≲ℎ^{𝜂−𝜃/2}‖𝑣‖_𝐻^𝜂_(Ω). (65)

Proof. Let 𝐼_ℎ be the projection operator from 𝐻^𝜂(Ω)∩𝐻₀^𝜃/2(Ω) to 𝑆_ℎ. From the definition of 𝐿₂-norm, we obtain ^𝑅₀𝐷_𝑥^𝜃/2(𝑅_ℎ𝑣−𝑣),^𝑅_𝑥𝐷_𝑋^𝜃/2(𝑅_ℎ𝑣−𝑣)_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝜃/2𝑅_ℎ𝑣−𝑣,^𝑅₀𝐷_𝑥^𝜃/2(𝑅_ℎ𝑣−𝑣)_𝐿2(Ω)=^𝑅₀𝐷_𝑥^𝜃/2(𝑅_ℎ𝑣−𝑣),^𝑅_𝑥𝐷_𝑋^𝜃/2(𝑅_ℎ𝑣−𝐼_ℎ𝑣)_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝜃/2(𝑅_ℎ𝑣−𝑣),^𝑅₀𝐷_𝑥^𝜃/2(𝑅_ℎ𝑣−𝐼_ℎ𝑣)_𝐿2(Ω)+^𝑅₀𝐷_𝑥^𝜃/2(𝑅_ℎ𝑣−𝑣),^𝑅_𝑥𝐷_𝑋^𝜃/2(𝐼_ℎ𝑣−𝑣)_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝜃/2(𝑅_ℎ𝑣−𝑣),^𝑅₀𝐷_𝑥^𝜃/2(𝐼_ℎ𝑣−𝑣)_𝐿2(Ω). (66) Let 𝜒=𝑅_ℎ𝑣−𝐼_ℎ𝑣∈𝑆_ℎ. From (64), we obtain ^𝑅₀𝐷_𝑥^𝜃/2(𝑅_ℎ𝑣−𝑣),^𝑅_𝑥𝐷_𝑋^𝜃/2(𝑅_ℎ𝑣−𝑣)_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝜃/2(𝑅_ℎ𝑣−𝑣),^𝑅₀𝐷_𝑥^𝜃/2(𝑅_ℎ𝑣−𝑣)_𝐿2(Ω)≤^𝑅₀𝐷_𝑥^𝜃/2(𝑅_ℎ𝑣−𝑣)_𝐿2(Ω)^𝑅_𝑥𝐷_𝑋^𝜃/2(𝐼_ℎ𝑣−𝑣)_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝜃/2(𝑅_ℎ𝑣−𝑣)_𝐿2(Ω)^𝑅₀𝐷_𝑥^𝜃/2(𝐼_ℎ𝑣−𝑣)_𝐿2(Ω)≤^𝑅₀𝐷_𝑥^𝜃/2(𝑅_ℎ𝑣−𝑣)_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝜃/2(𝑅_ℎ𝑣−𝑣)_𝐿2(Ω)×^𝑅_𝑥𝐷_𝑋^𝜃/2(𝐼_ℎ𝑣−𝑣)_𝐿2(Ω)+^𝑅₀𝐷_𝑥^𝜃/2(𝐼_ℎ𝑣−𝑣)_𝐿2(Ω). (67)
Note that ^𝑅₀𝐷_𝑥^𝜃/2(𝑅_ℎ𝑣−𝑣),^𝑅_𝑥𝐷_𝑋^𝜃/2(𝑅_ℎ𝑣−𝑣)_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝜃/2(𝑅_ℎ𝑣−𝑣),^𝑅₀𝐷_𝑥^𝜃/2(𝑅_ℎ𝑣−𝑣)_𝐿2(Ω)≳^𝑅₀𝐷_𝑥^𝜃/2𝑅_ℎ𝑣−𝑣_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝜃/2𝑅_ℎ𝑣−𝑣_𝐿2(Ω)². (68) Thus, combining (67) and (68), we obtain ^𝑅₀𝐷_𝑥^𝜃/2(𝑅_ℎ𝑣−𝑣)_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝜃/2(𝑅_ℎ𝑣−𝑣)_𝐿2(Ω)≲^𝑅₀𝐷_𝑥^𝜃/2(𝐼_ℎ𝑣−𝑣)_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝜃/2(𝐼_ℎ𝑣−𝑣)_𝐿2(Ω). (69) By (62), we see that 𝑣−𝐼_ℎ𝑣_𝐿2(Ω)+ℎ∇(𝑣−𝐼_ℎ𝑣)_𝐿2(Ω)≲ℎ^𝜂|𝑣|_𝐻^𝜂_(Ω). (70) By interpolation properties, we obtain ^𝑅₀𝐷_𝑥^𝜃/2(𝑣−𝐼_ℎ𝑣)_𝐿2(Ω)≤(𝑣−𝐼_ℎ𝑣)_{𝐻𝜃/2(Ω)}≲ℎ^{𝜂−𝜃/2}|𝑣|_𝐻^𝜂/2_(Ω). (71) Similarly, we have ^𝑅_𝑥𝐷_𝑋^𝜃/2(𝑣−𝐼_ℎ𝑣)_𝐿2(Ω)≤(𝑣−𝐼_ℎ𝑣)_{𝐻𝜃/2(Ω)}≲ℎ^{𝜂−𝜃/2}|𝑣|_𝐻^𝜂/2_(Ω). (72) Combining (71) and (72), we can obtain (65).

Theorem 11. For 0<𝛼_𝑠<⋯<𝛼₁<𝛼<1, let 𝑢_ℎ and 𝑢 be the solutions of (63) and (61), respectively. Then, it holds 𝑢(𝑡,𝑥)−𝑢_ℎ(𝑡,𝑥)_{𝐻𝛾/2(Ω)}=𝑂Δ𝑥^{𝜂−𝛾/2}. (73)

Proof. We write 𝑢_ℎ−𝑢=𝜗+𝜌, (74) where 𝜗=𝑢_ℎ−𝑅_ℎ𝑢, 𝜌=𝑅_ℎ𝑢−𝑢. The second term is easily bounded by Lemma 10 and has the obvious estimate 𝜌(𝑥)²_{𝐻𝛾/2(Ω)}≲ℎ^2𝜂−𝛾‖𝑢‖²_𝐻^𝜂_(Ω). (75)
In order to estimate 𝜗, for all 𝜒∈𝑆_ℎ, we get ^𝑅₀𝐷^𝛼_𝑡𝜗(𝑡),𝜒_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡𝜗(𝑡),𝜒_𝐿2(Ω)+𝑘_𝛽𝐶_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝜗(𝑡),^𝑅_𝑥𝐷_𝑋^𝛽/2𝜒_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝛽/2𝜗(𝑡),^𝑅₀𝐷_𝑥^𝛽/2𝜒_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝜗(𝑡),^𝑅_𝑥𝐷_𝑋^𝛾/2𝜒_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝛾/2𝜗(𝑡),^𝑅₀𝐷_𝑥^𝛾/2𝜒_𝐿2(Ω)=^𝑅₀𝐷^𝛼_𝑡𝑢_ℎ(𝑡)−𝑅_ℎ𝑢(𝑡),𝜒_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡𝑢_ℎ(𝑡)−𝑅_ℎ𝑢(𝑡),𝜒_𝐿2(Ω)+𝑘_𝛽𝐶_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝑢_ℎ(𝑡)−𝑅_ℎ𝑢(𝑡),^𝑅_𝑥𝐷_𝑋^𝛽/2𝜒_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝛽/2𝑢_ℎ(𝑡)−𝑅_ℎ𝑢(𝑡),^𝑅₀𝐷_𝑥^𝛽/2𝜒_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝑢_ℎ(𝑡)−𝑅_ℎ𝑢(𝑡),^𝑅_𝑥𝐷_𝑋^𝛾/2𝜒_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝛾/2𝑢_ℎ(𝑡)−𝑅_ℎ𝑢(𝑡),^𝑅₀𝐷_𝑥^𝛾/2𝜒_𝐿2(Ω)=^𝑅₀𝐷^𝛼_𝑡𝑢_ℎ(𝑡),𝜒_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡𝑢_ℎ(𝑡),𝜒_𝐿2(Ω)+𝑘_𝛽𝐶_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝑢_ℎ(𝑡),^𝑅_𝑥𝐷_𝑋^𝛽/2𝜒_𝐿2(Ω)+𝑘_𝛽𝐶_𝛽^𝑅_𝑥𝐷_𝑋^𝛽/2𝑢_ℎ(𝑡),^𝑅₀𝐷_𝑥^𝛽/2𝜒_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝑢_ℎ(𝑡),^𝑅_𝑥𝐷_𝑋^𝛾/2𝜒_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅_𝑥𝐷_𝑋^𝛾/2𝑢_ℎ(𝑡),^𝑅₀𝐷_𝑥^𝛾/2𝜒_𝐿2(Ω)−^𝑅₀𝐷^𝛼_𝑡𝑅_ℎ𝑢(𝑡),𝜒_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡𝑅_ℎ𝑢(𝑡),𝜒_𝐿2(Ω)+𝑘_𝛽𝐶_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝑅_ℎ𝑢(𝑡),^𝑅_𝑥𝐷_𝑋^𝛽/2𝜒_𝐿2(Ω)+𝑘_𝛽𝐶_𝛽^𝑅_𝑥𝐷_𝑋^𝛽/2𝑅_ℎ𝑢(𝑡),^𝑅₀𝐷_𝑥^𝛽/2𝜒_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝑅_ℎ𝑢(𝑡),^𝑅_𝑥𝐷_𝑋^𝛾/2𝜒_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅_𝑥𝐷_𝑋^𝛾/2𝑅_ℎ𝑢(𝑡),^𝑅₀𝐷_𝑥^𝛾/2𝜒_𝐿2(Ω)=𝑓(𝑡),𝜒_𝐿2(Ω)−^𝑅₀𝐷^𝛼_𝑡𝑅_ℎ𝑢(𝑡),𝜒_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡𝑅_ℎ𝑢(𝑡),𝜒_𝐿2(Ω)+𝑘_𝛽𝐶_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝑅_ℎ𝑢(𝑡),^𝑅_𝑥𝐷_𝑋^𝛽/2𝜒_𝐿2(Ω)+𝑘_𝛽𝐶_𝛽^𝑅_𝑥𝐷_𝑋^𝛽/2𝑅_ℎ𝑢(𝑡),^𝑅₀𝐷_𝑥^𝛽/2𝜒_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝑅_ℎ𝑢(𝑡),^𝑅_𝑥𝐷_𝑋^𝛾/2𝜒_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅_𝑥𝐷_𝑋^𝛾/2𝑅_ℎ𝑢(𝑡),^𝑅₀𝐷_𝑥^𝛾/2𝜒_𝐿2(Ω)=𝑃₀𝐷_𝑡𝑢(𝑡),𝜒_𝐿2(Ω)−𝑃₀𝐷_𝑡𝑅_ℎ𝑢(𝑡),𝜒_𝐿2(Ω)=𝐼−𝑅_ℎ𝑃₀𝐷_𝑡𝑢(𝑡),𝜒_𝐿2(Ω)=−𝑃₀𝐷_𝑡𝜌(𝑡),𝜒_𝐿2(Ω). (76)
Choosing 𝜒=𝜗(𝑡) and integrating on both sides with respect to 𝑡 on [0,𝑇], we obtain 𝑃₀𝐷_𝑡𝜗(𝑡),𝜗(𝑡)_𝐿2(Ω)+𝑘_𝛽𝐶_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝜗(𝑡),^𝑅_𝑥𝐷_𝑋^𝛽/2𝜗(𝑡)_𝐿2(Ω)+𝑘_𝛽𝐶_𝛽^𝑅_𝑥𝐷_𝑋^𝛽/2𝜗(𝑡),^𝑅₀𝐷_𝑥^𝛽/2𝜗(𝑡)_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝜗(𝑡),^𝑅_𝑥𝐷_𝑋^𝛾/2𝜗(𝑡)_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅_𝑥𝐷_𝑋^𝛾/2𝜗(𝑡),^𝑅₀𝐷_𝑥^𝛾/2𝜗(𝑡)_𝐿2(Ω)=−𝑃₀𝐷_𝑡𝜌(𝑡),𝜗(𝑡)_𝐿2(Ω). (77)
By Lemmas 2–4, for any small 𝜖>0, we have 𝐴_𝛼^𝑅₀𝐷_𝑡^𝛼/2𝜗(𝑡),^𝑅₀𝐷_𝑡^𝛼/2𝜗(𝑡)_𝐿2(Ω)+𝐴_𝛼𝑖^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝜗(𝑡),^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝜗(𝑡)_𝐿2(Ω)+𝑘_𝛽𝐶_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝜗(𝑡),^𝑅_𝑥𝐷_𝑋^𝛽/2𝜗(𝑡)_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝛽/2𝜗(𝑡),^𝑅₀𝐷_𝑥^𝛽/2𝜗(𝑡)_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝜗(𝑡),^𝑅_𝑥𝐷_𝑋^𝛾/2𝜗(𝑡)_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝛾/2𝜗(𝑡),^𝑅₀𝐷_𝑥^𝛾/2𝜗(𝑡)_𝐿2(Ω)=𝐴_𝛼^𝑅₀𝐷_𝑡^𝛼/2𝜌(𝑡),^𝑅₀𝐷_𝑡^𝛼/2𝜗(𝑡)_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖𝐴_𝛼𝑖^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝜌(𝑡),^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝜗(𝑡)_𝐿2(Ω)≤𝐶_𝜖𝐴_𝛼^𝑅₀𝐷_𝑡^𝛼/2𝜌(𝑡)²_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖𝐴_𝛼𝑖^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝜌(𝑡)²_𝐿2(Ω)+𝜖𝐴_𝛼^𝑅₀𝐷_𝑡^𝛼/2𝜗(𝑡)²_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖𝐴_𝛼𝑖^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝜗(𝑡)²_𝐿2(Ω), (78) where 𝐶_𝜖 is a constant respect to 𝜖, 𝐴_𝛼=cos(𝜋𝛼/2) and 𝐴_𝛼𝑖=cos((𝜋𝛼_𝑖)/2),𝑖=1,…,𝑠.
For sufficiently small 𝜖>0, by (75), we obtain 𝑃₀𝐷_𝑡𝜗(𝑡),𝜗(𝑡)_𝐿2(Ω)≤𝐶_𝜖^𝑅₀𝐷_𝑡^𝛼/2𝜌(𝑡)²_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝜌(𝑡)²_𝐿2(Ω)≤𝐶_𝜖ℎ^{𝜂−𝛾/22𝑅}₀𝐷_𝑡^𝛼/2𝑢(𝑡)²_𝐻𝜂(Ω)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝑢(𝑡)²_𝐻𝜂(Ω). (79) Combining (75) with (79), we obtain that 𝑢−𝑢_{ℎ𝐻𝛾/2(Ω)}=𝜌+𝜗_{𝐻𝛾/2(Ω)}=𝑂ℎ^{𝜂−𝛾/2}. (80)

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 find 𝑢_ℎ(𝑡,⋅)∈𝑆_ℎ such that (63) holds.

In terms of the basis {𝜓_𝑖}^𝑀−1_𝑖=1⊆𝑆_ℎ, write 𝑢_ℎ(𝑡,𝑥)=∑^𝑀−1_𝑗=1𝑢_𝑗(𝑡)𝜓_𝑗(𝑥), and insert it into (63). After the time discretization on 𝑡_𝑛, in Section 3, we get ^𝑀−1_𝑗=1Δ𝑡^𝑛−𝛼_𝑘=1𝜔^(𝛼)_𝑘𝑛𝑈_𝑗^𝑛−𝑘𝜓_𝑗,𝜓_{𝑖𝐿2(Ω)}+^𝑠_𝑖=1𝑎_𝑖^𝑀−1_𝑗=1Δ𝑡^{−𝛼𝑖𝑛}_𝑘=1𝜔^(𝛼𝑖)_𝑘𝑛𝑈_𝑗^𝑛−𝑘𝜓_𝑗,𝜓_{𝑖𝐿2(Ω)}+𝑘_𝛽𝐶_𝛽^𝑀−1_𝑗=1𝑈^𝑛_𝑗^𝑅₀𝐷_𝑥^𝛽/2𝜓_𝑗,^𝑅_𝑥𝐷_𝑋^𝛽/2𝜓_{𝑖𝐿2(Ω)}+^𝑅_𝑥𝐷_𝑋^𝛽/2𝜓_𝑗,^𝑅₀𝐷_𝑥^𝛽/2𝜓_{𝑖𝐿2(Ω)}+𝑘_𝛾𝐶_𝛾^𝑀−1_𝑗=1𝑈^𝑛_𝑗^𝑅₀𝐷_𝑥^𝛾/2𝜓_𝑗,^𝑅_𝑥𝐷_𝑋^𝛾/2𝜓_{𝑖𝐿2(Ω)}+^𝑅_𝑥𝐷_𝑋^𝛾/2𝜓_𝑗,^𝑅₀𝐷_𝑥^𝛾/2𝜓_{𝑖𝐿2(Ω)}=𝑓,𝜓_{𝑖𝐿2(Ω)},∀𝜓_𝑖∈𝑆_ℎ,𝑖=1,…,𝑀−1. (81)

To obtain the value of {𝑈^𝑛_𝑗}^𝑀−1_𝑗=1, Let 𝑈_𝑗=(𝑈¹_𝑗,𝑈²_𝑗,…,𝑈^𝑛_𝑗)^𝑇. From (81), we obtain a vector equation Δ𝑡^−𝛼𝜔^(𝛼)_0𝑛𝑀𝑈_𝑛+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖𝜔^(𝛼𝑖)_0𝑛𝑀𝑈_𝑛+𝑘_𝛽𝐶_𝛽𝑆₁^(𝛽/2)+𝑆₂^(𝛽/2)𝑈_𝑛+𝑘_𝛾𝐶_𝛾𝑆₁^(𝛾/2)+𝑆₂^(𝛾/2)𝑈_𝑛=𝑀𝐹_𝑛−Δ𝑡^𝑛−𝛼_𝑘=1𝜔^(𝛼)_𝑘𝑛𝑀𝑈_𝑛−𝑘−^𝑠_𝑖=1Δ𝑡^{−𝛼𝑖𝑛}_𝑘=1𝜔^(𝛼𝑖)_𝑘𝑛𝑀𝑈_𝑛−𝑘, (82) where 𝑀={(𝜓_𝑗,𝜓_𝑖)_𝐿2(Ω)}^𝑀−1_𝑗,𝑖=1 is the mass matrix; 𝑆₁^(𝜃) and 𝑆₂^(𝜃) are stiffness matrices with 𝜃=𝛽/2,𝛾/2, as follows: 𝑆₁^(𝜃)=^𝑅₀𝐷^𝜃_𝑥𝜓_𝑗,^𝑅_𝑥𝐷^𝜃_𝑋𝜓_{𝑖𝐿2(Ω)}^𝑀−1_𝑗,𝑖=1,𝑆₂^(𝜃)=^𝑅_𝑥𝐷^𝜃_𝑋𝜓_𝑗,^𝑅₀𝐷^𝜃_𝑥𝜓_{𝑖𝐿2(Ω)}^𝑀−1_𝑗,𝑖=1, (83) and 𝐹_𝑁=(𝑓₁,…,𝑓_𝑀−1)^𝑇 is a vector valued function. Then, from (82), we can obtain the solution 𝑈_𝑛.

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 𝜔^(𝛼)_𝑘𝑛 in Section 3, we can obtain the following lemma easily.

Lemma 12. For 0<𝛼<1, the coefficients 𝜔^(𝛼)_𝑘𝑛, (𝑘=1,…,𝑛) satisfy the following properties: (i)𝜔^(𝛼)_0𝑛>0 and 𝜔^(𝛼)_𝑘𝑛<0 for 𝑘=1,2,…,𝑛, (ii)Γ(2−𝛼)∑^𝑛_𝑘=1𝜔^(𝛼)_𝑘𝑛=(1−𝛼)𝑛^−𝛼+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. Let 𝑈^𝑛 denote the approximation to 𝑢_ℎ at 𝑡=𝑡_𝑛 and 𝜒(⋅)=𝑈^𝑛(⋅) and the right-hand side 𝑓=0. From (63), we have Δ𝑡^−𝛼1Γ()𝑈2−𝛼^𝑛,𝑈^𝑛_𝐿2(Ω)+^𝑛−1_𝑘=1𝜔^(𝛼)_𝑘𝑛𝑈^𝑛−𝑘,𝑈^𝑛_𝐿2(Ω)+𝜔^(𝛼)_𝑛𝑛𝑈⁰,𝑈^𝑛_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖1Γ2−𝛼_𝑖𝑈^𝑛,𝑈^𝑛_𝐿2(Ω)+^𝑛−1_𝑘=1𝜔^(𝛼𝑖)_𝑘𝑛𝑈^𝑛−𝑘,𝑈^𝑛_𝐿2(Ω)+𝜔^(𝛼𝑖)_𝑛𝑛𝑈⁰,𝑈^𝑛_𝐿2(Ω)+𝑘_𝛽𝐶_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝑈^𝑛,^𝑅_𝑥𝐷_𝑋^𝛽/2𝑈^𝑛_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝛽/2𝑈^𝑛,^𝑅₀𝐷_𝑥^𝛽/2𝑈^𝑛_𝐿2(Ω)+𝑘_𝛾𝐶_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝑈^𝑛,^𝑅_𝑥𝐷_𝑋^𝛾/2𝑈^𝑛_𝐿2(Ω)+^𝑅_𝑥𝐷_𝑋^𝛾/2𝑈^𝑛,^𝑅₀𝐷_𝑥^𝛾/2𝑈^𝑛_𝐿2(Ω)=0. (84) Using the Cauchy-Schwarz inequality, ±(𝑈^𝑛−𝑘,𝑈^𝑛)≤(1/2)(‖𝑈^𝑛−𝑘‖²_𝐿2(Ω)+‖𝑈^𝑛‖²_𝐿2(Ω)) for 𝑘=0,1,2,…,𝑛, by Definition 1 and Lemma 12, we obtain Δ𝑡^−𝛼()()𝑛2Γ2−𝛼1+1−𝛼^−𝛼+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖2Γ2−𝛼_𝑖1+1−𝛼_𝑖𝑛^−𝛼𝑖𝑈^𝑛2_𝐿2(Ω)+𝑘_𝛽^𝑅₀𝐷_𝑥^𝛽/2𝑈^𝑛2_𝐿2(Ω)+𝑘_𝛾^𝑅₀𝐷_𝑥^𝛾/2𝑈^𝑛2_𝐿2(Ω)≤Δ𝑡^−𝛼2−^𝑛−1_𝑘=1𝜔^(𝛼)_𝑘𝑛𝑈^𝑛−𝑘2_𝐿2(Ω)−𝜔^(𝛼)_𝑛𝑛𝑈⁰²_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖2−^𝑛−1_𝑘=1𝜔^(𝛼𝑖)_𝑘𝑛𝑈^𝑛−𝑘2_𝐿2(Ω)−𝜔^(𝛼𝑖)_𝑛𝑛𝑈⁰²_𝐿2(Ω). (85)
We prove the stability of (63) by induction. From the beginning, we have Δ𝑡^−𝛼()((+2Γ2−𝛼1+1−𝛼))_𝑠𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖2Γ2−𝛼_𝑖1+1−𝛼_𝑖𝑈¹²_𝐿2(Ω)≤Δ𝑡^−𝛼()((+2Γ2−𝛼1−1−𝛼))_𝑠𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖2Γ2−𝛼_𝑖1−1−𝛼_𝑖𝑈⁰²_𝐿2(Ω). (86) The induction basis ‖𝑈¹‖_𝐿2(Ω)≤‖𝑈⁰‖_𝐿2(Ω) is presupposed. For the induction step, we have ‖𝑈^𝑛‖_𝐿2(Ω)≤‖𝑈^𝑛−1‖_𝐿2(Ω)≤⋯≤‖𝑈⁰‖_𝐿2(Ω). Then, using this result, by Lemma 12, we obtain Δ𝑡^−𝛼1+(1−𝛼)(𝑛+1)^−𝛼()+2Γ2−𝛼_𝑠𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖1+1−𝛼_𝑖(𝑛+1)^−𝛼𝑖2Γ2−𝛼_𝑖𝑈^𝑛+12_𝐿2(Ω)≤Δ𝑡^−𝛼1−(1−𝛼)(𝑛+1)^−𝛼()+2Γ2−𝛼_𝑠𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖1−1−𝛼_𝑖(𝑛+1)^−𝛼𝑖2Γ2−𝛼_𝑖𝑈⁰²_𝐿2(Ω). (87) Here, 0<1−𝛼<1 and 0<1−𝛼_𝑖<1 for 𝑖=1,…,𝑠. After squaring at both sides of the above inequality, we obtain ‖𝑈^𝑛+1‖_𝐿2(Ω)≤‖𝑈⁰‖_𝐿2(Ω).

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 Δ𝑡 and space step size Δ𝑥, which have been shown in Theorems 9 and 11.

Example 1. Consider MT-TS-RFADE, in 𝑡∈[0,0.5],𝑥∈[0,1], ^𝐶₀𝐷^𝛼_𝑡𝑢(𝑡,𝑥)+^𝐶₀𝐷^𝛼1_𝑡𝑢(𝑡,𝑥)−^𝑅𝐷^𝛽_|𝑥|𝑢()−𝑡,𝑥_𝑅𝐷^𝛾_|𝑥|𝑢()(),𝑢()𝑡,𝑥=𝑓𝑡,𝑥0,𝑥=𝑥²(1−𝑥)²(),𝑢()()[].,𝑥∈0,1𝑡,0=𝑢𝑡,1=0,𝑡∈0,0.5 (88) The exact solution of (88) is 𝑢(𝑡,𝑥)=(2𝑡−1)²𝑥²(1−𝑥)². From the definition of the Riemann-Liouville differential operator, it holds ^𝑅_𝑎𝐷^𝛼_𝑡(𝑡−𝑎)^𝛽=Γ𝛽+1Γ()𝛽+1−𝛼𝑡−𝑎^𝛽−𝛼, (89) where 𝑎 is a positive constant. From (89), we can choose right-hand side function 𝑓(𝑡,𝑥) to satisfy (88).
Choosing 𝛼=0.9 and 𝛼₁=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 Δ𝑥=0.05 and time step size Δ𝑡=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 Example 1.

Example 2. Consider MT-TS-RFADE (88) with conditions as follows: 𝑢()[],𝑢()()[].0,𝑥=0,𝑥∈0,1𝑡,0=𝑢𝑡,1=0,𝑡∈0,0.5 (90) Let the exact solution is 𝑢(𝑡,𝑥)=𝑡^0.9𝑥²(1−𝑥)². We choose 𝛼₁=0.3,𝛼=0.4,𝛽=0.4,𝛾=1.7 and obtain the numerical solution and exact solution when 𝑡=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 Δ𝑥=0.1; time step size is Δ𝑡=0.01.

Figure 2 

Numerical solution and exact solution for Example 2.

Example 3. Consider MT-TS-RFADE (88) with the zero Dirichlet boundary conditions, for 𝑡∈[0,1], 𝑥∈[0,1]. We require that the exact solution is 𝑢(𝑡,𝑥)=sin(2𝜋𝑡)𝑥²(1−𝑥)².

For this example, in the first test, we obtain the numerical solution and exact solution when 𝑡=0.05, 0.1, 0.25, 0.4 in Figure 3, where we choose time step size Δ𝑡=0.01 and space step size Δ𝑥=0.01 with 𝛼₁=0.3, 𝛼=0.5, and 𝛽=0.9, 𝛾=1.9.

Figure 3 

Numerical solution and exact solution for Example 3.

In the second test, we check the convergence rates of numerical solutions with respect to the fractional orders 𝛼₁,𝛼,𝛽,  and 𝛾. We fix 𝛼₁=0.2,𝛼=0.8,𝛽=0.8, and 𝛾=1.8 and choose Δ𝑥=0.001 which is small enough such that the space discretization errors are negligible as compared with the time errors. Choosing step size Δ𝑡=1/2^𝑖(𝑖=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 Δ𝑡=0.001 and choosing Δ𝑥=1/2^𝑖(𝑖=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 𝐿₂-norm and 𝐻¹-norm of errors in Figure 4.

Table 1 

Convergence rate in time for Example 3.

Table 2 

Convergence rate in space for Example 3.

(a)

(b)

Figure 4 

𝐻1-norm and 𝐿2-norm of errors for Example 3; here, Δ𝑥=0.001 (a) and Δ𝑡=0.001 (b).

Acknowledgments

The authors are grateful to the referees for their valuable comments. This work is supported by the National Natural Science Foundation of China (11101109 and 11271102), the Natural Science Foundation of Hei-long-jiang Province of China (A201107), and SRF for ROCS, SEM.