# A Modified Equation Approach to Selecting a Nonstandard Finite Difference Scheme Applied to the Regularized Long Wave Equation

Authors: E. Momoniat
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

Two nonstandard finite difference schemes are derived to solve the regularized long wave equation. The criteria for choosing the “best” nonstandard approximation to the nonlinear term in the regularized long wave equation come from considering the modified equation. The two “best” nonstandard numerical schemes are shown to preserve conserved quantities when compared to an implicit scheme in which the nonlinear term is approximated in the usual way. Comparisons to the single solitary wave solution show significantly better results, measured in theL
2 and L
∞ norms, when compared to results obtained using a Petrov-Galerkin finite element method and a splitted quadratic B-spline collocation method. The growth in the error when simulating the single solitary wave solution using the two “best” nonstandard numerical schemes is shown to be linear implying the nonstandard finite difference schemes are conservative. The formation of an undular bore for both steep and shallow initial profiles is captured without the formation of numerical instabilities.


---

## Body

## 1. Introduction

In this paper we derive two nonstandard finite difference schemes to solve the regularized long wave (RLW) equation. The RLW equation has been derived by Peregrine [1] to model the propagation of a unidirectional weakly nonlinear and weakly dispersive water wave known as an undular bore. The model equation is given by the mixed space-time partial differential equation:
(1)
∂
u
∂
t
-
∂
3
u
∂
t
∂
x
2
+
∂
u
∂
x
+
u
∂
u
∂
x
=
0
,
x
∈
ℝ
,
t
>
0
,
where
(2)
u
(
x
,
0
)
=
g
(
x
)
,
x
∈
ℝ
.
The RLW equation (1) satisfies the Dirichlet boundary conditions
(3)
u
⟶
0
,
x
±
∞
.
Benjamin et al. [2] have proposed the RLW equation (1) as an alternative to the Korteweg de Vries (KdV) equation [3]
(4)
∂
u
∂
t
+
∂
u
∂
x
+
u
∂
u
∂
x
+
∂
3
u
∂
x
3
=
0
,
to model phenomena in wave propagation and other physical systems. The criteria for selecting the “best” nonstandard approximation to the nonlinear term in the RLW equation come from considering the effects of the modified equation on the single solitary wave solution admitted by the RLW equation (1). The modified equation is obtained by considering coefficients 𝒪
(
Δ
t
/
Δ
x
2
) of a Taylor series expansion of an implicit scheme and nonstandard finite difference schemes approximating the RLW equation. We show that the two “best” nonstandard approximations to the nonlinear term make the nonstandard finite difference schemes approximating the RLW equation (1) second-order accurate, that is, 𝒪
(
Δ
t
2
/
Δ
x
4
), where Δ
t is the time step length and Δ
x is the spatial step length. To the best of our knowledge this is the first paper in which the analysis of the modified equation has been used to justify the choice of nonstandard approximation to a nonlinear term.The nonstandard finite difference schemes are reduced to linear tridiagonal systems that are solved using a Thomas algorithm [4]. We show that the nonstandard numerical schemes preserve the conserved quantities (these quantities have been obtained by Olver [5]):
(5)
I
1
=
∫
-
∞
∞
‍
u
d
x
,
(6)
I
2
=
∫
-
∞
∞
‍
(
u
2
+
(
u
x
)
2
)
d
x
,
(7)
I
3
=
∫
-
∞
∞
‍
(
u
3
+
3
u
2
)
d
x
.
For completeness, we show that a standard implicit scheme in which the nonlinear term is approximated in the usual way does not preserve (5)–(7). Numerical solutions obtained from simulations of the nonstandard finite difference schemes modelling the single solitary wave solution [1] of the RLW equation (1) given by
(8)
u
(
x
,
t
)
=
3
v
sec
h
2
(
k
0
x
+
k
1
t
+
x
0
)
,
where
(9)
k
0
=
1
2
v
1
+
v
,
k
1
=
-
1
2
v
(
1
+
v
)
,
and x
0 is an arbitrary constant, show significantly better results when measured in the L
2 and L
∞ norms when compared to simulations of the single solitary wave solution using a Petrov-Galerkin finite element method [6] and a splitted quadratic B-spline collocation method [7]. We also show that a nonstandard finite difference scheme is able to capture the formation of an undular bore for both steep and shallow initial profiles without the formation of numerical instabilities on the solution profile.Various numerical approaches have been developed to solve the RLW equation (1). Peregrine [1] implemented a finite difference scheme to solve (1) that is first-order accurate in time and second-order accurate in space. Eilbeck and McGuire [8, 9] have derived a three-level second-order accurate space and time scheme. An implicit finite difference scheme has been implemented by Kutluay and Esen [10]. Jain et al. [11] have developed a numerical scheme that is second-order accurate in both space and time based on cubic spline interpolations. Bona et al. [12] review fourth-order finite difference schemes to solve (1). Numerical methods based on spline Galerkin techniques have also been implemented [7, 13–19]. An analysis of finite element techniques applied to the RLW equation has been undertaken by Luo and Liu [20]. Pseudospectral methods have been developed by Ben-yu and Manoranjan [21] and Sloan [22] with extensions to Chebyshev collocation methods investigated by Ali [23]. Avilez-Valente and Seabra-Santos [6] show that their approach using a Petrov-Galerkin finite element method outperforms numerical methods based on spline Galerkin as well as pseudospectral methods but is equivalent to the splitted quadratic B-spline collocation method introduced by Daǧ et al. [7]. Avilez-Valente and Seabra-Santos [6] use finite elements in both time and space. They introduce additional streamline upwind terms in the weight functions to ensure that dispersion correction and selective dissipation occur. Two versions of the numerical scheme are derived. A nonlinear version of the resulting predictor-corrector scheme is found to be third-order accurate. A linear version of the resulting scheme is found to be fourth-order accurate. Both versions of the numerical scheme are implicit and conditionally stable. Raslan [24] has employed a collocation method based on a cubic B-spline method to solve the RLW equation. Bhardwaj and Shankar [25] developed a numerical technique based on quintic splines. Cubic [26], quartic [27], and quintic [28, 29] B-spline collocation methods applied to a splitted RLW equation have been investigated to solve (1). Irk [30] extends the work of Bhardwaj and Shankar [25] by considering an Adams-Moulton time-integration scheme coupled with a quintic spline collocation method for the spatial variable. Araújo and Durán [31] have discussed the importance of using numerical schemes that conserve the invariants of (1). More specifically, they show that the error growth in schemes that are conservative is linear whilst nonconservative schemes have quadratic growth in the errors. Araújo and Durán [31] base their results on the fact that the RLW equation admits a Hamiltonian structure. Also, one of its invariants is generated by a symmetry group. Solitary wave solutions admitted by the RLW equation then occur as critical points of the Hamiltonian function [31]. This geometric interpretation of the solitary wave solution is used to determine the relationship between numerical schemes that conserve quantities (5)–(7) and the linear growth in error. Cai derives the Bridge's multisymplectic form of the RLW equation to derive a ten-point multisymplectic scheme to solve the RLW equation [32].Avilez-Valente and Seabra-Santos [6] perform a thorough comparison of their results with the work of Daǧ [14], Daǧ et al. [7], Gardner et al. [17], Gardner et al. [18], and Jain et al. [11]. Avilez-Valente and Seabra-Santos [6] conclude that their Petrov-Galerkin finite element formulation and the splitted quadratic B-spline collocation method implemented by Daǧ et al. [7] are the most effective approaches for solving the RLW equation. We compare our results to the results obtained by Avilez-Valente and Seabra-Santos [6] and Daǧ et al. [7]. We show the results obtained by Avilez-Valente and Seabra-Santos [6] and Daǧ et al. [7] in Tables 1 and 2, respectively. In Tables 1 and 2 the L
2 and L
∞ norms are given by
(10)
L
2
=
|
u
i
j
+
1
-
u
*
(
x
i
,
t
j
)
|
2
=
(
Δ
x
∑
i
=
0
n
‍
|
u
i
j
+
1
-
u
*
(
x
i
,
t
j
)
|
2
)
1
/
2
,
L
∞
=
|
u
i
j
+
1
-
u
*
(
x
i
,
t
j
)
|
∞
=
max
⁡
i
|
u
i
j
+
1
-
u
*
(
x
i
,
t
j
)
|
,
where u
*
(
x
i
,
t
j
) is the value of the single solitary wave solution (8) to the RLW equation at x
=
x
i and t
=
t
j for v
=
0.1 and x
0
=
0.Table 1
Table of values obtained by Avilez-Valente and Seabra-Santos [6] using a Petrov-Galerkin finite element method.
Time
I
1
I
2
I
3
L
2
×
1
0
3
L
∞
×
1
0
3
0
3.979927
0.8104625
2.579007
0.00000
0.00000
4
3.979944
0.8104625
2.579007
0.01549
0.00695
8
3.979947
0.8104625
2.579007
0.03000
0.01343
12
3.979946
0.8104625
2.579007
0.04306
0.01868
16
3.979939
0.8104625
2.579007
0.05464
0.02292
20
3.979909
0.8104625
2.579007
0.06493
0.02643Table 2
Table of values obtained by Daǧ et al. [7] using a splitted quadratic B-spline collocation method.
Time
I
1
I
2
I
3
L
2
×
1
0
3
L
∞
×
1
0
3
0
3.979926
0.8104646
2.579007
0.00000
0.00000
4
3.979930
0.8104626
2.578998
0.01757
0.00693
8
3.979924
0.8104606
2.578994
0.02249
0.00887
12
3.979926
0.8104621
2.578997
0.03355
0.01072
16
3.979918
0.8104624
2.578999
0.04075
0.01224
20
3.979890
0.8104625
2.578999
0.04315
0.01321The work undertaken in this paper is a major contribution to the field because of the accuracy of the results obtained for the low computational cost when compared to other methods. An example of the high computational cost that can be incurred to obtain accurate results can be seen in the work of Bona et al. [12]. Bona et al. [12] have implemented a Runge-Kutta approximation as a first step in a leap frog approximation to the time derivative. Bona et al. [12] approximate the spatial derivatives using a Gregory-Euler-Maclaurin (GEM) scheme which has a truncation error 𝒪
(
Δ
x
4
). Bona et al. [12] showed that by using a semidiscrete scheme they could obtain a tridiagonal system. Bona et al. [12] compare leap frog (LF), fourth-order Runge-Kutta (RK), and multistep (MS) predictor-corrector approximations to the time derivative. Fourth-order discretizations in space that lead to Störmer-Numerov (SN) type schemes [33] are also considered. Fourth-order Runge-Kutta approximations to the time derivative coupled with the SN approximation to the spatial derivative (SN-RK schemes) were simulated. The SN-RK scheme performed worse in the simulations undertaken by Bona et al. [12]. The Petrov-Galerkin finite element method of Avilez-Valente and Seabra-Santos [6] is also reduced to a tridiagonal system but has to be evaluated several times at each time step because of the high-order time integration method involved. The splitted quadratic B-spline collocation scheme implemented by Daǧ et al. [7] leads to a pentadiagonal system that also has to be evaluated multiple times at each time step because of the high-order time integration method that is chosen. As mentioned earlier, the nonstandard finite difference schemes are reduced to linear tridiagonal systems that we solve using a Thomas algorithm [4]. The implicit nature of the schemes means that at each time step the system matrix is recalculated and the system resolved. This is computationally not as expensive as the methods mentioned above.In this paper we investigate nonstandard approximations to the nonlinear term in the RLW equation to derive nonstandard finite difference schemes approximating (1). An introduction to nonstandard finite difference schemes and their applications to differential equations is given in the paper by Mickens [34]. Applications to ordinary differential equations are discussed in Mickens [35] and Anguelov and Lubuma [36]. Theoretical developments of nonstandard finite difference schemes and their applications to partial differential equations are investigated by Anguelov and Lubuma [37]. Mickens [38] considers the general application of nonstandard finite difference schemes to nonlinear reaction-advection equations. Mickens [39] extends these results by considering reaction-diffusion equations and derives positivity preserving nonstandard finite difference schemes. Jordan [40] derives a nonstandard finite difference scheme satisfying a positivity condition to model heat transfer with a quartic nonlinearity. Mickens [41] constructs nonstandard finite difference schemes to solve a Fisher partial differential equation that has application in ecology. Positivity preserving schemes for Maxwell-Cattaneo wave equations are derived by Mickens and Jordan [42].In this paper we select the “best” nonstandard approximation to the nonlinear term in the RLW equation by investigating the effects of the modified equation on the single solitary wave solution (8). Warming and Hyett [43] investigate the stability of simple linear partial differential equations by investigating the modified partial differential equation. In a series of papers Villatoro and Ramos [44–47] used the modified equation approach to analyze and improve the efficiency of an Euler forward difference method to solve linear and nonlinear differential equations. Junk and Yang [48] combine a modified equation and truncation error analysis in an asymptotic analysis of finite difference schemes. Rus and Villatoro [49] derive high-order finite difference schemes that preserve conserved quantities for Rosenau-Hyman equation based on the corresponding modified equations.The paper is divided up as follows: in Section2 we derive and analyze the nonstandard finite difference schemes. The numerical schemes are simulated in Section 3. Concluding remarks are made in Section 4.


