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

Abstract

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

1. Introduction

Many researchers have introduced various methods to solve the regularized long wave (RLW) of the form 𝑢_𝑡+𝑢_𝑥+𝜖𝑢𝑢_𝑥−𝜇𝑢_𝑥𝑥𝑡=0, (1) where 𝜖 and 𝜇 are positive constants.

The RLW equation was introduced to describe the behavior of the undular bore by Peregrine [1]. The RLW equation describes a lot of important physical phenomena, such as shallow water waves and ionic waves.

An analytical solution for the RLW equation was found under the restricted initial and boundary condition [2], but this solution is not very useful; the availability of numerical method is essential. In previous work [3–5] various numerical studies have been reported based on the finite difference. A recent study in [6] has shown that the meshless kernel-based method of lines provides a highly accurate and efficient method for the modified regularized long wave equation. In this study, the method of lines (MOL) is applied to obtain a numerical solution for RLW with a constant grid of central finite differences and fourth Runge-Kutta approximation in time. The performance of the method was tested on two known model problems.

2. The Method of Lines Solution of the RLW Equation

We shall consider (1) with the following initial and boundary conditions: 𝑢()(𝑥)𝑢()()𝑥,0=𝑔,𝑎≤𝑥≤𝑏,𝑎,𝑡=0,𝑢𝑏,𝑡=0, (2) where 𝑔(𝑥) is prescribed function.

The basic idea of the method of lines consists in a discretization of the space variable and in substitution of partial derivatives by finite-difference expressions; the time variable keeps a continuous form, and solutions above the lines parallel with the axis of this variable are computed.

Let us subdivide the solution domain of the RLW equation 𝑎≤𝑥≤𝑏 into uniform rectangular mesh by the lines 𝑥_𝑖=𝑖ℎ(𝑖=0,1,2,…,𝑁),whereℎ=𝑏−𝑎𝑁. (3) 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 𝜕𝑢_𝑖=𝑢𝜕𝑥^𝑖+1−𝑢_𝑖−1,𝜕2ℎ₂𝑢_𝑖𝜕𝑥²=𝑢_𝑖+1−2𝑢_𝑖+𝑢_𝑖−1ℎ². (4)(5) By substituting (4) and (5) into (1) and taking into account that 𝑢₀(𝑡)=𝑢_𝑁+1(𝑡)=0 the following system of ODEs is obtained: −ℎ²+2𝜇𝑑𝑢₁𝑑𝑡+𝜇𝑑𝑢₂=ℎ𝑑𝑡2𝑢₂+𝜖2𝑢²₂𝜇𝑑𝑢_𝑖−1−ℎ𝑑𝑡²+2𝜇𝑑𝑢_𝑖𝑑𝑡+𝜇𝑑𝑢_𝑖+1ℎ𝑑𝑡=−2𝑢_𝑖−1+𝜖2𝑢²_𝑖−1+ℎ2𝑢_𝑖+1+𝜖2𝑢²_𝑖+1𝜇𝑖=2,3,…,𝑁−1,𝑑𝑢_𝑁−1−ℎ𝑑𝑡²+2𝜇𝑑𝑢_𝑁ℎ𝑑𝑡=−2𝑢_𝑁−1+𝜖2𝑢²_𝑁−1. (6) This system can be written in the following form: 𝐀𝑑𝐔𝐔𝑑𝑡=𝐁𝐔,_𝑖(0)=𝑔𝑥_𝑖, (7) where 𝐔(𝑡)=𝑢₁(𝑡),𝑢₂(𝑡),…,𝑢_𝑛(𝑡)^𝑇,𝐔_𝑖(0)=𝑢₁(0),𝑢₂(0),…,𝑢_𝑛(0)^𝑇,𝐔(0)=𝑔𝑥₁,𝑔𝑥₂,…,𝑔𝑥_𝑛^𝑇, (8) where−ℎ𝐴=²𝜇ℎ+2𝜇𝜇−²𝜇ℎ+2𝜇𝜇−²⋯⋯−ℎ+2𝜇000000𝜇00⋮⋱⋮000...000²ℎ+2𝜇𝜇0𝜇−²𝜇ℎ+2𝜇𝜇−²,ℎ+2𝜇𝐵=2𝑢₂+𝜖2𝑢²₂,𝑢₃+𝜖2𝑢²₃−𝑢₁+𝜖2𝑢²₁,…,𝑢_𝑖+1+𝜖2𝑢²_𝑖+1−𝑢_𝑖−1+𝜖2𝑢²_𝑖−1,…,−𝑢_𝑁−1+𝜖2𝑢²_𝑁−1^𝑇. (9)It can be easily solved by computing 𝐴⁻¹; this yields a system of ordinary differential equations depending on the time variable 𝑡 in the form 𝑑𝐔_𝑖𝐔𝑑𝑡=𝑓_𝑖,𝑖=1,2,…,𝑁. (10) This system can be solved using some numerical method; in this paper the solutions are obtained using Runge-Kutta method.

