Stability and Convergence of an Effective Finite Element Method for Multiterm Fractional Partial Differential Equations

Abstract

A finite element method (FEM) for multiterm fractional partial differential equations (MT-FPDEs) is studied for obtaining a numerical solution effectively. The weak formulation for MT-FPDEs and the existence and uniqueness of the weak solutions are obtained by the well-known Lax-Milgram theorem. The Diethelm fractional backward difference method (DFBDM), based on quadrature for the time discretization, and FEM for the spatial discretization have been applied to MT-FPDEs. The stability and convergence for numerical methods are discussed. The numerical examples are given to match well with the main conclusions.

1. Introduction

In recent years, the numerical treatment and supporting analysis of fractional order differential equations has become an important research topic that offers great potential. The FEMs for fractional partial differential equations have been studied by many authors (see [1–3]). All of these papers only considered single-term fractional equations, where they only had one fractional differential operator. In this paper, we consider the MT-FPDEs, which include more than one fractional derivative. Some authors also considered solving linear problems with multiterm fractional derivatives (see [4, 5]). This motivates us to consider their effective numerical solutions for such MT-FPDEs, which have been proposed in [6, 7].

Let Ω=(0,𝑋)^𝑑, where 𝑑≥1 is the space dimension. We consider the MT-FPDEs with the Caputo time fractional derivatives as follows: 𝑃^𝐶𝐷_𝑡𝑢(𝑡,𝑥)−Δ_𝑥𝑢()()[]𝑢()𝑡,𝑥=𝑓𝑡,𝑥,𝑡∈0,𝑇,𝑥∈Ω,(1)0,𝑥=𝑢₀(𝑥)𝑢()[],𝑥∈Ω,(2)𝑡,𝑥=0,𝑡∈0,𝑇,𝑥∈𝜕Ω,(3) where the operator 𝑃(^𝐶𝐷_𝑡)𝑢(𝑡,𝑥) is defined as 𝑃^𝐶𝐷_𝑡𝑢(𝑡,𝑥)=^𝐶₀𝐷^𝛼_𝑡+^𝑠_𝑖=1𝑎_𝑖^𝐶₀𝐷^𝛼𝑖_𝑡𝑢(𝑡,𝑥),(4) with 0<𝛼_𝑠<𝛼_𝑠−1<⋯<𝛼₁<𝛼<1 and {𝑎_𝑖>0}^𝑠_𝑖=1. Here ^𝐶₀𝐷^𝛼_𝑡𝑢(𝑡,𝑥) denotes the left Caputo fractional derivative with respect to the time variable 𝑡 and Δ_𝑥 denotes the Laplace operator with respect to the space variable 𝑥.

Some numerical methods have been considered for solving the multiterm fractional differential equations. In [8], Liu et al. investigate some effective numerical methods for time fractional wave-diffusion and diffusion equations: ^𝐶₀𝐷^𝛼_𝑡𝑢(𝑡,𝑥)−𝑘Δ_𝑥𝑢(𝑡,𝑥)=𝑓(𝑡,𝑥),0<𝑥<𝐿,𝑡>0,(5) where 𝑘 and 𝐿 are arbitrary positive constants and 𝑓(𝑡,𝑥) is a sufficiently smooth function. The authors consider the implicit finite difference methods (FDMs) and prove that it is unconditionally stable. The error estimate of the FDM is 𝑂(Δ𝑡+Δ𝑡^2−𝛼+Δ𝑥), where Δ𝑡 and Δ𝑥 are the time and space step size, respectively. They also investigate the fractional predictor-corrector methods (FPCMs) of the Adams-Moulton methods for multiterm time fractional differential equations (1) with order {𝛼_𝑖}_{𝑖=1,…,𝑠} by solving the equivalent Volterra integral equations. The error estimate of the FPCM is 𝑂(Δ𝑡+Δ𝑡^{1+min{𝛼𝑙}}+Δ𝑥²). In recent years, there are some articles for the predictor-correction method for initial-value problems (see [9–14]). For the application of the FDMs, there have been many research articles as follows. In [15–20], Simos et al. investigate the numerical methods for solving the Schrödinger equation. In [21–24], the Runge-Kutta methods are considered and applied to get the numerical solution of orbital problems. For long-time integration, the Newton-Cotes formulae are considered in [25–27].

In [28], Badr investigate the FEM for linear multiterm fractional differential equations with one variable as follows: ^𝐶₀𝐷_𝑡^1+𝛼𝑢(𝑡)+^𝑠_𝑖=1𝐴_𝑖(𝑥)𝐷^𝛼𝑖𝑢(𝑥)(𝑡)𝛼=𝑓,𝛼≤𝑛,_𝑖<𝑛−1,0<𝑡<1,(6) where 𝐴_𝑖(𝑥) are known functions. The author gives the details of the modified Galerkin method for the above equations and makes the numerical example for checking the numerical method. In [29], Ford et al. consider the FEM for (5) with singular fractional order and obtain the error estimate 𝑂(Δ𝑡^2−𝛼+Δ𝑥²). In this paper, we follow the work in [29] and consider the FEM for solving MT-FPDEs (1)–(3). Then, we prove the stability and convergence of the FEM for MT-FPDEs and make the error estimate.

