The Lie-group method based on radial basis functions for solving nonlinear high dimensional generalized Benjamin-Bona-Mahony-Burgers equation in arbitrary domains

The aim of this paper is to introduce a new numerical method for solving the nonlinear generalized Benjamin-Bona-Mahony-Burgers (GBBMB) equation. This method is combination of group preserving scheme (GPS) with radial basis functions (RBFs), which takes advantage of two powerful methods, one as geometric numerical integration method and the other meshless method. Thus, we introduce this method as the Lie-group method based on radial basis functions (LG-RBFs). In this method, we use Kansas approach to approximate the spatial derivatives and then we apply GPS method to approximate first-order time derivative. One of the important advantages of the developed method is that it can be applied to problems on arbitrary geometry with high dimensions. To demonstrate this point, we solve nonlinear GBBMB equation on various geometric domains in one, two and three dimension spaces. The results of numerical experiments are compared with analytical solutions and the method presented in [6] to confirm the accuracy and efficiency of the presented method. Preprint submitted to Applied Mathematics and Computation 1 November 2017


Introduction
In mathematics and physics, nonlinear partial differential equations are partial differential equations with nonlinear terms. As said in [1,2], many phenomena in engineering and applied sciences are modeled by nonlinear evolution equations.They describe many different physical systems, ranging from gravitation to fluid dynamics, and have been used in mathematics to solve problems such as the Poincaré conjecture and the Calabi conjecture. They are hard to study: there are almost no generic methods that work for all such equations, and usually each particular equation has to be studied as a separate problem. The science of exact form solutions of the nonlinear partial differential equations helps the experiment of numerical solvers, aids in the stability analysis of solutions and guides to a better conception of nonlinear phenomena that these equations model [3]. But finding an exact solution for nonlinear PDEs is difficult. Therefore, many researchers are studying numerical methods for solving the nonlinear PDEs.
One particular class of such numerical methods that have been highly regarded in recent years is formed by geometric numerical integration methods [17,18]. The concept of a Lie group is very helpful in creating some trusty and potent numerical approaches to integrate ODEs, whilst preserving their invariant features. Indeed, preserving the Lie group structure under discretization plays an serious role in the improvement of a appropriate behavior in the error minimization [19]. So, by subscription the geometric structure and invariance of the original ODEs, new techniques can be devised, which are more strong, accurate and stable than the formal numerical approaches [16]. From the class of numerical Lie group methods, the so-called group preserving scheme (GPS), introduced by Liu [20], is a numerical approach which utilizes the Cayley transformation and the Padé approximations in the Minkowski space. Avoiding the spurious solutions and ghost fixed points are main benefits of this method [21][22][23][24][25][26].
In recent years, different techniques have been employed to apply the GPS for solving partial differential equations [27][28][29][30][31][32]. In most of these methods, using a technique known as the method of lines, all variables in the problem except one, are discrete and a system of ordinary differential equations is produced. Then, the GPS is applied for solving this system of ordinary differential equations. But some problems may arise in this process. The method is very sensitive in selecting the size of discrete step, this means that many criteria including stability, consistency and etc. should be considered. So, if the size of the discrete step selected is very small or large the consistency and stability of the method can be problematic and the solution will not be reliable. On the other hand, if for obtaining good results we select the length of discrete step too small, this will result in a large number of ordinary differential equations. Thus the method requires a lot of calculations. With these explanations, implementation of GPS on complex non-linear partial differential equations such as the nonlinear GBBMB equation is virtually impossible. For dealing with these difficulties, in this paper we propose a technique that combines the GPS with radial basis functions. We have named this scheme the Lie-group method based on radial basis functions.

