Numerical approximation of the generalized regularized long wave equation using Petrov-Galerkin finite element method

TL;DR

This paper develops a Petrov-Galerkin finite element method using cubic shape and quadratic B-spline weight functions to accurately and stably approximate solutions of the nonlinear generalized regularized long wave equation, validated through multiple test problems.

Contribution

It introduces a novel Petrov-Galerkin scheme with specific spline functions for the GRLW equation, providing stability analysis and demonstrating high accuracy in numerical simulations.

Findings

The scheme is unconditionally stable based on Fourier analysis.

Numerical results show high accuracy compared to existing methods.

The method effectively captures solitary wave interactions and evolution.

Abstract

The generalized regularized long wave (GRLW) equation has been developed to model a variety of physical phenomena such as ion-acoustic and magnetohydrodynamic waves in plasma, nonlinear transverse waves in shallow water and phonon packets in nonlinear crystals. This paper aims to develop and analyze a powerful numerical scheme for the nonlinear generalized regularized long wave (GRLW) equation by Petrov--Galerkin method in which the element shape functions are cubic and weight functions are quadratic B-splines. The suggested method is performed to three test problems involving propagation of the single solitary wave, interaction of two solitary waves and evolution of solitons with the Maxwellian initial condition. The variational formulation and semi-discrete Galerkin scheme of the equation are firstly constituted. We estimate accuracy of such a spatial approximation. Then Fourier stability analysis of the linearized scheme shows that it is unconditionally stable. To verify practicality and robustness of the new scheme error norms $L_{2}$, $L_{\infty }$ and three invariants $I_{1},I_{2}$ and $I_{3}$ are calculated. The obtained numerical results are compared with other published results and shown to be precise and effective.

Tables

Table 1. Table 1: Invariants and errors for single solitary wave with p = 2 , c = 1 , h = 0.2 , Δ ​ t = 0.025 , μ = 1 , x ∈ [ 0 , 100 ] . formulae-sequence 𝑝 2 formulae-sequence 𝑐 1 formulae-sequence ℎ 0.2 formulae-sequence Δ 𝑡 0.025 formulae-sequence 𝜇 1 𝑥 0 100 p=2,c=1,h=0.2,\Delta t=0.025,\mu=1,x\in\left[0,100\right].

$Time$	$I_1$	$I_2$	$I_3$	$L_2\times {10}^3$	$L_{\infty }\times {10}^3$
0	4.442866	3.299813	1.414214	0.000000	0.000000
2	4.442940	3.299938	1.414330	1.948707	1.190456
4	4.443005	3.300033	1.414425	2.362855	1.222540
6	4.443068	3.300124	1.414515	2.449792	1.198936
8	4.443129	3.300213	1.414604	2.448242	1.150862
10	4.443175	3.300302	1.414692	2.415468	1.079686

Table 2. Table 2: Comparisons of results for single solitary wave with p = 2 , 𝑝 2 p=2, c = 1 , 𝑐 1 c=1, h = 0.2 , ℎ 0.2 h=0.2, Δ ​ t = 0.025 , Δ 𝑡 0.025 \Delta t=0.025, μ = 1 , 𝜇 1 \mu=1, x ∈ [ 0 , 100 ] . 𝑥 0 100 x\in\left[0,100\right].

$Method$	$I_1$	$I_2$	$I_3$	$L_2\times {10}^3$	$L_{\infty }\times {10}^3$
Analytic	4.44288	3.29983	1.41421	0.000000	0.000000
Our Method	4.443175	3.300302	1.414692	2.415468	1.079686
Petrov-Galerkin[32]	4.44288	3.29981	1.41416	3.00533	1.68749
Septic Collocation First Scheme[34]	4.442866	3.299822	1.414204	2.632463	1.393064
Septic Collocation Second Scheme[34]	4.442866	3.299715	1.414312	2.571481	1.340210
Cubic Galerkin[35]	3.801670	2.888066	0.979294	13.291080	8.478107
Cubic B-spline coll-CN[39]	4.442	3.299	1.413	16.39	9.24
Cubic B-spline coll + PA-CN[39]	4.440	3.296	1.411	20.3	11.2
Cubic B-spline collocation[40]	4.44288	3.29983	1.41420	9.30196	5.43718

Table 3. Table 3: Invariants and errors for single solitary wave with p = 3 , 𝑝 3 p=3, c = 6 / 5 , 𝑐 6 5 c=6/5, h = 0.1 , ℎ 0.1 h=0.1, Δ ​ t = 0.025 , Δ 𝑡 0.025 \Delta t=0.025, μ = 1 , 𝜇 1 \mu=1, x ∈ [ 0 , 100 ] . 𝑥 0 100 x\in\left[0,100\right].

$Time$	$I_1$	$I_2$	$I_3$	$L_2\times {10}^3$	$L_{\infty }\times {10}^3$
0	3.797185	2.881250	0.972968	0.000000	0.000000
2	3.797187	2.881258	0.973414	1.700682	1.174285
4	3.797187	2.881257	0.973473	2.805942	1.797229
6	3.797187	2.881255	0.973486	3.899300	2.428864
8	3.797200	2.881254	0.973487	5.007404	3.073644
10	3.797282	2.881293	0.973446	6.128029	3.722138

Table 4. Table 4: Comparisons of results for single solitary wave with p = 3 , 𝑝 3 p=3, c = 6 / 5 , 𝑐 6 5 c=6/5, h = 0.1 , ℎ 0.1 h=0.1, Δ ​ t = 0.025 , Δ 𝑡 0.025 \Delta t=0.025, μ = 1 , 𝜇 1 \mu=1, x ∈ [ 0 , 100 ] . 𝑥 0 100 x\in\left[0,100\right].

$Method$	$I_1$	$I_2$	$I_3$	$L_2\times {10}^3$	$L_{\infty }\times {10}^3$
Our Method	3.797282	2.881293	0..973446	6.128029	3.722138
Petrov-Galerkin[32]	3.79713	2.88123	0.972243	7.76745	4.70875
Septic Collocation First Scheme[34]	3.797185	2.881252	0.973145	8.972983	5.175982
Septic Collocation Second Scheme[34]	3.797133	2.881089	0.973128	7.778169	4.441873
Cubic Galerkin[35]	3.801670	2.888066	0.979294	13.291080	8.478107

