High-order Compact Difference Schemes for the Modified Anomalous Subdiffusion Equation 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).

High-order Compact Difference Schemes for the Modified Anomalous Subdiffusion Equationthanks:  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).

Hengfei Ding,   Changpin Li
1. School of Mathematics and Statistics, Tianshui Normal University, Tianshui 741001, China
2. Department of Mathematics, Shanghai University, Shanghai 200444, China
Abstract

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

1 Introduction

The phenomenological diffusion equation can be derived by the following Fick’s first law [1] (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 [2].

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 [3]. The corresponding fractional Fick’s law has been proposed [4].

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 [5],

where , and are the anomalous diffusion coefficients. Thus, the modified fractional anomalous diffusion equation is obtained [6],

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 [18]. In [19], 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 [20], Mohebbi et al. considered an unconditionally stable difference scheme of order . Wang and Vong [21] 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 [22],

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:

and

then

and

hold.

Proof. In view of the following approximation scheme [23]

then one obtains

and

That is,

and

This completes the proof.   

Lemma 2 [24] 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 [25], Tuan and Gorenflo introduced the following fractional central difference operator

and proved that

where .

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,

Now set

then the numerical formula (8) becomes

where

Applying the Crank-Nicolson method to equation (3) yields

Setting

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):

and

It is obvious that the local truncation errors of difference schemes (13) and (14) are and , respectively.

3 Solvability Analysis

Denote

and

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 [26]. 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

and

respectively.

Note that and , . Thus

and

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)

Lemma 4 ([8, 27]). The coefficients satisfy

Lemma 5. The coefficients satisfy

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

and

respectively.

So, we can easily get the following roundoff error equation

Now, we define the grid functions

then can be expanded in a Fourier series

where

Let

By the Parseval equality