## 2. Analysis of Modified Equation

In this section we analyze the modified equation that we obtain from a finite difference approximation of the RLW equation (1). We define u
i
j
=
u
(
x
i
,
t
j
), where the spatial domain is discretized into n
+
1 points and x
i
=
i
Δ
x, where Δ
x is the spatial step length. The time t
j is defined by t
j
=
j
Δ
t, where Δ
t is the time step length. Central difference approximations to the spatial derivatives are given by
(11)
∂
u
∂
x
≈
u
i
+
1
j
-
u
i
-
1
j
2
Δ
x
+
𝒪
(
Δ
x
2
)
,
(12)
∂
2
u
∂
x
2
≈
u
i
+
1
j
-
2
u
i
j
+
u
i
-
1
j
Δ
x
2
+
𝒪
(
Δ
x
2
)
.
A forward difference approximation to the time derivative is given by
(13)
∂
u
∂
t
≈
u
i
j
+
1
-
u
i
j
Δ
t
+
𝒪
(
Δ
t
)
.
We write the RLW equation (1) as
(14)
∂
∂
t
(
u
-
∂
2
u
∂
x
2
)
+
∂
u
∂
x
+
1
2
∂
u
2
∂
x
=
0
.
Substituting (11), (12), and (13) with (14) we obtain an implicit numerical scheme given by
(15)
-
1
Δ
t
(
u
i
+
1
j
+
1
-
u
i
+
1
j
)
+
(
Δ
x
2
+
2
Δ
t
)
(
u
i
j
+
1
-
u
i
j
)
-
1
Δ
t
(
u
i
-
1
j
+
1
-
u
i
-
1
j
)
+
Δ
x
2
(
u
i
+
1
j
-
u
i
-
1
j
)
+
Δ
x
4
(
(
u
i
+
1
j
)
2
-
(
u
i
-
1
j
)
2
)
=
0
,
where we have approximated the partial derivative ∂
u
/
∂
x explicitly. To obtain the modified equation we substitute the Taylor series approximations
(16)
u
i
+
1
j
=
u
(
x
i
+
1
,
t
j
)
=
u
(
x
i
,
t
j
)
+
Δ
x
∂
u
∂
x
(
x
i
,
t
j
)
+
Δ
x
2
2
!
∂
2
u
∂
x
2
(
x
i
,
t
j
)
+
O
(
Δ
x
3
)
,
u
i
-
1
j
=
u
(
x
i
-
1
,
t
j
)
=
u
(
x
i
,
t
j
)
-
Δ
x
∂
u
∂
x
(
x
i
,
t
j
)
+
Δ
x
2
2
!
∂
2
u
∂
x
2
(
x
i
,
t
j
)
-
O
(
Δ
x
3
)
,
and similar approximations for u
i
+
1
j
+
1, u
i
j
+
1, and u
i
-
1
j
+
1 with (15) to obtain the partial differential equation
(17)
∂
∂
t
(
u
-
∂
2
u
∂
x
2
)
+
∂
u
∂
x
+
1
2
∂
u
2
∂
x
+
Δ
t
2
Δ
x
2
(
∂
2
u
∂
t
2
-
∂
4
u
∂
t
2
∂
x
2
)
=
0
.
We let the modified equation be defined by Δ
ϵ. Multiplying the modified equation Δ
ϵ by Δ
t
/
(
2
Δ
x
2
) we obtain
(18)
Δ
t
Δ
ϵ
2
Δ
x
2
=
∂
2
u
∂
t
2
-
∂
4
u
∂
t
2
∂
x
2
.
To determine the effects of the modified equation (18) on the single solitary wave solution of the RLW equation (1) we substitute the single solitary wave solution (8) with (18) and plot Δ
t
Δ
ϵ
/
(
2
Δ
x
2
) at different times in Figure 1.Figure 1
Plot comparing the modified equations obtained for explicit (18) (
-
-
-
) and implicit (20) (solid line) approximations to the derivative ∂
u
/
∂
x after substituting the single solitary wave solution to the RLW equation (8) with (18) and (20) taking v
=
0.1 and x
0
=
0.We improve the stability and accuracy of the numerical scheme (15) by approximating the derivative ∂
u
/
∂
x implicitly to obtain an implicit scheme to solve the RLW equation (14) given by
(19)
-
1
Δ
t
(
u
i
+
1
j
+
1
-
u
i
+
1
j
)
+
(
Δ
x
2
+
2
Δ
t
)
(
u
i
j
+
1
-
u
i
j
)
-
1
Δ
t
(
u
i
-
1
j
+
1
-
u
i
-
1
j
)
+
Δ
x
4
(
u
i
+
1
j
+
1
+
u
i
+
1
j
-
u
i
-
1
j
+
1
-
u
i
-
1
j
)
+
Δ
x
4
(
(
u
i
+
1
j
)
2
-
(
u
i
-
1
j
)
2
)
=
0
.
The corresponding modified equation is given by
(20)
Δ
t
Δ
ϵ
2
Δ
x
2
=
∂
2
u
∂
t
2
+
∂
2
u
∂
t
∂
x
-
∂
4
u
∂
t
2
∂
x
2
.
Once again we determine the effects of the modified equation (20) on the single solitary wave solution of the RLW equation (1) by substituting the single solitary wave solution (8) with (20) and plotting the resulting function in Figure 1.From the results indicated in Figure1 we note that the impact of the modified equation on the single solitary wave solution of the RLW equation (1) is smaller for the modified equation (20) derived from the finite difference scheme (19) where the derivative ∂
u
/
∂
x is approximated implicitly than for (18) derived from the finite difference scheme (15) where the derivative ∂
u
/
∂
x is approximated explicitly. From this result we can conclude that the finite difference scheme (19) will yield better results than (15) when used to determine solutions to the RLW equation (1). We note that the curves in Figure 1 propagate with the wave velocity of the single solitary wave solution (8) maintaining a constant wave amplitude.We firstly consider the nonstandard approximation [41]
(21)
(
u
i
j
)
2
⟶
2
(
u
i
j
)
2
-
u
i
j
+
1
u
i
j
.
The implicit scheme (19) reduces to
(22)
-
1
Δ
t
(
u
i
+
1
j
+
1
-
u
i
+
1
j
)
+
(
Δ
x
2
+
2
Δ
t
)
(
u
i
j
+
1
-
u
i
j
)
-
1
Δ
t
(
u
i
-
1
j
+
1
-
u
i
-
1
j
)
+
Δ
x
4
(
u
i
+
1
j
+
1
+
u
i
+
1
j
-
u
i
-
1
j
+
1
-
u
i
-
1
j
)
+
Δ
x
4
[
(
2
(
u
i
+
1
j
)
2
-
u
i
+
1
j
+
1
u
i
+
1
j
)
-
(
2
(
u
i
-
1
j
)
2
-
u
i
-
1
j
+
1
u
i
-
1
j
)
]
=
0
,
where the corresponding modified equation is given by
(23)
Δ
t
Δ
ϵ
2
Δ
x
2
=
∂
2
u
∂
t
2
-
∂
u
∂
t
∂
u
∂
x
+
(
1
-
u
)
∂
2
u
∂
t
∂
x
-
∂
4
u
∂
t
2
∂
x
2
.We next consider the nonstandard approximation [36]
(24)
(
u
i
-
1
j
)
2
⟶
u
¯
i
-
1
j
u
i
-
1
j
+
1
,
(
u
i
+
1
j
)
2
⟶
u
¯
i
+
1
j
u
i
+
1
j
+
1
,
where
(25)
u
¯
i
-
1
j
=
1
2
(
u
i
-
1
j
+
u
i
j
)
,
u
¯
i
+
1
j
=
1
2
(
u
i
j
+
u
i
+
1
j
)
.
The implicit scheme (19) reduces to
(26)
-
1
Δ
t
(
u
i
+
1
j
+
1
-
u
i
+
1
j
)
+
(
Δ
x
2
+
2
Δ
t
)
(
u
i
j
+
1
-
u
i
j
)
-
1
Δ
t
(
u
i
-
1
j
+
1
-
u
i
-
1
j
)
+
Δ
x
4
(
u
i
+
1
j
+
1
+
u
i
+
1
j
-
u
i
-
1
j
+
1
-
u
i
-
1
j
)
+
Δ
x
8
[
(
u
i
j
+
u
i
+
1
j
)
u
i
+
1
j
+
1
-
(
u
i
-
1
j
+
u
i
j
)
u
i
-
1
j
+
1
]
=
0
,
where the corresponding modified equation is given by
(27)
Δ
t
Δ
ϵ
2
Δ
x
2
=
∂
2
u
∂
t
2
+
1
2
∂
u
∂
t
∂
u
∂
x
+
(
1
+
u
)
∂
2
u
∂
t
∂
x
-
∂
4
u
∂
t
2
∂
x
2
.
We substitute the single solitary wave solution (8) with (23) and (27) and plot the resulting functions in Figure 2 to determine the effects of the modified equations (23) and (27) on the single solitary wave solution of the RLW equation (1).Figure 2
Plot comparing the behaviour of the modified equations obtained for the nonstandard approximation (21) (
-
-
-
) and the nonstandard approximation (24) (solid line) to the nonlinear term in the RLW equation after substituting the single solitary wave solution to the RLW equation (8) with (23) and (27) taking v
=
0.1 and x
0
=
0.From the results indicated in Figure2 we note that the impact of the modified equation (27) on the single solitary wave solution of the RLW equation (1) is smaller than that of the modified equation (23). We can therefore conclude that the nonstandard approximation (24) is a better approximation to use for the nonlinear term in the RLW equation than (21). We can also conclude that the finite difference scheme (26) will produce better results than (22) when used to determine numerical solutions to the RLW equation. In both cases the functions obtained after substituting the single solitary wave solution of the RLW equation (1) with the modified equations (23) and (27) propagate with the wave velocity of the single solitary wave solution.We next consider the approximation [41]
(28)
(
u
i
j
)
2
⟶
u
i
j
+
1
u
i
j
.
Substituting (28) with (19) we obtain the nonstandard numerical scheme
(29)
-
1
Δ
t
(
u
i
+
1
j
+
1
-
u
i
+
1
j
)
+
(
Δ
x
2
+
2
Δ
t
)
(
u
i
j
+
1
-
u
i
j
)
-
1
Δ
t
(
u
i
-
1
j
+
1
-
u
i
-
1
j
)
+
Δ
x
4
(
u
i
+
1
j
+
1
+
u
i
+
1
j
-
u
i
-
1
j
+
1
-
u
i
-
1
j
)
+
Δ
x
4
[
u
i
+
1
j
+
1
u
i
+
1
j
-
u
i
-
1
j
+
1
u
i
-
1
j
]
=
0
.
The modified equation corresponding to (29) is given by
(30)
Δ
t
Δ
ϵ
Δ
x
2
=
∂
∂
t
(
∂
u
∂
t
-
∂
3
u
∂
t
∂
x
2
+
∂
u
∂
x
+
u
∂
u
∂
x
)
=
0
.
We consider the next order in Taylor series expansions (16) to obtain the modified equation
(31)
Δ
t
2
Δ
ϵ
6
Δ
x
4
=
∂
3
u
∂
t
3
+
3
2
∂
2
u
∂
t
2
∂
u
∂
x
+
3
2
∂
3
u
∂
t
2
∂
x
(
1
+
u
)
-
∂
5
u
∂
t
3
∂
x
2
.We finally consider the approximation [39]
(32)
(
u
i
j
)
2
⟶
1
3
(
u
i
-
1
j
+
u
i
j
+
u
i
+
1
j
)
u
i
j
+
1
.
The implicit scheme (19) reduces to
(33)
-
1
Δ
t
(
u
i
+
1
j
+
1
-
u
i
+
1
j
)
+
(
Δ
x
2
+
2
Δ
t
)
(
u
i
j
+
1
-
u
i
j
)
-
1
Δ
t
(
u
i
-
1
j
+
1
-
u
i
-
1
j
)
+
Δ
x
4
(
u
i
+
1
j
+
1
+
u
i
+
1
j
-
u
i
-
1
j
+
1
-
u
i
-
1
j
)
+
Δ
x
12
[
(
u
i
j
+
u
i
+
1
j
+
u
i
+
2
j
)
u
i
+
1
j
+
1
-
(
u
i
-
2
j
+
u
i
-
1
j
+
u
i
j
)
u
i
-
1
j
+
1
]
=
0
.
The corresponding modified equation is given by (30) and hence (31), respectively. We substitute the single solitary wave solution (8) with (31) and plot the resulting function in Figure 3.Figure 3
Plot showing the behaviour of the second-order modified equation (31) after substituting the single solitary wave solution to the RLW equation (8) with (31) taking v
=
0.1 and x
0
=
0.From the results plotted in Figure3 we observe that the second-order modified equation has the same order error as the first-order modified equations plotted in Figure 2. The solutions plotted in Figure 3 propagate with the wave velocity of the solitary wave solution.Since the first-order modified equation (30) is identically satisfied we can conclude that the nonstandard finite difference schemes (29) and (33) are second-order accurate. Hence we can conclude that the nonstandard approximations (28) and (32) are the “best” approximations to use for the RLW equation. In addition, we obtain the CFL condition Δ
t
/
Δ
x
2
≪
1 for the stability of the implicit finite difference scheme (19) and Δ
t
2
/
Δ
x
4
≪
1 for the stability of the nonstandard finite difference schemes (29) and (33). In the next section we simulate the nonstandard finite difference schemes (29) and (33) to determine which scheme conserves (5)–(7) the “best.”


