Numerical solutions of the generalized equal width wave equation using Petrov Galerkin method

TL;DR

This paper develops a Petrov-Galerkin numerical scheme with quadratic shape functions and linear B-spline weights to solve the generalized equal width wave equation, analyzing its stability, error bounds, and effectiveness in modeling soliton propagation.

Contribution

The paper introduces a novel Petrov-Galerkin method for the GEW equation, providing theoretical error analysis and demonstrating superior performance over existing schemes.

Findings

The scheme is unconditionally stable.

It accurately models solitary wave propagation.

It outperforms other numerical methods in efficiency and accuracy.

Abstract

In this article we consider a generalized equal width wave (GEW) equation which is a significant nonlinear wave equation as it can be used to model many problems occurring in applied sciences. As the analytic solution of the (GEW) equation of this kind can be obtained hardly, developing numerical solutions for this type of equations is of enormous importance and interest. Here we are interested in a Petrov-Galerkin method, in which element shape functions are quadratic and weight functions are linear B-splines. We firstly investigate the existence and uniqueness of solutions of the weak form of the equation. Then we establish the theoretical bound of the error in the semi-discrete spatial scheme as well as of a full discrete scheme at $t=t^{n}$. Furthermore, a powerful Fourier analysis has been applied to show that the proposed scheme is unconditionally stable. Finally, propagation of single and double solitary waves and evolution of solitons are analyzed to demonstrate the efficiency and applicability of the proposed numerical scheme by calculating the error norms (in $L_{2}(\Omega)$ and $L_{\infty}(\Omega)$). The three invariants ($% I_{1},I_{2}$ and $I_{3})$ of motion have been commented to verify the conservation features of the proposed algorithms. Our proposed numerical scheme has been compared with other published schemes and demonstrated to be valid, effective and it outperforms the others.

Tables

Table 1. Table 1: Invariants and errors for single solitary wave with p = 2 , 𝑝 2 p=2, c = 0.5 , 𝑐 0.5 c=0.5, h = 0.1 , ℎ 0.1 h=0.1, ε = 3 , 𝜀 3 \varepsilon=3, Δ ​ t = 0.2 , Δ 𝑡 0.2 \Delta t=0.2, μ = 1 , 𝜇 1 \mu=1, x ∈ [ 0 , 80 ] . 𝑥 0 80 x\in\left[0,80\right].

$Time$	$I_1$	$I_2$	$I_3$	$L_2$	$L_{\infty }$
0	3.1415863	2.6682242	1.3333283	0.00000000	0.00000000
5	3.1415916	2.6682311	1.3333406	0.00395289	0.00294851
10	3.1415934	2.6682352	1.3333413	0.00704492	0.00473785
15	3.1415948	2.6682434	1.3333413	0.00995547	0.00651735
20	3.1415961	2.6682568	1.3333413	0.01286582	0.00831346

Table 2. Table 2: Comparisons of results for single solitary wave with p = 2 , 𝑝 2 p=2, c = 0.5 , 𝑐 0.5 c=0.5, h = 0.1 , ℎ 0.1 h=0.1, ε = 3 , 𝜀 3 \varepsilon=3, Δ ​ t = 0.2 , Δ 𝑡 0.2 \Delta t=0.2, μ = 1 , 𝜇 1 \mu=1, x ∈ [ 0 , 80 ] 𝑥 0 80 x\in\left[0,80\right] at t = 20 . 𝑡 20 t=20.

$Method$	$I_1$	$I_2$	$I_3$	$L_2$	$L_{\infty }$
Analytic	3.1415961	2.6666667	1.3333333	0.00000000	0.00000000
Our Method	3.1415916	2.6682568	1.3333413	0.01286582	0.00831346
Cubic Galerkin[8]	3.1589605	2.6902580	1.3570299	0.03803037	0.02629007
Quintic Collocation First Scheme[31]	3.1250343	2.6445829	1.3113394	0.05132106	0.03416753
Quintic Collocation Second Scheme[31]	3.1416722	2.6669051	1.3335718	0.01675092	0.01026391
Petrov-Galerkin[32]	3.14159	2.66673	1.33341	0.0123326	0.0086082

Table 3. Table 3: Invariants and errors for single solitary wave with p = 3 , 𝑝 3 p=3, c = 0.3 , 𝑐 0.3 c=0.3, h = 0.1 , ℎ 0.1 h=0.1, Δ ​ t = 0.2 , Δ 𝑡 0.2 \Delta t=0.2, ε = 3 , 𝜀 3 \varepsilon=3, μ = 1 , 𝜇 1 \mu=1, x ∈ [ 0 , 80 ] . 𝑥 0 80 x\in\left[0,80\right].

$Time$	$I_1$	$I_2$	$I_3$	$L_2$	$L_{\infty }$
0	2.8043580	2.4664883	0.9855618	0.00000000	0.00000000
5	2.8043723	2.4665080	0.9855942	0.00183258	0.00177948
10	2.8043747	2.4665108	0.9855973	0.00291958	0.00233283
15	2.8043753	2.4665119	0.9855973	0.00372417	0.00285444
20	2.8043758	2.4665135	0.9855973	0.00448357	0.00337609

Table 4. Table 4: Comparisons of results for single solitary wave with p = 3 , 𝑝 3 p=3, c = 0.3 , 𝑐 0.3 c=0.3, h = 0.1 , ℎ 0.1 h=0.1, Δ ​ t = 0.2 , Δ 𝑡 0.2 \Delta t=0.2, ε = 3 , 𝜀 3 \varepsilon=3, μ = 1 , 𝜇 1 \mu=1, x ∈ [ 0 , 80 ] 𝑥 0 80 x\in\left[0,80\right] at t = 20 . 𝑡 20 t=20.

