# Runge-Kutta Integration of the Equal Width Wave Equation Using the Method of Lines

Authors: M. A. Banaja; H. O. Bakodah
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The equal width (EW) equation governs nonlinear wave phenomena like waves in shallow water. Numerical solution of the (EW) equation is obtained by using the method of lines (MOL) based on Runge-Kutta integration. Using von Neumann stability analysis, the scheme is found to be unconditionally stable. Solitary wave motion and interaction of two solitary waves are studied using the proposed method. The three invariants of the motion are evaluated to determine the conservation properties of the generated scheme. Accuracy of the proposed method is discussed by computing theL2 and L∞ error norms. The results are found in good agreement with exact solution.


---

## Body

## 1. Introduction

Nonlinear equations are widely used to describe complex phenomena in various fields of science, such as fluid mechanics, plasma physics, solid-state physics, elastodynamic problems, nonlinear Schrodinger equations, elasticity problems, free vibration analysis, and optical fibers [1–16]. Benjamin et al. [17] advocated that the partial differential equation (PDE) modeled the same physical phenomena equally well as the (KdV) equation given the same assumptions and approximations that were originally used by Korteweg and de Vries [18]. This PDE of Benjamin et al. [17] is now often called the BBM equation, although it is also known as the regularized long wave (RLW) equation. Morrison et al. [19] proposed the one-dimensional PDE as an equally valid and accurate model for the same wave phenomena simulated by the KdV and RLW equations. This PDE is called the equal width (EW) equation because the solutions for solitary waves with a permanent form and speed, for a given value of the parameter μ, are waves with an equal width or wavelength for all wave amplitudes. The equal width (EW) equation, which was introduced by Morrison et al. [19], is an important special kind of a nonlinear dispersive wave equation. It is defined as (1)ut+εuux-μuxxt=0.For a smooth functionu=u(x,t) on a domain Ω×0,T with Ω⊆R. Except for a single travelling solitary wave solution, no analytic solutions are known, and therefore numerical methods have to be used. The equation has solution of the form (2)ux,t=3csech2⁡px-x0-ct,p=12μ.That represents a solution for a solitary wave traveling at constant speedc and vanishing at ∓∞. When the EW equation is used to model waves generated in shallow water channel, the variables are normalized so that the distance x and water elevation u are scaled to the water depth h, and time t is scaled to h/g, where g is the acceleration due to gravity. There is experimental evidence to suggest that both descriptions break down if the amplitude of any wave exceeds 0.28 [20]. Different numerical solution methods were applied to solve EW equation such as L. R. Gardner and G. A. Gardner [21] using Galerkin’s method-based cubic B-spline finite elements, Ramos [22] using finite difference methods, Raslan [23] using quartic B-spline method, Saka [24] using finite element method, Ali et al. [25] using He’s Exp-function method, and Cheng and Liew [26] using the improved element-free Galerkin’s method.In this paper, the method of lines (MOL) solution of the EW equation is presented. This method consists of converting the EW equation with auxiliary conditions into a system of ordinary differential equations with corresponding auxiliary conditions. Then the fourth-order Runge-Kutta method is used to solve the system of first-order ordinary differential equations instead of finite difference methods [22] that are accurate and efficient as shown in [27].


## 2. The Method of Lines Solution of the EW Equation

The method of lines is a general technique for solving partial differential equations (PDEs) by typically using finite difference relationships for the spatial derivatives and ordinary differential equations for the time derivative. For time-dependent partial differential equations, meshless kernel-based methods were similarly based on a fixed spatial interpolation, but now the coefficients are time-dependent, and one obtains a system of ordinary differential equations for these. This is the well-known method of lines, and it turned to be experimentally useful in various cases.To apply the method of lines for solving the EW equation, firstly we subdivide the solution domain into uniform rectangular meshes by the lines(3)xi=ihi=0,1,2,…,N,h=b-aNtj=jkj=0,1,2,…,M,T=Mk.In numerical calculation, Ω=R is usually replaced by Ω=[a,b] and with zero boundary condition at both ends. Hence, the solution domain of the EW equation is the rectangle defined as a≤x≤b, 0≤t≤T. Then, the partial derivatives depending on spatial variables, ∂u/∂x and ∂2u/∂x2, in (1) are replaced by the well-known finite difference approximation at point xi:(4)∂ui∂x=ui+1-ui-12h∂2ui∂x2=ui+1-2ui+ui-1h2.By substituting (4) into (1) and introducing the boundary conditions u0t=uN+1t=0, this yields a system of ordinary differential equations which depend on t in the following form: (5)AdUidt=FUi,i=1,2,…,N,where(6)A=-h2+2μμ000μ-h2+2μμ⋯000μ-h2+2μ…00⋮⋱⋮00μ…μ0...⋯0-h2+2μμ0000μ-h2+2μ,FUj=εh2U1U2U2U3-U2U1UiUi+1-UiUi-1UN-1UN-2.Thus, we have the system of differential equations of one independent variable t. This system can be solved using Runge-Kutta method. The resulting system of ordinary differential equations is integrated with respect to time.Many existing numerical schemes for evolutionary problems in partial differential equations (PDEs) can be viewed as method of lines (MOL) schemes. In 1984, Verwer and Sanz-Serna [28] treated the convergence of one-step MOL schemes. Their main purpose was to set up a general framework for a convergence analysis applicable to nonlinear problems.


