A numerical investigation of the GRLW equation using lumped Galerkin approach with cubic B-spline

In this work, we construct the lumped Galerkin approach based on cubic B-splines to obtain the numerical solution of the generalized regularized long wave equation. Applying the von Neumann approximation, it is shown that the linearized algorithm is unconditionally stable. The presented method is implemented to three test problems including single solitary wave, interaction of two solitary waves and development of an undular bore. To prove the performance of the numerical scheme, the error norms \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{2}$$\end{document}L2 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${L_{\infty}}$$\end{document}L∞ and the conservative quantities \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I_{1}}$$\end{document}I1, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I_{2}}$$\end{document}I2 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I_{3}}$$\end{document}I3 are computed and the computational data are compared with the earlier works. In addition, the motion of solitary waves is described at different time levels.

the GEW equation is written as and the GRLW equation has the following form: in which physical boundary conditions U → 0 as x → ±∞, the subscripts t and x represent time and spatial differentiation, ε and p is a positive integer, µ is positive constant. The boundary and initial conditions are taken where f(x) is a localized disturbance inside the interval [a, b] and it will be considered later. In the fluid problems, U implies the vertical displacement of the water surface or similar physical quantity. In the plasma applications, U is denoted as negative of the electrostatic potential. That's why, the solitary wave solution of Eqs. (1), (2) and (3) helps us to understand the many physical phenomena with weak nonlinearity and dispersion waves such as nonlinear transverse waves in shallow water, ion-acoustic and magnetohydrodynamic waves in plasma and phonon packets in nonlinear crystals.
The RLW equation is obtained by taking p = 1 in GRLW equation (3). Up to now, many numerical including finite elements and analytical solution techniques have been presented on the RLW equation. The RLW equation was investigated with the growth of an undular bore by Peregrine (1966). Morrison et al. (1984) proposed the approximate analytical technique for the scattering of solitary waves of the RLW equation. Galerkin approach with linear, quadratic and quintic B-spline was used by Doğan (2002), Gardner et al. (1995) and Dağ et al. (2006). Collocation method was set up by Raslan (2001) and Saka et al. (2011) with quadratic and both sextic and septic B-splines functions. Esen and Kutluay (2006) obtained the numerical solution of the RLW equation with lumped Galerkin method using quadratic B-spline. Galerkin method with extrapolation techniques has been implemented to the RLW equation by Mei and Chen (2012). Later on, the RLW equation has been solved numerically by using von Neumann technique based on parametric quintic splines (Lin 2014).
If p = 2 in Eq. (3), the obtained equation is called as the modified regularized long wave (MRLW) equation. Finite element methods based on quintic, cubic and septic collocation were used for obtaining the numerical solution of the MRLW equation by Gardner et al. (1997), Khalifa et al. (2008) and Karakoç et al. (2014). Collocation method based on quintic B-spline functions with Rubin and Graves linearization technique was investigated for solving the MRLW equation by Karakoç et al. (2013). The MRLW equation was solved numerically by Ali (2009) using mesh free collocation method. Galerkin approach with cubic B-spline has been applied to MRLW equation by Karakoç et al. (2015).
When we consider the GRLW equation discussed here, there are some exact and numerical solution techniques on its. Hamdi et al. (2004) presented the exact solution technique. Numerical methods based on decomposition scheme, finite difference scheme and element free kp-Ritz were introduced for GRLW equation by Kaya (2004), EL-Danaf et al. (2014) and Guo et al. (2014). An approximate quasilinearization scheme was used to solve the GRLW equation with initial condition on the formation of undular bore by Ramos (2007). Roshan (2012) and Mohammadi (2015) have got the numerical results of the GRLW equation using finite element method based on Petrov Galerkin and exponential B-spline collocation. Also, Galerkin and lumped Galerkin method used here have been implemented to the EW, KdVB, Coupled KdV and MEW equations by Doğan (2005), Saka and Dağ (2009), Kutluay and Uçar (2013) and Esen (2006).
Inspired by the results of the applied numerical methods to similar type equations, we can say that lumped Galerkin approach is an accurate and efficient numerical technique. So, in this work, we have constructed the lumped Galerkin approach with cubic B-splines to get the numerical results of the GRLW equation.

