Solitary Wave Solutions of the Generalized Rosenau-KdV-RLW Equation

: This paper investigates the solitary wave solutions of the generalized Rosenau–Korteweg-de Vries-regularized-long wave equation. This model is obtained by coupling the Rosenau–Korteweg-de Vries and Rosenau-regularized-long wave equations. The solution of the equation is approximated by a local meshless technique called radial basis function (RBF) and the ﬁnite-difference (FD) method. The association of the two techniques leads to a meshless algorithm that does not requires the linearization of the nonlinear terms. First, the partial differential equation is transformed into a system of ordinary differential equations (ODEs) using radial kernels. Then, the ODE system is solved by means of an ODE solver of higher-order. It is shown that the proposed method is stable. In order to illustrate the validity and the efﬁciency of the technique, ﬁve problems are tested and the results compared with those provided by other schemes.

Kaya and Aassila calculated the explicit solutions of the KdV equation with an initial condition by using the Adomian decomposition method [13]. Özer and Kutluay applied an analytical-numerical method to the KdV equation [14]. The RLW equation was developed by Peregrine, as an alternative to the classical KdV formulation in order to investigate the behavior of the solution [15,16]. Benjamin et al. proved the existence and uniqueness of the solution of the RLW model and determined its exact expression subject to restrictions in the initial and boundary conditions [2]. The RLW is also adopted in the modeling of long waves with small amplitudes on the water surface [17]. A noteworthy feature of the RLW problem is that the collision between two solitary waves results either in sinusoidal solutions, or in secondary solitary waves [18]. Since the KdV cannot describe wave-all and wave-wave interactions, another model, known as the Rosenau equation, was proposed by Rosenau to describe the dynamics behavior of dense discrete systems [7]. Zuo studied the solitary wave and periodic with the initial condition u(x, 0) = f (x) (2) and boundary conditions u(a, t) = u(b, t) = 0, u x (a, t) = u x (b, t) = 0, u(a, t) = u(b, t) = 0, u xx (a, t) = u xx (b, t) = 0, (3) where u = u(x, t), is a real-valued function, the real constants α, β, γ and µ are non-negative, p ≥ 2 is a positive integer, and f (x) is a given smooth function.

Lemma 1.
(See [25].) The following conservative properties for the initial value problem (1) hold and where Q(0) and E(0) are constants depending on the initial conditions.
In recent years, various analytical and numerical methods have been used to approximate the solution of the initial boundary value problem (1)-(3). Razborova et al. presented a theoretical approach based on the Ansatz method for the Rosenau-KdV-RLW equation [9]. Later, Razborova et al. used a semi-inverse Variational Principle to retrieve a single solitary wave solution [22]. Additionally, Razborova et al. and Sanchez et al. discussed the solutions of the perturbed Rosenau-KdV-RLW equation [23,24]. Wongsaijai et al. constructed a three-level weighted average implicit finite difference (FD) technique [19]. Pan et al. presented a C-N pseudo-compact conservative numerical scheme based on the FD technique [20]. Fernández and Ramos investigated a three-point compact method with fourth-second accuracy [21]. Wang et al. and Hu and Wang formulated FD schemes with linear three-level [31] and high-accuracy conservative [33] [42].
In this paper, we use the local meshless radial basis function (RBF) for solving the general Rosenau-KdV and the Rosenau-RLW models. Section 2 formulates and discusses the local meshless RBF based on the finite difference (RBF-FD) technique for discretizing Equation (1). Section 3 provides five numerical examples and compares the results with those of other schemes proposed in the literature. Finally, Section 4 presents the concluding remarks.