$Method$	$I_1$	$I_2$	$I_3$	$L_2$	$L_{\infty }$
Our Method	2.8043758	2.4665135	0..9855973	0.00448357	0.00337609
Cubic Galerkin[8]	2.8187398	2.4852249	1.0070200	0.01655637	0.01370453
Quintic Collocation First Scheme[31]	2.8043570	2.4639086	0.9855602	0.00801470	0.00538237
Quintic Collocation Second Scheme[31]	2.8042943	2.4637495	0.9854011	0.00708553	0.00480470
Petrov-Galerkin[32]	2.80436	2.46389	0.98556	0.00484271	0.00370926

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

$Time$	$I_1$	$I_2$	$I_3$	$L_2$	$L_{\infty }$
0	2.6220516	2.3598323	0.7853952	0.00000000	0.00000000
5	2.6220846	2.3598808	0.7854675	0.00125061	0.00141788
10	2.6220915	2.3598891	0.7854783	0.00178634	0.00147002
15	2.6220920	2.3598898	0.7854785	0.00193428	0.00139936
20	2.6220923	2.3598903	0.7854785	0.00196046	0.00133416

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

$Method$	$I_1$	$I_2$	$I_3$	$L_2$	$L_{\infty }$
Our Method	2.6220923	2.3598903	0.7854785	0.00196046	0.00133416
Cubic Galerkin[8]	2.6327833	2.3730032	0.8023383	0.00890617	0.00821991
Quintic Collocation First Scheme[31]	2.6220508	2.3561901	0.7853939	0.00421697	0.00297952
Quintic Collocation First Scheme[31]	2.6219284	2.3559327	0.7851364	0.00339086	0.00247031
Petrov-Galerkin[32]	2.62206	2.35615	0.78534	0.00230499	0.00188285

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

	$t$	$0$	$30$	$60$	$90$	$100$
	Our Method	4.20653	4.20657	4.20622	4.20502	4.20517
$I_1$	[8]	4.20653	4.20653	4.20616	4.20490	4.20503
	[31] First	4.20653	4.20653	4.20653	4.20653	4.20653
	[31] Second	4.20653	4.20653	4.20653	4.20653	4.20653
	[32]	4.20655	4.20655	4.20655	4.20655	4.20655
	Our Method	3.08311	3.08318	3.08309	3.08220	3.08251
$I_2$	[8]	3.07987	3.07991	3.07947	3.07777	3.07797
	[31] First	3.07988	3.07988	3.07988	3.07988	3.07988
	[31] Second	3.07988	3.07988	3.07988	3.07988	3.07988
	[32]	3.97977	3.07980	3.07987	3.07974	3.07972
	Our Method	1.01636	1.01644	1.01664	1.01632	1.01634
$I_3$	[8]	1.01636	1.01638	1.01654	1.01616	1.01616
	[31] First	1.01636	1.01636	1.01636	1.01636	1.01636
	[31] Second	1.01636	1.01636	1.01636	1.01636	1.01636
	[32]	1.01634	1.01634	1.01634	1.01633	1.01634

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

	$t$	$0$	$30$	$60$	$90$	120
	Our Method	3.93307	3.93311	3.93393	3.93229	3.93037
$I_1$	[8]	3.93307	3.93309	3.93388	3.93222	3.93026
	[31] First	3.93307	3.93307	3.93307	3.93307	3.93307
	[31] Second	3.93307	3.93307	3.93307	3.93307	3.93307
	[32]	3.93309	3.93309	3.93309	3.93309	3.93308
	Our Method	2.94979	2.94985	2.95122	2.94939	2.94801
$I_2$	[8]	2.94521	2.94527	2.94703	2.94436	2.94212
	[31] First	2.94524	2.94524	2.94524	2.94524	2.94524
	[31] Second	2.94524	2.94523	2.94523	2.94523	2.94523
	[32]	2.94512	2.94510	2.94505	2.94520	2.94511
	Our Method	0.79766	0.79775	0.79952	0.79824	0.79811
$I_3$	[8]	0.79766	0.79770	0.79942	0.79812	0.79794
	[31] First	0.79766	0.79766	0.79766	0.79766	0.79766
	[31] Second	0.79766	0.79766	0.79766	0.79766	0.79766
	[32]	0.79761	0.79761	0.79762	0.79761	0.79761

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	1.7724537	1.3792767	0.8862269	1.7724537	1.3792767	0.7926655	1.7724537	1.3792767	0.7236013
0.1	4	1.7724537	1.5760586	0.8862269	1.7724537	1.6168691	0.7926655	1.7724537	1.6360543	0.7236013
	8	1.7724537	1.5838481	0.8862269	1.7724537	1.6245008	0.7926655	1.7724537	1.6481131	0.7236013
	12	1.7724537	1.5920722	0.8862269	1.7724537	1.6325922	0.7926655	1.7724537	1.6531844	0.7236013
[31]	12	1.7724	1.3786	0.8862	1.7724	1.3786	0.7928	1.7725	1.3786	0.7243
[32]	12	1.7724	1.3785	0.8861	1.7724	1.3787	0.7926	1.7734	1.3836	0.7224
	0	1.7724537	1.3162954	0.8862269	1.7724537	1.3162954	0.7926655	1.7724537	1.3162954	0.7236013
	4	1.7724537	1.5406812	0.8862269	1.7724537	1.5766908	0.7926655	1.7724537	1.6243519	0.7236013
