A New Functional Iterative Algorithm for the Regularized Long-Wave Equation Using an Integral Equation Formalism

A fundamental question in computational nonlinear partial differential equations is raised to discover if one could construct a functional iterative algorithm for the regularized long-wave (RLW) equation (or the Benjamin–Bona–Mahony equation) based on an integral equation formalism? Here, the RLW equation is a third-order nonlinear partial differential equation, describing physically nonlinear dispersive waves in shallow water. For the question, the concept of pseudo-parameter, suggested by Jang (Commun Nonlinear Sci Numer Simul 43:118–138, 2017), is introduced and incorporated into the RLW equation. Thereby, dual nonlinear integral equations of second kind involving the parameter are formulated. The application of the fixed point theorem to the integral equations results in a new (semi-analytic and derivative-free) functional iteration algorithm (as required). The new algorithm allows the exploration of new regimes of pseudo-parameters, so that it can be valid for a much wider range (in the complex plane) of pseudo-parameter values than that of Jang (2017). Being fairly simple (or straightforward), the iteration algorithm is found to be not only stable but accurate. Specifically, a numerical experiment on a solitary wave is performed on the convergence and accuracy of the iteration for various complex values of the pseudo-parameters, further providing the regions of convergence subject to some constraints in the complex plane. Moreover, the algorithm yields a particularly relevant physical investigation of the nonlinear behavior near the front of a slowly varying wave train, in which, indeed, interesting nonlinear wave features are demonstrated. As a consequence, the preceding question may be answered.


Introduction
As a starting point of this paper, a fundamental (as well as theoretical) question is raised with regard to computational partial differential equations: is it possible to construct a functional iterative algorithm for the regularized long-wave (RLW) equation (also termed the Benjamin-Bona-Mahony equation) [1] through an integral equation formalism? Here, η and c 0 denote wave elevation and characteristic velocity (of the square root of the product of gravitational acceleration g and water depth h 0 ) c 0 = gh 0 , (1.2) respectively. Physically, the nonlinear partial differential equation (1.1) models a one-way wave propagation of water waves over a horizontal bottom valid for (weakly) nonlinear and fairly long waves, being widely studied and frequently used for simulating nonlinear dispersive waves in shallow water in applied sciences and coastal engineering. Being an alternative model to the Korteweg-de Vries (KdV) equation, the RLW (1.1) was firstly used by Peregrine for its advantages in numerical computations [2] and examined extensively by Benjamin et al. [3]. Here, the term "regularized" is used, because of the much better nature of the dispersion relation between frequency ω BBM and wave number k when compared to that of the KdV equation. This is the main reason why (1.1) is much more referred for a one-way propagation model.
There exists only few about the analytical studies on the nonlinear equation (1.1); e.g., Benjamin et al. [3] found the analytical solution for the equation just under the restricted initial and boundary conditions [3]. Thus, finding its numerical solutions is of practical importance and the availability of accurate and efficient numerical methods is essential [4]. Over the last few decades, significant efforts have been made for developing useful numerical procedures for (1.1). The finite difference approach would be the most typical numerical scheme, approximating (1.1) with difference equations. For instance, Eilbeck and McGuire [5] investigated numerical solitary wave solutions of the RLW equation to show that they exhibit true soliton behavior, being stable on collision with other solitary waves [5]. Bhardwaj and Shankar [6] employed quintic spline technique and splitting method to develop a new-finite difference method, which is used to model solitary wave motion and undular bore development [6]. Cai [7] derived a 6-point multisymplectic Preissman scheme from its Bridges' multisymplectic form, where backward error analysis is implemented and the performance and the efficiency of the new scheme are illustrated by solving several test examples [7]. However, finite element schemes also have been one of the most commonly used numerical methods, which are extensively employed in the area of computational partial differential equations. Challenging issues are found; for example, Gardner et al. [8] used a B-spline finite element method to solve the RLW equation numerically, which is involved with a Galerkin method with quadratic B-spline finite elements [8]. Dag et al. [9] applied Cubic B-spline functions to develop a collocation method for the nonlinear numerical solutions, which is used to model solitary waves, undular bore development and wave generation [9]. Now we return to the question at the beginning about how to construct a functional iterative algorithm of the RLW equation based on an integral equation formalism. For that, we will introduce the notion of pseudo-parameter proposed by Jang [10], which is incorporated into (1.1), for establishing a nonlinear integral formalism involving the parameter. This enables us to obtain dual nonlinear integral equations, equivalent to (1.1). Here, it is pointed out that the integral equations are not singular but regular, which inherently differ from not only the boundary integral equations resulted from boundary element methods but also the weak integral formulations from finite element methods. We then apply Banach fixed point theorem [10,11] to the integral equations and derive a (required) functional iteration algorithm for integrating (1.1), implying that the answer to the question may be found.
The iteration algorithm derived in this paper is semi-analytic and derivative-free, which provides a fast convergence speed for achieving highly accurate numerical solutions as is illustrated in the numerical experiment on moving solitary waves. The numerical results are also compared to those of the conventional excellent methods. Here, it should be emphasized that the derived algorithm works even for a much wider range of pseudo-parameter values than that of Jang [10]; e.g., the pseudo-parameter can be a complex number (rather than a real number), in the present study. Thus, a complex (numerical) calculation is naturally required for performing the iteration, which leads to a complex numerical solution. This point is in sharp contrast to the previous article by Jang [10], where pseudo-parameter values were only restricted to some part of the real numbers and thus a real (not complex) numerical calculation was required. This may be best illustrated by showing the regions of convergence (ROCs) subject to some constraints in the complex plane, presented in numerical experiment section. In particular, for a further relevant investigation for nonlinear dispersive water waves, the algorithm is also found to be useful at observing the nonlinear behavior near a wave train front. Here, the nonlinear feature of the front wave is confirmed; i.e., the front wave has more peaked crests and flatter troughs, as it becomes steep, and the wave moves faster with larger amplitude.

