Solution to a Convection-Diffusion Problem Using a New Variable Inverse-Multiquadric Parameter in a Collocation Meshfree Scheme

The lack of reliable judgment on the choice of shape parameter is acknowledged as a drawback in the methodology of collocation meshfree by radial basis function (RBF). Some attempts have been proposed and tested only with specific classes of PDEs. Moreover, while most of this focusses on multiquadric (MQ) RBF, there has not been work done on invers-multiquadric (IMQ) RBF despite its’ increasing popularity. As a consequence, our main tasks in this work are, firstly, to numerically investigate the quality of each adaptive/variable RBF-shape parameter approaches presented in literature by applying them to the same type of problem, convection-diffusion class. Secondly, we proposed a new form of shape-parameter scheme to be used with the inverse-multiquadric (IMQ) type of RBF. The Kansa meshless method is implemented and it is interestingly found that the proposed form produces good results quality in terms of both matrix condition number, and the accuracy.


INTRODUCTION
Convection-diffusion equation describes phenomena including both convection and diffusion effects, and appears in various fields of natural sciences, e.g., heat transfer, weather prediction and atmospheric radioactivity propagation.It may also be treated as a simplified model of the system of the Navier-Stokes equations, which are representative equations in fluid dynamics.For this kind of physical phenomena, numerical studying of it is known as one of the great challenges due to the interaction between convective and diffusive processes normally giving rise to oscillatory and nonphysical solutions particularly when Peclet numbers are increasing or convective forces become dominant.Amongst several attempts and remedies; the enlargement of the local node adaptive, upwind scheme, the biased domain, and the nodal refinement and the adaptive analysis, for this challenge have been proposed, revisited, and nicely numerically demonstrated by Gu and Liu [1].In this work, they concluded that with extra efforts included, the meshfree method has some attractive advantages over the traditional schemes such as finite element and finite volume.A modified finite-difference with the technique of elimination building up a new Parameter in a Collocation Meshfree Scheme iterative scheme to solve the implicit difference equation was wonderfully presented in Zhang and Wang [2].With compared to the Jacobi iteration, it was found that their new approach provides higher rate of convergence when dealing with moderate convectivediffusive problems.Recently, in 2013, Sun and Li [3] proposed a method which combines a 6th-order compact difference scheme and 2nd-order Crank-Nicolson scheme called the alternating direction implicit method (ADI).Their numerical examples have proven that ADI preserves the higher order accuracy for convection-dominated problem.Moreover, one of our previous works, Chanthawara et.al [4], done under the context of the boundary element method has confirmed the difficulty encountered and this is the main reason for focusing on this challenging problem in this investigation.
Initially invented and developed for multivariate data and function interpolation, radial basis function schemes have now become an attractive numerical tool for tackling engineering problems modelled in the form of partial differential equation (PDEs).While all the traditional numerical methods; finite volume (FV), finite element (FE), and finite difference (FD) are suffering from the grid/mesh generation process which consumes a great amount of effort, these radial basis function methods are classified as 'meshfree' or 'meshless' methods as no grid/mesh is not required in the discretization process.With this fact, the methods have gained their popularity for their simplicity to implement and, moreover, their capability of dealing with multi-dimensional problems.In early 1990s, Kansa [5] used a set of global approximation function to approximate the field variables on both the domain and the boundary when solving PDEs.This is of the following expression; ( ) Over some given scattered data nodes , Ω is the problem domain.
The Radial Basis Functions (RBF), ϕ , are commonly found as multivariate functions whose values are dependent only on the distance from the origin and commonly assumed to be strictly positive definite.This means that ( ) ( ) x  and r∈  ; or, in other words, on the distance from a point of a given set { } j x , and ( ) ( )  where it can normally be defined as follows; For some fixed points represents the Euclidean distance.The whole class of meshfree/meshless method can be categorized into three classes; weak forms, strong forms, and mixed, all nicely documented in [6].Each of these has its own advantages/disadvantages depending on several factors involved including domain geometry, governing equations, boundary/initial conditions, computer arithmetic, etc.For RBFcollocation meshfree methods, in particular, it is a common fact that the accuracy depends heavily on the behavior of the radial basis function (RBF), ϕ , used itself.Many popular choices of RBFs have been proposed, studied, and developed during the past decades.Some of these are listed in Table 1.
K is an order modified Bessel function In this work, nevertheless, the RBF type being focused on is known as inversemultiquadric form expressed as follows; ( ) ( , ) r r where ..., 3 / 2, 1/ 2, 1/ 2, 3 / 2, ... β = − − and ε is the so-called 'shape parameter' is normally given in a 'ad-hoc' manner.This kind of radial basis function has been proven to be strictly positive definite as nicely documented in [7] for which the distance matrix of the interpolation problem is invertible.Nevertheless, throughout this work we focus on the most popular form of inverse-multiquadric, i.e.
1/ 2 β = − , only.This shape parameter is known to have a huge effect on the solution quality but it is, very often, given in an 'ad-hoc' manner.The essential effect of this parameter has been regarded as the main key leading to the final result quality [8,9,10].Hardy [11] suggests that by fixing the shape at 1/ (0.815 ) , and i d is the distance from the node to its nearest neighbor, good results should be anticipated.Also, in the work of Franke [12] where the choice of a fixed shape of the form 0.8 N D ε = where D is the diameter of the smallest circle containing all data nodes can also be a good alternative.Rippa [10] studied the selection of optimal shape parameter for RBFs and oncluded that in order to obtain the socalled 'reliable shape parameter', it is crucial to take in to consideration many factors; the number and distribution of data points, radial basis function, condition number of coefficient matrix and precision of computation into account.For any particular interpolation problem, the radial basis function and precision of computation remains similar throughout the domain.However, if the distribution of data points is not uniform, the optimal value of shape parameter will differ for each data point in local RBFs and would depend upon the number and distribution of data points within its own influence domain.
Some recent attempts to pinpoint the optimal value of involve the work of Zhang et al. [13] where they demonstrated and concluded that the optimal shape parameter is problem dependent.In 2002, Wang and Liu [14] pointed out that by analyzing the condition number of the collocation matrix, a suitable range of derivable values of can be found.Later in 2003, Lee et.al. [15] suggested that the final numerical solutions obtained are found to be less affected by the method when the approximation in equation ( 1) is applied locally rather than globally.
Numerical studies referred to until now have been mostly on somewhat 'fixed value' of the parameter.Alternatively, some works have recently focused on developing this shape in adaptive-or variable-manner in which the detail is provided in Section (3) of this work.Our investigation carried out in this work aims to complete three main objectives.Firstly, it is to gather, numerically test, and compared against one another, some forms of variable shape parameters proposed in literature.Secondly, it is to propose another mathematical formula of variable shape aimed to apply with inverse-multiquadric RBF.Last but not least, in order to justify the effectiveness of our newly-proposed inverse-multiquadric shape parameter with those gathered from literature, we implement all forms to the same type of PDEs known as convection-diffusion via the RBF-collocation method.
The layout of the work is as follows.We firstly introduce in Section (2) the methodology of the collocation methods using radial basis functions for partial differential equations.It is followed by the implementation of the method toward convection-diffusion type of PDEs.In Section (3), we provide a brief review on numerical works done on variable shape parameters.We also, in this section, propose a new form of parameter that behaves also variably.The main numerical experiment is then detailed in Section ( 4) together with all the results obtained and general discussion.The validation of the solutions produced from every forms of shape parameters is done via its analytical form.Some important findings or conclusions found in this work are then withdrawn in Section (5).