## 3. Simulation of Numerical Schemes

For the first part of this section we take as the initial condition from the single solitary wave solution (8)
(34)
u
(
x
,
0
)
=
u
0
=
3
v
sec
h
2
(
k
0
x
)
,
where u
0
=
[
u
0
0
,
u
1
0
,
…
,
u
n
0
] and k
0
=
(
1
/
2
)
v
/
(
1
+
v
). We choose v
=
0.1 and x
0
=
0 and divide the interval x
=
[
-
50,70
] into n
+
1 equidistant points taking Δ
t
=
0.0001 and n
=
5000. We simulate the proposed nonstandard numerical schemes using MATHEMATICA on a 2.66 GHz quad core Windows machine with 8.00 GB RAM.We firstly consider the implicit scheme (19) that we write as
(35)
(
-
1
Δ
t
+
Δ
x
4
)
u
i
+
1
j
+
1
+
(
Δ
x
2
+
2
Δ
t
)
u
i
j
+
1
+
(
-
1
Δ
t
-
Δ
x
4
)
u
i
-
1
j
+
1
=
(
-
1
Δ
t
-
Δ
x
4
)
u
i
+
1
j
+
(
Δ
x
2
+
2
Δ
t
)
u
i
j
+
(
-
1
Δ
t
+
Δ
x
4
)
u
i
-
1
j
-
Δ
x
4
[
(
u
i
+
1
j
)
2
-
(
u
i
-
1
j
)
2
]
.
We then write the numerical scheme (35) as the linear system
(36)
A
u
j
+
1
=
d
,
where
(37)
A
=
[
1
0
0
0
⋯
0
0
0
0
a
1
b
1
c
1
0
⋯
⋮
⋮
⋮
⋮
0
a
2
b
2
c
2
⋯
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋯
a
n
-
2
b
n
-
2
c
n
-
2
0
⋮
⋮
⋮
⋮
⋯
0
a
n
-
1
b
n
-
1
c
n
-
1
0
0
0
0
⋯
0
0
0
1
]
,
(38)
d
=
[
0
(
-
1
Δ
t
-
Δ
x
4
)
u
2
j
+
(
Δ
x
2
+
2
Δ
t
)
u
1
j
+
(
-
1
Δ
t
+
Δ
x
4
)
u
0
j
-
Δ
x
4
[
(
u
2
j
)
2
-
(
u
0
j
)
2
]
(
-
1
Δ
t
-
Δ
x
4
)
u
3
j
+
(
Δ
x
2
+
2
Δ
t
)
u
2
j
+
(
-
1
Δ
t
+
Δ
x
4
)
u
1
j
-
Δ
x
4
[
(
u
3
j
)
2
-
(
u
1
j
)
2
]
⋮
(
-
1
Δ
t
-
Δ
x
4
)
u
n
-
1
j
+
(
Δ
x
2
+
2
Δ
t
)
u
n
-
2
j
+
(
-
1
Δ
t
+
Δ
x
4
)
u
n
-
3
j
-
Δ
x
4
[
(
u
n
-
1
j
)
2
-
(
u
n
-
3
j
)
2
]
(
-
1
Δ
t
-
Δ
x
4
)
u
n
j
+
(
Δ
x
2
+
2
Δ
t
)
u
n
-
1
j
+
(
-
1
Δ
t
+
Δ
x
4
)
u
n
-
2
j
-
Δ
x
4
[
(
u
n
j
)
2
-
(
u
n
-
2
j
)
2
]
0
]
.For the scheme (35)
(39)
a
i
=
-
1
Δ
t
-
Δ
x
4
,
b
i
=
Δ
x
2
+
2
Δ
t
,
c
i
=
-
1
Δ
t
+
Δ
x
4
.
The boundary conditions (3) determine the first and last rows of A and d. We use a Thomas algorithm [4] to solve the linear system of equations (36). We compare the single solitary wave solution (8) to the numerical solution obtained from the linear system (36) corresponding to the implicit scheme (35) using the L
2 and L
∞ norms as defined in (10). We also evaluate the conserved quantities (5)–(7) using a trapezoidal quadrature rule. We tabulate our results in Table 3. We find that the implicit scheme (35) does not preserve I
2 or I
3. At t
=
20 the L
2 and L
∞ norms are significantly better than the L
2 and L
∞ norms obtained by Avilez-Valente and Seabra-Santos [6] and Daǧ et al. [7] indicated in Tables 1 and 2, respectively.Table 3
Table comparing values obtained from the implicit scheme (35) to the single solitary wave solution (8) for x
0
=
0 and v
=
0.1, for x
∈
[
-
50,70
], and taking n
=
5000 and Δ
t
=
0.0001.
Time
I
1
I
2
I
3
L
2
×
1
0
3
L
∞
×
1
0
3
0
3.97995
0.810462
2.57901
0.00000
0.00000
4
3.97995
0.810463
2.57901
0.0044822
0.00190039
8
3.97995
0.810464
2.57901
0.00884828
0.00377625
12
3.97995
0.810466
2.57902
0.0130448
0.00550027
16
3.97995
0.810467
2.57902
0.0170392
0.00705772
20
3.97995
0.810468
2.57903
0.0208237
0.00846861We next consider the nonstandard numerical scheme (29). The nonstandard numerical scheme (29) can be written as
(40)
[
-
1
Δ
t
+
Δ
x
4
(
1
+
u
i
+
1
j
)
]
u
i
+
1
j
+
1
+
(
Δ
x
2
+
2
Δ
t
)
u
i
j
+
1
+
[
-
1
Δ
t
-
Δ
x
4
(
1
+
u
i
-
1
j
)
]
u
i
-
1
j
+
1
=
(
-
1
Δ
t
-
Δ
x
4
)
u
i
+
1
j
+
(
Δ
x
2
+
2
Δ
t
)
u
i
j
+
(
-
1
Δ
t
+
Δ
x
4
)
u
i
-
1
j
.
We write the numerical scheme (40) as the system (36) where A is given by (37) and where(41)
d
=
[
0
(
-
1
Δ
t
-
Δ
x
4
)
u
2
j
+
(
Δ
x
2
+
2
Δ
t
)
u
1
j
+
(
-
1
Δ
t
+
Δ
x
4
)
u
0
j
(
-
1
Δ
t
-
Δ
x
4
)
u
3
j
+
(
Δ
x
2
+
2
Δ
t
)
u
2
j
+
(
-
1
Δ
t
+
Δ
x
4
)
u
1
j
⋮
(
-
1
Δ
t
-
Δ
x
4
)
u
n
-
1
j
+
(
Δ
x
2
+
2
Δ
t
)
u
n
-
2
j
+
(
-
1
Δ
t
+
Δ
x
4
)
u
n
-
3
j
(
-
1
Δ
t
-
Δ
x
4
)
u
n
j
+
(
Δ
x
2
+
2
Δ
t
)
u
n
-
1
j
+
(
-
1
Δ
t
+
Δ
x
4
)
u
n
-
2
j
0
]
,such that
(42)
a
i
=
-
1
Δ
t
-
Δ
x
4
(
1
+
u
i
-
1
j
)
,
b
i
=
Δ
x
2
+
2
Δ
t
,
c
i
=
-
1
Δ
t
+
Δ
x
4
(
1
+
u
i
+
1
j
)
.
We tabulate our results in Table 4. We find that the nonstandard numerical scheme (40) shows no deviation in the evaluation of the conserved quantities (5)–(7). Also, the L
2 and L
∞ norms are once again significantly better than the L
2 and L
∞ norms obtained by Avilez-Valente and Seabra-Santos [6] and Daǧ et al. [7] indicated in Tables 1 and 2, respectively.Table 4
Table comparing values obtained from the nonstandard numerical scheme (40) to the single solitary wave solution (8) for x
0
=
0 and v
=
0.1, for x
∈
[
-
50,70
], and taking n
=
5000 and Δ
t
=
0.0001.
Time
I
1
I
2
I
3
L
2
×
1
0
3
L
∞
×
1
0
3
0
3.97995
0.810462
2.57901
0.00000
0.00000
4
3.97995
0.810462
2.57901
0.00424465
0.00167944
8
3.97995
0.810462
2.57901
0.00839433
0.00337615
12
3.97995
0.810462
2.57901
0.012402
0.00495651
16
3.97995
0.810462
2.57901
0.0162382
0.00639525
20
3.97995
0.810462
2.57901
0.0198962
0.00770612We finally consider the nonstandard numerical scheme (33). The nonstandard numerical scheme (33) can be written as
(43)
[
-
1
Δ
t
+
Δ
x
4
+
Δ
x
12
(
u
i
j
+
u
i
+
1
j
+
u
i
+
2
j
)
]
u
i
+
1
j
+
1
+
(
Δ
x
2
+
2
Δ
t
)
u
i
j
+
1
+
[
-
1
Δ
t
-
Δ
x
4
-
Δ
x
12
(
u
i
-
2
j
+
u
i
-
1
j
+
u
i
j
)
]
u
i
-
1
j
+
1
=
(
-
1
Δ
t
-
Δ
x
4
)
u
i
+
1
j
+
(
Δ
x
2
+
2
Δ
t
)
u
i
j
+
(
-
1
Δ
t
+
Δ
x
4
)
u
i
-
1
j
.
We write (43) as the linear system (36) where A and d are given by (37) and (41) such that
(44)
a
1
=
-
1
Δ
t
-
Δ
x
4
-
Δ
x
12
(
2
u
1
j
+
u
0
j
)
,
a
i
=
-
1
Δ
t
-
Δ
x
4
-
Δ
x
12
(
u
i
-
2
j
+
u
i
-
1
j
+
u
i
j
)
,
i
=
2
,
…
,
n
-
1
,
c
n
-
1
=
-
1
Δ
t
+
Δ
x
4
+
Δ
x
12
(
2
u
n
-
1
j
+
u
n
j
)
,
c
i
=
-
1
Δ
t
+
Δ
x
4
+
Δ
x
12
(
u
i
j
+
u
i
+
1
j
+
u
i
+
2
j
)
,
i
=
1
,
…
,
n
-
2
,
b
i
=
Δ
x
2
+
2
Δ
t
,
i
=
1
,
…
,
n
-
1
.
The unknown values of u in a
1 and c
n
-
1 are determined from Neumann boundary conditions obtained from the central difference approximation (11)
(45)
u
i
+
1
j
-
u
i
-
1
j
2
Δ
x
=
0
.
Evaluating (45) at i
=
0 and i
=
n we get that u
-
1
j
=
u
1
j and u
n
+
1
j
=
u
n
-
1
j. We tabulate the results obtained from the numerical scheme (43) in Table 5. We find that the nonstandard numerical scheme (43) conserves I
1, I
2, and I
3 obtained from (5)–(7). Once again the L
2 and L
∞ norms at t
=
20 are significantly better than the L
2 and L
∞ norms obtained by Avilez-Valente and Seabra-Santos [6] and Daǧ et al. [7] indicated in Tables 1 and 2, respectively.Table 5
Table comparing values obtained from the nonstandard numerical scheme (43) to the single solitary wave solution (8) for x
0
=
0 and v
=
0.1, for x
∈
[
-
50,70
], and taking n
=
5000 and Δ
t
=
0.0001.
Time
I
1
I
2
I
3
L
2
×
1
0
3
L
∞
×
1
0
3
0
3.97995
0.810462
2.57901
0.00000
0.00000
4
3.97995
0.810462
2.57901
0.00514014
0.00207598
8
3.97995
0.810462
2.57901
0.0101538
0.00416179
12
3.97995
0.810462
2.57901
0.0149734
0.0060774
16
3.97995
0.810462
2.57901
0.019566
0.00780266
20
3.97995
0.810462
2.57901
0.023929
0.00936437In Figure4 we plot the time evolution of the L
2 and L
∞ norms for the implicit scheme (35) and the nonstandard numerical schemes (40) and (43). We note from Figure 4 that the implicit scheme (35) performs better than the nonstandard numerical scheme (43) in terms of the growth in the error. The nonstandard numerical (40) is the best performer in terms of the growth in the error. From the work of Araújo and Durán [31] we can conclude that because the growth in the errors is linear and the nonstandard numerical schemes preserve the invariants (5)–(7) the schemes (40) and (43) are conservative.Figure 4
Plot showing the time evolution of theL
2 and L
∞ norms for the implicit scheme (35) (solid line), the nonstandard numerical schemes (40) (
-
-
-
) and (43) (
⋯
).In this section we have also shown that the implicit scheme (35), with no nonstandard approximation to the nonlinear term, produces results that do not preserve the invariants (5)–(7), even though at t
=
20 the L
2 and L
∞ norms are significantly better than the results obtained by Avilez-Valente and Seabra-Santos [6] and Daǧ et al. [7]. The two nonstandard numerical schemes (40) and (43) preserve the invariants (5)–(7) and produce L
2 and L
∞ norms at t
=
20 that are significantly better than the results of Avilez-Valente and Seabra-Santos [6] and Daǧ et al. [7].From the results simulated above and captured in Tables3, 4, and 5 we note that for the initial single solitary wave (34) the implicit scheme (35) and the nonstandard numerical schemes (40) and (43) preserve the conserved quantities (5)–(7).In order to test how effectively the implicit scheme (35) and the nonstandard numerical schemes (40) and (43) solve the RLWE equation, we consider the difference u
(
x
,
t
+
Δ
t
)
-
u
(
x
,
t
) between two solutions. We define
(46)
Φ
0,1
,
2
=
u
(
x
,
t
j
+
Δ
t
)
-
u
(
x
,
t
j
)
=
u
j
+
1
-
u
j
,
where the subscripts 0,1
,
2 correspond to the implicit scheme (35) and the nonstandard numerical schemes (40) and (43), respectively. From (36) we have that A
u
j
+
1
=
d. Equation (36) gives u
j
+
1
=
A
-
1
d leading us to write (46) as
(47)
Φ
-
0,1
,
2
=
A
Φ
0,1
,
2
=
d
-
A
u
j
.
The ith row of the vector Φ
-
0,1
,
2 is given by
(48)
[
Φ
-
0,1
,
2
]
i
j
=
d
i
j
-
(
a
i
u
i
-
1
j
+
b
i
u
i
j
+
c
i
u
i
+
1
j
)
.For the implicit scheme (35) where (38) holds for d and (39) holds for a
i, b
i, and c
i we find that
(49)
[
Φ
-
0
]
i
j
=
Δ
x
2
[
u
i
-
1
j
(
1
+
1
2
u
i
-
1
j
)
-
u
i
+
1
j
(
1
+
1
2
u
i
+
1
j
)
]
,
(50)
[
Φ
-
0
]
0
j
=
0
,
[
Φ
-
0
]
n
j
=
0
.For the nonstandard numerical scheme (40) the vector d is given by (41) and a
i, b
i, and c
i are given by (42). Substitution with (47) yields
(51)
[
Φ
-
1
]
i
j
=
Δ
x
2
[
u
i
-
1
j
(
1
+
1
2
u
i
-
1
j
)
-
u
i
+
1
j
(
1
+
1
2
u
i
+
1
j
)
]
,
(52)
[
Φ
-
1
]
0
j
=
0
,
[
Φ
-
1
]
n
j
=
0
.
This is the same as the expression (49).For the nonstandard numerical scheme (43) the vector d is given by (41) and a
i, b
i, and c
i are given by (44). Substitution with (47) yields
(53)
[
Φ
-
2
]
i
j
=
Δ
x
2
[
u
i
-
1
j
+
1
6
(
u
i
-
2
j
+
u
i
-
1
j
+
u
i
j
)
u
i
-
1
j
-
u
i
+
1
j
-
1
6
(
u
i
j
+
u
i
+
1
j
+
u
i
+
2
j
)
u
i
+
1
j
]
,
(54)
[
Φ
-
2
]
1
j
=
Δ
x
2
[
u
0
j
+
1
6
(
u
0
j
+
2
u
1
j
)
u
0
j
-
u
2
j
-
1
6
(
u
1
j
+
u
2
j
+
u
3
j
)
u
2
j
]
,
(55)
[
Φ
-
2
]
n
-
1
j
=
Δ
x
2
[
u
n
-
2
j
+
1
6
(
u
n
-
3
j
+
u
n
-
2
j
+
u
n
-
1
j
)
u
n
-
2
j
-
u
n
j
-
1
6
(
2
u
n
-
1
j
+
u
n
j
)
u
n
j
]
,
(56)
[
Φ
-
2
]
0
j
=
0
,
[
Φ
-
2
]
n
j
=
0
.
Conditions (50), (52), and (56) occur as a result of the Dirichlet boundary conditions.To test how accurately the numerical schemes solve the RLWE we substitute in the solution to the RLWE equation given by (8) that we write as
(57)
u
(
x
i
,
t
j
)
=
u
i
j
=
3
v
sec
h
2
(
k
0
x
i
+
k
1
t
j
+
x
0
)
,
where
(58)
k
0
=
1
2
v
1
+
v
,
k
1
=
-
1
2
v
(
1
+
v
)
,
with the respective differences (49) and (53)–(55). We plot log
⁡
(
max
⁡
Φ
-
0
j
) and log
⁡
(
max
⁡
Φ
-
2
j
) for each t
j for different step lengths Δ
x in Figure 5. We take x
0
=
0 and v
=
0.1 for x
∈
[
-
50,70
]. In Figure 6 we plot the difference log
⁡
(
max
⁡
Φ
-
0
j
)
-
log
⁡
(
max
⁡
Φ
-
2
j
).Figure 5
Figures showing the quantitieslog
⁡
(
max
⁡
Φ
-
0
j
) and log
⁡
(
max
⁡
Φ
-
2
j
) for each t
j, where t
j
∈
[
0,20
], x
0
=
0, and v
=
0.1 for x
∈
[
-
50,70
] for the step lengths Δ
x
=
0.5
,
Δ
x
=
0.1, and Δ
x
=
0.05.Figure 6
Figures showing the quantitieslog
⁡
(
max
⁡
Φ
-
0
j
)
-
log
⁡
(
max
⁡
Φ
-
2
j
) for each t
j, where t
j
∈
[
0,20
], x
0
=
0, and v
=
0.1 for x
∈
[
-
50,70
] for the step lengths Δ
x
=
0.5
,
Δ
x
=
0.1, and Δ
x
=
0.05.It is clear from the expressions for[
Φ
-
0,1
,
2
]
i
j that the amplitude of the oscillations in the trigonometric functions obtained by substituting (57) with (49)–(55), respectively, is proportional to the step length Δ
x. This is seen in the results plotted in Figure 5 but also more clearly in the results plotted in Figure 6. We note that as Δ
x
→
0 the maximum values of the differences plotted in both Figures 5 and 6 tend to zero.We can therefore conclude that for the implicit numerical scheme (35) and the nonstandard numerical schemes (40) and (43) the difference given by (46) is oscillatory where the maximum height of the oscillation is proportional to the spatial step length. We can also conclude that the difference (46) tends to zero as Δ
x
→
0. We can therefore conclude that the first integral (5) that can be written as (46) on differentiating (5) with respect to time is satisfied by the implicit numerical scheme (35) and the nonstandard numerical schemes (40) and (43).


