Distinctive Shape Functions of Fractional Differential Quadrature for Solving Two-Dimensional Space Fractional Diffusion Problems

: The aim of this study is to utilize a differential quadrature method with various kernels, such as Lagrange interpolation and discrete singular convolution, to tackle problems related to the Riesz fractional diffusion equation and the Riesz fractional advection–dispersion equation. The governing equation for convection and diffusion depends on both spatial and transient factors. By using the block marching technique, we transform these equations into an algebraic system using differential quadrature methods and the Caputo-type fractional operator. Next, we develop a MATLAB program that generates code capable of solving the fractional convection–diffusion equation in (1+2) dimensions for each shape function. Our goal is to ensure that our methods are reliable, accurate, efficient, and capable of convergence. To achieve this, we conduct two experiments, comparing the numerical and graphical results with both analytical and numerical solutions. Additionally, we evaluate the accuracy of our findings using the L ∞ error. Our tests show that the differential quad-rature method, which relies mainly on the discrete singular convolution shape function, is a highly effective numerical approach for fractional convective diffusion problems. It offers superior accuracy, faster convergence, and greater reliability than other techniques. Furthermore, we study the impact of fractional order derivatives, velocity, and positive diffusion parameters on the results


Introduction
Many researchers have been drawn to the study and use of fractional differential equations in recent decades. Classical differential equations are generalized by fractional order differential equations to explain complex problems such as viscoelasticity, biological engineering, fluid mechanics, engineering, stream through porous media [1,2], fractional operators [3], mathematical physics, and other fields. These equations are classified into three kinds based on the presence of fractional order derivatives in space, time, or space-time [4]. Currently, several various techniques are applied to solve fractional diffusion problems such as the Laplace, Mellin, and Fourier transforms. Mirza et al. [5] solved the time-fractionalized advection-diffusion equation by exploiting the Atangana-Baleanu fractional derivative operator. Khan et al. [6] performed a new modification of the Adomian decomposition method to solve the fractional convection-diffusion equation. Attar et al. [7] proposed an analytical method called Akbari-Ganji's technique for solving nonlinear fractional differential equations. Shah et al. [8] presented fractional-order diffusion equations via the Natural transform decomposition method.
Solving the fractional convection-diffusion equation numerically is a difficult task, especially for high-dimensional cases, due to the fractional derivative characteristic of the space fractional derivative differential operator. Various numerical methods have been proposed for approximating one-dimensional fractional convection-diffusion equations, including the finite difference method [9], Galerkin method [10], collocation method [11], homotopy analysis transform method [12], and finite volume element method [13]. Ding et al. [14] suggested an improved matrix transform numerical approach to solve the onedimensional space fractional advection-dispersion model, and an analytical solution was obtained using the Padé approximation. Recently, a shifted Grünwald-Letnikov difference operator for space and a Crank-Nicolson scheme for time were used to solve the space fractional convection-diffusion equation with variable coefficients, which achieved second-order convergence in both time and space with extrapolation [15]. Liu et al. [16] developed a radial basis functions finite difference (RBF-FD) for solving the time fractional convection-diffusion equation. Saadeh [17] implemented a space finite volume and finite difference approaches to solve the fractional convection-diffusion equation.
To handle two-dimensional space fractional diffusion problems, numerical approaches such as the Galerkin finite element [18], the alternating direction implicit (ADI) [19], the Kronecker product splitting [20], and the finite volume approaches [21] are utilized. Zhuang et al. [22] investigated the implicit difference approach of a time fractional diffusion equation. Pingyang et al. [23] investigated the explicit and implicit difference methods for time-space fractional reaction-diffusion equations, as well as their stability and convergence. Povstenko [24] and Zhang et al. [25] used the integral transform approach to achieve analytical findings for the time fractional diffusion-wave equation with a source term in cylindrical dimensions. The Jacobi collocation technique was employed to solve a particular case of the fractional advection-diffusion equation with a nonlinear source term by Parvizi et al. [11]. Bu et al. [26] introduced a finite element multi-grid technique for multi-term time fractional advection diffusion equations. Povstenko et al. [27] reported two approximation conclusions in two distinct modifications of the space-timefractional advection-diffusion equation. The Laplace and Fourier transforms were used to get the basic solutions, as well as the numerical results. Mohyud-Din Email et al. [28] discovered a fully implicit finite difference approach for the time fractional advection diffusion equation based on extended cubic B-splines. In Marin et al. [29], they dealt with the mixed initial boundary value problem for a dipolar structure in the context of thermo elastic theory as well as Hölder-type stability. It is commonly understood that fractional operators are nonlocal, and hence, independent of the numerical methods utilized. To achieve an acceptable computing cost, the Toeplitz or Toeplitz times-diagonal structure was used and explored [30]. Tuan et al. [31] utilized a finite difference discretization with Caputo derivative sense to solve a two-dimensional space fractional diffusion equation. Devshali and Arora [32] introduced the differential transform method and the differential quadrature technique to solve the fractional diffusion equation in two dimensions. Jannelli and Speciale [33] applied Lie symmetries and an implicit classical numerical method with the Caputo definition of fractional derivative to solve a two-dimensional time-fractional-diffusion-reaction equation. Zureigat et al. [34] offered two compact finite difference schemes with a fuzzy Caputo generalized Hukuhara derivative to obtain a certain solution for fractional convection-diffusion equation. Mustafa et al. [35,36] solved fractional-order chaotic systems with the Sinc shape function and analyzed drug diffusion via a thin membrane using the differential quadrature approach.
In our framework, we solve the fractional order differential equation in both time and space via a differential quadrature approach with distinct shape function and Caputo definition fractional derivative. The distinct shape functions are Lagrange interpolation function [37,38], Delta Lagrange kernel, and Regularized Shannon kernel [39], which have been effectively used in the (1+2) dimensional fractional convection-diffusion equation when paired with the block marching approach. Additionally, we design a MATLAB code for each approach to obtain a numerical solution for three problems to be examined. When compared to earlier analytical [40,41] and numerical [13,42,43] methodologies, the derived numerical results attain excellent efficiency and accuracy. Furthermore, we present some parametric studies to highlight the reliability of our methods with the influence of fractional order derivative, the velocity, and positive diffusion parameters on the results.
The organization of this work is as follows: The formulation of the problems used to explain the proposed schemes is detailed in Section 2. A detailed explanation of these methods can be found in Section 3. Section 4 presents the overall performance of the methods via three types of fractional partial differential problems. The findings obtained are also discussed in the same section before Section 5 summarizes the main conclusions.