The paper is organized as follows. In Section 2, the weak formulation of the MT-FPDEs is given and the existence and uniqueness results for such problems are proved. In Section 3, we consider the convergence rate of time discretization of MT-FPDEs, based on the Diethelm fractional backward difference method (DFBDM). In Section 4, we propose an FEM based on the weak formulation and carry out the error analysis. In Section 5, the stability of this method is proven. Finally, the numerical examples are considered for matching well with the main conclusions.

2. Existence and Uniqueness

Let Γ(⋅) denote the gamma function. For any positive integer 𝑛 and 𝑛−1<𝛼<𝑛, the Caputo derivative are the Riemann-Liouville derivative are, respectively, defined as follows [30].(i)The left Caputo derivatives: ^𝐶₀𝐷^𝛼_𝑡𝑣(𝑡)1∶=Γ(𝑛−𝛼)^𝑡₀1(𝑡−𝜏)^{𝛼−𝑛+1}𝑑^𝑛𝑑𝜏^𝑛𝑣(𝜏)𝑑𝜏.(7)(ii)The left Riemann-Liouville derivatives: ^𝑅₀𝐷^𝛼_𝑡𝑣(𝑡)1∶=Γ()𝑑𝑛−𝛼𝑑𝑡^𝑛𝑛𝑡₀𝑣(𝜏)(𝑡−𝜏)^{𝛼−𝑛+1}𝑑𝜏.(8)(iii)The right Riemann-Liouville derivatives: ^𝑅_𝑡𝐷^𝛼_𝑇𝑣(𝑡)()∶=−1^𝑛Γ()𝑑𝑛−𝛼𝑑𝑡^𝑛𝑛𝑇_𝑡𝑣(𝜏)(𝜏−𝑡)^{𝛼−𝑛+1}𝑑𝜏.(9)

Let 𝐶^∞(0,𝑇) denote the space of infinitely differentiable functions on (0,𝑇) and 𝐶^∞₀(0,𝑇) denote the space of infinitely differentiable functions with compact support in (0,𝑇). We use the expression 𝐴≲𝐵 to mean that 𝐴≤𝑐𝐵 when 𝑐 is a positive real number and use the expression 𝐴≅𝐵 to mean that 𝐴≲𝐵≲𝐴. Let 𝐿₂(𝒬) be the space of measurable functions whose square is the Lebesgue integrable in 𝒬 which may denote a domain 𝒬=𝐼×Ω, 𝐼 or Ω. Here time domain 𝐼∶=(0,𝑇) and space domain Ω∶=(0,𝑋). The inner product and norm of 𝐿₂(𝒬) are defined by (𝑢,𝑣)_𝐿2(𝒬)∶=_𝒬𝑢𝑣𝑑𝒬,‖𝑢‖_𝐿2(𝒬)∶=(𝑢,𝑢)_𝐿^1/2_2(𝒬),∀𝑢,𝑣∈𝐿₂(𝒬).(10)

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(𝒬).(11) In the usual Sobolev space 𝐻^𝜎₀(𝒬), we also have the definition ‖𝑣‖_𝐻^𝜎_0(𝒬)∶=‖𝑣‖²_𝐿2(𝒬)+|𝑣|²_𝐻^𝜎_0(𝒬)^1/2,|𝑣|²_𝐻^𝜎_0(𝒬)∶=^𝑅0𝐷^𝜎_𝑡𝑣,^𝑅_𝑡𝐷^𝜎_𝑇𝑣_𝐿2(𝒬)().cos𝜋𝜎(12)

From [3], for 𝜎>0, 𝜎≠𝑛−1/2, the spaces ^𝑙𝐻^𝜎₀(𝒬), ^𝑟𝐻^𝜎₀(𝒬), and 𝐻^𝜎₀(𝒬) are equal, and their seminorms are all equivalent to |⋅|_𝐻^𝜎_0(𝒬). We first recall the following results.

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

From [3], we define the following space: 𝐵^𝛼/2(𝐼×Ω)=𝐻^𝛼/2𝐼,𝐿₂(Ω)∩𝐻^𝛼1/2𝐼,𝐿₂(Ω)∩⋯∩𝐻^𝛼𝑠/2𝐼,𝐿₂(Ω)∩𝐿₂𝐼,𝐻¹₀(Ω)=𝐻^𝛼/2𝐼,𝐿₂(Ω)∩𝐿₂𝐼,𝐻¹₀(Ω).(14)

Here 𝐵^𝛼/2(𝐼×Ω) is a Banach space with respect to the following norm: ‖𝑣‖_𝐵^𝛼/2_(𝐼×Ω)=‖𝑣‖²_{𝐻𝛼/2𝐼,𝐿2(Ω)}+‖𝑣‖²_{𝐿2𝐼,𝐻10(Ω)}^1/2,(15) where 𝐻^𝛼/2(𝐼,𝐿₂(Ω))∶={𝑣;‖𝑣(𝑡,⋅)‖_𝐿2(Ω)∈𝐻^𝛼/2(𝐼)}, endowed with the norm ‖𝑣‖_𝐻^𝛼/2_{(𝐼,𝐿2(Ω))}∶=‖𝑣(𝑡,⋅)‖_{𝐿2(Ω)𝐻𝛼/2(𝐼)}.(16)

