Compact difference schemes for the modified anomalous fractional sub-diffusion equation and the fractional diffusion-wave equation
In this paper, compact finite difference schemes for the modified anomalous fractional sub-diffusion equation and fractional diffusion-wave equation are studied. Schemes proposed previously can at most achieve temporal accuracy of order which depends on the order of fractional derivatives in the equations and is usually less than two. Based on the idea of weighted and shifted Grünwald difference operator, we establish schemes with temporal and spatial accuracy order equal to two and four respectively.
Keywords: Modified anomalous fractional sub-diffusion equation, Fractional diffusion-wave equation, Compact difference scheme, Weighted and shifted Grünwald difference operator
Since fractional differential equations turn out to model many physical processes more accurately than the classical ones, in the past decades, increasing attentions have been made on these equations. Readers can refer to the books [1, 2] for theoretical results on fractional differential equations. This paper concerns with methods for obtaining accurate numerical approximations to the solutions of fractional sub-diffusion equations and fractional diffusion-wave equations. A fractional sub-diffusion equation is an integro-partial differential equation obtained from the classical diffusion equation by replacing the first-order time derivative by a fractional derivative of order between zero and one. When the time derivative is of order between one and two, we get a fractional diffusion-wave equation. Fractional derivatives of order between zero and one are widely used in describing anomalous diffusion processes , while fractional diffusion-wave equations have applications in modeling universal electromagnetic, acoustic, and mechanical responses [4, 5].
Numerical methods for the modified anomalous fractional sub-diffusion equation and the diffusion-wave equation have been considered by many authors, one may refer to – and the references therein. We point out here that one of the main tasks for developing accurate finite difference scheme of fractional differential equations is to discretize the fractional derivatives. Noticing that fractional derivatives are defined through integrals, one can approximate the derivatives by interpolating polynomials . Using this idea, Sousa and Li developed a second order discretization for the Riemann-Liouville fractional derivative .
We note that there are alternative ways to tackle the problem. In , a shifted Grünwald formula was proposed by Meerschaert and Tadjeran to approximate fractional derivatives of order for fractional advection-dispersion flow equations. It is also worth to mention that stability of forward-Euler scheme and weighted averaged difference scheme based on Grünwald-Letnikov approximation were analyzed in [21, 22] for time fractional diffusion equations. Very recently, accurate finite difference schemes based on weighted and shifted Grünwald difference operator were developed for solving space fractional diffusion equations in [23, 24].
Inspired by their work on the weighted and shifted Grünwald difference operator, in this paper, we establish high order schemes for the fractional diffusion-wave equation and the modified anomalous fractional sub-diffusion equation, which were proposed and studied in [8, 9] and –, respectively. We remark that the order of temporal accuracy of schemes proposed previously can at most be a fraction depending on the fractional derivatives in the equations and is usually less than two. We show that the schemes proposed in this paper are of order , where and are the temporal and spatial step sizes respectively.
This paper is organized as follows. Some preliminaries will be given in the next section. Compact schemes are proposed and studied in Section 3 and 4 for the modified anomalous fractional sub-diffusion equation and the fractional diffusion-wave equation respectively. In the last section, numerical tests are carried out to justify the theoretical results.
We first recall that the Caputo fractional derivative of order for a function is defined as
with being the gamma function and, for , the Riemann-Liouville fractional derivative of order for is defined as
Closely related to the fractional derivatives of a function is the Riemann-Liouville fractional integral which is given by
In order to develop second order approximation of the Riemann-Liouville fractional derivative, we consider the shifted Grünwald approximation  to the Riemann-Liouville fractional derivative given by
where for . Inspired by , we similarly introduce the shifted operator to the Riemann-Liouville fractional integral defined by
where for .
The following lemma was given in .
Suppose . The Fourier transform of the Riemann-Liouville fractional integral satisfy the following:
where denotes the Fourier transform of .
Lemma 2.1 can be used to obtain
(i) Let and its Fourier transform belong to , and define the weighted and shifted Grünwald difference operator by
then we have
for , where and are integers and .
(ii) Let and belong to . Define the weighted and shifted difference operator by
then we have
for , where and are integers and .
Proof. The proof of (i) can be found in . The proof of (ii) is similar to that of (i) but we give it here for the completeness of our presentation.
Referring to the definition of , we let
Taking Fourior transform on (2), we get
In order to achieve second order accuracy, we let the coefficients and satisfy the following system:
which implies that and .
We are now ready to establish our high order compact schemes in the next two sections.
3 A compact scheme for the fractional sub-diffusion equation
Consider the following modified anomalous fractional sub-diffusion equation
where , . We have used and to denote the Riemann-Liouville fractional operators and with respect to the time variable .
(i) We note that, in [25, 26, 27], the fractional derivatives in the equation are of order and with . We have changed the notations here in order to match the presentation for the two types of equations discussed in this paper.
(ii) Without loss of generality, we have assumed the initial condition . If , one may consider the equation for instead.
To develop a finite difference scheme for the problem (5)–(6), we let and be the spatial and temporal step sizes respectively, where and are some given positive integers. For and , denote . For any grid function , we introduce the following notations:
With this notations, we study the problem under the following inner product and norms:
It is easy to check that, if , the following identity holds
which plays important role in our analysis.
We can therefore consider an weighted Crank-Nicolson type discretization for equation (5) given by
In order to raise the accuracy in the spatial direction, we need the following lemma:
() Denote . If , then it holds that
Based on Lemma 3.1, we therefore propose the following compact scheme:
It is easy to see that at each time level, the difference scheme (3)–(9) is a linear tridiagonal system with strictly diagonal dominant coefficient matrix, thus the difference scheme has a unique solution.
The following lemmas are critical for establishing the convergence of the proposed scheme.
Let be defined as in (7), then for any positive integer and real vector , it holds that
Proof. For simplicity of presentation, in this proof, we denote without ambiguity. One can easily check that, to prove the above quadratic form is nonnegative is equivalent to proving the symmetric Toeplitz matrix T is positive semi-definite, where
Notice that the generating function (see ) of T is given by
As mentioned in , we only need to consider the principal value of on which gives
Let , then one can easily check that
Therefore is nondecreasing with respect to and , which implies that . The lemma now follows as a result of the Grenander-Szegö Theorem .
We note here that the function can not be proved to be nonnegative by considering differentiation with respect to as in .
(Grownall’s inequality ) Assume that and are nonnegative sequences, and the sequence satisfies
where . Then the sequence satisfies
Proof. We can easily get the error equation:
Denote , multiplying (10) by and summing in , we obtain
Summing up for and noticing that
we get, by Lemma 3.2, that
then the desired result follows by Lemma 3.3.
We can therefore obtain