### 3.1. Undular Bore

We next apply the implicit scheme (35) and the nonstandard finite difference schemes (40) and (43) to model the evolution of an undular bore where the initial profile is given by
(59)
u
(
x
,
0
)
=
1
2
u
0
(
1
-
tanh
(
x
-
x
0
α
)
)
.
When simulating the development of an undular bore, we set the first row of the vectord in (38) and (41) to d
0
=
u
(
x
0
,
0
)
=
u
0. We simulate the evolution of the undular bore using the implicit scheme (35) and the nonstandard finite difference schemes (40) and (43). We tabulate the maximum value of u at the times indicated in Tables 6 and 7 for α
=
5 and α
=
2. From the results indicated in Tables 6 and 7 we note that the maximum values of u for all three schemes compare well for α
=
5 and α
=
2. We plot a simulation of the evolution of the undular bore using the nonstandard finite difference scheme (40) in Figures 7 and 8 as this scheme produced the best results as tabulated in Table 4. We choose u
0
=
0.1, x
0
=
0, and α
=
5 for the curves in Figure 7 and α
=
2 for the curves in Figure 8.Table 6
Table showingmax
u
(
x
,
t
) for each of the schemes, where α
=
5.
Time
(35)
(40)
(43)
50
0.117140
0.117136
0.117139
100
0.127146
0.127139
0.127140
150
0.134316
0.134306
0.134306
200
0.140250
0.140236
0.140233
250
0.145455
0.145438
0.145436Table 7
Table showingmax
u
(
x
,
t
) for each of the schemes where α
=
2.
Time
(35)
(40)
(43)
50
0.131598
0.131592
0.131591
100
0.137537
0.137527
0.137526
150
0.142442
0.142429
0.142427
200
0.146848
0.146831
0.146827
250
0.150901
0.150881
0.150879Figure 7
Evolution of initial undular bore (59), where u
0
=
0.1, x
0
=
0, and α
=
5 for Δ
t
=
0.001and n
=
2500 on the interval x
∈
[
-
50,500
] using the nonstandard finite difference scheme (40).Figure 8
Evolution of initial undular bore (59), where u
0
=
0.1, x
0
=
0, and α
=
2 for Δ
t
=
0.001and n
=
2500 on the interval x
∈
[
-
50,500
] using the nonstandard finite difference scheme (40).We note that the nonstandard numerical scheme (40) has captured the development of the undular bore without the formation of numerical instabilities. Our results are comparable with those obtained in Figures 9 and 11 of the paper by Avilez-Valente and Seabra-Santos [6]. We note that the front of the bore and maximum amplitude of the bore in Figures 7 and 8 are the same as those in Figures 9 and 11 of Avilez-Valente and Seabra-Santos [6], respectively.


