New group iterative schemes in the numerical solution of the two-dimensional time fractional advection-diffusion equation

Abstarct: Numerical schemes based on small fixed-size grouping strategies have been successfully researched over the last few decades in solving various types of partial differential equations where they have been proven to possess the ability to increase the convergence rates of the iteration processes involved. The formulation of these strategies on fractional differential equations, however, is still at its infancy. Appropriate discretization formula will need to be derived and applied to the time and spatial fractional derivatives in order to reduce the computational complexity of the schemes. In this paper, the design of new group iterative schemes applied to the solution of the 2D time fractional advection-diffusion equation are presented and discussed in detail. The Caputo fractional derivative is used in the discretization of the fractional group schemes in combination with the Crank–Nicolson difference approximations on the standard grid. Numerical experiments are conducted to determine the effectiveness of the proposed group methods with regard to execution times, number of iterations, and computational complexity. The stability and convergence properties are also presented using a matrix method with mathematical induction.The numerical results will be proven to agree with the theoretical claims. *Corresponding author: Alla Tareq Balasim, School of Mathematical Sciences, Universiti Sains Malaysia, Penang 11800, Malaysia; Department of Mathematics, College of Basic Education, University of Al-mustansiria, Baghdad, Iraq E-mail: alkhazrejy@yahoo.com


PUBLIC INTEREST STATEMENT
Many phenomena in science and engineering can be modelled as two dimensional fractional advection-diffusion equations satisfying appropriate initial and boundary conditions. With the advent of computing technology, effective numerical methods have been extensively formulated in solving these equations due to their simplicity and accuracy. In particular, finite difference formulas, are oftenly used to numerically discretize the differential equations to produce sparse systems of linear equations which may be suitably solved by iterative solvers. However, iterative solvers have the disadvantage of being too slow to converge which may increase the computation timings. To overcome this problem, suitable small fixed-size grouping strategies are applied to the mesh points of the solution domain following the discretization of the differential equation which results in schemes with faster convergence and lower computational complexity with comparable accuracies. This is very useful to engineers and scientists who are involved in time consuming simulation modeling processes.