RBF-Collocation Scheme for PDEs
For the methodology of RBF-collocation meshless method for numerically solving PDEs, it begins with considering a linear elliptic partial differential equation with boundary conditions, where ( ) g x and ( ) f x are known.We seek ( ) where d Ω ∈  , ∂Ω denotes the boundary of domain Ω , L and M are the linear elliptic partial differential operators and operating on the domain Ω and boundary domain ∂Ω , respectively.For Kansa's method, it represents the approximate solution ( ) u x  by the interpolation, using an RBF interpolation as expressed in equation (1).We can see that N linear dependent equations are required for solving N unknowns of j c .Substituting ( ) u x  into equation (4) and equation ( 5), we obtain the system of equations as follows; , , x where 1, , Above equations, we choose N collocation points on both domain Ω and boundary domain ∂Ω , and divide it into I N interior points and B N boundary points ( )  6) and equation ( 7) are chosen as collocation points.The previous substituting yields a system of linear algebraic equations which can be solved for seeking coefficient ' s c by rewriting equation ( 6) and equation ( 7) equation in matrix form as; where , , , and ( ) ( ) 8), the coefficient ' s c are computed from the following system; Therefore, the matrix c is substituted into equation ( 1) and the approximate solution of ( ) The system is known to provide solution if and only if the matrix A is non-singular, its inverse exists.This aspect is related directly to its condition number, as can be defined as; In the RBF-meshless method, it is well-known that this condition number is strongly affected by the magnitude of the shape parameter ε and the number of nodes involved, N . For this, the solutions obtained in this work are also monitored in terms of this condition number ensuring the solvability of the collocation matrix.