Integral Formalism
Here, we will present an integral formalism which corresponds to the RLW equation (1.1). We first modify the governing equation (1.1), following the concept of pseudo-parameter [10]. To be specific, we introduce a pseudo-parameter α ∈ C in the present study, where C denotes the set of complex numbers, and add the product of α · η to both sides of (1.1), leading to The forcing term ϕ in (2.1) involves η as well as its derivative η x , i.e., It immediately follows that (2.1), combined with (2.2), is still equivalent to the original Eq. (1.1), regardless of complex values of α. The parameter α introduced above will play a crucial role in formulating an integral formalism for (1.1), which is necessary for obtaining a system of nonlinear integral equations, as shall be discussed later. We begin with applying two successive integral transforms.

Integral Transformations
Here, associated with (2.1), we aim to find the successive Fourier and Laplace transforms of η. First of all, we apply Fourier integral transform to both sides of the partial differential equation of (2.1) to find an ordinary differential equation. We denote byη the Fourier transform F of η with respect to space x as [10,12] for a parameter k. Then, the inverse Fourier transform F −1 is expressed as If we assume a localized wave motion of η in x, i.e., we have the (well known) differentiation properties of Fourier transform This allows (2.1) to be converted into an (first order) ordinary differential equation forη, Next, the Laplace transform L ofη is to be performed on (2.7) with respect to time t, denoted byη * , i.e., for a parameter s [10,12]. This results simply in an algebraic equation forη * (i.e., L [F(η)]), Here, the Laplace-transform property of differentiation is employed together with initial condition; in which η 1 designates an initial wave profile, i.e., η 1 (x) = η(x, 0). (2.9) is readily solved forη * to giveη With the definition of a (complex-valued) functionω BBM such that the expression (2.11) may be further shortened; (2.11) admits the partial fractions representation in the simple form in terms of theω BBM