Based on the relation equation between the left Caputo and the Riemann-Liouville derivative in [31], we can translate the Caputo problem to the Riemann-Liouville problem. Then, we consider the weak formulation of (1) as follows. For 𝑓∈𝐵^𝛼/2(𝐼×Ω)^′, find 𝑢(𝑡,𝑥)∈𝐵^𝛼/2(𝐼×Ω) such that 𝒜(𝑢,𝑣)=ℱ(𝑣),𝑣∈𝐵^𝛼/2(𝐼×Ω),(17) where the bilinear form is, by Lemma 1, 𝒜(𝑢,𝑣)∶=^𝑅₀𝐷_𝑡^𝛼/2𝑢,^𝑅_𝑡𝐷_𝑇^𝛼/2𝑣_{𝐿2(𝐼×Ω)}+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝑢,^𝑅_𝑡𝐷^𝛼𝑖_𝑇^/2𝑣_{𝐿2(𝐼×Ω)}+∇_𝑥𝑢,∇_𝑥𝑣_{𝐿2(𝐼×Ω)},(18) and the functional is ℱ(𝑣)∶=(𝑓,𝑣)_{𝐿2(𝐼×Ω)}, 𝑓(𝑡,𝑥)∶=𝑓(𝑡,𝑥)+(𝑢₀(𝑥)𝑡^−𝛼/Γ(1−𝛼))+∑^𝑠_𝑖=1𝑎_𝑖(𝑢₀(𝑥)𝑡^−𝛼𝑖/Γ(1−𝛼_𝑖)).

Based on the main results in Subsection 3.2 in [32], we can prove the following existence and uniqueness theorem.

Theorem 2. Assume that 0<𝛼<1 and 𝑓∈𝐵^𝛼/2(𝐼×Ω)^′. Then the system (17) has a unique solution in 𝐵^𝛼/2(𝐼×Ω). Furthermore, ‖𝑢‖_𝐵^𝛼/2_(𝐼×Ω)≲𝑓_{𝐵𝛼/2(𝐼×Ω)′}.(19)

Proof. The existence and uniqueness of the solution of (17) is guaranteed by the well-known Lax-Milgram theorem. The continuity of the bilinear form 𝒜 and the functional ℱ is obvious. Now we need to prove the coercivity of 𝒜 in the space 𝐵^𝛼/2(𝐼×Ω). From the equivalence of ^𝑙𝐻^𝛼₀(𝐼×Ω), ^𝑟𝐻^𝛼₀(𝐼×Ω) and 𝐻^𝛼₀(𝐼×Ω), for all 𝑢,𝑣∈𝐵^𝛼/2(𝐼×Ω), using the similar proof process in [32], we obtain 𝒜(𝑣,𝑣)≳^𝑅₀𝐷_𝑡^𝛼/2𝑣,^𝑅₀𝐷_𝑡^𝛼/2𝑣_{𝐿2(𝐼×Ω)}+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝑣,^𝑅₀𝐷^𝛼𝑖_𝑡^/2𝑣_{𝐿2(𝐼×Ω)}+∇_𝑥𝑢,∇_𝑥𝑣_{𝐿2(𝐼×Ω)}≳‖𝑣‖²_{𝐵𝛼/2()𝐼×Ω}.(20) Then we take 𝑣=𝑢 in (17) to get ‖𝑢‖²_{𝐵𝛼/2(𝐼×Ω)}≲(𝑓,𝑢)_{𝐿2(𝐼×Ω)} by the Schwarz inequality and the Poincaré inequality.

3. Time Discretization and Convergence

In this section, we consider DFBDM for the time discretization of (1)–(3), which is introduced in [33] for fractional ordinary differential equations. We can obtain the convergence order for the time discretization for the MT-FPDEs. Let 𝐴=−Δ_𝑥, 𝐷(𝐴)=𝐻¹₀(Ω)∩𝐻²(Ω). Let 𝑢(𝑡), 𝑓(𝑡), and 𝑢(0) denote the one-variable functions as 𝑢(𝑡,⋅), 𝑓(𝑡,⋅), and 𝑢(0,⋅), respectively. Then (1) can be written in the abstract form, for 0<𝑡<𝑇, 0<𝛼_𝑠<⋯<𝛼₁<𝛼<1, with initial value 𝑢(0)=𝑢₀. Now we have ^𝑅₀𝐷^𝛼_𝑡𝑢−𝑢₀(𝑡)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡𝑢−𝑢₀(𝑡)+𝐴𝑢(𝑡)=𝑓(𝑡).(21)

