Radial Basis Functions Based Differential Quadrature Method for One Dimensional Heat Equation

: In this paper, Radial basis functions based differential quadrature method has been presented for solving one-dimensional heat equation. First, the given solution domain is discretized using uniform discretization grid point in both direction and the derivative involving the spatial variable, x is replaced by the sum of the weighting coefficients times functional values at each grid points. Next by using properties of linear independence of vector space and Multiquadratic radial basis function we can find all waiting coefficient at each grid points of solution domain and we obtain first order linear system of ordinary differential equation with N by N square coefficient Matrices. Then, the resulting first order linear ordinary differential equation is solved by fifth-order Runge-Kutta method. To validate the applicability of the proposed method, one model example is considered and solved for different values of the shape parameter ‘c’ and mesh sizes in the direction of the temporal variable; t and fixed value of mesh size in the direction of spatial variable, x. Numerical results are presented in tables in terms of root mean square (E 2 ), maximum absolute error (E ∞ ) and condition number K (A) of the system matrix. The numerical results presented in tables and graphs confirm that the approximate solution is in a good agreement with the exact solution.


Introduction
Partial Differential Equations (PDEs) are mathematical equations that are significant in modeling physical phenomena that occur in nature. Applications of PDEs can be found in physics, engineering, mathematics, and finance. Examples include: modeling mechanical vibration, heat, sound vibration, elasticity, and fluid dynamics [7]. Although PDEs have a wide range of applications to real world problems in science and engineering, the majority of PDEs do not have analytical solutions.
Traditionally, mesh-dependent methods such as the finite difference method (FDM), finite element method (FEM), and boundary element method (BEM) have been used to solve PDEs [16]. Despite their great success in the past decades in many branches of science and engineering, these mesh-based methods require meshes or grids as the solution domain. The costs and difficulties in creating quality meshes, however, constitute one of the major bottlenecks in these methods [17]. However, may be the complications of these methods include a slow rate of convergence, spatial dependence, instability, low accuracy, and difficulty of implementation in complex geometries [16,18].
But scientists in the field of computational mathematics have been trying to develop numerical methods by using computers for further application [18]. One of those numerical methods is differential quadrature method. The differential quadrature method (DQM) was introduced by Richard Bellman and his associates in the early of 1970s following the idea of integral quadrature [19,14]. The DQM are able to produce the accurate result using only a small number of grid points. They are mesh-based methods. Furthermore, the distribution of the nodes has limitations, One Dimensional Heat Equation and they must be clustered near the boundary and hence limits its usefulness [16].
However, mesh-less approximation techniques using radial basis functions (RBFs) have been developed over the last several decades. These methods are easy to implement, highly accurate, and truly mesh less, which avoids troublesome mesh generation for high-dimensional problems [20]. These methods are called radial basis function based differential quadrature methods [20]. Shu et. al. Solved partial differential equations by a global radial basis functionbased differential quadrature method.
In recent decades, various RBFs-based methods have gained fast growing attention from abroad range of scientific computing and engineering applications such as multivariate scattered data processing, numerical solutions of PDEs, and machine learning [3]. Franke in 1982 [9] extensively tested 29 different algorithms on the typical benchmark function interpolation problems, and ranked the Multi-quadric Radial Basis Function [21] (MQ-RBF) and Thin Plate Spline Radial Basis Function (TPS-RBF) as two of the best candidates based on the following criteria: timing, storage, accuracy, and ease of implementation. Kansa (1990) first developed an RBF collocation scheme for solving PDEs of elliptic, parabolic, and hyperbolic types, in particular, using the MQ-RBFs. The methodology is now often called the Kansa method. The Kansa method is mesh less and has distinct advantages compared with the traditional methods: superior convergence, integration-free, and easy implementation. Mesh-less method are a class of numerical methods used for solving partial differential equations. In these methods, mesh generation on the spatial domain of the problem is not needed. This property is the main advantage of these techniques over the mesh dependent methods. Due to the wide range of the application of the one dimensional linear parabolic equation, several numerical methods have been developed.
Benyam Mebrate [1] presented Numerical solutions of a one-dimensional heat Equation together with initial condition and Dirichlet boundary conditions of the form: He presented computation of the numerical solutions by using two methods: Finite difference, and Finite element methods. The finite element methods are implemented by Crank-Nicolson method. This method does not always converge to the exact solutions for coarser step lengths. He got the accurate approximate solution with the length of time-step k=0.001 on0 < t <1. The method applied is not efficient as it requires longer CPU time and large memory size. Hooshmandasl et al., [10] used Chebyshev Wavelets Method to obtain a numerical solution of the onedimensional heat equation. They applied different wavelet families and the wavelet coefficients were calculated by the Galerkin collocation method. Generally, they developed Chebyshev wavelet method with operational matrices of integration for solution of one dimensional heat equation with Dirichlet boundary conditions which is fast, mathematically simple and guarantees the necessary accuracy for a relatively small number of grid points. This confirms that the method is not accurate for a relative large number of grid points and is difficult to apply for high dimension geometric spaces. In their work, Greengard and Li [12] also developed an explicit method for solving inhomogeneous heat equation in free space, following the time evolution of the solution in the Fourier domain. The error in the solution is simply the quadrature error in evaluating the solution and the solution is dependent on the full time-space history of the diffusion process with the time-step k=0.0025 and they got accurate solution for finer spatial step-length and time step. So, the method has required high cost regarding to storage capacity in the computational domain.
Kalyanil and Rao in [11] solved one-dimensional heat equation by using double interpolation. They used finite difference method for double interpolation method to solve onedimensional heat equation, but the method gives better accuracy only for small step length and is difficult to compute the solution in complex computational domain. Even though the accuracy of the aforementioned methods is promising, they require large memory and long computational time. Besides, the methods are not suitable for higher dimensional and problems involving complex geometries. So, the treatment of the mesh-size presents severe difficulties that have to be addressed to ensure accuracy of the solution, and efficiency of the method applied.
Tatari and Dehghan [15] solved heat equation using radial basis functions. They applied the Gaussian radial basis functions for obtaining solution of heat equation. The solution diverges when shape parameters are: 0.5, 1, and 2. As compared to the exact solution, the approximate solution needs further improvement. To this end, the aim of this paper is to construct efficient and accurate numerical method for solving one-dimensional heat equation.

