Finite Difference Method for Time-Space Fractional Advection–Diffusion Equations with Riesz Derivative

In this article, a numerical scheme is formulated and analysed to solve the time-space fractional advection–diffusion equation, where the Riesz derivative and the Caputo derivative are considered in spatial and temporal directions, respectively. The Riesz space derivative is approximated by the second-order fractional weighted and shifted Grünwald–Letnikov formula. Based on the equivalence between the fractional differential equation and the integral equation, we have transformed the fractional differential equation into an equivalent integral equation. Then, the integral is approximated by the trapezoidal formula. Further, the stability and convergence analysis are discussed rigorously. The resulting scheme is formally proved with the second order accuracy both in space and time. Numerical experiments are also presented to verify the theoretical analysis.


Introduction
The concepts of fractional calculus and entropy are becoming more popular for analyzing the dynamics of complex systems. The idea of entropy was introduced in the field of thermodynamics by Clausius (1862) and Boltzmann (1896) and was later applied by Shannon (1948) and Jaynes (1957) in information theory. Recently, more general entropy measures have being proposed for application in several types of complex systems due to the relaxation of the additivity axiom. The concept of entropy for analyzing the dynamics of multi-particle systems with integer and fractional order behavior was proposed in [1]. The entropy production rate for the fractional diffusion process was calculated in [2]. In [3] it has been shown that the total spectral entropy can be used as a measure of the information content in a fractional order model of anomalous diffusion.
Fractional calculus has been applied to almost every field of science, engineering, and mathematics during the last decades [4][5][6][7][8]. Particularly fractional calculus has significant impact in the fields of viscoelasticity and rheology, physics, electrical engineering, electrochemistry, signal and image processing, biology, biophysics and bioengineering, mechanics, mechatronics, and control theory.
Fractional calculus is indeed a worthwhile mathematical tool that can undertake more than integer calculus. The monographs authored by Samko, Kilbas, Marichev [9], Podlubny [10] and Kilbas, Srivastava, Trujillo [7] have been helpful in understanding the theory and applications of fractional differential equations.
Numerous numerical methods have been proposed for solving the time-fractional differential equations. In this paper, we convert the fractional differential equation into the equivalent integral equation. Then, fractional trapezoidal formula is used to approximate fractional integral which has second-order accuracy [11,12]. Early in 1993, Tang [13] presented a finite difference method for the numerical solution of the partial integro-differential equations with a weakly singular kernel based on the product trapezoidal formula. Chen et al. [14] proposed fractional trapezoidal rule (FTR) type difference scheme by combining the second order difference quotient for spatial discretization and the FTR alternating direction implicit method in the time stepping for a two-dimensional fractional evolution equation. Chen et al. [15] derived a fractional trapezoidal rule type difference scheme for fractional order integro-differential equation with second order accuracy both in temporal and spatial directions. Recently, a finite difference scheme has been proposed in [16] to solve time-space fractional diffusion equation of second-order accuracy in both time and space by employing trapezoidal rule. Numerical schemes for linear and nonlinear time-space fractional diffusion equations were constructed in [17] using the trapezoidal formula for temporal approximation and the centred difference approximation for the spatial Riesz fractional derivative. Several numerical schemes have been proposed to approximate Riesz fractional derivative based on numerical methods to approximate Riemann-Liouville derivative such as standard Grünwald-Letnikov formula (first-order accuracy), shifted Grünwald-Letnikov formula (first-order accuracy) [18], L-2 approximation method [19] (first-order accuracy), spline interpolation method [20] (second-order accuracy), weighted and shifted Grünwald-Letnikov formulas [21] (second and third-order accuracy), fractional average central difference formula [22] (second and fourth-order accuracy). It is worth mentioning that high-order algorithms for Riemann-Liouville derivatives were first proposed by Lubich [23], however, the high order algorithms for Riesz derivatives were constructed by Ding and Li [22,24,25].
Many researchers studied the fractional advection-dispersion equation (ADE) recently. Fractional ADE is used for the description of transport dynamics in the complex systems which are controlled by the anomalous diffusion and the non-exponential relaxation patterns [26]. The fractional ADE is also used in groundwater hydrology research to model the transport of passive tracers carried by the fluid flow in a porous medium [27]. Our aim is to investigate the time-space fractional ADE. Time nonlocality deals with memory effects, whereas space nonlocality describes the long-range interaction. The fundamental idea is that fractional order models convey more information about the underlying structure and dynamics of complex systems. Total Shannon spectral entropy for the case of anomalous diffusion governed by a fractional order diffusion equation generalized in space and in time is calculated in [3] as it can be used as a measure of the information content in a fractional order model of anomalous diffusion. This fractional order representation of the continuous time, random walk model of diffusion gives a spectral entropy minimum for normal (i.e., Gaussian) diffusion with surrounding values leading to greater values of spectral entropy. Povstenko et al. [28] examined the fundamental solutions to space-time fractional diffusion equation with mass absorption (mass release) in the case of axial symmetry. Liu et al. [29] considered time fractional ADE and the solution was obtained using variable transformation, Mellin and Laplace transforms, and H-functions. Povstenko and Kyrylych [30] discussed two different generalizations of the space-time fractional advection-diffusion equation. They studied the fundamental solutions to the corresponding Cauchy and source problems for one spatial variable using Laplace transform and Fourier transform with respect to time and spatial coordinate, respectively. Huang and Liu [31] also considered time-space fractional ADE and obtained the solution in terms of Green functions. Meerschaert et al. [18] proposed numerical methods to solve the one-dimensional space fractional ADE with variable coefficients on a finite domain. Tripathi et al. [32] investigated the approximate analytical solution of fractional order nonlinear diffusion equations by using the homotopy analysis method. Momani et al. [33] developed a reliable algorithm using the Adomian decomposition method to construct a numerical solution for the time-space fractional ADE. Liu et al. [34] proposed an approximation of the Lëvy-Feller advection-dispersion process by employing a random walk and finite difference methods. Finite difference methods [35], finite element methods [36], finite volume methods [37], homotopy perturbation methods [38] and spectral methods [39,40] are also employed to approximate the fractional ADE. Furthermore, recent advances in numerical linear algebra had a substantial impact on designing efficient methods for the solution of the resulting linear systems which are dense but whose computational cost can be essentially reduced to O(N log(N)) where N is the size of the underlying coefficient matrix (see [41][42][43][44] and references therein). In this article, we construct a numerical scheme for the time-space fractional ADE by transforming the fractional differential equations into equivalent Volterra integral equations. As it is known that numerical methods for an integral equation have better numerical stability over the schemes designed for equivalent differential equation. Also the numerical methods for an integral equation can be constructed based on the weaker smoothness requirement than that for the differential equation. To the best of our knowledge, all of the other higher order methods are proposed based on the discretizations for fractional derivative directly.
This paper is organized as follows. In Section 2, some useful notations and auxiliary lemmas are introduced. In Section 3, the fractional trapezoidal scheme is derived combined with the second-order fractional weighted and shifted Grünwald-Letnikov formula for the approximation of the Riesz derivative. Section 4 is devoted to the study of the stability and convergence of the proposed scheme. Some numerical experiments are presented to verify the efficiency of our theoretical results in Section 5. The last section concludes our work.