Let 0=𝑡₀<𝑡₁<⋯<𝑡_𝑁=𝑇 be a partition of [0,𝑇]. Then, for fixed 𝑡_𝑗, 𝑗=1,2,…,𝑁, we have ^𝑅₀𝐷^𝛼_𝑡𝑢−𝑢₀𝑡_𝑗=𝑡_𝑗^−𝛼Γ(−𝛼)¹₀𝑔(𝜏)𝜏^−1−𝛼𝑑𝜏,(22) where 𝑔(𝜏)=𝑢(𝑡_𝑗−𝑡_𝑗𝜏)−𝑢₀. Here, the integral is a Hadamard finite-part integral in [33] and [34].

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𝛼^(𝛼)_𝑘𝑗𝑔𝑘𝑗+𝑅_𝑗^(𝛼)𝑔,(23) where the weights 𝛼^(𝛼)_𝑘𝑗 are 𝛼(1−𝛼)𝑗^−𝛼𝛼^(𝛼)_𝑘𝑗=−1,for𝑘=0,2𝑘^1−𝛼−(𝑘−1)^1−𝛼−(𝑘+1)^1−𝛼()𝑘,for𝑘=1,2,…,𝑗−1,𝛼−1^−𝛼−(𝑘−1)^1−𝛼+𝑘^1−𝛼,for𝑘=𝑗,(24) and the remainder term 𝑅_𝑗^(𝛼)(𝑔) satisfies ‖𝑅_𝑗^(𝛼)(𝑔)‖≤𝛾_𝛼𝑗^𝛼−2sup_{0≤𝑡≤𝑇}‖𝑔^′′(𝑡)‖, where 𝛾_𝛼>0 is a constant.

Thus, for 𝜔^(𝛼)_𝑘𝑗=𝑗^−𝛼𝛼^(𝛼)_𝑘𝑗/Γ(−𝛼), we have ^𝑅₀𝐷^𝛼_𝑡𝑢−𝑢₀𝑡_𝑗=Δ𝑡^𝑗−𝛼_𝑘=0𝜔^(𝛼)_𝑘𝑗𝑢𝑡_𝑗−𝑡_𝑘(0)+𝑡−𝑢_𝑗−𝛼Γ()𝑅−𝛼^𝑗(𝛼)𝑔.(25)

Let 𝑡=𝑡_𝑗, we can write (21) as Δ𝑡^𝑗−𝛼_𝑘=0𝜔^(𝛼)_𝑘𝑗𝑢𝑡_𝑗−𝑡_𝑘(0)+−𝑢_𝑠𝑖=1𝑎_𝑖Δ𝑡^{−𝛼𝑖𝑗}_𝑘=0𝜔^(𝛼𝑖)_𝑘𝑗𝑢𝑡_𝑗−𝑡_𝑘−𝑢(0)+𝐴𝑢𝑡_𝑗=𝑓𝑡_𝑗−𝑡_𝑗^−𝛼Γ()𝑅−𝛼^𝑗(𝛼)𝑔−^𝑠_𝑖=1𝑎_𝑖𝑡^−𝛼𝑖_𝑗Γ−𝛼_𝑖𝑅^(𝛼𝑖)_𝑗𝑔,𝑗=1,2,3,….(26) Denote 𝑈^𝑗 as the approximation of 𝑢(𝑡_𝑗) and 𝑓_𝑗=𝑓(𝑡_𝑗). We obtain the following equation: Δ𝑡^𝑗−𝛼_𝑘=0𝜔^(𝛼)_𝑘𝑗𝑈^𝑗−𝑘−𝑈⁰+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{−𝛼𝑖𝑗}_𝑘=0𝜔^(𝛼𝑖)_𝑘𝑗𝑈^𝑗−𝑘−𝑈⁰+𝐴𝑈^𝑗=𝑓_𝑗.(27)

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

Let 𝑒^𝑗=𝑢(𝑡_𝑗)−𝑈^𝑗 denote the error in 𝑡_𝑗. Then we have the following error estimate.

Theorem 4. Let 𝑈^𝑗 and 𝑢(𝑡_𝑗) be the solutions of (27) and (21), respectively. Then one has ‖𝑈^𝑗−𝑢(𝑡_𝑗)‖≲Δ𝑡^2−𝛼.