## 3. Stability Analysis

By solving EW equation using the method of lines, firstly the spatial variables are discretized to obtain an ordinary differential equations system in the time variable that can be discretized by an ordinary differential equation (ODE) solver. The simplest ODE solver is the forward Euler method and it is used widely for analysis of the stability properties of the spatial discretization. However, while forward Euler is ideal for analysis of the stability properties of a given spatial discretization, it is only first-order accurate. In practice, high order time discretization that preserves all the stability properties of forward Euler is needed. In [29], high order strong stability preserving time discretization method for the semidiscrete method of lines approximations of PDE is developed. These methods are derived by assuming that the first-order forward Euler time discretization of the method of lines ODE is strongly stable under a certain norm when the time step Δt is suitably restricted, and then try to find a higher order time discretization (Runge-Kutta or multistep) that maintains strong stability for the same norm perhaps under a different time-step restriction.The stability analysis constitutes the essential study of the numerical solution of PDEs. In general, this is because such study provides the means by which the step size and the numerical integration scheme for the given differential equation could be selected so as to secure manageable numerical solution.The stability analysis is based on the Neumann theory in which the growth factor of the error in a typical mode of amplitudee^n is(7)Ein=e^nejKih,j=-1,ih=x,where K is a Fourier number and h, the finite difference size, is determined from a linearization of the numerical scheme. Assuming u in the nonlinear term as a constant γ=maxi⁡ui; this enables one to discuss the stability in the linearized sense. The numerical method of lines of EW equation (8)ut+εuux-μuxxt=0gives the system of ordinary differential equations (9)μdui-1dt-h2+2μduidt+μdui+1dt=hε2uiui+1-ui-1,i=1,2,…,N.A trial solution is assumed and substituted into (9). The trial solution must take into account the variation of u(x,t) with both x and t or i and t. So (10)ux,t=cψtφx.In accordance with a method proposed by von Neumann, the function φx can be of the following form: (11)φx=ejKx,j=-1,x=ih.Substituting (10) and (11) into (9) gives (12)dψdt=-jhεγsin⁡Kh4μsin2⁡Kh/2+h2ψ.Equation (12) shows the growth factor for the error of the following form: (13)g=-jhεγsin⁡Kh4μsin2⁡Kh/2+h2.In [30], Evans and Raslan proved that g≤1 and the scheme is virtually unconditionally stable, at least for any practical problem, where h is a small quantity and γ represents the single speed and will usually be around unity.Kreiss and Scherer in [31] derived the conditions of local stability of Runge-Kutta methods when applied to hyperbolic partial differential equations, as when the time in (12) is discretized by using a locally stable Runge-Kutta method, the resulting completely discretized method is stable provided that kg≤R with locally stable Runge-Kutta methods whose stability region contains a half circle. Consider (14)σw≤R,Realσ≤0,where k denotes the time step and for computational purposes the Runge-Kutta method is only useful if it is stable for sufficient small k/h.


## 4. Test Problems

In this section, the results of the numerical solution of the EW equation based on the MOL using Runge-Kutta time discretization are presented. The numerical simulation includes the propagation of a solitary wave, the interaction of two solitary waves, the development of an undular bore, and the temporal evaluation of a Maxwellian initial pulse.In order to confirm the accuracy and efficiency of the method,L2 and L∞ error norms are used and defined by (15)L2=h∑i=0Nuei-uai21/2L∞=maxi⁡uei-uai,where ue denotes the exact solution and ua denotes the numerical solution. According to [32], in many cases, there are the following time invariants:(16)C1=∫abux,tdxC2=∫abu2x,t+μux2x,tdxC3=∫abu3x,tdx.By sufficiently fine spatial resolution, there is no problem to maintain these invariants to reasonable accuracy. Wherever possible, this statement is supported by providing numerical results. The constants of the motion, for a solitary wave of amplitude 3c and depending on k as given by (2), may be evaluated analytically to give (17)C1=6cp,C2=12c2p+48pc2μ5,C3=144c35p.