0.05	8	1.7724537	1.6342604	0.8862269	1.7724537	1.6367952	0.7926655	1.7724537	1.6554614	0.7236013
	12	1.7724537	1.6835979	0.8862269	1.7724537	1.6372439	0.7926655	1.7724537	1.7079133	0.7236013
[31]	12	1.7724	1.3159	0.8864	1.7725	1.3160	0.7940	1.7735	1.3188	0.7345
[32]	12	1.7724	1.3160	0.8861	1.7724	1.3156	0.7922	1.7724	1.3177	0.7245

Equations

$$U_t+U_x+\epsilon U{U}_x-\mu U_{xxt}=0,$$

$$U_t+\epsilon U{U}_x+\mu U_{xxx}=0,$$

$$U_t+\epsilon U{U}_x-\mu U_{xxt}=0,$$

$$U_t+\epsilon U^p{U}_x-\mu U_{xxt}=0,$$

$$\begin{array}[]{l}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,\\ U(x,0)=f(x),~{}~{}\ \ \ \ a\leq x\leq b,\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(\Omega ),$$

$$\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$$

$$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 ).$$

$$\underset{\xi \in S_h}{\inf }\Vert u-\xi \Vert \le C{h}^3\Vert u{\Vert }_3.$$

$$\left(u_{ht},\xi \right)+\left(\nabla u_{ht},\nabla \xi \right)=-\left(\mathcal{F}(u_h),\nabla \xi \right),\forall \xi \in S_h,$$

$$\|u_{h}\|_{1}^{2}=\|u_{0,h}\|_{1}^{2},\ t\in\ (0,\ T],\$$

$$\Vert u_h{\Vert }_{L^{\infty }(L^{\infty }(\Omega ))}\le C\Vert u_{0,h}{\Vert }_1$$

$$\mathcal{A}(u,v)=(\nabla u,\nabla v),\forall u,v\in H_0^1(\Omega ),$$

$$|\mathcal{A}(u,v)|\le M\Vert u{\Vert }_1\Vert v{\Vert }_1,\forall u,v\in H_0^1(\Omega )$$

$$\mathcal{A}(u,u)\ge \alpha \Vert u{\Vert }_1,\forall u\in H_0^1(\Omega ),\text{for some }\alpha \in \mathbb{R}.$$

$$\mathcal{A}(u-\tilde{u},\xi )=0,\xi \in S_h.$$

$$\Vert u-u_h\Vert \le C{h}^3,$$

$$\alpha \Vert u-\tilde{u}{\Vert }_1^2$$

$$=\mathcal{A}(u-\tilde{u},u-\xi ),\xi \in S_h.$$

$$\Vert u-\tilde{u}{\Vert }_1\le \underset{\xi \in S_h}{\inf }\Vert u-\xi {\Vert }_1,$$

$$\Vert \psi {\Vert }_1\le C{h}^2\Vert u{\Vert }_3,\text{and so }\Vert \psi \Vert \le C{h}^3\Vert u{\Vert }_3.$$

$$\Vert {\psi }_t\Vert \le C{h}^3\Vert u_t{\Vert }_3.$$

$$({\theta }_t,\xi )+(\nabla {\theta }_t,\nabla \xi )=({\psi }_t,\xi )-(\mathcal{F}(u)-\mathcal{F}(u_h),\nabla \xi ).$$

$$\frac{1}{2}\frac{d}{dt}\Vert \theta {\Vert }_1^2\le \Vert {\psi }_t\Vert \Vert \theta \Vert +\Vert \mathcal{F}(u)-\mathcal{F}(u_h)\Vert \Vert \nabla \theta \Vert .$$

$$\Vert \mathcal{F}(u)-\mathcal{F}(u_h)\Vert \le C(\Vert \psi \Vert +\Vert \theta \Vert ),$$

$$\frac{d}{dt}\Vert \theta {\Vert }_1^2\le C\left(\Vert {\psi }_t{\Vert }^2+\Vert \psi {\Vert }^2+\Vert \theta {\Vert }^2+\Vert \nabla \theta {\Vert }^2\right).$$

Full text

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

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

1 Introduction

Nonlinear partial differential equations are extensively used to explain complex phenomena in different fields of science, such as plasma physics, fluid mechanics, hydrodynamics, applied mathematics, solid state physics and optical fibers. One of the important issues to nonlinear partial differential equations is to seek for exact solutions. Because of the complexity of nonlinear differential equations, exact solutions of these equations are commonly not derivable. Owing to the fact that only limited classes of these equations are solved by analytical means, numerical solutions of these nonlinear partial differential equations are very functional to examine physical phenomena. The regularized long wave (RLW) equation,

[TABLE]

is a symbolisation figure of nonlinear long wave and can define many important physical phenomena with weak nonlinearity and dispersion waves, including nonlinear transverse waves in shallow water, ion-acoustic and magneto hydrodynamic waves in plasma, elastic media, optical fibres, acoustic-gravity waves in compressible fluids, pressure waves in liquid–gas bubbles and phonon packets in nonlinear crystals [1]. The RLW equation was first suggested to describe the behavior of the undular bore by Peregrine [2, 3], who constructed the first numerical method of the equation using finite difference method. RLW equation is an alternative description of nonlinear dispersive waves to the more usual

[TABLE]