The RBF-FD Collocation Method
A mesh-free (or meshless) method adopts an algebraic system of equations for the complete domain without requiring a pre-defined mesh discretization of the domain and its boundary [43,44]. Mesh-free techniques are used to approximate scattered data, since generating meshes is one of the most laborious tasks of mesh-based numerical processes. Indeed, a mesh-free technique provides a low-cost alternative to schemes involving finite volume, finite difference, finite element, multivariate splines, and wavelets, all requiring node connectivity. Meshless techniques eliminate the mesh generation step and a collection of scattered data can be used. The RBF is one of the most widely used meshless techniques and reveals good performance in case of multidimensional scattered data interpolation [43,44].
Given a set of scattered node data X C = {x 1 , . . . , x N } ⊆ R n and the corresponding function values u i = u(x i ), i = 1, 2, . . . , N, the RBF interpolant is represented in the form . . , N, are RBF with shape parameter c, and the operation · 2 represents the Euclidean norm [44,45]. Some popular choices of RBF include the linear, Cubic, Multiquadric (MQ), Gaussian (GA), and thin-plate spline (TPS) versions with dependence r, r 3 , √ c 2 + r 2 , exp (−cr 2 ), and r 4 ln(r), respectively, where r = x − x j 2 . The coefficients {α j } N j=1 of Equation (6) are computed by imposing interpolation conditions S(x i ) = u i , i = 1, . . . , N. The relation (6) can be written in the following matrix form where The non-singularity of the associated linear system was proven in [46]. The main pros of the RBF collocation method when solving PDEs are its simplicity, easy application to different PDEs, and efficiency for solving problems involving complex domains. On the other hand, the major con of this method is related to the problem of full matrices. These matrices are strongly sensitive to the shape parameter c selected in the RBF and, therefore, they become difficult to solve in problems where we have too many unknowns. This problem arises from the fact that using the RBF interpolation increases the condition numbers of the related matrices for a large number of nodes. This occurs particularly when one selects inadequate data centers and uses basic functions that are infinitely smooth, such as the MQ, with extreme values of the shape parameter c [45].
The notation of local differentiation is popular in the RBF literature, particularly for time-dependent PDEs. The local radial basis function (RBF) generated by finite differences (RBF-FD), raised considerable interest owing to the structure of their differentiation and interpolation matrices [47,48]. It is possible to control the degree of sparsity of the differentiation and interpolation matrices produced by the local RBF. This sparsity can take advantage of parallelism and solve large problems [49,50]. The local RBFs have also been employed to reduce the model order. In some situations, researchers have found that the local RBF technique can produce the same degree of accuracy as the global RBF technique with a smaller mesh size [49][50][51][52][53]. Although small mesh sizes result in smaller ODE systems, the overall accuracy is maintained. Interested readers can find examples of the application of local RBFs to problems in the geosciences in [54,55]. Garshasbi et al. used the RBF collocation method for approximating the shallow water model named the Camassa-Holm equation [56]. Uddin connected the RBF to the pseudo-spectral method, known as RBF-PS method to approximate the equal width equation [57]. Nikan  implemented the local RBF-FD meshless method for generalized Korteweg-de Vries-Burgers [61] and Kawahara [62] equations.
Let us consider that In the local RBF-FD collocation method, the linear differential operator L at every point can be approximated only the stencil instead of applying the complete number of point, i.e., where x i = x i 1 is the center point of stencil I i . Figure 1 gives an example of a domain with 9 grid points and a stencil size of n i = 4. At the point x 3 , the n i − 1 = 3 nearest neighbors are used in the computation. Figure 2 shows the sparsity patterns for N = 50 for two stencil sizes n i = 11 and n i = 15. By deriving the RBF expansion of u(x) in Equation (8), the weighted differences of stencil node can be obtained from the system as: where Indeed, it is necessary to solve a small-sized linear system with a conditionally positive definite coefficient matrix in each stencil. The weighted differences of the stencil nodes w 1 , w 2 , . . . , w n i can be determined from the above system.
The first, second and third order derivatives can be approximated with the help of the function values at a set of n i nodes (including x i ) in the stencil of x i . That is, we can write where w x,l i,j represents the weighted differences of stencil node for the order derivatives l = {1, 2, 3}. We can obtain the following semi-discrete system by considering the notations above as The above equation can be represented as We must note that the matrices where (16) is of the form Equation (17) is an ODE with respect to u and it can be solved by means of an ODE solver in MATLAB such as ode113 or ode45. Let τ = T/M and t n = nτ, for n = 0, 1, . . . , M, so that the mesh {t n : n = 0, 1, . . . , M} is uniform. The initial solution u 0 is the starting vector. The package ode45 is an explicit Runge-Kutta of order 4(5) formula of the Dormand-Prince pairs [63]. The ode45 is a one-step solver that computes u t n given the solution at the preceding time point u t n−1 . On the other hand, the ode113 is a variable-order Predict-Evaluate-Correct-Evaluate solver of the Adams-Bashforth-Moulton type [64]. This solver might be more efficient than the ode45 for close tolerances and, in particular, when the ODE file function is particularly expensive to evaluate. A multi-step solver, such as the ode113, needs the solutions corresponding to more than one preceding time point for calculating the current solution. Hereafter, we calculate the differentiation matrices, expressed by W x , W xx and W xx , only one time outside the time-stepping operation. Additionally, merely matrix-vector multiplications are required within the time-stepping operation.

Stability Analysis
The method of lines represents the idea of using the FD method in the time direction t to solve a coupled system of ODEs. The numerical stability of the method of lines is investigated by a rule of tumb. The method of lines is stable if the eigenvalues of the (linearized) spatial discretization operator, scaled by τ, lie in the stability region of the time-discretization operator [57,65]. One defines the stability region as the portion of a multifaceted plane consisting of eigenvalues which result in the generation of bounded solutions. The coefficient matrix eigenvalues determine the stability of Equation (16) [66]. Hence, we need only to demonstrate that every eigenvalue Re(λ i ) belonging to the coefficient matrix has a non-positive real term Re(λ i ), where λ i , i = 1, 2, . . . , N, represents of the matrix eigenvalues. In other words, for all i = 1, 2, . . . , N, we must have Re(λ i ) ≤ 0 for obtaining stable solutions. The reader is referred to [66] for further details. In order to investigate the stable and unstable eigenvalue ranges of the Rosenau-KdV-RLW model, one must compute the eigenvalues belonging to the matrix W, which are scaled by τ.