Formulation of the Problem
To explain our idea, we consider three types of fractional partial differential equations on a finite domain that model advection and diffusion processes:

The One-Dimensional Riesz Space Fractional Advection Equation
In this section, we consider the following Riesz fractional advection equation: where 1 > > 2 and refers to a Riesz space fractional derivatives [44]. ( , ) is the solute concentration [45], " " represents the x-direction, and " " represents time. The expression "C x ≥ 0" represents the coefficient for the velocity parameter. Subject to the initial and boundary conditions are taken as follows [42]: Additionally, the exact solution is given in reference [42]:

The One-Dimensional Riesz Space Fractional Advection-Dispersion Equation
In this section, we consider the following one-dimensional Riesz space fractional advection-dispersion: where ≥ 0 and > 0 express the velocity parameter and positive diffusion coefficients, respectively; 0 < < 1 and 1 < < 2 are the Riesz space fractional order derivatives, defined a follows: ( ) are defined as the left-and right-side Riemann-Liouville derivatives.
The one-dimensional Riesz space fractional advection-dispersion Equation (5) has a physical meaning [46]. Physical considerations of a fractional advection-dispersion transport model restrict 0 < < 1, 1 < < 2, and we assume > 0 and ⩾ 0 so that the flow is from left to right. The physical meaning of applying homogeneous Dirichlet boundary conditions is that the boundary is placed sufficiently far away from an evolving plume that no significant concentrations reach it [47]. When = 1 and = 2, Equation (5) simplifies to the classical advection-dispersion equation. However, this paper exclusively focuses on fractional cases. Specifically, when = 0, Equation (5) transforms into the Riesz fractional advection equation as Equation (1) [48]; on the other hand, when ≠ 0, we obtain the Riesz fractional advection-dispersion equation [46] as Equation (5).
Subject to the initial and boundary conditions, they are given as follows [42]: In addition, the exact solution is offered in reference [42]:

The 2-D Fractional Diffusion Equation with Source Term [42]
Let us take a bounded domain as Ω = [0, ] × [0, ], Ω = [0, ] for our discretization of the problem. The aim is to find the numerical approximation of the 2-D fractional diffusion problem with zero Dirichlet boundary conditions over the finite domain Ω × Ω .
where 1 < , < 2 are fractional order derivatives, and and express the diffusion coefficients. With the source term: Subject to the initial and boundary conditions, they are given as follows [42]:

Method of Solution
Here, we seek to solve the fractional order differential equation in both time and space using a differential quadrature approach with distinct shape functions. These shape functions with block marching are used to solve (1+2) dimensions fractional convectiondiffusion problems.
Firstly, we start with the definition of a fractional derivative, where there exist many definitions, and the most accepted Caputo's definition is considered in this work.

Caputo's Fractional Derivative
Caputo proposed Caputo's fractional derivative, a novel definition of a fractional derivative based on the Riemann-Liouville Fractional Derivative [49]. Also, Weilbeer [50] showed this definition as follows: Let ∈ + , and if " " is a positive integer, then − 1 < < . The Riemann-Liouville Fractional Derivative of a function ( ) with order , expressing ( , ), is defined as follows: Caputo's Fractional Derivative of order is defined as follows: where ( ) is the fractional derivative of ( ), and " " is the integration lower limit. For = , the classical definition of integer order derivative is obtained. Secondly, we introduce the definition of differential quadrature method with the following different shape functions:

Lagrange Interpolation Polynomial (PDQM)
According to this shape function, the functional values for any unknown " " at a definite number of grid points " " can be expressed as follows [37,38]: Thus, the different derivatives of this unknown "υ" can be determined as follows: where ℜ ( ) is the weighting coefficient for the nth derivative. But the key to the accuracy of DQM is in determining the weighting coefficients. Therefore, they are different based on the shape function. Consequently, the weighting coefficients ℜ (1) for the first derivative and ℜ (2) for the second derivative can be found by differentiating Equation (17):

