Abstract
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.
1. 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–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–9] are of the order , where . Here, we propose a method for the time fractional differential heat equations with the accuracy of order .
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 , 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]: If is nonzero, then there are some problems about the existence of the solutions for the heat equation (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.
2. Discretization of the Problem
In this section we introduce the basic ideas for the numerical solution of the time fractional heat equation (1) by compact finite difference scheme.
For some positive integers and , the grid sizes in space and time for the finite difference algorithm are defined by and , respectively. The grid points in the space interval are the numbers , , and the grid points in the time interval are labeled , . The values of the functions and at the grid points are denoted by and , respectively.
As in the classical Crank-Nicholson difference scheme, we use the approximation [13] to the fractional derivative at , and then where , , and , for .
Definition 2. Define the average operator as follows: where is a grid function and is the space of the grid functions.
Lemma 3. Suppose , then
Proof (see [14]). We use the Taylor expansion of each term about , and then we obtain the following truncation error for any , where :
3. 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 ; The difference scheme above can be written in matrix form as follows: where , , , , , and .
Here and are the matrices of the form:
We note that the unspecified entries are zeros at the matrices above.
Using the idea of the modified Gauss-Elimination method, we can convert (10) into the following form:
This way, the two-step form of difference schemes in (10) is transformed to one-step method as in (12).
Now, we need to determine the matrices and satisfying the last equality. Since , we can select and . Combining the equalities , and and the matrix equation (10), we have Then, we write where .
So, we obtain the following pair of formulas [15]: where .
4. The Spectral Stability of the Method
The spectral stability analysis is done by using the analysis of the eigenvalues of the iteration matrix () of the scheme (12).
Let denote the spectral radius of a matrix , that is, the maximum of the absolute value of the eigenvalues of the matrix .
We will prove that , , by induction
since is a zero matrix .
Moreover, , and is a lower triangular matrix of the following form:
Therefore < .
Since and , we can write .
Now, assume . After some calculations we find thatand we already know that , and , nonzero eigenvalue of for . If , then we have two subcases. (a)If , (b)If , If , then we have two subcases.(a)If , (b)If ,
So, we have proved that whenever then it follows that . So for any , where .
Remark 4 (see [16]). It is well known that for any , as if and only if . 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 guarantees eventual decay of the errors, but does not control the intermediate growth of the errors. Then, it is easy to understand that is a necessary condition for stability but not always sufficient.
5. The Fourier Stability of the Method
We analyze the stability of the difference scheme by a Fourier analysis. Let be the approximate solution and define , , . Then, we write and obtain the following roundoff error equation for (9):
We now define the grid functions: and then can be expanded in a Fourier series as follows: where and we introduce the following norm: and applying the Parseval equality , we obtain
Based on the above analysis we can suppose that the solution of (24) has the following form , where and . Substituting the above expression into (24), we obtain
After simplifications, we write
Theorem 5. for , where is the solution of (30).
Proof. We will use mathematical induction for the proof.
We start with , and then
Then
and therefore .
Now, assume that , . We need to prove that . Indeed,
Theorem 6. The finite difference scheme (9) is unconditionally stable.
Proof. Using the last theorem and Parseval equality, we obtain which means the proposed difference scheme is unconditionally stable.
6. Numerical Analysis
Example 7. One has
Exact solution of this problem is . 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.
(a)
(b)
The errors in Table 1 are calculated by the following formula:
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.
7. 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.