Numerical Analysis of an 𝐻¹-Galerkin Mixed Finite Element Method for Time Fractional Telegraph Equation

Abstract

We discuss and analyze an 𝐻¹-Galerkin mixed finite element (𝐻¹-GMFE) method to look for the numerical solution of time fractional telegraph equation. We introduce an auxiliary variable to reduce the original equation into lower-order coupled equations and then formulate an 𝐻¹-GMFE scheme with two important variables. We discretize the Caputo time fractional derivatives using the finite difference methods and approximate the spatial direction by applying the 𝐻¹-GMFE method. Based on the discussion on the theoretical error analysis in 𝐿²-norm for the scalar unknown and its gradient in one dimensional case, we obtain the optimal order of convergence in space-time direction. Further, we also derive the optimal error results for the scalar unknown in 𝐻¹-norm. Moreover, we derive and analyze the stability of 𝐻¹-GMFE scheme and give the results of a priori error estimates in two- or three-dimensional cases. In order to verify our theoretical analysis, we give some results of numerical calculation by using the Matlab procedure.

1. Introduction

In this paper, our purpose is to present and discuss a mixed finite element method for the time fractional telegraph equation 𝜕^2𝛼_0,𝑡𝑢(𝑥,𝑡)+2𝜅𝜕^𝛼_0,𝑡𝑢()()()(),()𝑥,𝑡−Δ𝑢𝑥,𝑡+𝛽𝑢𝑥,𝑡=𝑓𝑥,𝑡𝑥,𝑡∈Ω×𝐽, (1) with boundary condition 𝑢(𝑥,𝑡)=0,(𝑥,𝑡)∈𝜕Ω×𝐽, (2) and initial conditions 𝑢(𝑥,0)=𝑢₀(𝑥),𝑢_𝑡(𝑥,0)=𝑢₁(𝑥),𝑥∈Ω, (3) where Ω(⊂𝑅^𝑑,𝑑=1,2,3) is a bounded domain with boundary 𝜕Ω and 𝐽=(0,𝑇] is the time interval with 0<𝑇<∞. The coefficients 𝜅>0 and 𝛽≥0 are two constants, 𝑓(𝑥,𝑡) is a given source function, 𝑢₀(𝑥) and 𝑢₁(𝑥) are two given initial functions, and the time Caputo fractional-order derivatives 𝜕^𝛼_0,𝑡𝑢(𝑥,𝑡) and 𝜕^2𝛼_0,𝑡𝑢(𝑥,𝑡) are defined, respectively, by 𝜕^𝛼_0,𝑡𝑢()=1𝑥,𝑡Γ(1−𝛼)^𝑡₀𝜕𝑢(𝑥,𝜏)𝜕𝜏𝑑𝜏(𝑡−𝜏)^𝛼,𝜕^2𝛼_0,𝑡𝑢()=1𝑥,𝑡Γ(2−2𝛼)^𝑡₀𝜕²𝑢(𝑥,𝜏)𝜕𝜏²𝑑𝜏(𝑡−𝜏)^2𝛼−1, (4) where 1/2<𝛼<1.

In the current literatures, we can see that some numerical methods for solving fractional partial differential equations (PDEs), which include finite element methods [1–8], mixed finite element methods [9], finite difference methods [10–24], finite volume methods [25, 26], spectral methods [27], and discontinuous Galerkin methods [28–31], have been considered and analyzed. In 2014, Liu et al. [9] gave some theoretical error analysis for a class of fractional PDE based on a nonstandard mixed method in spatial direction and a finite difference scheme in time direction. In [32], Zhao and Li discussed finite element method for the fractional telegraph equation. In 2014, Wei et al. [31] studied the numerical solution for time fractional telegraph equation based on the LDG method. But, we have not seen any related studies on mixed finite element methods for solving the fractional telegraph equation.

Recently, some people have made use of the method to obtain the numerical solution for some partial differential equations since Pani (in 1998) [33] proposed an 𝐻¹-GMFE method. This method includes some advantages, such as avoiding the LBB consistency condition, allowing different polynomial degrees of the finite element spaces, and obtaining the optimal a priori estimates in both 𝐻¹ and 𝐿²-norms. In [34], Pani and Fairweather discussed some detailed a priori error results of two numerical schemes based on the 𝐻¹-GMFE method for linear parabolic integrodifferential equations. In [35], Pani et al. gave some error analysis based on 𝐻¹-GMFE scheme for the partial differential equation of hyperbolic type. At the same time, a modified 𝐻¹-GMFE procedure was also proposed and analyzed for the case of two or three dimensions. Guo and Chen [36, 37] obtained some theoretical error analysis and numerical results of 𝐻¹-GMFE method for RLW equation and Sobolev equation. In 2013, Liu et al. [38] introduced another auxiliary variable, which is different from the one in [36], and then proposed and studied an explicit multistep 𝐻¹-GMFE scheme for a RLW equation. Liu and Li [39] and Zhou [40] gave some different discussions on 𝐻¹-GMFE method for pseudohyperbolic equations (heat transport equation), respectively. In [41], Che et al. studied the 𝐻¹-GMFE method for a nonlinear integrodifferential equation. Recently, Shi et al. [42, 43] proposed some nonconforming mixed scheme based on the 𝐻¹-GMFE method.

Based on the above review on 𝐻¹-GMFE method, we easily see that the method was studied based on the integer-order partial differential equations. However, the theoretical results of 𝐻¹-GMFE method for fractional telegraph equation have not been presented and analyzed.

In this paper, our aim is to give some detailed a priori error analysis and numerical results on the 𝐻¹-GMFE method for time fractional telegraph equation. We apply the difference schemes to approximate the time fractional derivatives and use the 𝐻¹-GMFE method to discretize spatial direction. We obtain some optimal a priori error results of one dimension for the scalar unknown in 𝐿² and 𝐻¹-norms. Moreover, we also get an a priori error result of the optimal 𝐿²-norm for the auxiliary variable. At the same time, we use the 𝐻¹-GMFE method to deal with the cases in several dimensions and analyze the stable results for the 𝐻¹-GMFE scheme. We calculate some numerical results to verify the theoretical analysis of 𝐻¹-GMFE method for fractional telegraph equation.

The layout of the paper is as follows. In Section 2, we formulate an 𝐻¹-GMFE scheme for time fractional telegraph equation (5). In Section 3, we introduce some lemmas of two important projections and two difference approximations for time fractional derivatives and then analyze some a priori error results. In Section 4, some theoretical results are given in the cases of two and three dimensions. In Section 5, we choose a numerical example to verify the theoretical analysis of our method. In Section 6, we give some remarks and extensions about the 𝐻¹-GMFE method for fractional PDEs.

For the need of study, we denote the natural inner product as (⋅,⋅) in (𝐿²(Ω))^𝑑, 𝑑=1,2,3. Further, we write the classical Sobolev spaces 𝑊^𝑚,2(Ω) as 𝐻^𝑚 with norm ‖𝑧‖_𝑚=[∑_{0≤|𝛽|≤𝑚}∫_Ω|𝐷^𝛽𝑧|²𝑑𝑥]^1/2. When 𝑚=0, we simply write the norm ‖⋅‖₀ as ‖⋅‖.

2. An 𝐻¹-GMFE Scheme in One Space Dimension

In this section, we first consider the 𝐻¹-GMFE method for the following time fractional telegraph equation in 1D case: 𝜕^2𝛼_0,𝑡𝑢(𝑥,𝑡)+2𝜅𝜕^𝛼_0,𝑡𝑢()−𝜕𝑥,𝑡²𝑢(𝑥,𝑡)𝜕𝑥²()(),()+𝛽𝑢𝑥,𝑡=𝑓𝑥,𝑡𝑥,𝑡∈Ω×𝐽, (5) with boundary condition 𝑢𝑥_𝐿,𝑡=𝑢𝑥_𝑅,𝑡=0,𝑡∈𝐽, (6) and initial conditions 𝑢(𝑥,0)=𝑢₀(𝑥),𝑢_𝑡(𝑥,0)=𝑢₁(𝑥),𝑥∈Ω, (7) where Ω=[𝑥_𝐿,𝑥_𝑅]⊂𝑅.

