Variable Shape Parameter Strategy in Local Radial Basis Functions Collocation Method for Solving the 2 D Nonlinear Coupled Burgers ’ Equations

This study aimed at investigating a local radial basis function collocation method (LRBFCM) in the reproducing kernel Hilbert space. This method was, in fact, a meshless one which applied the local sub-clusters of domain nodes for the approximation of the arbitrary field. For time-dependent partial differential equations (PDEs), it would be changed to a system of ordinary differential equations (ODEs). Here, we intended to decrease the error through utilizing variable shape parameter (VSP) strategies. This method was an appropriate way to solve the two-dimensional nonlinear coupled Burgers’ equations comprised of Dirichlet and mixed boundary conditions. Numerical examples indicated that the variable shape parameter strategies were more efficient than constant ones for various values of the Reynolds number.


Introduction
Contrary to conventional numerical methods in solving the partial differential equations (PDEs), meshless methods [1], it is not essential to utilize meshes.Collocation methods are, in fact, meshless and easy to program.In addition, they allow some kinds of approaches to solve the PDEs.Considering the translation of kernels as trial functions, meshless collocation in asymmetric and symmetric forms is described in [2][3][4].It is highly successful, since the arising linear systems are easy to produce, leading to such good accuracy with the beneficial range of computational expenses.Furthermore, it has been newly proven [5] that the symmetric collocation [2,3] utilizing kernel basis is optimal along all linear PDE solvers using the same input data.By following this utilization of kernels, we can solve the PDEs.An overview of kernel methods prior to the year 2006 is presented in [6], while their recent variations are in [7][8][9][10][11][12] and the related references.
The selection of a shape parameter in a radial or kernel basis function has been a challengeable topic for some decades (see [13] and references over there).There are many experimental observations on the treatment of kernel-based methods under scaled (shape) parameter.In addition, there are optimization methods attempting to provide harmony between bad conditions and small errors.
There are some strategies in the literature to select shape parameters.Variable shape parameter (VSP) strategies can be compared with constant shape parameter (CSP) strategies.Many mathematicians utilize the CSPs in the radial basis function (RBF) approximations [14,15] due to their easy analysis in comparison with the VSPs; however, there are many findings from a broad applications set [16][17][18], indicating the advantages of using the VSPs.A number of strategies to choose a VSP would be investigated in this work.
Meshless kernel-based approaches are on the basis of a fixed spatial interpolation for time-dependent PDEs.By the means of method of lines, they can be converted to a system of the ordinary differential equations (ODEs) in time [19,20].
In general, the global collocation methods consider the whole domain.Although this method is simple to implement, the obtained collocation matrix is ill-conditioned, especially for large-scale problems.
Therefore, various localized meshless methods have been recommended in the literature to solve this problem (see [21,22] and the references therein).A main idea behind the local RBF collocation method is the utilization of the local sub-clusters of domain nodes (see Figure 1), called local domains of influence, with the local RBFs for the approximation of fields.In other words, for the approximation of function in any nodes of domain, we consider the local sub-clusters of domain nodes containing that node and specified number of nearest neighboring nodes of domain nodes.With the chosen influence domain, an approximation function is considered as a sum of weighted local RBFs with VSP.After that, the collocation approach is utilized to determine weights.In the following, all essential differential operators can be made via any operator on the approximation function.The most important benefit of utilizing the local method would be that the overlapping influence domain leads to many small matrices for every center node, rather than a large collocation matrix.As a result, sparse global derivative matrices would be achieved.Therefore, less computer storage and flops are required by this method.In principle, it is probable to use both uniform and random nodal arrangements in the method implementation owing to the approach meshless feature, but the accuracy wise efficiency of the uniform nodal arrangement is better than random nodes.One of the advantages of using uniform nodal arrangement is that small spatial derivatives and collocation matrices (of the size of sub-domain that is 5 × 5 in the current case) corresponding to every stencil require to be computed just once.This stores a considerable amount of the CPU time and also memory.In this work, by considering the 2D nonlinear coupled Burgers' equations, we try to show the validation of the proposed method.Indeed, we consider the Dirichlet and mixed boundary conditions.The accuracy, stability and efficiency are considered for the high Reynolds number, Re.We can find the applications and some numerical methods for this equation in [23][24][25][26][27][28][29].In addition, the numerical investigation of the two-dimensional coupled Burgers' equations can be found in [30].The 2D nonlinear coupled Burgers' equations are considered by a implicit finite-difference scheme in [31], the element-free characteristic method in [32], the variational multiscale element-free Galerkin method in [33], the Chebyshev pseudospectral method in [34], the global RBF method in [24], and the local RBF collocation method in [35].
The rest of the paper is organized as follows.In Section 2, a useful summary of the kernel-based trial functions is provided.In Section 3, the recommended local reproducing kernel method is given.In Section 4, the model is solved numerically and conclusions are provided in the last section.

