Finite Difference Method for Two-Sided Two Dimensional Space Fractional Convection-Diffusion Problem with Source Term

: In this paper, we have considered a numerical difference approximation for solving two-dimensional Riesz space fractional convection-diffusion problem with source term over a ﬁnite domain. The convection and diffusion equation can depend on both spatial and temporal variables. Crank-Nicolson scheme for time combined with weighted and shifted Grünwald-Letnikov difference operator for space are implemented to get second order convergence both in space and time. Unconditional stability and convergence order analysis of the scheme are explained theoretically and experimentally. The numerical tests are indicated that the Crank-Nicolson scheme with weighted shifted Grünwald-Letnikov approximations are effective numerical methods for two dimensional two-sided space fractional convection-diffusion equation.


Introduction
Differential equation described fractional partial differential equations are appropriated to explain complex problems like viscoelasticity, electroanalytical chemistry, biology, fluid mechanics, engineering [1], physics [1,2], fractional operators [3] and flows in porous media [4][5][6][7][8]. Through the advection and dispersion processes, pollutants create a contaminant plume within an aquifer, the movement of which in an aquifer is described by transport model. One of the very rich transport model is advection-dispersion model, which is used to describe the transport phenomena in different fields of science. Solute transport is important to predict the solute concentration in aquifers, rivers, lakes and streams too.
The remain arrangement of this paper is organized as follows-in Section 2, we introduce some preliminary remarks, lemmas and definitions. We have shown the formulation of one dimensional Riesz space fractional convection-diffusion problem with Crank-Nicolson and weighted shifted Grünwald-Letnikov difference scheme in Section 3. In Section 4, we have described the formulation with discretization of two-dimensional Riesz space fractional convection-diffusion problem. In Section 5, unconditional stability and convergence order analysis of the scheme have done using CNADI-WSGD. In Section 6, numerical simulations are implemented to show the importance of our theoretical study and the conclusions are discussed in Section 7.

Preliminary Remarks
Definition 1. The Riesz differential operator which is given by analytic continuation in the whole range 0 < α ≤ 2 with α = 1 as: where the fractional derivative, with n ∈ N and coefficient K = −1 2cos(απ/2) , are the left and right Riemann-Liouville fractional derivatives. From this definition the fractional integral operators −∞ I α x u(x, t) and x I α ∞ u(x, t) are the left and right Weyl fractional integrals as defined in Reference [44]: Lemma 1 ([44,45]). Let α > 0 and Γ(.) represents gamma function, then the following are properties of binomial coefficients: 2.
Theorem 1. Let u(x) has n − 1 continuous derivatives on the closed interval [a, b] with the derivatives u (n) (x) are integrable for x ≥ a or x ≤ b, then for each α(n − 1 < α ≤ n), the left and right Riemann-Liouville fractional derivatives exist and coincide with the corresponding (left and right) Grünwald-Letnikov fractional derivatives.
Proof. The left standard Grünwald-Letnikov fractional derivative is given by the limit expression on [a, x], where x−a n = h = b−x n , n − 1 < α ≤ n. Here our aim is to evaluate the limit described in Equation (6). For the evaluation of the limit, we are assuming the function u(x) continuous on [a, x] and for α > 0, we have: where . We need to transform Equation (7) to the following form using the property 1 of Lemma 1.
Similarly, we have to apply property 1 of Lemma 1 repeatedly m times, after simplification we get: Now, we need to evaluate the limit of Equation (9).
, denote the second sum. Let us find the limit of p th -term of the first sum.
In order to evaluate the limit of second sum U h s , we have to follow the property of binomial coefficients of Lemma 1.
From property 4 of Lemma 1, we have By considering Equations (12) and (13) we have that: (14) Now by combining Equations (10) and (14), we have finalized the general limit evaluation as: By taking n = m + 1 or n − 1 = m with n − 1 < α ≤ n, the left Grünwald-Letnikov fractional derivative over the closed interval [a, x] is written as: Similarly, the right standard Grünwald-Letnikov fractional derivative on the closed interval Thus, for a → −∞, u (p) (x) approaches to zero and Equation (16) leads to have: which gives the left Riemann-Liouville fractional derivative as we are expected and it is also exists for n − 1 < α ≤ n. In a similar proof, we also have the right Riemann-Liouville fractional derivative as b → ∞: Left Riemann-Liouville fractional derivative: Right Riemann-Liouville fractional derivative: As it is discussed in Reference [46], the shifted Grünwald-Letnikov difference operator with first order ∆ (α) which is defined as, approximates the left and right Riemann-Liouville fractional derivatives. Here p is an integer that shifts the approximation p-shift to the right and g α k = (−1) k α k are the coefficients of the power series for the function (1 − z) α , for all |z| ≤ 1,with: Lemma 2 ([47]). The coefficients g (α) k satisfy the following properties for 0 < α < 1.
Applying the above Theorem 1, and weighted shifted Grünwald-Letnikov fractional derivative derivation from Reference [46] for 0 < α < 1, 1 < β ≤ 2, the left and right Riemann-Liouville fractional derivatives of u(x) over a bounded interval at each point x can be formulated as: where ω (α) 0 The properties of the weighted coefficients ω

