Numerical analysis of an H1-Galerkin mixed finite element method for time fractional telegraph equation.

We discuss and analyze an H(1)-Galerkin mixed finite element (H(1)-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 H(1)-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 H(1)-GMFE method. Based on the discussion on the theoretical error analysis in L(2)-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 H(1)-norm. Moreover, we derive and analyze the stability of H(1)-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 with boundary condition and initial conditions where Ω(⊂R , d = 1,2, 3) is a bounded domain with boundary ∂Ω and J = (0, T] is the time interval with 0 < T < ∞. The coefficients κ > 0 and β ≥ 0 are two constants, f(x, t) is a given source function, u 0(x) and u 1(x) are two given initial functions, and the time Caputo fractional-order derivatives ∂0, u(x, t) and ∂0, 2 u(x, t) are defined, respectively, by 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 H 1-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 H 1 and L 2-norms. In [34], Pani and Fairweather discussed some detailed a priori error results of two numerical schemes based on the H 1-GMFE method for linear parabolic integrodifferential equations. In [35], Pani et al. gave some error analysis based on H 1-GMFE scheme for the partial differential equation of hyperbolic type. At the same time, a modified H 1-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 H 1-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 H 1-GMFE scheme for a RLW equation. Liu and Li [39] and Zhou [40] gave some different discussions on H 1-GMFE method for pseudohyperbolic equations (heat transport equation), respectively. In [41], Che et al. studied the H 1-GMFE method for a nonlinear integrodifferential equation. Recently, Shi et al. [42, 43] proposed some nonconforming mixed scheme based on the H 1-GMFE method.

Based on the above review on H 1-GMFE method, we easily see that the method was studied based on the integer-order partial differential equations. However, the theoretical results of H 1-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 H 1-GMFE method for time fractional telegraph equation. We apply the difference schemes to approximate the time fractional derivatives and use the H 1-GMFE method to discretize spatial direction. We obtain some optimal a priori error results of one dimension for the scalar unknown in L 2 and H 1-norms. Moreover, we also get an a priori error result of the optimal L 2-norm for the auxiliary variable. At the same time, we use the H 1-GMFE method to deal with the cases in several dimensions and analyze the stable results for the H 1-GMFE scheme. We calculate some numerical results to verify the theoretical analysis of H 1-GMFE method for fractional telegraph equation.

The layout of the paper is as follows. In Section 2, we formulate an H 1-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 H 1-GMFE method for fractional PDEs.

For the need of study, we denote the natural inner product as (·, ·) in (L 2(Ω)), d = 1,2, 3. Further, we write the classical Sobolev spaces W (Ω) as H with norm ||z|| = [∑0≤|∫ | D z|2 dx]1/2. When m = 0, we simply write the norm ||·||0 as ||·||.

2. An H 1-GMFE Scheme in One Space Dimension

In this section, we first consider the H 1-GMFE method for the following time fractional telegraph equation in 1D case: with boundary condition and initial conditions where Ω = [x , x ] ⊂ R.

In order to get the H 1-GMFE formulation, we first introduce an auxiliary variable σ = ∂u(x, t)/∂x and split (5) into the following first-order system by

Multiplying (8) by −∂w/∂x, w ∈ H 1, integrating with respect to space from x to x , and using an integration by parts with ∂u(x , t)/∂t = ∂u(x , t)/∂t = 0 and ∂2 u(x , t)/∂t 2 = ∂2 u(x , t)/∂t 2 = 0, we easily get Multiply (9) by ∂v/∂x, v ∈ H 0 1, and integrate with respect to space from x to x to obtain For formulating finite element scheme, we now choose the finite element spaces V ⊂ H 0 1 and W ⊂ H 1, which satisfy the following approximation properties: for 1 ≤ p ≤ ∞ and k, r positive integers [33], Based on the chosen finite element spaces, the semidiscrete H 1-GMFE scheme is described by

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 P u ∈ V for the variable u by Then the following estimates hold, for j = 0,1:

Lemma 2 .

Further, one also defines an elliptic projection R σ ∈ W of σ as the solution of where B(σ, w) = (σ , w ) + λ(σ, w). Here λ > 0 is chosen to satisfy Then the following estimates are found: for j = 0,1,

3.2. Approximation of Time-Fractional Derivative

For formulating fully discrete scheme, let 0 = t 0 < t 1 < t 2 < ⋯

Lemma 3 (see [27]).

The time fractional derivative ∂0, σ(x, t) at t = t is approximated by, for 0 < α < 1, and then holds

Lemma 4 (see [1]).

The time fractional order derivative ∂0, 2 σ(x, t) at t = t is estimated by, for 1 < 2α < 2, and then holds

In the next analysis, we will derive and prove some a priori error results for u 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): where E 0 = E + E 2 .