Kernel-Based Trial Functions
Definition 1.Let Γ ⊂ R d be an arbitrary nonempty set.A function K : Γ × Γ → R is called (real) kernel on Γ. Definition 2. Let H be a real Hilbert space of function f : Ω → R. A function K : Ω × Ω → R is called reproducing kernel for H if [36] : ) holds for all x, y ∈ Γ.
for any finite set of points X = {x 1 , . . ., x n } ⊂ R d and any real numbers 1 , . . ., n .Furthermore, the function K is called positive definite on Γ if the quadratic form (1) is zero only for X = 0.
We consider a smooth symmetric positive definite kernel K : Γ × Γ → R on the spatial domain Γ.With each kernel, there is a reproducing "native" Hilbert space: of functions defined on Γ in the sense: where the inner product is related to the property of the kernel by: K(x, •), K(y, •) N K = K(x, y) for all x, y ∈ Γ.
For scattered nodes x 1 , . . ., x n ∈ R d , the translates K j (x) = K(x j , x) are the trial functions, and we intend to begin our work with them.Owing to the smoothness and explicitness of the kernel K, take derivatives with respect to both arguments.Therefore, we can achieve cheap derivatives of the K j .If the kernel is translation-invariant on R d , we have: One of the most important kernels with significant properties is the radial kernel [37].The radial kernel can be defined as: for a scalar function: The most important examples are the Whittle-Matern kernels , and K ν is the modified Bessel function of the second kind [38].

Constant Scaled Kernels
Kernels on R d can be scaled by a positive factor c by examining the new kernel: Moreover, the constant scaled translation-invariant and radial kernels on R d can be defined as: respectively.Large c increases the condition number of kernel matrices, and small c indicates sharp peaks, leading to approximating functions imperfectly.The selection of shape parameter c is regarded as a problem existing for over two decades (see [13,17,39] and the references therein), playing a crucial role in finding the numerical solution of the PDEs (see [40]).

