High-Order Compact Difference Scheme for the Numerical Solution of Time Fractional Heat Equations

A high-order finite difference scheme is proposed for solving time fractional heat equations. The time fractional derivative is described in the Riemann-Liouville sense. In the proposed scheme a new second-order discretization, which is based on Crank-Nicholson method, is applied for the time fractional part and fourth-order accuracy compact approximation is applied for the second-order space derivative. The spectral stability and the Fourier stability analysis of the difference scheme are shown. Finally a detailed numerical analysis, including tables, figures, and error comparison, is given to demonstrate the theoretical results and high accuracy of the proposed scheme.


Introduction
In the last decades, more and more attention has been placed on the development and research of fractional differential equations, because they can describe many phenomena, physical and chemical processes more accurately than classical integer order differential equations [1][2][3][4]. And the finite difference method is an efficient tool for solving fractional partial differential equations.
There are many different discretizations in time variable equipped with the compact difference scheme in spatial variable. The approximations given in [4][5][6][7][8][9] are of the order ( + ℎ 4 ), where 1 ≤ < 2. Here, we propose a method for the time fractional differential heat equations with the accuracy of order ( 2 + ℎ 4 ).
In this work, we consider the following time fractional heat equation: where the term ( , )/ denotes -order modifying Riemann-Liouville fractional derivative [10] given with the following formula: where Γ(⋅ ) is the gamma function.
Remark 1. If ( ) = 0, then the Riemann-Liouville and the modified Riemann-Liouville fractional derivatives are identical, since the Riemann-Liouville derivative is given by the following formula [11]:  (1). To rectify the situation two main approaches can be used; the modified Riemann-Liouville fractional derivative can be used [10] or the initial condition should be modified [12]. We chose the first approach in our work.
Proof (see [14]). We use the Taylor expansion of each term about , and then we obtain the following truncation error for any , where 1 ≤ ≤ − 1:

Compact Finite Difference Scheme
If we apply the operator to both side of (1), and use the approximation (4), then we obtain the following difference scheme which is accurate of order ( 2 + ℎ 4 ); The Scientific World Journal 3 The difference scheme above can be written in matrix form as follows: where ] . Here ( +1)×( +1) and ( +1)×( +1) are the matrices of the form: ] .
We note that the unspecified entries are zeros at the matrices above.

The Spectral Stability of the Method
The spectral stability analysis is done by using the analysis of the eigenvalues of the iteration matrix (1 ≤ ≤ ) of the scheme (12).

4
The Scientific World Journal Let ( ) denote the spectral radius of a matrix , that is, the maximum of the absolute value of the eigenvalues of the matrix .
Remark 4 (see [16]). It is well known that for any ∈ R × , → 0 as → ∞ if and only if ( ) ≤ 1. We note that if is normal, then ‖ ‖ = ( ) but when the matrix is not normal the spectral radius gives no indication of the magnitude of the roundoff error for finite . In this case a condition of the form ( ) ≤ 1 guarantees eventual decay of the errors, but does not control the intermediate growth of the errors. Then, it is easy to understand that ( ) ≤ 1 is a necessary condition for stability but not always sufficient.
Proof. We will use mathematical induction for the proof.

Numerical Analysis
Exact solution of this problem is ( , ) = ( 4 + 1) 5 (1 − ). The solution by the proposed scheme is given in Figure 1. The errors when solving this problem are listed in Table 1 for various values of time and space nodes.
The errors in Table 1 are calculated by the following formula: 8 The Scientific World Journal It can be concluded from Table 1 and Figure 1 above that when the time-step size is reduced by a factor of 1/4 and the spatial step size is reduced by a factor of 1/2, then the error decreases by about 1/16. The numerical results support the claim about the order of the convergence.

Conclusion
In this work, the compact difference scheme was successfully applied to solve the time fractional heat equations. The second order approximation for the Riemann-Liouville fractional derivative is equipped with the higher order compact difference schemes. The Fourier analysis and the spectral stability method are used to show that the proposed scheme is unconditionally stable. Numerical results are in good agreement with the theoretical results.