Numerical Methods for Solving Space Fractional Partial Differential Equations Using Hadamard Finite-Part Integral Approach

We introduce a novel numerical method for solving two-sided space fractional partial differential equations in two-dimensional case. The approximation of the space fractional Riemann–Liouville derivative is based on the approximation of the Hadamard finite-part integral which has the convergence order O(h3-α)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(h^{3- \alpha })$$\end{document}, where h is the space step size and α∈(1,2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha \in (1, 2)$$\end{document} is the order of Riemann–Liouville fractional derivative. Based on this scheme, we introduce a shifted finite difference method for solving space fractional partial differential equations. We obtained the error estimates with the convergence orders O(τ+h3-α+hβ)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(\tau +h^{3-\alpha }+ h^{\beta })$$\end{document}, where τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} is the time step size and β>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta >0$$\end{document} is a parameter which measures the smoothness of the fractional derivatives of the solution of the equation. Unlike the numerical methods for solving space fractional partial differential equations constructed using the standard shifted Grünwald–Letnikov formula or higher order Lubich’s methods which require the solution of the equation to satisfy the homogeneous Dirichlet boundary condition to get the first-order convergence, the numerical method for solving the space fractional partial differential equation constructed using the Hadamard finite-part integral approach does not require the solution of the equation to satisfy the Dirichlet homogeneous boundary condition. Numerical results show that the experimentally determined convergence order obtained using the Hadamard finite-part integral approach for solving the space fractional partial differential equation with non-homogeneous Dirichlet boundary conditions is indeed higher than the convergence order obtained using the numerical methods constructed with the standard shifted Grünwald–Letnikov formula or Lubich’s higher order approximation schemes.


Introduction
Consider the following space fractional partial differential equation, with 1 < < 2 , 0 ≤ x, y ≤ 1 , 0 < t < T: where f is a source/sink term and u 0 , 1 , 2 , 1 , 2 are well-defined initial and boundary values, respectively. Here the Riemann-Liouville left-sided fractional derivative R 0 D x f (x) is defined by Similarly, we may define the Riemann-Liouville right-sided fractional derivative as follows: There are several ways to approximate the Riemann-Liouville fractional derivative in the literature. Meerschaert and Tadjeran [30] used the Grünwald-Letnikov formula to obtain the first-order scheme O(h) to approximate the Riemann-Liouville fractional derivative. Lubich [27] introduced the higher order schemes with order O(h p ), p = 1, 2, ⋯ 6, to approximate the Riemann-Liouville fractional derivative. Diethelm [9,10] obtained the scheme to approximate the Riemann-Liouville fractional derivative with the convergence order O(h 2− ), 0 < < 2 using the Hadamard finite-part integral approach; see other higher order schemes to approximate the Riemann-Liouville fractional derivative in Li and Zeng [22].
Based on the different schemes for approximating the Riemann-Liouville fractional derivatives, many numerical methods are introduced for solving space fractional partial differential Eqs. (1)-(4): finite difference methods [4-6, 15-18, 21, 24, 25, 28, 29, 31-40], finite element methods [1-3, 7, 8, 11-14, 23, 26] and spectral methods [19,20]. Meerschaert and Tadjeran [30] introduced a shifted finite difference method based on the Grünwald-Letnikov formula for solving the two-sided space fractional partial differential equation in one-dimensional case and proved that the convergence order of the numerical method is O(h). Meerschaert and Tadjeran [29] also considered the finite difference method for solving the fractional advection-dispersion equation in one-dimensional case using the Grünwald-Letnikov formula. The second-order shifted finite difference methods for solving fractional partial differential equations based on the Grünwald-Letnikov formula are discussed in both oneand two-dimensional cases in Tadjeran et al. [35] and Tadjeran and Meerschaert [34]. Now (1) u t (t, x, y) = R 0 D x u(t, x, y) + R x D 1 u(t, x, y) + R 0 D y u(t, x, y) + R y D 1 u(t, x, y) + f (t, x, y), (2) u(t, 0, y) = 1 (t, y), u(t, 1, y) = 2 (t, y), (3) u(t, x, 0) = 1 (t, x), u(t, x, 1) = 2 (t, x), (4) u(0, x, y) = u 0 (x, y), we turn to the Lubich's higher order schemes. When Lubich's higher order schemes with no shifts are applied for solving space fractional partial differential equations, the obtained finite difference methods are unstable as for using the Grünwald-Letnikov formula. With shifted Lubich higher order methods, it shows that the corresponding numerical methods for solving space fractional partial differential equations have only first-order accuracy; see [4,5]. In [22,Section 2.2], Li and Zeng introduced other higher order schemes, for example, L2, L2C schemes, to approximate the Riemann-Liouville fractional derivative. However, to the best of our knowledge, there are no works available in the literature to use the L2 and L2C methods for solving space fractional partial differential equations. The numerical methods discussed in [22,Chapter 4] for solving space fractional partial differential equations are also based on the Grünwald-Letnikov formula and Lubich's higher order schemes; see other recent works for solving space fractional partial differential equations in [1-3, 5, 21, 24, 25, 39, 40]. All the numerical methods constructed using the Grünwald-Letnikov formula or Lubich's higher order methods for solving space fractional partial differential equations require the solution of the equation satisfies the homogeneous Dirichlet boundary condition. Otherwise, the experimentally determined convergence orders of such numerical methods are very low, e.g., see Table 4 in Example 2 in Sect. 3. Therefore, it is interesting to design some numerical methods which have the higher order convergence for solving the space fractional partial differential equation with respect to both homogeneous and non-homogeneous Dirichlet boundary conditions. The purpose of this paper is to introduce such finite difference methods for solving the space fractional partial differential equation.
Recently, Ford et al. [17] considered the finite difference method for solving the space fractional partial differential equation in one-dimensional case where the Riemann-Liouville fractional derivative is approximated using the Hadamard finite-part integral; see also [15,16,37]. The convergence order O( + h 3− + h ), ∈ (1, 2), > 0 of the numerical method in [17] is proved in the maximum norm for both homogeneous and non-homogeneous Dirichlet boundary conditions. In this paper, we will extend the method in Ford et al. [17] to solve space fractional partial differential equations in twodimensional case. The corresponding error estimates in this paper are proved using a completely different way from Ford et al. [17]. The error estimates with the convergence order O( + h 3− + h ), ∈ (1, 2), > 0 hold for both homogeneous and non-homogeneous Dirichlet boundary conditions.
The main contributions of this paper are as follows.
(i) A new finite difference method for solving space fractional partial differential equations in two-dimensional case is introduced and the convergence order is