Now, we formulate a fully discrete procedure: find (u , σ ) ∈ V × W (n = 1,…, M − 1) such that

Making a combination of (24)-(25) with two projections (14) and (16), we get the following error equations: where

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

Theorem 5 .

With σ 0 = u 0 and σ (0) = u 1, suppose that u ∈ H 0 1∩H , σ ∈ H , σ 0 = R σ(0), and σ (0) = L σ (0), where L is the L 2 projection defined by (w − L w, w ) = 0, w ∈ W . Then there exists a positive constant C(u, σ, α, β) free of space-time discrete parameters h and Δt such that, for 1/2 < α < 1, and for α → 1

Proof

For the need of error analysis, we first consider the L 2-norm ||ϑ || and the H 1-norm ||ϑ ||1. Taking v = ϑ in (26) and using Poincaré inequality based on the H 0 1-space and Cauchy-Schwarz inequality, we easily get In the next analysis, we will give the estimates of ||δ || and ||ϑ || in L 2-norm. Noting that then (27) may be rewritten as We take w = δ in (35) and multiply by to arrive at By the simple calculation, we get the following equality: By applying the similar process of calculation to (38), we use Cauchy-Schwarz inequality to get Noting (33), we use Cauchy-Schwarz inequality and Young inequality to have Substitute (38), (39), and (40) into (37) and use Cauchy-Schwarz inequality to arrive at By an application of Young inequality, we get

Noting that B 2 − B 2 > 0, we arrive at At the same time, noting that B − B > 0, we use the similar method to have By a combination of (43) and (44) with (18) and noting that B 0 = B 0 2 = 1, we arrive at Substitute (45) into (42) and note that μ 0||δ ||1 ≤ B(δ , δ )/||δ ||1 ≤ B(δ , δ )/||δ || to get In the following discussion, we apply mathematical induction to obtain the error result Take n = 1 in (46) and note that B 0 2 = 1; it is easy to find that the following inequality holds: Assuming that the inequalities ||δ || ≤ (C(u, σ, α, β)/B 2)(h + h + Δt 2 + Δt 1+2) hold for j = 1,2,…, n − 1, we now prove that the inequality (47) holds. Noting that 1/B < 1/B , j = 1,2, and combining (43) and (44) with (46), we have 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 n 1−2/B 2. We now use Taylor formula to arrive at where 0 < θ < 1 and 1/2 < α < 1. Noting that (nΔt)2 ≤ T 2 and (50), we have Substitute (51) into (33) and use (18) to get By combining (51) and (18) with triangle inequality, the L 2-norm estimate ||σ − σ || is got. Similarly, the estimate ||u − u || in the L 2-norm and the estimate ||u −u ||1 in the H 1-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 C(u, σ, α, β)/(2 − 2α) → ∞ as α → 1. So, we have to look for another error estimate's process. Noting the fact that nΔt ≤ T (n = 1,2,…, M), we can obtain the following error inequality:

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 H 1-GMFE Scheme for Several Spaces Variables

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

Let L 2(Ω) = (L 2(Ω)), (d = 2 or 3) with inner product and norm Further, let with norm Taking q = ∇u, we use a similar process to the system (10) and (11) to get the H 1-Galerkin mixed weak formulation {u, q}:[0, T] ↦ H 0 1 × W by The corresponding time semidiscrete system is defined by where R 0 = O(Δt), which can be estimated by a similar process to E 0 .

In order to get fully discrete mixed finite element scheme, we now choose the finite element spaces V ⊂ H 0 1 and W ⊂ W, which satisfy the following approximation properties: for 1 ≤ p ≤ ∞ and k, r positive integers [33], The fully discrete H 1-GMFE scheme is to find {u , q }:[0, T] ↦ V × W such that

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 q 0 = ∇u 0 and q (0) = ∇u 1, the following inequality holds:

By a similar discussion as in [32], we have Based on the given initial value conditions and (63), we arrive at

Theorem 7 .

The following stable inequality for the system (60) and (61) holds:

In (60), we take v = ∇u and use Cauchy-Schwarz inequality and Poincaré inequality to arrive at In (61), we choose w = q and use Cauchy-Schwarz inequality, Young inequality, and (66) to get By the simple simplification for (67), we easily get For the case n = 1 in (68), we multiply (36) to get easily Noting that B 0 = 1 = B 0 2 and Lemma 6, we use Cauchy-Schwarz inequality to get By the simple calculation for (70), we arrive at For the case n ≥ 1, using a similar process to the proof of Theorem 5 based on the mathematical induction, we can get 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 Further, let be the standard finite element interpolant of q.

Let , ; see [33, 44]; we obtain Now, we can get the following theorem of a priori error estimates based on the above contents.

Theorem 8 .

With q 0 = ∇u 0, q (0) = ∇u 1, u ∈ H 0 1∩H , and q ∈ H , there exists a positive constant C(u, q, α, β) free of space-time meshes h and Δt such that, for 1/2 < α < 1, and for α → 1

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 u(x, t) = t 2+sin(2πx), for all (x, t)∈[0,1]×[0,1], and the determined source term f(x, t) 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 H 1-GMFE method. Now, we divide interval [0,1] into equal-length intervals e = [x , x ], 0 ≤ i ≤ N − 1. Now we choose the piecewise linear spaces V = span⁡{φ 1, φ 2,…, φ } and W = span⁡{ψ 0, ψ 1,…, ψ } with index k = r = 1. Let Based on the above expressions (77), we can get the following equivalent algebraic equation of (25): where a = Δt −2/Γ(3 − 2α), b = Δt −/Γ(2 − α), , , ,

Based on the algebraic equation (78), we calculate and get the detailed numerical results listed in Tables 1–4.

Table 1

Spatial convergence results in L 2-norm with fixed Δt = 1/4000 and changed α.

L 2-norm	α	h_1=120	h_2=140	h_3=180	Order (h_1h_2) l l l l l	Order (h_2h_3) g g g g g
||u − u h||	0.6	5.163644E − 03	1.280561E − 03	3.129308E − 04	2.013469	2.039022
	0.7	5.051593E − 03	1.251162E − 03	3.044436E − 04	2.013469	2.039022
	0.8	4.916797E − 03	1.215362E − 03	2.936933E − 04	2.016333	2.049004
	0.9	4.758728E − 03	1.172117E − 03	2.795694E − 04	2.021459	2.067839
||σ − σ h||	0.6	4.396500E − 03	1.122373E − 03	3.219747E − 04	1.969804	1.801533
	0.7	5.141035E − 03	1.312040E − 03	3.759918E − 04	1.970247	1.803038
	0.8	6.036700E − 03	1.543010E − 03	4.444122E − 04	1.968013	1.795777
	0.9	7.087004E − 03	1.822010E − 03	5.343041E − 04	1.959645	1.769798

Table 4

Space-time convergence results in H 1-norm with h = 5Δt = 1/M and changed α.

H 1-norm	α	M 1 = 20	M 2 = 40	M 3 = 80	Order (M_1M_2)l^l^l^l^l	Order (M_2M_3)g_g_g_g_g
||u−u h||1	0.6	6.933582E − 01	3.484856E − 01	1.744222E − 01	0.992502	0.998515
	0.7	6.933669E − 01	3.484880E − 01	1.744230E − 01	0.992510	0.998519
	0.8	6.933819E − 01	3.484924E − 01	1.744244E − 01	0.992523	0.998525
	0.9	6.934089E − 01	3.485010E − 01	1.744275E − 01	0.992543	0.998536
||σ−σ h||1	0.6	4.375309E + 00	2.191391E + 00	1.096134E + 00	0.997538	0.999422
	0.7	4.375792E + 00	2.191479E + 00	1.096153E + 00	0.997639	0.999456
	0.8	4.376484E + 00	2.191614E + 00	1.096183E + 00	0.997779	0.999505
	0.9	4.377505E + 00	2.191830E + 00	1.096235E + 00	0.997973	0.999579

