{'fraction_non_alphanumeric': 0.05486855759252854, 'fraction_numerical': 0.02602905569007264, 'mean_word_length': 4.32926267281106, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A linearly implicit conservative difference scheme for the generalized Rosenau-Kawahara-RLW equation

This paper concerns the numerical study for the generalized Rosenau-Kawahara-RLW equation obtained by coupling the generalized Rosenau-RLW equation and the generalized Rosenau-Kawahara equation. We first derive the energy conservation law of the equation, and then develop a three-level linearly implicit difference scheme for solving the equation. We prove that the proposed scheme is energy-conserved, unconditionally stable and second-order accurate both in time and space variables. Finally, numerical experiments are carried out to confirm the energy conservation, the convergence rates of the scheme and effectiveness for long-time simulation.


Introduction
The nonlinear wave is one of the most important scientific research areas. During the past several decades, many scientists developed different mathematical models to explain the wave behavior, such as the KdV equation, the RLW equation, the Rosenau equation, and many others. In the following, we give a short review of these important wave models.
The well known KdV equation was first introduced by Boussinesq [1] in 1877 and rediscovered by Diederik Korteweg and Gustav de Vries [2] in 1895. Since then, there are a lot of studies on this equation and its variational form. Here we just mention some of the recent work. Kudryashov [3] reviewed the travelling wave solutions for the KdV and the KdV-Burgers equations proposed by [4], Biswas [5] studied the solitary wave solution for KdV equation with power-law nonlinearity and timedependent coefficients, while Wang et al. [6] investigated the solitons, shock waves for the potential KdV equation. In addition to the theoretical studies, readers can refer to [7,8] for the numerical simulations of the KdV equation and the generalized KdV equation. The regularized long-wave (RLW) equation (also known as Benjamin-Bona-Mahony equation) was first proposed as a model for small-amplitude long wave of water in a channel by Peregrine [9,10]. The regularized long-wave (RLW) equation and its different variational forms were well studied both theoretically and numerically in the literature. Readers can refer to [11][12][13][14][15][16][17] for numerical studies and [18][19][20] for theoretical studies. When study the compact discrete systems, the well-known KdV equation can not describe the wave-wave and wave-wall interactions. To overcome the shortcoming of KdV equation, Rosenau proposed the following so-called Rosenau equation [21,22]: The existence and uniqueness of the solution for the Rosenau equation were theoretically proved by [23]. Besides the theoretical analysis, numerical studies of equation (3) also exist in the literature, see [24][25][26][27][28] and references therein. For further consideration of nonlinear waves, the term −u xxt is included in the Rosenau equation. The resulting equation is usually called the Rosenau-RLW equation [29][30][31]: The above equation was further extended into the generalized Rosenau-RLW equation [32][33][34]: where m ≥ 1 is a positive integer. The initial boundary value problems for the Rosenau-RLW equation and the generalized Rosenau-RLW equation have been well studied numerically in the past years [29][30][31][32][33][34]. For example, Pan et al. [29] developed a three-level linearly implicit conservative scheme for the usual Rosenau-RLW equation (4), where the method achieves second-order accuracy both in time and space variables. Zuo et al. [32] proposed a secondorder Crank-Nicolson finite difference scheme for the generalized Rosenau-RLW equation (5), however, the method is a nonlinear scheme, iterations are needed at each time step. Later on, Pan et al. [33] used a three-level linearly implicit conservative scheme for the same generalized Rosenau-RLW equation (5), where the method is a linear scheme and is also second-order convergent both in time and space variables.
On the other hand, to consider another behavior of nonlinear waves, the viscous term u xxx needs to be included in the Rosenau equation (3). The resulting equation is usually called the Rosenau-KdV equation [35][36][37][38][39]: For numerical investigations, Hu et al. [39] proposed a second-order linear conservative finite difference method for the Rosenau-KdV equation. However, numerical methods for the initial-boundary value problem of the Rosenau-KdV equation have not been studied widely. By coupling the above Rosenau-RLW equation and Rosenau-KdV equation, one can obtain the following Rosenau-KdV-RLW equation [40][41][42][43][44], For numerical investigation, Wongsaijai et al. [40] proposed a three level implicit conservative finite difference method for the above Rosenau-KdV-RLW equation. For theoretical studies, the solitary wave, shock waves, conservation laws as well as the asymptotic behavior for the Rosenau-KdV-RLW equation with power law nonlinearity were studied by [41][42][43][44], where the power law nonlinearity means the the last term in the left-hand side of equation (7) is replaced by a general nonlinear term (u p ) x and p is any positive integer. In addition, the following Kawahara equation arose in the theory of shallow water waves with surface tension [45]. Equation (8) [46][47][48]. In addition to the theoretical studies, readers can refer to [49][50][51] for the numerical studies of the Kawahara equation and the modified Kawahara equation.
As one more step consideration of the nonlinear wave, Zuo [38] obtained the Rosenau-Kawahara equation by adding another viscous term −u xxxxx to the Rosenau-KdV equation (6), and studied the solitary solution and periodic solution of the Rosenau-Kawahara equation. The Rosenau-Kawahara equation is given as follows [38]: For numerical study, Hu et al. [52] proposed a two level nonlinear Crank-Nicolson scheme and another three-level implicit linear conservative finite difference scheme for the Rosenau-Kawahara equation, both methods are proved to be second-order convergent. And Biswas [53] investigated the solitary solution and the two invariance of following equation, where a, b, c, ν, α, λ are real constants, m is a positive integer, which indicates the power law nonlinearity. Equation (10) is referred as the generalized Rosenau-Kawahara equation [52]. As far as we know, besides the analytical solitary solution obtained by [53], there is no numerical study for this generalized Rosenau-Kawahara equation.
Besides all the literature mentioned in the above, there are also a lot of related references which used both analytical and numerical methods to study different nonlinear models. For example, Triki et al. [54] discussed the solitary solution for the Gear-Grimshaw model, Bhrawy et al. [55] provided solitons, cnoidal waves and snoisal waves for the Whitham-Broer-Kaup system, while Shokri et al. [58] numerically solved the two-dimensional complex Ginzburg-Landau equation by using the meshless radial basis functions method, and Dehghan et al. [59] applied a semi-analytical method for solving the Rosenau-Hyman equation. Readers can refer to [54][55][56][57] for more discussions on analytically finding solitons and other wave solutions for different nonlinear models, and refer to [58][59][60] for more discussions on numerical studies.
In this paper, we will study the following generalized Rosenau-Kawahara-RLW equation: which is obtained by coupling the generalized Rosenau-RLW equation (5) and the generalized Rosenau-Kawahara equation (10). Analytical solitary solutions of the above equation (11) were recently obtained by [61] through the sech ansatz method.
In this work, we shall consider the numerical simulation of above equation (11) with initial condition and boundary conditions where m is a positive integer, a, b, c, ν are all real constants, while α, λ are set to be positive constants, x l is a large negative number and x r is a large positive number. In this paper, we only discuss the solitary solution of equation (11). By solitary wave assumptions, the solitary solution and its derivatives have the following asymptotic values: u → 0 as x → ±∞, and ∂ n u ∂x n → 0 as x → ±∞, for n ≥ 1. Thus, the boundary conditions (13) are meaningful for the solitary solution of equation (11). And we assume that the wave peak is initially located at x = 0, and x l , x r , which are large numbers, are used to assure that the solitary wave peak always locates inside the domain [x l , x r ] during the time interval [0, T ]. Similar set-ups are used in [39,40,52].
Note that equation (11) reduces to the generalized Rosenau-RLW equation when c = 0, ν = 0. Specially, it reduces to the usual Rosenau-RLW equation (4) when a = m = α = λ = 1, b = 2, c = 0, ν = 0. And equation (11) reduces to the generalized Rosenau-Kawahara equation when α = 0. Specially, it reduces to the usual Rosenau-Kawahara equation For numerically solving differential equations, the total accuracy of a particular method is affected not only the order of accuracy of the method, but also other factors. The conservative property of the method is another factor that has the same or possibly even more impact on results. For example, one successful and active research is to construct structure-preserving schemes (or called symplectic schemes) for the ODE systems (see [62] and the references therein). Better solutions can be expected from numerical schemes which have effective conservative properties rather than the ones which have nonconservative properties [63]. And Li et al. [64] even pointed out that in some areas, the ability to preserve some invariant properties of the original differential equation is a criterion to judge the success of a numerical simulation.
For the initial boundary value problem (11)-(13), we will show that it satisfies a fundamental energy conservation law in the next section. In addition, the wave equation is highly nonlinear due to the third nonlinear term in the left-hand side of the equation (11). When considering the finite difference scheme for the equation (11), the usual Crank-Nicolson scheme will lead to a nonlinear scheme with heavy computation while other standard linearized discretizations for the nonlinear term, e.g., one step Newton's method or a second-order extrapolation method, will loss the the energy conservative property. An ideal scheme should have relative less computational cost, can preserve energy, be unconditionally stable and maintain second-order accuracy.
In this paper, a three-level linearly implicit finite difference method for the initial value problem (11)-(13) will be presented. The fundamental energy conservation is preserved by the presented numerical scheme. The existence and uniqueness of the numerical solution are also proved. Moreover, numerical analysis shows that the method is secondorder convergent both in time and space variables, and the method is unconditionally stable. Numerical results confirm well with the theoretical results.
The rest of the paper is organized as follows: Section 2 shows the energy conservation of the initial boundary value problem (11)- (13). Section 3 gives the detail description of the three-level linearly implicit finite difference method, the proof for the discrete conservative property, the existence and uniqueness as well as the convergence and stability of the numerical solution. Numerical results are shown in section 4. Conclusions are provided in the final section.
, then the solution of (11)-(13) satisfies the energy conservation law: for

Proof. Multiplying (11) by 2u and integrating over the interval
Using the integration by parts, one can easily obtain where the boundary conditions (13) are used. Thus, only the first, fifth and sixth term in the left hand side of (15) are nonzero, all other terms are zero. This yields, Therefore, for any t ∈ [0, T ]. This completes the proof.

Numerical scheme
In this section, we give a complete description of our numerical method for the initial value problem (11)- (13). We first describe our solution domain and its grid. The computational domain is defined as For convenience, the difference operators, inner product and norms are defined as follows: where U n = (U n −1 , · · · , U n M+1 ) is the numerical solution at time t n = nτ. The essential of our scheme is that the third term in the left-hand side of (11) is rewritten and discretized as and which is a second-order approximation around (x i = x l + ih, t n = nτ). And all other terms of left-hand side of (11) are discretizated by using the standard second-order central difference method around (x i = x l + ih, t n = nτ) .
The detailed numerical scheme is as follows: and Obviously, the above conditions (21) will give U j 1 = U j M−1 = 0 on two internal points and U j −1 = U j M+1 = 0 on two fictitious points, for any 0 ≤ j ≤ N. Thus, U j ∈ Z 0 h , for any 0 ≤ j ≤ N. In addition, we can see that (19) is a three-level linearly implicit scheme and the coefficient matrix of the linear system (19) is banded. Therefore, the resulting linear system of equations can be solved efficiently using a direct linear solver, such as the LU decomposition method.
Since the scheme is a three-level method, to start the computation, we need to give the method for computation of U 1 . The U 1 is computed through the following Crank-Nicolson scheme: where it is a nonlinear scheme and is second-order accurate both in time and space variables. We point out that this discretization (18) for the third term in the left-hand side of (11) is specially designed so that the numerical scheme (19)- (21) can guarantee the energy conservation (see the proof of theorem 2 in the next subsection). Readers can refer to references [29,40,52] for the similar kind of treatments for the nonlinear term uu x in their equations to achieve conservative schemes. In addition, our scheme (19)- (22) is only nonlinear at first time step when computing U 1 , and all other steps are linear. Thus, the computational cost is relatively cheap. In the following several sections of the paper, we will show that our scheme is second-order accurate both in time and space variables, can preserve the energy identity (14), and is unconditionally stable.
The following lemmas are well known results, which are essential for existence, uniqueness, convergence, and stability of our numerical solution. Lemma 1 and 2 can be verified through direct computation, lemma 3 is the discrete form of the Sobolev's embedding theorem, which can be found in lemma 1, page 110 of [65]. In the rest part of the paper, unless otherwise indicated, C is the notation referring to a general positive constant, which may have different values in different contexts.

Discrete conservation
Theorem 2. Suppose u 0 ∈ C 7 0 [x l , x r ], then the solution of finite difference scheme (19)-(21) satisfies U n ∞ ≤ C, U n x ∞ ≤ C, for any 0 ≤ n ≤ N, and also satisfies the following discrete energy conservation: for any 0 ≤ n ≤ N − 1, where E n is the discrete energy at time t = (n + 1 2 )τ. ≤ N) into account, and after computing the inner product of equation (19) withŪ n , i. e., U n+1 +U n−1 2 , we have