Differential Quadrature Method
Differential quadrature method is one of the most efficient numerical methods to solve partial differential equations [18]. Different differential quadrature methods were developed by Bellman R. and Castiin 1971 as cited in [19]. A variety of methods have been developed based on the DQ method, including the polynomial-based differential quadrature (PDQ) and the Fourier-expansion-based differential quadrature methods (FDQ) [19] The basic idea behind the DQ method is to determine the weighting coefficients for any order derivatives by using a weighted sum of functional values at a set of selected grid points [5]. PDQM and FDQM are highly efficient method by using a small number of grid points, they are not efficient when the number of grid-points is large and they are also sensitive to grid point distribution. While the PDQ and FDQ methods are able to obtain accurate results using only a small number of grid points, they are mesh-based methods [8,16].
To overcome the limitations for the applications of DQ method, a new class of methods called mesh-free methods has surfaced. Each class of methods offers numerous and, in many ways, complementary benefits. In the ideal case we seek a method defined on arbitrary geometries that behaves regularly in any dimension, and avoids the cost of mesh generation. The ability to locally refine areas of interest in a practical fashion is also desirable. Fortunately, mesh-free methods provide all of these properties: based wholly on a set of independent points in n-dimensional space, there is minimal cost for mesh generation, and refinement is as simple as adding new points where they are needed. A subset of mesh-free methods of particular interest to the numerical modeling community today revolves around Radial Basis Functions (RBFs). This method is called radial basis function differential quadrature method (RDQM). RBF methods are based on a superposition of translates of these radially symmetric functions, providing a linearly independent but non-orthogonal basis to interpolate between nodes in ndimensional space [2].  (Chenoweth, 2012). RBFs are a powerful tool in interpolating multivariable functions or approximation solution of partial or ordinary differential equations. The most popular radial basis functions are:
The global approach uses information from every center in the domain to approximate a function value or derivative at a single point. In contrast, the local method only uses a small subset of the available centers [4].
For given data ( . ) j j x f the smooth function p defined by enforcing that ( ) . The choice of the basic functions will determine which methods are available for solving system of interpolation and whether such a solution even exists. If the interpolation matrix is symmetric positive definite, then the linear system has a unique solution [14].

Description of the Method
Consider the one dimensional heat equation of the form: Subject to initial condition: and boundary conditions Here, 1 ( ) f x , 1 ( ) g t and 2 ( )