### 4.1. Single Solitary Wave

We first study the motion of a single solitary wave. This is derived from the initial condition. Consider(18)ux,0=3csech⁡2⁡px-x0,where p=1/(2μ) and c is a constant. This follows from the analytic solution (2) of the EW equation. We choose μ=1.0, c=1.0, x0=15, ε=1, h=0.2, and k=0.2 through the interval [0,80]. In this case the problem reduces to(19)ut+uux-uxxt=0(20)ux,0=3sech2⁡0.5x-15.To apply the method, let(21)xi=0.2ii=0,1,2,…,400tj=0.2jj=0,1,2,…,200.By substituting (4) into (19) and introducing the boundary conditions u0t=uN+1t=0, this yields a system of ordinary differential equations depending on t in the following form: (22)dUidt=A-1FUi,i=1,2,…,N,where A-1 is the inverse of the matrix A and (23)A=-2.0410001-2.041⋯0001-2.04…00⋮⋱⋮001…10...⋯-2.0410001-2.04.This system can be solved using Runge-Kutta method. The solution for times from t=0.0 and 40.0 is shown in Figure 1.Figure 1
Solitary wave amplitude 3.0 att=0.0 and t=40.In Table1, we examine various space/time step combinations and compare with simulation results given in [33]. The most accurate simulation is obtained in the case of h=k=0.04 for which L2 error has a value of 7.5×10-3.Table 1
Solitary wave amplitude 3.0,x0=15, 0≤x≤80, simulation results at t=40, and various space/time step combinations.
h
Δ
t
Present method
[26]
C
1
C
2
C
3
L
2
C
1
C
2
C
3
L
2
0.2
0.2
12.00000
28.7868
57.5944
0.386253
12.0474
28.9512
58.0881
0.218173
0.1
0.1
12.00000
28.7971
57.5998
0.0965104
12.0079
28.8247
57.6826
0.068558
0.05
0.05
12.00000
28.7993
57.6
0.0241263
11.9977
28.7842
57.5548
0.037365
0.04
0.04
12.00000
28.7995
57.6
0.0154409
11.9935
28.8177
57.6549
0.007512The simulation is run to timet=40 and the three invariants C1, C2, and C3 whose analytical values can be obtained as C1=12.0, C2=28.8, and C3=57.6 are listed for the duration of the simulation. We found that the MOL with Runge-Kutta integration is more accurate than Galerkin’s method with linear element in studying solitary wave of amplitude 3.0 for h=k=0.04.A solitary wave of amplitude 0.3 has also been modeled, and the results of the simulation are given in Table2. By the time t=80, L2 and L∞ error norms remain less than 1.8×10-3. The profiles of the solitary waves at time t=0 and t=80 are compared in Figure 2.Table 2
Solitary wave amplitude 0.3,x0=10, h=0.03, Δt=0.05, and 0≤x≤30.
Time
Present method
[26]
C
1
C
2
C
3
L
2
L
∞
C
1
C
2
C
3
L
2
L
∞
10
1.20001
0.287997
0.0576
0.0000332973
0.0000334152
1.20446
0.28953
0.05806
0.001647
0.001004
20
1.20004
0.287997
0.0576
0.0000565609
0.0000457087
1.20890
0.29102
0.05851
0.003556
0.002395
30
1.20005
0.287997
0.0576
0.0000777466
0.0000502313
1.21325
0.29244
0.05894
0.005813
0.004040
40
1.20005
0.287997
0.0576
0.0000987841
0.0000518949
1.21737
0.29377
0.05934
0.008496
0.005939
50
1.20005
0.287997
0.0576
0.000119954
0.0000613023
1.22173
0.29523
0.05978
0.011760
0.008165
60
1.20005
0.287997
0.0576
0.000141203
0.0000729713
1.22581
0.29655
0.06018
0.015532
0.010651
70
1.20005
0.287997
0.0576
0.000162472
0.0000844715
1.22988
0.29786
0.06058
0.019843
0.013401
80
1.20004
0.287997
0.0576
0.000183785
0.000095878
1.23387
0.29915
0.06097
0.024697
0.016425Figure 2
Solitary wave amplitude 3.0 att=0.0 and t=80.In further simulation of solitary waves of smaller magnitude, (2) is taken as initial condition with x0=10.0, c=0.03, and c=0.01 that the solitary waves have amplitudes 0.09 and 0.03. The simulations are run to time t=80 and L2 and L∞ error norms and the invariants C1, C2, and C3 are recorded throughout the simulation; see Tables 3 and 4.Table 3
Solitary wave amplitude 0.09,x0=10, h=0.05, Δt=0.05, and    0≤x≤30.
Time
Present method
[26]
C
1
C
2
C
3
L
2
L
∞
C
1
C
2
C
3
L
2
L
∞
10
0.359991
0.0259194
0.0015552
6.4777 × 10−6
4.02862 × 10−6
0.36081
0.02600
0.00156
0.000285
0.0001520
20
0.359998
0.0259194
0.0015552
0.0000126282
7.01317 × 10−6
0.36164
0.02608
0.00157
0.000579
0.000329
30
0.360002
0.0259194
0.0015552
0.0000185759
9.41599 × 10−6
0.36247
0.02616
0.00158
0.000884
0.000529
40
0.360006
0.0259194
0.0015552
0.0000243984
0.00001254
0.36330
0.02625
0.00159
0.001204
0.000750
50
0.360008
0.0259194
0.0015552
0.0000301437
0.0000156344
0.36413
0.02633
0.00159
0.001541
0.000993
60
0.36001
0.0259194
0.0015552
0.0000358406
0.0000186939
0.36497
0.02641
0.00160
0.001898
0.001255
70
0.360012
0.0259194
0.0015552
0.0000415067
0.0000217177
0.36581
0.02649
0.00161
0.002278
0.001536
80
0.360013
0.0259194
0.0015552
0.0000471519
0.0000247077
0.36665
0.02658
0.00162
0.002683
0.001836Table 4
Solitary wave amplitude 0.03,x0=10, h=0.05, Δt=0.05, and 0≤x≤30.
Time
Present method
[26]
C
1
C
2
C
3
L
2
L
∞
C
1
C
2
C
3
L
2
L
∞
10
0.119995
0.00287993
0.0000576
7.35307 × 10−7
4.93052 × 10−7
0.12009
0.00288
0.000058
0.000039
0.000030
20
0.119996
0.00287993
0.0000576
1.45427 × 10−6
9.39189 × 10−7
0.12021
0.00289
0.000058
0.000083
0.000061
30
0.119997
0.00287993
0.0000576
2.15926 × 10−6
1.34287 × 10−6
0.12033
0.00289
0.000058
0.000130
0.000091
40
0.119998
0.00287993
0.0000576
2.85232 × 10−6
1.70815 × 10−6
0.12045
0.00290
0.000058
0.000174
0.000119
50
0.119999
0.00287993
0.0000576
3.5352 × 10−6
2.03866 × 10−6
0.12056
0.00290
0.000058
0.000215
0.000144
60
0.119999
0.00287993
0.0000576
4.20941 × 10−6
2.33772 × 10−6
0.12067
0.00291
0.000058
0.000254
0.000167
70
0.120000
0.00287993
0.0000576
4.87624 × 10−6
2.60833 × 10−6
0.12078
0.00291
0.000058
0.000292
0.000188
80
0.120000
0.00287993
0.0000576
5.53677 × 10−6
2.85319 × 10−6
0.12088
0.00291
0.000059
0.000330
0.000206