Inverse Integral Transformations
We start with taking the inverse Laplace transform L −1 on the first term of the right hand side of (2.13), which is for a ∈ C from table [10,12]. Performing further the inverse Fourier transform F −1 on the above result, being denoted by I 1 , would lead to by Fourier convolution theorem, where the notation * means Fourier convolution. I 1 in (2.16) can be calculated as which describes physically the superposition of (complex) elementary wave solutions e −iθ(x,t;k) with a wave spectrumη 1 , i.e., the Fourier transform of the initial profile η 1 (x); Next, we move on to the inverse Laplace transform L −1 on the second term of (2.13), where the formula 1/(s − a) = L exp(at) for a ∈ C and Laplace convolution theorem are used [10,12]. And we take the inverse Fourier transform F −1 on the above result, giving a triple integral I 2 , by the use of Fourier formula (2.3), its inverse (2.4) and (2.12).

Remark 2
The integral I 2 of (2.23) can be also expressed as a single integral in the form Then, the integrand ζ physically represents the superposition of elementary wave solutions e −iθ(x,t−τ ;k) with time shift τ , which has a wave spectrumφ(k, τ )/(1+h 2 0 k 2 /6) parametrized by the τ :

Establishing Integral Formalism
This subsection is devoted to the realization of an integral formalism from the results in the previous subsection. This is simply achieved by recalling that η is recovered by taking the two successive inverse Fourier and Laplace transforms on both sides of (2.13), i.e., for a pseudo-parameter α ∈ C, η is symbolically represented as

Derivation of Functional Iteration Algorithm
In this section, we shall derive a new functional iteration algorithm for solving the RLW equation (1.1). To this end, (1.1) is transformed into a system of nonlinear integral equations by using the results of the previous section. And then, we employ Banach contraction theorem [10,11], applied to the system of integral equations, which makes us derive an (numerical) iterative strategy (i.e., a functional iteration algorithm for (1.1)).

Dual Integral Equations
With the integral formalism (2.26) established in Sect. 2.3, we will obtain nonlinear integral equations. For that, let us substitute (2.2) into (2.26) to find where ψ indicates a new variable, corresponding physically to wave slope η x . And, we differentiate the above with respect to x, giving The pair of (3.1) and (3.2) then constitutes two coupled nonlinear integral equations of second kind for the two unknown functions of η and ψ. The integral equations of the form (3.1) and (3.2) may be abbreviated, with some preliminaries of background material of function spaces. Denote by B T the collection of complex-valued continuous and bounded functions defined on R × [0, T ] for T > 0, where R is the set of real numbers. Analogously, let B 0 be the collection of real-valued continuous and bounded functions defined on R. Then, associated with (3.1) and (3.2), we can define an integral operator J : for any η 1 ∈ B 0 . In a similar way, we define K : for any η ∈ B T . By virtue of the integral operators J in (3.3) and K in (3.4), we can write (3.1) as The two operator equations of (3.5) and (3.6) can be considered as an abbreviated form for the integral equations described by (3.1) and (3.2). Here, the new operators J x and K x have the expressions and respectively.

Deriving Functional Iterative Algorithm
In this subsection, we present the derivation of a functional iteration algorithm to solve (1.1) by means of the integral equations obtained in the previous subsection. Let us define a (integral) nonlinear operator for any η, ψ ∈ B T , so that (3.5) can be more shortened into Then, the integral equations (3.1) and (3.2) are much more abbreviated as Here, the nonlinear operator T has the form (the superscript t denotes the transpose) in terms of a new variable η = (η, ψ) t , in which the operator T x stands for the right hand side of (3.6), i.e., Based on the fact that η is invariant under the mapping T in (3.11), we would like to propose a functional iteration algorithm for η n ≡ (η n , ψ n ) t , where η n , ψ n ∈ B T , with the initial guess by appealing to Banach fixed point theorem [10,11]. That is, (3.14) is a recurrence relation which recursively defines a functional sequence η n , n = 0, 1, 2, ....