Multiquadric Radial Basis Functions (MQ-RBFs)
Multiquadric RBF was proposed by Hardy as cited by Ding et al., in [6] for the interpolation of topographical surfaces. It has very simple mathematical form as where c is a shape parameter and c >0. For time dependent problems, the multiquadric radial basis function is given by where N is the number of grid points in the direction of the spatial variable x , t is the temporal variable.
Definition -1: [4] A function ϕ is completely monotonic MQ-RBF is an example of global infinitely differentiable and conditional positive definite function.
Theorem -2: [13]. For a set of N distinct data points to be approximated by the set of N approximation function φ and this function must be nonsingular. Therefore, by the above definition and theorems, the interpolation matrix resulting from MQ-RBFs is invertible, smooth, and continuously differentiable on its domain. Thus, the first and second order derivatives of multiquadric radial basis functions are given by:

Approximating the Partial Derivatives and Determining the Weighting Coefficients
From the primary idea of differential quadrature method, Here, ij w are the weighting coefficients. Since ( ) Substituting Eqs. (8) and (10) into Eq. (7), we get: From Eq. (12) we have the following system of equations: Since φ is global support RBFs it produces a dense matrix A. By theorems 1 and 2 the matrix A is conditional positive definite and non-singular. Therefore, the system in Eq. (13) has a unique solution. The condition number of A is calculated as For large number of ( ) A κ , the system in Eq. (13) is illconditioned. This is because of the presence of shape parameter" c "that affects both the condition number of the interpolation matrices A and the accuracy of the method. For a fixed number of interpolation points N, the condition number of A depends on the shape parameter c and support of the RBFs [22]. Also, the condition number grows with N for fixed value of the shape parameter c [15]. Now the weighting coefficients (w) can be obtained from linear system of equation given in Eq. (13) by using Gaussian-Elimination method. In thiscase, the augmented matrix is written as [ ] By row operation, the augmented matrix is changed to the upper triangular matrix.
Now, we can calculate the weighting coefficients using backward substitution.
The th N weighting coefficient N w can be obtained by using the formula: The ( ) 1 th N − weighting coefficient 1 N w − can also be obtained by using the formula: Generally, the th i weighting coefficient can be obtained by using the formula: can be used to approximate the second-order derivative in the x-direction for unknown smooth function at th i nodes. From the procedure of DQ approximation of derivatives, it can be observed that the weighting coefficients are only dependent on the selected RBFs and the distribution of the supporting points in the local support [6].

Results and Discussion
From Eqs. (1-3), (8), (10) and (12), the resulting system of initial-boundary value problem is given as: Subject to initial condition and boundary conditions

Criteria for Investigating the Accuracy of the Method
The root Mean Square (RMS) error ( 2 E ), maximum absolute error ( E ∞ ) are used to measure the accuracy of the method. The RMS error and maximum absolute error are calculated as follows in [15],

Numerical Experiments
Example1: Consider the classical heat equation considered by Dehghan &Tatari, in [15], The exact solution is given by: The numerical results are presented in tables in terms of 2 , as the means for measuring the accuracy the present method.

Discussion
As depicted in Table 1, the present method is able to generate convergent numerical solution when c = 0.5, 1, and 2 at which the method presented by Tatari and Dehghan, in [15] Table 1 confirms this issue. That is the smaller the value of E ∞ the less the effect of the condition number on the accuracy of the approximate solution. As it can also be seen from Tables 1-Table 3, the condition number of the system matrix is kept constant for same values of the shape parameter c in each table. This is because the condition number is depends only on the step length the spatial variable and the shape parameter c. Comparison among Table 1- Table 3 shows that more accurate result is generated when t ∆ = 0.001. This indicates that the present method is suitable for problems that require long time interval. The simulations presented in Figures 1 and 2 shows that the approximate solution obtained by the present method is in a good agreement with the exact solution.

Conclusion
In this thesis, RBFS-DQM is used to solve one dimensional Heat equation. First, the domain is discretized using the uniform step length and the derivative involving the spatial variable 'x' is replaced by the sum of weighting coefficient times the functional value at each grid point via Multiquadraic radial basis functions. Then the resulting first order linear system of ODE is solved by fifth order Runge-Kutta method.
To validate the applicability of the method, one model example is considered and solved by varying the value of shape parameter c and time-step k, and keeping the steplength h fixed. As it can be seen from the numerical results presented in tables and graphs, the present method is superior over the method developed by Tatari and Dehghan (2010) and approximates the exact solution very well. In a nutshell, the present method is conceptually simple, easy to use and readily adaptable for computer implementation for solving one-dimensional heat equation.