Radial Basis Functions Based Differential Quadrature Method

for One Dimensional Heat Equation

Kedir Aliyi a, Alemayehu Shiferaw b, Hailu Muleta b

a
Ambo University College of Natural and Computational sciences, Department of

mathematics , Ambo, Ethiopia

b
Jimma University College of Natural Sciences Department of Mathematics , Jimma,

Ethiopia

Corresponding author’s address: E-mail: kediraliyi39@gmail.com

Abstract

In this paper, the Radial basis functions based differential quadrature method has been

presented for solving the one-dimensional heat equation. First, the given solution domain

is discretized and the derivative involving the spatial variable is replaced by the sum of

the weighting coefficients times functional values at each grid point by Multiquadratic

radial basis function. Then, the resulting first-order linear ordinary differential equation

is solved by the 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.

Numerical results are presented in tables in terms of root mean square , maximum

absolute error , and condition number of the system matrix. The numerical

results presented in tables and graphs confirm that the approximate solution is in good

agreement with the exact solution.

1
Keywords: Heat equation, Shape parameter, Multiquadric radial basis functions,

Weighting coefficients, Mesh-free method.

1. 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. 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 [15]. 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 [16]. However, maybe the

complications of these methods include a slow rate of convergence, spatial dependence,

instability, low accuracy, and difficulty of implementation in complex geometries [15].

But scientists in the field of computational mathematics have been trying to develop

numerical methods by using computers for further application. One of those numerical

methods is the differential quadrature method. The differential quadrature method

(DQM) was introduced by Richard Bellman and his associates in the early 1970s

following the idea of integral quadrature [4]. The DQM can produce an accurate result

using only a small number of grid points. They are mesh-based methods. Furthermore,

2
the distribution of the nodes has limitations, and they must be clustered near the boundary

and hence limits its usefulness [15].

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 meshless, which avoids troublesome mesh generation for high-

dimensional problems [15]. These methods are called radial basis function-based

differential quadrature methods. Shu in [20] solved partial differential equations by a

global radial basis function-based 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 [8] extensively tested 29 different algorithms on the typical benchmark function

interpolation problems and ranked the Multi-quadric Radial Basis Function (MQ-RBF)

andThin 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 develope the first multiduadric RBF collocation scheme for solving PDEs of

elliptic, parabolic, and hyperbolic types, in particular, using the MQ-RBFs [19]. The

methodology is now often called the Kansa method. The Kansa method is meshless and

has distinct advantages compared with the traditional methods: superior convergence,

integration-free, and easy implementation. Mesh-less methods 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

3
of the application of the one-dimensional linear parabolic equation, several numerical

methods have been developed.

BenyamMebrate in [1] presented Numerical solutions of a one-dimensional heat Equation

together with the 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 the 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 on 0 < t <1. The method applied is not efficient as it

requires longer CPU time and large memory size. Hooshmandasl et al., in [9] used

Chebyshev Wavelets Method to obtain the anumerical solution of the one-dimensional

heat equation. They applied different wavelet families and the wavelet coefficients were

calculated by the Galerkin or collocation method. Generally, they developed the

Chebyshev wavelet method with operational matrices of integration for the solution of

the 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 relatively large

number of grid points and is difficult to apply for high dimension geometric spaces.

Greengard and Li in [11] 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 the accurate solution for finer spatial step-length and

4
time step. So, the method has required high cost regarding storage capacity in the

computational domain.

Kalyanil and Rao in [10] solved the one-dimensional heat equation by using double

interpolation. They used the finite difference method for double interpolation method to

solve the one-dimensional heat equation, but the method gives better accuracy only for

small step length and is difficult to compute the solution in a 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 the accuracy

of the solution, and efficiency of the method applied.

Tatari and Dehghan in [14] solved the heat equation using radial basis functions. They

applied the Gaussian radial basis functions for obtaining the solution of the 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, this paper

aims to construct an efficient and accurate numerical method for solving a one-

dimensional heat equation.

2. Preliminaries

2.1. Differential Quadrature Method

The differential quadrature method is one of the most efficient numerical methods to

solve partial differential equations. Different differential quadrature methods were

developed by Bellman R. and Castiin 1971 as cited in [17]. A variety of methods have

been developed based on the DQ method, including the polynomial-based differential

5
quadrature (PDQ) and the Fourier-expansion-based differential quadrature methods

(FDQ) [17]. 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 [4]. PDQ and FDQM are highly efficient methods 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

can obtain accurate results using only a small number of grid points, they are mesh-based

methods [7, 15].

To overcome the limitations for the applications of the 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 a 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 the 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 n-dimensional

space [2]..

6
2.2. Radial Basis Function

