A novel local Hermite radial basis function‐based differential quadrature method for solving two‐dimensional variable‐order time fractional advection–diffusion equation with Neumann boundary conditions

A novel Hermite radial basis function‐based differential quadrature (H‐RBF‐DQ) method is presented in this paper based on 2D variable order time fractional advection–diffusion equations with Neumann boundary conditions. The proposed method is designed to treat accurately for derivative boundary conditions, which considerably improve the approximation results and extend the range of applicability for the method of RBF‐DQ. The advantage of the present method is that the Hermite interpolation coefficients are only dependent of the point distribution yielding a substantially better imposition of boundary conditions, even for time evolution. The proposed algorithm is thoroughly validated and is demonstrated to handle the fractional calculus problems with both Dirichlet and Neumann boundaries very well.


INTRODUCTION
Nowadays, the field of fractional calculus has attracted interest of researchers in several areas including mathematics, physics, chemistry, engineering, and even finance and social sciences. Main reason for this development is that a fractional model could describe simply and elegantly the complex characteristics of physical models. Therefore, the classical integer order models, for example, Maxwell model and Oldroyd-B model in non-Newtonian fluids, have been updated to the fractional order ones for complex viscoelastic materials [1][2][3]. Due to the growing number of applications of fractional calculus in science and engineering, numerical methods are being developed to provide tools for solving such problems [4]. Nevertheless, the intrinsic complexity of fractional calculus, caused partially by non-local properties of fractional derivatives and integrals, makes it rather difficult to find efficient numerical methods in this field. Despite this fact, however, the literature exhibits a growing interest and improving achievements in numerical methods for fractional calculus in general [5][6][7][8][9][10][11][12][13][14].
It seems enough to mention here that currently the numerical methods commonly used in fractional calculus included finite difference method and finite element method. The finite element method has now become the most widely used numerical method in research and application, and has played a great role in scientific research and engineering analysis. The finite element method discretized the complex region into elements, and a node based interpolation function approximation was established on the element for each variable. The governing equation was established by the variational principle or the weighted residual method. It effectively overcame the limitation of the finite difference method for the complex geometrical domain, and can handle the different boundaries flexibly. However, when the finite element method was used to deal with the problems with large deformation, dynamic crack propagation, and fluid solid interaction, the grids can be distorted. Therefore, mesh reconstruction for large deformation was unavoidable during the process of numerical simulation, which not only limited the computational efficiency, but also led to the damage of the numerical precision. In order to get rid of the solution dependence on the grid, the meshless method, only based on nodes, has made great progress since 1990s [15][16][17][18]. The meshless method with the structure of the test function is built on a series of discrete nodes. The connection between the field point and the node is no longer realized by the element. Hence, it gets rid of the constraints of the grid or element, and no longer needs the process of the grid reconfiguration when it involves the grid distortion and the grid movement. In order to simplify the mesh generation for the complex computational domains, the meshless/meshfree methods attracted more attention in recent years. With the effort of scientists, many kinds of effective meshless methods appeared, including smooth particle hydrodynamics method (SPH), element-free Galerkin method (EFG), meshless local Petrov-Galerkin method (MLPG), reproducing kernel particle method (RKPM), radial bases functions method (RBFs), finite point method (FPM), moving least square method (MLS), and so on [16][17][18]. Their main difference lies in the use of different test functions or different equivalent forms of differential equations. Generally, at the present, meshless methods can be divided into two categories, that is, weak-form method and collocation method [16]. In the former, the meshless method gives the discrete form of the problem based on the weak form controlled by variational principle of partial differential equations (PDEs). The characteristic of this method is that the solution is of high accuracy and good stability, but the amount of calculation is high, and the numerical integration is demanded. In the later, the collocation method directly satisfies the differential equation or boundary condition at discrete points. In addition, the collocation method is simple and does not require numerical integration. Furthermore, the computational efficiency is higher than that of the weak form method.
In 1990, Kansa first applied the collocation method and RBFs method with application to computational fluid dynamics [19,20]. After that, there are lots of research works on the applications of the RBFs collocation methods (see [11,[21][22][23][24][25][26][27][28][29][30][31] and their references). The collocation method with RBFs is a real meshless method for solving PDEs, and it is independent of the spatial dimensions [17,32]. As we all know, for stationary grid, the coefficients of finite difference scheme for spatial derivative is always fixed. But the common collocation method with RBFs for temporal problem, the interpolation coefficients should be calculated in every step during the time advancing process. To enhance the computational efficiency and treat complex geometries, Shu et al. developed a meshfree local radial basis function-based differential quadrature (RBF-DQ) approach to solve the two-dimensional incompressible Navier-Stokes equations [21]. RBF-DQ method originated from the concept of differential quadrature (DQ) [33,34], however, it cannot directly be applied to the problems with irregular geometries. Due to the vast flexibility, RBF-DQ method has been successfully used to study many scientific and engineering problems with irregular geometries [22,32,[35][36][37][38][39][40][41]. In this method, the interpolation coefficients are only dependent on the distributions of the points and independent of the solution. However, the accuracy achieved by direct meshless collocation scheme is a bit poor especially on the boundaries [42,43]. Furthermore, the collocation scheme has difficulties in dealing with Neumann boundary conditions. In order to improve the accuracy on the derivative boundaries, recently, Hermite RBFs method has been developed by many authors [16,17,[43][44][45]. It should be noted that, for the most of the Hermite RBFs methods, the interpolation coefficients depend on the solutions. In this paper, we developed a novel Hermite RBF-DQ (H-RBF-DQ) method for solving PDEs based on a variable-order time fractional advection-diffusion equation with both Dirichlet and Neumann boundary conditions. The advantage of the present method is that the interpolation coefficients are only dependent of the point distributions, therefore, it can greatly reduce the cost of computation. The present method can be seen as an extension of local RBF-DQ meshless approach for Neumann boundary problems in [46], but the accuracy for approximating to the Neumann boundary conditions has been much improved. Without loss of generality, in this paper, the method is validated by a special variable-order time fractional advection-diffusion equation with Neumann boundary conditions.
There are some researches on meshless method for the numerical simulations of time fractional advection-diffusion equations. For example, MLS method is developed to solve the constant-order time fractional advection-diffusion equation in the papers [47,48]. For the variable-order time fractional advection-diffusion equation, it can be referred to the paper [49]. However, the shape function in MLS method is not satisfied with Kronecker function character and the boundary conditions cannot be directly enforced, which increase the difficulty for the treatment of boundaries. Furthermore, the above works were all on the problem with Dirichlet boundary conditions. In [11], a global Kansa method is developed to solve a variable-order time-fractional advection-diffusion equation arising in anomalous transport. And also, in [50], a novel global RBF-based meshless method with particular boundary treatment for solving time-fractional transport equations is presented. In [46], a meshless local RBF-DQ approach to solve a variable-order time fractional advection-diffusion equation subjected to Neumann boundary conditions. In [51][52][53], boundary collocation and localized collocation methods for anomalous heat conduction and diffusion on surfaces are presented. In order to improve the approximating accuracy on the Neumann boundary conditions, in this paper, a novel H-RBF-DQ method is proposed to deal with both Dirichlet and Neumann boundary conditions based on a two-dimensional variable order time fractional advection diffusion equation. The numerical results demonstrate the robustness and the versatility of the proposed method.
The paper is organized as follows. In Section 2, the modeling equations and their time discretization are presented. The review of RBF-DQ method is presented in Section 3. For the treatment of the Neumann boundary condition, in Section 3, based on the Hermite RBF interpolation, a novel H-RBF-DQ method is proposed to improve the approximating accuracy on the Neumann boundaries.
Numerical results with applications are presented in Section 4 to demonstrate the accuracy of the proposed method. Finally, the conclusion is presented in Section 5.

THE MODELING EQUATION AND TIME DISCRETIZATION APPROXIMATION
The variable-order time fractional advection-diffusion equation can be formulated as [49] c (1) subject to the following initial and boundary conditions where Δu =  . Ω is a bounded domain in R 2 , Ω is the boundary of Ω, (x, t) > 0 is the diffusion coefficient function, v(x,t) is advection velocity, and f (x, t), g(x), h(x, t) are given functions.
In Equation (1), the variable-order time fractional where Γ is the Gamma function.
In this study, we focus on the Hermite RBF-DQ method for the variable-order time fractional advection-diffusion equation on complex geometries with Neumann boundary conditions. First, in this section, we give the time discretization approximation for Equation (1). Using the time discretization method proposed in [49] for Caputo derivative with the notation where t k = kΔt for k = 0, 1, … , M and M = T/Δt. The truncation error R k+1 is subjected to [47,49,54] Substituting Equation (4) into Equation (1), we can obtain where and where C and C are constant.
Using the notations u k = u k (x) as the numerical approximation to u( (1) can be discretized as follows Equation (9) can be rearranged as

SPATIAL DISCRETIZATION METHOD
In order to obtain a final numerical scheme of Equation (10), the spatial discretization of variable-order time fractional advection-diffusion equation should be given. In this section, we first review RBF-based differential quadrature method for the spatial discretization, which has been reported in our recent paper [46]. Then the details of the new Hermite RBF-based differential quadrature method will be proposed to improve the accuracy of approximating solutions near the Neumann boundary.

Review of RBF-based differential quadrature method
The idea of differential quadrature method came from the numerical integral, that any integral over a closed domain can be approximated by a linear weighted sum of all the functional values in the integral domain. In differential quadrature method, similar to other numerical integral, the derivative value u (m) (x) at the center x i are approximated by a linear weighted sum of the function values at a set of nodes in a closed domain as for i = 1, · · ·, N, where N is the total number of the points inside the computational domain, and w (m) ij , w (m) ij are the weighting coefficients for derivatives of order m with respect to x and y, respectively. The weighting coefficients of w (m) ij , and w (m) ij are determined by the base functions as the test function in Equation (11), which will greatly enhance the computational efficiency [21,55]. Although, there are many different kind of RBFs available, due to the excellent performance for the scattered data interpolation, the multiquadrics (MQ) basis function is used here as the test function [21,22]. In the RBF-DQ method, the function in the region of Ω can be locally approximated by MQ RBFs as with shape parameter c j . To make the problem be well-posed, the equation is enforced. Substituting Equation (13) into Equation (12) gives where The number of unknowns in Equation (14) is N. As the setting in the paper [21], N+1 can be replaced by i , and Equation (14) can be rewritten as where g i (x, y) = 1 and g j (x, y), j = 1, · · ·, N but j ≠ i given by Equation (15), formed a base vector to cover the function space of h(x, y).
In RBF-DQ method, the most important contribution is that the weighting coefficients of w (m) ij , and w (m) ij are determined by all the base functions as the test function in Equation (11). Based on this idea, we can obtain The N equations of (17) formed a linear system with N unknowns for the given i, and the weighting coefficients w (m) ik can be solved by a numerical method. In a similar manner, the weighting coefficients w (m) ij of the y-derivatives can also be computed by Equation (17) with x substituted by y. Substitution of Equation (11) into time discretized Equation (10) at point In Equation (18), the following notations are used It should be noted that when the number of nodes, N, is large, the coefficient matrix of (17) may be ill-conditioned. This limits its application. Hence, we use the local RBF-DQ method to solve the variable-order time fractional advection-diffusion equation.
The key of local RBF-based differential quadrature method is that the m-order derivatives of u(x) at x i are approximated by the function values at a set of nodes in the neighborhood of x i with N i supporting nodes (including x i ). That is The corresponding coefficients w (m) ij and w (m) ij can be determined by Equation (17) with N i local support nodes in the neighbor of x i .
Substitution of Equation (19) into time discretized Equation (10) at point x i , by the local RBF-based differential quadrature method, yields for i = 1, 2, · · ·, N. Although the discretization of Equation (1) is simple, the RBF-DQ approximation of the function contains a shape parameter c that could be node-dependent and must be determined by the user. In this study, we use the method of normalization of supporting region suggested in the paper of [21], which can also be referred to [46].

Hermite RBF-DQ method
In order to improve the approximated accuracy near the boundary with Neumann boundary condition, for the unknown computational nodes with supporting points on the Neumann boundary, we suggest to use the Hermite method and include the information on the boundaries to obtain the spatial discretization of Equation (10). For convenience of illustration, a sample of mesh points distribution in H-RBF-DQ method on hybrid boundary conditions is shown in Figure 1, where the Neumann boundary conditions are enforced on the upper boundary and the Dirichlet boundary conditions are satisfied on the lower boundary. According to the position of the unknown points, we divide them into four categories. Type A denotes the point on the Neumann boundary conditions. Type B denotes the internal point with supporting points on the Dirichlet boundary conditions. Type C denotes the point with all internal supporting points. And Type D denotes the internal point with supporting points on the Neumann boundary conditions (see Figure 1). For the relationship structure of points, we use the algorithm in paper [56]. The approximation of a function u(x) (see [16,43]) can be written in a Hermite interpolation form by a linear combination of RBFs at all the N i supporting nodes (including the N b i Neumann boundary points) within the local support domain and the normal derivatives along the normal direction n = (n x , n y ) at the Neumann boundary points as are interpolation coefficients, and j , b l are radial basis function, respectively.
In this study, similar to the function g j (x, y) given in the previous section, the MQ basis function is used as the radial basis function, that is is the coordinate for the lth normal derivative boundary point. So we can get b l With the constrains of (13) and the notation of (15), Equation (21) can be reformulated as and where g i (x, y) = 1, g j (x, y), j = 1, · · ·, N i but j ≠ i and b l n l , l = 1, · · · , N b i are a base vector for the function space of h(x, y).
and the coefficients vectorâ can be obtained by the interpolations at N i points and the normal derivatives at N b i points on the Neumann boundary in Equations (24) and (25). The process can also be expressed by matrix formulation as followsû = Ψâ.
Thus, the unknown coefficients vectorâ Finally, the approximation form of function can be obtained like and the matrix of shape functions can be expressed as follows Hence, the function u(x) can be approximately expressed as From the interpolation formulation (33), we can directly get Then the mth order derivative value u (m) x m at the center x i can be approximated by a linear weighted sum of the function values at N i points and the normal derivatives at N b points on Neumann boundary as In order to get the coefficients in Equations (35), by using the idea of differential quadrature in Section 3, we take all the base functions in (31) as the test functions and solve the following system In a similar manner, the weighting coefficients w (m) ik , and s (m) ij of the mth order y-derivatives can also be computed by Equation (36) with x substituted by y.
With the notation of (35) for the spatial derivative, the numerical scheme at the node D (see Figure 1) with supporting points on the Neumann boundary will become where the source term F k+1 i denotes the impact on Neumann boundary conditions, defined by The boundary treatment for Dirichlet boundary condition in collocation method is trivial. However, the collocation scheme has difficulties in dealing with Neumann boundary conditions. For the present Hermite-type RBF-DQ method (37), the Neumann boundary condition at point ( (see the point A in Figure 1) can be formulated as Although, the derivatives in the Neumann condition can also be discretized by the local RBF-DQ method of (19) as the accuracy is much lower than the present H-RBF-DQ method, which will be shown in the next section.
It is worth pointing out that, in the present H-RBF-DQ method, the normal derivatives at the Neumann boundary nodes are added as additional DOFs, which are enforced by Equation (35). For an unknown collocation node at x i , if its local support nodes do not include the Neumann boundary node, the conventional differential quadrature formulas (19) are used to approximate the derivatives in Equation (10), and Equation (20) is used to discretize the variable-order time fractional advection-diffusion equation. If its support neighboring nodes include the Neumann boundary nodes, the formulation (35) is used. Hence, by the present H-RBF-DQ method, if the computational node with supporting points on the Neumann boundary, the scheme of (37) is used, otherwise the scheme (20) is used. The throughout procedure can be referred to the Algorithm 1. Although the above description is based on the variable-order time fractional advection-diffusion equation, the same treatment is still suitable for other Neumann boundary-value problems.
We must point out that the present Hermite RBF point interpolation method is different from ones in [16,43]. The method in [16,43] firstly required to obtain the matrix of shape functions in (30); then the differential of components of must be calculated to obtain Equation (34). Furthermore, for time related problems, their method demands that the calculation of the matrix of shape functions with their differential of components should be solved at every time step. In our method, however, the coefficients in (35) is only related to the support points set; and for time related problems, the system of (36) is just solved only once, therefore, the coefficients in (35) can be stored before the time evolution. (1) and (2) 2: generate the mesh pointes 3: set the time step Δt 4: set the label to distinguish the boundary nodes and inner nodes 5: ascertain the supporting points for every unknown node 6: calculate the weights in Equations (19) and (35)

1: input the known terms in Equations
with the exact solution u(x, t) = x 2 + y 2 + t 2 . Equation (41) was given in paper [49] to evaluate the MLS method to deal with Dirichlet boundary condition. In the present study, the Neumann boundary condition is considered by the present Hermite RBF-DQ method. The unit circle domain is considered to test the accuracy. The two constant orders ( = 0.5, 0.8) and a variable order are considered. The spatial convergence is evaluated by three different grids. For example, the point distribution with N = 214 is shown in Figure 2a. A contour plot of numerical solution is shown in Figure 2b. The comparison between the new H-RBF-DQ method (39) and the RBF-DQ method (40) (see the paper [46]) is presented in Table 1. The results show that the errors can be greatly reduced by the use of the present method. This demonstrates that the new developed H-RBF-DQ method is more accurate than the RBF-DQ method. In addition, the matrix condition numbers and time costs are also provided in Table 1. Although the condition numbers for the H-RBF-DQ method are slightly bigger than the RBF-DQ method, they are still in a reasonable region. In Table 2, the numerical comparisons at final time t = 2 are also given. It is very interesting that the L ∞ and RMS error accumulation are slowly increased. However, for the present special example, with the increase of time, the solution will increase, which will lead to the decrease of relative L 2 error. The grid convergence for variable-order (42) is provided in Table 3. The condition numbers and the time costs are also reported in Table 3. The time convergence for variable-order (42) is provided in Table 4.
Similar to the RBF-DQ method, the approximated accuracy for the H-RBF-DQ method is also affected by the shape parameter c. The impact of the shape parameter c is drawn in Figure 3, which shows the choice c 2 = 5 is suitable for the present example. For the choosing of c 2 , one can refer to papers [21,46,[57][58][59].

Example 2
We will consider a variable-order fractional equations, say     Abbreviation: RMS, root-mean-square. In order to show the versatility of the method, an irregular computational domain with Neumann boundary condition is considered again. The computational domain is formed by the polar coordinate method [60], and is defined by = n + 1 n 2 (n + 1 − cos n ), ∈ [0, 2 ]. As an example, the geometry for n = 6 with = 0.1 is considered. The distribution of grid points with total number 2125 points is shown in Figure 4. The numerical solution and the error are shown in Figure 5, which demonstrate the present method is effective for Gauss pulse. Furthermore, a complicated boundary condition, a hybrid Dirichlet and Neumann boundary condition is also studied here. For the sake of simplification, we rethink Equation (43) with hybrid boundary conditions, but the computational domain is formed by a pentagonal star and a circular arc with 4 which is shown in Figure 6. It is worth pointing out that such a geometry was first proposed in [46] for testing the RBF-DQ method. The point coordinates on the boundary of the pentagonal star are represented by the following formula. where the parameters R = 2 and r = 1 are considered. As a hybrid boundary condition, the Neumann boundary condition is enforced on the right of the shape labeled by red "⋆," and the Dirichlet boundary condition is enforced on the rest. Clearly, the computational domain is singular, which increases the difficulty of using numerical methods. The meshless distribution with 1364 points is shown in Figure 6. The total cost of computation is 140 s. It can be seen that although the boundaries with the Neumann boundary conditions include singular corner points, the final numerical result seems still very fine as shown in Figure 7, where is the error distribution plot. Hence, it can be concluded that the present method handles the Neumann boundary conditions very well.

Example 3
As the final example, we studied a fractional approximate modeling of air pollution, in which the process of transport and diffusion of pollutants in the atmosphere with obstacles is investigated. The equation represents a mass balance statement of pollutant materials when they are transported through the air, and is defined by where C(x, t) is the concentration of pollutant materials, v is the wind velocity and is the diffusion coefficient.
Although, the authors in [49] considered the modeling of air pollution, the homogeneous Dirichlet boundary condition is not realistic in the practical applications and the computational domain without  The numerical solution and the error of Example 2 enforced by the hybrid Dirichlet and Neumann boundary conditions obstacles is very simple. Because the penetrability cannot be satisfied near the computational boundaries, the method suggested in [49] is assumed that the air pollution does not touch the boundaries. In our recent research [46], a solution for this problem was provided from an application with the Neumann boundary conditions with obstacles by using RBF-DQ method, in which the geometry is complicated with obstacles. For this case, the wind velocity distribution affected by the obstacles is not constant, and is obtained by the method of computational fluid dynamics.
Here, we reconsidered the air pollution modeling with constant = 0.2, but with (x, t) = 0.55 and (x, t) = 0.55 + 0.45 sin(xyt), and a single source of pollutant with strength f = 5 is located in the domain of 2.5 ≤ x ≤ 2.7 and 2.5 ≤ y ≤ 2.7. In this application, a computational grid with 3145 points and the time step Δt = 0.02 are used. For the simplicity, the boundaries of holes in the Using the H-RBF-DQ method, the pollutant distributions at different time are shown in Figure 8. It is found that the present method captures the pollutant feature very well throughout the computational domain. As expected, the concentration of pollutant materials with the variable-order time model (x, t) = 0.55+0.45 sin(xyt) moves more slowly than that with the constant-order model (x, t) = 0.55. It should be noted that the present method guarantees the values of concentration are positive, and the air pollution is successfully transferred across the boundaries as shown in Figure 8. The results show that the proposed method can properly capture the distribution of pollutant on the domain with a given velocity field.

CONCLUSION
A local meshless H-RBF-DQ method on Neumann boundary condition is developed in this paper to simulate the fractional calculus problems with complex geometries. Although the method is presented based on the variable-order time fractional advection-diffusion equations, it is also suitable for other PDEs with the Neumann boundary conditions. It is worth noting that the present Hermite RBF method is very different from ones in [16,43]. The advantage for the use of the current method is that the coefficients of the shape functions are only related to the support points set, hence, for time related problems, the system of equations is just solved only once, therefore, the coefficients can be stored before the time evolution. Numerical examples show the present method is more accurate than the RBF-DQ method. The current method is successfully applied to a variable-order fractional model of air pollution. It is an evident that the numerical method can keep the values of the concentration are positive. Furthermore, our numerical examples demonstrate that the current method is suitable for complicated boundary value problems, including hybrid Dirichlet and Neumann boundary conditions. In the future work, the high order temporal discretization for other PDEs will be studied.