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

A method of lines approach to the numerical solution of nonlinear wave equations typified by the regularized long wave (RLW) is presented.Themethod 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 with L ∞ and L 2 error norms and the conservation properties of mass, energy, and momentum under the RLW equation.


Introduction
Many researchers have introduced various methods to solve the regularized long wave (RLW) of the form 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][4][5] various numerical studies have been reported based on the finite difference.A recent study in [6] has shown that the meshless kernelbased 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.
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 equation  ≤  ≤  into uniform rectangular mesh by the lines   = ℎ ( = 0, 1, 2, . . ., ) , where ℎ =  −   .The partial derivatives depending on spatial variable in RLW equation are replaced by known central finite difference approximations at point   .This leads to a tridiagonal system of ordinary differential equation By substituting (4) and ( 5) into (1) and taking into account that  0 () =  +1 () = 0 the following system of ODEs is obtained: This system can be written in the following form: where where It can be easily solved by computing  −1 ; this yields a system of ordinary differential equations depending on the time variable  in the form U   =  (U  ) ,  = 1, 2, . . ., .
This system can be solved using some numerical method; in this paper the solutions are obtained using Runge-Kutta method.

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 where  denotes the time step, and for computational purposes the Runge-Kutta method is only useful if it is stable for sufficiently small /ℎ.
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 Δ for stability.The numerical method of lines of RLW equation gives the system of ordinary differential equations First,   is assumed as constant in the nonlinear term, where  = max    .We can now inquire about the eigenvalue of the  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 (, ) with both  and  or  and .Therefore, a product solution is assumed: Further, in accordance with a method proposed by Von Neumann, the function () can be of the form where  is a Fourier number.Substituting of ( 16) and ( 17) in (15) gives Equation ( 18) written in terms of eigenvalue  becomes with  a complex scalar, and denote  =  for numerical discretization with step size  = Δ which can be written as The classical Runge-Kutta (RK4) method contains a segment of the imaginary axis in its stability region and is only mildly dissipative, specifically; |R()| is only a little less than 1 for Note that all  eigenvalues given by (18) have pure imaginary parts, which means that the system of  first-order ODEs are marginally stable.

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 the  2 and  ∞ error norms defined by where   denotes the exact solution and   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 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 (  )   and its first derivative ((  )  )   can be computed from (4).

Motion of the Solitary Waves
Case 1.The equation represents a single soliton with amplitude 3 and constant speed +1.For the purpose of comparing with the earlier work, the computations are done for the parameters  = 1,  = 1, and  0 = 0. Consider  = 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  = 18.The error norms at  = 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].
Case 2. In the second experiment, a smaller solitary wave of amplitude 0.09 ( = 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   error norms remain less than  2 = 0.637 × 10 −3 and  ∞ = 0.233 × 10 −3 , whereas they are  2 = 0.638 × 10 −3 and  ∞ = 0.233 × 10 −3 for the finite difference method [7] and  2 = 0.572 × 10 −3 and  ∞ = 0.364 × 10 −3 for the finite element method [9].Moreover, it can be seen that the error norms ( 2 and  ∞ ) are somewhat small compared with those quoted by others [7,9].
The invariants  1 ,  2 , and  3 for RLW at different values of time  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  1 = 2.109407,  2 = 0.127302, and  3 = 0.388806, it can be seen from Table 7 that the quantities  1 ,  2 , and  3 obtained using present method are very close to the corresponding exact values and compared with the results obtained by [8].

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 where   = 3  sec ℎ 2 [  ( −   )],  = 1, 2, and ), 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 (, ) = (, ) = 0.
Figure 3 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  (see Table 8).The computed values of the invariants  1 ,  2 , and  3 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.
Case 2. The interaction of two negative solitary waves with initial condition (25) was studied over the region 0 ≤  ≤ 120 with  1 = 0.6,  2 = 0.8,  1 = 82,  2 = 67,  =  = 1, and Δ = 0.2, Δ = 0.5.The invariants  1 ,  2 , and  3 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. are carried out for the evolution of solitary waves for  = 0.04.We take the smaller space step ℎ = 0.01 and time step Δ = 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  = 12.The values of  1 ,  2 , and  3 are given in Table 10; each value is satisfactorily constant.

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 norms  2 and  ∞ .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 which  2 and  ∞ cannot be evaluated, we verify the accuracy of the method by computing the conservation quantities  1 ,  2 , and  3 .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.

Figure 1 :
Figure 1: Graphs of initial solution and solution at time = 20 and amplitude = 0.3.

Figure 5 :
Figure 5: The motion of Maxwellian initial condition initial solution and solution at time = 0 and 12.
(12)298 × 10 −5 1.33000 × 10 −4 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)with locally stable Runge-Kutta methods.The stability region contains a half circle

Table 8 :
Invariant values for the interaction of two positive solitary waves over the region 0 ≤  ≤ 120.

Table 9 :
Invariant values for the interaction of two negative solitary waves over the region 0 ≤  ≤ 120.

Table 10 :
Invariant values for Maxwellian initial condition,  = 0.04. the solution range in order to maintain accuracy.The range here is doubled from −40 ≤  ≤ 60 to −80 ≤  ≤ 120.The error norms for  = 0.03 are recorded in Table6for different values of time .With the range −40 ≤  ≤ 60, Δ = 0.1, excellent results were obtained,  2 and  ∞