# The Method of Lines Solution of the Regularized Long-Wave Equation Using Runge-Kutta Time Discretization Method

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

---

## Abstract

A method of lines approach to the numerical solution of nonlinear wave equations typified by the regularized long wave (RLW) is presented. The method developed uses a finite differences discretization to the space. Solution of the resulting system was obtained by applying fourth Runge-Kutta time discretization method. Using Von Neumann stability analysis, it is shown that the proposed method is marginally stable. To test the accuracy of the method some numerical experiments on test problems are presented. Test problems including solitary wave motion, two-solitary wave interaction, and the temporal evaluation of a Maxwellian initial pulse are studied. The accuracy of the present method is tested withL∞ and L2 error norms and the conservation properties of mass, energy, and momentum under the RLW equation.


---

## Body

## 1. Introduction

Many researchers have introduced various methods to solve the regularized long wave (RLW) of the form(1)ut+ux+ϵuux-μuxxt=0,
where ϵ and μ are positive constants.The RLW equation was introduced to describe the behavior of the undular bore by Peregrine [1]. The RLW equation describes a lot of important physical phenomena, such as shallow water waves and ionic waves.An analytical solution for the RLW equation was found under the restricted initial and boundary condition [2], but this solution is not very useful; the availability of numerical method is essential. In previous work [3–5] various numerical studies have been reported based on the finite difference. A recent study in [6] has shown that the meshless kernel-based method of lines provides a highly accurate and efficient method for the modified regularized long wave equation. In this study, the method of lines (MOL) is applied to obtain a numerical solution for RLW with a constant grid of central finite differences and fourth Runge-Kutta approximation in time. The performance of the method was tested on two known model problems.


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

We shall consider (1) with the following initial and boundary conditions:
(2)u(x,0)=g(x),a≤x≤b,u(a,t)=0,u(b,t)=0,
where g(x) is prescribed function.The basic idea of the method of lines consists in a discretization of the space variable and in substitution of partial derivatives by finite-difference expressions; the time variable keeps a continuous form, and solutions above the lines parallel with the axis of this variable are computed.Let us subdivide the solution domain of the RLW equationa≤x≤b into uniform rectangular mesh by the lines
(3)xi=ih(i=0,1,2,…,N),whereh=b-aN.
The partial derivatives depending on spatial variable in RLW equation are replaced by known central finite difference approximations at point xi. This leads to a tridiagonal system of ordinary differential equation
(4)∂ui∂x=ui+1-ui-12h,(5)∂2ui∂x2=ui+1-2ui+ui-1h2.
By substituting (4) and (5) into (1) and taking into account that u0(t)=uN+1(t)=0 the following system of ODEs is obtained:
(6)-(h2+2μ)du1dt+μdu2dt=h2[u2+ϵ2u22]μdui-1dt-(h2+2μ)duidt+μdui+1dt=-h2[ui-1+ϵ2ui-12]+h2[ui+1+ϵ2ui+12]i=2,3,…,N-1,μduN-1dt-(h2+2μ)duNdt=-h2[uN-1+ϵ2uN-12].
This system can be written in the following form:
(7)AdUdt=BU,Ui(0)=g(xi),
where
(8)U(t)=[u1(t),u2(t),…,un(t)]T,Ui(0)=[u1(0),u2(0),…,un(0)]T,U(0)=[g(x1),g(x2),…,g(xn)]T,
where(9)A=[-(h2+2μ)μμ-(h2+2μ)μμ-(h2+2μ)⋯000000μ00⋮⋱⋮000...000⋯-(h2+2μ)μ0μ-(h2+2μ)μμ-(h2+2μ)],B=h2[u2+ϵ2u22,(u3+ϵ2u32)-(u1+ϵ2u12),…,(ui+1+ϵ2ui+12)-(ui-1+ϵ2ui-12),…,-(uN-1+ϵ2uN-12)]T.It can be easily solved by computing A-1; this yields a system of ordinary differential equations depending on the time variable t in the form
(10)dUidt=f(Ui),i=1,2,…,N.
This system can be solved using some numerical method; in this paper the solutions are obtained using Runge-Kutta method.


## 3. Stability Analysis of the Numerical Method

