An Efficient Approach to Numerical Study of the MRLW Equation with B-Spline Collocation Method

Abstract

A septic B-spline collocation method is implemented to find the numerical solution of the modified regularized long wave (MRLW) equation. Three test problems including the single soliton and interaction of two and three solitons are studied to validate the proposed method by calculating the error norms 𝐿₂ and 𝐿_∞ and the invariants 𝐼₁, 𝐼₂, and 𝐼₃. Also, we have studied the Maxwellian initial condition pulse. The numerical results obtained by the method show that the present method is accurate and efficient. Results are compared with some earlier results given in the literature. A linear stability analysis of the method is also investigated.

1. Introduction

The generalized regularized long wave equation is given by 𝑈_𝑡+𝑈_𝑥+𝛿𝑈^𝑝𝑈_𝑥−𝜇𝑈_𝑥𝑥𝑡=0, (1) where 𝑝 is a positive integer and 𝛿 and 𝜇 are positive constants. This equation is one of the most important nonlinear wave equation used a model for small amplitude long waves on the surface of water in a channel [1, 2]. A few authors solved the equation numerically: among others, Zhang [3] used a finite difference method for a Cauchy problem and Kaya [4] applied the Adomian decomposition method and a quasilinearization method based on finite differences was used by Ramos [5]. Roshan [6] implemented the Petrov-Galerkin method using a linear hat function as the trial function and a quintic B-spline function as the test function. A mesh-free technique for the numerical solution of the equation has been presented by Mokhtari and Mohammadi [7]. For   𝑝=1, 𝑈_𝑡+𝑈_𝑥+𝛿𝑈𝑈_𝑥−𝜇𝑈_𝑥𝑥𝑡=0. (2) Equation (2) is known as regularized long wave equation, originally introduced to describe the behavior of the undular bore by Peregrine [1] and later widely studied by Benjamin et al. [8]. The RLW equation has been solved numerically by finite element methods [9–22], finite difference methods [23–26], Fourier pseudospectral [27], and mesh-free method [28]. For  𝑝=2, 𝑈_𝑡+𝑈_𝑥+6𝑈²𝑈_𝑥−𝜇𝑈_𝑥𝑥𝑡=0. (3) Another particular case of (2) is called modified regularized long wave (MRLW) equation. Like RLW equation, the MRLW equation has been solved numerically by various methods. Among many others, a collocation solution to the equation using quintic B-spline finite element method is developed by Gardner et al. [29]. Khalifa et al. [30, 31] obtained the numerical solutions of the equation using finite difference method and cubic B-spline collocation finite element method. Solutions based on collocation method with quadratic B-spline finite elements and the central finite difference method for time are investigated by Raslan [32]. The equation was solved with a collocation finite element method using quadratic, cubic, quartic, and quintic B-splines to obtain the numerical solutions of the single solitary wave by Raslan and Hassan [33]. Haq et al. [34] have designed a numerical scheme based on quartic B-spline collocation method for the numerical solution of the equation. Ali [35] has formulated a classical radial basis functions (RBFs) collocation method for solving the equation. Karakoc et al. [36] have obtained a type of the quintic B-spline collocation procedure in which nonlinear term in the equation is linearized by using the form introduced by the Rubin and Graves [37] to solve the equation. A Petrov-Galerkin method using cubic B-spline function as trial function and a quadratic B-spline function as the test function is set up to solve the equation by Karakoc and Geyikli [38]. A homotopy analysis method has been employed to obtain approximate numerical solution of the modified regularized long wave (MRLW) equation with some specified initial conditions by Khan et al. [39].

In the present paper, a numerical scheme based on the septic B-spline collocation method has been set up for solving the MRLW equation with a variant of both initial and boundary conditions. This paper is set out as follows. In Section 2, septic B-spline collocation scheme is presented. Also stability analysis is considered. In Section 3, test problems including single, two, and three solitary waves and Maxwellian initial condition are discussed. Finally in Section 4, a summary is given at the end of the paper.

2. Septic B-Spline Finite Element Solution

Consider the MRLW Equation (3) given with the following boundary conditions, 𝑈()()𝑈𝑎,𝑡=0,𝑈𝑏,𝑡=0,_𝑥(𝑎,𝑡)=𝑈_𝑥()𝑈𝑏,𝑡=0,_𝑥𝑥(𝑎,𝑡)=𝑈_𝑥𝑥()𝑈𝑏,𝑡=0,_𝑥𝑥𝑥(𝑎,𝑡)=𝑈_𝑥𝑥𝑥(𝑏,𝑡)=0,𝑡>0, (4) and the initial condition 𝑈(𝑥,0)=𝑓(𝑥)𝑎≤𝑥≤𝑏. (5)

In order to be able to apply the numerical method, the solution region of the problem is restricted over an interval 𝑎≤𝑥≤𝑏. Space interval [𝑎,𝑏] is partitioned into uniformly sized finite elements of length ℎ by the nodes 𝑥_𝑚 such that 𝑎=𝑥₀<𝑥₁<⋯<𝑥_𝑁−1<𝑥_𝑁=𝑏 and ℎ=𝑥_𝑚+1−𝑥_𝑚,𝑚=−3,−2,…,𝑁+2,𝑁+3. The set of septic B-spline functions {𝜙₋₃(𝑥),𝜙₋₂(𝑥),𝜙₋₁(𝑥),…,𝜙_𝑁+1(𝑥),𝜙_𝑁+2(𝑥),𝜙_𝑁+3(𝑥)} forms a basis over the problem domain [𝑎,𝑏]. A global approximation 𝑈_𝑁(𝑥,𝑡) is expressed in terms of septic B-splines as 𝑈_𝑁(𝑥,𝑡)=^𝑁+3_𝑖=−3𝜙_𝑖(𝑥)𝛿_𝑖(𝑡), (6) where 𝛿_𝑖(𝑡)s are time dependent parameters to be determined from the initial, boundary, and collocation conditions.