Korteweg-de Vries (KdV) equation [4]. This equation was first generated by Korteweg and de Vries to symbolise the action of one dimensional shallow water solitary waves [5]. The equation has found numerous applications in physical sciences and engineering field such as fluid and quantum mechanics, plasma physics, nonlinear optics, waves in enharmonic crystals, bubble liquid mixtures, ion acoustic wave and magneto-hydrodynamic waves in a warm plasma as well as shallow water waves. The Equal Width (EW) wave equation

[TABLE]

which is less well recognised and was introduced by Morrison et al. [6] is a description alternative to the more common KdV and RLW equations. This equation is named equal width equation, because the solutions for solitary waves with a perpetual form and speed, for a given value of the parameter $\mu$, are waves with an equal width or vawelength for all wave amplitudes [7]. The solutions of this equation are sorts of solitary waves called as solitons whose figures are not changed after the collision. GEW equation, procured for long waves propagating in the positive $x$ direction takes the form

[TABLE]

where $p$ is a positive integer, $\varepsilon$ and $\mu$ are positive parameters, $t$ is time and $x$ is the space coordinate, $U(x,t)$ is the wave amplitude. Physical boundary conditions require $U\rightarrow 0$ as $\left|x\right|\rightarrow\infty$. For this work, boundary and initial conditions are chosen

[TABLE]

where $f(x)$ is a localized disturbance inside the considered interval and will be designated later. In the fluid problems as known, the quantity $U$  is associated with the vertical displacement of the water surface but in the plasma applications, $U$ is the negative of the electrostatic potential. That’s why, the solitary wave solution of Eq.$(\ref{gew})$ helps us to find out the a lot of physical phenomena with weak nonlinearity and dispersion waves such as nonlinear transverse waves in shallow water, ion-acoustic and magneto- hydrodynamic waves in plasma and phonon packets in nonlinear crystals [8]. The GEW equation which we tackle here is based on the EW equation and relevant to the both generalized regularized long wave (GRLW) equation [9, 10] and the generalized Korteweg-de Vries (GKdV) equation [11]. These general equations are nonlinear wave equations with $(p+1)$th nonlinearity and have solitary wave solutions, which are pulse-like. The investigate of GEW equation ensures the possibility of investigating the creation of secondary solitary waves and/or radiation to get insight into the corresponding processes of particle physics [12, 13]. This equation has many implementations in physical situations for example unidirectional waves propagating in a water channel, long waves in near-shore zones, and many others [14]. If $p=1$ is taken in Eq.$(\ref{gew})$ the EW equation [15-20] is obtained and if $p=2$ is taken in Eq.$(\ref{gew}),$ the obtained equation is named as the modified equal width wave (MEW) equation [21-27]. In recent years, various numerical methods have been improved for the solution of the GEW equation. Hamdi et al. [7] generated exact solitary wave solutions of the GEW equation. Evans and Raslan [28] investigated the GEW equation by using the collocation method based on quadratic B-splines to obtain the numerical solutions of the single solitary wave, interaction of solitary waves and birth of solitons. The GEW equation solved numerically by a B-spline collocation method by Raslan [29]. The homogeneous balance method was used to construct exact travelling wave solutions of generalized equal width equation by Taghizadeh et al. [30]. The equation is solved numerically by a meshless method based on a global collocation with standard types of radial basis functions (RBFs) by [14]. Quintic B-spline collocation method with two different linearization techniques and  a lumped Galerkin method based on B-spline functions were employed to obtain the numerical solutions of the GEW equation by Karakoc and Zeybek, [8, 31] respectively. Roshan [32], applied Petrov-Galerkin method using the linear hat function and quadratic B-spline functions as test and trial function respectively for the GEW equation.

In this study, we have constructed a lumped Petrov-Galerkin method for the GEW equation using quadratic B-spline function as element shape function and linear B-spline function as the weight function. Context of this work has been planned as follows:

A semi-discrete Galerkin finite element scheme of the equation along with the error bounds are demonstrated in Section 2.

• -

A full discrete Galerkin finite element scheme has been studied in Section 3.

• -

Section 4 is concerned with the construction and implementation of the Petrov-Galerkin finite element method to the GEW equation.

• -

Section 5 contains a linear stability analysis of the scheme.

• -

Section 6 includes analysis of the motion of single solitary wave, interaction of two solitary wave and evolution of solitons with different initial and boundary conditions.

• -

Finally, we conclude the study with some remarks on this study.

2 Variational formulation and its analysis

The higher order nonlinear initial boundary value problem (4) can be written as

[TABLE]

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

[TABLE]

and the boundary conditions

[TABLE]

To define the weak form of the solutions of (6) 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 considered as an usual normed space of real valued functions on $\Omega$ and

[TABLE]

where $D=\frac{\partial}{\partial x}$. We denote the norm on $H^{k}(\Omega)$ by $\|\cdot\|_{k}$ which is the well known usual $H^{k}$ norm, and when $k=0$, $\|\cdot\|_{0}=\|\cdot\|$ represents $L_{2}$ norm and $(\cdot,\cdot)$ represents the standard $L_{2}$ inner product [33, 34].

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

[TABLE]

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

[TABLE]

with $u(0)=u_{0}$. Here we state the uniqueness theorem without proof which can be well established following [33, 34].

Theorem 1**.**

If $u$ satisfies (9) then

[TABLE]

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

Theorem 2**.**

Assume that $u_{0}\in H_{0}^{1}(\Omega)$ and $T>0$. Then there exists one and only one $u\in H_{0}^{1}(\Omega)$ satisfying (9) for any $T>0$ such that

[TABLE]

2.1 Semi-discrete Galerkin Scheme