Implementation towards Convection-Diffusion Type of PDE
Consider the governing partial differential equation of convection diffusion problems expressed as; In this practice, its boundary condition is of Dirichlet type as follows: In order to implement the RBF-collocation meshless procedure, the above governing equation can be approximated using the summation and becomes; This system can be written in the form of; This system can be generated in matrix form and the approximate solutions then can be obtained by substituting the coefficients , j c obtained by solving the above matrix form, in equation (1).

THE PROPOSED IMQ VARIABLE SHAPE PARAMETER
It is known that the effectiveness of the RBF-collocation method can well be influenced by the choice of the shape parameter.Many researchers have agreed that using variable shape parameters provides superior in solution quality (See [16] and references herein).The result of allowing the shape parameter to vary locally is that each column of the interpolation matrix, matrix A in (8), no longer contains constant entries leading to lower condition number.Several strategies have been proposed for providing reliable numerical solution accuracy and they are revisited as shown in Table 2.
The variable abbreviated as VAR-1 is clearly in an exponential manner and was used in Kansa [1] before its further modified version was later invented in Kansa and Carlson [17], noted as VAR-2.In their work, it was demonstrated that if 2 min ε and 2 max ε varied by several orders of magnitude, then a very satisfactory result quality can well be expected.Later, a linear form of variable shape parameter was proposed and applied to both interpolation and some benchmark partial differential equations by Sarra [18] and Sarra and Sturgill [16], as noted in Table 2 by VAR-3 and VAR-4 respectively.Here, the command 'rand' is the MATLAB function that returns uniformly distributed pseudo-random numbers on the unit interval.It was proven in their work that the variable shape outperformed the fixed value of parameter especially when the scheme includes the information about the minimum distance of a center to its nearest neighbor, n h , with also a user input value µ .In terms or the condition number produced by equation (11), it was also found to be considerable smaller over most of the average shape range.

Abbreviation used in this work
Reference Formulation of ε for jth-element Kansa and Carlson [11] max min min ( ) Sarra and Sturgill [10] ( ) In this work, we proposed a new form of variable shape parameter where both linear and exponential manners are taken into consideration, expressed as in equation (18).
( ) This weight then sets j ε to become its exponential manner when j N = .This proposed variable shape, equation ( 18) is referred to as " VAR-5 " throughout the work.

