INTRODUCTION
The radial basis functions (RBFs), ϕ 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) = ϕ(r)∈R, the set of real numbers, with x∈R^{n} and r∈R or in other words, on the distance from a point of a given set {x_{j}} and ϕ(xx_{j}) = ϕ(r_{j})∈R where can normally defined as expressed in Eq. 1:
for some fixed points x∈R^{n} and r_{j} = xx_{j}_{2} is the Euclidean distance. Two classes of numerical tools that involve the use of RBFs are the socalled ‘meshfree method’ and ‘boundary element method’. There are many forms of RBFs invented, proposed, tested and documented nowadays under the context of finding solution to different type partial differential equations (PDEs). Below is the list of some wellknown and most widely used RBFs:
• 
Linear (LR): 1+r 
• 
Gaussian (GU): 
• 
Cubic (CU) : r^{3} 
• 
Polyharmonic (PY) in R^{3}: r^{2n1}, n∈ℕ 
• 
Polyharmonic (PY) in R^{2}: r^{2n} ln(r), n∈ℕ 
• 
Multiquadric (MQ): 
• 
Inverse multiquadric (IMQ): 
• 
Thinplate spline (TPS): r^{2} ln(r) 
• 
Matern/sobolev (MS): *[K_{V} (r) r^{V}] 
• 
*[K_{V} is an order Bessel function] 
In this study, nevertheless, the RBF type being focused on is the ‘Inversemultiquadric (IMQ)’ form expressed as follows:
where, β = ..., 3/2, 1/2, 1/2, 3/2,... and ε is the socalled ‘shape parameter’ and is normally given in an ‘ad hoc’ manner. Different values of ε, illustrated in Fig. 1, lead to different impact on the final results as well acknowledged. This kind of radial basis function has been proven to be strictly positive definite as nicely documented by Buhmann^{1}, meaning that the distance matrix of the interpolation problem is invertible. Throughout this study, it focused on the most popular form of the inverse multiquadric, i.e., β = 1/2, only.