3 2 The Finite Difference Method
In this section, we shall extend the method in Ford et al. [17] for solving the space fractional partial differential equation in one-dimensional case to solve space fractional partial differential equations (1)-(4) in two-dimensional case. For simplicity with the notations, we also assume that the boundary values are equal to 0, i.e., 1 = 2 = 1 = 2 = 0.

and
The remainder term R l (g) satisfies, for every g ∈ C 3 (0, 1), Similarly, we may consider the approximation of the right-sided Riemann-Liouville frac- Using the same argument as for the approximation of R To obtain a stable finite difference method, we will consider the following shifted equation: where the errors produced by the shifted terms are denoted by We assume that R 0 D x u(t, x, y), R 0 D y u(t, x, y) satisfy the following Hölder conditions.

Assumption 1
For any x * , x * * , y * , y * * ∈ ℝ , there exist constants C > 0 and > 0 such that We also assume that R x D 1 u(t, x, y), R y D 1 u(t, x, y) satisfy the following Hölder conditions.

Assumption 2
For any x * , x * * , y * , y * * ∈ ℝ , there exist constants C > 0 and > 0 such that Remark 1 To make Assumptions 1 and 2 hold, we need to assume that the solution u satisfies some regularity conditions. In some circumstances, such conditions are easy to check, for example, We now turn to the discretization scheme of (10). Discretizing u t at t = t n+1 using the backward Euler method and discretizing R , y = y m−1 , respectively, using the Diethem's finite difference method introduced in Lemma 1, we obtain where S n+1 l,m can be defined as S n+1 2j in [17, (27)] and the weights (12) are defined by (7) and (8).

3
Then we have Let U n l,m ≈ u(t n , x l , y m ) denote the approximate solution of u(t n , x l , y m ) . We define the following finite difference method for solving (1)- (4): where Q n+1 l,m is some approximation of S n+1 l,m , defined as in [17, (30)] which satisfies Now we come to our main theorem in this work.
Theorem 1 Assume that u(t n , x l , y m ) and U n l,m are the solutions of (12) and (16), respectively. Assume that Assumptions 1 and 2 hold. Then there exists a norm ‖ ⋅ ‖ such that Proof Let e n l,m = u(t n , x l , y m ) − U n l,m . Subtracting (16) from (12), we get the following error equation, with = ∕h : (e n+1 l,m − e n l,m ) −      (23) Ae n+1 = e n + R n+1 , Here where and, with E = I (M−1)×(M−1) , We shall show that there exists a norm ‖ ⋅ ‖ such that Assume (24) holds at the moment, we have, by (23), noting that n = t n ≤ T, where we use the fact e 0 = 0. It remains to show (24). It suffices to show all the eigenvalues of A are greater than or equal to 1, which implies that all the eigenvalues of A −1 are less than or equal to 1. If all the eigenvalues of A −1 are less than or equal to 1, then there exists some norm ‖ ⋅ ‖ such that ‖A −1 ‖ ≤ 1 [33]. To show all the eigenvalues of A are greater than or equal to 1, we may use the well-known Gershgorin lemma. Let We have which imply that |a l,k |.
By Lemma 2, we get which implies that all the eigenvalues of A satisfy, by Gershgorin lemma, that is, all the eigenvalues of A are greater than 1 which implies (24). Together these estimates complete the proof of Theorem 1.