Table 5. Table 5: Invariants and errors for single solitary wave with p = 4 , 𝑝 4 p=4, c = 4 / 3 , 𝑐 4 3 c=4/3, h = 0.1 , ℎ 0.1 h=0.1, Δ ​ t = 0.01 , Δ 𝑡 0.01 \Delta t=0.01, μ = 1 , 𝜇 1 \mu=1, x ∈ [ 0 , 100 ] . 𝑥 0 100 x\in\left[0,100\right].

$Time$	$I_1$	$I_2$	$I_3$	$L_2\times {10}^3$	$L_{\infty }\times {10}^3$
0	3.468709	2.671691	0.729204	0.000000	0.000000
2	3.468718	2.671714	0..729969	0.967786	0.708600
4	3.468719	2.671714	0.730017	1.040242	0.591250
6	3.468720	2.671714	0.730027	1.102854	0.611363
8	3.468731	2.671714	0.730028	1.183442	0.715175
10	3.468799	2.671742	0.730001	1.283420	0.821650

Table 6. Table 6: Comparisons of results for single solitary wave with p = 4 , 𝑝 4 p=4, c = 4 / 3 , 𝑐 4 3 c=4/3, h = 0.1 , ℎ 0.1 h=0.1, Δ ​ t = 0.01 , Δ 𝑡 0.01 \Delta t=0.01, μ = 1 , 𝜇 1 \mu=1, x ∈ [ 0 , 100 ] . 𝑥 0 100 x\in\left[0,100\right].

$Method$	$I_1$	$I_2$	$I_3$	$L_2\times {10}^3$	$L_{\infty }\times {10}^3$
Our Method	3.468799	2.671742	0..730001	1.283420	0.821650
Petrov-Galerkin[32]	3.46866	2.67168	0.728881	2.46065	1.56620
Septic Collocation First Scheme[34]	3.468709	2.671696	0.729258	3.351740	2.049733
Septic Collocation Second Scheme[34]	3.468671	2.671658	0.729237	2.698709	1.656002
Cubic Galerkin[35]	3.470439	2.674445	0.731987	1.511394	0.857585

Table 7. Table 7: Invariants for interaction of two solitary waves with p = 3 . 𝑝 3 p=3.

	$t$	$0$	$2$	$4$	$6$
	Our Method	9.690777	9.690777	9.690777	9.690777
	[32]	9.69075	9.69074	9.69074	9.69074
$I_1$	[34] First	9.690777	9.690777	9.690777	9.690778
	[34] Second	9.690777	9.688117	9.686015	9.683462
	[35]	9.6907	9.6906	9.6898	9.6901
	Our Method	12.944360	12.928161	12.957476	12.988509
	[32]	12.9444	12.9452	12.9453	12.9454
$I_2$	[34] First	12.944391	12.944392	12.944393	12.944394
	[34] Second	12.944391	12.939062	12.970312	13.002753
	[35]	12.9443	12.9440	12.9418	12.9426
	Our Method	17.018706	17.034905	17.005590	16.974557
	[32]	17.0184	16.9835	16.9261	16.9113
$I_3$	[34] First	17.018675	17.02567	16.981696	16.952024
	[34] Second	17.018675	17.02400	16.992754	16.960313
	[35]	17.0187	17.0324	16.9849	16.9557

Table 8. Table 8: Invariants for interaction of two solitary waves with p = 4 . 𝑝 4 p=4.

	$t$	$0$	$2$	$4$	$6$
	Our Method	8.834272	8.834272	8.834272	8.834272
	[32]	8.83427	8.84204	8.84209	8.83434
$I_1$	[34] First	8.834272	8.834160	8.834053	8.8339467
	[34] Second	8.834272	8.564186	8.435464	8.327161
	[35]	8.8342	8.7089	8.6518	8.6134
	Our Method	12.170706	11.339311	11.209384	15.812521
	[32]	12.1697	12.3700	12.5703	12.6103
$I_2$	[34] First	12.170887	12.170537	12.170205	12.169873
	[34] Second	12.170887	11.939598	11.977097	11.814722
	[35]	12.1707	11.7871	11.6179	11.4992
	Our Method	14.029604	14.860999	14.990927	10.387789
	[32]	14.0302	13.9607	13.9805	14.6974
$I_3$	[34] First	14.029423	14.413442	14.351624	14.292901
	[34] Second	14.029423	14.260712	14.223214	14.385588
	[35]	14.0296	12.9204	12.1972	11.9640

Table 9. Table 9: Maxwellian initial condition for different values of μ . 𝜇 \mu.

$\mu$	$t$	p=2			p=3			p=4
		$I_1$	$I_2$	$I_3$	$I_1$	$I_2$	$I_3$	$I_1$	$I_2$	$I_3$
	0.01	1.772481	1.378659	0.760911	1.772481	1.378655	0.760779	1.772422	1.378551	0.760310
0.1	0.03	1.772475	1.378639	0.760890	1.772435	1.378538	0.759621	1.772189	1.378049	0.755954
	0.05	1.772463	1.378599	0.760847	1.772328	1.378278	0.757292	1.771793	1.377149	0.747477
	0.01	1.772480	1.315994	0.823572	1.772480	1.315988	0.823384	1.772382	1.315812	0.822659
0.05	0.03	1.772469	1.315959	0.823526	1.772388	1.315772	0.821657	1.771994	1.314976	0.816132
	0.05	1.772445	1.315887	0.823429	1.772171	1.315282	0.818116	1.71289	1.313372	0.803134
	0.01	1.772480	1.284661	0.854901	1.772478	1.284651	0.854699	1.772340	1.284399	0.853806
0.025	0.03	1.772462	1.284610	0.854823	1.772317	1.284299	0.852665	1.771760	1.283154	0.846013
	0.05	1.772424	1.284502	0.854658	1.771933	1.283467	0.848301	1.770623	1.280410	0.829240

Equations

$$u_t+u_x+au{u}_x-b{u}_{xxt}=0,$$

$$u_t+u_x+p(p+1)u^p{u}_x-\mu u_{xxt}=0,$$