## 3.1. Undular Bore

We next apply the implicit scheme (35) and the nonstandard finite difference schemes (40) and (43) to model the evolution of an undular bore where the initial profile is given by
(59)
u
(
x
,
0
)
=
1
2
u
0
(
1
-
tanh
(
x
-
x
0
α
)
)
.
When simulating the development of an undular bore, we set the first row of the vectord in (38) and (41) to d
0
=
u
(
x
0
,
0
)
=
u
0. We simulate the evolution of the undular bore using the implicit scheme (35) and the nonstandard finite difference schemes (40) and (43). We tabulate the maximum value of u at the times indicated in Tables 6 and 7 for α
=
5 and α
=
2. From the results indicated in Tables 6 and 7 we note that the maximum values of u for all three schemes compare well for α
=
5 and α
=
2. We plot a simulation of the evolution of the undular bore using the nonstandard finite difference scheme (40) in Figures 7 and 8 as this scheme produced the best results as tabulated in Table 4. We choose u
0
=
0.1, x
0
=
0, and α
=
5 for the curves in Figure 7 and α
=
2 for the curves in Figure 8.Table 6
Table showingmax
u
(
x
,
t
) for each of the schemes, where α
=
5.
Time
(35)
(40)
(43)
50
0.117140
0.117136
0.117139
100
0.127146
0.127139
0.127140
150
0.134316
0.134306
0.134306
200
0.140250
0.140236
0.140233
250
0.145455
0.145438
0.145436Table 7
Table showingmax
u
(
x
,
t
) for each of the schemes where α
=
2.
Time
(35)
(40)
(43)
50
0.131598
0.131592
0.131591
100
0.137537
0.137527
0.137526
150
0.142442
0.142429
0.142427
200
0.146848
0.146831
0.146827
250
0.150901
0.150881
0.150879Figure 7
Evolution of initial undular bore (59), where u
0
=
0.1, x
0
=
0, and α
=
5 for Δ
t
=
0.001and n
=
2500 on the interval x
∈
[
-
50,500
] using the nonstandard finite difference scheme (40).Figure 8
Evolution of initial undular bore (59), where u
0
=
0.1, x
0
=
0, and α
=
2 for Δ
t
=
0.001and n
=
2500 on the interval x
∈
[
-
50,500
] using the nonstandard finite difference scheme (40).We note that the nonstandard numerical scheme (40) has captured the development of the undular bore without the formation of numerical instabilities. Our results are comparable with those obtained in Figures 9 and 11 of the paper by Avilez-Valente and Seabra-Santos [6]. We note that the front of the bore and maximum amplitude of the bore in Figures 7 and 8 are the same as those in Figures 9 and 11 of Avilez-Valente and Seabra-Santos [6], respectively.