Preliminaries
In this paper, we will consider the following time-space fractional ADE with the initial conditions and the Dirichlet boundary conditions where 0 < α < 1, 0 < β 1 < 1, 1 < β 2 ≤ 2, K β 1 ≥ 0, K β 2 > 0 and c D α t is the Caputo fractional derivative. In addition, ∂ β 1 ∂|x| β 1 and ∂ β 2 ∂|x| β 2 are the Riesz fractional derivatives of order β 1 and β 2 respectively, defined on the domain [0, L] as follows [19] where In the interval [a, b], let x j = jh, (j = 0, 1, ..., M) be mesh points in space, where h = b−a M is the uniform spatial step size. Meerschaert and Tadjeran [18] showed that the standard Grünwald-Letnikov difference formula was often unstable for time dependent problems and they proposed the shifted Grünwald difference operators to approximate the left and right Riemann-Liouville fractional derivatives that have the first order accuracy given by, where p, q are integers, and g (γ) In fact, the coefficients g for all |z| ≤ 1, and they can be evaluated recursively by using the following relation

Lemma 1 ([35]
). Suppose that 0 < β 1 < 1, then the coefficients g Tian et al. in [21] derived the following weighted shifted Grünwald difference operators based on the multi-step method and their Fourier transforms belong to L 1 (R), then the weighted and shifted Grünwald difference operators satisfy uniformly for x ∈ R, where p, q are integers and p = q.
Consider a function f (x) under the same assumptions as in Lemma 3 on the bounded interval [a, b], if f (a) = 0 or f (b) = 0, the function f (x) can be zero extended for x < a or x > b. In addition, then, the γ order left and right Riemann-Liouville fractional derivatives of f (x) at each point x can be approximated with the second order accuracy as follows then b α n = (n + 1) α − n α , (n = 0, 1, 2, ...) satisfy the following properties there exists a positive constant C > 0, such that τ ≤ Cb α n τ α , n = 1, 2, 3, ....

Fractional Trapezoid Formula
In this subsection and in the sequel, the symbol C denotes a generic constant, whose value may be different from one line to another. Integrating both sides of (17) with respect to the time t from t n to t n+1 , and using Lemmas 4, 7 and 9, we obtain where R 2 depends on h and by using Lemma 4, we have Here R 3 depends on τ and h as by Lemma 7, we get Hence, the solution for system (1)-(3) can be approximated by the following scheme: where f n j = f (x j , t n ). We can write this system into the following matrix form: where Let then we obtain the following numerical scheme based on fractional trapezoid formula with initial and boundary conditions

Stability
In this section, we analyze the stability and convergence for the scheme (26).
The spatial errors and convergence orders of the proposed scheme for solving (35) are shown in Tables 1 and 2 for different values of β 1 and β 2 , respectively. Fixing τ = 0.01, α = 0.9, β 2 = 1.6. The L 2 -norm is used to compute the error of the numerical solution at the last time step by Next, we fix the spatial step size h = 0.01 and vary the time step. Table 3 presents the errors and convergence order for various values of α at time T = 1. The numerical convergence order in the spatial and temporal direction is O(τ 2 + h 2 ), as in Theorem 2.

Conclusions
In this article, we have proposed a finite difference method for solving a class of time-space fractional advection-dispersion equation. We combined the trapezoidal formula, which is well known for the numerical integration of Riemann-Liouville integral, with the Grünwald-Letnikov discretization of the Riesz fractional derivative in space to obtain a numerical scheme. We proved that our proposed scheme is stable and convergent with the accuracy of O(τ 2 + h 2 ) under the sufficient conditions 2 α+1 ≤ 3 and √ 17−1 2 ≤ β 2 ≤ 2. However, our numerical experiments given in Tables 2 and 3 depict that when 2 α+1 > 3 and β 2 < √ 17−1 2 the presented numerical method is still stable and convergent for various temporal and spatial time steps. Finally, some numerical experiments for the fractional finite difference method are given that agree very well with our theoretical results.