$$u(x,0)=f_1(x),a\le x\le b,$$

$$\begin{array}[]{lll}u(a,t)=0,~{}~{}~{}~{}~{}u(b,t)=0,&&\\ u_{x}(a,t)=0,~{}~{}~{}~{}~{}u_{x}(b,t)=0,&&\\ u_{xx}(a,t)=0,~{}~{}~{}~{}~{}u_{xx}(b,t)=0,&t>0&\end{array}$$

$$u_t-\mu \Delta u_t=\nabla \mathcal{F}(u),$$

$$u(x,0)=f_1(x),a\le x\le b,$$

$$\begin{array}[]{lll}u(a,t)=0,~{}~{}~{}~{}~{}u(b,t)=0,&&\\ u_{x}(a,t)=0,~{}~{}~{}~{}~{}u_{x}(b,t)=0,&&\\ u_{xx}(a,t)=0,~{}~{}~{}~{}~{}u_{xx}(b,t)=0,&t>0.&\end{array}$$

$$H_0^k(\Omega )=\left\{v\in H^k(\Omega ):D^{i}v=0\text{on }\partial \Omega ,i=0,1,\cdots ,k-1\right\}$$

$$(u_t,\xi )-\mu (\Delta u_t,\xi )=(\nabla \mathcal{F}(u),\xi ).$$

$$\left(u_t,\xi \right)+\mu \left(\nabla u_t,\nabla \xi \right)=-\left(\mathcal{F}(u),\nabla \xi \right),\forall \xi \in H_0^1,$$

$$\Vert u(t){\Vert }_1=\Vert u_0{\Vert }_1,t\in (0,T],\text{and }\Vert u{\Vert }_{L^{\infty }(L^{\infty }(\Omega ))}\le C\Vert u_0{\Vert }_1$$

$$\left(u_t,u\right)+\mu \left(\nabla u_t,\nabla u\right)=-\left(\mathcal{F}(u),\nabla u\right)$$

$$\frac{1}{2}\frac{d}{dt}\left[\Vert u{\Vert }^2+\mu \Vert \nabla u{\Vert }^2\right]={\int }_{\Omega }u[\nabla \cdot \mathcal{F}(u)]dx.$$

$$u\nabla \cdot \mathcal{F}(u)=\nabla \cdot [\mathcal{F}(u)u]-\nabla \cdot [\mathcal{G}(u)],$$

$${\int }_{\Omega }u[\nabla \cdot \mathcal{F}(u)]dx={\int }_{\Omega }\nabla (u\mathcal{F}(u))dx=0.$$

$$\frac{1}{2}\frac{d}{dt}\left(\Vert u{\Vert }_1^2\right)=0,$$

$$\Vert u{\Vert }_1^2=\Vert u_0{\Vert }_1^2,$$

$$u\in L^{\infty }(0,T,H_0^1(\Omega ))\text{with }(u(x,0),\xi )=(u_0,\xi ),\xi \in H_0^1(\Omega ),$$

$$v^m=span\{w_1,w_2,\cdots ,w_m\}.$$

$$u^m(t)=\sum _{i=1}^m{c}_i(t)w_i,$$

$$\left(u_t^m,\xi \right)+\left(\nabla u_t^m,\nabla \xi \right)=-\left(\mathcal{F}(u^m),\nabla \xi \right),\forall \xi \in v^m,$$

$$u_{0,m}=u^m(0)=\sum _{i=1}^m{c}_i(0)w_i=P^m{u}_0.$$

$$\|u^{m}\|_{\infty}\leq C\|u_{0}\|_{1}\$$

$$\Vert \mathcal{F}(u^m){\Vert }^2\le C\Vert u_0{\Vert }_1^2$$

$$\left(u_t^m,u_t^m\right)+\left(\nabla u_t^m,\nabla u_t^m\right)=-\left(\mathcal{F}(u^m),\nabla u_t^m\right).$$

$$\Vert u_t^m{\Vert }_1^2=-\left(\mathcal{F}(u^m),\nabla u_t^m\right)$$

$$\Vert u_t^m{\Vert }_1\le C\Vert u_0{\Vert }_1.$$

$$\left(u_t^m,w_i\right)+\left(\nabla u_t^m,\nabla w_i\right)=-\left(\mathcal{F}(u^m),\nabla w_i\right).$$

$$\left(W_t^m,\xi \right)+\left(\nabla W_t^m,\nabla \xi \right)=-\left(\mathcal{F}(W^m),\nabla \xi \right).$$

$$\frac{d}{dt}\Vert W{\Vert }_1\le C\Vert W{\Vert }_1.$$

Full text

Keywords: GRLW equation; Petrov-Galerkin; Cubic B-splines; Solitary waves; Soliton.

AMS classification: 65N30, 65D07, 74S05,74J35, 76B25.

1 Introduction

The GRLW equation was originated by a famous nonlinear analyst Peregrine who first successfully introduced the regularized long wave equation as a perfect alternative to the famous KdV equation for studying soliton phenomena and as a mathematical model for small amplitude long waves on the surface of water  [1]. Nonlinear evolution equations play fundemental roles in various fields of science mostly in physics, applied mathematics and in engineering problems. Analytical solutions of these equations are commonly not derivable, particularly when the nonlinear terms are contained. Numerical solutions of these equations are very practical to analyze the physical phenomena due to the fact that analytical solutions of these equations are found for the restricted boundary and initial conditions. The regularized long wave (RLW) equation

[TABLE]