Fig. 1: 
Infinitely smooth inversemultiquadric RBF measured numerically at different shape values, ε 
Appearing as an alternative numerical method, finite element, finite volume and finite difference, over the last two decades, the boundary element method (BEM) is now known to be another important tool for solving a wide range of applied sciences and engineering that involve linear as well as certain types of nonlinear partial differential equations (PDEs). Like other numerical schemes, the method is not without difficulties or challenges when it faces problems with nonlinear, transient, coupled and nonhomogeneous form of PDEs.
Brebbia and Butterfield^{2} proposed an improved version of the scheme and they named it as ‘dual reciprocity boundary element method (DRBEM)’. In the process of DRBEM, the solution is divided into two parts: Complementary solutions of its homogeneous form and the particular solutions of the inhomogeneous counterpart. Since the particular solutions are not always available especially in complex problems, the inhomogeneous term of the PDE is approximated by a series of simple functions and transformed to the boundary integrals employing particular solutions of considered problem. The most widely used approximating functions in DRBEM are radial basis functions (RBFs) for which particular solutions can be easily determined^{3}.
When it comes to using a radial basis function, particularly those forms with the parameter, it is then down to the user’s decision to choose what they believe is best and very often this is done under rather ‘ad hoc’ manner. For multiquadric (MQ) type of RBF, there are some investigations done on providing information regarding choosing optimal shape value. Hardy^{4} suggests that by fixing the shape at ε = 1/(0.815d), where, and d_{i} is the distance from the node to its nearest neighbor, good results should be anticipated. Also, in the work of Franke and Schaback^{5} where the choice of a fixed shape of the form where, D is the diameter of the smallest circle containing all data nodes, can also be a good alternative.
Some recent attempts to pinpoint the optimal value of ε involve the work of Zhang et al.^{6} where they demonstrated and concluded that the optimal shape parameter is problem dependent. In 2002, Wang and Lui^{7} 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.^{8} suggested that the final numerical solutions obtained are found to be less affected by the method when the approximation is applied locally rather than globally.
At this point, nevertheless, there are several facts to be noticed at. Firstly, most of the proposed forms are concerned with either numerical aspect or fixed values only, meaning that no real physical aspects included or being under consideration. Secondly, all the works were done under the context of pure RBFcollocation methods, nothing is done under DRBEM. Thirdly, the shape parameter contained in the inversemultiquadric type of RBF has not, by far, been taken a look at. This prompts our initiative idea of this study.
One of the classical nonlinear parabolic equation system widely acknowledged as the Burgers’ equations, named after the great Physicist^{9}. The equations retain the nonlinear of the governing equation in a number of applications, flow through a shock wave traveling in viscous fluid, the phenomena of turbulence, air/water pollution or chemical compound carried by fluids, sedimentation of two kinds of particles in fluid suspensions under the effect of gravity (Burgers^{9}, Nee and Duan^{10}).
For large values of the Reynolds number (Re), one of the major difficulties is due to inviscid boundary layers produced by the steepening effect of the nonlinear advection term in Burgers’ equations. This difficulty is also encountered in an inviscid Navier–Stokes equation for a convection dominated flow and in fact, the famous coupledBurgers’ equations are one of the principle equations used to evaluate any newly proposed numerical methods^{11}.
This all leads to our main objectives of this investigation which are, firstly, it proposed a new form of shape parameter contained in the inversemultiquadric RBF. Secondly, the proposed shape was designed to take into consideration the local physical of the problem at hand governed by the dimensionless Reynolds number (Re). Thirdly, the methodology of DRBEM was studied, applied and computationally implemented to one of the most challenging types of PDEs known as ‘Burgers’ equation’, famous for its rich in transient, couple and nonlinear phenomena. The results obtained from this investigation were validated against other alternative numerical works in literature when available.
MATERIALS AND METHODS
DRBEM for Burgers’ equations: Let us consider the system of twodimensional unsteady, nonlinear and coupled Burgers’ equations, expressed as follows:
Subject to the initial conditions:
u(x, y, 0) = β_{1}(x, y) and v(x, y, 0) = β_{2} (x, y)
and the boundary conditions:
u(x, y, t) = γ_{1}(x, y, t) and v(x, y, t) = γ_{2} (x, y, t)
for (x, y)∈∂Ω, t>0, u(x, y, t) and v(x, y, t) are the velocity components, β_{1}, β_{2}, γ_{1} and γ_{2} are known functions and Re is the Reynolds number, described by Keannakham et al.^{12}
The mathematical construction of the dual reciprocity boundary element method can start with the Poisson equation as follows:
Which as its equivalent integral form, given by Brebbia and Butterfield^{2} as:
where, u* the fundamental solution and the term is defined as , where n is the unit outward normal to Γ and can be written as follows:
Next, the method was applied with N and L being the number of boundary and internal nodes, respectively, b can be now approximated by:
Or, can be written as follows:
Or:
where, the function f is the radial basis function which is, in this work, the inversemultiquadric type. With this radial basis function, then the following form was obtained:
for some particular solution, Applying Green’s theorem, the boundary element approximation to Eq. 9, then it becomes, at a node ith:
For i = 1, 2,..., N.
After introducing the interpolation function and integrating over each boundary elements, the above Eq. 10 can be rewritten in terms of nodal values as:
where, the definition of the terms H_{ik} and G_{ik} are described by Toutip^{13}. The index k is used for the boundary nodes which are the field points. After application to all boundary nodes, using a collocation technique, Eq. 11 can be compactly expressed in matrix form as follows:
By substituting α = F^{–1}b from Eq. 8, into Eq. 12 making the right hand side of Eq. 12 a known vector. Therefore, it can be rewritten as:
where, Applying boundary conditions to Eq. 13, then it can be seen as the simple form as follows:
where, x contains N unknown boundary values of u’s and q’s.
After Eq. 14 is solved using standard techniques such as Gaussian elimination, the values at any internal node can be calculated from Eq. 11, i.e., c_{i} = 1 as expressed in Eq. 15 where each one involving a separate multiplication of known vectors and matrices:
Firstly, substituting Eq. 12 into Eq. 11 to get the equation system matrix which expressed as:
It is then possible to approximate the terms, and as below:
From Eq. 16 and Eq. 17, the following was obtained:
and:
where, and:
In the case, if were considered to be generated by using the same redial basis function, then:
For the time derivatives, the forward difference method was used to approximate time derivative and Now, substituting Eq. 3235, 38 and 39 in Eq. 43 and 44, then the following was obtained:
and
by setting and then the final forms of DRBEM for this type of equations A = πr^{2} as follows:
and
Note that the elements of matrices H, G and A depend only on geometrical data. Thus, they can all be computed once and stored.
Table 1: 
Some choices proposed in literature 