Septic B-splines 𝜙_𝑚(𝑥), (𝑚=−3(1)𝑁+3), at the knots 𝑥_𝑚 are defined over the interval [𝑎,𝑏] by [40] 𝜙_𝑚(𝑥)=1ℎ⁷𝑥−𝑥_𝑚−4⁷𝑥_𝑚−4,𝑥_𝑚−3𝑥−𝑥_𝑚−4⁷−8𝑥−𝑥_𝑚−3⁷𝑥_𝑚−3,𝑥_𝑚−2𝑥−𝑥_𝑚−4⁷−8𝑥−𝑥_𝑚−3⁷+28𝑥−𝑥_𝑚−2⁷𝑥_𝑚−2,𝑥_𝑚−1𝑥−𝑥_𝑚−4⁷−8𝑥−𝑥_𝑚−3⁷+28𝑥−𝑥_𝑚−2⁷−56𝑥−𝑥_𝑚−1⁷𝑥_𝑚−1,𝑥_𝑚𝑥_𝑚+4−𝑥⁷−8𝑥_𝑚+3−𝑥⁷+28𝑥_𝑚+2−𝑥⁷−56𝑥_𝑚+1−𝑥⁷𝑥_𝑚,𝑥_𝑚+1𝑥_𝑚+4−𝑥⁷−8𝑥_𝑚+3−𝑥⁷+28𝑥_𝑚+2−𝑥⁷𝑥_𝑚+1,𝑥_𝑚+2𝑥_𝑚+4−𝑥⁷−8𝑥_𝑚+3−𝑥⁷𝑥_𝑚+2,𝑥_𝑚+3𝑥_𝑚+4−𝑥⁷𝑥_𝑚+3,𝑥_𝑚+40otherwise. (7)

Using expansion (6) and trial function (7), the nodal values (𝑈_𝑁)_𝑚 and their first, second, and third derivatives (𝑈_𝑁)_𝑥,(𝑈_𝑁)_𝑥𝑥,(𝑈_𝑁)_𝑥𝑥𝑥 can be calculated at the nodal points 𝑥_𝑚 in terms of nodal parameters 𝛿_𝑚 by the following set of equations: 𝑈_𝑁𝑚=𝛿_𝑚−3+120𝛿_𝑚−2+1191𝛿_𝑚−1+2416𝛿_𝑚+1191𝛿_𝑚+1+120𝛿_𝑚+2+𝛿_𝑚+3,𝑈_𝑁𝑥𝑚=7ℎ−𝛿_𝑚−3−56𝛿_𝑚−2−245𝛿_𝑚−1+245𝛿_𝑚+1+56𝛿_𝑚+2+𝛿_𝑚+3,𝑈_{𝑁𝑥𝑥𝑚}=42ℎ²𝛿_𝑚−3+24𝛿_𝑚−2+15𝛿_𝑚−1−80𝛿_𝑚+15𝛿_𝑚+1+24𝛿_𝑚+2+𝛿_𝑚+3,𝑈_{𝑁𝑥𝑥𝑥𝑚}=210ℎ³−𝛿_𝑚−3−8𝛿_𝑚−2+19𝛿_𝑚−1−19𝛿_𝑚+1+8𝛿_𝑚+2+𝛿_𝑚+3. (8) The splines 𝜙_𝑚(𝑥)  and their six principle derivatives vanish outside the interval [𝑥_𝑚−4,𝑥_𝑚+4].

To apply the proposed method, Crank-Nicolson approximation for the space derivatives 𝑈_𝑥 and 𝑈_𝑥𝑥 and usual first order forward difference formula for the time derivative of the 𝑈 in (3) have been used, which lead to 𝑈^𝑛+1−𝑈^𝑛+𝑈Δ𝑡^𝑥𝑛+1+𝑈^𝑛_𝑥2𝑈+6²𝑈_𝑥^𝑛+1+𝑈²𝑈_𝑥^𝑛2𝑈−𝜇^{𝑛+1𝑥𝑥}−𝑈^𝑛_𝑥𝑥Δ𝑡=0. (9)

In order to linearize the nonlinear term (𝑈²𝑈_𝑥)^𝑛+1, we can write the term as follows, 𝑈²𝑈_𝑥=𝑈𝑈𝑈_𝑥, (10) and apply the linearization form introduced by Rubin and Graves [37] 𝑈𝑈_𝑥^𝑛+1=𝑈^𝑛+1𝑈^𝑛_𝑥+𝑈^𝑛𝑈_𝑥^𝑛+1−𝑈^𝑛𝑈^𝑛_𝑥 (11) to (9),  and we get 𝑈^𝑛+1+Δ𝑡2𝑈_𝑥^𝑛+1+3Δ𝑡𝑈^𝑛+1×𝑈^𝑛+1𝑈^𝑛_𝑥+𝑈^𝑛𝑈_𝑥^𝑛+1−𝑈^𝑛𝑈^𝑛_𝑥−𝜇𝑈^𝑛+1_𝑥𝑥=𝑈^𝑛−Δ𝑡2𝑈^𝑛_𝑥−3Δ𝑡𝑈^𝑛𝑈𝑈_𝑥^𝑛−𝜇𝑈^𝑛_𝑥𝑥. (12)