Proof. Subtracting (27) from (26), we obtain the error equation Δ𝑡^𝑗−𝛼_𝑘=0𝜔^(𝛼)_𝑘𝑗𝑒^𝑗−𝑘−𝑒⁰+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{−𝛼𝑖𝑗}_𝑘=0𝜔^(𝛼𝑖)_𝑘𝑗𝑒^𝑗−𝑘−𝑒⁰+𝐴𝑒^𝑗𝑡=−^𝑗−𝛼Γ()𝑅−𝛼^𝑗(𝛼)𝑔−^𝑠_𝑖=1𝑎_𝑖𝑡^−𝛼𝑖_𝑗Γ−𝛼_𝑖𝑅^(𝛼𝑖)_𝑗𝑔.(28)
Note that 𝑒⁰=𝑢(0)−𝑈⁰=0. Denote 𝑒^𝑗=𝛼^(𝛼)_0𝑗+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{𝛼−𝛼𝑖}Γ(−𝛼)Γ−𝛼_𝑖𝛼^(𝛼𝑖)_0𝑗+𝐴𝑡^𝛼Γ(−𝛼)⁻¹×^𝑗_𝑘=1𝛼^(𝛼)_𝑘𝑗𝑒^𝑗−𝑘+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{𝛼−𝛼𝑖}Γ(−𝛼)Γ−𝛼_𝑖^𝑗_𝑘=1𝛼^(𝛼𝑖)_𝑘𝑗𝑒^𝑗−𝑘−𝑅_𝑗^(𝛼)𝑔−^𝑠_𝑖=1𝑎_𝑖Γ(−𝛼)Γ−𝛼_𝑖𝑡^{𝛼−𝛼𝑖}_𝑗𝑅^(𝛼𝑖)_𝑗𝑔.(29) Let ‖⋅‖ denote the 𝐿₂-norm, then we have 𝑒^𝑗≤𝛼^(𝛼)_0𝑗+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{𝛼−𝛼𝑖}Γ(−𝛼)Γ−𝛼_𝑖𝛼^(𝛼𝑖)_0𝑗+𝐴𝑡^𝛼Γ(−𝛼)⁻¹×^𝑗_𝑘=1𝛼^(𝛼)_𝑘𝑗𝑒^𝑗−𝑘+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{𝛼−𝛼𝑖}Γ(−𝛼)Γ−𝛼_𝑖^𝑗_𝑘=1𝛼^(𝛼𝑖)_𝑘𝑗𝑒^𝑗−𝑘+𝑅_𝑗^(𝛼)𝑔+^𝑠_𝑖=1𝑎_𝑖Γ(−𝛼)Γ−𝛼_𝑖𝑡^{𝛼−𝛼𝑖}_𝑗𝑅^(𝛼𝑖)_𝑗𝑔.(30)
Note that 𝐴 is a positive definite elliptic operator with all of eigenvalues 𝜆>0. Since 𝛼^(𝛼)_0𝑗<0 and Γ(−𝛼)<0, we have 𝛼^(𝛼)_0𝑗+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{𝛼−𝛼𝑖}Γ(−𝛼)Γ−𝛼_𝑖𝛼^(𝛼𝑖)_0𝑗+𝐴𝑡^𝛼Γ(−𝛼)⁻¹=sup_𝜆>0𝛼^(𝛼)_0𝑗+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{𝛼−𝛼𝑖}Γ(−𝛼)Γ−𝛼_𝑖𝛼^(𝛼𝑖)_0𝑗+𝜆𝑡^𝛼Γ(−𝛼)⁻¹≲−𝛼^(𝛼)_0𝑗−^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{𝛼−𝛼𝑖}Γ(−𝛼)Γ−𝛼_𝑖𝛼^(𝛼𝑖)_0𝑗⁻¹.(31) Hence, 𝑒^𝑗≤𝛼(1−𝛼)𝑗^𝑗−𝛼_𝑘=1𝛼^(𝛼)_𝑘𝑗𝑒^𝑗−𝑘+^𝑠_𝑖=1𝛼_𝑖1−𝛼_𝑖𝑗^{−𝛼𝑖𝑗}_𝑘=1𝛼^(𝛼𝑖)_𝑘𝑗𝑒^𝑗−𝑘+𝛼(1−𝛼)𝛾_𝛼𝑛⁻²sup_{0≤𝑡≤𝑇}𝑢^′′+^𝑠_𝑖=1𝛼_𝑖1−𝛼_𝑖𝛾_𝛼𝑖𝑛⁻²sup_{0≤𝑡≤𝑇}𝑢^′′.(32)
Denote 𝑑₁=1 and 𝑑_𝑗=1+𝛼(1−𝛼)𝑗^{−𝛼𝑗−1}_𝑘=1𝛼^(𝛼)_𝑘𝑗𝑑_𝑗−𝑘𝑑,𝑗=2,3,…,𝑛,_𝑖𝑗=1+𝛼_𝑖1−𝛼_𝑖𝑗^{−𝛼𝑖𝑗−1}_𝑘=1𝛼^(𝛼𝑖)_𝑘𝑗𝑑_𝑗−𝑘,𝑗=2,3,…,𝑛,(33) where 𝑖=1,2,…,𝑠. By induction and Lemma 3, then we have 𝑒^𝑗≲𝛼(1−𝛼)𝑛⁻²sup_{0≤𝑡≤𝑇}𝑢^′′(𝑡)⋅𝑑_𝑗+^𝑠_𝑖=1𝛼_𝑖1−𝛼_𝑖𝑛⁻²sup_{0≤𝑡≤𝑇}𝑢^′′(𝑡)⋅𝑑^𝑖_𝑗≲𝑛⁻²sup_{0≤𝑡≤𝑇}𝑢^′′(𝑡())sin𝜋𝛼𝜋𝑗^𝛼+^𝑠_𝑖=1𝑛⁻²sup_{0≤𝑡≤𝑇}𝑢^′′(𝑡)sin𝜋𝛼_𝑖𝜋𝑗^𝛼𝑖≲Δ𝑡^2−𝛼+^𝑠_𝑖=1Δ𝑡^2−𝛼𝑖.(34)

4. Space Discretization and Convergence

In this section, we will consider the space discretization for MT-FPDEs (1) and show the complete process and details of numerical scheme. The variational form of (1) is to find 𝑢(𝑡,⋅)∈𝐻¹₀(Ω), such that, for all 𝑣∈𝐻¹₀(Ω), ^𝑅₀𝐷^𝛼_𝑡𝑢(𝑡,𝑥),𝑣_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡𝑢(𝑡,𝑥),𝑣_𝐿2(Ω)+∇_𝑥𝑢,∇_𝑥𝑣=𝑓(𝑡,𝑥),𝑣_𝐿2(Ω).(35)

