A FINITE DIFFERENCE SOLUTION OF THE REGULARIZED LONG-WAVE EQUATION

A linearized implicit finite difference method to obtain numerical solution of the onedimensional regularized long-wave (RLW) equation is presented. The performance and the accuracy of the method are illustrated by solving three test examples of the problem: a single solitary wave, two positive solitary waves interaction, and an undular bore. The obtained results are presented and compared with earlier work.


Introduction
In this study, we will consider the one-dimensional RLW equation with the physical boundary conditions U → 0 as x → ±∞, where t is time, x is the space coordinate, U(x,t) is the wave amplitude, and ε and μ are positive parameters.The RLW equation (1.1) was first introduced by Peregrine [1] to describe the development of an undular bore.This equation is one of the most important nonlinear wave equations which can be used to model a large number of problems arising in various areas of applied sciences [2,3].The RLW equation has been solved analytically for a restricted set of boundary and initial conditions.Therefore, the numerical solution of the RLW equation has been the subject of many papers.Various numerical techniques particularly including finite difference [4][5][6][7][8], finite element [9][10][11][12][13][14][15][16][17][18][19], and spectral [20][21][22][23] methods have been used for the solution of the RLW equation.
In this paper, we have used a linearized implicit finite difference method to investigate the motion of a single solitary wave, development of two positive solitary waves interaction, and an undular bore for the RLW equation (1.1).

Method of solution
For the numerical treatment, the spatial variable x of the problem is restricted over an interval a ≤ x ≤ b.In this study, we consider the RLW equation (1.1) with the homogeneous boundary conditions U(a,t) = 0, t > 0, U(b,t) = 0, t > 0, (2.1) and the initial condition where f (x) is a prescribed function.
The solution domain a ≤ x ≤ b, t > 0 is divided into subintervals Δx in the direction of the spatial variable x and Δt in the direction of time t such that x i = iΔx, i = 0(1)N (NΔx = b − a); t j = jΔt, j = 0(1)J, and the numerical solution of U at the grid point (iΔx, jΔt) is denoted by U i, j .
In the finite difference method, the dependent variable and its derivatives are approximated by the finite difference approximation.This approximation will lead to either a single explicit equation or a system of difference equations.Applying the classical implicit finite difference method to nonlinear problems normally gives nonlinear system of equations which cannot be solved directly.
Equation (1.1) can be written as Using the forward difference approximation for ∂U/∂t, the Crank-Nicolson difference approximation for ∂U/∂x and ∂U 2 /∂x, and the central difference approximation for ∂ 2 U/∂x 2 at the point (i, j + 1),

∂U ∂t
respectively, (2.3) yields the system of algebraic equations S. Kutluay and A. Esen 3 for i = 1(1)N − 1 and j = 0(1)J with a truncation error of O(Δt) + O(Δx) 2 .The scheme is a nonlinear system of equations in U i, j+1 and it needs to use an iteration technique to evaluate the solution.
Using the central difference operator δ defined by δ x U i, j = U i+1, j − U i−1, j , (2.5) can be written as By Taylor expansion of U 2 i, j+1 about the point (i, j) we obtain Hence in terms of order Δt, U 2 i, j+1 ), and taking (2.6), with some manipulations, leads to (i = 1(1)N − 1) a system of linear equations for W i .This approximation is second order in both space and time as regards truncation error.Obviously, the solution at the ( j + 1)th time level is obtained from (2.8) as U i, j+1 = U i, j + W i .Since the stability parameter Δt/(Δx) 2 depends not only on the form of the finite difference scheme (2.9) but also generally upon the solution U(x,t) being obtained, the complications and difficulties may arise in the analysis of stability.In order to show how good the numerical solutions are in comparison with the exact ones, we will use the L 2 and L ∞ error norms defined by (2.10) 4 Regularized long-wave equation