### 4.2. Two Solitary Waves

As a second test problem for the EW equation, we chose the interaction of solitary waves, as mentioned in the used literature. The initial function is(24)ux,t=u1+u2,where (25)ui=3cisech2⁡0.5x-xi-ci,i=1,2,and solved the EW equation over the region 0≤x≤80 considering c1=1.5 and c2=0.75, x1=10, x2=25, h=0.1, and k=0.1. The shape of two solitary waves with given parameters is illustrated at the times t=0.0 and 30 in Figure 3. As seen in Figure 3, two solitary waves at the time t=0.0 are propagated to the right with velocities dependent upon their magnitudes and reached a stage where the larger solitary wave has passed through the smaller solitary wave and emerged in their original position. The values of the three invariants obtained by the present methods for this numerical experiment are recorded in Table 5. They are found to be very close to values given in [32] and the analytical ones given by C1=12c1+c2=27, C2=28.8c12+c22=81, and C3=57.6c13+c23=218.7. Additionally, we have observed that C3 remains almost constant, while C1 and C2 are affected more from the interaction of the two solitary waves during run of the algorithm.Table 5
Invariants for the interaction of two solitary waves.
Time
Present method
[32]
C
1
C
2
C
3
C
1
C
2
C
3
5
27.0002
80.9921
218.701
27.00001
81.00029
218.70226
10
27.0002
80.9892
218.679
27.00000
80.94427
218.66114
15
27.0002
80.9648
218.488
27.00055
80.94396
218.32641
20
27.0002
80.9874
218.67
27.00011
80.99415
218.65843
25
27.0002
80.9909
218.697
27.00044
80.00148
218.70712
30
27.0002
80.9913
218.698
27.00024
80.00140
218.70694Figure 3
Interaction of two solitary waves.