Discrete Singular Convolution (DSCDQM)
The singular convolution can be written as follows [51,52]: where ( − ) and ( ) are a singular kernel and a test function space element, respectively.
The shape function in this type depends on the choice of kernel type. But this shape function has numerous kernels, so we will use two of them to describe the functional values of the unknown " " and its derivatives at a certain number of grid points " ", as follows: Kernel (1)-the shape function of Delta Lagrange Kernel (DSCDQM-DLK) is described as follows: Consequently, ℜ (1) and ℜ (2) can be given by differentiating Equation (18) as follows [39,[52][53][54]: Also, in kernel (2), the shape function of the Regularized Shannon kernel (DSCDQM-RSK) is introduced as follows: where , , and Δ represent the Regularized Shannon factor, the computational parameter, and the step size, respectively.
Consequently, ℜ (1) and ℜ (2) can be determined by differentiating Equation (21) , This debate has revealed that the kernel type, regular grid points " ", and bandwidth (2 + 1) play an important role in obtaining convergence and accuracy solutions.
Equations (29) and (30) can be proved as follows: For ∈ (0,1) let Cconsequently, Further, Similarly, for , ∈ (1,2). Finally, to deal with time-dependent equations in the three problems and transform the governing equations to algebraic equations, we use the block marching method. This method provides more accuracy for DQM with three kernels. Therefore, we explain this technique as follows:

Block Marching Technique with Differential Quadrature Discretization
The governing Equations (1), (5), and (11) for the three problems are (1+2)-dimensional (x,y) and time-dependent (t) equations. The time-dependent models are solved using the block marching approach [56]. This approach divides the semi-infinite domain into numerous time intervals ( 1 , 2 , 3 , … ) in the t-direction. Each block consists of one interval , with the x-direction domain ranging from 0 to Lx and y-direction domain ranging from 0 to Ly.
The grid point distribution is equal for all blocks, and these weighting coefficients are applied to all blocks by producing 1 = 2 = 3 = ⋯. The following mesh sizes are utilized in the and directions in the nth block [57]: where ℝ is the number of blocks, N and M denote the x and y-mesh width, and L denotes the block's time level.

Numerical Results
In this section of our paper, the developed DQ techniques are examined on three test problems and demonstrate the accuracy and efficiency of the presented approaches. The developed techniques such as PDQM [37,38], DSCDQM-DLK, DSCDQM-RSK [51,52,58], and block marching with notion of Caputo are used for studying fractional order differential equations in both time and space for (1+2) dimensions fractional convection-diffusion equations. We carried out our calculations by creating MATLAB code for each method. The primary objective of our paper is to evaluate the performance, validity, efficiency, and accuracy of the developed techniques. We confirm this by comparing the computed results with existing numerical solutions [13,42,43,59] and analytical solutions [40,41,44]. To assess the convergence and accuracy of the developed methods, we use the error computation method outlined in [13,42,43,59] (27) and (28) in (1): Also, to deal with the boundary and initial conditions (3,4), this was conducted through adapting the governing equations. Now, we start to demonstrate the obtained results to know the reliability, stability, convergence, and performance of DQM based on three types of kernels and block marching with Caputo sense as follows: Table 1 shows the effect of applying uniform and non-uniform PDQM on computation of solute concentration ( , ) at = 1.85, = 2, = 8. It is noticed that the ∞ error is decreasing when the grid points are increasing, that is, the PDQM is stability in the x-direction, as well as the computed results via non-uniform grids are higher agree with earlier numerical solutions [13,42,43,59] than uniform ones. Also, ∞ error (1.0885 × 10 −6 ) and performance time (0.17 s) for non-uniform PDQM achieved the least. Therefore, the computed results' validity, efficiency, and accuracy to earlier numerical solutions via non-uniform PDQM are better than uniform PDQM. Increasing number of grids (30 × 30) leads to more accuracy at any number of blocks (ℝ), but 3-blocks with 8-levels give the fewest CPU time. Furthermore, non-uniform PDQM with × = 30 × 30, ℝ = 3, = 8 is better than uniform PDQM.  Tables 2 and 3 investigate the effect of values on the accuracy of DSCDQM-RSK and DSCDQM-DLK such as regularized Shannon factor ( ) , computational parameter ( ) , and step size (∆). Table 2 explains that DSCDQM-RSK is more accurate than DSCDQM-DLK for computing the solute concentration ( , ) , compared with earlier numerical [13,42,43,59] and exact [40,41,44] solutions. Also, the bandwidth (2 + 1 = 7) and ( = 1.45 × ∆) are the most suitable choices for numerical results for problem (4.1), which achieve more efficient results. In addition, Tables 2 and 3 exhibit that increasing the time levels (L) and decreasing the step size ( ) for each block ( = 0.1, = 12) lead to more efficiency at a lower bandwidth.    Table 4 presents a comparison between DSCDQM-DLK and DSCDQM-RSK at different grid points (N), times (T), and number of blocks (ℝ). The computed results proved that DSCDQM-RSK is the most efficient method compared with DSCDQM-DLK and PDQM. This preference belongs to low CPU time (0.12 s), ∞ error (7.0012 × 10 −6 ). Using the previously obtained conditions for = 12, ℝ = 3 , = 0.1, = 1.45 × ∆ at = 1.9, the numerical results of solute concentration ( , ) and ∞ error norms via DSCDQM-DLK and DSCDQM-RSK at different values of velocity parameter coefficient ( ) and time (T) are illustrated in Table 5. The results realized that increasing the values of the velocity parameter coefficient and time lead to more matching with exact values [40,41,44], with high accuracy with ∞ = 10 −7 . Table 5. Numerical results of solute concentration ( , )and ∞ error norms via DSCDQM-DLK and DSCDQM-RSK for problem (4.1) at different grid points (N) and times (T). = 1.9, = 12, ℝ = 3 , = 0. Now, we can use DSCDQM-RSK with the previous best conditions to introduce detailed parametric studies for problem (4.1). Figure 1 shows the variance of solute concentration with different fractions ∈ (1,2) and velocity parameter coefficients ( ) . Also, this figure demonstrates that when the value of ( ) increases, the solute concentration decreases, but it increases with increasing fraction power up to ( ≈ 1.75), then decreases.      (25) and (26) in (5): This equation is subjected to the boundary and initial conditions (8,9). Therefore, we deal with these conditions by adapting the governing equations.
In problem (4.2), we have solved this problem via DQM based on three types of kernels and block marching with Caputo sense. Therefore, we focus on determining the values whose influence is on the accuracy and convergence of the calculated results as uniform and non-uniform grid points, number of blocks, regularized Shannon factor ( ) , computational parameter ( ), and step size (∆). The values of these parameters are recorded in Tables 6-8. In Table 6, it is found that the computed results for solute concentration ( , ) at = 1, = 1.25, = 2, = 2, = 8 by PDQM with non-uniform grids are higher in agreement with earlier numerical solutions [13,42,43,59] than uniform ones at × = 30 × 30, ℝ = 3, = 8, ∞ error (3.7374 × 10 −6 ) and performance time (0.18 s).   Also, Tables 7 and 8 produce the best values for bandwidth (2 + 1 = 9), and = 1.3 × ∆ that achieves DSCDQM-RSK is more accurate than DSCDQM-DLK for computing the solute concentration ( , ) , compared with earlier numerical [13,42,43,59] and exact [40,41,44] solutions. And these tables explain that increasing the time levels (L) and decreasing the step size ( ) for each block ( = 0.1, = 12) lead to more efficiency at a lower bandwidth.    Exact [40,41,44] 0.54316 In addition, Table 9 demonstrates that DSCDQM-RSK is the most efficient method compared with DSCDQM-DLK and PDQM, which recorded low CPU time (0.13 s), ∞ error (3.0054× 10 -5 ).   Moreover, from all results in Tables 6-9, it is noted that the computed results are more accurate and converge at ℝ = 3 , = 0.1, = 12, = 1.3 × ∆ and = 0.5, = 1.9 via DSCDQM-RSK. Then, in Table 10, we present the numerical results of solute concentration ( , ) and ∞ error norms via DSCDQM-DLK and DSCDQM-RSK at the different values of the velocity parameter coefficient ( , ) and time (T). The results proved that increasing the values of ( , ) and time lead to more matching with exact values [40,41,44], with high accuracy with ∞ = 10 −7  Finally, we introduce the parametric studies for problem (4.2) via DSCDQM-RSK to show the flexibility and reliability of this method. Figure 4 explains the effects of fractions ∈ (1,2) and ∈ (0,1) on the solute concentration. This figure indicates that when the fraction (α) increased, the solute concentration decreases. But it increases with increasing the fraction power (β) up to x ≈ 1.75, then decreases. Also, Figure 5 exhibits the influence of the solute concentration with change in the time and fractions of ∈ (1,2) and ∈ (0,1). The solute concentration is directly proportional to ( , ) up to ≈ 1.75 , then decreases and is inversely proportional to time. And Figure 6 displays the relation between solute concentration and different and at different fractions of ∈ (1,2) and ∈ (0,1) . The solute concentration is inversely proportional to different coefficients ( and ). (x,t) x = =   (27) and (28) in (11): , , , , , , This equation is subjected to the boundary and initial conditions (13,14). Therefore, we deal with these conditions by adapting the governing equations.
In problem (4. Order of Convergence log where ‖ (∆ , ∆ , ∆ )‖ ∞ is the maximum error for a two-dimensional space fractional problem. Table 11 shows the stability and convergence of DSCDQM-RSK for solving a twodimensional space fractional problem. Also, the performance time of this method is calculated to ensure that DSCDQM-RSK is an effective method for such problems.  2) on the solute concentration ( , , ). These figures show that this method is flexible when altering the fraction order ( or ), which raises the solute concentration when they increase. According to Figure 8, when fractional orders are equal = = 2, the diffusion coefficients have a significant impact on the solute concentration results, and their effect increases as they increase.

Conclusions
This work presents effective numerical techniques for solving fractional advection and advection-dispersion equations in (1+2) dimensional Riesz space. The developed methods utilize unique shape functions, including Lagrange interpolation polynomial, delta Lagrange, and Regularized Shannon kernels, within the framework of the differential quadrature method (DQM). Additionally, the Caputo type is applied to treat the fractional differential operators. Also, to achieve more accurate and efficient results we used the block marching technique to deal with time dependency. Therefore, a MATLAB program was used to design code for solving three tested problems.
According to the presented techniques, we solved the considered problems in which the numerical evidence confirmed the earlier established numerical [13,[42][43][44]59,60] and exact [40,41,44] solutions. Further, we confirmed that the convergence of the offered techniques by calculating L∞ error norms. Therefore, the results computed via our methods are much more accurate, stable, and efficient. The best values of parameters, which control our methods, are = 12, ℝ = 3 , = 0.1, mesh size > (7 × 7), and = 1.45 × ∆ for problem (1) and = 1.3 × ∆ for problems (2)(3), and these values were achieved at ∞ = 10 −7 and computation time = 0.036 s. Furthermore, we investigated the effect of the velocity parameter, positive diffusion coefficients, and fractional order derivative at different locations and times on the solute concentration. The solute concentration is directly proportional to (α, β) up to (x ≈ 1.75), then decreases and is inversely proportional to time. Also, the solute concentration is inversely proportional to different coefficients (Cx and dx). If the fractional orders are the same, the diffusion coefficients will strongly influence the solute concentration results, and the impact will become greater as the coefficients increase. Additionally, it has been discovered that the suggested methods possess a distinct ability to solve fractional partial differential equations with both initial and boundary conditions. As a result, there is optimism that these techniques can be employed to tackle other nonlinear fractional problems that arise in a variety of fields within applied science. Funding: The researchers wish to extend their sincere gratitude to the Deanship of Scientific Research at the Islamic University of Madinah for the support provided to the Post-Publishing Program (2).

Data Availability Statement:
The data presented in this study are available in the article.