For any $0<h<1$ let $S_{h}$ of $H_{0}^{1}(\Omega)$ be a finite dimensional subspace such that for $u\in H_{0}^{1}(\Omega)\cap H^{3}(\Omega)$, $\exists$ a constant $C$ independent of $h$ [33, 34, 35] such that

[TABLE]

Here it is our moto to look for solutions of a semi-discrete finite element formulation of (6) $u_{h}:[0,\ T]\rightarrow S_{h}$ such that

[TABLE]

where $u_{h}(0)=u_{0,h}\in S_{h}$ approximates $u_{0}$. We start here first by stating a priori bound of the solution of (11) below before establishing the original convergence result.

Theorem 3**.**

Let $u_{h}\in S_{h}$ be a solution of (11). Then $u_{h}\in S_{h}$ satisfies

[TABLE]

*and ***

[TABLE]

holds where $C$ is a positive constant.

Proof.

The proof is trivial, it follows from [36]. ∎

Our next goal is to establish the theoretical estimate of the error in the semi-discrete scheme (11) of (9). To that end here we start by considering the following bilinear form

[TABLE]

which satisfies the boundedness property

[TABLE]

and coercivity property (on $\Omega$)

[TABLE]

Let $\tilde{u}$ be an auxiliary projection of $u$ [33, 35, 34], then $\mathcal{A}$ satisfies

[TABLE]

Now the rate of convergence (accuracy) in such a spatial approximation (11) of (9) is given by the following theorem.

Theorem 4**.**

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

[TABLE]

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

Proof.

Letting $\mathcal{E}=u-u_{h}=\psi+\theta,$ where $\psi=u-\tilde{u}$ and $\theta=\tilde{u}-u_{h}$, we write

[TABLE]

From (12), (14) and [34] it follows that

[TABLE]

and thus (10) and (15) confirms the following inequalities

[TABLE]

Now applying $\frac{\partial}{\partial t}$ on (14) and having some simplifications yields [34]

[TABLE]

Also we subtract (11) from (9) to obtain

[TABLE]

Now we substitute $\xi=\theta$ in (16), and then apply Cauchy-Schwarz inequality to obtain

[TABLE]

Here

[TABLE]

comes from Lipschitz conditions of $\mathcal{F}$ and boundedness of $u$ and $u_{h}$. Thus

[TABLE]

So

[TABLE]

Hence Gronwall’s lemma, bounds of $\psi$ and $\psi_{t}$ confirms

[TABLE]

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

3 Full discrete scheme

Here we aim to find solution of the semi-discrete problem (11) over $[0,T]$, $T>0$. Let $N$ be a positive full number and $\Delta t=\frac{T}{N}$ so that $t^{n}=n\Delta t$, $n=0,\ 1,\ 2,\ 3,\cdots,\ N.$ Here we consider

[TABLE]

Using the above notations we present a time discretized finite element Galerkin scheme by

[TABLE]

where $U^{0}=u_{0,h}$.

Theorem 5**.**

If $U^{n}$ satisfies (17) then

[TABLE]

and there exists a positive constant $C$ such that

[TABLE]

Proof.

Substituting $\xi=U^{n-1/2}$ in (17) it is easy to see that

[TABLE]

Thus the proof of the first part of the theorem follows from a sum from $n=1$ to $J$ and that of the second part follows from the Sovolev embedding theorem [34]. ∎

Now we focus on to establishing the theoretical upper bound of the error in such a full discrete approximation (18) at $t=t^{n}$.

Theorem 6**.**

Let $h$ and $\Delta t$ be sufficiently small, then

[TABLE]

where $C$ is independent of $h$ and $\Delta t$.

Proof.

Let

[TABLE]

where $\psi^{n}=u^{n}-\tilde{u^{n}}$, $\theta^{n}=\tilde{u^{n}}-U^{n}$, $u^{n}=u(t^{n})$, and $\tilde{u^{n}}=\tilde{u}(t^{n})$. From (9) and (17) along with auxiliary projection defined in the previous section the following equality holds

[TABLE]

where $\tau^{n}=u^{n-1/2}-\partial_{t}u^{n}$. Now substituting $\xi$ by $\theta^{n-1/2}$ in (19) yields

[TABLE]

Now

[TABLE]

and from boundedness of $\|U^{n}\|_{\infty}$ and $\|u^{n}\|_{\infty}$ it yields

[TABLE]

since $\mathcal{F}$ is a Lipschitz function. Thus from (20), (21) and (22) it follows that

[TABLE]

So (23) can be simplified as

[TABLE]

Choosing $\Delta t>0$ so that $1-C\Delta t\geq 0$ and summing over $n=1,(1),J$, and from the bounds of $\|\psi^{n}\|$ and $\|\partial_{t}\psi^{n}\|$ yields

[TABLE]

and the rest follows from the triangular inequality and sobolev embedding theorem [34, 35]. ∎

4 Construction and Implementation of the method

We take into account a uniformly spatially distributed set of knots $a=x_{0}<x_{1}<...<x_{N}=b$ over the solution interval $a\leq x\leq b$ and $h=x_{m+1}-x_{m},$ $m=0,1,2,...,N$. For this partition, we shall need the following quadratic B-splines $\phi_{m}(x)$ at the points $x_{m},$ $m=0,1,2,...,N.$ Prenter [37] identified following quadratic B-spline functions $\phi_{m}(x)$, (m= $-1(1)$ $N$), at the points $x_{m}$ which generate a basis over the interval $[a,b]$ by

[TABLE]

