Mathematical Problems in Engineering

Volume 2013 (2013), Article ID 804317, 8 pages

http://dx.doi.org/10.1155/2013/804317

Research Article

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

H. O. Bakodah and M. A. Banaja

Department of Mathematics, Sciences Faculty for Girls, King Abdulaziz University, Saudi Arabia

Received 23 November 2012; Revised 20 February 2013; Accepted 26 March 2013

Academic Editor: Slim Choura

Copyright © 2013 H. O. Bakodah and M. A. Banaja. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

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:

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

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 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 ; this yields a system of ordinary differential equations depending on the time variable in the form

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

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

with locally stable Runge-Kutta methods. The stability region contains a half circle

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 .

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; 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.

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

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).

It is well known that (1) has analytic solution of the form

where is the wave velocity, the width and initially centered at .

4.1. Motion of the Solitary Waves

Case 1. The equation represents a single soliton with amplitude and constant speed . For the purpose of comparing with the earlier work, the computations are done for the parameters , , and . Consider , 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 . The error norms at 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    and over the region .

4
0.69452 × 10⁻⁴
0.264411 × 10⁻⁴

8
0.137889 × 10⁻³
0.534151 × 10⁻⁴

12
0.204502 × 10⁻³
0.792604 × 10⁻⁴

16
0.268849 × 10⁻³
0.103317 × 10⁻³

20
0.33080 × 10⁻³
0.125611 × 10⁻³

Table 2: Absolute errors for (1) with    and over the region .

Present method
MOL [7]

−49.5
3.44444 × 10⁻⁹
0.000000000

−39.5
3.01068 × 10⁻⁸
0.000000000

−29.5
5.43620 × 10⁻⁹
0.000000000

−19.5
2.01194 × 10⁻⁷
1.00000 × 10⁻⁶

−9.5
3.91714 × 10⁻⁶
11.0000 × 10⁻⁵

0.5
4.39242 × 10⁻⁵
2.03000 × 10⁻⁴

10.5
2.31414 × 10⁻⁵
3.59000 × 10⁻³

20.5
1.43949 × 10⁻⁴
1.94700 × 10⁻³

30.5
1.9499 × 10⁻⁴
2.38500 × 10⁻³

40.4
4.67298 × 10⁻⁵
1.33000 × 10⁻⁴

Table 3: and errors for (1) with    and over the region .

Present method
Finite difference
Finite element

4
0.6945 × 10⁻⁴
0.26441 × 10⁻⁴
0.1200 × 10⁻³
0.05000 × 10⁻³
0.1176 × 10⁻³
0.5622 × 10⁻³

8
0.13788 × 10⁻³
0.53415 × 10⁻⁴
0.2300 × 10⁻³
0.09000 × 10⁻³
0.23108 × 10⁻³
0.1088 × 10⁻³

12
0.20450 × 10⁻³
0.79260 × 10⁻⁴
0.3400 × 10⁻³
0.1400 × 10⁻³
0.3374 × 10⁻³
0.1544 × 10⁻³

16
0.26884 × 10⁻³
0.10331 × 10⁻³
0.45000 × 10⁻³
0.18000 × 10⁻³
0.43784 × 10⁻³
0.19377 × 10⁻³

20
0.33080 × 10⁻³
0.12561 × 10⁻³
0.55000 × 10⁻³
0.21000 × 10⁻³
0.53223 × 10⁻³
0.2272 × 10⁻³

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 , , and [8].

Table 4: Invariant values for a single solitary wave: amplitude = 0.3, and over the region .

Present method
Cubic spline collocation method [8]

0
3.97990
0.810461
2.57901
3.97992
0.810462
2.57900

4
3.97995
0.810461
2.57901
3.97995
0.810276
2.57839

8
3.97997
0.810461
2.57901
3.97997
0.810276
2.57839

12
3.97998
0.810461
2.57901
3.97998
0.810276
2.57839

16
3.97999
0.810461
2.57901
3.97998
0.810276
2.57839

20
3.97996
0.810461
2.57901
3.97996
0.810276
2.57839

For , the relative changes in the invariants are very small (the quantities , , and change by less than , , and , respectively, by ). Initial solution and solitary wave profile at 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 () 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 to . The error norms for are recorded in Table 6 for different values of time . With the range , , excellent results were obtained, and error norms remain less than and , whereas they are and for the finite difference method [7] and and 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    and over the region .

4
0.51431 × 10⁻⁵
0.14805 × 10⁻⁵

8
0.10281 × 10⁻⁴
0.29961 × 10⁻⁵

12
0.15408 × 10⁻⁴
0.45367 × 10⁻⁵

16
0.20518 × 10⁻⁴
0.60908 × 10⁻⁵

20
0.25606 × 10⁻⁴
0.76472 × 10⁻⁵

Table 6: and errors for (1) with    and over the region .

Present method
Finite difference
Finite element

4
0.14956 × 10⁻³
0.12290 × 10⁻³
0.1500 × 10⁻³
0.1230 × 10⁻³
0.3983 × 10⁻³
0.22450 × 10⁻³

8
0.32000 × 10⁻³
0.16600 × 10⁻³
0.3200 × 10⁻³
0.1660 × 10⁻³
0.50528 × 10⁻³
0.21706 × 10⁻³

12
0.466159 × 10⁻³
0.17900 × 10⁻³
0.4670 × 10⁻³
0.1790 × 10⁻³
0.5357 × 10⁻³
0.20919 × 10⁻³