### 4.3. The Undular Bore

When a deeper stream of water flows into an area of the still water in a long horizontal channel, a bore is formed. To study development of an undular bore followed earlier by Peregrine [20], the following initial condition is used:(26)ux,0=0.5u01-tanh⁡x-x0d,where ux,0 denotes the evaluation of the water level above the equilibrium surface at time t=0. The change in the water level of magnitude u0 is centered on x=x0 and d measures the steepness of the change. We insistently choose the parameters d=5, μ=0.16666667, u0=0.1, and x0=0.0 in the region -20≤x≤50 to make a comparison with earlier works. The simulation is run until time t=800 and the values of the quantities C1, C2, and C3 are recorded in Table 6 with the values given in [34]. Figure 4 shows the undular bore profiles at times t=0, t=200, t=400, t=600, and t=800 for the gentle slope d=5. For the steep slope d=2, the results are compared with those given in [34] in Table 7.Table 6
Invariant for the undular bore withd=5.
Time
Present method
[34]
C
1
C
2
C
3
C
1
C
2
C
3
100
2.49294
0.24107
0.0236792
2.503585
0.242143
0.023787
200
2.99292
0.307831
0.0312025
3.003583
0.308783
0.031279
300
3.49292
0.374669
0.0387447
3.503573
0.375402
0.038765
400
3.99292
0.441516
0.0462886
4.003564
0.442020
0.046250
500
4.49292
0.508363
0.0538328
4.503559
0.508638
0.053735
600
4.99292
0.575211
0.0613769
5.003568
0.575259
0.061221
700
5.49292
0.64206
0.0689211
5.503578
0.641880
0.068707
800
5.99292
0.708908
0.0764653
6.003578
0.708498
0.076192Table 7
Invariant for the undular bore withd=2.
Time
Present method
[34]
C
1
C
2
C
3
C
1
C
2
C
3
100
2.49300
0.256324
0.025949
2.503504
0.257274
0.026029
200
2.99300
0.32317
0.0334927
3.003506
0.323895
0.033514
300
3.49300
0.390018
0.0410368
3.503511
0.390514
0.041000
400
3.99300
0.456867
0.048581
4.003510
0.457133
0.048485
500
4.49300
0.523715
0.0561252
4.503502
0.523750
0.055970
600
4.99300
0.590564
0.0636694
5.003488
0.590365
0.063455
700
5.49300
0.657413
0.0712136
5.503493
0.656985
0.070940
800
5.99298
0.724261
0.0787579
6.003478
0.723605
0.078426Figure 4
Undulation profiles ford=5 with different times t=200, t=400, t=600, and t=800.


### 4.4. The Maxwellian Initial Condition

The evaluation of an initial Maxwellian pulse into solitary waves is examined, using an initial condition in the following form [21]: (27)ux,0=exp⁡-x-72.For k=0.01, h=0.05, ε=1, and μ=0.001,0.04, and 0.2 through the interval [-20,20] and the values of the quantities C1, C2, and C3 are given in Table 8.Table 8
Invariant for the Maxwellian initial condition.
μ
Time
Present method
C
1
C
2
C
3
0.04
0.5
1.77245
1.30349
1.02348
1
1.77245
1.30375
1.02398
2.0
1.77245
1.30494
1.02627
3.0
1.77245
1.30681
1.02977
4.0
1.77245
1.30822
1.03238
0.2
0.5
1.77245
1.50384
1.02334
1
1.77245
1.50388
1.02340
2.0
1.77245
1.50401
1.02357
3.0
1.77245
1.50415
1.02374
4.0
1.77245
1.50422
1.02383
0.001
0.5
1.77245
1.25472
1.02373
1
1.77245
1.25574
1.02638
2.0
1.77245
1.30886
1.17932
3.0
1.77245
1.3841
1.37128
4.0
1.77245
1.42333
1.45041


## 4.1. Single Solitary Wave