Substituting the approximate solution 𝑈_𝑁 and putting the nodal values of 𝑈 and its derivatives given by (8) into (12) one obtains the following iterative system for 𝑚=0,1,…,𝑁∶𝛿^𝑛+1_𝑚−31+𝛼₁−𝛼₂−𝛼₃+𝛿^𝑛+1_𝑚−2120+120𝛼₁−56𝛼₂−24𝛼₃+𝛿^𝑛+1_𝑚−11191+1191𝛼₁−245𝛼₂−15𝛼₃+𝛿_𝑚^𝑛+12416+2416𝛼₁+80𝛼₃+𝛿^𝑛+1_𝑚+11191+1191𝛼₁+245𝛼₂−15𝛼₃+𝛿^𝑛+1_𝑚+2120+120𝛼₁+56𝛼₂−24𝛼₃+𝛿^𝑛+1_𝑚+31+𝛼₁+𝛼₂−𝛼₃=𝛿^𝑛_𝑚−31−𝛼₁+𝛼₂−𝛼₃+𝛿^𝑛_𝑚−2120−120𝛼₁+56𝛼₂−24𝛼₃+𝛿^𝑛_𝑚−11191−1191𝛼₁+245𝛼₂−15𝛼₃+𝛿^𝑛_𝑚2416−2416𝛼₁+80𝛼₃+𝛿^𝑛_𝑚+11191−1191𝛼₁−245𝛼₂−15𝛼₃+𝛿^𝑛_𝑚+2120−120𝛼₁−56𝛼₂−24𝛼₃+𝛿^𝑛_𝑚+31−𝛼₁−𝛼₂−𝛼₃, (13) where 𝛼₁=3Δ𝑡,𝛼₂=7Δ𝑡2ℎ,𝛼₃=42𝜇ℎ². (14) The newly obtained iterative system (13) consists of 𝑁+1 linear equation in 𝑁+7 unknowns (𝛿₋₃,𝛿₋₂,𝛿₋₁,𝛿₀,…,𝛿_𝑁,𝛿_𝑁+1,𝛿_𝑁+2,𝛿_𝑁+3)^𝑇. To obtain a unique solution of this system, six additional constraints are required. Applying the boundary conditions (4) and using the values of (8), these constraints are used and this enables us to eliminate the unknowns 𝛿₋₃,𝛿₋₂,𝛿₋₁,𝛿₀,…,𝛿_𝑁,𝛿_𝑁+1,𝛿_𝑁+2,𝛿_𝑁+3 from system (13). So system (13) is reduced to a septa-diagonal system of (𝑁+1) linear equations in (𝑁+1) unknowns given by 𝐴𝛿^𝑛+1=𝐵𝛿^𝑛, where 𝛿^𝑛=(𝛿₀,𝛿₁,𝛿₂,…,𝛿_𝑁−2,𝛿_𝑁−1,𝛿_𝑁)^𝑇. The coefficient matrixes are given by𝑎𝐴=¹¹𝑎₁₂𝑎₁₃𝑎₁₄𝑎₂₁𝑎₂₂𝑎₂₃𝑎₂₄𝑎₂₅𝑎₃₁𝑎₃₂𝑎₃₃𝑎₃₄𝑎₃₅𝑎₃₆⋱𝑎_{𝑚,𝑚−3}𝑎_{𝑚,𝑚−2}𝑎_{𝑚,𝑚−1}𝑎_𝑚,𝑚𝑎_𝑚,𝑚+1𝑎_𝑚,𝑚+2𝑎_𝑚,𝑚+3⋱𝑎_{𝑛−1,𝑛−4}𝑎_{𝑛−1,𝑛−3}𝑎_{𝑛−1,𝑛−2}𝑎_{𝑛−1,𝑛−1}𝑎_{𝑛−1,𝑛}𝑎_{𝑛−1,𝑛+1}𝑎_{𝑛,𝑛−3}𝑎_{𝑛,𝑛−2}𝑎_{𝑛,𝑛−1}𝑎_𝑛,𝑛𝑎_𝑛,𝑛+1𝑎_{𝑛+1,𝑛−2}𝑎_{𝑛+1,𝑛−1}𝑎_𝑛+1,𝑛𝑎_{𝑛+1,𝑛+1},𝑏𝐵=¹¹𝑏₁₂𝑏₁₃𝑏₁₄𝑏₂₁𝑏₂₂𝑏₂₃𝑏₂₄𝑏₂₅𝑏₃₁𝑏₃₂𝑏₃₃𝑏₃₄𝑏₃₅𝑏₃₆⋱𝑏_{𝑚,𝑚−3}𝑏_{𝑚,𝑚−2}𝑏_{𝑚,𝑚−1}𝑏_𝑚,𝑚𝑏_𝑚,𝑚+1𝑏_𝑚,𝑚+2𝑏_𝑚,𝑚+3⋱𝑏_{𝑛−1,𝑛−4}𝑏_{𝑛−1,𝑛−3}𝑏_{𝑛−1,𝑛−2}𝑏_{𝑛−1,𝑛−1}𝑏_{𝑛−1,𝑛}𝑏_{𝑛−1,𝑛+1}𝑏_{𝑛,𝑛−3}𝑏_{𝑛,𝑛−2}𝑏_{𝑛,𝑛−1}𝑏_𝑛,𝑛𝑏_𝑛,𝑛+1𝑏_{𝑛+1,𝑛−2}𝑏_{𝑛+1,𝑛−1}𝑏_𝑛+1,𝑛𝑏_{𝑛+1,𝑛+1}, (15)where 𝑎₁₁=−41472𝛼297³,𝑏₁₁=−41472𝛼297³,𝑎₁₂=−73224𝛼297³,𝑏₁₂=−73224𝛼297³,𝑎₁₃=−20736𝛼297³,𝑏₁₃=−20736𝛼297³,𝑎₁₄=−648𝛼297³,𝑏₁₄=−648𝛼297³,𝑎₂₁=215411+297215411𝛼297¹−11761𝛼297²+17893𝛼297³,𝑏₂₁=215411−297215411𝛼297¹+11761𝛼297²+17893𝛼297³,𝑎₂₂=616822.5+297616822.5𝛼297¹+43513.5𝛼297²+38665.5𝛼297³,𝑏₂₂=616822.5−297616822.5𝛼297¹−43513.5𝛼297²+38665.5𝛼297³,𝑎₂₃=340284+297340284𝛼297¹+78528𝛼297²−2532𝛼297³,𝑏₂₃=340284−297340284𝛼297¹−78528𝛼297²−2532𝛼297³,𝑎₂₄=35525.5+29735525.5𝛼297¹+16682.5𝛼297²−7109.5𝛼297³,𝑏₂₄=35525.5+29735525.5𝛼297¹+16682.5𝛼297²−7109.5𝛼297³,𝑎₂₅=1+𝛼₁+𝛼₂−𝛼₃,𝑏₂₅=1−𝛼₁−𝛼₂−𝛼₃,𝑎₃₁=34432+29734432𝛼297¹−15424𝛼297²−5920𝛼297³,𝑏₃₁=34432−29734432𝛼297¹+15424𝛼297²−5920𝛼297³,𝑎₃₂=352833+297352833𝛼297¹−71871𝛼297²−3561𝛼297³,𝑏₃₂=352833−297352833𝛼297¹+71871𝛼297²−3561𝛼297³,𝑎₃₃=717432+297717432𝛼297¹+120𝛼297²+23880𝛼297³,𝑏₃₃=717432−297717432𝛼297¹−120𝛼297²+23880𝛼297³,𝑎₃₄=353726+297353726𝛼297¹+72766𝛼297²−4454𝛼297³,𝑏₃₄=353726−297353726𝛼297¹−72766𝛼297²−4454𝛼297³,𝑎₃₅=120+120𝛼₁+56𝛼₂−24𝛼₃,𝑏₃₅=120−120𝛼₁−56𝛼₂−24𝛼₃,𝑎₃₆=1+𝛼₁+𝛼₂−𝛼₃,𝑏₃₆=1+𝛼₁+𝛼₂−𝛼₃,𝑎_{𝑚,𝑚−3}=1+𝛼₁−𝛼₂−𝛼₃,𝑏_{𝑚,𝑚−3}=1−𝛼₁+𝛼₂−𝛼₃,𝑎_{𝑚,𝑚−2}=120+120𝛼₁−56𝛼₂−24𝛼₃,𝑏_{𝑚,𝑚−2}=120−120𝛼₁+56𝛼₂−24𝛼₃,𝑎_{𝑚,𝑚−1}=1191+1191𝛼₁−245𝛼₂−15𝛼₃,𝑏_{𝑚,𝑚−1}=1191−1191𝛼₁+245𝛼₂−15𝛼₃,𝑎_𝑚,𝑚=2416+2416𝛼₁+80𝛼₃,𝑏_𝑚,𝑚=2416−2416𝛼₁+80𝛼₃,𝑎_𝑚,𝑚+1=1191+1191𝛼₁+245𝛼₂−15𝛼₃,𝑏_𝑚,𝑚+1=1191−1191𝛼₁−245𝛼₂−15𝛼₃,𝑎_𝑚,𝑚+2=120+120𝛼₁+56𝛼₂−24𝛼₃,𝑏_𝑚,𝑚+2=120−120𝛼₁−56𝛼₂−24𝛼₃,𝑎_𝑚,𝑚+3=1+𝛼₁+𝛼₂−𝛼₃,𝑏_𝑚,𝑚+3=1−𝛼₁−𝛼₂−𝛼₃,(1)𝑎𝑚=4𝑛−2,_{𝑛−1,𝑛−4}=1+𝛼₁−𝛼₂−𝛼₃,𝑏_{𝑛−1,𝑛−4}=1−𝛼₁+𝛼₂−𝛼₃,𝑎_{𝑛−1,𝑛−3}=120+120𝛼₁−56𝛼₂−24𝛼₃,𝑏_{𝑛−1,𝑛−3}=120−120𝛼₁+56𝛼₂−24𝛼₃,𝑎_{𝑛−1,𝑛−2}=353726+297353726𝛼297¹−72766𝛼297²−4454𝛼297³,𝑏_{𝑛−1,𝑛−2}=353726−297353726𝛼297¹+72766𝛼297²−4454𝛼297³,𝑎_{𝑛−1,𝑛−1}=717432+297717432𝛼297¹−120𝛼297²+23880𝛼297³,𝑏_{𝑛−1,𝑛−1}=717432−297717432𝛼297¹+120𝛼297²+23880𝛼297³,𝑎_{𝑛−1,𝑛}=352833+297352833𝛼297¹+71871𝛼297²−3561𝛼297³,𝑏_{𝑛−1,𝑛}=352833−297352833𝛼297¹−71871𝛼297²−3561𝛼297³,𝑎_{𝑛−1,𝑛+1}=34432+29734432𝛼297¹+15424𝛼297²−5920𝛼297³,𝑏_{𝑛−1,𝑛+1}=34432−29734432𝛼297¹−15424𝛼297²−5920𝛼297³,𝑎_{𝑛,𝑛−3}=1+𝛼₁−𝛼₂−𝛼₃,𝑏_{𝑛,𝑛−3}=1−𝛼₁+𝛼₂−𝛼₃,𝑎_{𝑛,𝑛−2}=35525.5+29735525.5𝛼297¹−16682.5𝛼297²−7109.5𝛼297³,𝑏_{𝑛,𝑛−2}=35525.5−29735525.5𝛼297¹+16682.5𝛼297²−7109.5𝛼297³,𝑎_{𝑛,𝑛−1}=340284+297340284𝛼297¹−78528𝛼297²−2532𝛼297³,𝑏_{𝑛,𝑛−1}=340284−297340284𝛼297¹+78528𝛼297²−2532𝛼297³,𝑎_𝑛,𝑛=616822.5+297616822.5𝛼297¹−43513.5𝛼297²+38665.5𝛼297³,𝑏_𝑛,𝑛=616822.5−297616822.5𝛼297¹+43513.5𝛼297²+38665.5𝛼297³,𝑎_𝑛,𝑛+1=215411+297215411𝛼297¹+11761𝛼297²+17893𝛼297³,𝑏_𝑛,𝑛+1=215411−297215411𝛼297¹−11761𝛼297²+17893𝛼297³,𝑎_{𝑛+1,𝑛−2}=−648𝛼297³,𝑏_{𝑛+1,𝑛−2}=−648𝛼297³,𝑎_{𝑛+1,𝑛−1}=−20736𝛼297³,𝑏_{𝑛+1,𝑛−1}=−20736𝛼297³,𝑎_𝑛+1,𝑛=−73224𝛼297³,𝑏_𝑛+1,𝑛=−73224𝛼297³,𝑎_{𝑛+1,𝑛+1}=−41472𝛼297³,𝑏_{𝑛+1,𝑛+1}=−41472𝛼297³. (16)

