On a Novel Numerical Scheme for Riesz Fractional Partial Differential Equations

: In this paper, we consider numerical solutions for Riesz space fractional partial differential equations with a second order time derivative. We propose a Galerkin ﬁnite element scheme for both the temporal and spatial discretizations. For the proposed numerical scheme, we derive sharp stability estimates as well as optimal a priori error estimates. Extensive numerical experiments are conducted to verify the promising features of the newly proposed method.


Introduction
The Riesz fractional partial differential equations (RPDEs) arise in various applications of practical importance and their numerical solutions have received considerable attentions in the literature in recent years; see [1][2][3][4][5][6][7] and references cited therein.In this article, we aim to develop a novel numerical scheme to solve linear space RPDEs of second order in time.
There are several numerical methods developed for solving RPDEs.For the space RPDE problems of first order in time, we would like to mention the works [3,5,[7][8][9][10], where the authors developed various numerical schemes that combine finite-difference and finite-element methods.For the linear space RPDE problems of second order in time, we propose a space-time bilinear finite element scheme for the numerical approximation.In the temporal direction, unlike these finite difference schemes in [7,10,11], we develop a C 0 -continuous Galerkin linear finite element method to descritize the time derivatives (cf.[4,[12][13][14][15]). Since the approximate solution function is a continuous piecewise linear polynomial in the whole temporal interval, the regularity assumptions on the exact solution in the error analysis can be relaxed [4,16,17] (see Remark 1 in what follows).In the spatial direction, we apply the linear finite element method to discretize the space fractional order derivative.We establish sharp stability estimates for the proposed numerical method (see Lemma 1).By introducing an unusual interpolation operator Ĩ satisfied by (14) and (15) (cf.[13,18]), we further derive optimal a-priori error estimates respectively in the L 2 norm and the energy norm (see Theorem 1).
The rest of the paper is organized as follows.In Section 2, we introduce the problem setup and construct the numerical scheme for our subsequent study.In Section 3, we establish the stability estimate of the proposed numerical method and derive some related theoretical results.Applying the technique of dual argument [14,17], we derive the optimal error estimates in Section 4. The numerical examples demonstrating the promising features of the approximate method are presented in Section 5.The paper is concluded in Section 6 with some relevant discussions.

Problem Setup and the Numerical Scheme
Consider the linear space RPDE problem as follows: where K γ u(x, t) = −q γ ∂ 2γ u(x, t) ∂|x| 2γ , 1/2 < γ ≤ 1, q γ > 0, the Riesz fractional derivative [19,20] is given as Assume that the force f ∈ L 2 (I × I), and the initial displacement and velocity u 0 , In what follows we shall establish a novel scheme to compute the numerical solution of problem (1), for this scheme the linear finite elements are used both in spatial and temporal directions.For positive integers N, L, let 0 partitions of I and I, respectively.Denote the temporal subintervals Similarly, we introduce the spatial subintervals = 0 , where P 1 (e) = {w; w is a polynomial of degree ≤ 1 on set e}, and then we obtain the following bilinear finite element space, In the sequel, we let S n hk be the restriction of S hk to I n .To solve problem (1) numerically, the bilinear finite element method can now be formulated: determine U ∈ S hk satisfying with where the elliptic projection operator [5,17] The test function ẇ in (2) may be discontinuous at point t n .Therefore if we set ẇ to vanish outside I n , (2) becomes: with 0 ≤ n ≤ N − 1, we deduce easily Substituting Equation (5) into the first equation of (4) and taking ẇ to be ψ on I n with ψ ∈ S h 0 , leads the explicit representation of discrete scheme (2): find {U n } N n=0 ∈ S h 0 and with 0 ≤ n ≤ N − 1.Therefore, we solve the Equation (6) to implement the scheme (2) in practical calculation.Once U n and Un − are available, we deduce Un+1 − by solving the first equation of ( 6), then derive the function U on I × I n by applying the second equation of ( 6), and hence U n+1 .In this way, we finally obtain the finite element solution U in given domain.

Stability for the Numerical Method and Some Theoretical Results
From now on, we use the usual symbols for Sobolev spaces [21].For real r > 0, we also write • H r (I) and • L 2 (I) as • r and • 0 , respectively.Let χ : [0, T] −→ S be a Lebesgue measurable function with Banach space S endowed norm • S .Denote [22] where C represents a positive constant independent of the functions and parameters involved, which are not necessarily the same at different occurrences.
Applying the following coercivity and boundedness of B γ (•, •) [2,23], and using (7) we deduce that Using the above inequality and noting that U is a piecewise constant in the whole temporal interval I, we obtain U ∞, 0 The proof is complete.
Applying Lemma 1, we can easily know that the numerical method (2) for problem (1) is uniquely solvable.
In the following, we introduce the dual problem of (1) with f = 0, which satisfies Then we rewrite (9) in the same expression as for (1) by applying the change of temporal variables t = T − t, and solve this problem numerically by using the method (2).Hence, using change of temporal variables t = T − t second time, we establish the approximate method for the problem (9): find Θ ∈ S hk such that where ϕ 0 , ϕ 1 ∈ S h 0 are any two functions.By summing (10) with respect to n, we obtain that Moreover, using Lemma 1 we can estimate Θ ∞, 0 ϕ 0 γ + ϕ 1 0 (12) with the function Θ ∈ S hk given by (10).
For the elliptic projection operator R h given by (3), we need some fundamental results described as follows.

Convergence Analysis for the Numerical Method
In this section, in order to obtain error bounds for the numerical scheme (2), we suppose the exact solution u of problem (1) satisfies which implies functions u and K γ u are C 1 -continuous in the temporal direction [22].which implies that u and L α u are C 1 -continuous in the time direction [22].
Introduce the interpolation Ĩ u = ũ ∈ P 1 ( Ī) of u, where P 1 ( Ī) is the set of continuous functions that are first order polynomials on each interval I n , i.e., with 0 ≤ m ≤ N. Using the interpolation definition ( 14) and ( 15), for 0 and then Recalling Taylor's formula with the integral form of remainder [24], we have Define from Equations ( 16)-( 18) we find and therefore Hence, For we conclude from ( 16) that With the help of ( 19), we can estimate It is clear that combining the last three estimates with (20) we finally obtain For w ∈ S n hk , by multiplying (1) by ẇ and by integrating on I n , we obtain the equality [7,20] Denote e = U − u, since u is C 1 -continuous in the time direction, from (2) we conclude that We introduce the following essential identity, which is vital to error analysis for approximate method (2).Lemma 4. Let u solve (1) and Θ solve (10).Let I denote the identity operator.Under the regularity assumptions (13) we have Proof.We replace w with Θ in equation( 22) and derive Now we shall simplify equality (24).From the decomposition η = u − R h u + (I − Ĩ )R h u and (3) we find that [4,7] Since Θ = Θn+1 − = Θn + is constant on I n , from ( 14) and (15) we see that Therefore, with the aid of ( 25) and ( 26) we can rewrite the right side term of equality (24) as Applying integration by parts and observing the equality we have that Summing in n and noting that From (11) we have that Applying the last two equalities, we express the left side term of equality (24) in the form Hence, we obtain the identity (23) by using ( 24), (27), and (28).
Theorem 1.Let u and U be the solutions of ( 1) and ( 2), respectively.Under Assumption 1, and assume that u satisfies the regularity conditions that u ∈ L ∞ I; H q (I) Then, for n = 1, 2, . . ., N, we have the error estimates where Proof.We only give the proof where n = N; the other cases can be treated analogously.In (10) we set Lemma 4,(12), Lemma 3 and the interpolation error estimate (21), we deduce that Using Lemma 2, the triangle inequality and ( 19), we can estimate Summing up the last three estimates we find that Thus, the desired result (29) follows from ( 31) and ( 32) with the aid of the triangle inequality.Again, we choose (10), and then conclude the following estimate analogously ζN Now the statement of (30) follows from the above inequality and Remark 1.For the computational scheme in [8], the convergence analysis requires u ttt ∈ L ∞ (I; H 2 (I)) if the linear finite element is applied in the spatial direction.Moreover, for the numerical method in [4], the error estimate requires u ttt ∈ L 1 (I; H 2γ (I)).Hence, in the temporal direction the regularity assumes for the exact solution that u is relaxed in convergence analysis for the present numerical scheme (2).

Numerical Tests
In order to assess the numerical scheme (2), we display some numerical experiments.Solve the Equation (1), equipped with u 0 2 is the exact solution of this problem [4].
Here, we only take T = 1 in our computation.It should be emphasized that since the constants in our estimates do not depend on time, the method (2) is still valid for a larger increase in T.
For χ ∈ H γ 0 (I), since the norms χ γ and χ * γ := are equivalent [2,5], by Theorem 1 we see that the predicted convergence orders for , and the predicted convergence orders for For simplicity, we use the numerical scheme (6) to compute this problem with a uniform partition ( . Nevertheless, it is emphasised that the method is also effective for the variable step-size implementation.Tables 1-4 show the error and numerical convergence order in the spatial direction with k = 1.25 × 10 −6 for different γ values.Tables 5 and 6 show the error E 0 , E M 0 and numerical convergence order in the temporal direction with h = 1/300 for different γ values.Then we compute the error E γ , E M γ and numerical convergence order in the temporal direction.To decrease the computational cost (avoid computing with a very small step-size h), we choose h ≈ k 1 2−γ (h fixed, k = 1/Fix(h γ−2 ), where Fix is the Floor function, i.e., Fix(x) := max{z ∈ Z : z ≤ x} ).Tables 7-9 show these numerical results with different γ values.It is clearly seen that these computational results verify the expected results in Theorem 1.

Conclusions
In this paper we proposed a novel numerical scheme to solve the linear space RPDE problems of second order in time.We applied the linear finite element method to numerically discretize in both the spatial and temporal directions.The stability and the optimal a priori error estimates of the newly proposed scheme are proven.Extensive numerical experiments verified the theoretical results and demonstrated the promising features of the proposed methods.In future study, we shall apply the numerical method developed for solving several practically important inverse problems associated with fractional PDEs [25,26]; see also [27][28][29][30][31][32] for the related backgrounds of those inverse problems.

Table 6 .
E M 0 and numerical convergence order in the temporal direction (h = 1/300).

Table 7 .
E γ , E M γ and numerical convergence order in the temporal direction (γ = 0.65, h ≈ k

Table 8 .
E γ , E M γ and numerical convergence order in the temporal direction (γ = 0.75, h ≈ k

Table 9 .
E γ , E M γ and numerical convergence order in the temporal direction (γ = 0.95, h ≈ k