The dual reciprocity boundary element method ( DRBEM ) with multiquadric radial basis function for coupled burgers ’ equations

Krittidej Chanthawara1, Sayan Kaennakham2 and Wattana Toutip3 1Department of Mathematics, Faculty of Science, Khonkaen University, Khonkaen, Thailand, 40002. krittidej.kku@gmail.com 2School of Mathematics, Institute of Science, Suranaree University of Technology, Nakhonratchasima, Thailand, 30000. sayan_kk@sut.ac.th 3Department of Mathematics, Faculty of Science, Khonkaen University, Khonkaen, Thailand, 40002. wattou@kku.ac.th


INTRODUCTION
Appearing as an alternative numerical method over the last two decades, the boundary element method (BEM) has become an important tool for solving a wide range of applied science and engineering that involve linear as well as certain types of nonlinear partial differential (PDE).The main feature of BEM is that the PDE problem at hand is converted to an equivalent boundary integral Equation (BIE) using Green's theorem and a fundamental solution [1].As a result, the discretization process takes place only on the domain boundary rather than the entire domain, as usually required in other classical numerical schemes i.e.Finite Element Method(FEM), Finite Difference Method(FDM) or Finite Volume Method(FVM), leading to a much smaller system of Equation and effectively, a considerable reduction of data required in the computing process while satisfactory results accuracy can still be expected.
Despite its attractive aspect mentioned, nevertheless, there are some difficulties in extending the method to non-homogeneous, non-linear and time dependent problems.The main drawback encountered in these cases is the need to discretize the domain into a series of internal cells to deal with the terms taken to the boundary by application of the fundamental solution.This additional discretization then inevitably destroys some of the attraction of the method.To remedy this undesirable drawback, one of the attempts that is successful for the task is the so-called 'Dual Reciprocity Boundary Element Method (DRBEM)' which was first proposed by Nardini and Brebbia [2] in 1982.Ever since, the idea has extensively been developed by many researchers such as the perturbation DRBEM [3], the separation of variables DRBEM [4], and the Laplace transform DRBEM [5][6][7].
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 [8].
A radial basis function is a one variable and continuous function defined for r ≥ 0 that has been radicalized by composition with the Euclidean norm on d and can contain a free parameter or the shape parameter, normally denoted by e.If we are given a set of N centers, , a radial basis function then takes the form; ( Where and are coefficients to be determined by interpolation.For the this task, it is crucial for the final solution that an appropriate interpolation function f j is adopted and there are a number of investigations carried out on this important task.In applications related to BEM, in particular, as pointed out by Golberg and Chen [9] that the theory of RBFs provides a mathematical basis for the choice of trial function in the application of BEM, clearly indicating that the choice of interpolating functions is the key for the accuracy , an early work of Partridge and Brebbia [10] conceded that a simple linear form denoted as f = 1 + r is usually sufficient.This simple choice seems to be also widely adopted by many authors [8,[11][12][13],while later the computational work of Golberg [14] and Zhu et.al.[15] indicated that a combination of the thin plate splines (TPS), and some augmented linear terms 1,x, y in 2 (or 1, x, y, z in 3 ) could be more reliable.The comparatively higher accuracy, however, comes with the price as it has been demonstrated by Madych and Nelson [16], see also [17,18] that this type of RBF is non-differentiable in 3 and its convergence is very slow.Alongside with many attempts in this optimal RBF searching task, Golberg et al. [19] have shown that the so-called Multiquadric (MQ) interpolation, denoted as the form , with is highly effective when applied to the DRBEM.In spite of MQ's excellent performance, its accuracy is greatly affected by the choice of the shape parameter C. It is well-known that large values of parameter e can result in higher error in solution and decreasing e can improve the accuracy significantly [20].However, as pointed out by Schaback [21] , the direct way of computing this type of RBF interpolatant suffers from severe mathematically ill-conditioning.This then makes the search for optimal choice of e a very challenging topic [22], and remains as one of the main objectives of this paper.
The second objective of this study is concerned with applying the DRBEM to one of the classical non-linear parabolic equation system as the Burgers' Equation, named after the great Physicist Johannes Martinus Burgers' (1895-1981), as expressed in two dimensions as follows [23]; (2) Where Re being the Reynolds number, the ratio of inertial forces to viscous forces u(x,y,t), v(x,y,t) are the velocity components to be determined.For a domain with its boundary G, the above system is subject to the following initial conditions: And the boundary conditions: Here b 1, b 2 g 1 and g 2 are known functions.
The equation 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, sedimentation of two kinds of particles in fluid suspensions under the effect of gravity (see [24][25][26] and the references therein).With using the Hopf-Cole transformation, an analytic solution of the system is provided in 1983 by Fletcher [27].Ever since, the coupled twodimensional Burgers' Equation have been attracting and challenging a number of numerical investigations.Some of the recent ones include the application of a fully implicit finitedifference by Bahadir [23], that by using Adomain-Pade technique by Dehghan et.al.[28], that by adopting variational iteration method by Soliman [29], that by applying a meshfree technique by Ali et.al.[30], and also the work of Zhu et.al.[31](see also the references therein).Standing as one of the most challenging problem, this study aims to tackle by the method of DRBEM in conjunction with several shape parameters for multiquadric RBF in order to shed more light on optimal choice of the parameter when it comes to coupled 2-D Burgers' type of Eqns.

MATHEMATICAL BACKGROUND 2.1 THE DUAL RECIPROCITY BOUNDARY ELEMENT METHOD (DRBEM)
The way of formulating DRBEM is to consider the Poisson Equation; (4) Where b = b(x,y).The solution to above eqn(4) can be expressed as the sum of the solution of homogenous and a particular solution as; (5) Where u h is the solution of the homogeneous equation and is a particular solution of the Poisson equation such that .We approximate b as a linear combination of interpolation functions, as elaborated in Section 2.2 , for each of which we can find a particular solution at points which are situated in the domain and on its boundary.
If there are N boundary nodes and L internal nodes, there will be N+L interpolationfunctions, f j , and consequently N+L particular solutions, .The approximation of b over domain W is written, for an i th -node, in the following form.
Then we have a matrix equation for the unknown coefficients a j :j=1,2,3,…, N+L as follows; (7) This can easily be expanded as shown below.
The particular solution, , and the interpolation function, f j , i.e. the Multiquadric radial basis function used in this work, are linked through the relation.(8) Substituting eqn(8) into eqn(6) leads; (9) Therefore, from eqn(4) and with the expression of b in eqn(9), we then obtain; (10) Multiplying both sides with the fundamental solution, , and integrating over the domain yields; (11) Integrating by parts the Laplacian terms produces the following integral equation for each source note i th [32].(12) b x y f x y ( , ) ( , ) The dual reciprocity boundary element method (DRBEM) with multiquadric radial basis function for coupled burgers' equations The term in eqn ( 12) is defined as , where n is the unit outward normal to G, q u n ˆĵ j = ∂ ∂ q ˆj and can be written as follows; (13) Since eqn (12) contains no domain integrals, the next step is to write eqn (12) in discretised form, with summations over the boundary elements replacing the integrals.This gives, for a source node i th , the expression.
After introducing the interpolation function and integrating over each boundary elements, the above eqn( 14) can be re-written in terms of nodal values as; (15) Where the definition of the terms H ik and G ik can be found in [33].The index k is used for the boundary nodes which are the field points.After application to all boundary nodes, using a collocation technique, eqn(15) can be compactly expressed in matrix form as follows; (16)

Hu Gq
Hu Gq = ( ˆĵ If each of the vector and is considered to be one column of the matrices and respectively, then eqn( 16) may be written without the summation to produce; (17) By substituting from eqn (7), into eqn( 17) making the right hand side of eqn(17) a known vector.Therefore, it can be rewritten as; (18) Where .Applying boundary condition(s) to eqn (18), then it can be seen as the simple form as follow; (19) Where x contains N unknown boundary values of u's and q's.
After eqn( 19) is solved using standard techniques such as Gaussian elimination, the values at any internal node can be calculated from eqn (15), i.e. c i = 1 as expressed in eqn (20) where each one involving a separate multiplication of known vectors and matrices. (20)

THE MULTIQUADRIC RADIAL BASIS FUNCTION
Multiquadric RBF was proved by Hardy and Nelson [34] for the interpolation of topographical surfaces.It has very simple mathematical from as follows; (21) With multiquadric solution, eqn(8), as follows; Here the parameter C is known to play an important role in determining the accuracy of final solution and to be varying based on the problem at hand with the Euclidian distance, r, in 2 .
As pointed out by many researches [19,20,35,36] and nicely documented by Ding at al. [37] that when the value of shape parameter, C, is small, the resulting multiquadric radial basis function has a cone-like shape.When increasing its value, the peak of the cone gradually becomes more and more flat.Moreover, while the larger the shape parameter normally results in the smaller the approximation error, in practice, the multiquadric RBF approximation suffers from a trade-off principle.In other words, one can adjust the parameter C to improve the accuracy of multiquadric approximation with the price of getting more and more ill-conditioned coefficient matrix causing problems in the solving inverse matrices for eqn (19).The trade-off principle was initially observed in the multivariate interpolation process for scattered data; however, it also occurs in the numerical solution of PDEs.With this all complication and uncertainty on the most appropriate choice of C, it then remains one of our aims finding the ultimate interval value of C if there is any.

COMPUTER IMPLEMENTATION OF THE DRBEM FOR 2D COUPLED BURGERS' EQUASIONS 3.1 COMPUTING MATRICES CONSTRUCTION
We first substitute eqn (17) into eqn (16) to get the equation system matrix which expressed as; (23) Setting (24) Then eqn (23) becomes (25) Consider the following system of two-dimensional Burgers' Equation ( 26) (27) Subject to the initial conditions: (28) (29) where And the boundary conditions: where and D is a domain of the problem, D is its boundary; u(x, y, t) and v(x, y, t) are the velocity components to be determined, b 1 , b 2 , g 1 and g 2 are known functions and Re is the Reynolds number, which described in [31].
The analytic solution of eqn (26) and eqn (27) was given by Fletcher using the Hopf-Cole transformation [27].The numerical solutions of this equation system have been studied by several authors.Jain and Holla [28] developed two algorithms based on cubic spline function technique.Fletcher [29] has discussed the comparison of a number of different numerical approaches.
As described in [26], and are approximated similarly.
From eqn (23) and eqn (25) we obtain (36) and (37) where and In the case, if we consider to be generated by using the same redial basis function, then we get; (40) Then eqn (36) and eqn(37) becomes ( 43) and (44) For the time derivatives, we use the forward difference method to approximate time derivative and .The forward difference is expressed as; (45) where and Substituting eqn (32) to eqn (35), eqn (38) and eqn (39) in eqn(43) and eqn(44), then Then eqn(47) and eqn(48), respectively, becomes; (50) and ( 51) Substituting eqn(45) and eqn(46) in above eqns(50-51).The following expressions are obtained. ( Note that the elements of matrices H, G and A depend only on geometrical data.Thus, they can all be computed once and stored.

NUMERICAL ALGORITHM
In the computing step, we apply the initial conditions to the right hand side of eqn(52) and eqn(53), and boundary conditions to the left hand side.In each time step, the iteration is utilized to get the suitable solution as described below; 1. Constructing matrix C in eqn(49) by setting it to zero for the first iterative loop for both U and V. 2. Setting all the values appearing in eqn(52) and eqn(53), 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.3. Using the iterative scheme briefly expressed, running with k = 1, 2, … as A (k) x (k) = y (k) , we get the solutions of U and V.With a 'ad-hoc' value of e, the iteration process stops when; (54) 4. 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.

NUMERICAL EXAMPLES AND GENERAL DISCUSSION
Recalling eqn(2) and eqn(3), this section demonstrates the performance of the DRBEM by numerically solving three problems described as coupled Burgers' Equations in two dimensions.For solution validation, all computed results are compared against both the exact ones and other researches available with also the absolute error, defined as Abs Error = |Exact -Numerical|, is also clearly stated in the Tables of results for each case.Before going any further, however, it is important to point out two remarks resulted from our numerical experiments.Firstly, as previously mentioned in Section (2.2) the choice of the parameter C varies depending heavily on the nature of problems at hand [20,35,36,38].For this work, a large number of simulations were carried out using different C's values in the wide interval of [0,5] where we found that the most suitable values of C reside mostly within Secondly, for the study of boundary grid convergence, a number of simulations were performed and it was found as clearly depicted in Fig ( 1) that beyond 80 boundary nodes, the simulations tended to fall in to the same level of accuracy.As a result, for each example, we present only simulations obtained from cases with 80 boundary nodes.

( )
The first example is the well-known form where its exact solutions for validation are provided by using a Hopf-Cole transformation nicely done by Fletcher [27] and expressed as follows; (55) Where both the initial and boundary conditions to be imposed to the equation system are generated directly from the above exact form over the domain Table 1-4 show the computed values of both U-and V-values for Re=1 and at t=0.01 , obtained from this work with 4 values of C's; 1.1, 1.2, 1.4 and 1.5 ,compared against both its' exact solutions and those provided by Biazar and Aminikhah [39].In this case, only 13 internal nodes were chosen and it can be seen that the method constructed in this work performed well and results agree nicely with both references.
When increasing the Reynolds number up to 80 and with the same number of boundary nodes, ) 32 The dual reciprocity boundary element method (DRBEM) with multiquadric radial basis function for coupled burgers' equations 0.0E+00

EXAMPLE 2
This example deals with the 2D Burgers' Equations with the initial conditions and boundary conditions are obtained directly from the exact solutions which are available in the work of Biazar and Aminikhah [39], expressed as; The computational domain under the investigation is domain Table 5-6 show a satisfactory results computed from the algorithm constructed in this work when compared to another relative study of Hongqing and Meiyu [31].The tables provide the U -and V -values for Re = 1 and at t = 0.1 for different chosen values of parameter C. It has to be also noted that because most of the results obtained with using C = 1.1, 1.2, 1.3, 1.4 and 1.5 are very similar and therefore, only a few cases are shown here.Also, a total of internal nodes of 16 were selected for the computing process but only 8 of which are shown in the Tables where the referencing work is available.The following example gives a confirmation that DRBEM can well be applied to problems described by the coupled and 2D Burgers' equations.

EXAMPLE 3
This problem is on the domain with the following initial and boundary conditions: For this last example, we adopted the work of Aminikhah [40] for solution validation.It can be seen that acceptable result accuracies are shown in Table 7 -8 for computed solutions obtained at Re = 100, t = 0.5 with boundary nodes N = 80 and internal nodes L= 12 with the same range of C's.With an increasing in the Reynolds number, up to 500, Table 9 and Table 10 reveal that the approximated solutions produced from DRBEM are still in good agreement with both the exacts and those obtained in the reference.Moreover, at Re = 1000 the solution calculated still nicely lies in the same pattern with the exact as illustrated in Fig 5-6.

CONCLUSIONS
The methodology namely Dual Reciprocity Boundary Element was numerically studied to solve problems expressed as Burgers' equations in two dimensions.A computer code was implemented and used to reveal the performance quality of the method.A total of 3 problems u y t u y t u x t u x t t (0, , ) (1, , ) ( ,0, ) ( ,1, ) 0, 0   described as Burgers' equations in a coupled form were chosen to demonstrate how well the method can be applied to this type of problem.The second concentration was on finding the optimum choice of the parameter appeared in the radial basis function utilized at hand, the so-called Multiquadric, where it is well-known that neither too high nor too low value of this parameter can always lead to satisfactory results.A large number of simulations were carried out to cover all aspects possible and only some of which were presented here.Several important findings obtained from the work can be listed as follows; • Clearly demonstrated by all the examples, the implemented DRBEM can well be applied to problems expressed as the coupled Burgers' equations in two dimensions.• The best range of C's value is between 1.1 to 1.5 where we found that with C< 1.1, the final results tended to diverge from the exact ones whereas with C>1.5, the algorithm severely suffered from the ill-condition problem.• With high Reynolds number, the algorithm constructed is still able to produce reasonable good results.From these conclusions, therefore, it is remaining as our next challenge to further explore other radial basis functions for this kind of problem as well as to study different forms of MQ radial basis function.
Fig 2-3 clearly illustrate a good agreement in the prediction of values of U and V at t = 0.1 with C = 1.3 for both cases.

Figure 2 :Figure 3 :
Figure 2: Comparison of V-values for, at t = 0.1 with Re = 80, between the exact (left) and the computed ones (right)

C = 1 . 4 C
increasing the Reynolds number up to 80 and with the same number of boundary nodes, Fig 4 clearly illustrates another reasonable good agreement in the prediction of values of U and V at t = 0.1 with C = 1.3.

Figure 4 :Figure 5 :
Figure 4: Comparison of V-value, at t = 0.1 with R e = 1, between the exact (left) and the computed ones (right).

Figure 6 :
Figure 6: Comparison of V-value at t = 1 with R e = 1,000, between the exact (left) and the computed ones (right)

Table 1 :
Comparison of computational values of Figure 1: Comparison of Relative Errors for boundary grid convergence test cases, with a fixed number of internal nodes of 13, taken from Example 1

Table 2 :
Comparison of computational values V -vaules for

Table 3 :
Comparison of computational values of

Table 4 :
Comparison of computational values of V -vaules for

Table 5 :
Comparison of computational values of

Table 6 :
Comparison of computational values of