is one of the important model in physics media on account of it defines phenomena with weak nonlinearity and dispersion waves, involving nonlinear transverse waves in shallow water, ion-acoustic waves in plasma, hydromagnetic wave in cold plasma, plasma, elastic media, optical fibres, acoustic-gravity waves in compressible fluids, pressure waves in liquid–gas bubbles and acoustic waves inharmonic crystals. The solutions of this equation are sorts of solitary waves called as solitons whose form are not affected by a collision. It was first alleged by Peregrine [1, 2] for studying soliton phenomena and as a sample for small-amplitude long-waves on the surface of water in a channel and widely studied by Benjamin et al. [3]. In physical situations such as unidirectional waves propagating in a water channel, long-crested waves in near-shore zones, and many others, the RLW equation serves as an alternative model to the Korteweg–de Vries equation (KdV equation) [4, 5]. An exact solution of the equation was obtained under the limited initial and boundary conditions in [6] for this reason it got fascinate from a numerical point of view. Therefore, numerical solutions of the RLW equation have been the matter of many papers. Various effective methods have been presented to solve the equation such as finite difference method [7–10], pseudo-spectral method [11] , meshfree method [12], Adomian decomposition method [13] and various forms of finite element methods in [14–24] . Indeed, the RLW equation is a special case of the generalized regularized long wave (GRLW) equation which is an alternative to the KdV equation for describing nonlinear dispersive waves and can be used to characterise phenomena with weak nonlinearity and dispersion waves. It is defined as

[TABLE]

subject to the initial condition

[TABLE]

and the boundary conditions

[TABLE]

where $p$ is a positive integer, $\mu$ is positive constant and physical boundary conditions require $u$ and $u_{x}$ $\rightarrow 0$ that $u\rightarrow 0$ for $x\rightarrow\pm\infty$. In equation $(\ref{grlw})$ $u$ indicates dimensionless surface elevation, $x$ distance and $t$ time. On the other hand, the GRLW equation has received much less attention, presumably because of its higher nonlinearity for $p>2$ and the fact that it possesses a finite number of conserved quantities and admits solitary waves as solutions, but, unlike other equations, the stability of its solutions depends on their velocity [25]. Some solitary wave solutions for GRLW equations have been obtained by Hamdi et al. [26] and Ramos [27] studied solitary wave interactions based on the separation of the temporal and spatial derivatives. Zhang [28] implemented finite difference method for a Cauchy problem while Kaya [29], Kaya and El-Sayed [30] indicated the numerical solution of the GRLW equation by using the Adomian decomposition method. Roshan [32] have procured numerical solutions of the GRLW equation by the application of Petrov–Galerkin method, which uses a linear hat function as the trial function and a quintic B-spline function as the test function. Wang et al. [33] offered a mesh-free method for the GRLW equation based on the moving least-squares approximation. Karakoç [34] and Zeybek [35] have obtained solitary-wave solutions of the GRLW equation by using septic B-spline collocation and cubic B-spline lumped Galerkin method.

Numerical solutions of the GRLW    equation have been obtained by Soliman [36] using He’s variational iteration method. Mokhtari and Mohammadi [37] suggested the Sinc-collocation method for this equation. A time-linearization method that uses a Crank–Nicolson procedure in time and three point, fourth-order accurate in space, compact difference equations, is presented and used to determine the solutions of the GRLW equation and a modified version thereof (mGRLW) by C.M. García-López and J. I. Ramos [38]. The another special case of the GRLW equation is the modified regularized long wave (MRLW) equation for $p=2$. MRLW equation was solved numerically by various methods [39, 40, 41, 42, 43, 44, 45].

Spline functions are a class of piecewise polynomials which provide continuity features being subject to the degree of the polynomials. They are spectacular mathematical instrument for numerical approximations because of their numerous popular specialities. One kind of splines, noted as B-splines, has been used in obtaining the numerical solution of the GRLW equation [34, 35, 40, 41]. Assemblies of B-splines are used as trial functions in the Petrov-Galerkin methods. Especially, cubic B-splines associated with finite element methods have been verified to give very smooth solutions, and use of the cubic B-splines as shape functions in the finite element method warranties continuity of the first and second-order derivatives of trial solutions at the mesh points [19].

In this study, we have designed a lumped Petrov-Galerkin method for the GRLW equation using cubic B-spline function as element shape function and quadratic B-spline function as the weight function. The plan of this paper is as follows:

• •

In Section 2, the governing equation and its variational formulation and newly established theorems are presented.

• •

A semi-discrete Galerkin scheme of the equation is notified in Section 3.

• •

In Section 4, a lumped Petrov-Galerkin finite element technique has been practiced to GRLW equation. Resulting system can be solved with a sort of the Thomas algorithm.

• •

Section 5, is dedicated to stability analysis of the method.

• •

The results of numerical examples are reported in Section 6. The last section is a brief conclusion.

2 Variational formulation and energy estimates

Here we are dedicated to write the initial-boundary value problem in a variational form, and use this weak form to derive some estimates for its solution. We start by proving existance and uniqueness of solutions by using this variational form. The higher order nonlinear initial boundary value problem (2) can be written as

[TABLE]

where $\mathcal{F}(u)=-u(1+pu^{p}),$ subject to the initial condition

[TABLE]

and the boundary conditions

[TABLE]

To define the weak form of the solutions of (5) and to investigate the existence and uniqueness of solutions of the weak form we define the following spaces.

Here $H^{k}(\Omega)$, $k\geq 0$ (integer) is an usual normed space of real valued functions on $\Omega$ and

[TABLE]

where $D=\frac{\partial}{\partial x}$. We denote the norm on this space by $\|\cdot\|_{k}$ which is the usual $H^{k}$ norm, and $k=0$ $\|\cdot\|_{0}=\|\cdot\|$ represents $L_{2}$ norm and $(\cdot,\cdot)$ represents $L_{2}$ inner product [31].

Multiplying (5) by $\xi\in H_{0}^{1}(\Omega)$, and then integrating over $\Omega$ we have

[TABLE]

Applying Green’s formula on the above inner products we aim to find $u(\cdot,t)\in H_{0}^{1}$ so that

[TABLE]

with $u(0)=u_{0}$.

Theorem 1**.**

If $u$ is a solution of (8) then

[TABLE]

holds if $u_{0}\in H_{0}^{1}$, and $C$ is a positive constant.

Proof.

Replacing $\xi\in H_{0}^{1}$ by $u\in H_{0}^{1}$ in (8) results

[TABLE]

with $u(0)=u_{0}$ which gives

[TABLE]

Now

[TABLE]

if $u\in H_{0}^{1}$ where $\mathcal{G}^{\prime}(u)=\mathcal{F}(u)$. For the simplicity of the analysis from now on in this section we fix $\mu=1$. The analysis for a general $\mu$ is the similar. Also, from the initial conditions in (9) we have $u=0$ on $\partial\Omega$ and so $\mathcal{G}(0)=0$, and then