NUMERICAL EXPERIMENT AND RESULTS
To demonstrate how effective the variables VAR-1, VAR-2, VAR-3, VAR-4 and VAR-5 can be, let us consider a 2D convection-diffusion problem in steady state, as given and studied in Gu and Liu [19].The governing equation is expressed as follows; The computational domain is taken to be ( and the coefficients are 0 , 0 In order to cover as wide aspect as possible, a large number of experiments were carried out but only some is presented here.Table 3 -Table 5 show the optimal shape parameter value obtained at different levels of diffusion coefficients ( γ ), the number of nodes ( N ), and the range of min max ( , ) ε ε .At relative high value of 10 γ = , it is found that the proposed form of variable shape yielded the same level of solutions quality as those obtained from other forms gathered from literatures.This is still the case even when the number of nodes and the min max ( , ) ε ε range are increasing.It is interesting to see that the proposed variable shape is still numerically compatible with the others even when the problem becomes more convective-dominated, , where the same order of error norms magnitudes are found as seen in Table 4 and Table 5.In terms of the condition number of the interpolation matrix, Figure 2 shows the plots where the condition number is seen to decrease around the optimal value of ε .The larger ε causes ( ) δ Λ A , as expected, to decrease meaning that the numerical process is farther from being affected by the 'singularity' problem.With high value of 1 γ > , solution profiles obtained from every variable shape under this investigation look similar, as shown in Figure 2. The effect of the number of collocation nodes is shown in Figure 4.As expected, the clearest effect takes place at the corner of the domain, 0.9 , 1.0 x y < < , where the boundary layers form and grow.The smoother solution profiles around this area as shown can be attributed to better support gained from local denser collocation nodes via. the collocation scheme used.
In order to compare the results obtained in this study to one of the benchmarks in literature namely Gu and Liu [19], and also to some of our pervious works, a new error indicator (Err), is adopted and defined as; Err ( ) Table 6 contains the error, defined in equation ( 23), comparing amongst those obtained from this study, from our previous work, and those from Gu and Liu [19].NAA* and IMQ* represent respectively, errors produced by the use of another popular numerical method called 'dual reciprocity boundary element method (DRBEM)' in which RBF plays also a crucial role.In this work, the optimal shape parameter was obtained purely by performing a large number of numerical experiments.It can be seen from Table 6 that when γ is comparatively high, i.e. 1 γ > , all solutions are in good agreement with the exact and close to one another.The error tends to get larger when the convective force becomes dominant, i.e. 1 γ < , and this is well-known as the problem is encountering the instability problem.
Nevertheless, it is interestingly to see that the variable shape proposed in this work produced a comparatively, slightly better quality than both Var-1, Var-3, and NAA*, IMQ*.

5.
In this investigation, three main objectives are aimed.Firstly, we have successfully applied the methodology of meshless method using radial basis function to PDE of convectiondiffusion type.Secondly, we have gathered variable shape parameters proposed in literature that are designed mostly for multiquadric type of RBF.We also performed numerical test on these shapes in conjunction with inverse multiquadric RBF.Last but not least, we have proposed in this work a new form of variable shape parameter that acts both linearly and exponentially, depending on the local distant norm.A large number of numerical experiments have been carried out in order to gather as much information as possible and some useful findings achieved in this work are as follows; • It is demonstrated that the RBF-collocation method can well be another alterative numerical tool for solving convection-diffusion type of PDEs.• All forms of shape parameters investigated here resulted in the same trend of condition number.
• Observed from all forms under investigation, the optimal values of shape parameter ε is seen to be between 0.9 -1.50 in most cases where the range of values of ε , i.e. min max ( , ) ε ε , is seen not to significantly influence the accuracy of the final results.
• The local node-density gives a notable effect where higher gradient in local problem parameters takes place.• The variable shape parameter proposed in this work is numerically proven to perform well when compared to those that are available in literature.A slightly better solution quality is also seen in the case where the problem becomes convective-dominated.
In the case where the problem is even more convective-dominated (i.e. 1 γ  ), nevertheless, the situation could be totally different, as encounter by many researchers.At this range it is well-known that the problem of instability can well ruin the whole simulation.From this, it is our further objective that some local feature such as the Peclet number should be taken into consideration when designing a new variable shape for inverse-multiquadric.This prompts the main objective of our future investigation.

Figure 1 :
Figure 1: Infinitely smooth inverse-multiquadric RBF; profile at different shape parameters denotes the set of collocation points, of boundary points.The centers jx used in equation (

.
The automatically self-adjusted parameter ζ is introduced and act as a weight function leading j ε to equal to the linear mode when 1 j = .

and 1 β
= in which γ is a given constant of diffusion coefficient.The boundary conditions are set with 0 u = on all four sides.The exact solution for this problem is given by;

Table 2 :
Variable shaper parameter schemes proposed in literature

Table 3 :
Solution quality comparison at different numbers of nodes (100 and 289), and the value of min

Table 4 :
Solution quality comparison at different numbers of nodes (100 and 289), and the value of min

Table 5 :
Solution quality comparison at different numbers of nodes (100 and 289), and the value of min