Numerical Examples
We shall consider in this section four numerical examples to illustrate that the numerical results are consistent with our theoretical results.
where It is easy to check that u(t, x, y) = 4e −t x 2 (2 − x) 2 y 2 (2 − y) 2 is the exact solution.
Note that the error estimate satisfies, by Theorem 1, with = min(3 − , ), a l,l − r l > 1, l = 1, 2, ⋯ , (M − 1) 2 , 1 < a l,l − r l < < a l,l + r l , In the numerical method (16), we simply ignore the errors n+1 l,m in (11) which are produced by the shifted terms. Of course, if we use the numerical methods (16) to calculate the approximate solutions, the spatial error should be O(h + h 3− ) . Since the exact solutions are given in our numerical examples, the errors n+1 l,m in (11) produced by the shifted terms can be calculated exactly. Thus, the convergence order should be O(h 3− ) if we include n+1 l,m in the numerical method (16). In general, we do not know the exact solutions of the equation. In such case, we may approximate n+1 l,m using the computed solutions U n to improve the convergence orders. In all our numerical simulations in this section, the numerical method (16) will include n+1 l,m defined by (11), which makes the experimentally determined order of convergence (EOC) independent of > 0.
We will observe the convergence orders with respect to the space step size. To see this, we shall choose sufficiently small time step size = 2 −10 and the different space step sizes h l = 2 −l , l = 2, 3, 4, 5, 6 such that the computational error is dominated by the space step size O(h 3− ), 1 < < 2 . Denote ‖e N l ‖ = ‖U N − u(t N )‖ the L 2 norm of the error at t N = 1 calculated with the step size h l . We then have which implies that Hence, the convergence order satisfies The experimentally determined orders of convergence (EOC) for the numerical method (16) are provided in Table 1 with respect to the different . We observe that the convergence order is indeed O(h 3− ) which is consistent with Theorem 1.
Next, we solve the equation in Example 1 using the finite difference method introduced in Meerschaert and Tadjeran [30] where the Riemann-Liouville fractional derivatives are approximated using the Grünwald-Letnikov formula which requires the solution of the equation satisfies the homogeneous Dirichlet boundary condition; see some other shifted and weighted Grünwald difference operator to approximate the Riemann-Liouville fractional derivative in [36]. The convergence order of the finite difference method in [30] is O(h) and we indeed observe this in Table 2 for solving (33)- (35). From now on, we call the finite difference method in Meerschaert and Tadjeran [30] as "the shifted Grünwald-Letnikov method ".

Example 2
In this example, we will consider the following space fractional partial differential equation with non-homogeneous Dirichlet boundary conditions, with 1 < < 2 , 0 ≤ x, y ≤ 2: where It is easy to see that u(t, x, y) = 4e −t x 2 (2 − x) 2 y 2 (2 − y) 2 + 5 is the exact solution of the equation.
In Table 3, we show the convergence orders using the numerical method (16). We see that for some , the convergence orders can reach O(h 3− ) and for some other the convergence orders are less than O(h 3− ) . But in most cases, the convergence orders of the numerical method (16) are greater than 1 for solving (30)- (32) with the non-homogeneous Dirichlet boundary conditions.
In Table 4, we use "the shifted Grünwald-Letnikov method" introduced in Meerschaert and Tadjeran [30] for solving (30)- (32). We observe that the convergence orders are very low because of the non-homogeneous boundary conditions. From Tables 3 and 4, we observe that the numerical method (16) introduced in this paper has higher order convergence than "the shifted Grünwald-Letnikov method" introduced in Meerschaert and Tadjeran [30] for solving space fractional partial differential equations with non-homogeneous boundary conditions. In the next example, we shall investigate the convergence orders of the numerical method (16) for solving space fractional partial differential equations where the solutions of the equations are not sufficiently smooth.
Here the exact solution has the form u(t, x, y) = e −t x 1 y 1 . We will consider two different 1 : the nonsmooth solution case with 1 = and the smooth solution case with 1 = 3.
In Table 5, we obtain the experimentally determined orders of convergence (EOC) for the different = 1.2, 1.4, 1.6, 1.8 . Since the solution is not sufficiently smooth, the convergence orders are less than O(h 3− ) as we expected.
For the smooth solution case with 1 = 3 , in Table 6, we observe that the convergence orders are almost 3 − as we expected.
In our final example, we consider a two-sided space fractional partial differential equation.
In Table 7, we observe that the convergence orders of the numerical method (16) for solving this equation are also O(h 3− ) as we expected.

Conclusions
In this paper, we construct a new and reliable finite difference method for solving the space fractional partial differential equations. The error estimates are proved and the convergence order of the numerical method depends on the smoothness of the solution of the equation. The convergence orders are proved for both homogeneous and non-homogeneous Dirichlet boundary conditions. Numerical examples show that the proposed numerical method in this paper has much higher convergence order than the shifted Grünwald-Letnikov method proposed in Meerschaert and Tadjeran [30] for solving space fractional partial differential equations with non-homogeneous boundary conditions. (38) u(0, x) = 4x 2 (2 − x) 2 y 2 (2 − y) 2 ,