Computational Results and Comparisons
This section considers five test problems assessing the effectiveness and accuracy of the proposed method for various values of h, τ and c. To measure the accuracy of method in comparison with the exact solution, we compute the following error norms: where u and u exact denote the numerical solution and exact solution, respectively. In addition, the invariants of motion are evaluated by It should be noted that the Gaussian function is used as a basis and the computations were performed in MATLAB R2016a with a computer system having a configuration including Intel(R) Core(TM) i5-2330 CPU 3.60 GHz and 8.00G RAM.  Table 1 lists the approximation errors in terms of L ∞ , L 2 and RMS with τ = 0.01 and n i = 489. Table 2 compares the obtained results with those provided by the techniques described in [38,41]. It is seen that the errors obtained by the proposed technique are inferior to the others. Figure 3 depicts the motion of the single solitary wave with h = τ = 0.125 over the spatial intervals x ∈ [−40, 60] (left) and x ∈ [−70, 100] (right) at final times T ∈ {0, 30, 40}. We verify that the single solitons move to the right at a constant speed and preserve their amplitude and shape with increasing time as anticipated. Figure 4 represents the absolute errors L ∞ at final times T ∈ {0, 30, 40}. Figure 5 portraits the eigenvalues of the linearized differentiation operator A and B (left and right panels, respectively) with N = 100. We observe that the eigenvalues calculated for A and B are zero or have negative values. The eigenvalues belonging to the linearized differentiation operators are real and negative or are complex with a negative real term. Hence, the stability of the proposed system for this case is proven. , k 23 = 5 + √ 34 10 · Table 3 reports the L ∞ , L 2 and RMS errors with τ = 0.01 and n i = 489. Table 4 compares the results with those obtained by the techniques described in [34,41]. It is clear that the results of the new method are considerably more accurate. Table 5 illustrates the conservative law of the discrete energy E. Figure 6 depicts the motion of single solitary wave with h = τ = 0.125 (left) and h = τ = 0.0625 (right) over the spatial interval x ∈ [−60, 90] at final times T ∈ {0, 10, 40}. The single solitons move to the right at a constant speed preserving their amplitude and shape. Figure 7 represents the absolute error L ∞ at final times T = {30, 40}. Figure 8 plots the eigenvalues of the matrices A and B (left and right panels, respectively) with N = 100. The eigenvalues calculated for A are negative values. For what concerns B, they are zero or have negative values. Therefore, the stability of the proposed system is confirmed.             Table 6 compares the results of the proposed method with those resulting from the schemes in [19,41]. The computational efficiency is clearly superior to the performance exhibited by the other schemes. Figure 9 plots the motion of single solitary wave with h = τ = 0.5, (left) h = τ = 0.25 (right) over the spatial interval x ∈ [−40, 100] at final times T ∈ {0, 30, 40}. The peak of the solitary waves remains the same during the simulation. Figure 10 shows the eigenvalues of the matrices A and B (left and right panel, respectively) with N = 100. The eigenvalues calculated for A and B have zero or negative values. Hence, the stability of the proposed system for this case is confirmed.     Table 7 compares the results of proposed method with those obtained with other schemes [19,25,33,42]. In this case, the accuracy of the method is slightly better than those achieved with the rest. Figure 11 depicts the motion of the single solitary wave with h = τ = 0.4 (left) and h = τ = 0.2 (right) over the spatial interval x ∈ [−50, 150] at final times T ∈ {8, 16, 24, 32}. The crest of the soliton clearly remains the same during the simulation. Figure 12 plots the eigenvalues of the matrices A and B (left and right panels, respectively) with N = 100. The eigenvalues calculated for A are negative values, while for B they have zero or negative values. Hence, the stability of the proposed system for this case is verified.
The initial boundary value problem (1)-(3) includes the following conservative quantities: related to mass and energy. The quantities I 1 and I 2 are applied to measure the conservation properties of the present method, calculated by means of the trapezoidal rule for the Rosenau-RLW equation. Tables 8 and 9 compare the results of the proposed method with those obtained from the schemes presented in [26,36,37,39]. It can be observed that the computational results are clearly better than the others and that the invariants I 1 and I 2 remain constant during the simulation. Figure 13 plots the motion of the single solitary wave for various p at T = {0, 30, 60} in the spatial interval x ∈ [−60, 120]. The single solitons move to the right at a constant speed and conserve their amplitudes and shapes. Figure 14 shows the eigenvalues of the linearized differentiation operator A and B (left and right panels, respectively) with N = 100. The eigenvalues calculated for A and B are zero, or have negative values. Therefore, the stability of the proposed system for this case is confirmed.   1.1875 × 10 −3 3.5725 × 10 −3 [36] 2.2547 × 10 −2 6.1304 × 10 −2

Conclusions
We adopted the local meshless RBF-FD to calculate the approximate numerical solutions of the general nonlinear Rosenau-RLW equation without performing any linearization or transformation of the equation. In order to demonstrate the accuracy of the proposed numerical technique, the error invariants and error norms were computed, and the results were compared with others available in the literature. The local RBF-FD technique was verified to be remarkably accurate. In conclusion, the method is sufficiently accurate and fast due to its limited computational load.
Author Contributions: All authors contributed equally to this paper. All authors read and approved the final paper.
Funding: This research received no external funding.