Numerical examples and results
All computations were executed on a Pentium 4 PC in the Fortran code using double precision arithmetic.The RLW equation (1.1) satisfies only three conservation laws given as which respectively correspond to mass, momentum, and energy [24].In the simulations the invariants I 1 , I 2 , and I 3 are monitored to check the conservation of the numerical scheme.For the computation of U x in (3.1), we used a central finite difference approximation.
To implement the performance of the method, three test problems will be considered: the motion of a single solitary wave, the interaction of two positive solitary waves, and the undular bore.

The motion of a single solitary wave.
We first consider (1.1) with the boundary conditions U → 0 as x → ±∞ and the initial condition The exact solution of this problem is This solution corresponds to the motion of a single solitary wave with amplitude 3c and width k, initially centered at x 0 , where v = 1 + εc is the wave velocity and k = (1/ 2)(εc/μv) 1/2 .This solution will also be used over an interval a ≤ x ≤ b.For this problem the theoretical values of the invariants are [14] which are recorded throughout the simulations.For the purpose of comparing with the earlier work, all computations are done for the parameters ε = 1, μ = 1, and x 0 = 0. Table 3.1 displays a comparison of the values of the invariants and error norms obtained by the present method with those obtained using the cubic finite difference method developed by Jain et al. [6] and implemented by Gardner et al. [10] for c = 0.1.As it is seen from the table, the numerical values of invariants obtained from (3.1) are in very good agreement with their analytical values obtained from (3.4).The quantities in the invariants remain almost constant during the computer run.For the proposed finite difference method at times t = 0 and t = 20, change in I 1 is 0.5 × 10 −4 , and I 2 and I 3 are exact up to the last recorded digit, whereas for the cubic finite difference method, they are 0.43227, 0.086883, and 0.2746, respectively.The error norms at each time obtained by the present method are smaller than those given in [6,10].For the present method at t = 20, the error norms are L 2 = 0.55 × 10 −3 and L ∞ = 0.21 × 10 −3 , whereas they are L 2 = 196.1 × 10 −3 and L ∞ = 67.35× 10 −3 for the cubic finite difference method.In Table 3.2 the time evolution of the invariants I 1 , I 2 , and I 3 , and of the error norms L 2 and L ∞ for c = 0.03, is compared with the cubic finite difference method [6,10].Again the present method produces good results.
The rates of convergence for the proposed numerical method in space sizes Δx m and time steps Δt m can be calculated by respectively [18].The convergence rates computed by the present method for values of space size Δx m and a fixed value of the time step Δt are recorded in Table 3.3.It is clearly seen that the scheme provides remarkable reductions in convergence rates for the smaller space sizes.
6 Regularized long-wave equation  The profiles of the solitary waves at times t = 0 and t = 20 and the error distributions of the analytical and numerical solutions at t = 20 for c = 0.1 with the range −40 ≤ x ≤ 60 and for c = 0.03 with the range −80 ≤ x ≤ 120, Δx = 0.125 and Δt = 0.1 are shown in Figure 3.1.For c = 0.1, the amplitude is 0.3 at time t = 0 while it is 0.299919 at time t = 20 (Figure 3.1(a)) and so the relative change in the amplitude is about 0.027%.It is seen that the maximum error is about between −4 × 10 −3 and 4 × 10 −3 (Figure 3.1(b)).For c = 0.03, the amplitude is 0.09 at time t = 0 while it is 0.089997 at time t = 20 (Figure 3.1(c)) and so the relative change in the amplitude is about 0.0033%.It is observed that the maximum error is about between −6 × 10 −4 and 6 × 10 −4 (Figure 3.1(d)).