Remark 3
Being semi-analytic and derivative-free iterative strategy, the iteration (3.14) yields a solution η of (1.1) as well as its derivative η x (= ψ) as a bonus, as the number of iteration n increases, provided that (3.14) converges.
We recall that the model (1.1) physically describes a one-way wave propagation of nonlinear dispersive waves η(x, t) over a horizontal bottom as discussed in Introduction, where η(x, t) is real-valued, if a real-valued initial condition η 1 (x) = η(x, 0) is taken into account. Suppose that η is a (real-valued) wave solution of the initial value problem for (1.1) with the real-valued initial condition. Then, the η should also satisfy the α-parameterized partial differential equations of (2.1) and (2.2) together with the initial condition, for any pseudoparameter α ∈ C, due to the fact that (2.1) and (2.2), for α ∈ C, are equivalent to (1.1) as mentioned at the beginning of Sect. 2. And, the η even further should satisfy the system of integral equations (3.1) and (3.2), or alternatively, the integral operator equation (3.11), because (3.1) and (3.2) describe an (equivalent) integral equation formulation for (1.1) combined with the initial condition.
If T in the integral operator equation (3.11) has a contractive nature, the η then may be expressed as a limit of a sequence {η n }, n = 0, 1, 2, ..., generated by (3.14) (recalling that η is invariant under T) and the sequence converges (See Banach fixed point theorem [11]). The sequence is a Cauchy sequence, which has a unique limit in a Banach space (such as our solution space), because every Banach space is topologically a Hausdorff space (Roman [11]); i.e., η n → η as n → ∞. (3.16) Remark 4 (3.16) indicates that a limit of a sequence {η n }, n = 0, 1, 2, ...by (3.14) is unique and real-valued, if a real-valued initial condition η 1 (x) = η(x, 0) is imposed. However, if n is finite, the output η n may be complex-valued, whose real part converges to the true solution of (1.1) as n → ∞, while its imaginary part goes to zero. That is, followed by the fact that η n = Re (η n ) + i · Im (η n ) → η as n → ∞ by (3.16) and the limit η is real-valued.

Iterative Algorithm for the Special Case α = 0
Notice that it is a complex value that the pseudo-parameter α of the iteration (3.14) involves. This demands that the output of the numerical iteration (3.14) should be complex-valued for a finite n. However, if we set α = 0 in (3.14), the output becomes real-valued. Suppose that the pseudo-parameter α ∈ C vanishes in (3.14). Then, I 1 in (2.17) or (2.20) reduces to where θ represents the phase functionθ in (2.19) when α = 0, i.e., Here, the kernel in (3.18) results from the Euler formula e −iθ = cos θ − i sin θ (3.20) and the anti-symmetric property because ω BBM and (thus) θ are odd in k. Similarly, from (2.23) because of the anti-symmetry (in k) (3.23) Therefore, (3.14) becomes for n = 0, 1, 2, ..., respectively, due to (3.18) and (3.22), when α = 0; with emphasis on that I 1 in (2.17) and I 2 in (2.23) remain pure real. It is clear that the functional iteration consisting of (3.24) and (3.25) derived above is a recurrence relation which does not involve a pseudo-parameter α ∈ C at all (i.e., it is pseudo-parameter free.) and thus both η n and ψ n are real-valued. That is, the output of the iteration of (3.24) and (3.25) is real-valued, as mentioned above in this subsection.

Remark 5
As a special case of (3.14), the pair of (3.24) and (3.25), which corresponds to α = 0, is the simplest (or the most efficient) recurrence relation among our α-parameter family of (3.14); in fact, the pair is an iterative strategy whose output is real-valued. Even though the pair is the simplest, it produces the highest levels of accuracy in the present iterative strategy (3.14), as can be seen later (see Sect. 4).