Numerical Approximation for One Dimensional Two-Sided Convection-Diffusion Problem with Source Term
We have considered the one-dimensional two-sided space fractional convection-diffusion equation, with initial condition: and with zero Dirichlet boundary conditions: The analytic solution for Riesz space fractional convection-diffusion equation is developed in Reference [49] using the spectral representation on a finite interval [0, L]. Reference [50] used Laplace transform and Fourier transform method for finding analytical solution of Riesz space fractional convection-diffusion problem with initial and zero Dirichlet boundary conditions. Here our discretization is based on the finite interval [0, L] into a uniform mesh with the space step h = L/N x and the time step τ = T/N t , where N x , N t are positive integers and the set of grid points is denoted by x m = mh and t n = nτ for 0 ≤ m ≤ N x and 0 ≤ n ≤ N t . Let We have used the following notations for our formulation: The Riesz space fractional convection-diffusion equation for 0 < α < 1, 1 < β < 2 can be written with following expression.
Theorem 1 allows us to use the Riemann-Liouville fractional derivative definition for the formulation of the problem. The weighted shifted Grünwald-Letnikov derivative formula for approximating the two-sided fractional derivative is derived in References [46,48] for space fractional derivative and Crank-Nicolson scheme for time are used. where Assume U n m be the numerical approximation of the solution u n m , then the CN-WSGD formulation for RSFCDEs become: By denoting, Therefore, the system of equations takes the form: where I is the (N x − 1) × (N t − 1) identity matrix with A m,j as the matrix coefficients. These matrix coefficients for m = 1, 2, 3, ..., For the convenience of implementation, using the matrix form of the grid functions,