A brief review of the GBBMB equation
In this paper, we consider the nonlinear GBBMB equation of the following form with initial condition and boundary condition where F (u) is a nonlinear function, f is a function of space and time variables, Ω ⊆ R n , ∆ and ∇ are the Laplacian and gradient operators, respectively. By choosing F (u) = 1 2 u 2 , Eq. (1) is called the Benjamin-Bona-Mahony-Burgers (BBMB) equation. If remove the term ∆u, we have the generalized Benjamin-Bona-Mahony (GBBM) equation. If we remove the term ∆u and also put F (u) = 1 2 u 2 , Eq. (1) is called the Benjamin-Bona-Mahony (BBM) equation.
Among the applications of BBM equation we can mention the following [4][5][6]: • Analysis of the surface waves of long wavelength in liquids; • hydromagnetic waves in cold plasma; • acoustic-gravity waves in compressible fluids; • acoustic waves in harmonic crystals.
Different analytical and numerical methods have been used for solving the BBMB equation in the literature. Here, we review some of them. The authors in [4] use the first integral method for getting an analytic solution of the modified BBM equation.
The BBM equation is solved using a generalization of the well-known tanh-coth method by Gómez and co-authors [7]. The nonlinear stability of planar viscous shock profiles of the GBBMB equation in two dimensions is analyzed in [8]. Guo [9] and his co-workers studied the optimal decay rates of solutions for the GBBM equation in multi-dimensional space n ≥ 3. Authors of [10] showed an exponential decay rate of the solutions to the Cauchy problem of the GBBMB equation by employing the space-time weighted energy method. The BBM and modified BBM equations are solved using the exp-function method to obtain some new soliton solutions by Noor and his co-workers [11]. Authors of [12] studied the convergence rate of the global solutions of the GBBMB equation to the corresponding degenerate boundary layer solutions in the half-space. A Crank-Nicolson-type finite difference scheme with the proof of the second-order convergence in the discrete H 1 -norm is discussed in [13] for solving the BBM equation. Omrani [14] presented a technique based on the standard Galerkin method and Crank-Nicolson formula for space and time variables, respectively, and he proved the convergence of a linearized Galerkin modified scheme. The interested readers are referred to [15,5,6].
Liu [20] has derived a Lie group transformation for the augmented dynamical system on the future cone, and developed the group preserving scheme for an effective numerical solution of nonlinear ODEs. Consider a system of n ordinary differential equations: where u(x) is an n-dimensional vector, x is an independent variable and f is a vector-valued function of u and x (the vector field). We can transform Eq. (4) into the following (n + 1)-dimensional augmented dynamical system: Here we assume ∥u∥ > 0, and hence the above system is well-defined.
Obviously, the first equation in Eq. where is a Minkowski metric, I n is the identity matrix, and the superscript T stands for the transpose. In terms of (u, ∥u∥), Eq. (6) becomes where the dot between two n-dimensional vectors denotes the Euclidean inner product. The cone condition is thus the most natural constraint that we can impose on the dynamical system (5).
Consequently, we have a n + 1-dimensional augmented differential system, with a constraint (6), where satisfying is an element of the Lie algebra so(n, 1) of the proper orthochronous Lorentz group SO 0 (n, 1). Notice that, although not stated explicitly, the matrix A depends indeed on U through u and ∥u∥.
Once the problem is formulated in SO 0 (n, 1), it is crucial that any numerical scheme used to solve (9) preserve the Lie group structure. Otherwise, the cone condition (6) is not satisfied and the numerical approximation is no longer consistent. This can be done if one is able to generate approximations U k by where U k denotes the numerical value of U at the discrete x k , and G(k) ∈ SO 0 (n, 1) verifies where G 0 0 > 0 is the 00th component of G. A first-order approximation is obtained by taking propose instead the Lie-group version of the simple forward Euler method, namely Since A(k) ∈ so(n, 1) and U k ∈ SO 0 (n, 1), then U k+1 ∈ SO 0 (n, 1) by construction.
The exponential mapping exp[∆xA(k)] admits the closed-form representation where and For notational convenience, we have used f k = f(u k , x k ). Substituting the above exp[∆xA(k)] for G(k) into Eq. (12), we obtain This scheme therefore preserves group properties for all ∆x > 0. In the practical numerical calculation, we only need Eq. (21) such that by knowing the initial value u 0 can obtain u k , k = 1, 2, 3, · · · . Theorem 1. The cone condition is preserved by Eqs. (21) and (22) for every ∆x.
Proof. As a matter of fact, this is a trivial consequence of the fact that U T k gU k = 0 and (12), with G(k) = exp[∆xA(k)]. Alternatively, one can show that as follows. From Eq. (21), we have , (24) after doing the above inner product and possible simplifications, we obtain Therefore, according to definition α k and β k from Eqs. (19) and (20) we obtain where the last equal is obtaied from Eq. (22). Thus, the proof is finished.
Theorem 2. Scheme (21) unconditionally preserves the fixed points of the differential equation (4) and their stability properties.
Therefore, it is obvious that This means that u k is a fixed point of the discretized mapping (21) if and only if the point u k is an equilibrium (critical, fixed) point of the system (4). We next check the property of the fixed point. The Jacobian of the map (21) is At the fixed point f k = 0, we thus have Recalling that η k > 0, the property of the fixed point is not altered by the map (21). More precisely, the map (21) has the same type of stability that the system (4) has.
According to the preceding theorems, the long-term behavior of the original system can be described very well by this numerical scheme.