We search the approximation $U_{N}(x,t)$ to the solution $U(x,t),$ which use these splines as the trial functions

[TABLE]

in which unknown parameters $\delta_{j}(t)$ will be computed by using the boundary and weighted residual conditions. In each element, using $h\eta=x-x_{m}$ $(0\leq\eta\leq 1)$ local coordinate transformation for the finite element $[x_{m},x_{m+1}],$ quadratic B-spline shape functions $(\ref{3})$ in terms of $\eta$ over the interval $[0,1]$ can be reformulated as

[TABLE]

All quadratic B-splines, except that $\phi_{m-1}(x),\phi_{m}(x)$ and $\phi_{m+1}(x)$ are zero over the interval $[x_{m},x_{m+1}].$ Therefore approximation function $(\ref{4})$ over this element can be given in terms of the basis functions $(\ref{5})$ as

[TABLE]

Using quadratic B-splines $(\ref{5})$ and the approximation function $(\ref{6}),$ the nodal values $U_{m}$ and $U_{m}^{\prime}$ at the knot are found in terms of element parameters $\delta_{m}$ as follows:

[TABLE]

Here weight functions $L_{m}$ are used as linear B-splines. The linear B-splines $L_{m}$ at the knots $x_{m}$ are identified as [37]:

[TABLE]

A characteristic finite interval $[x_{m},x_{m+1}]$ is turned into the interval $[0,1]$ by local coordinates $\eta$ concerned with the global coordinates using $h\eta=x-x_{m}$ $(0\leq\eta\leq 1).$ So linear B-splines $L_{m}$ are given as

[TABLE]

Using the Petrov-Galerkin method to Eq.$(\ref{gew}),$ we obtain the weak form of Eq.$(\ref{gew})$ as

[TABLE]

Applying the change of variable $x\rightarrow\eta$ into Eq.$(\ref{8})$ gives rise to

[TABLE]

where $\hat{U}$ is got to be constant over an element to make the integral easier. Integrating Eq.$(\ref{9})$ by parts and using Eq.$(\ref{gew})$ leads to

[TABLE]

where $\lambda=\frac{\varepsilon\hat{U}^{p}}{h}$ and $\beta=\frac{\mu}{h^{2}}.$ Choosing the weight functions $L_{m}$ with linear B-spline shape functions given by $(\ref{lin1})$ and replacing approximation $(\ref{6})$ into Eq.$(\ref{100})$ over the element $[0,1]$ produces

[TABLE]

which can be obtained in matrix form as

[TABLE]

In the above equations and overall the article, the dot denotes differentiation according to $t$ and $\delta^{e}=(\delta_{m-1},\delta_{m},\delta_{m+1},\delta_{m+2})^{T}$ are the element parameters. $A_{ij}^{e},B_{ij}^{e},C_{ij}^{e}$ and $D_{ij}^{e}$ are the $2\times 3$ rectangular element matrices represented by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $i$ takes $m,m+1$ and $j$ takes $m-1,m,m+1$ for the typical element $[x_{m},x_{m+1}].$ A lumped value for $U$ is attained from ($\frac{U_{m}+U_{m+1}}{2}$)p as

[TABLE]

Formally aggregating together contributions from all elements leads to the matrix equation

[TABLE]

where global element parameters are $\delta=(\delta_{-1},\delta_{0},...,\delta_{N},\delta_{N+1})^{T}$ and the $A$, $B,C$ and $\lambda D$ matrices are derived from the corresponding element matrices $A_{ij}^{e},B_{ij}^{e},C_{ij}^{e}$ and $D_{ij}^{e}.$ Row $m$ of each matrices has the following form;

[TABLE]

where

[TABLE]

Implementing the Crank-Nicholson approach $\delta=\frac{1}{2}(\delta^{n}+\delta^{n+1})$ and the forward finite difference $\dot{\delta}=\frac{\delta^{n+1}-\delta^{n}}{\Delta t}$ in Eq.$(\ref{120})$ we get the following matrix system:

[TABLE]

where $\Delta t$ is time step. Implementing the boundary conditions ($\ref{Bou.Con.})$ to the system $(\ref{14})$, we make the matrix equation square. This system is efficaciously 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 performed at each time step to cope with the nonlinearity. As a result, a typical member of the matrix system $(\ref{14})$  may be written in terms of the nodal parameters $\delta^{n}$ and $\delta^{n+1}$ as:

[TABLE]

where

[TABLE]

To start the iteration for computing the unknown parameters, the initial unknown vector $\delta^{0}$ is calculated by using Eqs.($\ref{Bou.Con.}).$ Therefore, 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$ related with a variant of the Thomas algorithm, the initial vector $\delta^{0}$ is easily obtained from the following matrix form

[TABLE]

5 Stability analysis

In this section, to show the stability analysis of the numerical method, we have used Fourier method based on Von-Neumann theory and presume that the quantity $U^{p}$ in the nonlinear term $U^{p}U_{x}$ of the equation $(\ref{gew})$ is locally constant. Substituting the Fourier mode $\delta_{j}^{n}=g^{n}e^{ijkh}$ where $k$ is mode number and $h$ is element size, into scheme ($\ref{15})$

[TABLE]

is obtained and where

[TABLE]

$|g|$ is found $1$ so our linearized scheme is unconditionally stable.

6 Computational results and discussions

The objective of this section is to investigate the deduced algorithm using different test problems relevant to the dispersion of single solitary waves, interaction of two solitary waves and the evolution of solitons. For the test problems, we have calculated the numerical solution of the GEW equation for $p=2,3$ and $4$ using the homogenous boundary conditions and different initial conditions. The $\mathit{L}_{2}$