Definition 2.2.1 A function is called radial if there exists a univariate

function such that was and is a norm in

( is typically the Euclidean norm).

A radial basis function is a one-variable, continuous function defined for that

has been radicalized by composition with the Euclidean norm on . RBFs may have a

free parameter, called the shape parameter, denoted by c. RBFs are a class of radially

symmetric functions symmetric about a point called the center.

It is given in the form of

Here, the value of the univariate function is a function of the Euclidean distance from

the center of point and is given by: in one-

dimensional space[18]. RBFs are a powerful tool in interpolating multivariable functions

or approximation solutions of partial or ordinary differential equations.

The most popular radial basis functions are:

i. Multiquadric RBF:

ii. Inverse Multiquadric RBF:

iii. Gaussian RBF:

iv. Thin plate spline RBF:

7
These radial basis functions (except TPS) are parametric functions with shape (control)

parameter c. They are smooth and continuously differentiable functions. Except for

theMultiquadraic RBF, interpolation matrices of these radial basis functions are positive

definite, but the interpolation matrices for Multiquadraic is conditional positive definite

[18]. The RBFs are usually divided into global and local support RBFs. It is said to be

global support if the and local supporting if the [6]. 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 [18].

Forgiven data the smooth function p defined by enforcing that .

Therefore, for any scattered set of point N with center the radial basis interpolation

function is given as:

By enforcing the interpolation condition , we determine the coefficients

from the matrix for . 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 [13]..

8
 3. Description of the Method, Results, and Discussion

3.1 Description of the Method

Consider the one-dimensional heat equation of the form:

, (3.1)

subject to initial condition:

, (3.2)

and boundary conditions

, (3.3)

Here, , and are continuous and differentiable functions. The

computational domain is partitioned as:

, where and are mesh-size of and , respectively.

3.2. Multiquadric Radial Basis Functions (MQ-RBFs)

Multiquadric RBF was proposed by Hardy as cited by Ding et al. in [5] for the

interpolation of topographical surfaces. It has a very simple mathematical form as

, where are a shape parameter and >0. For time-dependent problems, the

multiquadricradial basis function is given by

, (3.4)

9
, where is the number of grid points in the direction of the spatial variable , is the

temporal variable.

Definition 3.2.1. (Chenoweth, 2012).

A function is completely monotonic on if and

where and

Theorem 3.2.1 (Chenoweth, 2012)

Let and for Let be completely monotonic

on , then for any set of the distinct center , the matrix B

with an entry is invertible. Such function is said to be a conditional

positive definite.

MQ-RBF is an example of a global infinitely differentiable and conditionally positive

definite function.

Theorem 3.2.2 (Micchell’s, in [12])

For a set of distinct data points to be approximated by the set of 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:

(3.5)

10
 (3.6)

3.3 Approximating the Partial Derivatives and Determining the Weighting

Coefficients

From the primary idea of the differential quadrature method, Eq. (3.1) can be

expressed as a linear combination of functional values at each grid point. It is given as:

, (3.7)

Here, are the weighting coefficients. Since is a basis function for one-

dimensional space, every function defined on one-dimensional space can be expressed as

a linear combination of these basis functions using the properties of linear independence

of vector in vector space. Now for any constant and at least one of such that

, (3.8)

Since = the first and second-

order derivatives of Eq. (3.8) become

, (3.9)

, (3.10)

Substituting Eqs.(4.8) and (3.10) into Eq. (3.7), we get:

(3.11)

11
 (3.12)

From Eq. (3.12) we have the following system of equations:

,(3.13)

, and

Since is global support of RBFs it produces a dense matrix A. By theorems 3.2.1 and

3.2.2 the matrix A is conditional positive definite and non-singular. Therefore, the system

in Eq. (3.13) has a unique solution. The condition number of A is calculated as

(3.14)

For a large number of , the system in Eq. (3.13) is ill-conditioned. This is because of

the presence of shape parameter “ ” that affects both the condition number of the

interpolation matrices and the accuracy of the method. For a fixed number of

interpolation points N, the condition number depends on the shape parameter and

12
support of the RBFs. Also, the condition number grows with for fixed value of the

shape parameter c [14].

Now the weighting coefficients (w) can be obtained from the linear system of the

equation given in Eq. (3.13) by using the Gaussian-Elimination method. In this case, the

augmented matrix is written as

(3.15)

, where and for

By row operation, the augmented matrix is changed to the upper triangular matrix.

(3.16)

Now, we can calculate the weighting coefficients using backward substitution.

The weighting coefficient can be obtained by using the formula:

(3.17)

The weighting coefficient can also be obtained by using the formula:

(3.18)

Generally, the weighting coefficient can be obtained by using the formula:

13
 (3.19)

for