Proposed variable IMQ shape: Regarding the studies in the search for optimal choice of the shape parameter, many outstanding and wellknown forms proposed in the past are listed in Table 1. In this table, it can be seen that one of the first attempts in searching for variable form of shape parameter was done in 1990 by Kansa^{14 }which the experimental manner was looked at. This form was then modified by Kansa and Carlson^{15}. In their work, and were observed that the better quality of results can be achieved if these two values vary by several orders of magnitude. Sarra^{16 }and Sarra and Sturgill^{17 }later invented a linear form of variable shape parameter and applied to both interpolation and some benchmark partial differential equations. The command ‘rand’ is the MATLAB function that returns uniformly distributed pseudorandom numbers on the unit interval. In their study, the variable form of shape parameter was shown to be superior to the fixed form particularly when the system includes some information about the minimum distance of a center to its nearest neighbor h_{n}, with also a user input value μ. It is, nevertheless, to be noticed that all the attempts mentioned so far are under the context of the collocation meshfree method where there is no, by far, work carried out under the DRBEM. This leads to one of our main objectives of this investigation.
In this study, it proposed a new form of variable shape parameter where both linear and exponential manners are taken into consideration, expressed as in Eq. 2730:
Where:
and
and ζ is set to correspond to local Reynolds number (Re_{ij}) and is defined as follows:
and the dimensionless number called The Reynolds Number, (Re), which is the ratio of the inertia to the viscous or friction force and is defined, in a global sense, as:
where, L is the characteristic length scale, U is the velocity and v is the kinematic viscosity. In this investigation, this number is locally defined as:
where, v^{t} is defined in a local sense as:
where, N+L is the total number of nodes (both boundary and internal) involved in the system, the superscript t indicates the time level. The subscript i and j mean at ith and jth node, respectively.
Numerical setup: The computation process is done following the below algorithm:
• 
Initializing the matrix C by setting it to zero for both U and V 
• 
Setting all the values appearing in Eq. 25 and 26 and then applying the initial conditions to the right hand side and boundary condition to the left hand side, reducing to the form Ax = y 
• 
Using the iterative scheme as A^{(k)}x^{(k)} = y^{(k)}, where, k = 1, 2, 3.... The solutions of for both U and Vcan then be achieved when the algorithm stops i.e. when it reaches the following condition: 
with an user’s given value of a tolerance ω.
• 
The solutions obtained from step 3 will be used as the initial conditions in next timestep and the computing process continues until the time level given is reached 
All numerical simulations under this investigation were carried out on an Intel(R) Core(TM) i75500U CPU @ 2.4 GHz with 8.00 GB of RAM.
RESULTS
To demonstrate the effectiveness of the proposed shape parameter, the methodology was applied to two benchmark examples governed by the Burgers’ form of equations. The results obtained were validated against both the corresponding exact solutions and those published in literature. In some cases, the root mean square (RMS) error norm was used and it was defined as below:
It was noted here that the Reynolds number written as ‘Re’ from now on is the global Reynolds number described and used in Eq. 3133.
Experiment 1: The first example was the wellknown form where its exact solutions for validation were provided by using a HopfCole transformation nicely done by Fletcher^{18 }and was expressed as follows:
and:
where, and both the initial and boundary conditions to be imposed to the equation system was generated directly from the above exact form over the domain Ω = {(x, y):0<x<1, 0<y<1}.
Table 29 showed the computed values of both U and Vvalues for Re = 1 at t = 0.01 and Re = 10 at t = 0.5, respectively, obtained from this work compared against both its’ exact solutions and those provided by Biazar and Aminikhah^{19}. In this case, only 13 internal nodes were chosen and it can be seen that the method constructed in this study performed well and results agreed nicely with both references.
Table 2: 
Comparison of Uvelocity for Re = 1, t = 0.01 with boundary node N = 80 and internal nodes L = 13, for (ε_{min}, ε_{max}) = (0.1, 1) 