[TABLE]

and $\mathit{L}_{\infty}$

[TABLE]

error norms are considered to measure the efficiency and accuracy of the present algorithm and to compare our results with both exact values, Eq.$(\ref{exct})$, as well as other results in the literature whenever available. The exact solution of the GEW equation is taken [28, 31] to be

[TABLE]

which corresponds to a solitary wave of amplitude $\sqrt[p]{\frac{c(p+1)(p+2)}{2\varepsilon},}$ the speed of the wave traveling in the positive direction of the $x$-axis is $c$, width $\frac{p}{2\sqrt{\mu}}$ and $x_{0}$ is arbitrary constant. With the homogenous boundary conditions, solutions of GEW equation possess three invariants of the motion introduced by

[TABLE]

related to mass, momentum and energy, respectively.

6.1 Propagation of single solitary waves

For the numerical study in this case, we firstly select $p=2$, $c=0.5$, $h=0.1$, $\Delta t$ $=0.2$, $\mu=1$, $\varepsilon=3$ and $x_{0}$ $=30$ through the interval $[0,80]$ to match up with that of previous papers [8, 31, 32]. These parameters represent the motion of a single solitary wave with amplitude $1.0$ and the program is performed to time $t=20$ over the solution interval. The analytical values of conservation quantities are $I_{1}$ $=3.1415927$, $I_{2}$ $=2.6666667$ and $I_{3}=1.3333333.$ Values of the three invariants as well as $\mathit{L}_{2}$ and $\ \mathit{L}_{\infty}$-error norms from our method have been found and noted in Table $(\ref{400})$. Referring to Table$(\ref{400}),$ the error norms $\mathit{L}_{2}$ and $\mathit{L}_{\infty}$ remain less than $1.286582\times 10^{-2}$, $8.31346\times 10^{-3}$ and they are still small when the time is increased up to $t=20$. The invariants $I_{1},I_{2}$, $I_{3}$ change from their initial values by less than $9.8\times\ 10^{-6},$ $3.2\times\ 10^{-5}$ and $1.3\times\ 10^{-5},$ respectively, throughout the simulation. Also, this table confirms that the changes of the invariants are in agreement with their exact values. So we conclude that our method is sensibly conservative. Comparisons with our results with exact solution as well as the calculated values in [8, 31, 32] have been made and showed in Table$(\ref{401})$ at $t=20$. This table clearly shows that the error norms got by our method are marginally less than the others. The numerical solutions at different time levels are depicted in Fig. $(\ref{500}).$ This figure shows that single soliton travels to the right at a constant speed and conserves its amplitude and shape with increasing time unsurprisingly. Initially, the amplitude of solitary wave is $1.00000$ and its top position is pinpionted at $x=30$. At $t=20,$ its amplitude is noted as $0.999416$ with center $x=40$. Thereby the absolute difference in amplitudes over the time interval $[0,20]$ are observed as $5.84\times 10^{-4}$. The quantile of error at discoint times are depicted in Fig.$(\ref{501})$ . The error aberration varies from $-8\times 10^{-2}$ to $1\times 10^{-2}$ and the maximum errors happen around the central position of the solitary wave.

For our second experiment, we take the parameters $p=3,$ $c=0.3,$ $h=0.1,$ $\Delta t=0.2,$ $\varepsilon=3,$ $\mu=1$, $x_{0}$ $=30$ with interval $[0,80]$ to coincide with that of previous papers [8, 31, 32]. Thus the solitary wave has amplitude $1.0$ and the computations are carried out for times up to $t=20.$ The values of the error norms $L_{2},$ $L_{\infty}$ and conservation quantities $I_{1},I_{2}$,$I_{3}$ are found and tabulated in Table $(\ref{402})$. According to Table$(\ref{402})$ the error norms $\mathit{L}_{2}$ and $\mathit{L}_{\infty}$ remain less than $4.48357\times 10^{-3}$, $3.37609\times 10^{-3}$ and they are still small when the time is increased up to $t=20$ and the invariants $I_{1},I_{2}$,$I_{3}$ change from their initial values by less than $1.78\times\ 10^{-5},$ $2.52\times\ 10^{-5}$, $3.55\times\ 10^{-5},$ respectively. Therefore we can say our method is satisfactorily conservative. In Table$(\ref{403})$ the performance of the our new method is compared with other methods [8, 31, 32] at $t=20$. It is observed that errors of the method [8, 31, 32] are considerably larger than those obtained with the present scheme. The motion of solitary wave using our scheme is graphed at time $t=0,10,20$ in Fig.$(\ref{502}).$ As seen, single solitons move to the right at a constant speed and preserves its amplitude and shape with increasing time as anticipated. The amplitude is $1.00000$ at $t=0$ and located at $x=30$, while it is $0.999522$ at $t=20$ and located at $x=36$. Therefore the absolute difference in amplitudes over the time interval $[0,20]$ are found as $4.78\times 10^{-4}$. The aberration of error at discrete times are drawn in Fig.$(\ref{503}).$ The error deviation varies from $-3\times 10^{-3}$ to $4\times 10^{-3}$ and the maximum errors arise around the central position of the solitary wave.