Numerical Performance
This subsection concerns numerical experiments performed on solitary wave propagations to examine whether the functional iteration algorithm (3.14) works. Thereby, we will investigate the effect of the pseudo-parameter α ∈ C introduced in Sect. 2 on the numerical performance of the iteration (e.g., the iteration's convergence speed and accuracy). We first introduce three appropriate dimensionless variables, defined by Then, the (original) physical partial differential equation (1.1) is transformed into a dimensionless RLW equation [13] (see "Appendix A") It is well known that (4.2) allows a solitary wave (or soliton) moving to the right in the form where 3d indicates dimensionless amplitude, v wave velocity(= d + 1), and 2π/κ wave width: In original physical variables, (4.3) may be re-expressed as By iterating the recursion of (3.14) together with the (zero) initial guess (3.15), we now numerically simulate the solitary wave of (4.3) or (4.5) using the initial condition imposed by Recall that (3.14) involves the pseudo-parameter α which is in general complex. For this reason, a complex domain calculation is generally demanded for the numerical simulation as was discussed before. However, we start with the simplest case of (3.24) and (3.25) (when α = 0).

Simulations for α = 0
We iterate the pair of (3.24) and (3.25), corresponding to case α = 0, with the initial condition (4.6) or (4.7) to obtain the solitary wave solution of (4.3) or (4.5). For the simulation, the dimensionless amplitude 3d is chosen as 0.3, and gravitational acceleration g and water depth h 0 as 9.80665 m/s 2 and 1.0 m, respectively. Here, the usual Simpson integration rule is used to estimate the integrals appearing in the pair (but, with the trapezoidal rule in time), where we discretize dimensionless temporal interval 0 < t * < T into (N − 1) equal length-parts t * = T /(N −1) and dimensionless spatial interval a < x * < b into (M −1) equally spaced panels x * = (b − a)/(M − 1), giving x * i (these points will be used for error estimation later) for positive integers i (1 ≤ i ≤ M),    Figure 4 illustrates that the discrete root mean square norm Err 2 and the maximum error Err ∞ [14] decay rapidly as the number of iterations n increases;  Table 2 Numerical values given by the proposed iteration and other conventional excellent methods for Err 2 in (4.9), Err ∞ in (4.10) and invariants in (4.11-4.13) at specified dimensionless time t * = 20: α = 0 The proposed The conventional x * = 1/8, t * = 0. and 10) where two norms · 2 and · ∞ are taken over the x * variable only, with dimensionless time t * being kept fixed (t * = 10, 20). The values of errors are tabulated in Table 1 for various specified dimensionless times. Note that the order of accuracy is of 10 −6 , being accurate. They are also compared to other results of the conventional excellent methods as is Table 2.
It may be worth estimating numerical values for the following three conservative (invariant) quantities corresponding to mass, momentum and energy [14,15]: (4.12) and where, we review that (1.1) is known to have exactly the three conservation laws (4.11-4.13). Being compared to the analytical results, the numerical values of the invariants of C 1 , C 2 , C 3 in (4.11-4.13) by the proposed iteration of (3.24) and (3.25) are shown for various dimensionless times in Table 1. For other two types of increments, e.g., x * = t * = 0.50 and x * = t * = 0.25, the corresponding errors and invariants are presented as well in Table 2.