Then, the coefficient vector can be used to approximate the second-order derivative

in the x-direction for unknown smooth function at 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 [5].

3.4. Results and Discussion

From Eqs. (3.1-3.3), (3.8), (3.10), and (3.12), the resulting system of initial-boundary

value problem is given as:

, (3.20)

subject to the initial condition

, (3.21)

and boundary conditions

(3.22)

3.5. Criteria for Investigating the Accuracy of the Method

The Root Mean Square (RMS) error ( ), maximum absolute error ( ) are used to

measure the accuracy of the method. The RMS error and maximum absolute error are

calculated as followed by Tatari and Dehghan,in [14].

14
 , (3.23)

(3.24)

Here, and are the exact and approximation solutions of Eqs. (3.1),

(3.2), and (3.3), respectively.

3.6 Numerical Experiments

Example1: Consider the classical heat equation considered by Dehghan&Tatari, in [14]

, , ,

with initial condition

and boundary conditions

, ,

The exact solution is given by:

The numerical results are presented in tables in terms of and as the means

for measuring the accuracy of the present method.

Table 1.Point wise absolute error, root mean square error, and Condition Number of

example1 when h=0.1, and 0.01.

Our Method

Shape parameter c

0.5 3.2888E-01 3.3052E-0 4.0762E+07

15
 1 6.2787E-02 6.3100E-01 1.0457E+12

2 1.5974E-02 1.6054E-01 8.0464E+16

3 7.2063E-03 7.2422E-02 2.6283E+17

4 4.0799E-03 4.1003E-02 6.4077E+16

5 2.6197E-03 2.6327E-02 2.9761E+17

6 1.8226E-03 2.6327E-02 3.5931E+17

7 1.3405E-03 1.3472E-03 7.3104E+17

Tatari, M. and Dehghan,M.,2010

Shape parameter c

0.5 9.0651E+39 7.0000E+73 1.9676E+19

1 3.3136E+17 2.0000E+29 2.0434E+17

2 5.0227E+1 5.0001E+0 9.3751E+17

3 1.0196E-01 7.1266E-03 6.6616E+17

4 1.0384E-01 7.4346E-03 4.0811E+17

5 1.0400E-03 7.6029E-03 1.1529E+16

6 1.0150E-03 7.5571E-03 1.4832E+16

7 9.5715E-03 7.3241E-03 3.0267E+15

Table 2. Pointwise absolute error, root mean square error, and Condition number of

example1 when h=0.1, and 0.025

Our Method

Shape parameter c

0.5 1.7625E+01 1.1286E+02 4.0762E+07

16
 1 3.3389E-01 3.3556E-00 1.0457E+12

2 7.5336E-02 3.3389E-01 8.0464E+16

3 3.0695E-02 1.9655E-01 2.6283E+17

4 1.6767E-02 1.0736E-01 6.4077E+16

5 1.0588E-02 6.7799E-02 2.9761E+17

6 7.3003E-03 4.6745E-02 3.5931E+17

7 5.3403E-03 3.4195E-02 7.3104E+17

Table 3. Pointwise absolute error, root mean square error, and Condition number of

example1 when h=0.1, and 0.001.

Our Method

Shape parameter c

0.5 2.9601E-03 9.3653E-02 4.0762E+07

1 1.3346E-03 4.2226E-02 1.0457E+12

2 4.4049E-04 1.3937E-02 8.0464E+16

3 2.0921E-04 6.6190E-03 2.6283E+17

4 1.2064E-04 3.8168E-03 6.4077E+16

5 8.1757E-05 2.5867E-03 2.9761E+17

6 5.7186E-05 1.8093E-03 3.5931E+17

7 4.2198E-05 1.3351E-03 7.3104E+17

17
Figure 1: Graph of the Approximate and Exact solution at &

a. Surface sketch of Exact solution b. Surface sketch of Approximate solution

1 1

0.9 0.9

0.8 0.8

0.7 0.7

0.6 0.6
u(x,t)-->

0.5 0.5
u(x,t)

0.4 0.4

0.3 0.3

0.2 0.2

0.1 0.1

0 0
12 12
10 10
120 120
8 100 8 100
6 80 6 80
4 60 4 60
40 40
2 2
20 20
x 0 0
t x 0 0
t

Figure 2: Surface graphs of example 1 showing the physical behavior of the one-

dimensional heat equation when &

3.6.1. Discussion

18
As depicted in Table 1, the present method can generate a convergent numerical solution

when c = 0.5, 1, and 2 at which the method presented by Tatari and Dehghan, 2010 fails

to produce a convergent solution. The condition number of the system matrix of the

present method is in the range 4.0762E+07 7.3104E+17 whereas the condition