A lumped Galerkin method
Firstly, the solution domain limited to a finite interval [a, b] Prenter (1975) described the cubic B-spline functions φ m (x), ( m= −1(1) N + 1), at the nodes x m which form a basis over the interval [a, b] by Each cubic B-spline φ m covers four finite intervals, hence each finite interval [x m , x m+1 ] is covered by four splines. The approximate solution U N (x, t) is denoted in terms of the cubic B-splines by in which the unknown time-dependent quantities δ j (t) will be computed by using the boundary and weighted residual conditions. Using the equality hη = x − x m such that 0 ≤ η ≤ 1, the finite interval [x m , x m+1 ] is converted to more easily workable interval [0, 1]. So, the cubic B-splines (5) depending on variable η over the gap [0, 1] are reorganized with Here we should mention that except for φ m−1 (x), φ m (x), φ m+1 (x) and φ m+2 (x), all cubic B-spline functions are null over the finite element [0, 1]. Thus, approximation function (6) in terms of element parameters δ m−1 , δ m , δ m+1 , δ m+2 and B-spline element shape functions φ m−1 , φ m , φ m+1 , φ m+2 can be expressed over the interval [0, 1] by The nodal values of U , U ′ , U ′′ with respect to the time parameters δ m are derived from B-splines (7) and trial function (8) as follows: where the superscript ′ and ′′ symbolize first and second derivative to η, respectively. When applying the Galerkin's approach with weight function W(x) to Eq. (3), the weak form of Eq. (3) is obtained as Implementing the change of variable x → η to integral (10), which yields where Ů is considered to be a constant over an element to simplify the integral. Applying partial integration once to (11), this leads to the following equality: in which = p(p + 1)Ů p and β = µ h 2 . Substituting cubic B-splines (7) instead of the weight function W(x) and trial function (8) into integral equation (12) forms where δ e = (δ m−1 , δ m , δ m+1 , δ m+2 ) T and the dot states differentiation to t, which can be written in matrix form by The element matrices are By considering together contributions from all elements, the matrix equation (14) takes the form where δ = (δ −1 , δ 0 , ..., δ N , δ N +1 ) T is a nodal parameters. The A, B, C and D are septadiagonal matrices and their line of m is where Implementing the forward finite difference δ = δ n+1 −δ n �t and Crank-Nicolson approach δ = 1 2 (δ n + δ n+1 ) to Eq. (15), we obtain the matrix system Using the boundary conditions given by Eq. (4), the (N + 3) × (N + 3) system (16) is reduced to (N + 1) × (N + 1) matrix system. Since the row m of A, B, C and D has seven elements, the system (16) comprises of the diagonal matrix with seven columns element (known as septa-diagonal matrix). The septa-diagonal matrix system can be solved by using Thomas algorithm (see subsection ). In this solution procedure, we need to two or three inner iterations δ n * = δ n + 1 2 δ n − δ n−1 at each time step to minimize the nonlinearity. After all of these processes, we can easily achieve the recurrence relationship between two time steps n and n + 1 which is an ordinary member of the matrix system (16) where To initiate the iteration, the initial vector δ 0 must be calculated by using the initial and boundary conditions. Also, using the relations at the knots U N (x m , 0) = U (x m , 0), m = 0, 1, . . . , N and derivative condition U ′ N (x 0 , 0) = U ′ (x N , 0) = 0 together with a variant of the Thomas algorithm, the initial vector δ 0 can be easily computed from the following matrix form

The solution of septa-diagonal matrix system with Thomas algorithm
As used in Fortran program and given by Zaki (2000), the solution method of septadiagonal matrix system with Thomas algorithm is expressed as follows: The septa-diagonal system can be written by In the first step, the parameters are organized with and As a second step, we calculate the following parameters Now we obtain the solution

Stability analysis
In order to determine the linear stability analysis of the numerical algorithm, we use the Fourier method and assume that the quantity U p in the non-linear term U p U x of GRLW equation is locally constant. Substituting the Fourier mode δ n m = g n e imkh where k is mode number, h is the element size and i = √ −1, into the scheme (17) The modulus of |g| is 1, so the linearized scheme is unconditionally stable.

Numerical examples and results
In this section, we have applied the lumped Galerkin method to three test problems including single solitary wave, interaction of two solitary waves and development of an undular bore. These three examples are formed by using different values of initial condition. To demonstrate the efficiency and accuracy of the presented numerical scheme, the L 2 and L ∞ error norms are calculated by using the solitary wave solution in Eq. (22) and the following equalities: Furthermore, so as to indicate that the numerical approach keeps the properties related to mass, momentum and energy, we observe the changes of the invariants The exact solution of GRLW equation given in Gardner et al. (1997) and Roshan (2012) has the form  where p c(p+2) 2p is amplitude, c + 1 is the speed of the wave traveling in the positive direction of the x-axis, x 0 is arbitrary constant.