[TABLE]

Thus (10) takes the form

[TABLE]

and so

[TABLE]

completes the proof of the first part. The second part follows from Sobolev embedding theorem  [31, 48, 47]. ∎

Theorem 2**.**

There exists a unique solution of (8) for any $T>0$ such that

[TABLE]

if $u_{0}\in H_{0}^{1}$ for any $T>0$.

Proof.

In order to prove the uniqueness of the solution of (8) we consider an orthogonal basis $\{w_{i}\}_{i=1}^{\infty}$ for $H_{0}^{1}(\Omega)$ and

[TABLE]

Now we define

[TABLE]

for each $t>0$ to satisfy

[TABLE]

with $u^{m}(0)=u_{0,m}$ where

[TABLE]

Here $P^{m}$ is an orthogonal projection onto finite dimensional space $v^{m}$, and $u_{0,m}\rightarrow u_{0}\in H_{0}^{1}(\Omega)$ [31, 47]. Hence the weak form (11) can be written as a system of first order nonlinear ordinary differential equation and there exist a positive time $t_{m}>0$ such that the nonlinear system of differential equations has a unique solution $u^{m}$ over $(0,\ t_{m})$.

Also from Theorem 1 it is easy to see that

[TABLE]

and

[TABLE]

which shows that $\mathcal{F}(u^{m})$ is bounded in $L^{\infty}(0,T,L_{2}(\Omega))$. Now by setting $\xi=u_{t}^{m}$ in (11)

[TABLE]

Thus

[TABLE]

which yields

[TABLE]

Hence $\{u^{m}\}$ and $\{u_{t}^{m}\}$ are uniformly bounded in $L^{\infty}(0,T,H_{0}^{1}(\Omega))$.

By setting $\xi=w_{i}$ in (11) we have

[TABLE]

Thus the existence of solutions of the problem follows from the denseness of $\{w_{i}\}$ in $H_{0}^{1}(\Omega)$.

Considering $u$ and $v$ as two solutions of (8) with $u(0)=0$ and $v(0)=0$, we define $W=u-v$. Then $W(0)=0$. Also

[TABLE]

Replacing $\xi$ by $W$ in the above equation and following the boundedness of $u$ and $v$ one obtains [31, 48]

[TABLE]

Integrating the above inequality over $[0,\ t]$ yields

[TABLE]

Now applying Gronwall’s Lemma it is easy to see that

[TABLE]

which confirms $W=0$ completes the proof [31, 48]. ∎

3 The semidiscrete Galerkin B-spline finite element method

Consider $0<h<1$. A finite dimensional subspace $S_{h}$ of $H_{0}^{1}(\Omega)$ is considered such that for $u\in H_{0}^{1}(\Omega)\cap H^{4}(\Omega)$, there exists a constant $C$ independent of $h$ [31, 47, 48] such that

[TABLE]

Here we aim to find solutions of a semi-discrete finite element formulation of (5) $u_{h}:[0,\ T]\rightarrow S_{h}$ such that

[TABLE]

with $u_{h}(0)=u_{0,h}\in S_{h}$ is an approximation of $u_{0}$. Before establishing the original convergence result we first prove a priori bound of the solution of (13) below.

Theorem 3**.**

The solution $u_{h}\in S_{h}$ of (13) satisfies

[TABLE]

*and ***

[TABLE]

holds where $C$ is a positive constant.

Proof.

The proof is trivial from our discussion in the previous section (Theorem 1). ∎

Now we move onto establish the theoretical bound of the error in the semi-discrete scheme (13) of (8).

To that end we consider the following bilinear form

[TABLE]

which satisfies the boundedness property

[TABLE]

and coercivity property (on $\Omega$)

[TABLE]

Here $A$ satisfies

[TABLE]

where $\tilde{u}$ is an auxiliary projection of $u$ [31, 47, 48]. Now the accuracy result in such a semi-discrete approximation (13) of (8) can be established by the following theorem.

Theorem 4**.**

Let $u_{h}\in S_{h}$ be a solution of (13) and $u\in H_{0}^{1}(\Omega)$ be that of (8), then the following inequality holds

[TABLE]

where $C>0$ if $\|u(0)-u_{0,h}\|\leq Ch^{4}$ holds.

Proof.

Considering $e=u-u_{h}$ we write

[TABLE]

Here

[TABLE]

Also It follows from (14) and (16) and [48] that

[TABLE]

So (12) and (17) confirms the following inequalities

[TABLE]

Applying $\frac{\partial}{\partial t}$ on (16) and having some simplifications yields [48]

[TABLE]

Also subtracting (13) from (8) it is easy to see that

[TABLE]

Now substituting $\xi=\theta$ in (18), and then applying Cauchy-Schwarz inequality one gets

[TABLE]

Here

[TABLE]

comes from Lipschitz conditions and boundedness of $u$ and $u_{h}$ and thus

[TABLE]

So

[TABLE]

Hence Gronwall’s lemma, bounds of $\nu$ and $\nu_{t}$ confirms

[TABLE]

if $\theta(0)=0$, completes the proof [48, 47]. ∎

4 Numerical implementations of the scheme

Let us consider the solution domain is limited to a finite interval $a\leq x\leq b.$ Partition the interval $[a,b]$ at points by $x_{m}$ where $a=x_{0}<x_{1}<...<x_{N}=b$ and let $h=\frac{b-a}{N},$ $m=0,1,2,...,N$. On this partition, we shall need the following cubic B-splines $\phi_{m}(x)$ at the points $x_{m},$ $m=0,1,2,...,N.$ The cubic B-spline functions $\phi_{m}(x)$, (m= $-1(1)$ $N+1$) are identified as follows [46]

[TABLE]

We search the approximation solution $u_{N}(x,t)$ to the exact solution $u(x,t)$ which uses these cubic B-splines as trial functions

[TABLE]

where $\delta_{j}(t)$ are time depended quantities or the nodal parameters to be detected from boundary and weighted residual conditions. Applying the following transformation

[TABLE]