3. Stability Analysis of the Numerical Method

The stability of the method of lines for partial differential equations represents the most important factor for their solution; hence, it is a critical factor that should be handled carefully. The importance lies in its unique ability of judging acceptable solution for the given equation. The criticality depends on the nature of the eigenvalues of the matrix representation in connection with their values. Kreiss and Scherer in [10] derived the conditions of local stability of Runge-Kutta methods when applied to hyperbolic partial differential equations.

Theorem 1. Consider a well-posed initial boundary value problem for a system of partial differential equations. The space derivatives are discretized in such a way that the resulting system of ODEs 𝑑𝑢𝑑𝑡=𝑄𝑢 (11) 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 |𝜎(𝑤)|≤𝑅,Real𝜎≤0, (13) 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 𝑢_𝑡+𝑢_𝑥+𝜖𝑢𝑢_𝑥−𝜇𝑢_𝑥𝑥𝑡=0 (14) gives the system of ordinary differential equations 𝜇𝑑𝑢_𝑖−1−ℎ𝑑𝑡²+2𝜇𝑑𝑢_𝑖𝑑𝑡+𝜇𝑑𝑢_𝑖+1=ℎ𝑑𝑡2𝑢_𝑖+1+𝑢_𝑖−1+𝜖𝑢_𝑖𝑢_𝑖+1−𝑢_𝑖−1,𝑖=1,2,3,…,𝑁. (15) 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: 𝑢(𝑥,𝑡)=𝑐𝜓(𝑡)𝜑(𝑥). (16) Further, in accordance with a method proposed by Von Neumann, the function 𝜑(𝑥) can be of the form 𝜑(𝑥)=𝑒^𝑗𝑘𝑥,𝑗=√−1, (17) where 𝑘 is a Fourier number.

Substituting of (16) and (17) in (15) gives 𝑐𝜇𝑑𝜓𝑒𝑑𝑡^{𝑗𝑘𝑥𝑖+1}−ℎ²+2𝜇𝑑𝜓𝑒𝑑𝑡^{𝑗𝑘𝑥𝑖}+𝜇𝑑𝜓𝑒𝑑𝑡^{𝑗𝑘𝑥𝑖−1}=ℎ1+𝜖𝛾2𝜓𝑒^{𝑗𝑘𝑥𝑖+1}−𝜓𝑒^{𝑗𝑘𝑥𝑖−1}. (18) Equation (18) written in terms of eigenvalue 𝜆 becomes 𝑑𝜓=ℎ()𝑑𝑡1+𝜖𝛾𝑖sin𝑘Δ𝑥−4𝜇sin²(𝑘Δ𝑥/2)+ℎ²𝜓,ℎ()𝜆=1+𝜖𝛾sin𝑘Δ𝑥−4𝜇sin²(𝑘Δ𝑥/2)+ℎ²𝑖 (19)(20) with 𝜆 a complex scalar, and denote 𝑧=𝑝𝜆 for numerical discretization with step size 𝑝=Δ𝑡 which can be written as 𝜓_𝑛+1=ℛ(𝑧)𝜓_𝑛,𝑛=0,1,2,…,𝑁. (21) The classical Runge-Kutta (RK4) method contains a segment of the imaginary axis in its stability region and is only mildly dissipative, specifically; |ℛ(𝑧)| is only a little less than 1 for |𝐼𝑚(𝑧)|≤1.2.

Note that all 𝑁 eigenvalues given by (18) have pure imaginary parts, which means that the system of 𝑁 first-order ODEs are marginally stable.

4. Numerical Results