We first study the motion of a single solitary wave. This is derived from the initial condition. Consider(18)ux,0=3csech⁡2⁡px-x0,where p=1/(2μ) and c is a constant. This follows from the analytic solution (2) of the EW equation. We choose μ=1.0, c=1.0, x0=15, ε=1, h=0.2, and k=0.2 through the interval [0,80]. In this case the problem reduces to(19)ut+uux-uxxt=0(20)ux,0=3sech2⁡0.5x-15.To apply the method, let(21)xi=0.2ii=0,1,2,…,400tj=0.2jj=0,1,2,…,200.By substituting (4) into (19) and introducing the boundary conditions u0t=uN+1t=0, this yields a system of ordinary differential equations depending on t in the following form: (22)dUidt=A-1FUi,i=1,2,…,N,where A-1 is the inverse of the matrix A and (23)A=-2.0410001-2.041⋯0001-2.04…00⋮⋱⋮001…10...⋯-2.0410001-2.04.This system can be solved using Runge-Kutta method. The solution for times from t=0.0 and 40.0 is shown in Figure 1.Figure 1
Solitary wave amplitude 3.0 att=0.0 and t=40.In Table1, we examine various space/time step combinations and compare with simulation results given in [33]. The most accurate simulation is obtained in the case of h=k=0.04 for which L2 error has a value of 7.5×10-3.Table 1
Solitary wave amplitude 3.0,x0=15, 0≤x≤80, simulation results at t=40, and various space/time step combinations.
h
Δ
t
Present method
[26]
C
1
C
2
C
3
L
2
C
1
C
2
C
3
L
2
0.2
0.2
12.00000
28.7868
57.5944
0.386253
12.0474
28.9512
58.0881
0.218173
0.1
0.1
12.00000
28.7971
57.5998
0.0965104
12.0079
28.8247
57.6826
0.068558
0.05
0.05
12.00000
28.7993
57.6
0.0241263
11.9977
28.7842
57.5548
0.037365
0.04
0.04
12.00000
28.7995
57.6
0.0154409
11.9935
28.8177
57.6549
0.007512The simulation is run to timet=40 and the three invariants C1, C2, and C3 whose analytical values can be obtained as C1=12.0, C2=28.8, and C3=57.6 are listed for the duration of the simulation. We found that the MOL with Runge-Kutta integration is more accurate than Galerkin’s method with linear element in studying solitary wave of amplitude 3.0 for h=k=0.04.A solitary wave of amplitude 0.3 has also been modeled, and the results of the simulation are given in Table2. By the time t=80, L2 and L∞ error norms remain less than 1.8×10-3. The profiles of the solitary waves at time t=0 and t=80 are compared in Figure 2.Table 2
Solitary wave amplitude 0.3,x0=10, h=0.03, Δt=0.05, and 0≤x≤30.
Time
Present method
[26]
C
1
C
2
C
3
L
2
L
∞
C
1
C
2
C
3
L
2
L
∞
10
1.20001
0.287997
0.0576
0.0000332973
0.0000334152
1.20446
0.28953
0.05806
0.001647
0.001004
20
1.20004
0.287997
0.0576
0.0000565609
0.0000457087
1.20890
0.29102
0.05851
0.003556
0.002395
30
1.20005
0.287997
0.0576
0.0000777466
0.0000502313
1.21325
0.29244
0.05894
0.005813
0.004040
40
1.20005
0.287997
0.0576
0.0000987841
0.0000518949
1.21737
0.29377
0.05934
0.008496
0.005939
50
1.20005
0.287997
0.0576
0.000119954
0.0000613023
1.22173
0.29523
0.05978
0.011760
0.008165
60
1.20005
0.287997
0.0576
0.000141203
0.0000729713
1.22581
0.29655
0.06018
0.015532
0.010651
70
1.20005
0.287997
0.0576
0.000162472
0.0000844715
1.22988
0.29786
0.06058
0.019843
0.013401
80
1.20004
0.287997
0.0576
0.000183785
0.000095878
1.23387
0.29915
0.06097
0.024697
0.016425Figure 2
Solitary wave amplitude 3.0 att=0.0 and t=80.In further simulation of solitary waves of smaller magnitude, (2) is taken as initial condition with x0=10.0, c=0.03, and c=0.01 that the solitary waves have amplitudes 0.09 and 0.03. The simulations are run to time t=80 and L2 and L∞ error norms and the invariants C1, C2, and C3 are recorded throughout the simulation; see Tables 3 and 4.Table 3
Solitary wave amplitude 0.09,x0=10, h=0.05, Δt=0.05, and    0≤x≤30.
Time
Present method
[26]
C
1
C
2
C
3
L
2
L
∞
C
1
C
2
C
3
L
2
L
∞
10
0.359991
0.0259194
0.0015552
6.4777 × 10−6
4.02862 × 10−6
0.36081
0.02600
0.00156
0.000285
0.0001520
20
0.359998
0.0259194
0.0015552
0.0000126282
7.01317 × 10−6
0.36164
0.02608
0.00157
0.000579
0.000329
30
0.360002
0.0259194
0.0015552
0.0000185759
9.41599 × 10−6
0.36247
0.02616
0.00158
0.000884
0.000529
40
0.360006
0.0259194
0.0015552
0.0000243984
0.00001254
0.36330
0.02625
0.00159
0.001204
0.000750
50
0.360008
0.0259194
0.0015552
0.0000301437
0.0000156344
0.36413
0.02633
0.00159
0.001541
0.000993
60
0.36001
0.0259194
0.0015552
0.0000358406
0.0000186939
0.36497
0.02641
0.00160
0.001898
0.001255
70
0.360012
0.0259194
0.0015552
0.0000415067
0.0000217177
0.36581
0.02649
0.00161
0.002278
0.001536
80
0.360013
0.0259194
0.0015552
0.0000471519
0.0000247077
0.36665
0.02658
0.00162
0.002683
0.001836Table 4
Solitary wave amplitude 0.03,x0=10, h=0.05, Δt=0.05, and 0≤x≤30.
Time
Present method
[26]
C
1
C
2
C
3
L
2
L
∞
C
1
C
2
C
3
L
2
L
∞
10
0.119995
0.00287993
0.0000576
7.35307 × 10−7
4.93052 × 10−7
0.12009
0.00288
0.000058
0.000039
0.000030
20
0.119996
0.00287993
0.0000576
1.45427 × 10−6
9.39189 × 10−7
0.12021
0.00289
0.000058
0.000083
0.000061
30
0.119997
0.00287993
0.0000576
2.15926 × 10−6
1.34287 × 10−6
0.12033
0.00289
0.000058
0.000130
0.000091
40
0.119998
0.00287993
0.0000576
2.85232 × 10−6
1.70815 × 10−6
0.12045
0.00290
0.000058
0.000174
0.000119
50
0.119999
0.00287993
0.0000576
3.5352 × 10−6
2.03866 × 10−6
0.12056
0.00290
0.000058
0.000215
0.000144
60
0.119999
0.00287993
0.0000576
4.20941 × 10−6
2.33772 × 10−6
0.12067
0.00291
0.000058
0.000254
0.000167
70
0.120000
0.00287993
0.0000576
4.87624 × 10−6
2.60833 × 10−6
0.12078
0.00291
0.000058
0.000292
0.000188
80
0.120000
0.00287993
0.0000576
5.53677 × 10−6
2.85319 × 10−6
0.12088
0.00291
0.000059
0.000330
0.000206