the finite interval $[x_{m},x_{m+1}]$ is turned into more easily practicable interval $[0,1]$. Therefore cubic B-spline shape functions $(\ref{3})$ depending on variable $\eta$ on the region $[0,1]$ rearranged with

[TABLE]

All splines, apart from $\phi_{m-1}(x),\phi_{m}(x),\phi_{m+1}(x)$,$\phi_{m+2}(x)$ and their four principal derivatives are null over the region $[0,1].$ Thereby variation of the function $u(x,t)$ over element $[0,1]$ is approximated by

[TABLE]

where $\delta_{m-1},\delta_{m},\delta_{m+1},\delta_{m+2}$ and B-spline element functions $\phi_{m-1},\phi_{m},\phi_{m+1},\phi_{m+2}$ as element shape functions. The nodal values $u$ and its derivatives up to second order at the knots $x_{m}$ are given in terms of the parameters $\delta_{m}$ from the use of the B-splines $(\ref{5})$ and and the trial solution $(\ref{6})$:

[TABLE]

We take the weight functions $\Phi_{m}$ as quadratic B-splines. The quadratic B-splines $\Phi_{m}$ at the knots $x_{m}$ are defined as [46]:

[TABLE]

When we take into consideration of the transformation $(\ref{250})$, quadratic B-splines $\Phi_{m}$ are written as

[TABLE]

Performing the Petrov-Galerkin method to Eq.$(\ref{grlw}),$ the weak form of Eq.$(\ref{grlw})$ is attained as

[TABLE]

For a unique element $[x_{m},x_{m+1}]$ using transformation $(\ref{250})$ into Eq.$(\ref{8}),$ we obtain the following integral equation:

[TABLE]

where $\hat{u}$ is accepted to be constant over an element to ease the integral. Integrating Eq.$(\ref{9})$ by parts and using Eq.$(\ref{grlw})$ which yields:

[TABLE]

where $\lambda=\frac{1+p(p+1)\hat{u}^{p}}{h}$ and $\beta=\frac{\mu}{h^{2}}.$ Assuming the weight function $\Phi_{i}$ with quadratic B-spline shape functions given by Eq.$(\ref{quad})$ and substituting approximation $(\ref{6})$ into integral Eq.$(\ref{100})$, we get the element contributions in the form:

[TABLE]

where $\delta^{e}=(\delta_{m-1},\delta_{m},\delta_{m+1},\delta_{m+2})^{T}$ are the element parameters and dot states differentiation to $t$ which can be written in matrix form as follows:

[TABLE]

The element matrices $A_{ij}^{e},B_{ij}^{e},C_{ij}^{e}$ and $D_{ij}^{e}$ are rectangular $3\times 4$ given by the following integrals;

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $i$ takes the values $1,2,3$ and the $j$ takes the values $m-1,m,m+1,m+2$ for the typical element $[x_{m},x_{m+1}].$ A lumped value for $u$ is obtained from $\left(u_{m}+u_{m+1}\right)^{2}/4$ as

[TABLE]

Combining contributions from all elements induces to the following matrix equations

[TABLE]

where $\delta=(\delta_{-1},\delta_{0},...,\delta_{N},\delta_{N+1})^{T}$ global element parameters. The $A$, $B$, $C$ and $\lambda D$ matrices are rectangular and row $m$ of each has the following form:

[TABLE]

where

[TABLE]

Replacing the time derivative $\dot{\delta}$ by the forward difference approximation $\dot{\delta}=\frac{\delta^{n+1}-\delta^{n}}{\Delta t}$ and the parameter $\delta$ by the Crank-Nicolson formulation $\delta=\frac{1}{2}(\delta^{n}+\delta^{n+1}),$ then Eq. $(\ref{13})$ reduce to the following matrix system:

[TABLE]

where $t$ is time step. Applying the boundary conditions ($\ref{bndry})$ to the system $(\ref{14})$, $(N+1)\times(N+1)$ matrix system is obtained. This last system is actively solved with a variant of the Thomas algorithm but in solution process, two or three inner iterations $\delta^{n\ast}=\delta^{n}+\frac{1}{2}(\delta^{n}-\delta^{n-1})$ are also practiced at each time step to overcome the non-linearity. Ultimately, a typical member of the matrix system $(\ref{14})$  is written in terms of the nodal parameters $\delta^{n}$ and $\delta^{n+1}$ as:

[TABLE]

where

[TABLE]

To start the evolution of the vector of parameters $\delta^{n},$ $\delta^{0}$ must be calculated by using the periodic boundary condition and initial condition $u(x,0).$ So, using the relations at the knots $u_{N}(x_{m},0)=u(x_{m},0)$, $m=0,1,2,...,N$ and $u_{N}^{{}^{\prime}}(x_{0},0)=u^{{}^{\prime}}(x_{N},0)=0$ together with a variant of the Thomas algorithm, the initial vector $\delta^{0}$ is easily got from the following matrix equation

[TABLE]

5 Stability analysis

In this section, like other authors [24, 32, 40, 41] our stability analysis is based on the Von Neumann theory in which the growth factor of a typical Fourier model defined as

[TABLE]

where $k$ is mode number and $h$ is element greatness. To implement the Fourier stability analysis, Eq. $(\ref{grlw})$ needs to be linearized by assuming that the quantity $u^{p}$ in the nonlinear term $u^{p}u_{x}$ is locally constant. Substituting the Fourier mode ($\ref{260})$ into linearized scheme of ($\ref{15})$, we get

[TABLE]

where

[TABLE]

so that $|g|$ is $1$ and our linearized scheme is neutrally stable.

6 Computer implementations and illustrations

In this part, we introduce the results of the numerical experiments of our algorithm for the solution of the GRLW Eq. $(\ref{intl})$–$(\ref{bndry})$ for a single solitary wave and an interaction of two solitary waves. We also display the development of the Maxwellian initial condition into solitary waves. In order to demonstrate how favorable our numerical algorithm foresees the position and amplitude of the solution as the simulation progresses, we provides for the  following error norms:

[TABLE]

and

[TABLE]

With the boundary condition $u\rightarrow 0$ for $x\rightarrow\pm\infty$ the exact solution of the GRLW equation is [38]

[TABLE]