Simulations for α Other than Zero
In Sect. 4.1.1, we examined the special case of the iterative algorithm (3.14), i.e., the pair of (3.24) and (3.25) which corresponds to (3.14) when α = 0, for the numerical experiment on a solitary wave. However, in the present subsection, we extend the pseudo-parameter α to complex numbers (of course, which include real numbers) to explore its effect on the numerical results. Using the same integration numerical rules as was employed in Sect. 4.1.1, we perform (3.14) together with the initial condition (4.6) or (4.7) to simulate again the solitary wave solution of (4.3) or (4.5). Here, the dimensionless computational domains are chosen to be −35 < x * < 35 and 0 < t * < 4 respectively, but with the larger increments of x * = 0.5 and t * = 0.2 (than before). It should be noted that, for non-zero complex numbers α, the iterative solutions η n (x, t) for finite integers n produced by (3.14) are in general complex-valued (as was pointed out in Sect. 3.2), because the (two) concerning operators J in (3.3) and K in (3.4) are complexvalued. So, generally, the real parts of iterative solutions converge to a true solution as the number of iteration n increases, while the imaginary parts of iterative solutions converge to zero (See Remark 4).
However, we find it interesting to observe an exceptional case, where the iterative solutions η n (x, t) for finite integers n produced by (3.14) are proved to be all real-valued (not complexvalued), especially when α is pure real (because the iteration (3.14) becomes essentially a real-valued (recurrence) equation if α is real. See "Appendix B"). Bearing in mind this, the iterative algorithm (3.14) yields Fig. 5, which exhibits only convergence behaviors of the real parts of η n (x, t), when the pseudo-parameter is selected to be pure real; α = 1. Here, the real parts of the iterative solutions converge to a true solution as the number of iteration n increases. Figures 6 and 7 present the convergence characteristics of (3.14) for the pure imaginary case, α = i. The convergence behavior (real part) of Fig. 6 is analogous to that of the  Fig. 5, and the imaginary parts of the iterative solutions are seen to converge to zero in Fig. 7. Figures 8 and 9 account for the nature of solution convergence for α =1+i. Figures 10, 11 and 12 plot the iterative solitary wave profiles with dimensionless time kept fixed (t * =4). Here, note again that Fig. 10 only shows convergence behaviors of the real parts, because the iterative solutions are real-valued with α pure real, as was discussed above.
As was discussed earlier, the pseudo-parameter α was artificially introduced at the beginning of this study and it may be crucial in that it affects the nature of convergence and accuracy of the iteration (3.14) derived (or proposed) in this paper. Indeed, Fig. 13 reveals the characteristics of convergence and accuracy depending on the parameter α for a fixed dimensionless time t * = 4, where the real and imaginary parts of errors are defined as 14) Here, η * ex is the exact solution of the moving solitary wave of (4.3). It is interesting to see that the smaller the norm of α tends to contribute to the faster iteration's convergence and more accurate solutions.

Region of α ∈ C Satisfying Requirement
We have, so far, performed the numerical simulation of a solitary wave to find the numerical performances of the algorithm, when given a specified α. However, in this subsection, we  which is an inequality of a performance requirement with a tolerance ε. The two regions (ε = 0.1), plotted in the complex plane C as in Fig. 14a, b, respectively, both look symmetric with respect to the real axis R and similar. In addition, the former is bigger than the latter. This is because the algorithm produces a less accurate solution with higher dimensionless amplitude, when given a computational (discretization) condition; therefore, ultimately, there Here, it is also observed that the α ∈ C, considered for the numerical experiment in the above subsections, are inside the two regions (i.e., α ∈ 1 , 2 ), which, in fact, are located near the origin in the complex plane C (see the enlarged images in Fig. 14). The regions can be also considered as the ROCs subject to the constraint (4.16), due to the fact that they clearly show parts of ROCs of the present iterative algorithm. If α were limited to real numbers as was done in Jang [10], then, the corresponding ROCs subject to the constraint (4.16) would be the intersections of 1 and 2 with the real axis R. Because

