Finite Difference Approximation Method for a Space Fractional Convection–Diffusion Equation with Variable Coefficients

Space non-integer order convection–diffusion descriptions are generalized form of integer order convection–diffusion problems expressing super diffusive and convective transport processes. In this article, we propose finite difference approximation for space fractional convection–diffusion model having space variable coefficients on the given bounded domain over time and space. It is shown that the Crank–Nicolson difference scheme based on the right shifted Grünwald–Letnikov difference formula is unconditionally stable and it is also of second order consistency both in temporal and spatial terms with extrapolation to the limit approach. Numerical experiments are tested to verify the efficiency of our theoretical analysis and confirm order of convergence.


Introduction
Fractional differential equations (FDE) have attracted the attention of many researchers and scientists due to their importance in different fields of study such as viscoelasticity, fluid mechanics, physics, biology, engineering, and flows in porous media (see [1][2][3][4][5][6] and the references cited therein). As different experiments and implementations have shown, non-integer space derivatives have been used to develop anomalous diffusion to which a particle spreads at a rate inconsistent with the integer Brownian motion problem in the direction of both time and space. When non-integer order is replaced by the second order derivative in a diffusion equation, it acts to enhance the process which we call super-diffusion [7][8][9][10][11][12]. Laboratory experiments and field-scale tracer dispersion breakthrough curves (BTCs) are suitable for exhibiting early time arrivals that are not captured by the integer order derivatives and these non-Fickian phenomena can be controlled by non-classical order convection-diffusion and dispersion equations (FCDE) as it was explained in [13]. To increase the number of applications, there should be significant interest in constructing numerical schemes to solve a well known space fractional convection-diffusion model that has space variable coefficients. In most cases, non-integer order differential problems have no exact solution, so various iterative and numerical approximations [3,9,14] must be pointed out in advance. In general, these kinds of approaches have become important in finding the approximate solutions of fractional differential equations, so extensive numerical methods have been developed for space fractional convection-diffusion equations such as the spectral method [15], finite volume method [16,17], finite difference method [2,9,14,[18][19][20][21][22][23][24][25][26], finite element method [27][28][29][30] and collocation method [31,32].
When the discretization of domain over the region (which belongs to the geometry) is not complex, finite difference approximations are easier and faster than other methods (see [16,33] for further details) to get numerical solutions. In [34], the author used an unconditional stable difference method for time-space fractional convection-diffusion problems with space variable coefficients with first order convergence both in time and space. The Crank-Nicolson finite difference method for one-sided space fractional diffusion equations using an extrapolation method to get second order convergence was studied in [23]. In [9], the explicit and implicit finite difference methods are discussed for a one-sided space fractional convection-diffusion equation with first order convergence in both time and space. A first-order implicit finite difference discretization method for a two-sided space fractional diffusion equation (SFDE) is also applied in [10]. Recently, an unconditionally stable second order accurate difference method for a two-sided time-space fractional convection-diffusion equation was constructed in [35] using the weighted and Shifted Grünwald-Letnikov difference approximation. It is not suitable to apply the weighted combined with shifted Grünwald-Letnikov difference approximation for one-sided Riemann-Liouville fractional derivative to have second order accurate in space. To deal with such issues, it is important to develop a numerical scheme that leads to evaluate a one-sided space fractional convection-diffusion problem. Thus, the main focus of our study is to have temporal and spatial second order convergence estimates for one-sided space fractional convection-diffusion equations based on a stable finite difference method and using spatial extrapolation to the limit approach. The scheme has been treated using the Crank-Nicholson method with the novel Shifted Grünwald-Letnikov difference approximation and the algorithm has been examined both theoretically and experimentally.
Let us consider space-fractional convection-diffusion equation with variable coefficients: with the given initial condition: and homogeneous Dirichlet boundary conditions: where c(x), d(x) and g(x) are continuous functions on [L, R] and p(x, t) is continuous function on is the fluid variable velocity which means the system is evolving in space due to a velocity field and p(x, t) is sink term so that the fluid transport is from left to right. For the case of integer order (α = 2), Equation (1) gives to the classical convection-diffusion equation(CDE). In this study, we have only considered the fractional derivative case which describes a physical meaning in [36] and it involves only a left-sided fractional order derivative. We have assumed that this one-dimensional space fractional convection-diffusion problem has sufficiently smooth and unique enough solutions. The structure of this paper is arranged as follows. In Section 2, we introduce some preliminary remarks, lemmas and definitions and we show the formulation of the new Crank-Nicolson with right Shifted Grünwald-Letnikov difference scheme in Section 3. In Section 4, we describe the unconditional stability using Gerschgorin Theorem and convergence order analysis of the scheme. In Section 5, numerical tests are implemented to show the relevance of our theoretical study and the conclusions are put in Section 6.

Definition 2.
The left hand side and the right hand side fractional order derivatives, respectively, in Equation (1) are the Riemann-Liouville fractional derivatives with order α which are given by:

Definition 3 ([3]
). Let u be given on . The standard Grünwald-Letnikov estimate for 1 < α ≤ 2 with positive order α is defined by the formula, we also define the Grünwald-Letnikov difference operator as: where is called Grünwald-Letnikov coefficient which is the Taylor series expansion ω(z) = (1 − z) α which is the generating function. We can expressed the coefficients by the following recursive relations.
Lemma 1 ([37]). Assume that 1 < α ≤ 2, then Grünwald-Letnikov coefficients ω (α) k satisfy: The Shifted Grünwald-Letnikov difference operator expression is suitable for our purpose because, it allows us to estimate (D α * u) (x), which is defined in Equation (2), numerically in an accurate way. According to [14], right shifted Grünwald-Letnikov difference operator with p shifts for α th order Left R-L fractional derivative of u(x, t), x ∈ [L, R] at x = x m can be expressed as: where Lemma 2 ([38,39]). Let u ∈ C 2n ( ) that has a finite degree of smoothness with (D α + u)(x) which is approximated by h −α ∆ α h u (x) possesses an asymptotic expansion in integer powers of the step-length h, then an expansion in even powers of h for the Shifted operator can be written in the form: Lemma 3 ([39]). Let u ∈ C n+3 ( ) all derivative of u up to the order n + 4 belong to L 1 ( ). Then the Fourier transform of the Grünwald-Letnikov difference operator defined in Equation (5), iŝ Theorem 1. Let u ∈ C 2n+3 ( ) with all derivatives of u up to order 2n + 3 belong to L 1 ( ). For p ≥ 0 define the shifted Grünwald-Letnikov operator: (2), for any computable coefficient a 2k , which is independent of h, u and x, we have uniformly in x ∈ .
Proof of Theorem 1. We closely follow the result described in [9,10] for the unshifted Grünwald-Letnikov formula and also in [23] for the shifted Grünwald-Letnikov formula. We can see that with the Riemann-Lebesgue lemma, the assumptions on u indicates for real positive constant C 1 and from the condition which is imposed on u, we have From Lemma 3 for all t ∈ the Fourier transform for u(x) of the Grünwald-Letnikov approximation is From the definition of Fourier transform, we have observed that for a constant a ∈ , we have: have the Taylor expansion where a 2k =(−1) 2k ( α a 2k )=( α a 2k ), converges absolutely for |z| ≤ 1 since the function ω α,p (z) is bounded on . The shifted Grünwald difference approximation (∆ h,p )u(x) ∈ L 1 ( ). Thus, we have since ω α,p (−ith) is analytic around the origin, we express it as an even power expansions which absolutely convergent for all |z| ≤ R. For this a bounded function ω α,p (z) on , there exist a real positive constant C 2 which satisfy: is bounded uniformly in x ∈ . For any value |x| ≤ R , we have which is bounded on . For the other case |x| > R also, we have where Then, this implies that Equation (15) holds for all x ∈ . From Equation (17), we can write and with the conditions imposed on u, we can say that (1 + |x| 2n+3ũ (t) is bounded on . Thus, |t| 2α−3 |ũ(t)| ∈ L 1 (R). This implies that, for k ∈ with C = C 1 C 2 . Therefore using the Fourier inversion transform, we have At last, we have (10), it can be seen that for p = α/2, the error takes its minimum value and a second order convergence is achieved. We need the grid points x m − (k − p)h to find an optimal positive integer p that makes p − α/2 is minimum. It is numerically proved in [3] that for the value 0 < α ≤ 1, p = 0 is acceptable; while for 1 < α ≤ 2, p = 1 is optimal.

Remark 1. From Equation
Remark 2. Theorem 1 is the base of Extrapolation to the limit. Therefore one can apply it the Shifted Grünwald-Letnikov difference operator to obtain the convergence rate with arbitrary high order h k , k = 1, 2, 3, ..., n such that

Problem Formulation of the Scheme
Consider the following one-dimensional space fractional convection-diffusion problem: which is based on shifted Grünwald-Letnikov difference method with 1 < α ≤ 2 on a finite domain L < x < R.

Crank-Nicolson Scheme for Time and Shifted Grünwald Difference Scheme for Space Discretization
We partition the finite interval [L, R] with a uniform mesh in the space size step h = (R − L)/N x and the time step τ = T/N t , in which N x , N t are non-negative integers and the set of grid size points is symbolized by x m = mh and t n = nτ for 0 We use the following notations: . Applying the C-N technique for the time discretization of Equation (20) gives to In space discretization we have used the central finite difference method for the convection term and the Shifted Grünwald-Letnikov operator for the space fractional derivative with the approach of spatial Extrapolation to the limit, respectively. See the full discretization of the scheme: Multiplying Equation (22) by τ the discretization equation, we have The above equation is used to predict the values of u(x, t) at time n + 1, so all the values of u at time n are assumed to be known. For simplification Both the convection and diffusion variable coefficients are (N x − 1) × (N x − 1) diagonal matrices which are defined by These discretization together with Dirichlet boundary conditions which results in a linear system of equations for which the coefficient matrix is the sum of lower triangular and upper-diagonal matrices. The above discretization can be re-arranged to yield: Denoting U n m as the numerical approximation of u n m , we can construct the C-N scheme for Equation (20) 1 I is the (N x − 1) × (N t − 1) identity matrix with A m,n as the matrix coefficients. These coefficients,for m = 1, 2, 3, ..., N x − 1, n = 1, 2, ..., N t − 1 are given by: The finite difference scheme (24) and (26) defines a linear system of equations as

Theorem 2.
Suppose that 1 < α ≤ 2, the coefficient matrix defined in Equations (24)-(27), then the diagonal matrix and the coefficient matrix satisfy: Proof of Theorem 2. As we have seen from the coefficient matrix defined in Equation (27), This implies that ∑ N x −1 n=0,m =1 |A m,n | < A m,m . Therefore, the diagonal matrix is strictly dominant.

Theoretical Analysis of Finite Difference Scheme
In general for analyzing convergence and stability, we consider the following description.

Boundedness of the Fractional Scheme
The Classical Crank-Nicolson scheme combines the stability of an implicit finite difference method with its accuracy which produce second order convergence in both space and time.
Theorem 3. Crank-Nicolson scheme for solving space fractional convection-diffusion equations given by the following problem: which is based on shifted Grünwald-Letnikov difference approximation scheme is bounded for 1 < α ≤ 2.
Proof of Theorem 3. Consider C-N scheme for the space-fractional convection-diffusion problem for Here, we have shown the convergence and boundedness of the scheme by taking the smaller time-step in terms of Lax-Richtmyer stability analysis that uses a weaker bound (see [40]). Our matrix A has an eigenvalues of λ that have positive real parts, and, we also have found a strictly dominant matrix. These eigenvalues which are centered in the disks at each diagonal entries as: From the Gerschgorin Theorem in [41], the radius of the matrix can be expressed as Since from the Grünwald coefficients we have ω For a bounded ratio of time-step τ and space-step h with nτ ≤ T, we have

From the relation of Parseval's Theorem, [40]
which shows that the scheme is bounded.

Stability Analysis
Theorem 4. Let U n m be the numerical approximation of the exact solution u n m , then the C-N finite difference scheme (28) is unconditionally stable.

Proof of Theorem 4.
Consider the matrix coefficient of the difference approximation for the problem (20) can be written as described above Let e n = e n 1 , e n 2 , e n 3 , ..., e n N x −1 , and take the relation between the error e n+1 in U n+1 and the error e n in U n which is given by the linear system First of all, we must show that the (non-real valued) eigenvalues of the coefficient matrices A have positive real parts. For ω (α) 1 = −α with fractional order 1 < α < 2 and k = 1; we have ω (α) k > 0. In addition to this, −ω α 1 = α ≥ ∑ N k=0,k =1 ω α k for the value N > 1. As stated in Gerschgorin Theorem ( [41], pp. 136-139), the eigenvalues of the given matrix A are inside the disks centered at each diagonal entry.
These Gerschgorin disks are belong to the right half of the complex plane. Thus, the eigenvalue of the coefficient matrix A has positive real part which implies that A has an eigenvalue λ if and only if (I − A) has an eigenvalue (1 − λ) if and only if (I + A) −1 (I − A) has an eigenvalue 1−λ 1+λ . From the first part of this sentence, we have seen that all the eigenvalues of the matrix given by (I + A) have a radius larger than unity which implies the matrix is invertible. Now we can see from the above description the real part of λ is non-negative which we can conclude that Thus, the spectral radius of the system matrix (I + A) −1 (I − A) is strictly less than unity which implies that the difference scheme is unconditionally stable.

Convergence Analysis
First of all we have given the Truncation error of the C-N scheme. It is obvious to conclude that: From the above Extrapolation to the limit Theorem for n = 1, we got Therefore the local truncation error of (20) is given by T n+1 m = O(τ 2 + τh) Theorem 5. Let u n m be the exact solution of problem (20), and U n m be the solution of the finite difference scheme (26), then for all 1 ≤ n ≤ N t , we have the estimate where u n m − U n m ∞ =max 1≤m≤N x |u n m − U n m |, c is a non-negative constant independent of h and τ with ||.|| stands for the discrete L 2 -norm.

Remark 3.
The Crank-Nicolson scheme, for classical convection-diffusion equation, provides stable C-N finite difference method that is second order convergence in time and space. Also a study based on C-N finite difference method with the spatial extrapolation to the limit method, see Theorem 1, is used to get temporal and spatial second order for one-sided SFCDEs with space variable coefficients.

Numerical Tests
Problem test 1 1. Consider the space-fractional diffusion type of problem: homogeneous Dirichlet boundary condition The exact solution is All numerical experiments are implemented using Theorem 1 and C-N scheme with the space domain, 0 < x < 1 and time domain, 0 < t < T. Figure 1 shows the maximum error produced by C-N scheme for large enough time domain and numerical solution is close enough to the exact solution using C-N scheme with α = 1.5 in Figure 2. The maximum error and second order convergence for the fractional diffusion and fractional convection-diffusion equation with variable coefficients are given in Tables 1-3.
homogeneous Dirichlet boundary condition with variable convection-diffusion coefficients respectively, and source term The exact solution is u(x, t) = e −2t (x α − x) Figures 3 and 4 show the numerical and exact solutions for fractional diffusion and fractional convection-diffusion problems with large enough time domain in example 1 and 2, respectively.The exact and numerical solution of fractional convection-diffusion equation by C-N scheme is also given in Figure 5. In Table 4, the maximum error and first order convergence in space is obtained using C-N scheme without extrapolation to the limit approach by fixing the time step.    The exact solution is Problem test 4 is experimented with the grid size reduction extrapolation approach stated in [23]. We have smooth enough numerical and exact solutions by using C-N scheme in Figure 6, and Table 5 shows the maximum error with the error rate is given for space fractional convection-diffusion problem with a grid size reduction extrapolation method.

Conclusions
The one dimension space fractional diffusion and fractional convection-diffusion problem with space variable coefficients is solved by the fractional C-N scheme based on the Extrapolation to the limit approach of right shifted Grünwald-Letnikov approximation. The fractional C-N method, for the fractional diffusion problem and fractional convection-diffusion equation with space variable coefficients, is consistent and unconditionally stable with second order convergence. Numerical examples confirmed that the C-N method is suitable for the space fractional convection-diffusion problem even for a large value of time domain.
Author Contributions: The authors contributed equally to the writing and approved the final manuscript of this paper. All authors have read and agreed to the published version of the manuscript.