Let ℎ denote the maximal length of intervals in Ω and let 𝑟 be any nonnegative integer. We denote the norm in 𝐻^𝑟(Ω) by ‖⋅‖_𝐻^𝑟_(Ω). Let 𝑆_ℎ⊂𝐻^𝑟₀ be a family of finite element spaces with the accuracy of order 𝑟≥2, that is, 𝑆_ℎ 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 ℎ, 𝑣∈𝐻^𝑏(Ω)∩𝐻¹₀(Ω), inf_{𝜒∈𝑆ℎ}𝑣−𝜒_𝐿2(Ω)+ℎ∇_𝑥𝑣−𝜒_𝐿2(Ω)≤𝐶ℎ^𝑏‖𝑣‖_𝐻^𝑏_(Ω),1≤𝑏≤𝑟.(36)

The semidiscrete problem of (1) is to find the approximate solution 𝑢_ℎ(𝑡)=𝑢_ℎ(𝑡,⋅)∈𝑆_ℎ and 𝑓(𝑡)=𝑓(𝑡,⋅) for each 𝑡 such that ^𝑅₀𝐷^𝛼_𝑡𝑢_ℎ(𝑡),𝜒_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖^𝑅₀𝐷^𝛼𝑖_𝑡𝑢_ℎ(𝑡),𝜒_𝐿2(Ω)+∇_𝑥𝑢_ℎ(𝑡),∇_𝑥𝜒_𝐿2(Ω)=𝑓(𝑡),𝜒_𝐿2(Ω),∀𝜒∈𝑆_ℎ.(37)

Let 𝑈^𝑁=𝑢_ℎ(𝑡_𝑁,𝑥). After the time discretization, we have Δ𝑡^𝑁−𝛼_𝑘=0𝜔^(𝛼)_𝑘𝑁𝑈^𝑁,𝜒_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{−𝛼𝑖𝑁}_𝑘=0𝜔^(𝛼𝑖)_𝑘𝑁𝑈^𝑁,𝜒_𝐿2(Ω)+∇_𝑥𝑈^𝑁,∇_𝑥𝜒_𝐿2(Ω)=𝑓(𝑡),𝜒_𝐿2(Ω),∀𝜒∈𝑆_ℎ.(38)

In terms of the basis {𝜓_𝑚}^𝑀−1_𝑚=1⊆𝑆_ℎ, choosing 𝜒=𝜓_𝑚, writing 𝑢_ℎ𝑡_𝑁,𝑥=^𝑀−1_𝑗=1𝑈^𝑁_𝑗𝜓_𝑗(𝑥),(39) and inserting it into (38), one obtains ^𝑀−1_𝑗=1Δ𝑡^𝑁−𝛼_𝑘=1𝜔^(𝛼)_𝑘𝑁𝑈_𝑗^𝑁−𝑘𝜓_𝑗,𝜓_{𝑚𝐿2(Ω)}+^𝑠_𝑖=1𝑎_𝑖^𝑀−1_𝑗=1Δ𝑡^{−𝛼𝑖𝑁}_𝑘=1𝜔^(𝛼𝑖)_𝑘𝑁𝑈_𝑗^𝑁−𝑘𝜓_𝑗,𝜓_{𝑚𝐿2(Ω)}+^𝑀−1_𝑗=1𝑈^𝑁_𝑗∇_𝑥𝜓_𝑗,∇_𝑥𝜓_{𝑚𝐿2(Ω)}=𝑓,𝜓_{𝑚𝐿2(Ω)},𝑚=1,2,…,𝑀−1.(40)

Let 𝑈_𝑁=(𝑈^𝑁₁,𝑈^𝑁₂,…,𝑈^𝑁_𝑀−1)^𝑇. From (40), we obtain a vector equation Ψ₁Δ𝑡^−𝛼𝜔^(𝛼)_0𝑁𝑈_𝑁+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖𝜔^(𝛼𝑖)_0𝑁𝑈_𝑁+Ψ₂𝑈_𝑁=Ψ₁𝐹_𝑁−Δ𝑡^𝑁−𝛼_𝑘=1𝜔^(𝛼)_𝑘𝑁𝑈_𝑁−𝑘−^𝑠_𝑖=1𝑎_𝑖Δ𝑡^{−𝛼𝑖𝑁}_𝑘=1𝜔^(𝛼𝑖)_𝑘𝑁𝑈_𝑁−𝑘,(41) where initial condition is 𝑈₀=𝑢(0,𝑥), Ψ₁∶={(𝜓_𝑗,𝜓_𝑚)_𝐿2(Ω)}^𝑀−1_𝑗,𝑚=1 is the mass matrix, Ψ₂ is stiffness matrix as Ψ₂∶={(∇_𝑥𝜓_𝑗,∇_𝑥𝜓_𝑚)_𝐿2(Ω)}^𝑀−1_𝑗,𝑚=1, and 𝐹_𝑁∶=(𝑓₁,…,𝑓_𝑀−1)^𝑇 is a vector valued function. Then, we can obtain the solution 𝑈_𝑁 at 𝑡=𝑡_𝑁.

Let 𝑅_ℎ∶𝐻¹(Ω)→𝑆_ℎ be the elliptic projection, defined by (∇_𝑥𝑅_ℎ𝑢,∇_𝑥𝜒)_𝐿2(Ω)=(∇_𝑥𝑢,∇_𝑥𝜒)_𝐿2(Ω), for all 𝜒∈𝑆_ℎ.