## 4.2. Two Solitary Waves

As a second test problem for the EW equation, we chose the interaction of solitary waves, as mentioned in the used literature. The initial function is(24)ux,t=u1+u2,where (25)ui=3cisech2⁡0.5x-xi-ci,i=1,2,and solved the EW equation over the region 0≤x≤80 considering c1=1.5 and c2=0.75, x1=10, x2=25, h=0.1, and k=0.1. The shape of two solitary waves with given parameters is illustrated at the times t=0.0 and 30 in Figure 3. As seen in Figure 3, two solitary waves at the time t=0.0 are propagated to the right with velocities dependent upon their magnitudes and reached a stage where the larger solitary wave has passed through the smaller solitary wave and emerged in their original position. The values of the three invariants obtained by the present methods for this numerical experiment are recorded in Table 5. They are found to be very close to values given in [32] and the analytical ones given by C1=12c1+c2=27, C2=28.8c12+c22=81, and C3=57.6c13+c23=218.7. Additionally, we have observed that C3 remains almost constant, while C1 and C2 are affected more from the interaction of the two solitary waves during run of the algorithm.Table 5
Invariants for the interaction of two solitary waves.
Time
Present method
[32]
C
1
C
2
C
3
C
1
C
2
C
3
5
27.0002
80.9921
218.701
27.00001
81.00029
218.70226
10
27.0002
80.9892
218.679
27.00000
80.94427
218.66114
15
27.0002
80.9648
218.488
27.00055
80.94396
218.32641
20
27.0002
80.9874
218.67
27.00011
80.99415
218.65843
25
27.0002
80.9909
218.697
27.00044
80.00148
218.70712
30
27.0002
80.9913
218.698
27.00024
80.00140
218.70694Figure 3
Interaction of two solitary waves.


## 4.3. The Undular Bore