Before starting the solution process, initial parameters 𝛿⁰ must be determined by using the initial condition and the following derivatives at the boundaries: 𝑈_𝑁(𝑥,0)=𝑈𝑥_𝑚𝑈,0𝑚=0,1,2,…,𝑁,_𝑁𝑥(𝑎,0)=0,𝑈_𝑁𝑥()𝑈𝑏,0=0,_𝑁𝑥𝑥(𝑎,0)=0,𝑈_𝑁𝑥𝑥()𝑈𝑏,0=0,_{𝑁𝑥𝑥𝑥}(𝑎,0)=0,𝑈_{𝑁𝑥𝑥𝑥}(𝑏,0)=0. (17) So we have the following matrix form for the initial vector 𝛿⁰: 𝑀𝛿⁰=𝑘, (18) where𝑀=15362712768248273181210568.5811047968110063.518196008196597811957688196474⋱8112011120119124161191120111209647481195768819659781960018110063.58110479681210568.58182731,812476827121536 (19)𝛿⁰=(𝛿₀,𝛿₁,𝛿₂,…,𝛿_𝑁−2,𝛿_𝑁−1,𝛿_𝑁)^𝑇,and 𝑘=(𝑈(𝑥₀,0),𝑈(𝑥₁,0),…,𝑈(𝑥_𝑁−1,0),𝑈(𝑥_𝑁,0))^𝑇.

