High-Order Algorithms for Riesz Derivative and Their Applications ( I )

and Applied Analysis 3 Therefore, the modified first-order approximation scheme is constructed as follows: ∂ α u (x m , t) ∂|x| α = − Ψ α h ( m+1


Introduction
In recent decades, fractional differential equations are undergoing rapid development due to their wide application in physics, engineering, economics, and some other research realms [1][2][3]. Several methods have been introduced for some special fractional differential equations to seek their analytical solutions, such as the integral transformation method (Laplace's transform, Fourier's transform, and Mellin's transform), Adomian decomposition method, the method of separating variables, and so on [4][5][6][7]. However, the exact solutions of most fractional differential equations cannot be obtained. So it becomes important to develop numerical methods for these equations. It is worth mentioning that the fractional linear multistep (also high order) methods for the Riemann-Liouville integrals and derivatives were firstly proposed in [8]. High-order numerical methods for Caputo derivatives were firstly constructed in [9]. Other fractional differential equations include Caputo or (and) Riemann-Liouville derivatives; there also existed some numerical approaches [10][11][12][13]. So, in this paper, our main aim is to construct high-order numerical methods for Riesz fractional derivative and Riesz fractional differential equations.
The Riesz fractional derivative was derived from the kinetics of chaotic dynamics [14,15]. For the Riesz fractional differential equations, there have existed several analytical and numerical methods. Zhang and Liu [16] studied the analytical solutions of space Riesz and time Caputo fractional partial differential equations. Later on, they [17] used Galerkin finite element approximation and a backward difference technique in Riesz fractional advection-dispersion equation. By using the Laplace and Fourier transform, Chen et al. [18] obtained the analytical solution for the space Riesz fractional reaction dispersion equation, and they also constructed an explicit finite difference scheme for it. Shen et al. [19] presented the explicit and implicit difference schemes for the Riesz fractional advection-dispersion equation and the space Riesz time Caputo fractional advection-dispersion equation. Recently, they [20] proposed a novel numerical method for the space Riesz fractional advection-dispersion equation based on fractional central difference operator. Ozdemir et al. [21] were concerned with the numerical solutions of a two-dimensional problem. Yang et al. [22] obtained some numerical solutions for the space Riesz fractional diffusion and advection-dispersion equations on a finite domain by using three numerical methods. Ç elik and Duman [23] ( , ) = 1 Γ (2 − ) ∫ ( − ) in which Γ(⋅) indicates Euler's Gamma function. The structure of the paper is outlined as follows. In Section 2, we list existing numerical approximation of the Riesz fractional derivative just for reference. Then, we focus on constructing one second-order scheme and two fourthorder numerical schemes. In Section 3, we only present one fourth-order numerical algorithm for the space Riesz fractional diffusion equation on a finite domain. In Sections 4 and 5, we discuss the stability and convergence of the difference scheme, respectively. In Section 6, some numerical examples are shown which support the theoretical results derived in the above sections. Finally, the last section concludes this paper.
(I) By the Standard Grünwald-Letnikov Formula. Based on the above assumption and (3), we can obtain the first-order approximation formula (II) By the Shifted Grünwald-Letnikov Formula. In [25], Meerschaert and Tadjeran show that the above standard Grünwald-Letnikov formula is often unstable for time dependent problems. Hence, they proposed the following shifted Grünwald-Letnikov formulas for the left and right Riemann-Liouville derivatives in order to overcome the instability: RL , Therefore, the modified first-order approximation scheme is constructed as follows: (III) By the L2 Approximation Method. Note that the left and right Riemann-Liouville derivatives can be rewritten as (1 < < 2): RL , Hence, we can obtain a first-order scheme for the left and right Riemann-Liouville fractional derivatives [22]: RL , where ( ) = ( + 1) 2− − 2− , = 0, 1, . . . , − 1, or = 0, 1, . . . , − − 1. Therefore, applying the above two formulas and (3) gives in which ( ) is defined as above. where in which in which with = − 1, , + 1.
Combining (14), (16), and (3) gives where ( ) (V) By the Fractional Centered Difference Method. In [27], Ortigueira introduced a symmetrical fractional centered difference operator as follows: Later on, Ç elik and Duman [23] proved that the above symmetrical fractional centered difference operator for the Riesz fractional derivative has the following relation: RL , where ℓ 1 and ℓ 2 are two arbitrary integers and in which ℓ 1 , ℓ 2 , and ℓ 3 are three arbitrary integers and ( Naturally, we can obtain the following second-order and third order numerical formulas for the Riesz fractional derivative: respectively.
Here, we construct another second-order scheme and two fourth-order numerical schemes for the Riesz fractional derivative. In order to construct the new computational schemes, we introduce the following theorem.
Lemma 1 (see [29]). Let > 0, ( , ) ∈ ∞ 0 (R) with respect to . The Fourier transforms of the left and right Riemann-Liouville fractional derivatives with respect to are wherê( , ) denotes the Fourier transform of the function ( , ) with respect to ; that is, In [30], Tuan and Gorenflo introduced the following left fractional central difference operator: Similarly, we define the following right fractional central difference operator: Analogous to the integer-order finite difference formula, we define the following fractional average operator: Then we can get the following fractional left and right average central difference operators based on (30), (31), and (32), respectively: 6 Abstract and Applied Analysis Here, we always assume that ±ℎ can commute with the infinite summation.
For the fractional left and right average central difference operators defined in (33) and (34), one has the following result.
uniformly holds for ∈ R.
Next, we construct fourth-order difference scheme for the left and right Riemann-Liouville derivatives based on (33) and (34) by the following theorem.
uniformly holds for ∈ R, where 2 denotes second-order central difference operator with respect to and is defined by Proof. It is similar to Theorem 2.
Combining (3) and Theorems 2 and 3, we can get the following high-order difference schemes for the Riesz derivative: Finally, we derive another fourth-order numerical method for the Riesz fractional derivative which is presented in the following theorem. 8

Abstract and Applied Analysis
Proof. Here, we use the Fourier transform method to prove it [23]. From [27], we know that the generating function with coefficients ( ) satisfies From (3) and Lemma 1, we get the Fourier transform of the Riesz fractional derivative as follows: Applying the Fourier transform to the difference operator with respect to and using (51) give Since ( , ) ∈ 7 (R) and its partial derivatives up to order seven with respect to belong to 1 (R), there exists a positive constant̃1 such that So, using (56) and (57) leads tô wherẽ3 =̃1̃2. At this moment, taking the inverse Fourier transformation in both sides of (55) and noting (52) give In view of (58), we have Abstract and Applied Analysis 9 wherẽ=̃3/(2 − ) ; that is to say, This finishes the proof.
Furthermore, (61) can be rewritten as

The Numerical Method for the Space Riesz Fractional Diffusion Equation
In this section, we only develop the fourth-order numerical method for the Riesz fractional diffusion equation (1). The second-order methods derived in the preceding section for (1) are easy so they are omitted here. By Taylor expansion, one has ( , −1 ) = ( , ) − ( , ) Let be the approximation solution of ( , ). Substituting (63) into (68) and removing the high-order terms, one has the following finite difference scheme for (1): where 1 = /48ℎ and 2 = ( /2ℎ )(1 + ( /12)).

Stability Analysis
Now we perform the detailed stability analysis for the difference scheme (71). Firstly, we introduce some lemmas and a definition for the following discussion.
Lemma 5 (see [31]). Let A be an − 1 order positive definite matrix. Then for any parameter ≥ 0, the following two inequalities: hold.
Definition 6 (see [32]). Let Toeplitz matrix Q have the following form: If the diagonals { } −1 =− +1 are the Fourier coefficients of function ( , ), that is, then the function ( , ) is called the generating function of Q .
where min and max are the minimum and maximum values of ( , ). Moreover, if min < max , then all eigenvalues of Q satisfy for all > 0. And furthermore if min ≥ 0, then Q is positive definite. respectively.

11
Therefore the generating function of the matrix is that is, It immediately follows that is a positive definite matrix in view of Lemma 7.
Proof. Let and̃be the exact and numerical solutions of difference equation (71). Since the matrix ( + ) −1 is invertible, then we can obtain the following error equation: where E =̃− and = ( + ) −1 ( − ). By (82), we have Using (83) and Lemma 5, one has which means difference scheme (71) is unconditionally stable, and so is the difference scheme (69).

Convergence Analysis
In this section, we give the local truncation error analysis of difference scheme (69). Proof. We define the local truncation error at the point ( , ) as Utilizing the Taylor expansion, we have Substituting (1) into (87), we get the following local truncation error: Next, we study the convergence of difference scheme (69).
Theorem 11. The order of the convergence of difference scheme (69) is O( 2 + ℎ 4 ); that is to say, there exists a positive constant such that Proof. Let = ( , ) − . Then it follows from (69) and (85) that then (90) can be written in a matrix form Furthermore, we can obtain where = ( /2)( + ) −1 . According to Lemma 5 and Theorem 10, we get Combining (93) and (94) gives Noting that ≤ = , we easily get where = 1 . This completes the proof.

Numerical Examples
In the present section, some numerical examples are presented to demonstrate the theoretical analysis.
The Riesz fractional derivative of the function ( ) is   Example 3. We consider the equation where ( , ) = (2 + 1) 2 4 (1 − ) 4 + 2 +1 sec ( 2 ) together with the initial condition ( , 0) = 0 and homogeneous boundary value conditions. The analytical solution of (99) is ( , ) = 2 +1 4 (1 − ) 4 . Table 3 lists the maximum error and time and space convergence orders for different . Numerical results show that the convergence order of the difference scheme (69) is ( 2 + ℎ 4 ), which are in line with our theoretical analysis. Figures 1 and 2 show the comparison of the analytical and numerical solutions with = 1.1 at = 0.2 and = 1.9 at = 0.8, respectively. It can be seen that the numerical solutions are in excellent agreement with the analytical solutions.     By Figures 3-6, we find that the corresponding two groups of figures are almost the same, respectively.

Conclusions
In this paper, we construct one second-order method and two fourth-order methods for the Riesz fractional derivative. Later on, a fourth-order numerical formula is used for the Riesz space fractional diffusion equation. The stability and convergence of this numerical method are analyzed by the matrix method in detail. Finally, some numerical results are given to demonstrate the effectiveness of numerical method. From the above results, it is possible to claim that the methods and techniques discussed in this paper are useful for solving some other fractional differential equations with Riesz derivatives.

Numerical solution
Analytical solution