Introduction
Fractional calculus, which is the calculus of integrals and derivatives in random order, dates back as far as the more popular integer calculus, and has been gaining significant interest over the past few years with its history and development being explored in detail by Oldham and Spanier (1974), Miller and Ross (1993), Samko, Kilbas, and Marichev (1993) and Podlubny (1998). Fractional differential equations (FDEs) can be used to model many problems in a wide field of applications. They are defined as equations that utilize fractional derivatives and considered as powerful tools that can describe the memory and hereditary characteristics of various materials. Several researchers have explored the use of FDEs in the fields of chemistry (Gorenflo, Mainardi, Moretti, Pagnini, & Paradisi, 2002;Seki, Wojcik, & Tachiya, 2003), physics (Henry & Wearne, 2000;Metzler, Barkai, & Klafter, 1999;Metzler & Klafter, 2000;Wyss, 1986) and other scientific and engineering spheres (Baeumer, Benson, Meerschaert, & Wheatcraft, 2001;Benson, Wheatcraft, & Meerschaert, 2000;Cushman & Ginn, 2000;Mehdinejadiani, Naseri, Jafari, Ghanbarzadeh, & Baleanu, 2013). Since there are at most no exact solutions to the majority of fractional differential equations, it is necessary to resort to approximation and numerical methods (Abdelkawy, Zaky, Bhrawy, & Baleanu, 2015;Balasim & Ali, 2015;Baleanu, Agheli, & Al Qurashi, 2016;Bhrawy & Baleanu, 2013). Over the past decade, there has been an influx of numerical methods development for solving various types of FDEs (Agrawal, 2008;Ali, Abdullah, & Mohyud-Din, 2017;Chen, Deng, & Wu, 2013;Chen, Liu, Anh, & Turner, 2011;Chen, Liu, & Burrage, 2008;Leonenko, Meerschaert, & Sikorskii, 2013;Li, Zeng, & Liu, 2012;Liu, Zhuang, Anh, Turner, & Burrage, 2007;Shen, Liu, & Anh, 2011;Sousa & Li, 2015;Su, Wang, & Wang, 2011;Uddin & Haq, 2011;Zhang, Huang, Feng, & Wei, 2013;Zheng, Li, & Zhao, 2010;Zhuang, Gu, Liu, Turner, & Yarlagadda, 2011;Zhuang, Liu, Anh, & Turner, 2009).  used implicit and explicit difference techniques to solve time fractional reaction-diffusion equations, while Liu et al. (2007) proposed for these techniques to be employed in solving the space-time fractional advection-dispersion equation by replacing the first-order time derivative by the Caputo fractional derivative, and the first-order and second-order space derivatives by the Riemman-Liouville fractional derivatives. The use of radial basis function (RBF) approximation method was discussed in Uddin and Haq (2011) to solve the time fractional advection-dispersion equation in a bounded domain. The application of finite element method was also seen in Zheng et al. (2010), in solving the space fractional advection-diffusion equation under non-homogeneous initial boundary conditions. In Chen et al. (2011), a numerical method with first-order temporal accuracy and second-order spatial accuracy was established for solving the variable order Galilei advection-diffusion equation with a nonlinear source term, while Zhuang et al. (2009) used the explicit and implicit Euler approximation to solve a variable order fractional advection-diffusion equation with a nonlinear source term. A new numerical solution for the 2D fractional advection-dispersion equation with variable coefficients in a finite field was also introduced by . In 2011, Shen et al. (2011) proposed the explicit and implicit finite difference approximations for the space-time Riesz-Caputo fractional advection-diffusion equation where the implicit scheme was proven to be unconditionally stable, while the explicit scheme exhibits a conditionally stable property. The development of an implicit meshless method was also seen in Zhuang et al. (2011) in solving the time-dependent fractional advection-diffusion equation where a discretized system of equations was obtained using the moving least squares (MLS) meshless shape functions.
In general, finding numerical solutions to FDEs using iterative finite difference schemes is not a straightforward process and can be a challenging task due to several reasons. Firstly, all the earlier solutions have to be saved if the current solution is to be computed, making the calculations even more complex and costly in terms of CPU usage time in cases where traditional point implicit difference approaches are used. Secondly, although the point implicit techniques are stable, a significant amount of CPU usage time is required when many unknowns are involved. In recent decades, grouping strategies have been proven to possess characteristics that are able to reduce the spectral radius of the generated matrix resulting from the finite difference discretization of the differential equation, and therefore increase the convergence rates of the iterative algorithms. They have been shown to reduce the computational timings compared to their point wise counterparts in solving several types of partial differential equation (Ali & Kew, 2012;Evans & Yousif, 1986;Kew & Ali, 2015;Ng & Ali, 2008;Othman & Abdullah, 2000;Tan, Ali, & Lai, 2012;Yousif & Evans, 1995). These group methods are also easy to implement and suitable to be implemented on parallel computers due to their explicit nature. However, till date these strategies have not been tested on solving FDEs, particularly for the 2D cases. Therefore, in this study, the formulation of new group iterative methods is presented in solving the following 2Dl time fractional advection diffusion equation where 0 < < 1, a x , a y , b x , b y are positive constants and f(x, y, t) is nonhomogeneous term subjected to the following initial and Dirichlet boundary conditions This equation plays an important role in describing transport dynamics in complex systems which are governed by anomalous diffusion and non-exponential relaxation patterns (Zhuang, Gu, Liu, Turner, & Yarlagadda, 2011). This paper is outlined as follows. Section 2 presents the proposed group iterative methods obtained from the Crank-Nicolson difference approximation followed by the stability and convergence analysis of the difference schemes in Sections 3 and 4, respectively. Section 5 presents the discussion on the computational effort involved in solving Equation (1) using the proposed methods with regard to the arithmetic operations for each iteration. Finally, the results of the numerical experiments are presented and discussed in Section 6.