Nonlinear Behavior Near Wave Train Front
In this subsection, we will explore the nonlinear behavior near the front of a slowly varying wave train due to a local disturbance by using the present algorithm. This is a fascinating matter in water waves, in fact, which has been significantly examined in the past (e.g., in relation with tsunami). However, the classical (theoretical) study was restricted essentially to the linear approximation; cf. see Kajiura [21] who found the leading wave (linear) solution in the form of an Airy function [21].
For simulating a wave train, we first consider an initial value problem of the governing equation (1.1) and the initial conditions 17) or, alternatively, for constants A i > 0 and B > 0, which are known as the Maxwellian conditions [22]. And then, the algorithm (3.24) and (3.25), which corresponds to α = 0, is iterated to solve the above initial value problem for three different A i > 0, i.e., A 1 = 0.01, A 2 = 0.20 and A 3 = 0.27, where the concerning computational domains are selected as 20m < x < 140 m (M = 481, x = 1/3 m) and 0sec< t < 35sec (N = 211, t = 1/6 sec); equivalently, -49< x * < 343( x * ≈ 0.82), 0< t * < T = 268( t * ≈ 1.28) in dimensionless form. Figure 16 indicates a typical convergence behavior of the algorithm, say for A 3 = 0.27, in which, similarly as in the previous subsection (e.g., in Fig. 2), qualitatively reasonable iterative solutions η * n are found to be achieved at the early stage of iteration. Furthermore, it is clearly seen that the Maxwellian profile (4.17) propagates into space and time in a dispersive manner.
We will next determine where to stop the (iterative) algorithm (3.24) and (3.25), based on the norm difference e between η * n+1 and η * n in this study; we would like to stop the iteration with a (tolerance) bound δ and an end time T in numerical experiment (= 268, see the (just) above). Figure 17 uncovers plots of e against iteration n for A i > 0, i = 1, 2, 3, which illustrates that smaller amplitudes A i make e decay much more rapidly, as n increases. If we set, say δ = O(10 −9 ) in (4.19), we then stop the iteration algorithm at n = 6 for A 1 = 0.01, because e has order O(10 −9 ) at n = 6 for A 1 = 0.01 as shown in Fig. 17. Similarly, we stop the iteration at n = 25 for A 2 = 0.20 and at n = 32 for A 3 = 0.27, respectively, where e = O(10 −9 ) for both cases A 2 and A 3 (see Fig. 17). These results imply that our numerically converged solutions become η * 6 , η * 25 and η * 32 for A 1 = 0.01, for A 2 = 0.20 and A 3 = 0.27, respectively, if δ = O(10 −9 ) is chosen for (4.19).

Remark 6
No analytical solution of the present initial value problem has been found, so there may be no way to estimate quantitative errors of the iterative solutions such as (4.9) and (4.10). Figure 18 shows the numerically converged solutions of the η * 6 , η * 25 and η * 32 for three different amplitudes, A i > 0, i = 1, 2, 3. These describe physically the development of nonlinear wave train fronts. Here, we clearly observe the nonlinear behaviors near wave train It is found that the nonlinear behavior of the leading wave (or the wave front) is in a good agreement with the linear one, when the smallest amplitude, i.e., A 1 = 0.01, of (4.17) is considered, as expected. However, there is a sensible discrepancy between them for the other amplitudes. Especially, the nonlinear feature of the leading wave (or the wave front within blue-colored rectangles) is most revealed for the case A 3 = 0.27, where we can find steep waves with peaked crests and flat troughs like a cnoidal wave. And we further discover the interesting nonlinear phenomenon that high-amplitude waves move faster than waves of low amplitude.

Concluding Remarks
The basic idea underlying this paper may be characterized by the question which concerns how to construct a functional iterative algorithm for the RLW equation via an integral equation formalism. For resolving the question, we have introduced the pseudo-parameter concept [10], merged into the original RLW equation. This provides the (pseudo-parameterized) dual regular integral equations, making us derive a recurrence relation using fixed point iteration. It is a (required) new functional iterative algorithm for integrating the RLW equation.
The functional iterative algorithm derived in this paper is semi-analytic and derivativefree, which is straightforward to apply in that the only thing required is just the use of regular numerical integration in every iterative process (because the concerning integral equations are regular). For checking the algorithm's performance, we have made a numerical experiment with solitary wave propagations to examine not only the convergence properties and accuracy characteristics but also the ROCs subject to some constraints in the complex plane. Notably, the algorithm is also shown to be good at investigating the evolution of a slowly varyng nonlinear wave train and the nonlinear behavior near the front of the wave train, which demonstrates fundamental and interesting physical phenomena of nonlinear dispersive water waves.
It is noted that the pseudo-parameter, artificially introduced here, may be regarded as an important parameter, which has several roles. For the first, it connects the original RLW equation with the dual integral equations. Second, it is related with the iteration's performance, e.g., the convergence speed and accuracy of the iteration algorithm. Finally, we would like to close this paper by stressing the point that the values of the pseudo-parameter in this paper extend over a much wider range than those in the previous study in Jang [10]. This may be clearly confirmed by presenting diagrams showing the ROCs subject to some constraints plotted in the complex plane.