2.1. A Linear Stability Analysis

We have investigated stability analysis by applying the von-Neumann approach in which the growth factor of typical Fourier mode is given by 𝛿^𝑛_𝑗=𝜉^𝑛𝑒^{𝑖𝑗𝑘ℎ},𝑖=√−1, (20) where 𝑘 is a mode number and ℎ is the element size. To apply this method, we have linearized the nonlinear term 𝑈²𝑈_𝑥 by considering 𝑈² as a constant such as 𝑍_𝑖 in (9). If we substitute (20) into the iterative system (13) we obtain the following equation: 𝜉^𝑛+1=𝑔𝜉^𝑛, (21) where 𝑔 is the growth factor. We have identified the collocation points with the nodes and used (8) to evaluate 𝑈_𝑚 and its space derivatives in (3). This leads to a set of ordinary differential equations in the following form: 𝛿̇_𝑚−3+120𝛿̇_𝑚−2+1191𝛿̇_𝑚−1+2416𝛿̇_𝑚+1191𝛿̇_𝑚+1+120𝛿̇_𝑚+2+𝛿̇_𝑚+3+7ℎ1+6𝑍_𝑚×−𝛿_𝑚−3−56𝛿_𝑚−2−245𝛿_𝑚−1+245𝛿_𝑚+1+56𝛿_𝑚+2+𝛿_𝑚+3−42𝜇ℎ²𝛿̇_𝑚−3+24𝛿̇_𝑚−2+15𝛿̇_𝑚−1−80𝛿̇_𝑚+15𝛿̇_𝑚+1+24𝛿̇_𝑚+2+𝛿̇_𝑚+3=0, (22) where 𝑍_𝑚=(𝛿_𝑚−3+120𝛿_𝑚−2+1191𝛿_𝑚−1+2416𝛿_𝑚+1191𝛿_𝑚+1+120𝛿_𝑚+2+𝛿_𝑚+3)². Here ⋅ denotes derivative with respect to time. If the parameters 𝛿_𝑖’s and their time derivatives in (22) are discretized by the Crank-Nicolson formula and usual forward finite difference approximation, respectively, 𝛿_𝑖=𝛿^𝑛+𝛿^𝑛+12,𝛿̇_𝑖=𝛿^𝑛+1−𝛿^𝑛,Δ𝑡 (23) we obtain a recurrence relationship between two time levels 𝑛 and 𝑛+1 relating two unknown parameters 𝛿_𝑖^𝑛+1, 𝛿^𝑛_𝑖 for 𝑖=𝑚−3,𝑚−2,…,𝑚+2,𝑚+3, 𝛾₁𝛿^𝑛+1_𝑚−3+𝛾₂𝛿^𝑛+1_𝑚−2+𝛾₃𝛿^𝑛+1_𝑚−1+𝛾₄𝛿_𝑚^𝑛+1+𝛾₅𝛿^𝑛+1_𝑚+1+𝛾₆𝛿^𝑛+1_𝑚+2+𝛾₇𝛿^𝑛+1_𝑚+3=𝛾₇𝛿^𝑛_𝑚−3+𝛾₆𝛿^𝑛_𝑚−2+𝛾₅𝛿^𝑛_𝑚−1+𝛾₄𝛿^𝑛_𝑚+𝛾₃𝛿^𝑛_𝑚+1+𝛾₂𝛿^𝑛_𝑚+2+𝛾₁𝛿^𝑛_𝑚+3, (24) where 𝛾₁=1−𝐸−𝑀,𝛾₂𝛾=120−56𝐸−24𝑀,₃=1191−245𝐸−15𝑀,𝛾₄𝛾=2416+80𝑀,₅=1191+245𝐸−15𝑀,𝛾₆𝛾=120+56𝐸−24𝑀,₇=1+𝐸−𝑀,𝑚=0,1,…,𝑁,𝐸=7Δ𝑡2ℎ1+6𝑍_𝑚,𝑀=42𝜇ℎ². (25)