Lemma 5 (see [35]). Assume that (36) holds, then with 𝑅_ℎ and 𝑣∈𝐻^𝑏(Ω)∩𝐻¹₀(Ω), we have ‖𝑅_ℎ𝑣−𝑣‖_𝐿2(Ω)+ℎ‖∇_𝑥(𝑅_ℎ𝑣−𝑣)‖_𝐿2(Ω)≤𝐶ℎ^𝑏‖𝑣‖_𝐻^𝑏_(Ω) for 1≤𝑏≤𝑟.

In virtue of the standard error estimate for the FEM of MT-FPDEs, one has the following theorem which can be proved easily by Lemma 5 and the similar proof in [35].

Theorem 6. For 0<𝛼_𝑠<⋯<𝛼₁<𝛼<1, let 𝑢_ℎ∈𝑆_ℎ and 𝑢(𝑡,⋅)∈𝐻¹₀(Ω) be, respectively, the solutions of (37) and (1), then ‖𝑢−𝑢_ℎ‖_𝐿2(Ω)≲ℎ²‖𝑢‖_𝐿2(Ω).

5. Stability of the Numerical Method

In this section, we analyze the stability of the FEM for MT-FPEDs (1)–(3). Now we do some preparations before proving the stability of the method. Based on the definition of coefficients 𝜔^(𝛼)_𝑘𝑗 in Section 3, we can obtain the following lemma easily.

Lemma 7. 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-FPDEs in this section as follows.

Theorem 8. The FEM defined as in (38) is unconditionally stable.

Proof. In (38), let 𝜒(⋅)=𝑈^𝑗(⋅) at 𝑡=𝑡_𝑗 and the right hand 𝑓=0. We haveΔ𝑡^−𝛼1Γ()𝑈2−𝛼^𝑗,𝑈^𝑗_𝐿2(Ω)+^𝑗−1_𝑘=1𝜔^(𝛼)_𝑘𝑗𝑈^𝑗−𝑘,𝑈^𝑗_𝐿2(Ω)+𝜔^(𝛼)_𝑗𝑗𝑈⁰,𝑈^𝑗_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖1Γ2−𝛼_𝑖𝑈^𝑗,𝑈^𝑗_𝐿2(Ω)+^𝑗−1_𝑘=1𝜔^(𝛼𝑖)_𝑘𝑗𝑈^𝑗−𝑘,𝑈^𝑗_𝐿2(Ω)+𝜔^(𝛼𝑖)_𝑗𝑗𝑈⁰,𝑈^𝑗_𝐿2(Ω)+∇_𝑥𝑈^𝑗,∇_𝑥𝑈^𝑗_𝐿2(Ω)=0.(42)
Using, Cauchy-Schwarz inequality, ±(𝑈^𝑗−𝑘,𝑈^𝑗)≤(1/2)(‖𝑈^𝑗−𝑘‖²_𝐿2(Ω)+‖𝑈^𝑗‖²_𝐿2(Ω)) for 𝑘=0,1,2,…,𝑁 and Lemma 7, we get Δ𝑡^−𝛼()()𝑗2Γ2−𝛼1+1−𝛼^−𝛼+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖2Γ2−𝛼_𝑖1+1−𝛼_𝑖𝑗^−𝛼𝑖𝑈^𝑗2_𝐿2(Ω)+∇_𝑥𝑈^𝑗2_𝐿2(Ω)≤Δ𝑡^−𝛼2−^𝑗−1_𝑘=1𝜔^(𝛼)_𝑘𝑗𝑈^𝑗−𝑘2_𝐿2(Ω)−𝜔^(𝛼)_𝑗𝑗𝑈⁰²_𝐿2(Ω)+^𝑠_𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖2−^𝑗−1_𝑘=1𝜔^(𝛼𝑖)_𝑘𝑗𝑈^𝑗−𝑘2_𝐿2(Ω)−𝜔^(𝛼𝑖)_𝑗𝑗𝑈⁰²_𝐿2(Ω).(43)
We prove the stability of (37) by induction. Since when 𝑗=1, we have Δ𝑡^−𝛼()((+2Γ2−𝛼1+1−𝛼))_𝑠𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖2Γ2−𝛼_𝑖1+1−𝛼_𝑖𝑈¹²_𝐿2(Ω)≤Δ𝑡^−𝛼()((+2Γ2−𝛼1−1−𝛼))_𝑠𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖2Γ2−𝛼_𝑖1−1−𝛼_𝑖𝑈⁰_𝐿2(Ω).(44) The induction basis ‖𝑈¹‖_𝐿2(Ω)≤‖𝑈⁰‖_𝐿2(Ω) is presupposed. For the induction step, we have ‖𝑈^𝑗‖_𝐿2(Ω)≤‖𝑈^𝑗−1‖_𝐿2(Ω)≤⋯≤‖𝑈⁰‖_𝐿2(Ω). Then using this result, by Lemma 7, we obtain Δ𝑡^−𝛼1+(1−𝛼)𝑗+1^−𝛼()+2Γ2−𝛼_𝑠𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖2Γ2−𝛼_𝑖1+1−𝛼_𝑖𝑗+1^−𝛼𝑖𝑈^𝑗+12_𝐿2(Ω)≤Δ𝑡^−𝛼1−(1−𝛼)𝑗+1^−𝛼()+2Γ2−𝛼_𝑠𝑖=1𝑎_𝑖Δ𝑡^−𝛼𝑖2Γ2−𝛼_𝑖1−1−𝛼_𝑖𝑗+1^−𝛼𝑖𝑈⁰²_𝐿2(Ω).(45)
Here 0<1−𝛼<1. After squaring at both sides of the above inequality, we obtain ‖𝑈^𝑗+1‖_𝐿2(Ω)≤‖𝑈⁰‖_𝐿2(Ω).

