# A novel Hermite RBF-based differential quadrature method for solving two-dimensional variable-order time fractional advection-diffusion equation with Neumann boundary condition

###### Abstract

In this paper, a novel Hermite radial basis function-based differential quadrature method (H-RBF-DQ) is presented. This new method is designed to treat derivative boundary conditions accurately. The developed method is very different from the original Hermite RBF method. In order to illustrate the specific process of this method, although the method can be used to study most of partial differential equations, the numerical simulation of two-dimensional variable-order time fractional advection-diffusion equation is chosen as an example. For the general case of irregular geometry, the meshless local form of RBF-DQ was used and the multiquadric type of radial basis functions are selected for the computations. The method is validated by the documented test examples involving variable-order fractional modeling of air pollution. The numerical results demonstrate the robustness and the versatility of the proposed approach.

Jianming Liu and Xinkai Li^{1}^{1}1Corresponding author. xkl@dmu.ac.uk

School of Mathematics and Statistics, Jiangsu Normal University, Xuzhou 221116, China

Faculty of Technology, De Montfort University, Leicester LE1 9BH, England

Keywords: Hermite RBF-DQ; Differential quadrature method; Variable-order time fractional; Neumann boundary condition; Radial basis function

## 1 Introduction

Currently, numerical methods commonly used in scientific and engineering calculations include 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 discrete the complex region into elements, and a node based interpolation function approximation is established on the element for the solving variable. The control equation is established by the variational principle or the weighted residual method. It effectively overcomes the limitation of the finite difference method for the shape of the region, and can handle the boundary of various shapes flexibly. However, when the finite element method is used to deal with the problems with characters of large deformation, dynamic crack propagation, fluid solid interaction and so on, the grid can be distorted sometimes. Therefore, mesh reconstruction is unavoidable during the process of numerical simulation, which not only reduces the computational efficiency, but also leads to the damage of the calculation precision. In order to get rid of the dependence on the grid, the meshless method, only based on nodes, has made great progress since 1990s. The structure of the test function of the meshless method is built on a series of discrete nodes. The connection between the field point and the node is no longer realized by the element, so it gets rid of the constraints of the grid or element, and no longer needs the grid reconfiguration when it involves the grid distortion and the grid movement [3].

To treat the complex computational domains and avoid the mesh generation, meshless techniques have become very popular in recent years. There are many kinds of meshless methods at present, such as, 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, 7]. Their main difference is the use of different try functions or equivalent forms of differential equations. At present, meshless methods can be divided into two categories: weak-form method and collocation method. From the weak form variational principle of partial differential equation, the weak form meshless method gives the discrete form of the problem. 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. The collocation method directly satisfies the differential equation or boundary condition at discrete points. This method is simple and does not require numerical integration. The computation efficiency is higher than that of the weak form method.

In 1990, Kansa combined the collocation method with RBFs for computational fluid dynamics [11, 12]. After that, there are lots of research on the applications on the RBFs collocation methods in the papers[22, 21, 25, 19, 4, 6, 1] and their references. RBFs collocation method is a real meshless method for solving partial differential equations, and it has nothing to do with the spatial dimension [7, 5]. To treat any geometries, Shu et al. developed a meshfree local RBF-DQ approach to solve the two-dimensional incompressible Navier-Stokes equations [22]. Radial basis function-based differential quadrature method (RBF-DQ) is a meshless method which has originated from the concept of differential quadrature. Classical differential quadrature (DQ) method began from the idea of conventional integral quadrature [20], but it cannot directly be applied to problems with irregular geometries. Due to the vast flexibility, RBF-DQ method has been successfully used to study many scientific and engineering problems [21, 10, 2, 5, 9, 8]. However, the accuracy achieved by direct meshless collocation scheme is a bit poor especially on boundary [15, 17]. Furthermore, the collocation scheme has difficulties in dealing with Neumann boundary conditions. In order to improve the accuracy near the derivative boundary, recently, Hermte RBFs methods have been developed by many authors [17, 16, 7, 23, 13]. In this paper, we developed a novel Hermite RBF-based differential quadrature method for solving partial differential equations (PDEs). The present method can decrease greatly the cost of computation. The present method can be seen as an extension of meshless local RBF-DQ approach in [22], but the accuracy approximating the Neumann boundary conditions can be improved greatly. To show the method, we chose a special two-dimensional variable-order time fractional advection-diffusion equation with Neumann boundary condition to demonstrate.