Substituting the Fourier mode (20) into (24) leads to the growth factor of the form 𝜉=𝑎−𝑖𝑏,𝑎+𝑖𝑏 (26) where ()[]+()[]+()[],[][][],𝑎=1208+40𝑀+1191−15𝑀cosℎ𝑘120−24𝑀cos2ℎ𝑘1−𝑀cos3ℎ𝑘𝑏=245𝐸sinℎ𝑘+56𝐸sin2ℎ𝑘+𝐸sin3ℎ𝑘 (27)|𝜉|²=1; therefore the linearized scheme is unconditionally stable.

3. Numerical Examples and Results

In this section, we have obtained numerical solution of the MRLW equation for motion of single solitary wave, interaction of two and three solitary waves, and development of the Maxwellian initial condition into solitary waves. Accuracy of the method is measured by using the following error norms: 𝐿₂=𝑈^exact−𝑈_𝑁2≃ℎ^𝑁_𝐽=1𝑈^exact_𝑗−𝑈_𝑁𝑗²,𝐿_∞=𝑈^exact−𝑈_𝑁∞≃max_𝑗𝑈^exact_𝑗−𝑈_𝑁𝑗. (28) The discrete conservation properties of the MRLW equation corresponding to mass, momentum, and energy are determined by finding the following three invariants [41]: 𝐼₁=^𝑏_𝑎𝑈𝑑𝑥≃ℎ^𝑁_𝐽=1𝑈^𝑛_𝑗,𝐼₂=^𝑏_𝑎𝑈²+𝜇𝑈_𝑥²𝑑𝑥≃ℎ^𝑁_𝐽=1𝑈^𝑛_𝑗²+𝜇𝑈_𝑥^𝑛_𝑗,𝐼₃=^𝑏_𝑎𝑈⁴−𝜇𝑈²_𝑥𝑑𝑥≃ℎ^𝑁_𝐽=1𝑈^𝑛_𝑗⁴−𝜇𝑈_𝑥^𝑛_𝑗. (29)

3.1. The Motion of Single Solitary Wave

As a first problem, MRLW equation (3) is considered with the boundary conditions 𝑈→0 as 𝑥→±∞ and the initial condition 𝑈(𝑥,0)=√𝑐secℎ𝑝𝑥−𝑥₀. (30) Single solitary wave solution of the MRLW equation has an analytical solution of the form 𝑈(𝑥,𝑡)=√𝑐secℎ𝑝𝑥−(𝑐+1)𝑡−𝑥₀, (31) where 𝑝=𝑐/𝜇(𝑐+1) and 𝑥₀ and 𝑐 are arbitrary constants. This solution corresponds to motion of single solitary wave with amplitude √𝑐, initially centered at 𝑥₀ and with wave velocity 1+𝜇𝑐. For this problem the analytical values of the invariants are [29] 𝐼₁=𝜋√𝑐𝑝,𝐼₂=2𝑐𝑝+2𝜇𝑝𝑐3,𝐼₃=4𝑐²−3𝑝2𝜇𝑝𝑐3. (32)