To illustrate the effectiveness of the method, numerical results portraying a single soliton solution for the RLW equation are presented. In order to show how good the numerical solutions are in comparison with exact ones, we will use the 𝐿₂ and 𝐿_∞ error norms defined by 𝐿₂=ℎ^𝑁_𝑖=0𝑢_𝑒𝑖−𝑢_𝑎𝑖^21/2,𝐿_∞=max_𝑖𝑢_𝑒𝑖−𝑢_𝑎𝑖, (22) 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 𝐶₁=^𝑏_𝑎𝑢𝑑𝑥≅ℎ^𝑁_𝑖=1𝑢^𝑛_𝑖,𝐶₂=^𝑏_𝑎𝑢²+𝜇𝑢_𝑥²𝑑𝑥,≅ℎ^𝑁_𝑖=1𝑢_𝑁^𝑛_𝑖²+𝜇𝑢_𝑥^𝑛_𝑖²𝐶₃=^𝑏_𝑎𝑢³+3(𝑢)²𝑑𝑥,≅ℎ^𝑁_𝑖=1𝑢_𝑁^𝑛_𝑖³+3(𝑢)^𝑛_𝑖². (23) 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).

It is well known that (1) has analytic solution of the form 𝑢(𝑥,𝑡)=3𝑐secℎ²𝑘𝑥−𝜈𝑡−𝑥₀, (24) where 𝜈=1+𝜀𝑐 is the wave velocity, 𝑘=(1/2√𝜇)𝜀/(1+𝜀𝑐) the width and initially centered at 𝑥₀.

4.1. 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. 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].

Table 1 

𝐿₂ and 𝐿_∞ errors for (1) with 𝑐=0.1  Δ𝑥=0.1,Δ𝑡=0.1 and over the region −40≤𝑥≤60.

Table 2 

Absolute errors for (1) with 𝑐=0.1  𝑛=199 and over the region −50≤𝑥≤50,𝑡=18.

Table 3 

𝐿₂ and 𝐿_∞ errors for (1) with 𝑐=0.1  Δ𝑥=0.1,Δ𝑡=0.1 and over the region −40≤𝑥≤60.

The error norms at each time are smaller than those given in [7, 9]. The invariant values 𝐶₁, 𝐶₂, and 𝐶₃ are given in Table 4 and compared with the results obtained by [8]. From this table, it can be seen that the conservation properties (using the present method) are very good in comparison with the analytical values for the invariants given by 𝐶₁=6𝑐/𝑘=3.9799497, 𝐶₂=(12𝑐²/𝑘)+(48𝑐²𝜇/5)=0.81046249, and 𝐶₃=(36𝑐²/𝑘)+(144𝑐³/5𝑘)=2.579007 [8].

Table 4 

Invariant values for a single solitary wave: amplitude = 0.3, Δ𝑥=0.125,Δ𝑡=0.1 and over the region −40≤𝑥≤60.

For 𝑐=0.1, the relative changes in the invariants are very small (the quantities 𝐶₁, 𝐶₂, and 𝐶₃ change by less than 1.03×10⁻⁵, 1.49×10⁻⁶, and 3×10⁻⁶, respectively, by 𝑡=20). Initial solution and solitary wave profile at 𝑡=20 are depicted in Figure 1.

Figure 1 

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

Case 2. In the second experiment, a smaller solitary wave of amplitude 0.09 (𝑐=0.03) has been modeled, and the results of simulation are given in Table 5. As the amplitude of the solitary wave is reduced the pulse broadens, and it may be necessary to increase the solution range in order to maintain accuracy. The range here is doubled from −40≤𝑥≤60 to −80≤𝑥≤120. The error norms for 𝑐=0.03 are recorded in Table 6 for different values of time 𝑡. With the range −40≤𝑥≤60, Δ𝑡=0.1, excellent results were obtained, 𝐿₂ and 𝐿_∞ error norms remain less than 𝐿₂=0.637×10⁻³ and 𝐿_∞=0.233×10⁻³, whereas they are 𝐿₂=0.638×10⁻³ and 𝐿_∞=0.233×10⁻³ for the finite difference method [7] and 𝐿₂=0.572×10⁻³ and 𝐿_∞=0.364×10⁻³ for the finite element method [9]. Moreover, it can be seen that the error norms (𝐿₂ and 𝐿_∞) are somewhat small compared with those quoted by others [7, 9].