Multiplying (22) both sides by
and using the similar techniques as above, one can obtain Thus, (31) can be rewritten as Since u 0 ∈ C 7 0 [x l , x r ] and the initial condition (20) are used in the numerical method, the right-side of (33) is bounded. By assumptions, α, λ are positive constants, therefore, By using lemma 3, we have U n ∞ ≤ C. In addition, through direct computation, one can verify that U n xx = U n xx , for any 0 ≤ n ≤ N. Thus, Again by using lemma 3, we have U n x ∞ ≤ C. This completes the proof.

Existence and uniqueness Theorem 3. The finite difference scheme (19)-(21) has a unique solution.
Proof. To prove the theorem, we proceed by the mathematical induction. Suppose U 1 , · · · , U n (1 ≤ n ≤ N − 1) are solved uniquely, we now consider the equation (19) for U n+1 . Assume that U n+1,1 , U n+1,2 are two solutions of (19) and let W n+1 = U n+1,1 − U n+1,2 , then it is easy to verify that W n+1 satisfies the following equation: Taking the inner product of (37) with W n+1 , we have where are used. The first five identities of (39) are directly from lemma 1 and 2, and the last one can be obtained as follows: From (38) and the definition of the · -norm, one can see that (38) has only a trivial solution. Thus, (19) determines U n+1 uniquely. This completes the proof.