For the computational work, two sets of parameters have been chosen and discussed. First of all, we have taken the parameters 𝑐=1,𝜇=1,ℎ=0.2,𝑥₀=40,𝑘=0.025 over the interval [0,100] to compare our results with [6, 29, 30, 35, 36]. Thus, the solitary wave has an amplitude 1.0 and the computations are done up to time 𝑡=10 to obtain the invariants and error norms 𝐿₂ and 𝐿_∞. Values of the three invariants and error norms are reported in Table 1. It is clearly seen from the table that the error norms are satisfactorily small enough and the computed values of invariants are in good agreement with their analytical values 𝐼₁=4.4428829,𝐼₂=3.2998316, and 𝐼₃=1.4142135. The percentage of the relative error of the conserved quantities 𝐼₁, 𝐼₂, and 𝐼₃ is calculated with respect to the conserved quantities at 𝑡=0. Percentage of relative changes of 𝐼₁, 𝐼₂, and 𝐼₃ is found to be 0.001×10⁻³%,3.261×10⁻³%, and 7.609×10⁻³%, respectively. Thus, the quantities 𝐼₁, 𝐼₂, and 𝐼₃ remain constant during the computer run. For this case in Table 2, we compare the values of the invariants and error norms obtained by using the present method and some earlier methods [6, 29, 30, 35, 36]. From the table, we observed that the error norms obtained by the present method are less than those of other methods [6, 29, 30, 35, 36]. Figure 1 shows the motion of solitary wave with 𝑐=1,ℎ=0.2,𝑘=0.025 at different time levels. It is observed that the solitary wave moves to the right with constant velocity and amplitude. At 𝑡=0, the amplitude is 1.0 which is located at 𝑥=40, while it is 0.999950 at 𝑡=10 located at 𝑥=60. The absolute difference in amplitudes at times 𝑡=0 and 𝑡=10 is found to be 5×10⁻⁵, so there is a little change between the amplitudes.

Table 1 

Invariants and error norms for single solitary wave with 𝑐=1, ℎ=0.2, 𝑘=0.025, 0≤𝑥≤100.

Table 2 

Errors and invariants for single solitary wave with 𝑐=1, ℎ=0.2, 𝑘=0.025, 0≤𝑥≤100, at 𝑡=10.

Figure 1 

Single solitary wave with 𝑐=1,ℎ=0.2,Δ𝑡=0.025,0≤𝑥≤100 at times 𝑡=0,2,4,6,8, and 10.

Secondly, we have taken the parameters 𝜇=1,𝑐=0.3,ℎ=0.1,𝑘=0.01 and 𝑥₀=40 with range [0,100] to enable comparison with [6, 30, 32, 34–36]. So the solitary wave has amplitude 0.547723. Simulations are run up to time 𝑡=20. Error norms 𝐿₂ and 𝐿_∞ and conserved quantities are tabulated in Table 3, together with the results obtained in [6, 30, 32, 34–36]. These results show high degree of accuracy and efficiency of the method. The invariants 𝐼₁, 𝐼₂, and 𝐼₃ have changed by less than 0.001×10⁻³%,0.189×10⁻³%, and 1.656×10⁻³% percent, respectively. Moreover, the variation of the invariants 𝐼₂ and 𝐼₃ from the initial variants is less than 2.5×10⁻⁶ and 2.6×10⁻⁶, respectively, whereas the changes in the invariant 𝐼₁ approach to zero throughout the interval. The perspective views of traveling solitary wave at different time levels are shown in Figure 2. The distributions of the errors at time 𝑡=10 and 𝑡=20 are shown graphically for solitary waves amplitudes 1 and 0.3 in Figure 3. It is seen that the maximum errors are close to the tip of the solitary waves and between −6×10⁻³ and 6×10⁻³, −2×10⁻⁴, and 2×10⁻⁴, respectively. The CPU times for an Intel(R) Core i5, 2.53 GHz, are also given in Tables 4 and 5.

Table 3 

Invariants and error norms for single solitary wave with 𝑐=0.3, ℎ=0.1, 𝑘=0.01, 0≤𝑥≤100.

Table 4 

Invariants and error norms for single solitary wave with 𝑐=0.6, ℎ=0.1, 𝑘=0.1, −40≤𝑥≤60.

Table 5 

Invariants and error norms for single solitary wave with 𝑐=0.18, ℎ=0.1, 𝑘=0.1, −80≤𝑥≤120.

Figure 2 

Single solitary wave with 𝑐=0.3,ℎ=0.1,Δ𝑡=0.01,0≤𝑥≤100 at times 𝑡=0,5,10,15, and 20.

(a)

(b)

Figure 3 

Error with (a) 𝑐=1,ℎ=0.2,Δ𝑡=0.025,𝑡=10,0≤𝑥≤100 and (b) 𝑐=0.3,ℎ=0.1,Δ𝑡=0.01,𝑡=20,0≤𝑥≤100.

3.2. Interaction of Two Solitary Waves

As a second problem, interaction of two well separated solitary waves having different amplitudes and traveling in the same direction is considered by using the initial condition 𝑈(𝑥,0)=²_𝑗=1𝐴_𝑗secℎ𝑝_𝑗𝑥−𝑥_𝑗, (33) where 𝐴_𝑗=𝑐_𝑗, 𝑝_𝑗=𝑐_𝑗/𝜇(𝑐_𝑗+1), 𝑗=1,2, and 𝑐_𝑗 and 𝑥_𝑗 are arbitrary constants. The analytical values of the invariants can be found as [29] 𝐼₁=²_𝑗=1𝜋𝑐_𝑗𝑝_𝑗𝐼=11.467698,₂=²_𝑗=12𝑐_𝑗𝑝_𝑗+2𝜇𝑝_𝑗𝑐_𝑗3𝐼=14.629243,₃=²_𝑗=14𝑐²_𝑗3𝑝_𝑗−2𝜇𝑝_𝑗𝑐_𝑗3=22.880466. (34)