Variable Shape Parameter
In many cases, it has been mentioned that VSP strategies would generate more valid results in comparison with CSP.A negative outcome of utilizing a variable shape is that the system matrix is not symmetric anymore.Up to now, the choice of the optimal value of the shape parameter remains an open question; no mathematical theory has been developed up to now to determine its optimal value.
If the kernel K is radial, i.e., K(x, y) = φ r 2 , then the variably scaled radial kernels on R d can be defined as: where c j is the shape parameter corresponding to the jth center.For example, the Whittle-Matern kernel with VSP is given by: In this paper, we have proposed some strategies to choose the VSP: strategy 1 : strategy 5 : Note that, c min and c max are positive parameters and denote the minimum and maximum of c j , respectively.
In [41,42], the linear VSP is in the following form : Here, we modified the linear shape parameter formula, and considered c j as a new VSP by strategies 1 and 2. Strategies 3 and 5 are known as exponential and trigonometric VSPs, respectively [41,43].
In this work, new VSP c j is proposed in strategy 4.
The approximation function via radial kernel can be written as follows: for a set of distinct nodal points {x 1 , x 2 , . . ., Then, the unknown parameter λ j , which is independent of r j , can be found by putting each node x i , i = 1, 2, . . ., n into Equation ( 2) and solving the following system of linear algebraic equations: AΛ = F, where:

Numerical Method
In this section, the proposed local reproducing kernel method is introduced.Considering the time-dependent PDE: with the initial condition: and Dirichlet or Neumann boundary conditions: Suppose L is a differential operator, I is a linear operator, u ∈ H, f ∈ F, where H and F are Hilbert spaces of functions on Γ and we assume that the problems (3)-( 6) are well-posed.Let us consider the discretization points x i , 1 ≤ i ≤ n, and a symmetric positive definite kernel K : Γ × Γ → R. Let us reorder the points successively into points {X I ∪ X D ∪ X N }, where X I = {x 1 , . . ., x z 1 } is the set of interior points, X D = {x z 1 +1 , . . .x z 1 +z 2 } is the set of Dirichlet boundary points, and k=1 , which contains the center x i and its n s − 1 nearest neighboring points and forms the radial basis n s corresponding to these points.To approximate the solution u(x, t) over S i , we use a linear combination as follows: In matrix form, we have: which gives: where: 2 , t , . . ., u [i] x The derivative of the approximate solution can also be approximated at the center locations by applying a linear differential operator L to the local interpolation (7).Depending on the problem at hand, L will be either a single derivative operator or a linear combination of derivative operators.The equation: x ∈ S i , gives: where: . Now, we write the PDE (3) at a point x i , i = 1, . . ., z 1 as follows: Lu(x i , t) = f (x i , t), then: Thus: Now, let: By the Dirichlet boundary (5), we have: and by the Neumann boundary (6): Hence: in which: Now, let I i be a vector that contains the indices of center x i and its n s − 1 nearest neighboring points.
We consider the (n − z 1 − z 2 ) × n sparse matrix W as follows: Therefore, Equation ( 11) leads to: The unknown vector U N will be considered in terms of the unknown vector U I by solving the following equations: Furthermore, the initial condition (4) leads to: Let D be an z 1 × n sparse matrix as follows: where: .
Then, we have: By considering Equations ( 8), ( 10) and ( 12), Equation ( 9) leads to the following ODEs: with the initial conditions (13), where L is a linear differential operator and F is an operator, which may have some global sparse matrices from local contribution.

Validation of the Method
Let us consider the following system: where X = (x, y) ∈ Γ, Re is the Reynolds number, V 0 , U 0 , f D , g D , f N , and g N are known functions, Γ ⊂ R 2 is the domain set, ∂Γ is the boundary of the domain set Γ, ∆ is the Laplace operator, U and V are unknown functions.We consider the set of points {X I ∪ X D ∪ X N } for discretization equations, where X I = {X 1 , . . ., X z 1 } is the set of interior points, X D = {X z 1 +1 , . . ., X z 1 +z 2 } is the set of Dirichlet boundary points, and X N = {X z 1 +z 2 +1 , . . ., X n } is the set of Neumann boundary points, z 1 and z 2 are the number of interior and Dirichlet boundary points.For each X i (i = 1, . . ., n), we consider a stencil S i = {X [i] k } n s k=1 , which contains the center X i and its n s − 1 nearest neighboring points and form the radial basis n s corresponding to these points.Let us consider the symmetric positive definite kernel K : Γ × Γ → R, and try to find the functions U(•, t) | [0, T] → N K , and V(•, t) | [0, T] → N K , where: Here, N j (X) are the RBFs corresponding to the kernel K, which is reproducing in the native Hilbert space N K .Thus, we have: where: 1≤k≤n s ,1≤j≤n s , , We now write the PDE ( 14) at the interior points X i (i = 1, . . ., n i ) as follows: With the aid of Label ( 17), we have: where: Let I i be a vector that contains the indices of center X i and its n s − 1 nearest neighboring points.We consider the z 1 × n sparse matrices D 1 , D 2 and D 3 as follows: Then, Equation ( 14) leads to: where: and .* denotes the pointwise product between two matrices or vectors.The Dirichlet boundary conditions imply that: With the Neumann boundary conditions, we have: where: and: . Suppose: hence the vector U N and V N are as: By substituting Labels ( 22) and ( 23) in (21), we obtain the system of the ODEs with the initial conditions: where:

Numerical Results
The results of our scheme for the numerical solution of two problems with Dirichlet and mixed boundary conditions have been presented.We take the Matern kernel due to strong convergence rate with the RBF parameter ν = m − d/2 = 2, and RBF scale c, i.e., we work with the kernel: which is reproducing in the Hilbert space W 3 2 (R 2 ).In addition, we consider c min = 1 and c max = 10 in variable shape strategies for all experiments in examples.
In the implementation of this technique, we have also used Legendre points: belonging to the interval [−1, 1] that can be easily transferred to the interval [a, b] by the transformation y = b−a 2 x + a+b 2 , and uniform points: The ODE solver ode113 of MATLAB (R2014b, MathWorks, Natick, MA, USA) is used to solve the final ODE system (21).The accuracy of the numerical results, is measured by the maximum absolute error defined as : where u and ũ represent the exact and approximate solutions, respectively.Example 1.Take problems ( 14)- (16).The initial and Dirichlet boundary conditions can be chosen by considering the following solutions, which are stated in [35]: .
We will compare the obtained solution in our proposed method with the results described in [11,35] because there is not the exact solution.Fake oscillations have been observed by using finite element method (FEM), finite-difference method (FDM), element free Galerkin method [32] and Galerkin-reproducing kernel method [11].In [35], an adaptive upwind technique has been innovated to avoid wiggles in the recommended local RBF collocation method.As mentioned in [11], even a very fine grid cannot get rid of the oscillatory behavior caused by a sharp gradient.To omit the wiggles and defeat instabilities for Re = 1000, we have used RBFs with CSP as shown in Figure 18.
The solution is smooth and near the front we do not have any instability for small Re and, therefore, no VSP is required.In our case, for Re = 10, 000, Figure 19 shows the oscillatory behavior.To avoid these wiggles, we have used VSP (strategy 2) in order to stabilize the solution near the sharp front.

Conclusions
Six shape parameter strategies were compared in this study.One strategy used a constant shape, while the other five used a different value of the shape parameter at each center.Results show that VSPs can improve the condition number and the solution accuracy.The approach is successfully applied to solve the 2D nonlinear coupled Burgers' equations with Dirichlet boundary conditions for high Re.Owing to Neumann boundary conditions, instabilities appear near the sharp gradient without any special filtering technique for Re = 1000, 10,000 for mixed boundary conditions.These fake oscillations were also observed using the FDM, FEM and element free Galerkin methods.In [35], an adaptive upwind technique was devised to avoid these wiggles for Re = 1000.To get rid of the oscillatory behavior, we proposed constant and VSP for Re = 1000, 10,000, respectively.

Figure 1 .
Figure 1.The uniform node arrangement and the schematics of the local domains of influence in the interior, near boundary and corner points using ns = 5, 1D (left), and 2D (right), where × denotes the center node in 1D.
show absolute error distributions at time T = 2 for Re = 10 with strategy 3 for different number of points.These figures show, if we increase nodes, absolute error decreases.The absolute error distributions at time T = 2 for Re=100 with CSP are shown in Figures 5-7 with a different number of points showing the absolute error distributions at time T = 2 for Re=100.It can be seen that the error decreases with an increasing the number of points.Figures 8-12 present absolute error distributions at time T = 2 for Re=100 with different shape strategies.Figures 13 presents absolute error distributions at time T = 2 for Re=1000 by CSP.Figures 14-17 show absolute error distributions at time T = 2 for Re=1000 with different shape strategies.

Figure 18 .Figure 19 .
Figure 18.Numerical results at time T = 0.4, x = 0.5, with Re = 1000, c = 10, n = 441, on the left U and on the right V. (test problem 2) x n } is the set of Neumann boundary points, z 1 and z 2 are the number of interior and Dirichlet boundary points.For each x i

Table 1 .
Numerical results for U at time T = 2 with Re = 100, n = 441.

Table 2 .
Numerical results for V at time T = 2 with Re = 100, n = 441.