Formulation and Discretization of Two-Dimensional Riesz Space Fractional Convection Diffusion Equation with CNADI-WSGD Scheme
The analytic solution for two-dimensional Riesz space fractional anomalous diffusion equation is obtained by using the Fourier series expansion with homogeneous Dirichlet boundary condition. Let us take a bounded domain as Here our aim is to find the full numerical approximation of the two-dimensional Riesz space fractional convection-diffusion problem with zero Dirichlet boundary condition over a finite domain Ω × Ω t .
Consider the two-dimensional two-sided space fractional convection-diffusion problem with constant coefficients as: where 0 < α 1 , α 2 < 1, 1 < β 1 , β 2 < 2, c x , c y ≥ 0 and d x , d y > 0 express the velocity parameter and positive diffusion coefficients. Here the function u(x, y, t) is specified as solute concentration under the groundwater. The Riesz space fractional-order derivative is defined as: where and also from the coincides Theorem 1, we have the following left and right Riemann-Liouville fractional derivative definition for two dimension space fractional derivative.
To simplify our formulation, it is possible to symbolize the following operator as: By grouping like terms from Equations (43) and (44), we have: where T n m,j represent truncation error that can satisfy T n m,j ≤k τ 2 + h 2 1 + h 2 2 . Let us define the operators: with these operator definitions, the CNADI-WSGD scheme for the 2D-RSFCDE with homogeneous Dirichlet boundary conditions can be defined as an operator form: An alternating direction implicit Peacemann-Rachford is reduced a two-dimensional problem in to a one dimensional problem with a better computational efficient. For CNADI the operator can be expressed in the product form as: which produce an additional perturbation error in the form of τ 2 4 ∆ α x ∆ β y u n+1 m,j − u n m,j that has Taylor expansion as: As compared to the approximation errors, the additional perturbation errors is insignificant and the scheme defined in Equation (45) has second order accuracy in both space and time which is O τ 2 + h 2 1 + h 2 2 . The problem defined by Equation (47) can be simulated by the following efficient Peacemann-Rachford ADI approximation as it was presented in Reference [46] by considering u * m,j as an intermediate solution to make a numerical solution u n m,j at time t n to the numerical solution u n+1 m,j at time t n+1 . The corresponding iterative algorithms are: Algorithm 1: The first step is to solve the problem in the x-direction for each fixed y j to find an intermediate solution u * m,j in the form: Algorithm 2: The next step is to solve the problem in y-direction for each fixed x m as: Algorithm 3: We need to apply the homogeneous Dirichlet boundary conditions: Therefore now compute the boundary condition for the intermediate solution u * m,j which can be derived from subtracting Equation (50) from (49) to get: Therefore, the boundary conditions for u * m,j needed to solve each set of equations.
By setting U n m,j be the numerical approximation to exact solution u n m,j , we get the finite difference approximation for Equation (47):