We have studied the interaction of two positive solitary waves having the parameters 𝜇=1,ℎ=0.2,𝑘=0.025,𝑐₁=4,𝑐₂=1,𝑥₁=25,𝑥₂=55 through the interval 0≤𝑥≤250 to coincide with those used by [6, 30, 34–36]. The simulations are maintained up to 𝑡=20. Constant values 𝐼₁,𝐼₂, and 𝐼₃ at various time steps together with equivalent results for the previous methods are shown in Table 6. It is seen that the numerical values of the invariants remain almost constant during the computer run. The interaction of two solitary waves is shown in Figure 4. It can be seen from the figure that, at 𝑡=0, the wave with larger amplitude is to the left of the second wave with smaller amplitude. Since the taller wave moves faster than the shorter one, it catches up and collides with the shorter one at 𝑡=8 and then moves away from the shorter one as time increases. At 𝑡=20, the amplitude of larger waves is 2.001102 at the point 𝑥=127.4, whereas the amplitude of the smaller one is 0.996403 at the point 𝑥=92. It is found that the absolute difference in amplitude is 3.59×10⁻³ for the smaller wave and 1.10×10⁻³ for the larger wave for this case.

Table 6 

Comparison of the invariants for the interaction of two solitary waves with results from [30] with ℎ=0.2, 𝑘=0.025 in the region 0≤𝑥≤250.

Figure 4 

Interaction of two solitary waves with 𝑡=0,4,8,10,14, and  20.

3.3. Interaction of Three Solitary Waves

As a third problem, interaction of three solitary waves having different amplitudes and travelling in the same direction is studied. We consider (3) with initial conditions given by 𝑈(𝑥,0)=³_𝑗=1𝐴_𝑗secℎ𝑝_𝑗𝑥−𝑥_𝑗, (35) where 𝐴_𝑗=𝑐_𝑗, 𝑝_𝑗=𝑐_𝑗/𝜇(𝑐_𝑗+1), 𝑗=1,2,3, 𝑐_𝑗 and 𝑥_𝑗 are arbitrary constants. For this problem the analytical values of the invariants are found from (32) as 𝐼₁=³_𝑗=1𝜋𝑐_𝑗𝑝_𝑗𝐼=14.9801,₂=³_𝑗=12𝑐_𝑗𝑝_𝑗+2𝜇𝑝_𝑗𝑐_𝑗3𝐼=15.8218,₃=³_𝑗=14𝑐²_𝑗3𝑝_𝑗−2𝜇𝑝_𝑗𝑐_𝑗3=22.9923. (36)

In order to be able to compare with the previous works, computations are performed for the parameters 𝜇=1,ℎ=0.2,𝑘=0.025,𝑐₁=4,𝑐₂=1,𝑐₃=0.25,𝑥₁=15,𝑥₂=45,𝑥₃=60 over the interval 0≤𝑥≤250. The simulation is run from 𝑡=0 to 𝑡=45. Table 7 shows a comparison of the values of the invariants obtained by the present method with those obtained in [30, 34–36]. It is evident from the table that the obtained values of the invariants remain almost constant during the computer run which are all in good agreement with their analytical values given by (36). Figure 5 depicts the interaction of three solitary waves at different times. It is observed from Figure 5 that interaction started about time 𝑡=10, overlapping processes occurred between time 𝑡=15 and 𝑡=40, and waves started to resume their original shapes after the time 𝑡=40.

Table 7 

Comparison of invariants for the interaction of three solitary waves with results from [30] with ℎ=0.2, 𝑘=0.025 in the region 0≤𝑥≤250.

Figure 5 

Interaction of three solitary waves with 𝑡=0,5,8,15,20,and40.

3.4. The Maxwellian Initial Condition

As our last problem, the development of the Maxwellian initial condition, 𝑈(𝑥,0)=exp−(𝑥−40)², (37) into a train of solitary waves is considered. For the Maxwellian initial condition, behavior of the solution depends on the values of 𝜇. So we take 𝜇=0.1,𝜇=0.04,𝜇=0.015, and 𝜇=0.01. The numerical values of the invariants quantities during the simulations are given in Table 8. For 𝜇=0.1, only a single soliton occurred as depicted in Figure 6(a). When 𝜇=0.015 and 𝜇=0.01, two and three stable solitons occurred, respectively, as depicted in Figures 6(b) and 6(c). For 𝜇=0.04, Maxwellian initial condition has decayed into four solitary waves as depicted in Figure 6(d). As is seen, when 𝜇 is reduced, more solitary waves occurred. All figures are drawn up at time 𝑡=14.5. The peaks of the well developed wave lie on a straight line so that their velocities are linearly dependent on their amplitudes and also we observe a small oscillating tail appearing behind the last wave in all Maxwellian figures.

Table 8 

Invariants of MRLW equation using the Maxwellian initial condition.

(a)

(b)

(c)

(d)

Figure 6 

Maxwellian initial condition at 𝑡=14.5 with (a) 𝜇=0.1, (b) 𝜇=0.015, (c) 𝜇=0.01, and (d) 𝜇=0.04.

4. Conclusion

In this paper, a numerical treatment of the MRLW equation has been introduced using septic B-spline collocation finite element method. The nonlinear term in the equation has been linearized by using a form given in the paper [37]. To examine the performance of the scheme, four test problems have been studied. The performance and accuracy of the method have been tested by calculating the error norms 𝐿₂ and 𝐿_∞ and the invariant quantities 𝐼₁,𝐼₂, and 𝐼₃. Linear stability analysis proved that the present scheme is unconditionally stable. The experimental results of the algorithm are much satisfactory in comparison with the previous results. Our method can successfully be used to model the motion and interaction of the solitary waves. Thus, we can assert that our scheme is efficient and reliable for obtaining the numerical solutions of the other physically important nonlinear partial differential equations.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

The authors wish to thank the reviewers for their careful reading and making several useful comments improving the presented paper.