The interaction of two positive solitary waves.
We secondly consider (1.1) with the boundary conditions U → 0 as x → ±∞ and the initial condition [17] where A j = 4k 2 j /(1 − 4k 2 j ) (j = 1,2).For the simulation, all computations are done for the parameters k 1 = 0.4, 3, and Δt = 0.1 over the region 0 ≤ x ≤ 120.The experiment was run from t = 0 to t = 25 to allow the interaction to take place.Figure 3.2 shows the interaction of two positive solitary waves.As it is seen from the figure, at t = 0 a solitary wave with larger amplitude is on the left of the other solitary wave with smaller amplitude.The larger wave catches up with the smaller one as the time increases.At t = 0, the amplitude of the larger solitary wave is 5.33338 while the amplitude of the smaller one is 1.68598, whereas at t = 25, the amplitude of the larger solitary wave is 5.30235 at the point x = 86.7 while the amplitude of the smaller one is 1.67157 at the point x = 69.9.An oscillation of small amplitude trailing behind the solitary waves was observed.In order to see this oscillation occurring behind the waves in Figure 3.2 at time t = 25, the scale of the figure is magnified as in Figure 3.3.It is clearly seen that an oscillation of amplitude ∼ 2.2 × 10 −2 is trailing behind the solitary waves.Table 3.5 displays a comparison of the values of the invariants obtained by the present method with those obtained in [17].It is observed that the obtained values of the invariants remain almost constant during the computer run.At times t = 0 and t = 25, the relative changes in the invariants I 1 , I 2 , and I 3 for the present method are respectively 2.558 × 10 −3 %, 6.647 × 10 −3 %, and 9.797 × 10 −3 % whereas they are 0.352%, 0.570%, and 2.237% for the cubic B-spline collocation finite element method given in [17].It is clearly seen that each of the conserved quantities obtained by the present method is very well preserved.

The undular bore.
As our last test problem, we consider (1.1) with the physical boundary conditions U → 0 as x → ∞ and U → U 0 as x → −∞, and the initial condition where U(x,0), denotes the elevation of the water surface above the equilibrium level at time t = 0, U 0 represents the magnitude of the change in water level which is centered on x = x 0 , and d measures the steepness of the change.Under the above physical boundary conditions, the invariants I 1 , I 2 , I 3 are not constant but increase linearly throughout the simulation at the following rates [14]: respectively.
For the simulation, all computations are done for the parameters ε = 1.5, μ = 1/6, U 0 = 0.1, x 0 = 0, Δx = 0.24, Δt = 0.1, and d = 2,5 in the region −36 ≤ x ≤ 300.The simulation is run until time t = 250, and the values of the quantities I 1 , I 2 , I 3 with the position and amplitude of the leading undulation for the steep slope d = 2 and the gentle slope d = 5 are recorded in Table 3

Conclusion
A linearized implicit finite difference method was presented to obtain numerical solutions of the RLW equation.The efficiency of the method was tested on three numerical experiments of wave propagation: the motion of a single solitary wave, the development of two positive solitary waves interaction, and an undular bore, and its accuracy was examined by the error norms L 2 and L ∞ .The obtained results show that the error norms are reasonably small and the conservation properties are all very good.The results also suggest that the present method whose application is easier than many other numerical techniques such as finite element and spectral methods can be applied to a large number of physically important nonlinear wave problems with success.

. 6 .
The numerical values of variations in quantities I 1 , I 2 , I 3 are obtained as M 1 = 0.107500, M 2 = 0.010992, M 3 = 0.034096 for d = 2 and M 1 = 0.107500, M 2 = 0.010992, M 3 = 0.034101 for d = 5 which are in good agreement with the theoretical values M 1 = 0.105000, M 2 = 0.010667, M 3 = 0.033075 obtained from (3.8).The values of I 1 , I 2 , and I 3 increase linearly according to the values of M 1 , M 2 , and M 3 , respectively.The amplitudes of the leading undulation for d = 5 and d = 2 are 0.17710 and 0.18158, respectively.S. Kutluay and A. Esen 11

Figure 3 .
Figure 3.4 illustrates the undular bore profiles at t = 100 and t = 250 for the gentle slope d = 5 and the steep slope d = 2.As it can be seen that from the figure, the number of undulations formed increases with the decrease of d from d = 5 to d = 2.The number of undulations also increases with the increase of t, as expected.

Table 3 .
4displays the computed convergence rates for various values of time step Δt j and a fixed value of the space size Δx.Again a noticeable decrease in convergence rates is observed when the time step decreases.S. Kutluay and A. Esen 7

Table 3 .
5. Invariants for the interaction of two positive solitary waves.