In order to get the 𝐻¹-GMFE formulation, we first introduce an auxiliary variable 𝜎=𝜕𝑢(𝑥,𝑡)/𝜕𝑥 and split (5) into the following first-order system by 𝜕^2𝛼_0,𝑡𝑢(𝑥,𝑡)+2𝜅𝜕^𝛼_0,𝑡𝑢()−()𝑥,𝑡𝜕𝜎𝑥,𝑡()(),𝜕𝑥+𝛽𝑢𝑥,𝑡=𝑓𝑥,𝑡()𝜎−𝜕𝑢𝑥,𝑡𝜕𝑥=0. (8)(9)

Multiplying (8) by −𝜕𝑤/𝜕𝑥, 𝑤∈𝐻¹, integrating with respect to space from 𝑥_𝐿 to 𝑥_𝑅, and using an integration by parts with 𝜕𝑢(𝑥_𝐿,𝑡)/𝜕𝑡=𝜕𝑢(𝑥_𝑅,𝑡)/𝜕𝑡=0 and 𝜕²𝑢(𝑥_𝐿,𝑡)/𝜕𝑡²=𝜕²𝑢(𝑥_𝑅,𝑡)/𝜕𝑡²=0, we easily get 𝜕^2𝛼_0,𝑡𝜎,𝑤+2𝜅𝜕^𝛼_0,𝑡𝜎,𝑤+𝜕𝜎,𝜕𝑥𝜕𝑤𝜕𝑥=𝛽𝑢,𝜕𝑤−𝜕𝑥𝑓,𝜕𝑤.𝜕𝑥 (10) Multiply (9) by 𝜕𝑣/𝜕𝑥, 𝑣∈𝐻¹₀, and integrate with respect to space from 𝑥_𝐿 to 𝑥_𝑅 to obtain 𝜕𝑢,𝜕𝑥𝜕𝑣=𝜕𝑥𝜎,𝜕𝑣𝜕𝑥,∀𝑣∈𝐻¹₀. (11) For formulating finite element scheme, we now choose the finite element spaces 𝑉_ℎ⊂𝐻¹₀ and 𝑊_ℎ⊂𝐻¹, which satisfy the following approximation properties: for 1≤𝑝≤∞ and 𝑘, 𝑟 positive integers [33], inf_{𝑣ℎ∈𝑉ℎ}𝑣−𝑣_ℎ𝐿𝑝+ℎ𝑣−𝑣_{ℎ𝑊1,𝑝}≤𝐶ℎ^𝑘+1‖𝑣‖_𝑊^𝑘+1,𝑝,𝑣∈𝐻¹₀∩𝑊^𝑘+1,𝑝,inf_{𝑤ℎ∈𝑊ℎ}𝑤−𝑤_ℎ𝐿𝑝+ℎ𝑤−𝑤_{ℎ𝑊1,𝑝}≤𝐶ℎ^𝑟+1‖𝑤‖_𝑊^𝑟+1,𝑝,𝑤∈𝑊^𝑟+1,𝑝. (12) Based on the chosen finite element spaces, the semidiscrete 𝐻¹-GMFE scheme is described by 𝜕𝑢_ℎ,𝜕𝑥𝜕𝑣_ℎ=𝜎𝜕𝑥^ℎ,𝜕𝑣_ℎ𝜕𝑥,∀𝑣_ℎ∈𝑉_ℎ,𝜕^2𝛼_0,𝑡𝜎_ℎ,𝑤_ℎ+2𝜅𝜕^𝛼_0,𝑡𝜎_ℎ,𝑤_ℎ+𝜕𝜎_ℎ,𝜕𝑥𝜕𝑤_ℎ𝑢𝜕𝑥=𝛽_ℎ,𝜕𝑤_ℎ−𝜕𝑥𝑓,𝜕𝑤_ℎ𝜕𝑥,∀𝑤_ℎ∈𝑊_ℎ. (13)

3. Full Discrete Scheme and A Priori Error Estimates

3.1. Two Projection Lemmas

For a priori error estimates for fully discrete scheme, we introduce two projection operators [33, 44] in Lemmas 1 and 2.

Lemma 1. One defines a Ritz projection 𝑃_ℎ𝑢∈𝑉_ℎ for the variable 𝑢 by 𝑢_𝑥−𝑃_ℎ𝑢_𝑥,𝑣_ℎ𝑥=0,𝑣_ℎ∈𝑉_ℎ. (14) Then the following estimates hold, for 𝑗=0,1: 𝑢−𝑃_ℎ𝑢_𝑗≤𝐶_⋆ℎ^{𝑘+1−𝑗}‖𝑢‖_𝑘+1. (15)

Lemma 2. Further, one also defines an elliptic projection 𝑅_ℎ𝜎∈𝑊_ℎ of 𝜎 as the solution of 𝔅𝜎−𝑅_ℎ𝜎,𝑤_ℎ=0,𝑤_ℎ∈𝑊_ℎ, (16) where 𝔅(𝜎,𝑤)=(𝜎_𝑥,𝑤_𝑥)+𝜆(𝜎,𝑤). Here 𝜆>0 is chosen to satisfy 𝔅(𝑤,𝑤)≥𝜇₀‖𝑤‖²₁,𝑤∈𝐻¹,𝜇₀>0. (17) Then the following estimates are found: for 𝑗=0,1, 𝜎−𝑅_ℎ𝜎_𝑗≤𝐶_∗ℎ^{𝑟+1−𝑗}‖𝜎‖_𝑟+1. (18)

3.2. Approximation of Time-Fractional Derivative

For formulating fully discrete scheme, let 0=𝑡₀<𝑡₁<𝑡₂<⋯<𝑡_𝑀=𝑇 be a given partition of the time interval [0,𝑇] with step length Δ𝑡=𝑇/𝑀 and nodes 𝑡_𝑛=𝑛Δ𝑡, for some positive integer 𝑀. For a smooth function 𝜙 on [0,𝑇], define 𝜙^𝑛=𝜙(𝑡_𝑛). In the following analysis, for deriving the convenience of theoretical process, we now denote 𝐵^𝑗𝛼_𝑛−𝑘=(𝑛−𝑘+1)^{𝑗−𝑗𝛼}−(𝑛−𝑘)^{𝑗−𝑗𝛼}𝐷,𝑗=1,2_𝑡𝜎^𝑘=𝜎^𝑘−𝜎^𝑘−1Δ𝑡,𝐷_𝑡𝑡𝜎^𝑘−1=𝜎^𝑘−2𝜎^𝑘−1+𝜎^𝑘−2Δ𝑡². (19) Now, we will introduce two lemmas on the approximations of time fractional derivatives.

Lemma 3 (see [27]). The time fractional derivative 𝜕^𝛼_0,𝑡𝜎(𝑥,𝑡) at 𝑡=𝑡_𝑛 is approximated by, for 0<𝛼<1, 𝜕^𝛼_0,𝑡𝜎𝑥,𝑡_𝑛=Δ𝑡^1−𝛼Γ(2−𝛼)^𝑛_𝑘=1𝐵^𝛼_𝑛−𝑘𝐷_𝑡𝜎^𝑘+𝐸^𝑛_𝛼, (20) and then holds 𝐸^𝑛_𝛼≤𝐶₀Δ𝑡^2−𝛼. (21)