For our final treatment, we put the parameters $p=4,$ $c=0.2,$ $h=0.1,$ $\Delta t=0.2,$ $\varepsilon=3,$ $\mu=1$, $\ x_{0}$ $=30$ over the interval $[0,80]$ to make possible comparisons with those of earlier papers [8, 31, 32]. So solitary wave has amplitude $1.0$ and the simulations are executed to time $t=20$ 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 listed in Table$(\ref{4040})$ for several values of the time level $t$. It can be referred from Table$(\ref{4040}),$ the error norms $\mathit{L}_{2}$ and $\mathit{L}_{\infty}$ remain less than $1.96046\times 10^{-3}$, $1.33416\times 10^{-3}$ and they are still small when the time is increased up to $t=20$ and the invariants $I_{1},I_{2}$, $I_{3}$ change from their initial values by less than $4.07\times\ 10^{-5},$ $5.80\times\ 10^{-5}$ and $6.32\times\ 10^{-5},$ respectively, throughout the simulation. Hence we can say our method is sensibly conservative. The comparison between the results obtained by the current method with those in the other papers [8, 31, 32] is also documented in Table$(\ref{4050})$. It is noticeably seen from the table that errors of the current method are radically less than those obtained with the earlier methods [8, 31, 32]. For visual representation, the simulations of single soliton for values $p=4,c=0.2,h=0.1,\Delta t=0.2$ at times $t=0,10$ and $20$ 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=30$, while it is $0.999475$ at $t=20$ and located at $x=34$. The absolute difference in amplitudes at times $t=0$ and $t=10$ is $5.25\times 10^{-4}$ so that there is a little change between amplitudes. Error distributions at time $t=20$ are shown graphically in Figure$(\ref{505})$. As it is seen, the maximum errors are between $-1.5\times 10^{-3}$ to $1.5\times 10^{-3}$ and occur around the central position of the solitary wave.

6.2 Interaction of two solitary waves

Our second test problem pertains to the interaction of two solitary wave solutions of GEW equation having different amplitudes and traveling in the same direction. We tackle GEW equation with initial conditions given by the linear sum of two well separated solitary waves of various amplitudes as follows

[TABLE]

where $c_{j}$ and $x_{j}$, $\ j=1,2$ are arbitrary constants. For the computational work, two sets of parameters are considered by taking different values of $p,c_{i}$ and the same values of $h=0.1$, $\Delta t=0.025,$ $\varepsilon=3,$ $\mu=1$ over the interval $0\leq x\leq 80.$ We firstly take $p=3,$ $c_{1}=0.3,$ $c_{2}=0.0375.$ So the amplitudes of the two solitary waves are in the ratio $2:1.$ Calculations are done up to $t=100$. The three invariants in this case are tabulated in Table$(\ref{4051})$ . It is clear that the quantities are satisfactorily constant and very closed with the methods [8, 31, 32] during the computer run. Fig. $(\ref{40510})$ illustrates the behavior of the interaction of two positive solitary waves. At $t=100$, the magnitude of the smaller wave is $0.510619$ on reaching position $x=31.8$, and of the larger wave $0.999364$ having the position $x=46.7$, so that the difference in amplitudes is $0.010619$ for the smaller wave and $0.000636$ for the larger wave. For the second case, we have studied the interaction of two solitary waves with the parameters  $p=4,c_{1}=0.2,$ $c_{2}=1/80.$ So the amplitudes of the two solitary waves are in the ratio $2:1$.  For this case the experiment is run until time $t=120$. The three invariants in this case are recorded in Table$(\ref{40530}).$ The results in this table indicate that the numerical values of the invariants are good agreement with those of methods[8, 31, 32] during the computer run. Fig.$(\ref{4052})$ shows the development of the solitary wave interaction.

6.3 Evolution of solitons

Finally, another attracting initial value problem for the GEW equation is evolution of the solitons that is used as the Gaussian initial condition in solitary waves given by

[TABLE]

Since the behavior of the solution depends on values of $\mu$, we choose different values of $\mu=0.1$ and $\mu=0.05$ for $p=2,3,4$. The numerical computations are done up to $t=12$. Calculated numerical invariants at different values of $t$ are documented in Table$(\ref{4053}).$ From this table, we can easily see that as the value of $\mu$ increases, the variations of the invariants become smaller and it is seen that calculated invariant values are satisfactorily constant. The development of the evolution of solitons is presented in Figs.$(\ref{4055}),(\ref{4056})$ and $(\ref{4057})$. It is clearly seen in these figures that when the value of $\mu$ decreases, the number of the stable solitary wave increases.

7 Concluding remarks

. Solitary-wave solutions of the GEW equation by using Petrov-Galerkin method based on linear B-spline weight functions and quadratic B-spline trial functions, have been successfully obtained.

. Existence and uniqueness of solutions of the weak of the given problem as well as the proof of convergence has been proposed.

. Solutions of a semi-discrete finite element formulation of the equation and the theoretical bound of the error in the semi-discrete scheme are demonstrated.

. The theoretical upper bound of the error in such a full discrete approximation at $t=t^{n}$ has been proved.

. Our numerical algorithm has been tested by implementing three test problems involving a single solitary wave in which analytic solution is known and expanded it to investigate the interaction of two solitary waves and evolution of solitons where the analytic solutions are generally unknown during the interaction.

. The proffered method has been shown to be unconditionally stable.

. For single soliton the $L_{2}$ and $L_{\infty}$ error norms and for the three test problems the invariant quantities $I_{1}$, $I_{2}$ and $I_{3}$ have been computed. From the obtained results it is obviously clear that the error norms are sufficiently small and the invariants are marginally constant in all computer run. We can also see that our algorithm for the GEW equation is more accurate than the other earlier algorithms in the literature.

. Our method is an effective and a productive method to study behaviors of the dispersive shallow water waves.