Stability and Convergence Analysis of CNADI-WSGD Scheme
For discussing the stability and convergence of the scheme, we need to write our problem in matrix form. Thus, the Equation (49) can be put as: In a similar way, Equation (50) can be given in matrix form: where,ū n+1 q = u n+1 q,1 , u n+1 q,2 , ..., u n+1 Theorem 2. Assume that 0 < α 1 , α 2 < 1, 1 < β 1 , β 2 ≤ 2 , the coefficient matrices defined in Equations (55) and (57), then the diagonal matrix and coefficient matrix satisfy: tells us that A and B which are defined in Equations (54) and (56) are strictly diagonally dominant.
Proof. First we will consider the diagonal dominance of the coefficient matrix a m,j . Since K α 1 = 1 2cos(πα 1 /2) > 0 and K β 1 = 1 2cos(πβ 1 From Lemmas 4 and 5, we have: Sincec x > 0 andd x < 0, then we have: By looking Lemmas 4 and 5, we have seen that ω (α 1 ) 1 > 0 and ω As we have shown from Lemma 4, when k ≥ 3, ω a m,j , m = 1, 2, 3, ..., N x − 1, which means matrix A defined by the coefficient matrix a m,j , is strictly diagonally dominant. In the same way, the diagonally dominant result for matrix B can also be found as matrix A.

Stability Analysis of the CNADI-WSGD Method
In order to study the stability and convergence analysis for the CNADI-WSGD scheme , we are focused on the following description.
Let χ h = ν : ν = ν m,j : x m = mh 1 ; y j = jh 2 N x m=0 N y j=0 be the mesh grid function. For any ν = ν m,j ∈ χ h , we define our point-wise maximum norm as: and the discrete L 2 -norm Our next aim is to show the stability of CNADI-WSGD method which is defined as in the matrix form: where the matrices A and B define the operator τ 2 ∆ (α) x and y , respectively. The vector T n+1 absorbs the source term p n+1/2 m,j and the Dirichlet boundary condition in the formulated problem. Theorem 3. Let U n m,j be the numerical solution of the exact solution u n m,j , then CNADI-WSGD finite difference method (53) is unconditionally stable for 0 < α 1 , α 2 < 1 with 1 < β 1 , β 2 ≤ 2.
Proof. The matrices A and B are of size (N x − 1) N y − 1 × (N x − 1) N y − 1 . The commutative property defined in Reference [51], allows us to obtain the unconditional stability of CNADI-WSGD method. The matrix A which is N y − 1 × N y − 1 block diagonal matrix whose blocks are (N x − 1 × N x − 1) square super triangular matrices which is expressed as A = diag A 1 , A 2 , ..., A N y −1 . In the same way, the matrix B is a block matrix with (N x − 1) × (N x − 1) square diagonal matrices. The matrix B can be written as is the (m, j) th entry of the matrix B defined above. As we have seen from Theorem 2 , matrix A is diagonal dominant with entry a m,m > 0. The sum of the absolute value of the off-diagonal entries on the row m of matrix A is: implies that, Next we need to show that the eigenvalue of matrix A is negative real parts. For 0 < α 1 < 1, 1 < β 1 < 2, we can see that, We have noticed that, According to Greschgorin Theorem [52], the given eigenvalue of matrix A have non-positive real parts. Here we have noted that matrix A has an eigenvalue of λ 1 if and only if (I − A) has an eigenvalue of (1 − λ 1 ) if and only if (I − A) −1 (I + A) has an eigenvalue of (1 + λ 1 )/(1 − λ 1 ). From the first part of this statement, we can concluded that all eigenvalues of the matrix (I − A) have a spectral radius which is larger than unity indicates the matrix is invertible. Thus, every eigenvalue of the (I − A) −1 (I + A) has a spectral radius which is less than 1. Similarly, we can show that matrix B also satisfy the same property as matrix A. From the scheme (53), we can express the error e n+1 in U n+1 at time t n+1 and the error e n in U n at time t n as: where the identity matrix I is (N x − 1) × N y − 1 square. Hence, Equation (53) can be put in the form: Letting λ 1 and λ 2 be an eigenvalue of matrices A and B respectively, then it results from Equation (69) that the real parts of λ 1 and λ 2 are both negative. The spectral radius of each matrix is less than unity, which has followed that (I − A) −1 (I + A) n and (I − A) −1 (I + A) n which converges to null matrix (see Reference [46]). Therefore, we have concluded the scheme defined in Equation (53), is unconditionally stable.

Convergence Analysis of CNADI-WSGD Scheme
First of all we can express the truncation error of CNADI-WSGD difference method. So, it is easy to conclude that: where 0 < α 1 , α 2 < 1. It is the same to have a truncation error of O(τ 2 ) and O h 2 1 + h 2 2 for 1 < β 1 , β 2 < 2.
Therefore, the truncation error from Equation (43) is given by: Theorem 4. Assume u n m,j be the analytic solution, and let U n m,j be the approximation solution of the finite difference method (65), then for all 1 ≤ n ≤ N t , we have the estimate: where ||u n m,j − U n m,j || ∞ = max 1≤m≤N x ,1≤j≤N y |u n m,j − U n m,j | = |e nm ,ĵ |, C is a positive constant independent of h 1 , h 2 and τ with ||.|| stands for the discrete L 2 -norm.

Numerical Simulations
1. Consider the one dimensional RSFCDEs over a bounded domain with initial and Dirichlet boundary conditions: with the source term: The exact solution is u(x, t) = t β e αt x 2 (1 − x) 2 .
All the numerical simulations are done based on the finite space domain Ω × Ω t where Ω = [0, 1] × [0, 1] and Ω t = [0, 1]. The order of convergence both in space and time are calculated using the formula: (2) , where order 1 is the rate of convergence for one-dimensional two-sided space fractional convection-diffusion equation and order 2 is rate convergence of two dimensional two-sided space fractional equation. ||E(h, τ)|| ∞ is the maximum error for one dimensional space fractional problem and ||E(h x , h y , τ)|| ∞ for two dimensional space fractional problem, denoted as Max − Error. As we have seen in Table 1, the second order convergence and the maximum error are confirmed at each grid size for convection-dominance (i.e, c x > d x ) for one dimensional two-sided space fractional convection-diffusion equation with different space fractional order. As we have refined the grid size, the suitable maximum error is obtained. The convergence order and maximum error for a diffusion-dominance (i.e, c x < d x ) one dimensional two-sided space fractional convection-diffusion problem are shown in Table 2. Figure 1 shows the good agreement of exact and numerical solution of one-dimensional convection-diffusion equation with the coefficients c x = 0.5, d x = 1.5 and with fractional orders α = 0.75, β = 1.85 at N x = N t = 100 grid points.  2. Consider two-dimensional diffusion problem.(c x = c y = 0) with the source term: Table 3 shows that the maximum error and order of convergence for two-dimensional two-sided space fractional diffusion equation with different space fractional orders by taking c x = 0 = c y . For this numerical simulation, we have used same step-size for space and time (i.e, h x = h y = τ). The maximum time domain that used to obtain all the numerical results is T = 1 and the diffusion coefficients are d x = 2 = d y . The surface plot of u(x, y, t) with the diffusion coefficients d x = 2.5, d y = 1.5, β 1 = 1.25, β 2 = 1.85 at the mesh points h 1 = h 2 = τ = 0.01 is given in Figure 2. Table 3. Convergence rate and maximum error produced for example 2 at T = 1, d x = 2 = d y , h x = h y = τ. 3. Let us consider the two-dimensional Riesz space fractional convection-diffusion problem with bounded domain: with the source term: The exact solution is, u(x, y, t) = t β e αt x 2 (1 − x) 2 y 2 (1 − y) 2 .
In Table 4, we have found a numerical results that produce second order convergence rate and maximum error for two sided two dimensional space fractional convection-diffusion equation with diffusion-dominance (c x = 0.25 = c y , d x = d y = 2) phenomena. For this simulation we have taken a fixed value for β 2 (β 2 = 1.75) and for α 2 (α 2 = 0.5) with different values for α 1 , β 1 . Similarly in Table 5, we have considered the convection-dominance (c x = 2 = c y , d x = d y = 0.25) two-sided two-dimensional space fractional convection-diffusion problem with fixed β 2 (β 2 = 1.75) and for fixed α 2 (α 2 = 0.5). The order of convergence and maximum errors are calculated using the formula expressed in Equation (79). In Figure 3 the surface plot of exact and numerical solutions for two-dimensional convection-diffusion equations are investigated by considering the coefficients c x = c y = 2.5, d x = d y = 1.5 with orders α 1 = 0.75, α 2 = 0.75, β 1 = 1.85, β 2 = 1.85 at N x = N y = N t = 100 mesh grid points. Table 4. Convergence order produced with diffusion-dominance for example 3 at T = 1, c x = 0.25 = c y , d x = 2 = d y , h x = h y = τ.

Conclusions
In our study, we have developed an algorithm for two-dimensional two-sided space fractional convection-diffusion problem using the CNADI difference method for time discretization combined with WSGD scheme for the approximation of space fractional derivative. We have used a shifted category of standard Grünwald-Letnikov difference method and weighted version of the shifted Grünwald-Letnikov difference approximation with CNADI scheme to have unconditionally stable and second order convergence both in space and time without extrapolation. Moreover, unconditional stability and second order convergence is justified for convection-dominance two-sided two dimension space fractional convection-diffusion equation. Our theoretical study and analysis has been confirmed by our numerical simulation in Section 6. We will consider the space fractional reaction convection-diffusion equation in our near further research.
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 final version of the manuscript.