Lemma 4 (see [1]). The time fractional order derivative 𝜕^2𝛼_0,𝑡𝜎(𝑥,𝑡) at 𝑡=𝑡_𝑛 is estimated by, for 1<2𝛼<2, 𝜕^2𝛼_0,𝑡𝜎𝑥,𝑡_𝑛=Δ𝑡^2−2𝛼Γ(3−2𝛼)^𝑛_𝑘=1𝐵^2𝛼_𝑛−𝑘𝐷_𝑡𝑡𝜎^𝑘−1+𝐸^𝑛_2𝛼, (22) and then holds 𝐸^𝑛_2𝛼≤𝐶₀Δ𝑡. (23)

In the next analysis, we will derive and prove some a priori error results for 𝑢 and 𝜎.

3.3. Error Estimates for Fully Discrete Scheme

Based on the approximation formulas (20) and (22) of time-fractional derivatives, we obtain the time semidiscrete scheme of (10) and (11): 𝜕𝑢^𝑛,𝜕𝑥𝜕𝑣=𝜎𝜕𝑥^𝑛,𝜕𝑣𝜕𝑥,∀𝑣∈𝐻¹₀,Δ𝑡^2−2𝛼Γ(3−2𝛼)^𝑛_𝑘=1𝐵^2𝛼_𝑛−𝑘𝐷_𝑡𝑡𝜎^𝑘−1+,𝑤2𝜅Δ𝑡^1−𝛼Γ(2−𝛼)^𝑛_𝑘=1𝐵^𝛼_𝑛−𝑘𝐷_𝑡𝜎^𝑘,𝑤+𝜕𝜎^𝑛,𝜕𝑥𝜕𝑤𝑢𝜕𝑥=𝛽^𝑛,𝜕𝑤−𝑓𝜕𝑥^𝑛,𝜕𝑤+𝐸𝜕𝑥^𝑛0,𝑤,∀𝑤∈𝐻¹, (24) where 𝐸^𝑛₀=𝐸^𝑛_𝛼+𝐸^𝑛_2𝛼.

Now, we formulate a fully discrete procedure: find (𝑢^𝑛_ℎ,𝜎^𝑛_ℎ)∈𝑉_ℎ×𝑊_ℎ (𝑛=1,…,𝑀−1) such that 𝜕𝑢^𝑛_ℎ,𝜕𝑥𝜕𝑣_ℎ=𝜎𝜕𝑥^𝑛ℎ,𝜕𝑣_ℎ𝜕𝑥,∀𝑣_ℎ∈𝑉_ℎ,Δ𝑡^2−2𝛼Γ(3−2𝛼)^𝑛_𝑘=1𝐵^2𝛼_𝑛−𝑘𝐷_𝑡𝑡𝜎_ℎ^𝑘−1,𝑤_ℎ+2𝜅Δ𝑡^1−𝛼Γ(2−𝛼)^𝑛_𝑘=1𝐵^𝛼_𝑛−𝑘𝐷_𝑡𝜎^𝑘_ℎ,𝑤_ℎ+𝜕𝜎^𝑛_ℎ,𝜕𝑥𝜕𝑤_ℎ𝑢𝜕𝑥=𝛽^𝑛_ℎ,𝜕𝑤_ℎ−𝑓𝜕𝑥^𝑛,𝜕𝑤_ℎ𝜕𝑥,∀𝑤_ℎ∈𝑊_ℎ. (25)

Making a combination of (24)-(25) with two projections (14) and (16), we get the following error equations: 𝜕𝜗^𝑛,𝜕𝑥𝜕𝑣_ℎ=𝛿𝜕𝑥^𝑛,𝜕𝑣_ℎ+𝜚𝜕𝑥^𝑛,𝜕𝑣_ℎ𝜕𝑥,∀𝑣_ℎ∈𝑉_ℎ,Δ𝑡^2−2𝛼Γ(3−2𝛼)^𝑛_𝑘=1𝐵^2𝛼_𝑛−𝑘𝐷_𝑡𝑡𝛿^𝑘−1,𝑤_ℎ+2𝜅Δ𝑡^1−𝛼Γ(2−𝛼)^𝑛_𝑘=1𝐵^𝛼_𝑛−𝑘𝐷_𝑡𝛿^𝑘,𝑤_ℎ+𝔅𝛿^𝑛,𝑤_ℎ=−Δ𝑡^2−2𝛼Γ(3−2𝛼)^𝑛_𝑘=1𝐵^2𝛼_𝑛−𝑘𝐷_𝑡𝑡𝜚^𝑘−1,𝑤_ℎ−2𝜅Δ𝑡^1−𝛼Γ(2−𝛼)^𝑛_𝑘=1𝐵^𝛼_𝑛−𝑘𝐷_𝑡𝜚^𝑘,𝑤_ℎ+𝜆𝜚^𝑛,𝑤_ℎ+𝐸^𝑛₀,𝑤_ℎ+𝛽𝜙^𝑛+𝜗^𝑛,𝜕𝑤_ℎ,𝜕𝑥∀𝑤_ℎ∈𝑊_ℎ, (26)(27) where 𝜙^𝑛+𝜗^𝑛≐𝑢𝑡_𝑛−𝑃_ℎ𝑢^𝑛+𝑃_ℎ𝑢^𝑛−𝑢^𝑛_ℎ=𝑢𝑡_𝑛−𝑢^𝑛_ℎ;𝜚^𝑛+𝛿^𝑛≐𝜎𝑡_𝑛−𝑅_ℎ𝜎^𝑛+𝑅_ℎ𝜎^𝑛−𝜎^𝑛_ℎ=𝜎𝑡_𝑛−𝜎^𝑛_ℎ. (28)

In the following discussion, we will derive the proof for the fully discrete a priori error estimates.

Theorem 5. With 𝜎₀=𝑢_0𝑥 and 𝜎_𝑡(0)=𝑢_1𝑥, suppose that 𝑢^𝑛∈𝐻¹₀∩𝐻^𝑘+1, 𝜎^𝑛∈𝐻^𝑟+1, 𝜎⁰_ℎ=𝑅_ℎ𝜎(0), and 𝜎_ℎ𝑡(0)=𝐿_ℎ𝜎_𝑡(0), where 𝐿_ℎ is the 𝐿² projection defined by (𝑤−𝐿_ℎ𝑤,𝑤_ℎ)=0, 𝑤_ℎ∈𝑊_ℎ. Then there exists a positive constant 𝐶(𝑢,𝜎,𝛼,𝛽) free of space-time discrete parameters ℎ and Δ𝑡 such that, for 1/2<𝛼<1, 𝜎^𝑛−𝜎^𝑛_ℎ≤𝐶𝑢,𝜎,𝛼,𝛽𝑇^2𝛼−12−2𝛼Δ𝑡^1−2𝛼ℎ^{min{𝑘+1,𝑟+1}}+Δ𝑡^3−2𝛼+𝐶(𝜎)ℎ^𝑟+1,𝑢^𝑛−𝑢^𝑛_ℎ𝑗≤𝐶𝑢,𝜎,𝛼,𝛽𝑇^2𝛼−12−2𝛼Δ𝑡^1−2𝛼ℎ^{min{𝑘+1,𝑟+1}}+Δ𝑡^3−2𝛼+𝐶(𝜎)ℎ^𝑟+1+𝐶(𝑢)ℎ^{𝑘+1−𝑗},𝑗=0,1, (29)(30) and for 𝛼→1𝜎^𝑛−𝜎^𝑛_ℎ≤𝐶𝑢,𝜎,𝛽Δ𝑡⁻¹ℎ^{min{𝑘+1,𝑟+1}}(𝜎)ℎ+Δ𝑡+𝐶^𝑟+1,𝑢^𝑛−𝑢^𝑛_ℎ𝑗≤𝐶𝑢,𝜎,𝛽Δ𝑡⁻¹ℎ^{min{𝑘+1,𝑟+1}}(𝜎)ℎ+Δ𝑡+𝐶^𝑟+1+𝐶(𝑢)ℎ^{𝑘+1−𝑗},𝑗=0,1. (31)(32)