There are some researches on meshless method for the numerical simulations of time fractional advection-diffusion equation. MLS method is presented to solve the constant-order time fractional advection-diffusion equation in the paper [27, 18]. For the variable-order time fractional advection-diffusion equation, we can refer to the paper [24]. However, the shape function in MLS method is not satisfied with Kronecker function character and the boundary condition can not be directly enforced, which increase the difficulty for the treatment of boundary. Furthermore, the above works are all on the problem with Dirichlet boundary conditions. In recent work, we develop a meshless local RBF-DQ approach to solve the variable-order time fractional advection-diffusion equation. In order to improve the approximating accuracy near the Neumann boundary, in this paper, we develop a novel Hermite RBF-DQ method.

The paper is organized as follows. In Section 2, the modeling equation 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 4, based on the Hermite RBF interpolation, a novel Hermite RBF-DQ method is described to improve the approximating accuracy. In Section 5, computational results are presented based on two-dimensional variable-order time fractional advection-diffusion equation.

## 2 The modeling equation and time discretization approximation

The variable-order time fractional advection-diffusion equation [24] can be formulated as

(1) |

subject to the following general initial and boundary conditions

(2) |

where , , and denotes the operator with Dirichlet or Neumann boundary conditions. is a bounded domain in , is the boundary of , is the diffusion coefficient function, v(x,t) is advection velocity, and , , are given functions.

In Equation (1), the variable-order time fractional derivative with in the Caputo definition as

(3) |

In this study, we focus on the RBF-DQ method for the variable-order time fractional advection-diffusion equation on complex geometries with different boundary conditions. First, in this section, we give the time discretization approximation. Use the time discretization method proposed in the paper [24] for Caputo definition, with the notation we have

(4) |

where for and . The truncation error [14, 24, 27] is subjected to

(5) |

Substituting Equation (4) into Equation (1), and with the -weighted scheme (), we can obtain

(6) |

where

(7) |

and

(8) |

## 3 RBF-based differential quadrature method

### 3.1 Basic 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, as the numerical integral, the derivative value at the centre are approximated by a linear weighted sum of the function values at a set of nodes in a closed domain as

(11) |

for where are the weighting coefficients for derivatives of order with respect to and , respectively. In RBF-DQ method, the weighting coefficients of and are determined by all the base functions as test function in Equation (11) [26, 22].

There are many RBFs available. In this study, due to the better performance for the interpolation of 2D scattered data, multiquadrics (MQ) basis function is used as the test function. The function in the region of can be locally approximated by MQ RBFs as

(12) |

with shape parameter . To make the problem be well-posed, the equation

(13) |

is enforced. Substituting Equation (13) into equation (12) gives

(14) |

where

(15) |

The number of unknowns in Equation (14) is . As the setting in the paper [22], can be replaced by , and Equation (14) can be written as

(16) |

Where and but given by Equation (15) are a base vector for the function space of .

In RBF-DQ method, the weighting coefficients of and are determined by all the base functions as the test function in Equation (11), and we can obtain

(17a) | |||

(17b) |

### 3.2 Local RBF-based differential quadrature method

When the number of knots, , is large, the coefficient matrix of (17) may be ill-conditioned. This limits its application. Hence, we mainly 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 -order derivatives of at are approximated by the function values at a set of nodes in the neighborhood of with nodes (including ). That is

(20) |

The corresponding coefficients and can be determined by Equation (17) with local support nodes in the neighbor of .

Substitution of Equation (20) into time discretized equation (10) at point , by the local RBF-based differential quadrature method, yields

(21) |

for .

As shown in the previous subsection, the RBF-DQ approximation of the function contains a shape parameter that could be knot-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 [22] and set .

## 4 Boundary treatment

The boundary treatment for Dirichlet boundary condition in collocation method is trivial. However, the collocation scheme has difficulties in dealing with Neumann boundary conditions. In this paper, for the variable-order time fractional advection-diffusion equation, we present a new Hermite-type RBF-DQ method for the Neumann boundary condition.

### 4.1 Hermite RBF interpolation

The approximation of a function (see the papers [17, 16]) can be written in a linear combination of RBFs at all the nodes (including the Neumann boundary points) within the local support domain and the normal derivatives along the normal direction at the Neumann boundary points as

(22) |

where are interpolation coefficients, and are radial basis function, respectively.

In this study, as the Section 3, the MQ basis function is used as the radial basis function, that is

(23) |

where is the coordinate for the th normal derivative boundary point. So we can get

(24) |