The motion of single solitary wave
For this problem, we use the initial condition obtained by taking t = 0 in Eq. (22). To coincide with papers Dağ et al. (2006), Gardner et al. (1997), Khalifa et al. (2008), Ali (2009), Karakoç et al. (2013), Roshan (2012) and Mohammadi (2015),  Table 1 that the changes of the invariants are less than 0.04, 0.05 and 0.05 %, respectively. In Table 2, three  In the second case, we take the parameters p = 3, c = 1.2, h = 0.1, t = 0.025 and p = 3, c = 0.3, h = 0.1, t = 0.01. These produce the amplitude = 1 and amplitude = 0.6. The calculated quantities are presented in Tables 3 and 4. As can be seen in Table 3, the changes of the invariants are less than 0.5, 0.7 and 0.7 % . Table 4 shows that three invariants are almost constant as the time increases. Also, we observe that the quantities of the error norms L 2 and L ∞ are reasonably small, as expected.
Thirdly, if p = 4, c = 4/3, h = 0.1, t = 0.01 and p = 4, c = 0.3, h = 0.1, t = 0.01, the solitary wave has amplitude = 1 and 0.6. The obtained results are reported in Tables 5 and 6. Table 5 denotes that the changes of the invariants are less than 0.2, 0.3 and 0.3 % . On the other hand, this change is too little in Table 6. As in the parameters of p = 2, 3, the quantities of the error norms L 2 and L ∞ are sensibly small.
Finally, we study the parameters p = 2, 3, 4, 6, 8, 10 with c = 0.03 and c = 0.1, h = 0.1, t = 0.01. The calculated values are listed in Table 7 which clearly shows that the error norms are sufficiently small and remain less than 5.2 × 10 −3 with increasing time, p and c. In addition, the motion of single solitary wave is displayed at different times and the values of p in Fig. 1. From this figure, we can see that the solitary wave moves to the right at constant velocity and remains its shape and amplitude. When the values of p are increased, the peak position of single solitary wave rises.    In Table 8, we compare the quantity of invariants and error norms obtained by presented scheme with the ones given by earlier methods. From the table, we can conclude that three invariants are to be close to each other. The magnitude of our error norms is smaller than the ones given by Gardner et al. (1997), Khalifa et al. (2008), Ali (2009) and Roshan (2012) for p = 2 and it is almost same with the paper (Roshan 2012) for p = 3, 4.

The interaction of two solitary waves
In the second test problem, we have worked on which provides two positive solitary waves having different amplitudes of magnitudes 2 and 1 at the same direction, where c i and x i , i = 1, 2 are arbitrary constants.
The parameters are chosen to be first values p = 2, c 1 = 4, c 2 = 1, x 1 = 25, x 2 = 55 , h = 0.2, t = 0.025, µ = 1, x ∈ [0, 250]; second values p = 3, c 1 = 48/5, c 2 = 6/5 , x 1 = 20, x 2 = 50, h = 0.1, t = 0.01, µ = 1, x ∈ [0, 120] and third values p = 4, c 1 = 64/3, c 2 = 4/3, x 1 = 20, x 2 = 80, h = 0.125, t = 0.01, µ = 1, x ∈ [0, 200]. The numerical computations are given in Tables 9 and 10. The results in Tables show that the changes of the invariants from their initial state are as small as required and good agreement with those of Roshan (2012). The motion of two solitary waves is simulated   at different time levels in Figs. 2 and 3. These figures show that the initial position of the wave with larger amplitude is on the left of the second wave with smaller amplitude. As the time processes, the large wave catches up with the smaller one and overlapping process occurs. After a while, waves start to resume their original forms.

The development of an undular bore
As a last test problem, we have focused on the development of an undular bore given by which indicates the elevation of the water above the equilibrium surface at time zero. The change in water level of magnitude Eq. (24) is centered on x = x c . We study with the parameters U 0 = 0.1, µ = 1/6, x c = 0, d = 5, h = 0.1, �t = 0.1, x ∈ [−36, 300] to be consistent with earlier works (Peregrine 1966;Esen and Kutluay 2006;Mei and Chen 2012;Doğan 2005). The conservative quantities are recorded in Table 11. In this table, the changes of the invariants remain less than 1.1 × 10 −2 , 1.0 × 10 −3 and 2.0 × 10 −3 , respectively. The undulation profiles are depicted at time t = 50 and t = 200 when p = 2, 3, 4 in Figs. 4, 5 and 6. It is understood that the magnitude of the waves increases with rising the value of x. Later, undulations take the peak position and disappear.

Conclusion
The solitary-wave solutions of the GRLW equation have been successfully obtained by using lumped Galerkin method based on cubic B-spline functions. Also, the linearized scheme has been found to be unconditonally stable. The error norms L 2 , L ∞ and three conservative quantities I 1 , I 2 and I 3 have been computed for single solitary wave, interaction of two solitary waves and development of an undular bore. These computations demonstrate that our error norms are as small as required and they are smaller than the most of existing numerical calculations or too close to the best result in literature. The numerical algorithm conserves the properties related to mass, momentum and energy and the numerical values of them have been found to be in good agreement with earlier studies. In addition, the profiles of the solitary wave are similar to those of references. As a result, we can say that lumped Galerkin method is more practical, accurate and productive numerical approximation technique for GRLW equation and it can be reliably used to solve the similar type non-linear problems.