High-order Compact Difference Schemes for the Modified Anomalous Subdiffusion Equation††thanks: The work was partially supported by the National Natural Science Foundation of China under Grant No. 11372170, Key Program of Shanghai Municipal Education Commission under Grant No. 12ZZ084, the grant of “The First-class Discipline of Universities in Shanghai”, the Scientific Research Program for Young Teachers of Tianshui Normal University under Grant No. TSA1405, and Tianshui Normal University Key Construction Subject Project (Big data processing in dynamic image).
In this paper, two kinds of high-order compact finite difference
schemes for second-order derivative are developed. Then a second-order numerical scheme for a Riemann-Liouvile
derivative is established based on a fractional centered difference operator. We apply
these methods to a fractional anomalous subdiffusion
equation to construct two kinds of novel numerical schemes. The solvability,
stability and convergence analysis of these difference schemes are
studied by using Fourier method. The convergence orders of
these numerical schemes are and , respectively.
Finally, numerical experiments are displayed which are in line with the theoretical analysis.
Key words: Modified anomalous subdiffusion equation; High-order compact difference schemes; Fourier method; Riemann-Liouville derivative; Grünwald-Letnikov derivative
The phenomenological diffusion equation can be derived by the following Fick’s first law  (which describes steady-state diffusion),
Combing the following conservation law of energy
one can obtain the diffusion equation below (also known as Fick’s second law or the heat equation)
This equation well characterizes the classic diffusion phenomenon .
However, if the diffusion is abnormal, that is to say, it follows non-Gaussian statistics or can be interpreted as the Lévy stable densities, then the above equation can not well describe such anomalous diffusion. Generally speaking, the fractional differential equations can well describe and model these anomalous diffusion phenomena . The corresponding fractional Fick’s law has been proposed .
Combination of this equation with equation (2) gives
where , and is the anomalous diffusion coefficient. If , it is just the normal diffusion equation. Here is the Riemann-Liouville operator, which is defined as follows:
where is the Gamma function.
Recently, a modified fractional Fick’s law has been used to describe processes that become less anomalous as time progresses by the inclusion of a secondary fractional time derivative acting on a diffusion operator ,
where , and are the anomalous diffusion coefficients. Thus, the modified fractional anomalous diffusion equation is obtained ,
Till now, various kinds of anomalous diffusion equations have been studied numerically, see [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17] and many the references cited therein. However, it seems that only a few numerical studies are available for the two-term subdiffusions of the above form.
In the present paper, we aim to study the following modified anomalous diffusion equation with a source term
subject to the initial and Dirichlet boundary value conditions
where , , and are suitably smooth.
Jiang and Chen proposed a collocation method based on reproducing kernels to solve a modified anomalous subdiffusion equation (3) with a linear source term on a finite domain . In , Liu et al., constructed a conditionally stable difference scheme for equation (3) with a nonlinear source term, and they proved that the convergence order is by the energy method. In , Mohebbi et al. considered an unconditionally stable difference scheme of order . Wang and Vong  presented a compact method for the numerical simulation of the modified anomalous subdiffusion equation (3), and they achieved the convergence order . The aim of this paper is to propose much higher order numerical methods for equation (3). We construct two kinds of high-order compact difference schemes and provide a detailed study of the stability and convergence of the proposed methods by using the Fourier method. We demonstrate that the convergence orders are and , respectively. One of advantages of compact difference schemes are that they can produce highly accurate numerical solutions but involves the less number of grid points. Thus, compact schemes result in matrices that have smaller band-width compared with non-compact schemes. For example, a sixth-order finite difference scheme involves seven grid points, while sixth-order compact difference scheme only needs five grid points. Another additional advantage of the compact high order methods is that the methods described here leads to diagonal linear systems, thus allowing the use of fast diagonal solvers. Different from the typical differential equations, even if we use the lower order methods for solving the fractional differential equations, we still need more calculations and strong spaces. If we use the higher order methods for fractional differential equations, the calculations and memory capacities can not remarkably increase. In this sense, the higher order numerical methods for fractional calculus and fractional differential equations attract more and more interest.
The rest of this article is organized as follows. In Section 2, we firstly develop a sixth-order and an eight-order difference scheme for second-order derivative, next a second-order numerical scheme for the Riemann-Liouville derivative is proposed. Applications of these methods to equation (3) give two effective finite difference schemes. The solvability, stability and convergence of the numerical methods are discussed in Sections 3, 4 and 5, respectively. The numerical experiments are performed for equation (3) with the methods developed in this paper are given in Section 6, which support the theoretical analysis. Finally, concluding remarks are drawn in the last section.
2 Numerical Schemes
Let and where the grid sizes in time and space are defined by and , respectively.
Define the following centered difference operator as
then we have
It is well known that a second-order approximation for the derivative is given by the following second-order centered difference scheme
A fourth-order compact difference scheme has also been constructed ,
Next, we develop two high-order compact difference schemes for the second-order spatial derivative by the following lemma.
Lemma 1. Define the following two operators:
Proof. In view of the following approximation scheme 
then one obtains
This completes the proof.
Lemma 2  For the suitably smooth function with respect to , arbitrary different numbers , and , one has
Next, we develop a second order numerical scheme for the Riemann-Liouvile derivative at nongrid points .
In , Tuan and Gorenflo introduced the following fractional central difference operator
and proved that
Accordingly, we obtain the following form at point in view of equation (6),
Letting , and gives the following second-order numerical formula by using equation (7) and Lemma 2,
then the numerical formula (8) becomes
Applying the Crank-Nicolson method to equation (3) yields
and substituting (9) into (10) leads to
Let be the approximation solution of . Noting equation (11) and substituting (4) and (5) into (12) give the following two finite difference schemes for equation (3):
It is obvious that the local truncation errors of difference schemes (13) and (14) are and , respectively.
3 Solvability Analysis
Then we obtain the matrix form of difference scheme (13)
where , , matrices are given in the Appendix I.
Similarly, the matrix form of the difference scheme (14) is given by
where matrices are also given in the Appendix I.
Remark 1. In difference schemes (13) and (14), there are some points , , and outside of the interval , denoted as ghost-points, that are generally approximated using extrapolation formulas, see Appendix II for more details.
Lemma 3 . A circulant matrix is a Toeplitz matrix in the form
where each row is a cyclic shift of the preceding row, then matrix has eigenvector
and the corresponding eigenvalue
Theorem 1. The difference equations (15) and (16) are both uniquely solvable.
Proof. From Lemma 3, we know that the eigenvalues of the matrices and are
Note that and , . Thus
Therefore, the above two matrices are both nonsingular. The difference equations (13) and (14) are uniquely solvable. The proof is complete.
4 Stability Analysis
In this section, we analyze the stability of the difference schemes (13) and (14) by using the Fourier method.
4.1 Stability Analysis of Numerical Scheme (13)
Proof. (i) From the above analysis, we easily obtain the expressions of , , and
One has for if .
(ii) In view of Lemma 4, it is not difficult to obtain these relations by direct computations.
Let be the approximate solution of (13) and define
So, we can easily get the following roundoff error equation
Now, we define the grid functions
then can be expanded in a Fourier series
By the Parseval equality