With the constrains of (13) and the notation of (15), the equation of (22) can be reformulated as

(25) |

and

(26) |

where , but and are a base vector for the function space of .

Let

(27) |

and the coefficients vector

(28) |

can be obtained by the interpolations at points and the normal derivatives at points on Neumann boundary in Equations (25) and (26). The process can also be expressed by matrix formulation as follows

(29) |

Thus, the unknown coefficients vector

(30) |

Finally, the approximation form of function can be obtained like

(31) |

where

(32) |

and the matrix of shape functions can be expressed as follows

(33) |

Hence, the function can be approximately expressed as

(34) |

### 4.2 Hermite RBF-DQ method

From the interpolation formulation (34), we can directily get

(35) |

Then the -th order derivative value at the centre can be approximated by a linear weighted sum of the function values at points and the normal derivatives at points on Neumann boundary as

(36) |

In order to get the coefficients in Equations (36), by the idea of differential quadrature in Section 3, we use all the base functions in (32) as the test functions and solve the following system

(37a) | |||

(37b) | |||

(37c) |

In a similar manner, the weighting coefficients of the -th order -derivatives can also be computed by Equation (37) with substituted by . Then the Neumann boundary condition at point can be formulated as

(38) |

We must mention that the present Hermite RBF point interpolation method is very different from the ones in [17, 16]. The method in [17, 16] firstly required to obtain the matrix of shape functions in (31); then the differential of components of must be calculated to obtain Equation (35). Furthermore, for time related problems, their method demand that the calculation of the matrix of shape functions with their differential of components should be solved at every time step. In our method, the coefficients in (36) is only related to the support points set; and for time related problems, we just solve the system of (37) once. So the coefficients in (36) can be stored before the time evolution.

In the present Hermite RBF-DQ method, the normal derivatives at the Neumann boundary nodes are added as additional DOFs and Equations (38) are enforced. For an internal collocation node at , if its local support nodes do not include Neumann boundary node, the conventional differential quadrature formulas (20) are used to approximate the derivatives in Equation (10), and the Equation (21) is used to discretize the variable-order time fractional advection-diffusion equation. If its support neighboring nodes include Neumann boundary nodes, the formulation (35) is used.

Although, the derivatives in the Neumann condition can also be discretized by the local RBF-DQ method of (20) as

(39) |

the accuracy is far lower than the present Hermite RBF-DQ method, which will be shown in the next section.

## 5 Numerical examples

In this section, numerical experiments are carried out to demonstrate the effectiveness of RBF-DQ method developed in this paper for variable-order time fractional advection-diffusion equation with Dirichlet or Neumann boundary condition. Three different test problems are chosen to show the capability and accuracy of the proposed method. These test problems come from the paper [24], but the boundary conditions with Neumann condition and different solution domains are considered. Although Equation (21) is valid for any value of , we use as the famous implicit scheme in all numerical examples.

To show the accuracy of the proposed method, the errors and Root-Mean-Square (RMS) of errors are measured using the following definitions:

where is the numerical solution of .

Example 1. Consider the following variable-order time fractional advection-diffusion equation

(40) |

## 6 Conclusion

## References