The stability of the method of lines for partial differential equations represents the most important factor for their solution; hence, it is a critical factor that should be handled carefully. The importance lies in its unique ability of judging acceptable solution for the given equation. The criticality depends on the nature of the eigenvalues of the matrix representation in connection with their values. Kreiss and Scherer in [10] derived the conditions of local stability of Runge-Kutta methods when applied to hyperbolic partial differential equations.Theorem 1.
Consider a well-posed initial boundary value problem for a system of partial differential equations. The space derivatives are discretized in such a way that the resulting system of ODEs(11)dudt=Qu
is stable. Then, the time is discretized by using a locally stable (implicit or explicit) Runge-Kutta method or a multistep scheme. The resulting completely discretized method is stable, provided that
(12)∥kQ∥<R
with locally stable Runge-Kutta methods. The stability region contains a half circle
(13)|σ(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 sufficiently small k/h.In an attempt to gain some insight into the stability of method (7), the standard Fourier analysis is used to determine the condition to be imposed on the time step Δt for stability. The numerical method of lines of RLW equation
(14)ut+ux+ϵuux-μuxxt=0
gives the system of ordinary differential equations
(15)μdui-1dt-(h2+2μ)duidt+μdui+1dt=h2[ui+1+ui-1+ϵui(ui+1-ui-1)],i=1,2,3,…,N.
First,  ui is assumed as constant in the nonlinear term, where  γ=maxiui.We can now inquire about the eigenvalue of theN ordinary differential equation (15). Thus, a trial solution is assumed and substituted into (15). However, the trial solution must take into account the variation of u(x,t) with both x and t or i and t. Therefore, a product solution is assumed:
(16)u(x,t)=cψ(t)φ(x).
Further, in accordance with a method proposed by Von Neumann, the function φ(x) can be of the form
(17)φ(x)=ejkx,j=-1,
where k is a Fourier number.Substituting of (16) and (17) in (15) gives
(18)c[μdψdtejkxi+1-(h2+2μ)dψdtejkxi+μdψdtejkxi-1]=h(1+ϵγ)2(ψejkxi+1-ψejkxi-1).
Equation (18) written in terms of eigenvalue λ becomes
(19)dψdt=h(1+ϵγ)isin(kΔx)-(4μsin2(kΔx/2)+h2)ψ,(20)λ=h(1+ϵγ)sin(kΔx)-(4μsin2(kΔx/2)+h2)i
with λ a complex scalar, and denote z=pλ for numerical discretization with step size p=Δt which can be written as
(21)ψn+1=ℛ(z)ψn,n=0,1,2,…,N.
The classical Runge-Kutta (RK4) method contains a segment of the imaginary axis in its stability region and is only mildly dissipative, specifically; |ℛ(z)| is only a little less than 1 for |Im(z)|≤1.2.Note that allN eigenvalues given by (18) have pure imaginary parts, which means that the system of N first-order ODEs are marginally stable.


## 4. Numerical Results

To illustrate the effectiveness of the method, numerical results portraying a single soliton solution for the RLW equation are presented. In order to show how good the numerical solutions are in comparison with exact ones, we will use theL2 and L∞ error norms defined by
(22)L2=[h∑i=0N|(ue)i-(ua)i|2]1/2,L∞=maxi|(ue)i-(ua)i|,
where ue denotes the exact solution and ua denotes the exact numerical solution. The conservation properties of the RLW equation that will be validated by computing quantities corresponding to mass, momentum, and energy [11], respectively, are
(23)C1=∫abudx≅h∑i=1Nuin,C2=∫ab(u2+μ(ux)2)dx,≅h∑i=1N‍((uN)in)2+μ(((ux)in)2)C3=∫ab(u3+3(u)2)dx,≅h∑i=1N‍((uN)in)3+3(((u)in)2).
In the following test problems, the numerical solutions must control these conservation laws during propagation. Therefore, these quantities are used to measure the accuracy of the present method. Integrals for the conservation laws (23) were computed approximately with the trapezium rule. Nodal values (uN)in and its first derivative ((uN)x)in can be computed from (4).It is well known that (1) has analytic solution of the form
(24)u(x,t)=3csech2[k(x-νt-x0)],
where ν=1+εc is the wave velocity, k=(1/2μ)ε/(1+εc) the width and initially centered at x0.


### 4.1. Motion of the Solitary Waves

Case 1.
The equation represents a single soliton with amplitude3c and constant speed c+1. For the purpose of comparing with the earlier work, the computations are done for the parameters ε=1, μ=1, and x0=0. Consider c=0.1, as it seen from Table 1, the numerical results of the present method are in very good agreement with their analytical values obtained from the exact solution. Moreover, Table 2 displays a comparison of the approximate solution by the present method with those obtained using MOL developed by Griffiths and Schiesser in [12] at the time t=18. The error norms at t=18 obtained by the present method are smaller than those given in [12]. Also, Table 3 displays a comparison of error norms obtained by the present method with those obtained by using finite difference solution [7] and finite element [9].Table 1
L
2 and L∞ errors for (1) with c=0.1  Δx=0.1,Δt=0.1 and over the region -40≤x≤60.
t
L
2
L
∞
4
0.69452 × 10−4
0.264411 × 10−4
8
0.137889 × 10−3
0.534151 × 10−4
12
0.204502 × 10−3
0.792604 × 10−4
16
0.268849 × 10−3
0.103317 × 10−3
20
0.33080 × 10−3
0.125611 × 10−3Table 2
Absolute errors  for (1) with c=0.1  n=199 and over the region -50≤x≤50,t=18.
x
Present method
MOL [7]
−49.5
3.44444 × 10−9
0.000000000
−39.5
3.01068 × 10−8
0.000000000
−29.5
5.43620 × 10−9
0.000000000
−19.5
2.01194 × 10−7
1.00000 × 10−6
−9.5
3.91714 × 10−6
11.0000 × 10−5
0.5
4.39242 × 10−5
2.03000 × 10−4
10.5
2.31414 × 10−5
3.59000 × 10−3
20.5
1.43949 × 10−4
1.94700 × 10−3
30.5
1.9499 × 10−4
2.38500 × 10−3
40.4
4.67298 × 10−5
1.33000 × 10−4Table 3
L
2 and L∞ errors for (1) with c=0.1  Δx=0.1,Δt=0.1 and over the region -40≤x≤60.
t
Present method
Finite difference
Finite element
L
2
L
∞
L
2
L
∞
L
2
L
∞
4
0.6945 × 10−4
0.26441 × 10−4
0.1200 × 10−3
0.05000 × 10−3
0.1176 × 10−3
0.5622 × 10−3
8
0.13788 × 10−3
0.53415 × 10−4
0.2300 × 10−3
0.09000 × 10−3
0.23108 × 10−3
0.1088 × 10−3
12
0.20450 × 10−3
0.79260 × 10−4
0.3400 × 10−3
0.1400 × 10−3
0.3374 × 10−3
0.1544 × 10−3
16
0.26884 × 10−3
0.10331 × 10−3
0.45000 × 10−3
0.18000 × 10−3
0.43784 × 10−3
0.19377 × 10−3
20
0.33080 × 10−3
0.12561 × 10−3
0.55000 × 10−3
0.21000 × 10−3
0.53223 × 10−3
0.2272 × 10−3The error norms at each time are smaller than those given in [7, 9]. The invariant values C1, C2, and C3 are given in Table 4 and compared with the results obtained by [8]. From this table, it can be seen that the conservation properties (using the present method) are very good in comparison with the analytical values for the invariants given by C1=6c/k=3.9799497, C2=(12c2/k)+(48c2μ/5)=0.81046249, and C3=(36c2/k)+(144c3/5k)=2.579007 [8].Table 4
Invariant values for a single solitary wave: amplitude = 0.3,Δx=0.125,Δt=0.1 and over the region  -40≤x≤60.
t
Present method
Cubic spline collocation method [8]
C
1
C
2
C
3
C
1
C
2
C
3
0
3.97990
0.810461
2.57901
3.97992
0.810462
2.57900
4
3.97995
0.810461
2.57901
3.97995
0.810276
2.57839
8
3.97997
0.810461
2.57901
3.97997
0.810276
2.57839
12
3.97998
0.810461
2.57901
3.97998
0.810276
2.57839
16
3.97999
0.810461
2.57901
3.97998
0.810276
2.57839
20
3.97996
0.810461
2.57901
3.97996
0.810276
2.57839Forc=0.1, the relative changes in the invariants are very small (the quantities C1, C2, and C3 change by less than 1.03×10-5, 1.49×10-6, and 3×10-6, respectively, by t=20). Initial solution and solitary wave profile at t=20 are depicted in Figure 1.Figure 1
Graphs of initial solution and solution at time = 20 and amplitude = 0.3.Case 2.
In the second experiment, a smaller solitary wave of amplitude 0.09 (c=0.03) has been modeled, and the results of simulation are given in Table 5. As the amplitude of the solitary wave is reduced the pulse broadens, and it may be necessary to increase the solution range in order to maintain accuracy. The range here is doubled from -40≤x≤60 to -80≤x≤120. The error norms for c=0.03 are recorded in Table 6 for different values of time t. With the range -40≤x≤60, Δt=0.1, excellent results were obtained, L2 and L∞ error norms remain less than L2=0.637×10-3 and L∞=0.233×10-3, whereas they are L2=0.638×10-3 and L∞=0.233×10-3 for the finite difference method [7] and L2=0.572×10-3 and L∞=0.364×10-3 for the finite element method [9]. Moreover, it can be seen that the error norms (L2 and L∞) are somewhat small compared with those quoted by others [7, 9].Table 5
L
2 and L∞ errors for (1) with c=0.03  Δx=0.1,Δt=0.1 and over the region -80≤x≤120.
t
L
2
L
∞
4
0.51431 × 10−5
0.14805 × 10−5
8
0.10281 × 10−4
0.29961 × 10−5
12
0.15408 × 10−4
0.45367 × 10−5
16
0.20518 × 10−4
0.60908 × 10−5
20
0.25606 × 10−4
0.76472 × 10−5Table 6
L
2 and L∞ errors for (1) with c=0.03  Δx=0.1,Δt=0.1 and over the region -40≤x≤60.
t
Present method
Finite difference
Finite element
L
2
L
∞
L
2
L
∞
L
2
L
∞
4
0.14956 × 10−3
0.12290 × 10−3
0.1500 × 10−3
0.1230 × 10−3
0.3983 × 10−3
0.22450 × 10−3
8
0.32000 × 10−3
0.16600 × 10−3
0.3200 × 10−3
0.1660 × 10−3
0.50528 × 10−3
0.21706 × 10−3
12
0.466159 × 10−3
0.17900 × 10−3
0.4670 × 10−3
0.1790 × 10−3
0.5357 × 10−3
0.20919 × 10−3
16
0.56635 × 10−3
0.18454 × 10−3
0.5670 × 10−3
0.1850 × 10−3
0.54934 × 10−3
0.20165 × 10−3
20
0.63748 × 10−3
0.23300 × 10−3
0.6380 × 10−3
0.2330 × 10−3
0.57247 × 10−3
0.36498 × 10−3The invariantsC1, C2,  and C3 for RLW at different values of time t are also recorded in Table 7, where it can be seen that the invariants remain constant with respect to time. In comparison with the analytical values C1=2.109407, C2=0.127302, and  C3=0.388806, it can be seen from Table 7 that the quantities C1, C2, and C3 obtained using present method are very close to the corresponding exact values and compared with the results obtained by [8].Table 7
Invariant values for a single solitary wave: amplitude = 0.09,Δx=0.125,Δt=0.1 and over the region  -40≤x≤60.
t
Present method
Cubic spline collocation method [8]
C
1
C
2
C
3
C
1
C
2
C
3
0
2.10700
0.127302
0.38806
2.10702
0.127302
0.388804
4
2.10836
0.127302
0.38807
2.10838
0.127294
0.388781
8
2.10927
0.127302
0.38807
2.10928
0.127294
0.388782
12
2.10980
0.127302
0.38807
2.10981
0.127294
0.388782
16
2.10981
0.127302
0.38807
2.10981
0.127294
0.388782
20
2.10899
0.127302
0.38806
2.10900
0.127294
0.388781Forc=0.03, the relative changes in the quantities C1, C2, and C3 are changed by less than 4.1×10-4, 1.5×10-7, and 3.1×10-7, respectively, by t=20. Initial solution and solitary wave profile at t=20 are depicted in Figure 2.Figure 2
Graphs of initial solution and solution at time = 20 and amplitude = 0.09.


### 4.2. The Interaction of Solitary Wave

The interaction of two positive solitary waves is studied by using the initial condition given by the linear sum of two separate solitary waves of various amplitudes(25)u(x.0)=u1+u2,
where ui=3cisech2[ki(x-xi)], i=1,2, and ci=4ki2/(1-4ki2), with ε=μ=1.Case 1.
In this case, the numerical example consists in the interaction of two positive solitary waves defined by the initial condition (25), and the boundary conditions u(a,t)=u(b,t)=0.The RLW parameters used in the numerical simulation arek1=0.4, k2=0.3, x1=15, x2=35, ε=μ=1, and Δt=0.01. This wave moves across the space interval 0≤x≤120 and the time interval 0≤t≤25.Figure3 shows the interaction of two soliton solutions of RLW equation. From this figure, it can be seen that the faster pulse interacts with, and emerges ahead of, the slower pulse, with the shape and velocity of each soliton retained. Numerical check on the conservation of mass, momentum, and energy shows that the three quantities remain constant with respect to time t (see Table 8). The computed values of the invariants C1, C2, and C3 are satisfactorily constant compared with the corresponding invariant values obtained by [9]. The results of the present method at different times are shown in Figure 3, where it is possible to observe that the higher amplitude solitary wave passes through the smaller wave with no change in its waveform.Table 8
Invariant values for the interaction of two positive solitary waves over the region0≤x≤120.
t
Present method
[9]
C
1
C
2
C
3
C
1
C
2
C
3
0
37.9165
120.480
744.081
37.916520
120.52280
744.0815
4
37.9170
120.471
743.984
37.911700
121.16020
736.9443
8
37.9172
120.445
743.676
37.859750
119.73170
728.51730
12
37.9173
120.203
740.702
37.696670
119.63340
725.72360
16
37.9174
120.230
741.050
37.529160
119.41850
725.83940
20
37.9174
120.441
743.654
37.647300
119.80410
727.19480
25
37.9175
120.462
743.918
38.050100
119.83550
727.439200Figure 3
The motion of two positive solitary waves before and after the interaction at time = 0 and 25.Case 2.
The interaction of two negative solitary waves with initial condition (25) was studied over the region 0≤x≤120 with k1=0.6, k2=0.8, x1=82, x2=67, ε=μ=1, and Δt=0.2, Δx=0.5. The invariants C1, C2, and C3 recorded in Table 9; the computed values are satisfactorily constant compared with the corresponding invariant values obtained by [9]. The two negative solitary waves before and after the interaction are shown in Figure 4.Table 9
Invariant values for the interaction of two negative solitary waves over the region0≤x≤120.
t
Present method
[9]
C
1
C
2
C
3
C
1
C
2
C
3
0
−45.035
336.232
−1077.82
−45.034960
336.96020
−1077.824
4
−45.035
336.328
−1077.97
−45.035000
334.25310
−1048.6640
8
−45.035
336.014
−1068.16
−45.035080
332.98740
−1041.3070
12
−45.035
336.172
−1077.67
−45.035100
334.45920
−1048.8220
16
−45.035
336.295
−1077.80
−45.035110
334.23630
−1048.8480
20
−45.049
336.284
−1077.67
−45.035060
334.52330
−1048.7200Figure 4
The motion of two negative solitary waves before and after the interaction at time = 0 and 20.


### 4.3. Maxwellian Initial Condition

The numerical solutions of the RLW equation with the Maxwellian initial condition(26)u(x.0)=exp(-(x-7)2)
and boundary conditions
(27)u(0,t)=u(50,t)=0
are carried out for the evolution of solitary waves for μ=0.04. We take the smaller space step h=0.01 and time step Δt=0.01. For μ=0.04 the Maxwellian develops into a single solitary wave plus a small developed oscillating tail as shown in Figure 5 at time t=12. The values of C1, C2, and C3 are given in Table 10; each value is satisfactorily constant.Table 10
Invariant values for Maxwellian initial condition,μ=0.04.
t
Present method
[9]
C
1
C
2
C
3
C
1
C
2
C
3
0
1.77245
1.30332
4.78327
1.772454
1.303446
4.783268
3
1.77245
1.30466
4.79137
1.772500
1.341671
4.918245
6
1.77245
1.30457
4.79124
1.790807
1.343856
4.922796
9
1.77245
1.30436
4.79039
1.775987
1.343681
4.917124
12
1.77245
1.30417
4.78957
1.767234
1.377576
4.913779The motion of Maxwellian initial condition initial solution and solution at time = 0 and 12.
(a)
(b)


## 4.1. Motion of the Solitary Waves

Case 1.
The equation represents a single soliton with amplitude3c and constant speed c+1. For the purpose of comparing with the earlier work, the computations are done for the parameters ε=1, μ=1, and x0=0. Consider c=0.1, as it seen from Table 1, the numerical results of the present method are in very good agreement with their analytical values obtained from the exact solution. Moreover, Table 2 displays a comparison of the approximate solution by the present method with those obtained using MOL developed by Griffiths and Schiesser in [12] at the time t=18. The error norms at t=18 obtained by the present method are smaller than those given in [12]. Also, Table 3 displays a comparison of error norms obtained by the present method with those obtained by using finite difference solution [7] and finite element [9].Table 1
L
2 and L∞ errors for (1) with c=0.1  Δx=0.1,Δt=0.1 and over the region -40≤x≤60.
t
L
2
L
∞
4
0.69452 × 10−4
0.264411 × 10−4
8
0.137889 × 10−3
0.534151 × 10−4
12
0.204502 × 10−3
0.792604 × 10−4
16
0.268849 × 10−3
0.103317 × 10−3
20
0.33080 × 10−3
0.125611 × 10−3Table 2
Absolute errors  for (1) with c=0.1  n=199 and over the region -50≤x≤50,t=18.
x
Present method
MOL [7]
−49.5
3.44444 × 10−9
0.000000000
−39.5
3.01068 × 10−8
0.000000000
−29.5
5.43620 × 10−9
0.000000000
−19.5
2.01194 × 10−7
1.00000 × 10−6
−9.5
3.91714 × 10−6
11.0000 × 10−5
0.5
4.39242 × 10−5
2.03000 × 10−4
10.5
2.31414 × 10−5
3.59000 × 10−3
20.5
1.43949 × 10−4
1.94700 × 10−3
30.5
1.9499 × 10−4
2.38500 × 10−3
40.4
4.67298 × 10−5
1.33000 × 10−4Table 3
L
2 and L∞ errors for (1) with c=0.1  Δx=0.1,Δt=0.1 and over the region -40≤x≤60.
t
Present method
Finite difference
Finite element
L
2
L
∞
L
2
L
∞
L
2
L
∞
4
0.6945 × 10−4
0.26441 × 10−4
0.1200 × 10−3
0.05000 × 10−3
0.1176 × 10−3
0.5622 × 10−3
8
0.13788 × 10−3
0.53415 × 10−4
0.2300 × 10−3
0.09000 × 10−3
0.23108 × 10−3
0.1088 × 10−3
12
0.20450 × 10−3
0.79260 × 10−4
0.3400 × 10−3
0.1400 × 10−3
0.3374 × 10−3
0.1544 × 10−3
16
0.26884 × 10−3
0.10331 × 10−3
0.45000 × 10−3
0.18000 × 10−3
0.43784 × 10−3
0.19377 × 10−3
20
0.33080 × 10−3
0.12561 × 10−3
0.55000 × 10−3
0.21000 × 10−3
0.53223 × 10−3
0.2272 × 10−3The error norms at each time are smaller than those given in [7, 9]. The invariant values C1, C2, and C3 are given in Table 4 and compared with the results obtained by [8]. From this table, it can be seen that the conservation properties (using the present method) are very good in comparison with the analytical values for the invariants given by C1=6c/k=3.9799497, C2=(12c2/k)+(48c2μ/5)=0.81046249, and C3=(36c2/k)+(144c3/5k)=2.579007 [8].Table 4
Invariant values for a single solitary wave: amplitude = 0.3,Δx=0.125,Δt=0.1 and over the region  -40≤x≤60.
t
Present method
Cubic spline collocation method [8]
C
1
C
2
C
3
C
1
C
2
C
3
0
3.97990
0.810461
2.57901
3.97992
0.810462
2.57900
4
3.97995
0.810461
2.57901
3.97995
0.810276
2.57839
8
3.97997
0.810461
2.57901
3.97997
0.810276
2.57839
12
3.97998
0.810461
2.57901
3.97998
0.810276
2.57839
16
3.97999
0.810461
2.57901
3.97998
0.810276
2.57839
20
3.97996
0.810461
2.57901
3.97996
0.810276
2.57839Forc=0.1, the relative changes in the invariants are very small (the quantities C1, C2, and C3 change by less than 1.03×10-5, 1.49×10-6, and 3×10-6, respectively, by t=20). Initial solution and solitary wave profile at t=20 are depicted in Figure 1.Figure 1
Graphs of initial solution and solution at time = 20 and amplitude = 0.3.Case 2.
In the second experiment, a smaller solitary wave of amplitude 0.09 (c=0.03) has been modeled, and the results of simulation are given in Table 5. As the amplitude of the solitary wave is reduced the pulse broadens, and it may be necessary to increase the solution range in order to maintain accuracy. The range here is doubled from -40≤x≤60 to -80≤x≤120. The error norms for c=0.03 are recorded in Table 6 for different values of time t. With the range -40≤x≤60, Δt=0.1, excellent results were obtained, L2 and L∞ error norms remain less than L2=0.637×10-3 and L∞=0.233×10-3, whereas they are L2=0.638×10-3 and L∞=0.233×10-3 for the finite difference method [7] and L2=0.572×10-3 and L∞=0.364×10-3 for the finite element method [9]. Moreover, it can be seen that the error norms (L2 and L∞) are somewhat small compared with those quoted by others [7, 9].Table 5
L
2 and L∞ errors for (1) with c=0.03  Δx=0.1,Δt=0.1 and over the region -80≤x≤120.
t
L
2
L
∞
4
0.51431 × 10−5
0.14805 × 10−5
8
0.10281 × 10−4
0.29961 × 10−5
12
0.15408 × 10−4
0.45367 × 10−5
16
0.20518 × 10−4
0.60908 × 10−5
20
0.25606 × 10−4
0.76472 × 10−5Table 6
L
2 and L∞ errors for (1) with c=0.03  Δx=0.1,Δt=0.1 and over the region -40≤x≤60.
t
Present method
Finite difference
Finite element
L
2
L
∞
L
2
L
∞
L
2
L
∞
4
0.14956 × 10−3
0.12290 × 10−3
0.1500 × 10−3
0.1230 × 10−3
0.3983 × 10−3
0.22450 × 10−3
8
0.32000 × 10−3
0.16600 × 10−3
0.3200 × 10−3
0.1660 × 10−3
0.50528 × 10−3
0.21706 × 10−3
12
0.466159 × 10−3
0.17900 × 10−3
0.4670 × 10−3
0.1790 × 10−3
0.5357 × 10−3
0.20919 × 10−3
16
0.56635 × 10−3
0.18454 × 10−3
0.5670 × 10−3
0.1850 × 10−3
0.54934 × 10−3
0.20165 × 10−3
20
0.63748 × 10−3
0.23300 × 10−3
0.6380 × 10−3
0.2330 × 10−3
0.57247 × 10−3
0.36498 × 10−3The invariantsC1, C2,  and C3 for RLW at different values of time t are also recorded in Table 7, where it can be seen that the invariants remain constant with respect to time. In comparison with the analytical values C1=2.109407, C2=0.127302, and  C3=0.388806, it can be seen from Table 7 that the quantities C1, C2, and C3 obtained using present method are very close to the corresponding exact values and compared with the results obtained by [8].Table 7
Invariant values for a single solitary wave: amplitude = 0.09,Δx=0.125,Δt=0.1 and over the region  -40≤x≤60.
t
Present method
Cubic spline collocation method [8]
C
1
C
2
C
3
C
1
C
2
C
3
0
2.10700
0.127302
0.38806
2.10702
0.127302
0.388804
4
2.10836
0.127302
0.38807
2.10838
0.127294
0.388781
8
2.10927
0.127302
0.38807
2.10928
0.127294
0.388782
12
2.10980
0.127302
0.38807
2.10981
0.127294
0.388782
16
2.10981
0.127302
0.38807
2.10981
0.127294
0.388782
20
2.10899
0.127302
0.38806
2.10900
0.127294
0.388781Forc=0.03, the relative changes in the quantities C1, C2, and C3 are changed by less than 4.1×10-4, 1.5×10-7, and 3.1×10-7, respectively, by t=20. Initial solution and solitary wave profile at t=20 are depicted in Figure 2.Figure 2
Graphs of initial solution and solution at time = 20 and amplitude = 0.09.


## 4.2. The Interaction of Solitary Wave

The interaction of two positive solitary waves is studied by using the initial condition given by the linear sum of two separate solitary waves of various amplitudes(25)u(x.0)=u1+u2,
where ui=3cisech2[ki(x-xi)], i=1,2, and ci=4ki2/(1-4ki2), with ε=μ=1.Case 1.
In this case, the numerical example consists in the interaction of two positive solitary waves defined by the initial condition (25), and the boundary conditions u(a,t)=u(b,t)=0.The RLW parameters used in the numerical simulation arek1=0.4, k2=0.3, x1=15, x2=35, ε=μ=1, and Δt=0.01. This wave moves across the space interval 0≤x≤120 and the time interval 0≤t≤25.Figure3 shows the interaction of two soliton solutions of RLW equation. From this figure, it can be seen that the faster pulse interacts with, and emerges ahead of, the slower pulse, with the shape and velocity of each soliton retained. Numerical check on the conservation of mass, momentum, and energy shows that the three quantities remain constant with respect to time t (see Table 8). The computed values of the invariants C1, C2, and C3 are satisfactorily constant compared with the corresponding invariant values obtained by [9]. The results of the present method at different times are shown in Figure 3, where it is possible to observe that the higher amplitude solitary wave passes through the smaller wave with no change in its waveform.Table 8
Invariant values for the interaction of two positive solitary waves over the region0≤x≤120.
t
Present method
[9]
C
1
C
2
C
3
C
1
C
2
C
3
0
37.9165
120.480
744.081
37.916520
120.52280
744.0815
4
37.9170
120.471
743.984
37.911700
121.16020
736.9443
8
37.9172
120.445
743.676
37.859750
119.73170
728.51730
12
37.9173
120.203
740.702
37.696670
119.63340
725.72360
16
37.9174
120.230
741.050
37.529160
119.41850
725.83940
20
37.9174
120.441
743.654
37.647300
119.80410
727.19480
25
37.9175
120.462
743.918
38.050100
119.83550
727.439200Figure 3
The motion of two positive solitary waves before and after the interaction at time = 0 and 25.Case 2.
The interaction of two negative solitary waves with initial condition (25) was studied over the region 0≤x≤120 with k1=0.6, k2=0.8, x1=82, x2=67, ε=μ=1, and Δt=0.2, Δx=0.5. The invariants C1, C2, and C3 recorded in Table 9; the computed values are satisfactorily constant compared with the corresponding invariant values obtained by [9]. The two negative solitary waves before and after the interaction are shown in Figure 4.Table 9
Invariant values for the interaction of two negative solitary waves over the region0≤x≤120.
t
Present method
[9]
C
1
C
2
C
3
C
1
C
2
C
3
0
−45.035
336.232
−1077.82
−45.034960
336.96020
−1077.824
4
−45.035
336.328
−1077.97
−45.035000
334.25310
−1048.6640
8
−45.035
336.014
−1068.16
−45.035080
332.98740
−1041.3070
12
−45.035
336.172
−1077.67
−45.035100
334.45920
−1048.8220
16
−45.035
336.295
−1077.80
−45.035110
334.23630
−1048.8480
20
−45.049
336.284
−1077.67
−45.035060
334.52330
−1048.7200Figure 4
The motion of two negative solitary waves before and after the interaction at time = 0 and 20.


## 4.3. Maxwellian Initial Condition

The numerical solutions of the RLW equation with the Maxwellian initial condition(26)u(x.0)=exp(-(x-7)2)
and boundary conditions
(27)u(0,t)=u(50,t)=0
are carried out for the evolution of solitary waves for μ=0.04. We take the smaller space step h=0.01 and time step Δt=0.01. For μ=0.04 the Maxwellian develops into a single solitary wave plus a small developed oscillating tail as shown in Figure 5 at time t=12. The values of C1, C2, and C3 are given in Table 10; each value is satisfactorily constant.Table 10
Invariant values for Maxwellian initial condition,μ=0.04.
t
Present method
[9]
C
1
C
2
C
3
C
1
C
2
C
3
0
1.77245
1.30332
4.78327
1.772454
1.303446
4.783268
3
1.77245
1.30466
4.79137
1.772500
1.341671
4.918245
6
1.77245
1.30457
4.79124
1.790807
1.343856
4.922796
9
1.77245
1.30436
4.79039
1.775987
1.343681
4.917124
12
1.77245
1.30417
4.78957
1.767234
1.377576
4.913779The motion of Maxwellian initial condition initial solution and solution at time = 0 and 12.
(a)
(b)


## 5. Conclusions

MOL is developed to obtain numerical solutions of the RLW equation. The efficiency of the proposed method was tested on problems from the literature, and its accuracy was examined by the error normsL2 and L∞. Numerical experiments for the RLW equation were reported for a single soliton, interaction of two solitary waves, and Maxwellian initial condition. The obtained results show that the error norms are reasonably small. The results also suggest that the present method’s application is easier than many other numerical techniques such as finite difference and finite element. An important advantage to be gained from the use of this method is to produce very accurate results. Using the Von Neumann stability analysis, it is shown that the proposed method is marginally stable. The error norms computed by the present algorithm with different amplitudes compared to the previous results were found to be smaller.In cases in whichL2 and L∞ cannot be evaluated, we verify the accuracy of the method by computing the conservation quantities C1, C2, and C3. The numerical test problems are in agreement with related literature [7, 9, 12]. So we deduce that this algorithm is more accurate, and it will also be useful for solving similar nonlinear partial differential equations.



---

*Source: 804317-2013-04-17.xml*