Convergence and stability
Let u(x, t) be the solution of problem (11)-(13), U n i be the solution of the numerical schemes (19)- (21), and u n i = u(x i , t n ), e n i = u n i − U n i , then the truncation error of the scheme (19)- (21) can be obtained as follows: whereē n = e n+1 +e n−1 2 , 2 ≤ i ≤ M − 2 and 1 ≤ n ≤ N − 1. Since all terms in (19) are the second-order approximations of the corresponding terms in left-hand side of (11) around (x i = x l + ih, t n = nτ), by Taylor expansion, it can be easily obtained that r n i = O(τ 2 + h 2 ) if h, τ → 0 and u(x, t) ∈ C 7,3 , where C 7,3 is the set of functions which are seventh order continuous differentiable in space and third order continuous differentiable in time. This following lemma is a well known result.

Lemma 4. (Discrete Gronwall's inequality). Suppose that w(k) and ρ(k) are nonnegative functions while ρ(k) is a non-decreasing function. If
, and u(x, t) ∈ C 7,3 , then the numerical solution U n of the finite difference scheme (19)- (21) converges to the solution of the problem (11)- (13) in the sense of · ∞ , and the convergence rate is O(τ 2 + h 2 ), i.e., Proof. Taking the inner product of (41) with 2ē n , we have e n+1 2 − e n−1 2 +α( e n+1 x 2 − e n−1 where By using lemma 2, we obtain From the notations introduced at the beginning of section 3, for any 0 ≤ n ≤ N, we have the following inequality Thus, for any 0 ≤ n ≤ N, we have In addition, we have where theorem 2, Cauchy-Schwarz inequality and inequality (46) are used. Similarly, we have where theorem 2, Cauchy-Schwarz inequality and inequality (46) are used again. Furthermore, we have (r n , 2ē n ) ≤ 2 r n ē n ≤ r n 2 + ē n 2 ≤ r n 2 + 1 2 e n+1 2 + e n−1 2 , where Cauchy-Schwarz inequality are used. Substituting where C ′ = max( C α , C) and C is the positive constant in the above inequality. Replacing C ′ in the above inequality by the general positive constant notation C, we have then (51) can be rewritten as follows: which is equivalent to If τ, which is sufficiently small, satisfies τ < 1 3C (C is the positive constant in the inequality (52)), then 1 − Cτ > 0 and (52) gives where C ′′ = max(3C, 3) and we have used 2 1−Cτ < 3 since τ < 1 3C . Replacing C ′′ in the above inequality by the general positive constant notation C, we have Summing (53) from 1 to n − 1, we get where Since e 0 i = 0 and the Crank-Nicolson scheme (22) is used to compute U 1 , we have D 1 = O(τ 2 + h 2 ) followed by a simple analysis for the scheme (22). Therefore Using lemma 4, we obtain Thus, By using lemma 3, we have i.e., This completes the proof.
, then the solution U n of the finite difference scheme (19)- (21) is unconditionally stable with the · ∞ norm.
The proof of this theorem is similar as the above theorem.

Numerical results
In order to preform numerical accuracy test of the method proposed in this paper, one needs to know the exact solitary solutions beforehand. Exact solitary solutions of the equation (11) were previously obtained by [61] using the sech ansatz method. In this paper, we provide a different derivation through the sine-cosine method, which is listed in the appendix A. Example 1. We present the numerical results for the case m = 2, a = 1, b = 1, c = 2, α = 1, λ = 1, ν = 1. From [61] and appendix A, one can obtain the exact solitary solution as follows: and the initial condition is set as We first carry out the numerical convergence studies. For the spatial convergence, we set τ = 0.005 as the fixed time step and use 4 different spatial meshes: h = 0.8, 0.4, 0.2, 0.1, where τ is sufficient small such that the temporal error is negligible comparing to the spatial error (Here the time step is τ = 0.005 while the smallest spatial size is h = 0.1, thus h 2 >> τ 2 , therefore, the dominant errors are the spatial errors). The final time T is set to be 10, and x l = −40, x r = 200. Table 1 gives the errors between numerical solutions and exact solutions. We can see that the error decreases when the spatial mesh is refined and the convergence rate is two. Thus, the method is second-order convergent in space, which is consistent with theoretical results in the above section. For the temporal convergence, we set h = 0.005 as the fixed spatial mesh and use 4 different temporal meshes: τ = 0.8, 0.4, 0.2, 0.1, where h is sufficient small such that the spatial error is negligible comparing to the temporal error (Here the smallest time step is τ = 0.1 while the spatial size is h = 0.005, thus h 2 << τ 2 , therefore, the dominant errors are the temporal errors). The final time T is set to be 10, and x l = −40, x r = 200. Table 2 gives the errors between numerical solutions and exact solutions. Again, we can see that the error decreases when the temporal mesh is refined, and the convergence rate is also two. Thus, the method is second-order convergent in time, which is again consistent with theoretical results in the above section.
In order to show that the numerical scheme has the energy conservative property (23), we carry out another computation, where T = 100, x l = −40, x r = 240, h = 0.1, τ = 0.1 are used. Table 3 gives the quantities of E n at several time stages. We can see that E n is conserved exactly (up to 8 decimals) during the time evolution of the solitary wave.
Example 2. We present the numerical results for the case m = 4, a = 1, b = 1, c = 2, α = 1, λ = 1, ν = 1. From [61] and appendix A, one can obtain the exact solitary solution as follows: and the initial condition is set as Again, we carry out the spatial and temporal convergence. Table 4 and 5 give the errors between numerical solutions and exact solutions for spatial and temporal convergence, respectively. Once again, we can see that the method is second-order convergent both in time and space variables.
Additionally, table 6 provides the quantities of E n with T = 100, x l = −40, x r = 240, h = 0.1, τ = 0.1. Once again, we can see that E n is conserved and the method can be well used to study the solitary wave at long time.

Conclusions
In this paper, a three-level linearly implicit finite difference method for the initial value problem of the generalized Rosenau-Kawahara-RLW equation is developed. The fundamental energy conservative property is preserved by the current numerical scheme. The existence and uniqueness of the numerical solution are also proved. The method is shown to be second-order convergent both in time and space variables, and it is unconditionally stable. Moreover, exact solitary solutions are derived through sine-cosine method which is used for the numerical tests. Numerical results confirm well with the theoretical results.