- [1] Victor Bayona, Natasha Flyer, Bengt Fornberg, and Gregory A. Barnett. On the role of polynomials in rbf-fd approximations: Ii. numerical solution of elliptic pdes. Journal of Computational Physics, 332:257 – 273, 2017.
- [2] Y.L. Chan, L.H. Shen, C.T. Wu, and D.L. Young. A novel upwind-based local radial basis function differential quadrature method for convection-dominated flows. Computers & Fluids, 89:157 – 166, 2014.
- [3] Yumin Cheng. Meshfree Methods. Science Press, Beijing, 2015.
- [4] Mehdi Dehghan, Mostafa Abbaszadeh, and Akbar Mohebbi. A meshless technique based on the local radial basis functions collocation method for solving parabolic-parabolic Patlak-Keller-Segel chemotaxis model. Engineering Analysis with Boundary Elements, 56:129 – 144, 2015.
- [5] Mehdi Dehghan and Vahid Mohammadi. The numerical solution of Cahn-Hilliard (CH) equation in one, two and three-dimensions via globally radial basis functions (GRBFs) and RBFs-differential quadrature (RBFs-DQ) methods. Engineering Analysis with Boundary Elements, 51:74 – 100, 2015.
- [6] Mehdi Dehghan and Vahid Mohammadi. A numerical scheme based on radial basis function finite difference (RBF-FD) technique for solving the high-dimensional nonlinear Schrödinger equations using an explicit time discretization: Rungeâkutta method. Computer Physics Communications, 217:23 – 34, 2017.
- [7] G. Fasshauer. Meshfree Approximation Methods with MATLAB. World Scientific Publishers, Singapore, 2007.
- [8] Ahmad Golbabai and Ehsan Mohebianfar. A new method for evaluating options based on multiquadric RBF-FD method. Applied Mathematics and Computation, 308:130 – 141, 2017.
- [9] Ahmad Golbabai and Ahmad Nikpour. Computing a numerical solution of two dimensional non-linear Schrödinger equation on complexly shaped domains by RBF based differential quadrature method. Journal of Computational Physics, 322:586 – 602, 2016.
- [10] M.R. Hashemi and F. Hatam. Unsteady seepage analysis using local radial basis function-based differential quadrature method. Applied Mathematical Modelling, 35(10):4934 – 4950, 2011.
- [11] E.J. Kansa. Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-I surface approximations and partial derivative estimates. Computers & Mathematics with Applications, 19(8):127 – 145, 1990.
- [12] E.J. Kansa. Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-II solutions to parabolic, hyperbolic and elliptic partial differential equations. Computers & Mathematics with Applications, 19(8):147 – 161, 1990.
- [13] Artur Krowiak. Hermite type radial basis function-based differential quadrature method for higher order equations. Applied Mathematical Modelling, 40(3):2421 – 2430, 2016.
- [14] Yumin Lin and Chuanju Xu. Finite difference/spectral approximations for the time-fractional diffusion equation. Journal of Computational Physics, 225(2):1533 – 1552, 2007.
- [15] T.J. Liszka, C.A.M. Duarte, and W.W. Tworzydlo. hp-Meshless cloud method. Computer Methods in Applied Mechanics and Engineering, 139(1):263 – 288, 1996.
- [16] G. R. Liu and Y. T. Gu. An Introduction to Meshfree Methods and Their Programming. Springer Netherlands, 2005.
- [17] X. Liu, G. R. Liu, K. Tai, and K. Y. Lam. Radial point interpolation collocation method (RPICM) for partial differential equations. Computers & Mathematics with Applications, 50(8-9):1425–1442, 2005.
- [18] A. Mardani, M.R. Hooshmandasl, M.H. Heydari, and C. Cattani. A meshless method for solving the time fractional advection-diffusion equation with variable coefficients. Computers & Mathematics with Applications, 75(1):122 – 133, 2018.
- [19] C.M.C. Roque, D. Cunha, C. Shu, and A.J.M. Ferreira. A local radial basis functionsâfinite differences technique for the analysis of composite plates. Engineering Analysis with Boundary Elements, 35(3):363 – 374, 2011.
- [20] C. Shu. Differential Quadrature and Its Application in Engineering. Springer, London, 2000.
- [21] C. Shu, H. Ding, H.Q. Chen, and T.G. Wang. An upwind local RBF-DQ method for simulation of inviscid compressible flows. Computer Methods in Applied Mechanics and Engineering, 194(18):2001 – 2017, 2005.
- [22] C. Shu, H. Ding, and K. S Yeo. Local radial basis function-based differential quadrature method and its application to solve two-dimensional incompressible Navier-Stokes equations. Computer Methods in Applied Mechanics & Engineering, 192(7):941–954, 2003.
- [23] David Stevens, Henry Power, Michael Lees, and Herve Morvan. The use of PDE centres in the local RBF Hermitian method for 3D convective-diffusion problems. Journal of Computational Physics, 228(12):4606 – 4624, 2009.
- [24] A. Tayebi, Y. Shekari, and M.H. Heydari. A meshless method for solving two-dimensional variable-order time fractional advection diffusion equation. Journal of Computational Physics, 340:655–669, 2017.
- [25] Grady B. Wright and Bengt Fornberg. Scattered node compact finite difference-type formulas generated from radial basis functions. Journal of Computational Physics, 212(1):99 – 123, 2006.
- [26] Y. L. Wu and C. Shu. Development of RBF-DQ method for derivative approximation and its application to simulate natural convection in concentric annuli. Computational Mechanics, 29(6):477–485, 2002.
- [27] P. Zhuang, Y. T. Gu, F. Liu, I. Turner, and P. K. D. V. Yarlagadda. Time-dependent fractional advection-diffusion equations by an implicit MLS meshless method. International Journal for Numerical Methods in Engineering, 88(13):1346–1362, 2011.