where $\sqrt[p]{\frac{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. There are three conserved quantities

[TABLE]

for the GRLW equation. These are correspond to mass, momentum and energy respectively.

6.1 Dispersion of a single solitary wave

In our computational work for the first set, we prefer the parameters $p=2$, $c=1$, $h=0.2$, $\Delta t$ $=0.025$, $\mu=1$, $x_{0}$ $=40$ with interval $[0,100]$ to match up with that of previous papers [32, 34, 35, 39, 40]. These values yield the amplitude $1.0$ and the run of the algorithm is continued up to time $t=10$ over the solution region. Analytical values of the invariants are $I_{1}$ $=4.442883$, $I_{2}$ $=3.299832$ and $I_{3}=1.414214.$ Values of the three invariants as well as $\mathit{L}_{2}$ and $\ \mathit{L}_{\infty}$-error norms from our method have been computed and tabulated in Table $(\ref{400})$. Referring to Table $(\ref{400}),$ the error norms $\mathit{L}_{2}$ and $\mathit{L}_{\infty}$ remain less than $2.4154685\times 10^{-3}$ and $1.07968675\times 10^{-3},$ the invariants $I_{1},I_{2}$ and $I_{3}$ change from their initial values by less than $3.10\times\ 10^{-4},$ $4.89\times\ 10^{-4}$ and $4.79\times\ 10^{-4},$ respectively, throughout the simulation. Also, our invariants are almost stable as time increases and the agreement between numerical and analytic solutions is perfect. Hence our method is acceptedly conservative. Comparisons with our results with exact solution as well as the obtained values in [32, 34, 35, 39, 40] have been made and listed in Table $(\ref{401})$ at $t=10$. This table evidentially indicates that the error norms got by our method are marginally less than the others. The motion of solitary wave using our scheme is plotted at time $t=0,5,10$ in Fig. $(\ref{500}).$ It is obvious from the figure that the suggested method performs the motion of propagation of a solitary wave admissibly, which moved to the right with the preserved amplitude and shape. Initially, the amplitude of solitary wave is $1.00000$ and its top position is pinpionted at $x=40$. At $t=10,$ its amplitude is recorded as $0.99928$ with center $x=60$. Thereby the absolute difference in amplitudes over the time interval $[0,10]$ are observed as $7.16\times 10^{-3}$. The quantile of error at discoint times are depicted in Fig.$(\ref{501})$ . The error aberration varies from $-1\times 10^{-3}$ to $1\times 10^{-3}.$

For the second set, we select the parameters $p=3,c=6/5,h=0.1,\Delta t=0.025,\mu=1$, $x_{0}$ $=40$ with interval $[0,100]$ to coincide with that of previous papers [32, 34, 35]. These parameters produce the amplitude $1.0$ and the computations are carried out for times up to $t=10.$ The error norms $L_{2},$ $L_{\infty}$ and conservation quantities $I_{1},I_{2}$ and $I_{3}$ are computed, which are recorded in Table $(\ref{402})$. According to Table $(\ref{402})$ the error norms $\mathit{L}_{2}$ and $\mathit{L}_{\infty}$ remain less than $6.12802937\times 10^{-3}$ and $3.72213891\times 10^{-3},$ the invariants $I_{1},I_{2}$ and $I_{3}$ change from their initial values by less than $9.75\times\ 10^{-5},$ $4.32\times\ 10^{-5}$ and $4.78\times\ 10^{-4},$ respectively, during the simulation. Also, our invariants are almost constant as time increases. Therefore our method is satisfactorily conservative. In Table $(\ref{403})$ the performance of the our new method is compared with other methods [32, 34, 35] at $t=10$. It is observed that errors of the method [32, 34, 35] are considerably larger than those obtained with the present scheme. Perspective views of the traveling solitons are graphed at time $t=0,5,10$ in Fig.$(\ref{502}).$ It is clear from the figure that the single soliton moved to the right with the preserved amplitude and shape. The amplitude is $1.00000$ at $t=0$ and located at $x=40$, while it is $0.99958$ at $t=10$ and located at $x=62$. The absolute difference in amplitudes over the time interval $[0,10]$ are found as $4.2\times 10^{-4}$. The aberration of error at discrete times are modelled in Fig.$(\ref{503}).$ The error deviation varies from $-4\times 10^{-3}$ to $4\times 10^{-3}$.

Finally, we choose the parameters $p=4,c=4/3,h=0.1,\Delta t=0.01,\mu=1$, $x_{0}$ $=40$ over the region $[0,100]$ to compare with those of earlier papers [32, 34, 35]. These parameters lead to amplitude $1.0$ and the simulations are executed to time $t=10$ to invent the error norms $L_{2}$ and $L_{\infty}$ and the numerical invariants $I_{1},I_{2}$ and $I_{3}.$ For these values of the parameters, the conservation properties and the $L_{2}$-error as well as the $L_{\infty}$-error norms have been given in Table$(\ref{4040})$ for various values of the time level $t$. It can be noted from Table $(\ref{4040}),$ the error norms $\mathit{L}_{2}$ and $\mathit{L}_{\infty}$ remain less than $1.28342020\times 10^{-3}$ and $0.82165081\times 10^{-3},$ the invariants $I_{1},I_{2}$ and $I_{3}$ change from their initial values by less than $9.02\times\ 10^{-5},$ $5.08\times\ 10^{-5}$ and $7.97\times\ 10^{-4},$ respectively, throughout the simulation. Also, our invariants are almost constant as time increases. Therefore we can say our method is sensibly conservative. The comparison between the results obtained by the present method with those in the other studies [32, 34, 35] is also documented in Table $(\ref{4050})$. It is noticeably seen from the table that errors of the present method are radically less than those obtained with the earlier schemes [32, 34, 35]. For visual representation, the simulations of single soliton for values $p=4,c=4/3,h=0.1,\Delta t=0.01$ at times $t=0,5$ and $10$ are illustrated in Figure $(\ref{504})$. It is understood from this figure that the numerical scheme performs the motion of propagation of a single solitary wave, which moves to the right at nearly unchanged speed and conserves its amplitude and shape with increasing time. The amplitude is $1.00000$ at $t=0$ and located at $x=40$, while it is $0.99892$ at $t=10$ and located at $x=63.3$. The absolute difference in amplitudes at times $t=0$ and $t=10$ is $1.08\times 10^{-3}$ so that there is a little change between amplitudes. Error distributions at time $t=10$ are shown graphically in Figure $(\ref{505})$. As it is seen, the maximum errors are between $-6\times 10^{-4}$ to $1\times 10^{-3}$ and occur around the central position of the solitary wave.

6.2 Interaction of two solitary waves

As a second problem, we have focused on the behavior of the interaction of two solitary waves having different amplitudes and moving in the same direction. We provide for the GRLW equation with initial  conditions given by the linear sum of two well separated solitary waves of various amplitudes

[TABLE]

where $c_{j}$ and $x_{j}$, $j=1,2$ are arbitrary constants. For the simulation, we firstly choose $p=3,c_{1}=48/5,c_{2}=6/5,$ $h=0.1$, $\Delta t=0.01,\mu=1$ over the interval $0\leq x\leq 120.$ The amplitudes are in the ratio $2:1$. Calculations are performed to time $t=6$. The results are listed in Table $(\ref{4051})$. Referring to this table, it is noticed that the numerical values of the invariants are very closed with the methods [32, 34, 35] during the computer run. The initial function was placed with the larger wave to the left of the smaller one as seen in the Fig. $(\ref{40510})$a. Both waves move to the right with velocities dependent upon their magnitudes. According to Fig. $(\ref{40510})$, the larger wave catches up with the smaller wave at about $t=3$, the overlapping process continues until $t=4$, then two solitary waves emerge from the interaction and resume their former shapes and amplitudes. At $t=6$, the magnitude of the smaller wave is $1.00029$ on reaching position $x=60.0$, and of the larger wave $1.99213$ having the position $x=85.3$, so that the difference in amplitudes is $0.00029$ for the smaller wave and $0.00787$ for the larger wave. The changes of the invariants for this case are satisfactorily small. Secondly, to ensure an interaction of two solitary waves take place, calculation is carried out with the parameters  $p=4,c_{1}=64/3,c_{2}=4/3,$ $h=0.125$, $\Delta t=0.01,\mu=1$ over the interval $0\leq x\leq 200.$ The parameters give solitary waves of different amplitudes $2$ and $1$ having centers at $x_{1}=20$ and $x_{2}=80.$ The results are given in Table $(\ref{40530}).$ According to the this table, it is realized that the numerical values of the invariants are very closed with the methods [32, 34, 35] during the computer run. The initial function was placed with the larger wave to the left of the smaller one as seen in the Fig. $(\ref{4052})$a. Both waves move to the right with velocities dependent upon their magnitudes. According to Fig. $(\ref{4052})$, the larger wave catches up with the smaller wave at about $t=3$, the overlapping process continues until $t=5$, then two solitary waves emerge from the interaction and resume their former shapes and amplitudes.

6.3 The Maxwellian initial condition

Finally, we have examined the evolution of an initial Maxwellian pulse into solitary waves, arising as initial condition

[TABLE]

For this problem, the behavior of the solution depends on the value of $\mu$ [13, 32]. Therefore, we chose the values of $\mu=0.1,$ $\mu=0.05$ and $\mu=0.025$ for $p=2,3,4$. The numerical computations are done up to $t=0.05$. Calculated numerical invariants at different values of $t$ are shown in Table $(\ref{4053})$ and it is seen that calculated invariant values are satisfactorily constant. For $p=2$ and $\mu=0.1$; the variation of invariants $I_{1},I_{2}$ and $I_{3}$ from initial variants changes less than $1.02\times 10^{-3},4.48\times 10^{-3}$ and $8.69\times 10^{-3}$percent, respectively, and for $\mu=0.05;$ $2.01\times 10^{-3},8.35\times 10^{-3}$ and $17.92\times 10^{-3}$percent, respectively, and for $\mu=0.025;$ $3.19\times 10^{-3},12.72\times 10^{-3}$ and $29.28\times 10^{-3}$percent, respectively, for $p=3$ and $\mu=0.1$; the variation of invariants $I_{1},I_{2}$ and $I_{3}$ from initial variants changes less than $8.60\times 10^{-3},2.77\times 10^{-2}$ and $4.75\times 10^{-1}$percent, respectively, and for $\mu=0.05;$ $17.45\times 10^{-3},54.32\times 10^{-3}$ and $6.62\times 10^{-1}$percent, respectively, and for $\mu=0.025;$ $30.92\times 10^{-3},93.27\times 10^{-3}$ and $7.72\times 10^{-1}$percent, respectively, for $p=4$ and $\mu=0.1$; the variation of invariants $I_{1},I_{2}$ and $I_{3}$ from initial variants changes less than $38.82\times 10^{-3},1.09\times 10^{-3}$ and $1.76$ percent, respectively, and for $\mu=0.05;$ $67.23\times 10^{-3},1.99\times 10^{-1}$ and $2.48$ percent, respectively, and for $\mu=0.025;$ $1.04\times 10^{-1},3.31\times 10^{-1}$ and $3.0$ percent, respectively. The development of the Maxwellian initial condition into solitary waves is shown in Fig. $(\ref{4055})$.

7 Conclusion

In this work, a numerical technique based on a Petrov-Galerkin method using quadratic weight functions and cubic B-spline finite elements has been proffered to get numerical solutions of GRLW equation. We experimented our algorithm along with single solitary wave in which the exact solution is known and broadened it to examine the interaction of two solitary waves and Maxwellian initial condition where the exact solutions are unknown during the interaction. Variational formulation and semi-discrete Galerkin scheme of the equation are generated. Stability analysis have been done and the linearized numerical scheme have been obtained unconditionally stable. The accuracy of the method is investigated both $L_{2}$ and $L_{\infty}$ error norms and the invariant quantities $I_{1}$, $I_{2}$ and $I_{3}$. The obtained numerical results indicate that the error norms are satisfactorily small and the conservation laws are marginally constant in all computer program run. We can see that our numerical scheme for the equation is more accurate than the other earlier schemes found in the literature. Therefore, our numerical technique is suitable for getting numerical solutions of partial differential equations.