Standard approximation schemes for fractional advection-diffusion equations
The Caputo fractional derivative, D , of the order-is expressed as follows (Li, Qian, & Chen, (2011): where Γ(⋅) is the Euler Gamma function. Further details on the definitions and properties of fractional derivatives are available in Podlubny (1998).
We need to apply appropriate finite difference approximations to the time and space derivaties of (1), let h > 0 be the space step and k > 0 be the time step. The domain is assumed to be uniform in both x and y directions. Define x i = ih, y j = jh, {i, j = 0, 1,..., n}, and a mesh size of h = 1 n , where n is an arbitrary positive integer and t k = k , k = 0, 1,..., l. various approximations formulas could be obtained for (1) at the point of (x i , y j , t k ). By taking the average of the central difference approximations to the left side of (1) at the points (i, j, k + 1) and (i, j, k), the Caputo time fractional approximation (2) can be transformed to the following form (Karatay, Kale, & Bayramoglu, 2013) Using (3) in combination with the second-order Crank-Nicolson difference scheme for the right side of (1) will result in the following fractional standard point (FSP) formula Equation (4) can be simplified to become as follows 1-w1-sx-sy  Figure 1 shows the computational molecule for Equation (5).

Fractional explicit group method
In formulating the fractional explicit group (FEG) method, (5) is applied to any group of four points in the solution domain to generate a 4 × 4 system of equation as follows Mathematical software can be used to easily invert (6) to obtain the FEG formula where Figure 2 shows the construction of blocks of four points in the solution domain for the case n = 10. Note that if n is even, there will be ungrouped points near the upper and right sides of the boundary. The FEG method proceeds with the iterative evaluation of solutions in these blocks of four points using Equation (7) throughout the whole solution domain until convergence is achieved. For the case of even n, the solutions at the ungrouped points near the boundaries are computed using (5).
(5) The steps below are performed directly once for points ✷ and △ in Figure 3 at the time level k + 1: (a) For type ✷ points, the rotated h − spaced five-point approximation formula derived by rotating the x − y axis clockwise at 45 o was used (Tan, Ali, & Lai, 2012). This approximation formula was applied to the right side of Equation (1), in combination with (3) being applied to the left side of (1), to obtain the following rotated C-N formula: (b) For type △ points, the typical h − spaced formula (5) is used. Note that formula (5) involve points of type △ only.

Stability analysis
A scheme is considered to be stable if the errors cease to increase with the passing of time, and gradually become inconsequential as the computation progresses. Even though the spacing is different, Equation (5) gives rise to both the FEG and FMEG methods. Therefore, the stability of both methods can be analyzed in similar ways. Here, we will show the stability of the FMEG scheme using eigenvalues of the generated matrices with mathematical induction.
The following lemma is important to prove the stability.

Computational complexity
This section presents an analysis of the computational complexity with regard to the three techniques described for solving (1), which is the number of arithmetic operations per iteration. For simplicity, we assume s x = s y , c x = c y . Suppose m 2 internal points exist within the solution, with m = n − 1, where n has an even mesh size, then the ungrouped points will be found close to the right/upper boundaries. Internal mesh points have two main types, namely, the iteration points that participate in the iteration process and the direct points that are directly calculated once using the rotated and standard difference formulas following the attainment of the iteration convergence. Table 1 lists the number of internal mesh points for the three earlier methods, whereas Table 2 provides a summary of the number of arithmetic operations that are needed for each iteration and the direct solution following the convergence, not only for the explicit group methods, but also for the fractional standard point (FSP) method.

Numerical experiments and discussion of results
Several numerical examples are presented in this section to prove the effectiveness of the fractional explicit group methods in solving the 2-D TFADE (1) with a Dirichlet boundary condition. A computer with Windows 7 Professional and Mathematica software having a Core i7 GHZ and 4 GB of RAM was used for conducting the experiments.