1-dimensional case
Clearly, Ω is now an arbitrary interval in R. We consider the following expansion of u(x, t) [33][34][35][36][37][38][39][40][41]: in which where r j = {∥x − x j ∥, x j ∈ Ω} denotes the Euclidean distance between x and x j ({x j } M j=1 are usually called center nodes). As mentioned in [42] the accuracy of many schemes for interpolating scattered data with radial basis functions depends on a shape parameter c of the radial basis function. Most authors use the trial and error method for obtaining a good shape parameter that results in best accuracy. In this paper we also use from this method.
In Eq. (31) we have M + 2 nunknown variables. Besides the M equations resulting , two additional equations are required that add two condition as follows: where c n j ≡ c j (t n ). For the use of Kansas method, we let i=2 are interior points and x 1 and x M are boundary points.
By substituting Eq. (31) for interior points in Eq. (1) we obtain where ∇φ(r ij ) = {φ ′ (r j )} | x=x i and ∆φ(r ij ) = {φ ′′ (r j )} | x=x i . Now by defining the following vectors Eq. (35) is converted in u ′ (t) = f(t, u(t)), i.e., an equation of the form (4). Now, we use the GPS (21) and solve Eq. (35) respect to t. We define where τ = T /N is the step size of time variable and T is the final time. By using Eq. (21) for n = 0 we get the following form: and η 0 is calculated from Eq.
, and η n is calculated from Eq. (21). Using Eqs. (33), (38), (39) and condition boundary (3) we have the following matrix form: where also we have in which and also we obtain and finally we can write b n+1 After solving the algebraic system of equations Ac n+1 = B n+1 at each time step, we can construct the solution using Eq. (31). We use the LU decomposition method for solving linear algebraic system of equations Ac n+1 = B n+1 .

n-dimensional case (n ≥ 2)
Now Ω is an arbitrary interval in R n , n ≥ 2. For x = (x 1 , x 2 , ..., x n ) ∈ R n , we construct an approximate expansion of u(x, t) of the form in which where r j = {∥x − x j ∥, x j ∈ Ω} denotes the Euclidean distance between x and x j .
For the use of Kansas method, we let By substituting Eq. (46) for interior points in Eq. (1) we obtain where ∇φ(r ij ) = {∇φ(r j )} | x=x i and ∆φ(r ij ) = {∆φ(r j )} | x=x i . Now by definition the following vectors Eq. (49) has the form of Eq. (4). We use the group preserving scheme (21) and solve Eq. (49) respect to t. We define t k = nτ, n = 0, 1, 2, ..., N, where τ = T /N is the step size of time variable and T is the final time. For n = 0 we get the following form: and η 0 is calculated from Eq. (21). where , and η n is calculated from Eq. (21). Using Eqs. (52), (53) and condition boundary (3) we have the following matrix form: in which is M × M matrix, the operator Υ is where and also we obtain and finally we can write After solving the algebraic system of equations Ac n+1 = B n+1 at each time step with LU decomposition method, we can construct the solution using Eq. (46).
In this section we present the numerical results obtained with the proposed method on four test problems. We test the accuracy and stability of the method described in this paper by applying it with different values of h (distance points in space) and τ (length of time step). We performed our computations using Mathematica 10 software on a Pentium V, 2800 MHz CPU machine with 2 GB of memory. In particular we compute the following error norms: where RMS means root-mean-square, T is the final time, M is the total number of points and N is the total number of time discretization.