## 4. Concluding Remarks

In this paper we have shown the following:(i)
that the modified equation can be used as a criterion for selecting the “best” nonstandard finite difference scheme for solving the RLW equation (1);
(ii)
that the nonstandard numerical scheme we have chosen does indeed perform well by preserving the three conserved quantities (5)–(7) and capturing the single solitary wave solution (8) accurately;
(iii)
that the nonstandard finite difference scheme captures the development of an undular bore for both steep and shallow initial profiles.The implicit nature of the nonstandard numerical schemes (40) and (43) derived in this paper means that at each time step the system matrix A in (36) is recalculated and the system resolved. The numerical results obtained in this paper have been compared to the results obtained by Avilez-Valente and Seabra-Santos [6] who implemented a Petrov-Galerkin scheme and Daǧ et al. [7] who implemented a splitted quadratic B-spline collocation scheme. In the Petrov-Galerkin scheme implemented by Avilez-Valente and Seabra-Santos [6] finite elements are used in both the space and time dimensions. The time-integration is performed by use of a predictor-corrector method. Avilez-Valente and Seabra-Santos [6] summarize the computational efficiency of their method by indicating that their FEM method leads to a tridiagonal system of n
+
1 equations, where n is the number of spatial finite elements. The system matrix remains unchanged during the computation but three inner iterations are required at each time step as part of the predictor-corrector method. Avilez-Valente and Seabra-Santos [6] also comment on the fact that the splitted quadratic B-spline collocation scheme implemented by Daǧ et al. [7] leads to a pentadiagonal system of 2
n
+
4 equations. At each time step there are multiple inner iterations that need to be evaluated depending on the time-integration scheme that has been chosen. Unlike the Petrov-Galerkin method the system matrix must be reassembled after each time step. This makes the splitted quadratic B-spline collocation scheme implemented by Daǧ et al. [7] computationally more expensive than the Petrov-Galerkin scheme implemented by Avilez-Valente and Seabra-Santos [6]. The nonstandard numerical schemes derived in this paper are easy to implement, highly accurate, and computationally a lot cheaper than either the Petrov-Galerkin numerical scheme or the splitted quadratic B-spline collocation scheme.



---

*Source: 754543-2014-04-29.xml*