When a deeper stream of water flows into an area of the still water in a long horizontal channel, a bore is formed. To study development of an undular bore followed earlier by Peregrine [20], the following initial condition is used:(26)ux,0=0.5u01-tanh⁡x-x0d,where ux,0 denotes the evaluation of the water level above the equilibrium surface at time t=0. The change in the water level of magnitude u0 is centered on x=x0 and d measures the steepness of the change. We insistently choose the parameters d=5, μ=0.16666667, u0=0.1, and x0=0.0 in the region -20≤x≤50 to make a comparison with earlier works. The simulation is run until time t=800 and the values of the quantities C1, C2, and C3 are recorded in Table 6 with the values given in [34]. Figure 4 shows the undular bore profiles at times t=0, t=200, t=400, t=600, and t=800 for the gentle slope d=5. For the steep slope d=2, the results are compared with those given in [34] in Table 7.Table 6
Invariant for the undular bore withd=5.
Time
Present method
[34]
C
1
C
2
C
3
C
1
C
2
C
3
100
2.49294
0.24107
0.0236792
2.503585
0.242143
0.023787
200
2.99292
0.307831
0.0312025
3.003583
0.308783
0.031279
300
3.49292
0.374669
0.0387447
3.503573
0.375402
0.038765
400
3.99292
0.441516
0.0462886
4.003564
0.442020
0.046250
500
4.49292
0.508363
0.0538328
4.503559
0.508638
0.053735
600
4.99292
0.575211
0.0613769
5.003568
0.575259
0.061221
700
5.49292
0.64206
0.0689211
5.503578
0.641880
0.068707
800
5.99292
0.708908
0.0764653
6.003578
0.708498
0.076192Table 7
Invariant for the undular bore withd=2.
Time
Present method
[34]
C
1
C
2
C
3
C
1
C
2
C
3
100
2.49300
0.256324
0.025949
2.503504
0.257274
0.026029
200
2.99300
0.32317
0.0334927
3.003506
0.323895
0.033514
300
3.49300
0.390018
0.0410368
3.503511
0.390514
0.041000
400
3.99300
0.456867
0.048581
4.003510
0.457133
0.048485
500
4.49300
0.523715
0.0561252
4.503502
0.523750
0.055970
600
4.99300
0.590564
0.0636694
5.003488
0.590365
0.063455
700
5.49300
0.657413
0.0712136
5.503493
0.656985
0.070940
800
5.99298
0.724261
0.0787579
6.003478
0.723605
0.078426Figure 4
Undulation profiles ford=5 with different times t=200, t=400, t=600, and t=800.


## 4.4. The Maxwellian Initial Condition

The evaluation of an initial Maxwellian pulse into solitary waves is examined, using an initial condition in the following form [21]: (27)ux,0=exp⁡-x-72.For k=0.01, h=0.05, ε=1, and μ=0.001,0.04, and 0.2 through the interval [-20,20] and the values of the quantities C1, C2, and C3 are given in Table 8.Table 8
Invariant for the Maxwellian initial condition.
μ
Time
Present method
C
1
C
2
C
3
0.04
0.5
1.77245
1.30349
1.02348
1
1.77245
1.30375
1.02398
2.0
1.77245
1.30494
1.02627
3.0
1.77245
1.30681
1.02977
4.0
1.77245
1.30822
1.03238
0.2
0.5
1.77245
1.50384
1.02334
1
1.77245
1.50388
1.02340
2.0
1.77245
1.50401
1.02357
3.0
1.77245
1.50415
1.02374
4.0
1.77245
1.50422
1.02383
0.001
0.5
1.77245
1.25472
1.02373
1
1.77245
1.25574
1.02638
2.0
1.77245
1.30886
1.17932
3.0
1.77245
1.3841
1.37128
4.0
1.77245
1.42333
1.45041


## 5. Conclusion

The EW equation is numerically solved using MOL with Runge-Kutta integration. The numerical solution leads to an unconditionally stable algorithm with which accurate simulations of the motion of a solitary wave are found over an extended time scale. Single solitary waves of amplitude 3.0 forh=k=0.04 are found. MOL with Runge-Kutta integration is more accurate than Galerkin’s method with linear elements for a single solitary wave; it gives better conservation than that given in [33]. Moreover, for a smaller solitary wave amplitude = 0.03, excellent results are obtained.The performance of the method has been examined by studying the propagation of interaction of two solitary waves. Using the proposed method, the magnitude, profile, and position of solitary waves are faithfully retrieved. To make a comparison between results of MOL and some earlier published works, the same steph is used.The MOL based on Runge-Kutta integration was efficiently applied to an undular bore. The obtained results are good and accurate as compared to the results in [34]. In addition, the Maxwellian initial condition is simulated.The obtained results indicate that the present method is remarkably successful numerical technique for solving the EW equation. The results suggest that proposed method, whose application is easier than many other numerical schemed methods such as finite element methods and Galerkin’s method, can be applied to this type of the nonlinear problems with success. The method can be also used efficiently for solving a large number of physically important nonlinear problems.



---

*Source: 274579-2015-07-05.xml*