Test problem 1
We consider the 1-dimensional nonlinear GBBMB equation and boundary conditions

Test Problem 2
Next we consider the two dimensional nonlinear GBBMB equation with initial condition and boundary conditions with f (x, y, t) = sin(x + y) (

+ 2t − 2t 2 cos(x + y)
) . (67) Then the exact solution of this test problem is u(x, y, t) = 1 + t sin(x + y). We  Fig. 3. Also, we solve this problem on irregular domain that its criterion is r = 1 n 2 [1 + 2n + n 2 − (n + 1)cos(nθ)]. We plot this irregular domains for n = 4 and n = 8 in Fig. 4 which are marked with Ω 4 and Ω 5 , respectively. Tables 3 and 5 show errors obtained for the present method with h = 1/10, T = 5 on rectangular domain and with h = 1/10, T = 10 on triangular domain for different value of τ , respectively. From this tables we can see the method is convergent when τ is smaller. In Table 4 we compare the our results with method of [6] for final times T = 5 and T = 10. We can see in two cases the results of our method is better. In

Test Problem 3
For our next problem we take again a 2-dimensional nonlinear GBBMB equation, ∂u(x, y, t) ∂x + ∂u(x, y, t) ∂y = cos(u(x, y, t)) ∂u(x, y, t) ∂x + cos(u(x, y, t)) ∂u(x, y, t) ∂y but now with initial condition u(x, y, 0) = sech 2 (x + y − r) + sech 2 (x + y + r), x, y ∈ Ω, and boundary conditions We have chosen f (x, y, t) so that the exact solution of this test problem is u(x, y, t) = exp(t) ( sech 2 (x + y − r) + sech 2 (x + y + r) ) . We solve this problem with the new method presented in this article with several values of h, τ and r at different final time on the rectangular, triangular, and circular domains. Also, we solve this problem on irregular domains that their criterion is Ω 6 = { 1 n 2 (1 + 2n + n 2 − (n + 1)cos(nθ)) , (n = 16)} and Ω 7 = {exp(sin(θ))sin 2 (2θ) + exp(cos(θ))cos 2 (2θ)} which we plot this irregular domains in Fig 4. To show the accuracy of the method from Table 6 we can see the comparison between the present method and method of [6]. From Figs.
11-18 we can show the numerical results for this problem on different domains for various r.

Test Problem 4
In the last example, we consider the three-dimensional nonlinear GBBMB equation with initial condition u(x, y, z, 0) = 1, x, y, z ∈ Ω, and boundary conditions with f (x, y, z, t) = (4 + 3t)sin(x + y + z) Then the exact solution of this test problem is u(x, y, t) = 1 + tsin(x + y + z). We solve this test problem with LG-RBFs in the current article with several values of h and τ on the cubic and spherical domains that are shown in Fig 19. From Table   7 we can see the numerical results of our method is more accurate than method of [6] at different points. From Tables 8, 9 we see the errors obtained for this example in different times with h = 1/10 and τ = 1/100 on cubic and spherical domains, respectively. Also, we fixed z = 0 and T = 1 and plot the approximation solution and absolute error using the present method for this example on the cubic and spherical domains in Figs. 20 and 21, respectively.

Conclusions
The nonlinear GBBMB is an equation which is widely used in various sciences.
In this paper, we present a technique for solving this equation in one, two and three dimensions and different types of domains by combining a first-order Liegroup method with radial basis functions. Numerical results show the efficiency and accuracy of the new technique and its feasibility for solving different types of partial differential equations, in particular on domains with complicated geometry and with high dimensions. Extensions to higher orders in time will be the subject of a subsequent study.  Table 2 Comparison between errors of present method and method of [6] for example 1 with      Table 4 Comparison between errors of present method and method of [6] Table 6 Comparison between errors of present method and method of [6] for example 3 with r = 0, Method of [6] Present method Method of [6] Present method         The boundary points are red and interior points are blue. Table 7 Comparison between errors of present method and method of [6] for example 4 with h = 1/10 and T = 1 at different points of domain on cubic domain.