In Table 1, we calculate the numerical results of a priori error results and orders of convergence in L 2-norm for both u and σ with a fixed time step length Δt = 1/4000 and the changed spatial meshes h 1 = 1/20, h 2 = 1/40, and h 3 = 1/80. From the calculated data, we can clearly see that the orders of convergence in L 2-norm keep unchanged with the changed α = 0.6,0.7,0.8,0.9, which confirm the optimal second-order convergence results of H 1-Galerkin mixed finite element method. Similarly, we obtain optimal first-order rates of convergence of H 1-norm for u and σ in Table 2.

Table 2

Spatial convergence results in H 1-norm with fixed Δt = 1/4000 and changed α.

H 1-norm	α	h_1=120	h_2=140	h_3=180	Order (h_1h_2)l^l^l^l^l	Order (h_2h_3)g_g_g_g_g
	0.6	6.933461E − 01	3.484812E − 01	1.744206E − 01	0.992495	0.998510
||u−u h||1	0.7	6.933489E − 01	3.484816E − 01	1.744207E − 01	0.992499	0.998511
	0.8	6.933531E − 01	3.484823E − 01	1.744208E − 01	0.992505	0.998513
	0.9	6.933591E − 01	3.484833E − 01	1.744210E − 01	0.992513	0.998516
	0.6	4.374396E + 00	2.191164E + 00	1.096078E + 00	0.997386	0.999347
||σ−σ h||1	0.7	4.374658E + 00	2.191198E + 00	1.096083E + 00	0.997450	0.999363
	0.8	4.374979E + 00	2.191241E + 00	1.096089E + 00	0.997528	0.999383
	0.9	4.375365E + 00	2.191294E + 00	1.096098E + 00	0.997620	0.999406

In Table 3, we calculate the numerical results of a priori error results and orders of convergence in L 2-norm for both u and σ with different space-time step length h = 5Δt = 1/M  (M = 20,40,80). From the calculated data, we can find that the rates of convergence, which are higher than the results O(Δt 3−2 + h 2) of theory, gradually decrease with the increased α (which is taken from 0.6 to 0.9 with interval 0.1). In Table 4, some first-order convergence results for both u and σ in H 1-norm, which are unchanged with the changed α = 0.6,0.7,0.8,0.9, are given.

Table 3

Space-time convergence results in L 2-norm with h = 5Δt = 1/M and changed α.

L 2-norm	α	M 1 = 20	M 2 = 40	M 3 = 80	Order (M_1M_2)l^l^l^l^l	Order (M_2M_3)g_g_g_g_g
||u − u h||	0.6	4.781196E − 03	1.095478E − 03	2.262375E − 04	2.125811	2.275651
	0.7	4.587916E − 03	1.028885E − 03	2.011905E − 04	2.156757	2.354448
	0.8	4.319590E − 03	9.315674E − 04	1.630376E − 04	2.213162	2.514455
	0.9	3.941100E − 03	7.841911E − 04	2.154685E − 04	2.329321	1.863729
||σ − σ h||	0.6	6.937717E − 03	2.316455E − 03	8.737355E − 04	1.582542	1.406650
	0.7	8.221980E − 03	2.746089E − 03	1.033148E − 03	1.582108	1.410332
	0.8	1.000490E − 02	3.373945E − 03	1.275972E − 03	1.568199	1.402840
	0.9	1.251982E − 02	4.324761E − 03	1.669151E − 03	1.533521	1.373505

In view of the above analysis on the numerical results, we now announce that the time fractional telegraph equation can be well solved by the H 1-GMFE method.

6. Some Concluding Remarks and Extensions

As far as we know, more and more people have proposed and analyzed a lot of numerical methods for fractional partial differential equations. However, the discussions on mixed finite element methods for solving fractional partial differential equations are fairly limited. The a priori error analysis of H 1-GMFE method for time fractional telegraph equation, especially, has not been made and discussed. In this paper, we give the detailed proof's process of the error analysis on H 1-GMFE method for time fractional telegraph equation. Further, we provide a numerical procedure to verify the theoretical results of the studied method.

In the near future, our aim is to study an H 1-Galerkin moving mixed finite element method, which is based on a combination of H 1-GMFE methods and moving finite element methods.