16
0.56635 × 10⁻³
0.18454 × 10⁻³
0.5670 × 10⁻³
0.1850 × 10⁻³
0.54934 × 10⁻³
0.20165 × 10⁻³

20
0.63748 × 10⁻³
0.23300 × 10⁻³
0.6380 × 10⁻³
0.2330 × 10⁻³
0.57247 × 10⁻³
0.36498 × 10⁻³

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 , , and , 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, and over the region .

Present method
Cubic spline collocation method [8]

0
2.10700
0.127302
0.38806
2.10702
0.127302
0.388804

4
2.10836
0.127302
0.38807
2.10838
0.127294
0.388781

8
2.10927
0.127302
0.38807
2.10928
0.127294
0.388782

12
2.10980
0.127302
0.38807
2.10981
0.127294
0.388782

16
2.10981
0.127302
0.38807
2.10981
0.127294
0.388782

20
2.10899
0.127302
0.38806
2.10900
0.127294
0.388781

For , the relative changes in the quantities , , and are changed by less than , , and , respectively, by . Initial solution and solitary wave profile at 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

where , , and , with .

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 .

The RLW parameters used in the numerical simulation are , , , , , and . This wave moves across the space interval and the time interval .

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 .

Present method
[9]

0
37.9165
120.480
744.081
37.916520
120.52280
744.0815

4
37.9170
120.471
743.984
37.911700
121.16020
736.9443

8
37.9172
120.445
743.676
37.859750
119.73170
728.51730

12
37.9173
120.203
740.702
37.696670
119.63340
725.72360

16
37.9174
120.230
741.050
37.529160
119.41850
725.83940

20
37.9174
120.441
743.654
37.647300
119.80410
727.19480

25
37.9175
120.462
743.918
38.050100
119.83550
727.439200

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 with , , , , , and , . 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 .

Present method
[9]

0
−45.035
336.232
−1077.82
−45.034960
336.96020
−1077.824

4
−45.035
336.328
−1077.97
−45.035000
334.25310
−1048.6640

8
−45.035
336.014
−1068.16
−45.035080
332.98740
−1041.3070

12
−45.035
336.172
−1077.67
−45.035100
334.45920
−1048.8220

16
−45.035
336.295
−1077.80
−45.035110
334.23630
−1048.8480

20
−45.049
336.284
−1077.67
−45.035060
334.52330
−1048.7200

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

and boundary conditions

are carried out for the evolution of solitary waves for . We take the smaller space step and time step . For the Maxwellian develops into a single solitary wave plus a small developed oscillating tail as shown in Figure 5 at time . The values of , , and are given in Table 10; each value is satisfactorily constant.

Table 10: Invariant values for Maxwellian initial condition, .

Present method
[9]

0
1.77245
1.30332
4.78327
1.772454
1.303446
4.783268

3
1.77245
1.30466
4.79137
1.772500
1.341671
4.918245

6
1.77245
1.30457
4.79124
1.790807
1.343856
4.922796

9
1.77245
1.30436
4.79039
1.775987
1.343681
4.917124

12
1.77245
1.30417
4.78957
1.767234
1.377576
4.913779

(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.

References

D. H. Peregrine, “Calculation of the development of an undular bore,” Journal of Fluid Mechanics, vol. 25, no. 2, pp. 321–330, 1966.

T. B. Benjamin, J. L. Bona, and J. J. Mahony, “Model equations for long waves in nonlinear dispersive systems,” Philosophical Transactions of the Royal Society of London A, vol. 272, no. 1220, pp. 47–78, 1972.

J. C. Eilbeck and G. R. McGuire, “Numerical study of the regularized long-wave equation. I. Numerical methods,” Journal of Computational Physics, vol. 19, no. 1, pp. 43–57, 1975.

P. C. Jain and L. Iskandar, “Numerical solutions of the regularized long-wave equation,” Computer Methods in Applied Mechanics and Engineering, vol. 20, no. 2, pp. 195–201, 1979.

C. Qianshun, G. B. Wang, and B. L. Guo, “Conservative scheme for a model of nonlinear dispersive waves and its solitary waves induced by boundary motion,” Journal of Computational Physics, vol. 93, no. 2, pp. 360–375, 1991.

Y. Dereli, “Numerical solutions of the MRLW equation using meshless kernel based method of lines,” International Journal of Nonlinear Science, vol. 13, no. 1, pp. 28–38, 2012.

S. Kutluay and A. Esen, “A finite difference solutions of the regularized long wave equation,” Mathematical Problems in Engineering, vol. 2006, Article ID 85743, 14 pages, 2006.

D. Irk, I. Dag, and A. Dogan, “Numerical integration of the RLW equation using cubic splines,” The ANZIAM Journal, vol. 47, no. 1, pp. 131–142, 2005.

K. R. Raslan, “A computational method for the regularized long wave (RLW) equation,” Applied Mathematics and Computation, vol. 167, no. 2, pp. 1101–1118, 2005.

H.-O. Kreiss and G. Scherer, “Method of lines for hyperbolic differential equations,” SIAM Journal on Numerical Analysis, vol. 29, no. 3, pp. 640–646, 1992.

P. J. Olver, “Euler operators and conservation laws of the BBM equation,” Mathematical Proceedings of the Cambridge Philosophical Society, vol. 85, no. 1, pp. 143–160, 1979.

G. Griffiths and W. Schiesser, Traveling Wave Analysis of Partial Differential Equations, Academic Press, New York, NY, USA, 2012.