Proof. For the need of error analysis, we first consider the 𝐿²-norm ‖𝜗^𝑛‖ and the 𝐻¹-norm ‖𝜗^𝑛‖₁. Taking 𝑣_ℎ=𝜗^𝑛 in (26) and using Poincaré inequality based on the 𝐻¹₀-space and Cauchy-Schwarz inequality, we easily get 𝜗^𝑛≤𝜗^𝑛₁≤𝐶𝜕𝜗^𝑛𝜚𝜕𝑥≤𝐶^𝑛+𝛿^𝑛. (33) In the next analysis, we will give the estimates of ‖𝛿^𝑛‖ and ‖𝜗^𝑛‖ in 𝐿²-norm. Noting that ^𝑛_𝑘=1𝐵^𝛼_𝑛−𝑘𝐷_𝑡𝛿^𝑘=^𝑛−1_𝑘=0𝐵^𝛼_𝑘𝐷_𝑡𝛿^𝑛−𝑘,_𝑛𝑘=1𝐵^2𝛼_𝑛−𝑘𝐷_𝑡𝑡𝛿^𝑘−1=^𝑛−1_𝑘=0𝐵_𝑘^2𝛼𝐷_𝑡𝑡𝛿^{𝑛−𝑘−1}, (34) then (27) may be rewritten as Δ𝑡^2−2𝛼Γ(3−2𝛼)^𝑛−1_𝑘=0𝐵_𝑘^2𝛼𝐷_𝑡𝑡𝛿^{𝑛−𝑘−1},𝑤_ℎ+2𝜅Δ𝑡^1−𝛼Γ(2−𝛼)^𝑛−1_𝑘=0𝐵^𝛼_𝑘𝐷_𝑡𝛿^𝑛−𝑘,𝑤_ℎ+𝔅𝛿^𝑛,𝑤_ℎ=−Δ𝑡^2−2𝛼Γ(3−2𝛼)^𝑛_𝑘=1𝐵^2𝛼_𝑛−𝑘𝐷_𝑡𝑡𝜚^𝑘−1,𝑤_ℎ−2𝜅Δ𝑡^1−𝛼Γ(2−𝛼)^𝑛_𝑘=1𝐵^𝛼_𝑛−𝑘𝐷_𝑡𝜚^𝑘,𝑤_ℎ+𝜆𝜚^𝑛,𝑤_ℎ+𝐸^𝑛₀,𝑤_ℎ+𝛽𝜙^𝑛+𝜗^𝑛,𝜕𝑤_ℎ.𝜕𝑥 (35) We take 𝑤_ℎ=𝛿^𝑛 in (35) and multiply by Γ_𝛼≐Γ(3−2𝛼)Γ(2−𝛼)Δ𝑡^2𝛼, (36) to arrive at Γ(2−𝛼)^𝑛−1_𝑘=0𝐵_𝑘^2𝛼𝛿^𝑛−𝑘−2𝛿^{𝑛−𝑘−1}+𝛿^{𝑛−𝑘−2},𝛿^𝑛+2𝜅Γ(3−2𝛼)Δ𝑡^𝛼𝑛−1_𝑘=0𝐵^𝛼_𝑘𝛿^𝑛−𝑘−𝛿^{𝑛−𝑘−1},𝛿^𝑛+Γ_𝛼𝔅𝛿^𝑛,𝛿^𝑛=−Γ(2−𝛼)^𝑛−1_𝑘=0𝐵_𝑘^2𝛼𝜚^𝑛−𝑘−2𝜚^{𝑛−𝑘−1}+𝜚^{𝑛−𝑘−2},𝛿^𝑛−2𝜅Γ(3−2𝛼)Δ𝑡^𝛼𝑛−1_𝑘=0𝐵^𝛼_𝑘𝜚^𝑛−𝑘−𝜚^{𝑛−𝑘−1},𝛿^𝑛+𝜆Γ_𝛼𝜚^𝑛,𝛿^𝑛+Γ_𝛼𝐸^𝑛₀,𝛿^𝑛+𝛽Γ_𝛼𝜙^𝑛+𝜗^𝑛,𝜕𝜗^𝑛.𝜕𝑥 (37) By the simple calculation, we get the following equality: Γ(2−𝛼)^𝑛−1_𝑘=0𝐵_𝑘^2𝛼𝛿^𝑛−𝑘−2𝛿^{𝑛−𝑘−1}+𝛿^{𝑛−𝑘−2},𝛿^𝑛+2𝜅Γ(3−2𝛼)Δ𝑡^𝛼𝑛−1_𝑘=0𝐵^𝛼_𝑘𝛿^𝑛−𝑘−𝛿^{𝑛−𝑘−1},𝛿^𝑛=Γ(2−𝛼)+2𝜅Γ(3−2𝛼)Δ𝑡^𝛼𝛿^𝑛2()×𝐵+Γ2−𝛼_12𝛼−2𝐵₀^2𝛼𝛿^𝑛−1,𝛿^𝑛−𝐵_𝑛^2𝛼𝛿⁰,𝛿^𝑛+𝐵^2𝛼_𝑛−1𝛿⁻¹,𝛿^𝑛+^𝑛−1_𝑘=1𝐵^2𝛼_𝑘+1−2𝐵_𝑘^2𝛼+𝐵^2𝛼_𝑘−1𝛿^{𝑛−𝑘−1},𝛿^𝑛+2𝜅Γ(3−2𝛼)Δ𝑡^𝛼×^𝑛−1_𝑘=1𝐵^𝛼_𝑘−𝐵^𝛼_𝑘−1𝛿^𝑛−𝑘,𝛿^𝑛−𝐵^𝛼_𝑛−1𝛿⁰,𝛿^𝑛. (38) By applying the similar process of calculation to (38), we use Cauchy-Schwarz inequality to get −Γ(2−𝛼)^𝑛−1_𝑘=0𝐵_𝑘^2𝛼𝜚^𝑛−𝑘−2𝜚^{𝑛−𝑘−1}+𝜚^{𝑛−𝑘−2},𝛿^𝑛−2𝜅Γ(3−2𝛼)Δ𝑡^𝛼𝑛−1_𝑘=0𝐵^𝛼_𝑘𝜚^𝑛−𝑘−𝜚^{𝑛−𝑘−1},𝛿^𝑛≤Γ(2−𝛼)+2𝜅Γ(3−2𝛼)Δ𝑡^𝛼𝜚^𝑛()×+Γ2−𝛼2𝐵₀^2𝛼−𝐵₁^2𝛼𝜚^𝑛−1+𝐵_𝑛^2𝛼𝜚⁰+𝐵^2𝛼_𝑛−1𝜚⁻¹+^𝑛−1_𝑘=1𝐵^2𝛼_𝑘+1−2𝐵_𝑘^2𝛼+𝐵^2𝛼_𝑘−1𝜚^{𝑛−𝑘−1}𝛿^𝑛+2𝜅Γ(3−2𝛼)Δ𝑡^𝛼×^𝑛−1_𝑘=1𝐵^𝛼_𝑘−𝐵^𝛼_𝑘−1𝜚^𝑛−𝑘+𝐵^𝛼_𝑛−1𝜚⁰𝛿^𝑛. (39) Noting (33), we use Cauchy-Schwarz inequality and Young inequality to have 𝛽Γ_𝛼𝜙^𝑛+𝜗^𝑛,𝜕𝜗^𝑛𝜕𝑥≤𝛽Γ_𝛼𝜙^𝑛𝜕𝜗^𝑛𝜕𝑥+𝛽Γ_𝛼𝜗^𝑛𝜕𝜗^𝑛𝛽Γ𝜕𝑥≤𝐶_𝛼𝜙^𝑛𝜚^𝑛+𝛿^𝑛+𝐶𝛽Γ_𝛼𝜚^𝑛+𝛿^𝑛2≤𝐶𝛽Γ_𝛼𝜙^𝑛2+𝜚^𝑛2+𝐶𝛽Γ_𝛼𝛿^𝑛2. (40) Substitute (38), (39), and (40) into (37) and use Cauchy-Schwarz inequality to arrive at Γ(2−𝛼)+2𝜅Γ(3−2𝛼)Δ𝑡^𝛼−𝐶𝛽Γ_𝛼𝛿^𝑛2+Γ_𝛼𝔅𝛿^𝑛,𝛿^𝑛≤Γ(2−𝛼)+2𝜅Γ(3−2𝛼)Δ𝑡^𝛼𝜚^𝑛()×+Γ2−𝛼2𝐵₀^2𝛼−𝐵₁^2𝛼𝛿^𝑛−1+𝜚^𝑛−1+𝐵_𝑛^2𝛼𝛿⁰+𝜚⁰+𝐵^2𝛼_𝑛−1𝛿⁻¹+𝜚⁻¹+^𝑛−1_𝑘=1𝐵^2𝛼_𝑘+1−2𝐵_𝑘^2𝛼+𝐵^2𝛼_𝑘−1×𝛿^{𝑛−𝑘−1}+𝜚^{𝑛−𝑘−1}+𝜆Γ_𝛼𝜚^𝑛+Γ_𝛼𝐸^𝑛₀𝛿^𝑛+2𝜅Γ(3−2𝛼)Δ𝑡^𝛼×^𝑛−1_𝑘=1𝐵^𝛼_𝑘−𝐵^𝛼_𝑘−1𝛿^𝑛−𝑘+𝜚^𝑛−𝑘+𝐵^𝛼_𝑛−1𝛿⁰+𝜚⁰𝛿^𝑛+𝐶𝛽Γ_𝛼𝜙^𝑛2+𝜚^𝑛2. (41) By an application of Young inequality, we get Γ(2−𝛼)𝛿^𝑛+Γ_𝛼𝔅𝛿^𝑛,𝛿^𝑛≤1Γ()×Γ()()2−𝛼2−𝛼+2𝜅Γ3−2𝛼Δ𝑡^𝛼𝜚^𝑛()×+Γ2−𝛼2𝐵₀^2𝛼−𝐵₁^2𝛼𝛿^𝑛−1+𝜚^𝑛−1+𝐵_𝑛^2𝛼𝛿⁰+𝜚⁰+𝐵^2𝛼_𝑛−1𝛿⁻¹+𝜚⁻¹+^𝑛−1_𝑘=1𝐵^2𝛼_𝑘+1−2𝐵_𝑘^2𝛼+𝐵^2𝛼_𝑘−1×𝛿^{𝑛−𝑘−1}+𝜚^{𝑛−𝑘−1}+𝜆Γ_𝛼𝜚^𝑛+Γ_𝛼𝐸^𝑛₀+2𝜅Γ(3−2𝛼)Δ𝑡^𝛼Γ()×2−𝛼_{𝑛−1𝑘=1}𝐵^𝛼_𝑘−𝐵^𝛼_𝑘−1𝛿^𝑛−𝑘+𝜚^𝑛−𝑘+𝐵^𝛼_𝑛−1𝛿⁰+𝜚⁰+𝐶𝛽Γ_𝛼𝜙^𝑛+𝜚^𝑛. (42)
Noting that 𝐵_𝑘^2𝛼−𝐵^2𝛼_𝑘+1>0, we arrive at ^𝑛−1_𝑘=1𝐵^2𝛼_𝑘+1−2𝐵_𝑘^2𝛼+𝐵^2𝛼_𝑘−1≤^𝑛−1_𝑘=1𝐵_𝑘^2𝛼−𝐵^2𝛼_𝑘+1+^𝑛−1_𝑘=1𝐵^2𝛼_𝑘−1−𝐵_𝑘^2𝛼=𝐵₁^2𝛼−𝐵_𝑛^2𝛼+𝐵₀^2𝛼−𝐵^2𝛼_𝑛−1. (43) At the same time, noting that 𝐵^𝛼_𝑘−𝐵^𝛼_𝑘+1>0, we use the similar method to have ^𝑛−1_𝑘=1𝐵^𝛼_𝑘−𝐵^𝛼_𝑘−1=𝐵^𝛼₀−𝐵^𝛼_𝑛−1. (44) By a combination of (43) and (44) with (18) and noting that 𝐵^𝛼₀=𝐵₀^2𝛼=1, we arrive at 1Γ()Γ()()2−𝛼2−𝛼+2𝜅Γ3−2𝛼Δ𝑡^𝛼𝜚^𝑛()×+Γ2−𝛼2𝐵₀^2𝛼−𝐵₁^2𝛼𝜚^𝑛−1+𝐵_𝑛^2𝛼𝜚⁰+𝐵^2𝛼_𝑛−1𝜚⁻¹+^𝑛−1_𝑘=1𝐵^2𝛼_𝑘+1−2𝐵_𝑘^2𝛼+𝐵^2𝛼_𝑘−1𝜚^{𝑛−𝑘−1}+2𝜅Γ(3−2𝛼)Δ𝑡^𝛼×^𝑛−1_𝑘=1𝐵^𝛼_𝑘−𝐵^𝛼_𝑘−1𝜚^𝑛−𝑘+𝐵^𝛼_𝑛−1𝜚⁰+𝜆Γ_𝛼𝜚^𝑛+Γ_𝛼𝐸^𝑛₀+𝐶𝛽Γ_𝛼𝜙^𝑛+𝜚^𝑛≤𝐶Γ()Γ()()2−𝛼2−𝛼+2𝜅Γ3−2𝛼Δ𝑡^𝛼()×+Γ2−𝛼2𝐵₀^2𝛼−𝐵₁^2𝛼+𝐵_𝑛^2𝛼+𝐵^2𝛼_𝑛−1+𝐵₁^2𝛼−𝐵_𝑛^2𝛼+𝐵₀^2𝛼−𝐵^2𝛼_𝑛−1+2𝜅Γ(3−2𝛼)Δ𝑡^𝛼×𝐵^𝛼₀−𝐵^𝛼_𝑛−1+𝐵^𝛼_𝑛−1+𝜆Γ_𝛼ℎ^𝑟+1+𝐶Γ(3−2𝛼)Γ(2−𝛼)Δ𝑡^1+2𝛼=𝐶Γ()Γ()()2−𝛼2−𝛼+2𝜅Γ3−2𝛼Δ𝑡^𝛼ℎ^𝑟+1+𝐶𝑢,𝜎,𝛽Γ_𝛼ℎ^𝑟+1+ℎ^𝑘+1+𝐶Γ(3−2𝛼)Γ(2−𝛼)Δ𝑡^1+2𝛼. (45) Substitute (45) into (42) and note that 𝜇₀‖𝛿^𝑛‖₁≤𝔅(𝛿^𝑛,𝛿^𝑛)/‖𝛿^𝑛‖₁≤𝔅(𝛿^𝑛,𝛿^𝑛)/‖𝛿^𝑛‖ to get Γ(2−𝛼)𝛿^𝑛+𝜇₀Γ_𝛼𝛿^𝑛₁≤Γ(2−𝛼)2𝐵₀^2𝛼−𝐵₁^2𝛼𝛿^𝑛−1+𝐵_𝑛^2𝛼𝛿⁰+𝐵^2𝛼_𝑛−1𝛿⁻¹+^𝑛−1_𝑘=1𝐵^2𝛼_𝑘+1−2𝐵_𝑘^2𝛼+𝐵^2𝛼_𝑘−1𝛿^{𝑛−𝑘−1}+2𝜅Γ(3−2𝛼)Δ𝑡^𝛼Γ()×2−𝛼_{𝑛−1𝑘=1}𝐵^𝛼_𝑘−𝐵^𝛼_𝑘−1𝛿^𝑛−𝑘+𝐵^𝛼_𝑛−1𝛿⁰+𝐶Γ()Γ()()2−𝛼2−𝛼+2𝜅Γ3−2𝛼Δ𝑡^𝛼ℎ^𝑟+1+𝐶𝑢,𝜎,𝛽Γ_𝛼ℎ^𝑟+1+ℎ^𝑘+1+𝐶Γ(3−2𝛼)Γ(2−𝛼)Δ𝑡^1+2𝛼. (46) In the following discussion, we apply mathematical induction to obtain the error result 𝛿^𝑛≤𝐶𝑢,𝜎,𝛼,𝛽𝐵^2𝛼_𝑛−1ℎ^𝑘+1+ℎ^𝑟+1+Δ𝑡²+Δ𝑡^1+2𝛼. (47) Take 𝑛=1 in (46) and note that 𝐵₀^2𝛼=1; it is easy to find that the following inequality holds: 𝛿¹+Δ𝑡^𝛼𝛿¹₁≤𝐶𝑢,𝜎,𝛼,𝛽ℎ^𝑘+1+ℎ^𝑟+1+Δ𝑡²+Δ𝑡^1+2𝛼≤𝐶𝑢,𝜎,𝛼,𝛽𝐵₀^2𝛼ℎ^𝑘+1+ℎ^𝑟+1+Δ𝑡²+Δ𝑡^1+2𝛼. (48) Assuming that the inequalities ‖𝛿^𝑗‖≤(𝐶(𝑢,𝜎,𝛼,𝛽)/𝐵^2𝛼_𝑗−1)(ℎ^𝑘+1+ℎ^𝑟+1+Δ𝑡²+Δ𝑡^1+2𝛼) hold for 𝑗=1,2,…,𝑛−1, we now prove that the inequality (47) holds. Noting that 1/𝐵_𝑘^𝑗𝛼<1/𝐵^𝑗𝛼_𝑘+1, 𝑗=1,2, and combining (43) and (44) with (46), we have 𝛿^𝑛+Δ𝑡^𝛼𝛿^𝑛₁≤Γ(2−𝛼)2𝐵₀^2𝛼−𝐵₁^2𝛼𝐶𝑢,𝜎,𝛼,𝛽𝐵^2𝛼_𝑛−1×ℎ^𝑟+1+ℎ^𝑘+1+Δ𝑡²+Δ𝑡^1+2𝛼+𝐵_𝑛^2𝛼𝛿⁰+𝐵^2𝛼_𝑛−1𝛿⁻¹+^𝑛−1_𝑘=1𝐵^2𝛼_𝑘+1−2𝐵_𝑘^2𝛼+𝐵^2𝛼_𝑘−1×𝐶𝑢,𝜎,𝛼,𝛽𝐵^2𝛼_{𝑛−𝑘−2}×ℎ^𝑟+1+ℎ^𝑘+1+Δ𝑡²+Δ𝑡^1+2𝛼+2𝜅Γ(3−2𝛼)Δ𝑡^𝛼Γ()×2−𝛼_{𝑛−1𝑘=1}𝐵^𝛼_𝑘−𝐵^𝛼_𝑘−1𝐶𝑢,𝜎,𝛼,𝛽𝐵^2𝛼_{𝑛−𝑘−1}×ℎ^𝑟+1+ℎ^𝑘+1+Δ𝑡²+Δ𝑡^1+2𝛼+𝐵^𝛼_𝑛−1𝛿⁰+𝐶Γ()Γ()()2−𝛼2−𝛼+2𝜅Γ3−2𝛼Δ𝑡^𝛼ℎ^𝑟+1+𝐶𝑢,𝜎,𝛽Γ_𝛼ℎ^𝑟+1+ℎ^𝑘+1+𝐶Γ(3−2𝛼)Γ(2−𝛼)Δ𝑡^1+2𝛼≤𝐶𝑢,𝜎,𝛼,𝛽𝐵^2𝛼_𝑛−1ℎ^𝑟+1+ℎ^𝑘+1+Δ𝑡²+Δ𝑡^1+2𝛼. (49) Based on the process of mathematical induction, we claim that (47) holds.
For the need of the next proof, we have to estimate the term 𝑛^1−2𝛼/𝐵^2𝛼_𝑛−1. We now use Taylor formula to arrive at 𝑛^1−2𝛼𝐵^2𝛼_𝑛−1=𝑛^1−2𝛼𝑛^2−2𝛼−(𝑛−1)^2−2𝛼=1𝑛−(𝑛−1)(1+1/(𝑛−1))^2𝛼−1=()×()1𝑛−𝑛−11+2𝛼−1+()𝑛−12𝛼−1)(2𝛼−2×12!1+𝜃𝑛−1^2𝛼−31𝑛−1²⁻¹=()2−2𝛼+2𝛼−1)(2−2𝛼×12!1+𝜃𝑛−1^2𝛼−31𝑛−1⁻¹<1,2−2𝛼 (50) where 0<𝜃<1 and 1/2<𝛼<1. Noting that (𝑛Δ𝑡)^2𝛼−1≤𝑇^2𝛼−1 and (50), we have 𝛿^𝑛+Δ𝑡^𝛼𝛿^𝑛₁≤𝐶𝑢,𝜎,𝛼,𝛽𝐵^2𝛼_𝑛−1ℎ^𝑟+1+ℎ^𝑘+1+Δ𝑡²+Δ𝑡^1+2𝛼=𝑛^2𝛼−1Δ𝑡^2𝛼−1Δ𝑡^1−2𝛼𝑛^1−2𝛼𝐶𝑢,𝜎,𝛼,𝛽𝐵^2𝛼_𝑛−1×ℎ^𝑟+1+ℎ^𝑘+1+Δ𝑡²+Δ𝑡^1+2𝛼≤𝑇^2𝛼−1Δ𝑡^1−2𝛼𝐶2−2𝛼𝑢,𝜎,𝛼,𝛽𝐵^2𝛼_𝑛−1×ℎ^𝑟+1+ℎ^𝑘+1+Δ𝑡²+Δ𝑡^1+2𝛼=𝐶𝑢,𝜎,𝛼,𝛽𝑇^2𝛼−1×2−2𝛼Δ𝑡^1−2𝛼ℎ^𝑟+1+ℎ^𝑘+1+Δ𝑡²+Δ𝑡^3−2𝛼. (51) Substitute (51) into (33) and use (18) to get 𝜗^𝑛+𝜗^𝑛₁≤𝐶𝜚^𝑛+𝛿^𝑛≤𝐶𝑢,𝜎,𝛼,𝛽𝑇^2𝛼−1×2−2𝛼Δ𝑡^1−2𝛼ℎ^𝑟+1+ℎ^𝑘+1+Δ𝑡²+Δ𝑡^3−2𝛼+𝐶(𝜎)ℎ^𝑟+1. (52) By combining (51) and (18) with triangle inequality, the 𝐿²-norm estimate ‖𝜎^𝑛−𝜎^𝑛_ℎ‖ is got. Similarly, the estimate ‖𝑢^𝑛−𝑢^𝑛_ℎ‖ in the 𝐿²-norm and the estimate ‖𝑢^𝑛−𝑢^𝑛_ℎ‖₁ in the 𝐻¹-norm also are obtained by a combination of (52) and (15) with triangle inequality.
Now we mainly analyze the case 𝛼→1. We can find that the error inequality (29) in this case has no meaning since the coefficient 𝐶(𝑢,𝜎,𝛼,𝛽)/(2−2𝛼)→∞ as 𝛼→1. So, we have to look for another error estimate’s process. Noting the fact that 𝑛Δ𝑡≤𝑇 (𝑛=1,2,…,𝑀), we can obtain the following error inequality: 𝛿^𝑛≤𝐶𝑢,𝜎,𝛼,𝛽𝑛Δ𝑡Δ𝑡^1−2𝛼ℎ^𝑟+1+ℎ^𝑘+1+Δ𝑡^3−2𝛼. (53)
Now we can use induction to prove the inequality (53). The detailed proof is similar to the above process of analysis, so we do not give the detailed proof again. Making a combination of (53) and (18) with triangle inequality, we can get the estimate (31). A similar discussion for (32) can also be made.

4. An 𝐻¹-GMFE Scheme for Several Spaces Variables

In this section, we consider (1) with two and three spaces variables.

Let 𝐋²(Ω)=(𝐿²(Ω))^𝑑, (𝑑=2 or 3) with inner product and norm (𝐪,𝐰)=^𝑑_𝑖=1𝑞_𝑖,𝑤_𝑖,‖𝐰‖=^𝑑_𝑖=1𝑤_𝑖^21/2. (54) Further, let 𝐖=𝐇(𝐝𝐢𝐯;Ω)=𝐰∈𝐋²(Ω)∣∇⋅𝐰∈𝐿²(Ω) (55) with norm ‖𝐰‖_{𝐇(𝐝𝐢𝐯;Ω)}=‖𝐰‖²+‖∇⋅𝐰‖^21/2. (56) Taking 𝐪=∇𝑢, we use a similar process to the system (10) and (11) to get the 𝐻¹-Galerkin mixed weak formulation {𝑢,𝐪}∶[0,𝑇]↦𝐻¹₀×𝐖 by (∇𝑢,∇𝑣)=(𝐪,∇𝑣),∀𝑣∈𝐻¹₀,𝜕^2𝛼_0,𝑡𝐪,𝐰+2𝜅𝜕^𝛼_0,𝑡+()()−𝐪,𝐰∇⋅𝐪,∇⋅𝐰=𝛽𝑢,∇⋅𝐰𝑓,∇⋅𝐰,∀𝐰∈𝐖. (57) The corresponding time semidiscrete system is defined by ∇𝑢^𝑛,∇𝑣=𝐪^𝑛,∇𝑣,∀𝑣∈𝐻¹₀,Δ𝑡^2−2𝛼Γ(3−2𝛼)^𝑛_𝑘=1𝐵^2𝛼_𝑛−𝑘𝐷_𝑡𝑡𝐪^𝑘−1+,𝐰2𝜅Δ𝑡^1−𝛼Γ(2−𝛼)^𝑛_𝑘=1𝐵^𝛼_𝑛−𝑘𝐷_𝑡𝐪^𝑘,𝐰+∇⋅𝐪^𝑛𝑢,∇⋅𝐰=𝛽^𝑛,∇⋅𝐰−𝑓^𝑛,∇⋅𝐰+𝐑^𝑛₀,𝐰,∀𝐰∈𝐖, (58) where 𝐑^𝑛₀=𝑂(Δ𝑡), which can be estimated by a similar process to 𝐸^𝑛₀.

In order to get fully discrete mixed finite element scheme, we now choose the finite element spaces 𝑉_ℎ⊂𝐻¹₀ and 𝑊_ℎ⊂𝐖, which satisfy the following approximation properties: for 1≤𝑝≤∞ and 𝑘, 𝑟 positive integers [33], inf_{𝑣ℎ∈𝑉ℎ}𝑣−𝑣_ℎ+ℎ𝑣−𝑣_ℎ1≤𝐶ℎ^𝑘+1‖𝑣‖_𝑘+1,𝑣∈𝐻¹₀∩𝐻^𝑘+1,inf_{𝐰ℎ∈𝐖ℎ}𝐰−𝐰_ℎ+ℎ𝐰−𝐰_{ℎ𝐇(𝐝𝐢𝐯;Ω)}≤𝐶ℎ^𝑟+1‖𝐰‖_𝑟+1,𝐰∈𝐇^𝑟+1. (59) The fully discrete 𝐻¹-GMFE scheme is to find {𝑢^𝑛_ℎ,𝐪^𝑛_ℎ}∶[0,𝑇]↦𝑉_ℎ×𝐖_ℎ such that ∇𝑢^𝑛_ℎ,∇𝑣_ℎ=𝐪^𝑛_ℎ,∇𝑣_ℎ,∀𝑣_ℎ∈𝑉_ℎ,Δ𝑡^2−2𝛼Γ(3−2𝛼)^𝑛_𝑘=1𝐵^2𝛼_𝑛−𝑘𝐷_𝑡𝑡𝐪_ℎ^𝑘−1,𝐰_ℎ+2𝜅Δ𝑡^1−𝛼Γ(2−𝛼)^𝑛_𝑘=1𝐵^𝛼_𝑛−𝑘𝐷_𝑡𝐪^𝑘_ℎ,𝐰_ℎ+∇⋅𝐪^𝑛_ℎ,∇⋅𝐰_ℎ=𝛽𝑢^𝑛_ℎ,∇⋅𝐰_ℎ−𝑓^𝑛,∇⋅𝐰_ℎ,∀𝐰∈𝐖_ℎ. (60)(61)

4.1. The Analysis of Stability

In the following discussion, we will give the stability for the system (60) and (61). First, we need to obtain an important lemma on two initial value conditions.

Lemma 6. With 𝐪₀=∇𝑢₀ and 𝐪_𝑡(0)=∇𝑢₁, the following inequality holds: 𝐪_ℎ⁻¹≤𝐶𝑢₀₁+Δ𝑡𝑢₁₁. (62)

Proof. By a similar discussion as in [32], we have 𝐪_ℎ⁻¹≤𝐶𝐪₀+Δ𝑡𝐪_𝑡(0). (63) Based on the given initial value conditions and (63), we arrive at 𝐪_ℎ⁻¹≤𝐶∇𝑢₀+Δ𝑡∇𝑢₁≤𝐶𝑢₀₁+Δ𝑡𝑢₁₁. (64)

Theorem 7. The following stable inequality for the system (60) and (61) holds: 𝑢^𝑛_ℎ1+𝐪^𝑛_ℎ+𝜅Γ_𝛼∇⋅𝐪^𝑛_ℎ≤𝐶(𝛼)Γ_𝛼max_{0≤𝑗≤𝑛}𝑓^𝑗+𝑢₀₁+Δ𝑡𝑢₁₁. (65)

Proof. In (60), we take 𝑣_ℎ=∇𝑢^𝑛_ℎ and use Cauchy-Schwarz inequality and Poincaré inequality to arrive at 𝑢^𝑛_ℎ≤𝑢^𝑛_ℎ1≤𝐶∇𝑢^𝑛_ℎ≤𝐶𝐪^𝑛_ℎ. (66) In (61), we choose 𝐰_ℎ=𝐪^𝑛_ℎ and use Cauchy-Schwarz inequality, Young inequality, and (66) to get Δ𝑡^2−2𝛼Γ(3−2𝛼)^𝑛_𝑘=1𝐵^2𝛼_𝑛−𝑘𝐷_𝑡𝑡𝐪_ℎ^𝑘−1,𝐪^𝑛_ℎ+2𝜅Δ𝑡^1−𝛼Γ(2−𝛼)^𝑛_𝑘=1𝐵^𝛼_𝑛−𝑘𝐷_𝑡𝐪^𝑘_ℎ,𝐪^𝑛_ℎ+∇⋅𝐪^𝑛_ℎ²≤𝛽𝑢^𝑛_ℎ∇⋅𝐪^𝑛_ℎ+𝑓^𝑛∇⋅𝐪^𝑛_ℎ≤𝐶𝐪^𝑛_ℎ∇⋅𝐪^𝑛_ℎ+𝑓^𝑛∇⋅𝐪^𝑛_ℎ≤𝐶𝐪^𝑛_ℎ²+𝑓^𝑛2+12∇⋅𝐪^𝑛_ℎ². (67) By the simple simplification for (67), we easily get Δ𝑡^2−2𝛼Γ(3−2𝛼)^𝑛_𝑘=1𝐵^2𝛼_𝑛−𝑘𝐷_𝑡𝑡𝐪_ℎ^𝑘−1,𝐪^𝑛_ℎ+2𝜅Δ𝑡^1−𝛼Γ(2−𝛼)^𝑛_𝑘=1𝐵^𝛼_𝑛−𝑘𝐷_𝑡𝐪^𝑘_ℎ,𝐪^𝑛_ℎ+12∇⋅𝐪^𝑛_ℎ²≤𝐶𝐪^𝑛_ℎ²+𝑓^𝑛2. (68) For the case 𝑛=1 in (68), we multiply (36) to get easily Γ(2−𝛼)𝐵₀^2𝛼𝐪¹_ℎ−2𝐪⁰_ℎ+𝐪_ℎ⁻¹,𝐪¹_ℎ+2𝜅Δ𝑡^𝛼Γ(3−2𝛼)𝐵^𝛼₀𝐪¹_ℎ−𝐪⁰_ℎ,𝐪¹_ℎ+𝜅Γ_𝛼∇⋅𝐪¹_ℎ²≤𝐶Γ_𝛼𝐪¹_ℎ²+𝑓¹². (69) Noting that 𝐵^𝛼₀=1=𝐵₀^2𝛼 and Lemma 6, we use Cauchy-Schwarz inequality to get 12Γ(2−𝛼)+𝜅Δ𝑡^𝛼Γ(3−2𝛼)−𝐶Γ_𝛼𝐪¹_ℎ²+𝜅Γ_𝛼∇⋅𝐪¹_ℎ²≤𝐶Γ_𝛼𝑓¹²+𝐶(𝛼)𝐪_ℎ⁻¹²+𝐪⁰_ℎ²≤Γ_𝛼𝑓¹²+𝐶(𝛼)𝑢₀²₁+Δ𝑡²𝑢₁²₁≤𝐶(𝛼)Γ_𝛼max_{0≤𝑗≤𝑛}𝑓^𝑗+𝑢₀₁+Δ𝑡𝑢₁₁². (70) By the simple calculation for (70), we arrive at 𝐪¹_ℎ+𝜅Γ_𝛼∇⋅𝐪¹_ℎ≤𝐶(𝛼)Γ_𝛼max_{0≤𝑗≤𝑛}𝑓^𝑗+𝑢₀₁+Δ𝑡𝑢₁₁. (71) For the case 𝑛≥1, using a similar process to the proof of Theorem 5 based on the mathematical induction, we can get 𝐪^𝑛_ℎ+𝜅Γ_𝛼∇⋅𝐪^𝑛_ℎ≤𝐶(𝛼)Γ_𝛼max_{0≤𝑗≤𝑛}𝑓^𝑗+𝑢₀₁+Δ𝑡𝑢₁₁. (72) By a combination of (72) and (66), we get the conclusion of Theorem 7.

4.2. A Priori Error Results

For deriving the a priori error analysis, we define the Ritz projection 𝑢_ℎ∈𝑉_ℎ by ∇𝑢−𝑢_ℎ,∇𝑣_ℎ=0,𝑣_ℎ∈𝑉_ℎ. (73) Further, let 𝐪_ℎ∈𝑊_ℎ be the standard finite element interpolant of 𝐪.

Let 𝜌=𝐪−𝐪_ℎ, 𝜂=𝑢−𝑢_ℎ; see [33, 44]; we obtain 𝜂_𝑗≤𝐶ℎ^{𝑘+1−𝑗}‖𝑢‖_𝑘+1𝜌𝜌,𝑗=0,1,+ℎ_{𝐇(𝐝𝐢𝐯;Ω)}≤𝐶ℎ^𝑟+1‖𝐪‖_𝑟+1. (74) Now, we can get the following theorem of a priori error estimates based on the above contents.

Theorem 8. With 𝐪₀=∇𝑢₀, 𝐪_𝑡(0)=∇𝑢₁, 𝑢^𝑛∈𝐻¹₀∩𝐻^𝑘+1, and 𝐪^𝑛∈𝐇^𝑟+1, there exists a positive constant 𝐶(𝑢,𝐪,𝛼,𝛽) free of space-time meshes ℎ and Δ𝑡 such that, for 1/2<𝛼<1, 𝐪^𝑛−𝐪^𝑛_ℎ≤𝐶𝑢,𝐪,𝛼,𝛽𝑇^2𝛼−12−2𝛼Δ𝑡^1−2𝛼ℎ^{min{𝑘+1,𝑟}}+Δ𝑡^3−2𝛼+𝐶(𝐪)ℎ^𝑟+1,𝑢^𝑛−𝑢^𝑛_ℎ𝑗≤𝐶𝑢,𝐪,𝛼,𝛽𝑇^2𝛼−12−2𝛼Δ𝑡^1−2𝛼ℎ^{min{𝑘+1,𝑟}}+Δ𝑡^3−2𝛼+𝐶(𝐪)ℎ^𝑟+1+𝐶(𝑢)ℎ^{𝑘+1−𝑗},𝑗=0,1, (75) and for 𝛼→1𝐪^𝑛−𝐪^𝑛_ℎ≤𝐶𝑢,𝐪,𝛽Δ𝑡⁻¹ℎ^{min{𝑘+1,𝑟}}+Δ𝑡+𝐶(𝐪)ℎ^𝑟+1,𝑢^𝑛−𝑢^𝑛_ℎ𝑗≤𝐶𝑢,𝐪,𝛽Δ𝑡⁻¹ℎ^{min{𝑘+1,𝑟}}(𝐪)ℎ+Δ𝑡+𝐶^𝑟+1+𝐶(𝑢)ℎ^{𝑘+1−𝑗},𝑗=0,1. (76)

Proof. We can use a similar proof as in Theorem 5 to get the conclusion of Theorem 8, so we do not discuss that again.

5. Some Numerical Results

Here, in order to show the numerical performances on the rate of convergence and a priori error estimates, we now choose 𝜅=1/2 and 𝛽=0, the exact solution 𝑢(𝑥,𝑡)=𝑡^2+𝛼sin(2𝜋𝑥), for all (𝑥,𝑡)∈[0,1]×[0,1], and the determined source term 𝑓(𝑥,𝑡) by the exact solution in (5) and then calculate some numerical results by using Matlab procedure. In the numerical calculation, the time direction is approximated by finite difference schemes and the spatial direction is discretized by the 𝐻¹-GMFE method. Now, we divide interval [0,1] into equal-length intervals 𝑒_𝑖=[𝑥_𝑖,𝑥_𝑖+1], 0≤𝑖≤𝑁−1. Now we choose the piecewise linear spaces 𝑉_ℎ=span{𝜑₁,𝜑₂,…,𝜑_𝑁−1} and 𝑊_ℎ=span{𝜓₀,𝜓₁,…,𝜓_𝑁} with index 𝑘=𝑟=1. Let 𝑢^𝑛_ℎ=^𝑁−1_𝑘=1𝑢_𝑘𝜑_𝑘,𝜎^𝑛_ℎ=^𝑁_𝑘=0𝜎_𝑘𝜓_𝑘. (77) Based on the above expressions (77), we can get the following equivalent algebraic equation of (25): [(𝑎+𝑏)𝐄+𝐅]𝜎^𝑛=𝑎^𝑛−1_𝑘=1−𝐵^2𝛼_{𝑛−𝑘+1}+2𝐵^2𝛼_𝑛−𝑘−𝐵^2𝛼_{𝑛−𝑘−1}𝐄𝜎^𝑘−1+𝑎2𝐵₀^2𝛼−𝐵₁^2𝛼𝐄𝜎^𝑛−1+𝑏^𝑛_𝑘=1−𝐵^𝛼_{𝑛−𝑘+1}+𝐵^𝛼_𝑛−𝑘𝐄𝜎^𝑘−1−𝑎𝐵^2𝛼_𝑛−1𝐄𝜎⁻¹−𝑓^𝑛,𝐆𝑢^𝑛=𝐇𝜎^𝑛, (78) where 𝑎=Δ𝑡^−2𝛼