Table 3: 
Comparison of Uvelocity for Re = 1, t = 0.01 with boundary nodes N = 80 and internal nodes L = 13, for (ε_{min}, ε_{max}) = (0.1, 10) 


Table 4: 
Comparison of Vvelocity for Re = 1, t = 0.01 with boundary nodes N = 80 and internal nodes L = 13, for (ε_{min}, ε_{max}) = (0.1, 1) 


Table 5: 
Comparison of Vvelocity for Re = 1, t = 0.01 with boundary nodes N = 80 and internal nodes L = 13, for (ε_{min}, ε_{max}) = (0.1, 10) 


Table 6: 
Comparison of Uvelocity for Re = 10, t = 0.5 with boundary nodes N = 80 and internal nodes L = 13 for (ε_{min}, ε_{max}) = (0.1, 1) 


Table 7: 
Comparison of Uvelocity for Re = 10, t = 0.5 with boundary nodes N = 80 and internal nodes L = 13, for (ε_{min}, ε_{max}) = (0.1, 10) 


Table 8: 
Comparison of Vvalues for Re = 10, t = 0.5 with boundary nodes N = 80 and internal nodes L = 13, for (ε_{min}, ε_{max}) = (0.1, 1) 


Table 9: 
Comparison of Vvalues for Re = 10, t = 0.5 with boundary nodes N = 80 and internal nodes L = 13, for (ε_{min}, ε_{max}) = (0.1, 10) 



Fig. 2(ab): 
Comparison of Vvalues for (ε_{min}, ε_{max}) = (0.1, 10) at t = 0.1 with Re = 80 between the exact (left) and the computed ones (right) 
Solution profiles of both U and Vvelocity, at moderate Reynolds number, were illustrated in Fig. 24 showed the root mean square (RMS) error.
Experiment 2: For this example, it was adopted from the work of Aminikhah^{20} for solution validation. This problem is on the domain Ω = {(x, y): 0<x<1, 0<y<1)} with the following initial and boundary conditions:

Fig. 3(ab): 
Comparison of Uvalues for (ε_{min}, ε_{max}) = (0.1, 10) at t = 0.1 with Re = 80 between the exact (left) and the computed ones (right) 

Fig. 4(ab): 
Root mean square (RMS) error comparison at two intervals of (ε_{min}, ε_{max})and that obtained using a fixed value of 3, computed with boundary nodes N = 120 and internal nodes L = 24 at t = 1.0 
This leads to the exact solutions as the following forms:
Table 10: 
Comparison of Uvelocity for Re = 100, t = 0.5 with boundary nodes N = 80 and internal nodes L = 12, for (ε_{min}, ε_{max}) = (0.1, 1) 


Table 11: 
Comparison of Vvalues for Re = 100, t = 0.5 with boundary nodes N = 80 and internal nodes L = 12, for (ε_{min}, ε_{max}) = (0.1, 1) 


Table 12: 
Comparison of Uvelocity for Re = 500, t = 0 with boundary nodes N = 80 and internal nodes L = 12, for (ε_{min}, ε_{max}) = (0.1, 1) 

and
where, = exp (5π^{2} t/Re).

