A note on the stability of a second order finite difference scheme for space fractional diffusion equations

The unconditional stability of a second order finite 
difference scheme for space fractional diffusion equations is proved 
theoretically for a class of variable diffusion coefficients. In 
particular, the scheme is unconditionally stable for all one-sided 
problems and problems with Riesz fractional derivative. For problems 
with general smooth diffusion coefficients, numerical experiments 
show that the scheme is still stable if the space step is small 
enough.

1. Introduction. In this paper, we consider the following two-sided one dimensional fractional diffusion equation (FDE) ∂u(x, t) ∂t where the diffusion coefficients d ± (x) ≥ 0. Here a D α x u(x, t) and x D α b u(x, t) denote the α-order left and right Riemann-Liouville (RL) fractional derivatives of u, respectively, and they are defined as scheme based on the GL formula for time dependent FDE is unstable [8]. To overcome this problem, Meerschaert and Tadjeran propose to use the shifted GL approximation. For solving FDE (1), backward Euler scheme with shifted GL difference (BE-SGD scheme) is proposed in [9]. The BE-SGD scheme is proved to be unconditionally stable by using the Gershgorin theorem since the differentiation matrix for shifted GL discretization is diagonally dominant. Some efficient solvers for the BE-SGD scheme are provided in [7,17,18,19,20,10]. However, the BE-SGD scheme is only first order accurate.
To achieve second order accuracy, Crank-Nicolson scheme with shifted GL discretization is proposed in [14,15]. With Richardson's extrapolation, a numerical solution with second order accuracy in both time and space can be obtained. However, the computational cost of this approach is high since an extra numerical solution on fine grid is required. Recently, second order discretizations for RL fractional derivatives are being studied. Sousa and Li [12], and Tian et al. [16] present two different second order discretizations which possess an advantage that their differentiation matrices preserve the Toeplitz structure as the one for the shifted GL discretization. Nevertheless, they are not diagonally dominant for all α ∈ (1, 2). Due to this reason, the unconditional stability of those second order finite difference schemes proposed in [12,16,5] cannot be obtained analogously. In fact, those schemes are only proved to be unconditionally stable for the case with constant diffusion coefficients.
In this paper, stability of the Crank-Nicolson scheme with Sousa's second order discretization is studied. The scheme is proved to be unconditionally stable if (i) the diffusion coefficients d ± (x) = k ± d(x) where k ± ≥ 0 and d(x) ≥ 0; or (ii) the order of the fractional derivative α ∈ (α 0 , 2) where α 0 ≈ 1.5545. In particular, the scheme is unconditionally stable for all one-sided problems, and problems with Riesz derivative [1] instead of the RL derivative. For FDE (1) with smooth diffusion coefficients, numerical examples reveal that the scheme is still stable when the number of mesh points in space is large enough even though no restriction is imposed on the time step.
3. Stability and convergence of the finite difference scheme. Suppose u m and v m are two solutions of the finite difference scheme (5), and denote and furthermore, we have η m = H m η 0 . If the spectral radius of H is not larger than 1, then there exists a constant C such that ||H m || ≤ ||H|| m ≤ C < ∞ for some norm, and hence we have which means that the scheme is stable in some sense.
Recall that a matrix is said to be negative semi stable if each of its eigenvalue has nonpositive real part. Now, by showing that the matrix B is negative semi stable, two theorems are presented for the stability of the finite difference scheme (4). To prove the two theorems, we need the following lemmas.  [12,13]) For any α ∈ (1, 2), the sequence {q (3) satisfies the following properties: where α 0 ≈ 1.5545 is the root of the equation (6) is negative semi stable, and the symmetric matrix Q + Q T is negative semi definite. Lemma 3.3), we get that Σ + Σ T is negative semi definite and hence, by Lemma 3.2, the matrix Σ is negative semi stable. Now, since D is a diagonal matrix with nonnegative entries, by the Corollary 1.7.7 in [6], the matrix B = DΣ is also negative semi stable. Remark 1. We want to emphasize that our proof here works for the case that D has zero entries on its main diagonal. If all diagonal entries of D are strictly positive, D − 1 2 exists. Then we can prove that B = DΣ is negative semi stable by the following way. Note that B is similar to , and the matrix C is negative semi stable since is negative semi definite. Therefore B is negative semi stable. A similar proof is also given in the Theorem 4.1 in [2]. However, if D has zero entries on its main diagonal, the above proof cannot work since D − 1 2 does not exists.
Remark 2. The result in Theorem 3.4 reveals that the finite difference scheme when applied to one-sided FDE problems (k + = 0 or k − = 0), and FDE problems with Riesz fractional derivative (k + = k − ) is unconditionally stable for any diffusion coefficients.
By noting that Q is diagonally dominant if α ∈ (α 0 , 2), we have another theorem for the unconditional stability of the finite difference scheme.
Proof. For α ∈ (α 0 , 2), from Lemma 3.1, we have q 3 , · · · ≥ 0, and according to the Gerschgorin theorem, all eigenvalues of B are in the disks centered at < 0, i = 1, . . . , N , with radius and . . , N − 1. Therefore, every disk is contained in the left-half complex plan, and hence B is negative semi stable.
Remark 3. For the case α ∈ (1, α 0 ] with general diffusion coefficients d ± (x), the matrix B may not be negative semi stable. For example, consider N = 3 and α = 1.01, with Now, one can prove the second order convergence of the finite difference scheme (4) by using the similar technique given in [2,3].    Table 1. Maximum errors and observed orders for Example 1 with τ = h = 1/(N + 1), t = 1, and different α.
The exact solution of this example is u(x) = e −t x 2 (1 − x) 2 . From Figure 1, it is clear that ϑ ≤ 0 for all N , and hence the matrix B is negative semi stable. This supports the theoretical analysis in Theorem 3.4. Also, Table 1 shows the second order convergence of the finite difference scheme for α = 1.01, 1.1, 1.5, and 1.9. Now, we present other two examples in which the diffusion coefficients are not in the form of (7).
The exact solution of Example 2 is u( The exact solution of this example is u( In Figures 2 and 3, we plot the curves of ϑ defined in (8) versus N for α = 1.01, 1.1, 1.5 and 1.9. One can see that ϑ ≤ 0 for all N when α = 1.9 ∈ (α 0 , 2). This supports the theoretical result of Theorem 3.5. For other values of α being tested, we cannot guarantee that the matrix B is negative semi stable theoretically. Numerically, as shown in the left of Figures 2 and 3, when α = 1.01 and 1.1, it really occurs that ϑ > 0 for some N . However, in the right of Figures 2 and 3, it is clear that when N is large enough (N ≥ 100), we have ϑ ≤ 0, i.e., B is negative semi stable. Therefore, we guess that, for FDEs with smooth diffusion coefficients, the finite difference scheme (4) is stable provided that the number of node points in space is large enough even though no restriction is imposed on the time step τ . Tables 2 and 3 also show the second order convergence of the scheme.  Table 3. Maximum errors and observed orders for Example 3 with τ = h = 1/(N + 1), t = 1, and different α.

Conclusion.
A second order finite difference scheme for one-dimensional FDE (1) is studied in this paper. The scheme is proved to be unconditionally stable and convergent if (i) d ± (x) = k ± d(x) where k ± ≥ 0 and d(x) ≥ 0; or (ii) α ∈ (α 0 , 2) where α 0 ≈ 1.5545. Numerical examples are shown to support the theoretical results. Moreover, for FDEs with general smooth diffusion coefficients, numerical experiments reveal that the scheme is still stable if the number of node points in space is large enough (≥ 100 in our numerical tests), even though no restriction is imposed on the time step. Besides, we mention that our stability analysis for scheme (4) also works if the second order discretization for fractional derivatives are replaced by the one proposed in [16]. Similar theoretical and numerical results on the stability and convergence can be obtained analogously. Furthermore, we remark that the stability of finite difference method for a fractional-order differential linear complementarity problem arising in pricing option was studied in [4] recently. However, the approach in our paper is different from that used in [4].