6. Numerical Experiments

In this section, we present the numerical examples of MT-FPDEs to demonstrate the effectiveness of our theoretical analysis. The main purpose is to check the convergence behavior of numerical solutions with respect to Δ𝑡 and Δ𝑥, which have been shown in Theorem 4 and Theorem 6. It is noted that the method in [29] is a special case of the method in our paper for fractional partial differential equation with single fractional order. So, we just need to compare FEM in our paper with other existing methods in [8, 28].

Example 9. For 𝑡∈[0,𝑇], 𝑥∈(0,1), consider the MT-FPDEs with two variables as follows: ^𝐶₀𝐷^𝛼_𝑡𝑢(𝑡,𝑥)+^𝐶₀𝐷^𝛽_𝑡𝑢(𝑡,𝑥)−𝜕²_𝑥𝑢()(),𝑢()(),𝑢()()[],𝑡,𝑥=𝑓𝑡,𝑥0,𝑥=0,𝑥∈0,1𝑡,0=𝑢𝑡,1=0,𝑡∈0,𝑇(46) where the right-side function 𝑓(𝑡,𝑥)=(2𝑡^2−𝛼/Γ(3−𝛼))sin(2𝜋𝑥)+(2𝑡^2−𝛽/Γ(3−𝛽))sin(2𝜋𝑥)+4𝜋²sin(2𝜋𝑥)𝑡². The exact solution is 𝑢(𝑡,𝑥)=𝑡²sin(2𝜋𝑥).
We use this example to check the convergence rate (c. rate) and CPU time (CPUT) of numerical solutions with respect to the fractional orders 𝛼 and 𝛽.
In the first test, we fix 𝑇=1, 𝛼=0.9 and 𝛽=0.5 and choose Δ𝑥=0.001 which is small enough such that the space discretization errors are negligible as compared with the time errors. Choosing Δ𝑡=1/2^𝑖 (𝑖=2,4,…,7), we report that the convergence rate of FDM in time is nearly 1.15 in Table 1, which matches well with the result of Theorem 4. On the other hand, Table 2 shows that an approximate convergence rate is 2, by fixing Δ𝑡=0.001 and choosing Δ𝑥=1/2^𝑖 (𝑖=2,…,6), which matches well with the result of Theorem 6. In the second test, we give the convergence rate when 𝛼=0.5, 𝛽=0.25 for Δ𝑡 in Table 3, and Δ𝑥 in Table 4, respectively. We also report the 𝐿₂-norm and 𝐻¹-norm of errors in Figures 1 and 2, respectively.
Fixing Δ𝑥=0.001, 𝛼=0.9, and 𝛽=0.3 in (46), we compare the error and CPUT calculated by the FEM in this paper with the FDM in [8] and the FPCM in [8]. From Table 5, it can be seen that the FEM in this paper is computationally effective.

Table 1 

Convergence rate in time for (46) with 𝛼=0.9 and 𝛽=0.5.

Table 2 

Convergence rate in space for (46) with 𝛼=0.9 and 𝛽=0.5.

Table 3 

Convergence rate in time for (46) with 𝛼=0.5 and 𝛽=0.25.

Table 4 

Convergence rate in space for (46) with 𝛼=0.5 and 𝛽=0.25.

Table 5 

Comparison of error and CPUT for (46) with 𝛼=0.9 and 𝛽=0.3.

(a)

(b)

Figure 1 

𝐻1-norm and 𝐿2-norm of errors for (46) with 𝛼=0.9, 𝛽=0.5, Δ𝑥=0.001 (a), and Δ𝑡=0.001 (b).

(a)

(b)

Figure 2 

𝐻1-norm and 𝐿2-norm of errors for (46) with 𝛼=0.5, 𝛽=0.25, Δ𝑥=0.001 (a), and Δ𝑡=0.001 (b).

Example 10. Consider the following multierm fractional differential problem: ^𝐶₀𝐷^𝛼_𝑡𝑢(𝑡)+𝑡₀^−0.2𝐶𝐷^𝛽_𝑡𝑢(𝑡)=𝑓(𝑡),𝑢(0)=0,(47) where 𝑓(𝑡)=(4𝑡^1.5/Γ(2.5))+(12𝑡²/Γ(3.5)). For 𝛼=0.5 and 𝛽=0.3, the exact solution is 𝑢(𝑡)=𝑡²+𝑡^2.5.
For the problem (47), our method in this paper is just the DFBDM in Section 3. Therefore, we only need to compare M1 with the FEM in [28] (FEM2). In Table 6, although the convergence rate of FEM2 is higher than that of DFBDM, the error and CPUT of DFBDM are smaller than those of FEM2.

Table 6 

Comparison of error, convergence rate, and CPUT for (47) with 𝛼=0.5 and 𝛽=0.3.

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.