Fig. 5(ab): 
Comparison of Vvelocity for (ε_{min}, ε_{max}) = (0.1, 20) at t = 1 with Re = 1,000 between the exact (left) and the computed ones (right) 
Table 13: 
Comparison of Vvelocity for Re = 500, t = 0.5 with boundary nodes N = 80 and internal nodes L = 12, for (ε_{min}, ε_{max}) = (0.1, 1) 

Table 1013 provided results comparison amongst those produced in this study and both the exact solutions and other numerical investigation, Aminikhah^{20}. The solution profiles obtained at high Reynolds number were shown in Fig. 56 where Fig. 7 provided the comparison of RMS at different values of both Reynolds numbers and ranges of (ε_{min}, ε_{max}).

Fig. 6(ab): 
Comparison of Uvelocity for (ε_{min}, ε_{max}) = (0.1, 20) at t = 1 with Re = 1,000 between the exact (left) and the computed ones (right) 

Fig. 7(ab): 
Root mean square (RMS) error comparison at two intervals of (ε_{min}, ε_{max}) and that obtained using a fixed value of 1.5, computed with boundary nodes N = 120 and internal nodes L = 24 at t = 0.75 
DISCUSSION
Regarding the use of inverse multiquadric (IMQ) type of radial basis functions (RBFs), it is not surprising to realize that only a few works available in literature. Most of these works, moreover, were carried out under the context of meshfree/meshless method. Some recent works that adopted this type of RBF are Chantawara^{21} and Chuathong^{22}, where some fixed values of shape parameters were adopted. To validate results obtained in this paper, their suggested values of shape parameter were adopted.
The computational results produced in this work for the two benchmark examples are promising in several aspects. Firstly, by providing only two thresholds or bounding values of (ε_{min}, ε_{max}), searching for optimal choice becomes simpler. This feature can be seen more clearly when compared with those obtained by using a fixed value of ε, as shown in all Tables. The effect of the wide of the range (ε_{min}, ε_{max}) was seen not to have any significant effect on the final results quality ensuring the simplicity of the proposed shape form.
Another aspect that is desirable about the main findings of this study is the capability of the method to handle the instability of the problem when the Reynolds number increases. It is commonly found and is known to be a great challenge for any numerical procedure to tackle Burgers’ equations at high Reynolds number (>500), the flow becomes more turbulent. One of our previous works and many other works, Keannakham et al.^{12} and references herein, have shown that beyond this value of Reynolds number, the error norm is growing dramatically. It is shown clearly in this work that the proposed adaptive parameter is able to take into consideration the local physics that effects both the computing collocation matrix and the numerical algorithm itself.
CONCLUSION
In this study, the effectiveness of shape parameter contained in the inverse multiquadric type of RBF in conjunction with the method of DRBEM was investigated numerically. The investigation began with applying DRBEM to one of the most complicated PDEs namely Burgers’ equations. Then a new form of shape parameter that behaves locallyadaptive, i.e., it varies accordingly to the local change of the physics of the problem which, here, is the local Reynolds number(Re) was proposed. The proposed variable shape form also contains both linear and exponential aspects based on the distance between the center node i and pointed node j in the RBFcollocation numerical method adopted. Some important conclusions can now be withdrawn from the investigation and they are as follows:
• 
DRBEM has successfully been applied to Burgers’ equations using inverse multiquadric radial basis function at relatively high Reynolds number (Re) 
• 
It is found from all the results obtained in this work that the proposed shape parameter can outperform the fixed ones and certainly deserves further investigation 
• 
When compared with other numerical works, the accuracy lies in an acceptable level and moreover, gets better when the problem become advective dominated, at high Reynolds number 
SIGNIFICANCE STATEMENT
This study discovers the benefit of using the local physical aspect of the problem as the indicator for the radial basis function’s shape parameter which is designed to be locally adaptive. This study will help the researcher to uncover the critical area of choosing the optimal shape parameter when dealing with this particular form of equations. Thus, a new theory on these radial basis function’s adaptive shape parameters and possibly other related fields, may be arrived at.
ACKNOWLEDGMENTS
The corresponding author would like to express his sincere gratitude and appreciation to the Center of Excellence in Mathematics, Thailand for their financial support.