number of the system matrix presented by Tatari and Dehghan, 20103.0267E+15

1.9676E+19. Thus, the effect of the condition number on the accuracy of the numerical

solution is more significant on the method presented by Tatari and Dehghan, 2010 than

on the numerical solution of the present method. The value in Table 1 confirms this

issue. That is the smaller the value of 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 the same values of the

shape parameter c in each table. This is because the condition number depends only on

the step length the spatial variable and the shape parameter c. Comparison among Table

1-Table 3 shows that a more accurate result is generated when 0.001. This indicates

that the present method is suitable for problems that require a long time interval.

The simulations presented in Figures 1 and 2 show that the approximate solution obtained

by the present method is in good agreement with the exact solution.

3.6.2. Conclusion

In this thesis, RBFS-DQM is used to solve the 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

19
value at each grid point via Multiquadraic radial basis functions. Then the resulting first-

order linear system of ODE is solved by the 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 step-length h

fixed. As can be seen from the numerical results presented in tables and graphs, the

present method is superior to 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 the one-dimensional heat equation.

Acknowledgment

The authors are grateful to Jimma University for fully funding the project work so that it

has become reality.

REFERENCES

[1] BenyamMabrate (2015). Numerical solution of a one dimensional Heat Equation with

Dirichlet boundary conditions.American Journal of Applied Mathematics,

3(4):305-311.

[2] Boling, E.F. (2013). Multi-GPU solutions of geophysical PDEs with radial basis

Function-generated finite differences.PhD thesis Florida State University.

[3] Chen, W., Fu, Z. J., and Chen, C. S. (2014). Recent advances in radial basis function

collocation methods, Springer Briefs in Applied Sciences and Technology.

[4] Cheong, H.T. (2015). Parallel Localized differential quadrature approach in boundary

value problem.Journal of math. Analysis and appl. 31(2):127-134.

20
[5] Ding, H., Shu, C., and Tang, D.B. (2005). Error estimates of local multiquadric-based

differential quadrature (LMQDQ) methodthrough numerical experiments.

International journal for numerical method in Engin. 63(11):1513-1529.

[6] Erfanian, M., and Kosari, S. (2017). Using Chebyshev polynomials zeros as point grid

for numerical solution of nonlinear PDEs by differential quadrature-based radial basis

functions. Int.Journ.of Maths. Modelling and Computations, 7(1):67-77.

[7] Fasshauer, G. (2007). Meshfree Application Method with Matlab.Interdisciplinary

Mathematical Sciences.

[8] Franke, R. (1982).Scattered data interpolation: tests of some methods. Mathematicsof

computation, 38(157):181-200.

[9] Hooshmandasl M. R Haidari, M., and MaalekGhaini F.M. (2012).Numerical solution of one

dimensional heat equation by using Chebyshev Wavelets method.JACM an open access

journal,1(6), 1-7.

[10] Kalyanil, P., and Rao (2013).Numerical solution of heat equation through double

Interpolation.IOSR Journ.of math.6(6): 58-62.

[11] Li, J. R., and Greengard, L. (2007). On the numerical solution of the heat equation I:

Fast solvers in free space. Journal of Computational Physics, 226(2):1891-

1901.

[12] Micchelli, C.A (1984). Interpolation of scattered data: distance matrices and

conditionally positive definite functions. In Approximation theory and spline

functions:143-145. Springer, Dordrecht.

[13] Mongillo, M. (2011). Choosing basis functions and shape parameters for radial basis

function methods.SIAM undergraduate research online, 4(190-209):2-6.

[14] Tatari, M. and Dehghan,M.,(2010).A method for solving partial differential equations via

21
 radial basis functions:Application to the heat equation Engineering Analysis with

Boundary element 34(3),206-212.

[15] Waston, D. (2017). Radial Basis Function Differential Quadrature Method for the

numerical solution of partial differential equations. The Aquila Digital

Community, Dissertations.

[16] Yao,G. and Yu ,Z. (2012).A Localized Meshless Approach for Modeling Spatial temporal

Calcium Dynamics in Ventricular Myocytes.International Journal for numerical method in

biomedical engineering.28 (2),187-204.

[17] Zhang, Y., and Zong, Z. (2009).Advanced differential quadrature methods.Chapman

and Hall/CRC.USA.

[18] Chenoweth, M.N. (2012). A local radial basis function method for the numerical

Solution of partial differential equations.PhD thesis, Marshall University.

[19] Kansa, E. J. (1990). 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-9), 127-145.
[20] Shu,C., Ding, H., and Yeo, K. S. (2004). The solution of partial differential
equations by a global radial basis function-based differential quadrature
method. Engineering Analysis with Boundary Elements,28(10),1217-1226.

22