Table 5 

𝐿₂ and 𝐿_∞ errors for (1) with 𝑐=0.03  Δ𝑥=0.1,Δ𝑡=0.1 and over the region −80≤𝑥≤120.

Table 6 

𝐿₂ and 𝐿_∞ errors for (1) with 𝑐=0.03  Δ𝑥=0.1,Δ𝑡=0.1 and over the region −40≤𝑥≤60.

The invariants 𝐶₁, 𝐶₂, and 𝐶₃ 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 𝐶₁=2.109407, 𝐶₂=0.127302, and 𝐶₃=0.388806, it can be seen from Table 7 that the quantities 𝐶₁, 𝐶₂, and 𝐶₃ obtained using present method are very close to the corresponding exact values and compared with the results obtained by [8].

Table 7 

Invariant values for a single solitary wave: amplitude = 0.09, Δ𝑥=0.125,Δ𝑡=0.1 and over the region −40≤𝑥≤60.

For 𝑐=0.03, the relative changes in the quantities 𝐶₁, 𝐶₂, and 𝐶₃ are changed by less than 4.1×10⁻⁴, 1.5×10⁻⁷, and 3.1×10⁻⁷, respectively, by 𝑡=20. Initial solution and solitary wave profile at 𝑡=20 are depicted in Figure 2.

Figure 2 

Graphs of initial solution and solution at time = 20 and amplitude = 0.09.

4.2. The Interaction of Solitary Wave

The interaction of two positive solitary waves is studied by using the initial condition given by the linear sum of two separate solitary waves of various amplitudes 𝑢(𝑥.0)=𝑢₁+𝑢₂, (25) where 𝑢_𝑖=3𝑐_𝑖secℎ²[𝑘_𝑖(𝑥−𝑥_𝑖)], 𝑖=1,2, and 𝑐_𝑖=4𝑘²_𝑖/(1−4𝑘²_𝑖), 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.

The RLW parameters used in the numerical simulation are 𝑘₁=0.4, 𝑘₂=0.3, 𝑥₁=15, 𝑥₂=35, 𝜀=𝜇=1, and Δ𝑡=0.01. This wave moves across the space interval 0≤𝑥≤120 and the time interval 0≤𝑡≤25.

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 𝐶₁, 𝐶₂, and 𝐶₃ are satisfactorily constant compared with the corresponding invariant values obtained by [9]. The results of the present method at different times are shown in Figure 3, where it is possible to observe that the higher amplitude solitary wave passes through the smaller wave with no change in its waveform.

Table 8 

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

Figure 3 

The motion of two positive solitary waves before and after the interaction at time = 0 and 25.

Case 2. The interaction of two negative solitary waves with initial condition (25) was studied over the region 0≤𝑥≤120 with 𝑘₁=0.6, 𝑘₂=0.8, 𝑥₁=82, 𝑥₂=67, 𝜀=𝜇=1, and Δ𝑡=0.2, Δ𝑥=0.5. The invariants 𝐶₁, 𝐶₂, and 𝐶₃ recorded in Table 9; the computed values are satisfactorily constant compared with the corresponding invariant values obtained by [9]. The two negative solitary waves before and after the interaction are shown in Figure 4.

Table 9 

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

Figure 4 

The motion of two negative solitary waves before and after the interaction at time = 0 and 20.

4.3. Maxwellian Initial Condition

The numerical solutions of the RLW equation with the Maxwellian initial condition 𝑢(𝑥.0)=exp−(𝑥−7)² (26) and boundary conditions 𝑢(0,𝑡)=𝑢(50,𝑡)=0 (27) 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 𝐶₁, 𝐶₂, and 𝐶₃ are given in Table 10; each value is satisfactorily constant.

Table 10 

Invariant values for Maxwellian initial condition, 𝜇=0.04.

(a)

(b)

Figure 5 

The motion of Maxwellian initial condition initial solution and solution at time = 0 and 12.

5. Conclusions

MOL is developed to obtain numerical solutions of the RLW equation. The efficiency of the proposed method was tested on problems from the literature, and its accuracy was examined by the error norms 𝐿₂ 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 𝐿₂ and 𝐿_∞ cannot be evaluated, we verify the accuracy of the method by computing the conservation quantities 𝐶₁, 𝐶₂, and 𝐶₃. 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.