Explicit finite difference methods for the EW and RLW equations

Applied Mathematics and Computation 179 (2006) 622–638 www.elsevier.com/locate/amc

Explicit finite difference methods for the EW and RLW equations J.I. Ramos Room I-320-D, E.T.S. Ingenieros Industriales, Universidad de Ma´laga, Plaza El Ejido, s/n, 29013-Ma´laga, Spain

Abstract An extensive assessment of the accuracy of explicit finite difference methods for the solution of the equal-width (EW) and regularized long-wave (RLW) equations is reported. Such an assessment is based on the three invariants of these equations as well as on the magnitude of the errors of the numerical solution and has been performed as a function of the time step and grid spacing. Two of the methods presented here make use of three-point, fourth-order accurate, finite difference formulae for the first- and second-order spatial derivatives. Two methods are based on the analytical solution of secondorder ordinary differential equations which have locally exponential solutions, and the fourth technique is a standard finite difference scheme. A linear stability analysis of the four methods is presented. It is shown that, for the EW and RLW equations, a compact operator method is more accurate than locally exponential techniques that make use of compact operator approximations. The latter are reported to be more accurate than exponential techniques that employ second-order accurate approximations, and, these, in turn, are more accurate than the standard explicit method. Ó 2005 Elsevier Inc. All rights reserved. Keywords: Explicit finite difference methods; Compact operators; Locally exponential techniques; Equal-width equation; Regularized long-wave equation

1. Introduction The equal width (EW) and regularized long-wave (RLW) equations are special cases of the (damped) generalized regularized long-wave (GRLW) equation which can be written as ut þ aux þ ðup Þx ¼ luxx þ dutxx ;

1 < x < 1;

t > 0;

ð1Þ

where u, t and x denote the amplitude, time and spatial coordinate, respectively, a, , l and d are greater than or equal to zero and p P 2 is a constant. Eq. (1) includes the (first-order) linear wave equation for  = l = d = 0 and a 5 0, the (first-order) nonlinear wave equation for a = l = d = 0 and  5 0, the onedimensional heat transfer equation for a =  = d = 0 and l 5 0, the one-dimensional linear advection– diffusion equation for  = d = 0, a 5 0 and l 5 0, the one-dimensional nonlinear advection–diffusion equation for a = d = 0,  5 0 and l 5 0, the one-dimensional nonlinear Burgers’ equation for a = d = 0,  5 0 E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.12.003

J.I. Ramos / Applied Mathematics and Computation 179 (2006) 622–638

623

and l 5 0, the EW equation for a = l = 0,  = 1, d 5 0 and p = 2, and the RLW or the Benjamin–Bona– Mahony [1] for a = 1, p = 2,  5 0, d 5 0 and l = 0. On the other hand, the GRLW equation corresponds to a = 1, l = 0, a 5 0,  5 0, d 5 0 and p > 2. Eq. (1) also includes the dissipative EW, RLW and GRLW if l 5 0 in Eq. (1). It must be noted that the RLW equation can be transformed into the EW one by means of the transformation v = a + 2u. The RLW equation was formulated by Peregrine [2,3] as an alternative to the Korteweg-de Vries (KdV) equation for studying soliton phenomena and as a model for small-amplitude long waves on the surface of water. A rather interesting property of the RLW equation is that the collision of two solitary waves may result in the creation of secondary solitary waves or sinusoidal solutions. This phenomenon is somewhat analogous to what happens in subatomic physics where the collisions of particles create another particles and/or radiation. Therefore, a study of both the EW and RLW equations provide the opportunity of investigating the creation of secondary solitary waves and/or radiation to get insight into the corresponding processes of particle physics [4,5]. The collisions of solitary waves of the EW and RLW equations differ from the collision of solitons of the KdV equation in that, in the latter, solitons collide, interact and emerge unchanged from the collision except for a phase shift. By way of contrast, the numerical solution of the RLW equation indicates that two solitary waves can effectively pass through each other but with diminished amplitude, and create secondary waves. In addition, a notch in the negative amplitude wave may appear and this notch develops into a secondary negative amplitude wave, positive and negative amplitude waves may move in opposite directions, multiple secondary waves of positive and negative amplitude may be formed, and annihilation of positive and negative waves and the creation of a sinusoidal residual referred to as radiation, may occur. Previous numerical studies of the RLW equation include the use of quartic-B spline collocation methods for the spatial discretization together with an explicit fourth-order accurate Runge–Kutta technique for time integration [6], least-squares finite element techniques with linear space–time elements and cubic B-spline finite element techniques that use the Galerkin method [7,8], Fourier-pseudospectral methods for periodic boundary conditions [9–12], mixed finite element methods based on the introduction of an intermediate variable [13], operator-splitting techniques with cubic [14,15], quintic [16] B-splines and Crank–Nicolson time discretization, Adomian’s decomposition method [17,18], and perturbation (WKBJ) techniques [19]. In this paper, we present four explicit finite difference methods for the numerical solution of the EW and RLW equations. These methods are all based on the introduction of a new dependent variable and the separation of the space and time derivatives. The latter are discretized by means of first-order accurate, explicit differences, whereas spatial discretization is based on the use of three-point, fourth-order accurate, compact differences [20], piecewise local integration that results in locally exponential techniques [21,22], and standard second-order accurate discretizations. The methods presented here can also be used to obtain the numerical solution of the damped or undamped GRLW equation, as well as the other one-dimensional equations described above. The paper has been arranged as follows. In the next section, four explicit methods for the solution of the EW and RLW equations are presented and their linear stability is analyzed. In Section 3, an assessment of the explicit methods presented in the paper is reported in terms of the conservation of the first three invariants and the L2-norm errors of the numerical solution, as a function of the time step and grid spacing for initial conditions corresponding to the exact solutions of the EW and RLW equations. Since 3 also includes some results illustrating the effects of the parameters of the EW and RLW equations on the solution of these equations. A summary of the main findings of the paper is presented in the last section of the paper. 2. Numerical methods The explicit numerical methods considered in this paper are based on the separation of the time and space derivatives, as follows. Eq. (1) (with p = 2) can be written as the following system of partial differential equations: v  ut ; v  dvxx ¼ luxx  aux  uux ;

ð2Þ ð3Þ

624

J.I. Ramos / Applied Mathematics and Computation 179 (2006) 622–638

where Eq. (3) does not contain time derivatives and the third term of the right-hand side of this equation has been written in non-conservation form with the appropriate redefinition of . Eqs. (2) and (3) are to be solved in the interval (a, b) subject to u(t, a) = u(t, b) = 0, i.e., homogeneous Dirichlet boundary conditions in a finite interval, and u(0, x) = u0(x). A first-order accurate (Euler’s forward) method for Eq. (2) yields ¼ uni þ kvni ; unþ1 i

ð4Þ

where i and n denote the ith grid point, i.e., xi = a + (i  1)h, and the nth time level, i.e., tn = nk, respectively, k ba is the time step, n ¼ 0; 1; 2; . . . ; h ¼ NP1 is the grid spacing, and NP is the number of grid points, i = 1, 2, . . . , NP. In order to determine the value of vni that appears in Eq. (4), we have used the techniques presented in the following subsections. 2.1. Compact operator method By defining F  ux and G  uxx, Eq. (3) can be written as v  dvxx ¼ P ;

ð5Þ

where P  lG  bF, b  a + u, and then discretize Eq. (5) by means of a three-point, compact, fourth2 order accurate finite difference formula, i.e., ðvxx Þi ¼ h12 1þdd2 vvi=12 þ Oðh4 Þ, where d2vi = vi+1  2vi + vi1, which i yields 1 1 1 ðviþ1 þ 10vi þ vi1 Þ  2 d2 vi ¼ ðP iþ1 þ 10P i þ P i1 Þ; 12 12 h

ð6Þ

which provides a (diagonally dominant) tri-diagonal matrix for vi, with v1 = vNP = 0. Due to the dependence of F and G and, therefore, P on u, the following three-point, fourth-order accurate, compact difference equations were used to determine Fi and Gi, and, therefore, Pi: 1 d2 ui ðGiþ1 þ 10Gi þ Gi1 Þ ¼ 2 ; 12 h

ð7Þ

1 uiþ1  ui1 . ðF iþ1 þ 4F i þ F i1 Þ ¼ 6 2h

ð8Þ

and

Eqs. (7) and (8) can be solved to determine F ni and Gni subject to F1 = FNP = G1 = GNP = 0, and, therefore, ¼ lGni  bni F ni which can then be substituted into Eq. (6) to determine vni ; this value can then be substituted into Eq. (4) to calculate unþ1 . i The fact that boundary conditions on F and G are required is a consequence of the introduction of these two variables and the use of three-point, fourth-order accurate, compact discretizations for ux and uxx, while Eq. (1) only demands two boundary conditions for u. A linear Fourier-von Neumann stability analysis of Eqs. (4) and (6)–(8) with uni ¼ /n expðIiKhÞ; vni ¼ V n expðIiKhÞ; F ni ¼ fn expðIiKhÞ and Gni ¼ gn expðIiKhÞ, where I2 = 1 and K is the wave number, yields, after some lengthy algebra, the following amplification factor: P ni

/nþ1 ¼ R þ IS; /n

ð9Þ

where sin2 u  ;   1 sin2 u 1 þ 4D 3 sinð2uÞ   S ¼ CQ    ; 2 2   1 sin2 u 1  3 sin u 1 þ 4D 3 R ¼ 1  4F



ð10Þ ð11Þ

J.I. Ramos / Applied Mathematics and Computation 179 (2006) 622–638

625

 ¼ kb and D  ¼ d2 are the (diffusion) Fourier, convection (Courant) and diswhere Q ¼ 1  13 sin2 u; F ¼ kl ; C h h h2 Kh persion numbers, respectively, u ¼ 2 , and b is a constant which may be taken as the supremum of a + u. The amplification factor is here a complex number the modulus of which much be less than or equal to one  and D.  In this paper, we shall refer to the compact operator for linear stability and depends on u; F ; C method presented in this section as CE to indicate that it is a compact (C) explicit (E) technique. 2.2. Locally analytical method 1 This technique is based on the piecewise local integration of Eq. (5) as follows. Consider the interval xi1 6 x 6 xi+1 and assume that, in this interval, Eq. (5) may be approximated by the following linear one: v  dvxx ¼ P i ; which can be integrated analytically subject to v(xi1) = vi1, v(xi) = vi and v(xi+1) = vi+1 to yield       h h p ffiffi ffi p ffiffi ffi viþ1  2vi cosh þ vi1 ¼ 2P i cosh 1 ; d d

ð12Þ

ð13Þ

which results in a system of linear algebraic equations for vi with a tri-diagonal structure where v1 = vNP = 0. If, in Eq. (13), the first- and second-order spatial derivatives of u that appear in P are discretized by means of second-order accurate finite difference formulae, i.e. Pi ¼

l b ðuiþ1  2ui þ ui1 Þ  i ðuiþ1  ui1 Þ; 2h h2

ð14Þ

then, one could solve Eq. (13) (using uni ) to determine vni and this requires the inversion of a tri-diagonal matrix. Once vni is known, unþ1 can be easily determined from Eq. (4). Note that, since Eq. (12), is a second-order, i linear ordinary differential equation with constant coefficients that has (locally) exponential solutions, we shall refer to the method just described as EXE to indicate its exponential (EX) explicit (E) character. A linear Fourier-von Neumann stability analysis of Eqs. (4), (13) and (14) with uni ¼ /n expðIiKhÞ and n vi ¼ V n expðIiKhÞ, yields, after some lengthy algebra, an amplification factor of the same form as Eq. (9) but with Q R ¼ 1  4F sin2 u; ð15Þ Z  Q sinð2uÞ; ð16Þ S ¼ C Z pffiffiffiffi pffiffiffiffi   1 and Z ¼ cosh D   cosð2uÞ. The amplification factor is here a complex number the where Q ¼ cosh D  and D.  modulus of which much be less than or equal to one for linear stability and depends on u; F ; C 2.3. Locally analytical method 2 Instead of approximating Pi in Eq. (13) as in Eq. (14), one can use a compact operator approximation analogous to that presented above, as follows. Eqs. (7) and (8) are used to determine Gni and F ni , respectively, and these values are employed to determine P ni ¼ lGni  bni F ni which can be substituted into Eq. (13) to determine vni . Then, Eq. (4) can be employed to calculate uinþ1 . In this paper, we shall refer to the technique resulting from these approximations as EXCE to indicate its locally exponential (EX) character, the use of compact (C) operators to determine f, g and P, and its explicitness (E). A linear Fourier-von Neumann stability analysis of Eqs. (4), (7), (8) and (13) with uni ¼ /n expðIiKhÞ; vni ¼ V n expðIiKhÞ; F ni ¼ fn expðIiKhÞ and Gni ¼ gn expðIiKhÞ, yields, after some lengthy algebra, an amplification factor analogous to Eq. (9) but with Q sin2 u R ¼ 1  4F ; Z 1  13 sin2 u  Q sinð2uÞ ; S ¼ C Z 1  23 sin2 u

ð17Þ ð18Þ

626

J.I. Ramos / Applied Mathematics and Computation 179 (2006) 622–638

pffiffiffiffi pffiffiffiffi   1 and Z ¼ cosh D   cosð2uÞ. The amplification factor is here a complex number the where Q ¼ cosh D  and D.  modulus of which much be less than or equal to one for linear stability and depends on u; F ; C Since EXCE makes uses of three-point, fourth-order accurate, finite difference formulae for F and G in the determination of u and requires that F1 = FNP = G1 = GNP = 0, the method requires that u, ux and uxx be zero at the upstream and downstream boundaries. By way of contrast, Eq. (1) only requires two boundary conditions. 2.4. Explicit method In this technique, Eq. (3) is discretized by means of second-order accurate central difference formulae to yield vi 

d l b ðviþ1  2vi þ vi1 Þ ¼ 2 ðuiþ1  2ui þ ui1 Þ  i ðuiþ1  ui1 Þ; 2 2h h h

ð19Þ

subject to v1 = vNP = 0, which can be used to determine vni from uni by inverting a tri-diagonal matrix. vni can then be substituted into Eq. (4) to determine uinþ1 . In this paper, we shall refer to the explicit finite difference technique just described as E to indicate its explicit (E) character. A linear Fourier-von Neumann stability analysis of Eqs. (4) and (19) with uni ¼ /n expðIiKhÞ and n vi ¼ V n expðIiKhÞ, yields an amplification factor identical to Eq. (9) but with R ¼ 1  4F  S ¼ C

sin2 u ;  sin2 u 1 þ 4D

ð20Þ

sinð2uÞ .  sin2 u 1 þ 4D

ð21Þ

Since for linear stability the modulus of the amplification factor must be less than or equal to one, and it is  2 6 2F and F 6 1 ð1 þ DÞ,  which for strictly positive, it can be easily shown that these conditions require that C 2  the standard advection–diffusion equation that results upon setting d = 0 in Eq. (1), i.e., for D ¼ 0, reduce to  2 6 1. the well-known conditions F 6 12 and, therefore, C  It must be noted that, since C depends on b = a + u, one should use the supremum of b to determine the  2 6 2F . It should also be noted that the four explicit methods presented in this paper time step that satisfies C require the inversion of tri-diagonal matrices to determine vni and/or F ni and Gni .   1, a two-term Taylor’s series expansion of Eq. (13) yields E. Therefore, one It must be noted that, for D   1. should expect that, in the absence of round-off errors, EXE and E produce nearly identical results for D As pointed out by some authors, e.g., [11], numerical methods for equations that preserve some of the conserved quantities, i.e., invariants, of the original differential equations may exhibit a more accurate behavior in time than methods that do not preserve them. For l = 0, the EW and RLW equation have the following invariants: I1  I2  I3 

Z

þ1

1 Z þ1 1 Z þ1 1

uðt; xÞ dx ¼

Z

þ1

uð0; xÞ dx;

ð22Þ

1



 u2 ðt; xÞ þ du2x ðt; xÞ dx;

 a  3 u2 ðt; xÞ þ u3 ðt; xÞ dx; 

ð23Þ ð24Þ

provided that u(t, ±1) = ux(t, ±1) = uxx(t, ±1) = 0. If l 5 0, I1 is preserved, I3 is not conserved, and Z þ1 Z þ1  2  dI 2 d  u ðt; xÞ þ du2x ðt; xÞ dx ¼ 2l u2x ðt; xÞ dx; ð25Þ dt 1 dt 1 indicates that I2 decreases with time owing to dissipation or viscosity, i.e., l 5 0.

J.I. Ramos / Applied Mathematics and Computation 179 (2006) 622–638

627

3. Presentation of results In this section, we assess the accuracy of the four explicit methods described in the previous section as a function of the time step and grid spacing or number of grid points. Since as indicated previously, only the linear stability of E provides explicit stability conditions, these conditions are used as a guide for the time step and grid spacing selection in the other three techniques. In addition, it should be note that the explicit methods presented here can use variable time steps according to the evolution of the solution, i.e., according to the supremum of (a + ui) at each time step. 3.1. The EW equation As indicated previously, the EW equation depends on , l and d; the usual EW equation is characterized by a = l = 0 and has the following analytical exact solution: ue ðt; xÞ ¼

3c ; cosh ðcðx  xo  ctÞÞ

ð26Þ

2

where c ¼ 2p1 ffiffid, which corresponds to a solitary wave of amplitude A = 3c, speed c and width c located initially at x0. For a = 0, b = 100, x0 = 40, c = 0.3, l = a = 0,  = d = 1, and initial conditions given by the exact solution, calculations were performed with the finite difference methods described in the previous section for different time steps and grid spacings, and some sample results are presented in Tables 1–4. Table 1 shows that E preserves very accurately the first invariant, whereas the second and third invariants decrease slightly as the time step is decreased for NP = 401. The L2-norm errors increase slightly as the time step is decreased for NP = 401, presumably because of the increase in round-off errors. Table 1 also indicates that the second and third invariants slightly increase and decrease, respectively, as the number of grid points is increased whenever hk is kept constant; however, the L2-norm errors decrease substantially as NP is increased while keeping hk constant. Analogous results to those exhibited in Table 1 have been found for CE as indicated in Table 2, except that this method results in lower L2-norm errors than E owing to its fourth-order spatial accuracy, although CE requires longer computational times, performs more operations, and may be subject to larger round-off errors than E. In addition, while the L2-norm errors were found to increase as k was decreased for NP = 401, those of CE have been found to decrease with the decrease of the time step. The results presented in Table 3 were obtained with EXE and indicate that this method behaves in an analogous manner to E in regards to the decrease of the second and third invariants and the increase in the L2norm errors as the time step is decreased for NP = 401. However, the L2-norm errors are smaller than those of E for all the cases considered here except for that corresponding to NP = 201. As indicated in Table 4, EXCE is more accurate than EXE and E and its errors decrease as the time step or the grid spacing is decreased. 1 ð0Þ 2 ð0Þ 3 ð0Þ Figs. 1 and 2 show the normalized three invariants, i.e., RI 1 ¼ I 1 ðtÞI ; RI 2 ¼ I 2 ðtÞI ; RI 3 ¼ I 3 ðtÞI , and I 1 ð0Þ I 2 ð0Þ I 3 ð0Þ the L2-norm errors obtained with CE as functions of the time step and grid spacing, respectively. Fig. 1 shows that for the first four cases of Tables 1–4, the normalized first invariant is well-preserved regardless of the time Table 1 Values of the invariants I1, I2 and I3 and L2-norm errors at t = 40 as functions of the time step and grid spacing for the solution of the EW equation obtained with E (NP, k)

I1

I2

I3

L2-norm error

(401, 0.01000) (401, 0.00500) (401, 0.00250) (401, 0.00125) (201, 0.00125) (801, 0.0006250) (1601, 0.0003125)

3.600000 3.599999 3.599999 3.600000 3.600000 3.599999 3.600000

2.611728 2.598592 2.592086 2.588848 2.569942 2.592049 2.592429

1.578784 1.566885 1.560998 1.558070 1.557098 1.556687 1.555948

0.430810E02 0.488855E02 0.519549E02 0.535215E02 0.215065E01 0.130103E02 0.305161E03

628

J.I. Ramos / Applied Mathematics and Computation 179 (2006) 622–638

Table 2 Values of the invariants I1, I2 and I3 and L2-norm errors at t = 40 as functions of the time step and grid spacing for the solution of the EW equation obtained with CE (NP, k)

I1

I2

I3

L2-norm error

(401, 0.01000) (401, 0.00500) (401, 0.00250) (401, 0.00125) (201, 0.00125) (801, 0.0006250) (1601, 0.0003125)

3.599999 3.600000 3.600000 3.600001 3.600000 3.599999 3.600000

2.612544 2.599007 2.592304 2.588970 2.570723 2.592064 2.592432

1.579549 1.567285 1.561220 1.558205 1.558199 1.556701 1.555950

0.161169E02 0.795278E03 0.389948E03 0.188069E03 0.154111E03 0.999448E04 0.503415E04

Table 3 Values of the invariants I1, I2 and I3 and L2-norm errors at t = 40 as functions of the time step and grid spacing for the solution of the EW equation obtained with EXE (NP, k)

I1

I2

I3

L2-norm error

(401, 0.01000) (401, 0.00500) (401, 0.00250) (401, 0.00125) (201, 0.00125) (801, 0.0006250) (1601, 0.0003125)

3.599999 3.600000 3.600000 3.600000 3.600001 3.599998 3.599998

2.611949 2.598730 2.592185 2.588926 2.570921 2.592056 2.592430

1.579040 1.567065 1.561141 1.558195 1.558820 1.556696 1.555949

0.336729E02 0.393829E02 0.424675E02 0.440513E02 0.180214E01 0.106005E02 0.244791E03

Table 4 Values of the invariants I1, I2 and I3 and L2-norm errors at t = 40 as functions of the time step and grid spacing for the solution of the EW equation obtained with EXCE (NP, k)

I1

I2

I3

L2-norm error

(401, 0.01000) (401, 0.00500) (401, 0.00250) (401, 0.00125) (201, 0.00125) (801, 0.0006250) (1601, 0.0003125)

3.600001 3.599999 3.600000 3.600000 3.600000 3.599999 3.599999

2.612880 2.599115 2.592301 2.588911 2.569241 2.592063 2.592430

1.579960 1.567488 1.561321 1.558256 1.558527 1.556707 1.555951

0.339077E02 0.263893E02 0.228759E02 0.212119E02 0.789379E02 0.571546E03 0.164143E03

step; however, the second and third invariants and the L2-norm errors increase almost linearly with time, and increase as the time step is increased. Although not shown here, the first, second and third invariants obtained with E and EXE exhibit similar trends to those of Fig. 1, whereas the L2-norm errors of these methods increase almost linearly with time; they also increase slightly as the time step is decreased. On the other hand, the errors of EXCE have been found to decrease as the time step is decreased and seem to tend to asymptotic constant values for large times for cases 2–4 of Tables 1–4. Fig. 2 indicates that the first invariant is well preserved, i.e., it is nearly independent of the number of grid points, whereas the second and third invariants increase almost linearly with time at a rate that decreases as the time step is decreased, but they are nearly identical for cases 4 and 5 of Tables 1–4. The L2-norm errors increase in an almost quadratic fashion with time for cases 4, 6 and 7 of Tables 1–4, and decrease as the time step is decreased for these cases. However, for case 5, the error increases parabolically but tends to an almost constant asymptotic value for large times. Although not shown here, the second and third invariants obtained with E and EXE increase with time at a rate that decreases as the grid spacing is decreased, and the errors for case 5 are larger than those for cases 4, 6

J.I. Ramos / Applied Mathematics and Computation 179 (2006) 622–638 x 10

2

−7

0.012 0.01

0

0.008 RI

2

−2 RI1

629

−4

0.004

−6 −8 0

0.006

0.002 10

20

30

0

40

0

t

10

0.02

2

20

30

40

20

30

40

t

−3

x 10

EL error

0.01

1

2

RI

3

0.015

0.005

0

0

10

20 t

30

40

0 0

10 t

Fig. 1. Normalized first (top left), second (top right) and third (bottom left) invariants and L2-norm errors (bottom right) as functions of the time step for the solution of the EW equation obtained with CE: (NP, k) = (401, 0.01000), s; (401, 0.00500), ; (401, 0.00250), h; (401, 0.00125), v.

and 7 of Tables 1–4. The L2-norm errors of E also increase linearly with time. On the other hand, the errors of EXCE exhibit similar trends to those of E but increase in an almost parabolic fashion with time, and tend to asymptotic constant values for large times. Tables 5–8 show the first three invariants and the L2-norm errors of the four explicit methods considered in this paper as functions of the time step for a fixed number of grid points and as functions of the grid spacing for a fixed time step. Once gain, the four methods preserve very well the first invariant. The second and third invariants decrease slightly as the time step is decreased for NP = 401, whereas they increase slightly as the grid spacing is decreased for k = 0.1 for E, CE and EXCE. However, EXE predicts a third invariant that increases slightly as the grid spacing is decreased. The L2-norm errors shown in Tables 5–8 are larger than those of Tables 1–4 and decrease as the time step and the grid spacing are decreased. Tables 5–8 also show large errors in the second and third invariants, especially for NP = 401 and k = 0.2. Fig. 3 illustrates the effects of the parameters , l and d on the solutions of the EW equation for a = 0, NP = 801, h = 1/2, a = 0, b = 150, x0 = 40, k = 0.000625, initial conditions corresponding to the exact solution and the values of the parameters shown in Table 9. This figure indicates that the amplitude of the wave increases as c increases, while its width increases as d is increased. For l = 0, the solitary wave moves along a straight line, whereas the presence of viscosity/dissipation causes the wave to curve its trajectory as indicated in the bottom right graph of Fig. 3. In the presence of dissipation, the amplitude of the wave decreases whereas its width increases with time. 3.2. The RLW equation As indicated previously, the RLW equation depends on a, , l and d; the usual RLW equation is characterized by l = 0 and has the following analytical solution for a = 1: ue ðt; xÞ ¼

3c ; cosh ðcðx  xo  ce tÞÞ 2

ð27Þ

630

J.I. Ramos / Applied Mathematics and Computation 179 (2006) 622–638 −7

6

x 10

1.4

1 RI2

1

2 RI

0

−2

0.8 0.6 0.4

−4

0.2 0

10

20 t

30

0

40

−3

2

−3

1.2

4

−6

x 10

x 10

2

0 x 10

10

20 t

30

40

10

20 t

30

40

−4

RI

3

EL2 error

1.5

1

1

0.5

0

0

10

20 t

30

40

0

0

Fig. 2. Normalized first (top left), second (top right) and third (bottom left) invariants and L2-norm errors (bottom right) as functions of the grid spacing for the the solution of the EW equation obtained with CE: (NP, k) = (401, 0.00500), s; (201, 0.00125), ; (801, 0.0006250), h; (1601, 0.0003125), v.

Table 5 Values of the invariants I1, I2 and I3 and L2-norm errors at t = 40 as functions of the time step and grid spacing for the solution of the EW equation obtained with E (NP, k)

I1

I2

I3

L2-norm error

(401, 0.200) (401, 0.100) (401, 0.050) (401, 0.025) (201, 0.100) (801, 0.100) (1601, 0.100)

3.600000 3.599999 3.600000 3.600000 3.600000 3.600000 3.600000

3.270090 2.880173 2.723143 2.652160 2.830034 2.893544 2.896941

2.199281 1.826159 1.680490 1.615533 1.798870 1.833380 1.835213

0.326587E01 0.125228E01 0.453740E02 0.310718E02 0.132492E01 0.163510E01 0.173992E01

Table 6 Values of the invariants I1, I2 and I3 and L2-norm errors at t = 40 as functions of the time step and grid spacing for the solution of the EW equation obtained with CE (NP, k)

I1

I2

I3

L2-norm error

(401, 0.200) (401, 0.100) (401, 0.050) (401, 0.025) (201, 0.100) (801, 0.100) (1601, 0.100)

3.600000 3.599999 3.600000 3.600000 3.600001 3.600000 3.599999

3.298118 2.890414 2.727589 2.654247 2.867881 2.896155 2.897597

2.226784 1.835786 1.684602 1.617459 1.835094 1.835824 1.835826

0.393769E01 0.177183E01 0.842693E02 0.410724E02 0.174116E01 0.177460E01 0.177525E01

J.I. Ramos / Applied Mathematics and Computation 179 (2006) 622–638

631

Table 7 Values of the invariants I1, I2 and I3 and L2-norm errors at t = 40 as functions of the time step and grid spacing for the solution of the EW equation obtained with EXE (NP, k)

I1

I2

I3

L2-norm error

(401, 0.200) (401, 0.100) (401, 0.050) (401, 0.025) (201, 0.100) (801, 0.100) (1601, 0.100)

3.600000 3.599999 3.600000 3.599999 3.600000 3.599998 3.600000

2.882298 2.724104 2.652641 2.618623 2.837960 2.894085 2.897077

1.828214 1.681432 1.616028 1.585093 1.807345 1.833890 1.835341

0.133293E01 0.490421E02 0.247959E02 0.311449E02 0.105221E01 0.165895E01 0.174605E01

Table 8 Values of the invariants I1, I2 and I3 and L2-norm errors at t = 40 as functions of the time step and grid spacing for the solution of the EW equation obtained with EXCE I1

I2

I3

L2-norm error

(401, 0.200) (401, 0.100) (401, 0.050) (401, 0.025) (201, 0.100) (801, 0.100) (1601, 0.100)

3.600000 3.600000 3.600000 3.600000 3.599999 3.599999 3.599999

2.896153 2.729992 2.655304 2.619832 2.889894 2.897603 2.897961

1.841294 1.686929 1.618533 1.586267 1.858179 1.837184 1.836165

0.196983E01 0.101997E01 0.583430E02 0.378284E02 0.258587E01 0.182295E01 0.178727E01

4

10

2

5

u

u

(NP, k)

0 0

0 0 20

10 0

40

t

80 150

50

30

100 x

t

10

5

5

u

10 u

0

20

50

60

0 0

40 150

100 x

0 0 10

20 0

20 50

30 t

40 150

100 x

0

40 50

60 t

80 150

100 x

Fig. 3. Solution of the EW equation for the values of the parameters of Table 9: (case 1) top left; (case 2) top right; (case 3) bottom left; (case 4) bottom right.

pffiffiffiffiffiffi c , which corresponds to a solitary wave of amplitude A = 3c, speed c and width where ce = 1 + c and c ¼ 4dc e c located initially at x0. For a = 0, b = 100, x0 = 20, l = 0, c = 0.3, a =  = d = 1, and initial conditions given by the exact solution, calculations have been performed with the finite difference methods described in the previous section for dif-

632

J.I. Ramos / Applied Mathematics and Computation 179 (2006) 622–638

Table 9 Parameters of the EW equation Case

c

d

l

1 2 3 4

1 2 2 2

1 1 5 1

0.0 0.0 0.0 0.1

ferent time steps and grid spacings, and some sample results are presented in Tables 10–13. These tables indicate that the four explicit methods considered in this paper preserve very accurately the first invariant, whereas the second and third ones and the L2-norm errors decrease as k is decreased for NP = 401. The first invariant increases slightly and the second and third invariants decrease slightly as the number of grid points is increased when hk is kept constant. On the other hand, the L2-norm errors decrease substantially as the grid spacing is decreased while keeping kh constant. Tables 12 and 13 also show that the locally analytical methods EXE and EXCE are more accurate than E, except for the first two cases shown in these tables. For the RLW equation, EXE has been found to be more accurate (in terms of the L2-norm errors) than EXCE despite Table 10 Values of the invariants I1, I2 and I3 and L2-norm errors at t = 40 as functions of the time step and grid spacing for the solution of the RLW equation obtained with E (NP, k)

I1

I2

I3

L2-norm error

(401, 0.01000) (401, 0.00500) (401, 0.00250) (401, 0.00125) (201, 0.00125) (801, 0.0006250) (1601, 0.0003125)

7.503245 7.503246 7.503247 7.503246 7.503219 7.503259 7.503265

4.872303 4.785781 4.744000 4.723468 4.719965 4.713958 4.708991

17.387243 17.048864 16.885653 16.805496 16.798145 16.766430 16.746525

0.496251E02 0.263859E02 0.238826E02 0.261931E02 0.114336E01 0.613489E03 0.220363E03

Table 11 Values of the invariants I1, I2 and I3 and L2-norm errors at t = 40 as functions of the time step and grid spacing for the solution of the RLW equation obtained with CE (NP, k)

I1

I2

I3

L2-norm error

(401, 0.01000) (401, 0.00500) (401, 0.00250) (401, 0.00125) (201, 0.00125) (801, 0.0006250) (1601, 0.0003125)

7.503252 7.503248 7.503249 7.503244 7.503221 7.503261 7.503266

4.874431 4.786813 4.744528 4.723752 4.721616 4.713989 4.708992

17.395731 17.053043 16.887861 16.806753 16.806738 16.766556 16.746544

0.641792E02 0.315681E02 0.156986E02 0.791509E03 0.774218E03 0.415737E03 0.243777E03

Table 12 Values of the invariants I1, I2 and I3 and L2-norm errors at t = 40 as functions of the time step and grid spacing for the solution of the RLW equation obtained with EXE (NP, k)

I1

I2

I3

L2-norm error

(401, 0.01000) (401, 0.00500) (401, 0.00250) (401, 0.00125) (201, 0.00125) (801, 0.0006250) (1601, 0.0003125)

7.503247 7.503244 7.503246 7.503247 7.503219 7.503263 7.503270

4.873154 4.786196 4.744212 4.723582 4.720620 4.713973 4.708988

17.390678 17.050577 16.886585 16.806051 16.802126 16.766478 16.746532

0.535719E02 0.241640E02 0.147752E02 0.147934E02 0.686968E02 0.394810E03 0.210029E03

J.I. Ramos / Applied Mathematics and Computation 179 (2006) 622–638

633

Table 13 Values of the invariants I1, I2 and I3 and L2-norm errors at t = 40 as functions of the time step and grid spacing for the solution of the RLW equation obtained with EXCE (NP, k)

I1

I2

I3

L2-norm error

(401, 0.01000) (401, 0.00500) (401, 0.00250) (401, 0.00125) (201, 0.00125) (801, 0.0006250) (1601, 0.0003125)

7.503249 7.503248 7.503251 7.503249 7.503223 7.503263 7.503270

4.875628 4.787379 4.744805 4.723887 4.722147 4.714008 4.708995

17.400421 17.055262 16.888935 16.807283 16.808908 16.766624 16.746550

0.755209E02 0.432128E02 0.280965E02 0.211775E02 0.665990E02 0.711015E03 0.302923E03

the fact that the latter uses three-point, fourth-order accurate finite difference formulae to determine F, G and P, whereas the former employs second-order accurate finite differences for these terms. This may be due to the larger number of operations required by EXCE to obtain the solution. Tables 10, 12 and 13 also show that E yields much larger errors than EXE except for NP = 401 and k = 0.0125, and smaller errors than EXCE except for the fourth and fifth cases of the tables, i.e., except for NP = 401 and k = 0.0125 and NP = 201 and k = 0.00125, and this result is consistent with the analysis per  1 and this condition is not met for formed in Section 2 where it was shown that EXE coincides with E for D NP = 201. Tables 10–13 indicate that CE is more accurate (in terms of the L2-norm errors) than the other three techniques presented in this paper for the third, fourth and fifth cases of these tables. For cases 6 and 7, however, EXE is the most accurate method, whereas E is the most accurate technique for the first case. Figs. 4 and 5 show the normalized three invariants and the L2-norm errors obtained with CE as functions of the time step and grid spacing, respectively, for the cases reported in Tables 10–13. Fig. 4 shows that, for the

1.4

x 10

−3

0.04

1.2

2

0.8

RI

RI

1

1

0.6

0.02

0.4 0.2 0

0

10

20 t

30

0

40

0

10

20 t

30

40

10

20 t

30

40

−3

0.06

7

x 10

6 5

2

RI3

EL error

0.04

0.02

4 3 2 1

0

0

10

20 t

30

40

0

0

Fig. 4. Normalized first (top left), second (top right) and third (bottom left) invariants and L2-norm errors (bottom right) as functions of the time step for the solution of the RLW equation obtained with CE: (NP, k) = (401, 0.01000), s; (401, 0.00500), ; (401, 0.00250), h; (401, 0.00125), v.

634

J.I. Ramos / Applied Mathematics and Computation 179 (2006) 622–638

1.4

x 10

−3

−3

5

1.2

x 10

4

1 3 RI2

RI1

0.8 0.6

2

0.4 1

0.2 0

0

10

20 t

30

0

40

−3

5

x 10

8

4

x 10

10

20 t

30

40

10

20 t

30

40

−4

EL 2 error

6

RI3

3 2

4

2

1 0

0

0

10

20 t

30

40

0

0

Fig. 5. Normalized first (top left), second (top right) and third (bottom left) invariants and L2-norm errors (bottom right) as functions of the grid spacing for the the solution of the EW equation obtained with CE: (NP, k) = (401, 0.00500), s; (201, 0.00125), ; (801, 0.0006250), h; (1601, 0.0003125), v.

first four cases of Tables 10–13, the normalized first invariant increases linearly with time, but it is nearly independent of the time step; however, the second and third invariants and the L2-norm errors increase almost linearly with time, and increase as the time step is increased. Although not shown here, the first, second and third invariants obtained with E exhibit similar trends to those of Fig. 4 as functions of the time step, whereas the L2-norm errors of this method increase almost quadratically with time, but tend to the same asymptotic constant value for cases 2, 3 and 4 of Tables 10–13. On the other hand, the errors of EXE EXCE exhibit similar trends to those of CE and decrease as the time step is decreased. Fig. 5 indicates that the first invariant is nearly independent of the number of grid points but it increases with time, whereas the second and third invariants increase almost linearly with time at a rate that decreases as the time step is decreased, but they are nearly identical for cases 4 and 5 of Tables 10–13. The L2-norm errors also increase in an almost linear fashion with time for cases 6 and 7 of Tables 10–13, and decrease as the time step is decreased for these cases. However, for cases 4 and 5, the error is nearly the same. The first, second and third invariants predicted by E, EXE and EXCE exhibit similar trends to those shown in Fig. 5, but the L2-norm errors were found to increase linearly with time for cases 4, 6 and 7 of Tables 10–13, and decrease as the grid spacing is decreased. Case 4 resulted in the largest errors which exhibit a parabolic shape as a function of time and tend towards constant asymptotic values for large times. Tables 14–17 clearly indicate that the methods presented here preserve rather well the first invariant of the RLW equation. However, in spite that the last three cases of these tables show a third invariant in accord with that of Tables 10–13, the second invariant is predicted poorly. The fact that the second invariant of Tables 14– 17 is much larger than that of Tables 10–13 while the first invariant is preserved very well indicates that the reasons for the discrepancies between these two sets of tables is due to the term ux which appears in the second invariant (cf. Eq. (23)) and reflects that the numerical solutions corresponding to the values of NP and k reported in Tables 14–17 must exhibit oscillations.

J.I. Ramos / Applied Mathematics and Computation 179 (2006) 622–638

635

Table 14 Values of the invariants I1, I2 and I3 and L2-norm errors at t = 40 as functions of the time step and grid spacing for the solution of the RLW equation obtained with E (NP, k)

I1

I2

I3

L2-norm error

(401, 0.200) (401, 0.100) (401, 0.050) (401, 0.025) (201, 0.100) (801, 0.100) (1601, 0.100)

7.503259 7.503253 7.503252 7.503249 7.503226 7.503264 7.503273

20.017279 7.528773 5.740211 5.157949 7.348070 7.577495 7.589907

83.434036 28.031851 20.811701 18.508327 27.313272 28.224897 28.274067

0.335583E+00 0.878740E01 0.343900E01 0.147885E01 0.754583E01 0.913745E01 0.922864E01

Table 15 Values of the invariants I1, I2 and I3 and L2-norm errors at t = 40 as functions of the time step and grid spacing for the solution of the RLW equation obtained with CE (NP, k)

I1

I2

I3

L2-norm error

(401, 0.200) (401, 0.100) (401, 0.050) (401, 0.025) (201, 0.100) (801, 0.100) (1601, 0.100)

7.503259 7.503254 7.503251 7.503249 7.503227 7.503266 7.503275

20.961472 7.591734 5.756601 5.164124 7.584011 7.593503 7.593934

87.922615 28.290277 20.877106 18.532789 28.286205 28.290527 28.290535

0.353824E+00 0.924924E01 0.371744E01 0.169133E01 0.923192E01 0.925536E01 0.925828E01

Table 16 Values of the invariants I1, I2 and I3 and L2-norm errors at t = 40 as functions of the time step and grid spacing for the solution of the RLW equation obtained with EXE (NP, k)

I1

I2

I3

L2-norm error

(401, 0.200) (401, 0.100) (401, 0.050) (401, 0.025) (201, 0.100) (801, 0.100) (1601, 0.100)

7.503253 7.503249 7.503248 7.503245 7.503225 7.503267 7.503272

7.552972 5.746656 5.160400 4.918195 7.434186 7.583730 7.591487

28.131227 20.837481 18.518082 17.567059 27.669451 28.250456 28.280487

0.896014E01 0.354367E01 0.155602E01 0.695963E02 0.813055E01 0.918226E01 0.923995E01

Table 17 Values of the invariants I1, I2 and I3 and L2-norm errors at t = 40 as functions of the time step and grid spacing for the solution of the RLW equation obtained with EXCE (NP, k)

I1

I2

I3

L2-norm error

(401, 0.200) (401, 0.100) (401, 0.050) (401, 0.025) (201, 0.100) (801, 0.100) (1601, 0.100)

7.503253 7.503252 7.503252 7.503251 7.503227 7.503269 7.503274

7.628969 5.766107 5.167660 4.921364 7.737756 7.602737 7.596239

28.443010 20.914938 18.546722 17.579538 28.919468 28.328367 28.299967

0.949608E01 0.386701E01 0.181327E01 0.923009E02 0.102522E+00 0.931656E01 0.927354E01

J.I. Ramos / Applied Mathematics and Computation 179 (2006) 622–638 0.4

0.4

0.2

0.2

u

u

636

0 0

0 0 20

20 0

40

80 150

50

60

100 x

t

0.4

0.4

0.2

0.2

u

u

t

0

40

50

60

80 150

100 x

0 0

0 0 20

20 0

40 60 t

80 150

0

40

50

50

60

100 x

t

80 150

100 x

Fig. 6. Solution of the RLW equation for the values of the parameters of Table 18: (case 2) top left; (case 3) top right; (case 4) bottom left; (case 5) bottom right.

Table 18 Parameters of the RLW equation Case

c

d



a

l

1 2 3 4 5 6

0.03 0.10 0.10 0.10 0.10 0.10

1 1 2 1 1 1

1 1 1 1 2 1

0 0 0 0 0 1

0.0 0.0 0.0 0.1 0.0 0.0

The results presented in Tables 5–7, 8, 15–17 clearly indicate that, for the time steps and grid spacings considered in these tables, the explicit methods presented in this paper perform very well for the EW equation but rather poorly for the RLW equation, and this is a consequence of the linear advective term, i.e., a = 1, in the latter that increases the supremum of (a + u) and, therefore, the Courant number, and, therefore, requires a smaller time step for the RLW equation than for the EW one. Fig. 6 illustrates the effects of the parameters a, , l and d on the solutions of the RLW equation for NP = 801, h = 1/2, a = 0, b = 150, x0 = 40, k = 0.000625, initial conditions corresponding to the exact solution and the values of the parameters of Table 18. In case 1, it was observed that the solitary wave did not undergo any noticeable change as it propagated towards the downstream boundary, because the nonlinear terms are small owing to the small wave speed of this case. However, cases 2–5 yield the results presented in Fig. 6 which clearly show that the initial profile undergoes steepening and narrowing processes that result in the creation of secondary solitary waves in cases 2, 3 and 5. The creation of these secondary waves is accelerated by a decrease in d and an increase in  as indicated for cases 2, 3 and 5. Fig. 3 also shows that, even though the presence of the viscosity/dissipation decreases the solitary wave amplitude and increases the wave width, the nonlinear terms and those associated with the mixed third-order derivative try to create a secondary wave as illustrated in the third graph of Fig. 3 which corresponds to case 4. Although not shown here, case 6 exhibited a solitary wave that propagated at much faster speed owing to a and did not create a secondary wave.

J.I. Ramos / Applied Mathematics and Computation 179 (2006) 622–638

637

4. Conclusions Four explicit finite difference methods based on the introduction of a new dependent variable that results in the separation of the time and space derivatives for the numerical solution of the equal-width (EW) and regularized long-wave (RLW) equations are presented. The methods are based on the use of three-point, fourthorder accurate, compact operator techniques, analytical solutions of second-order, linear, ordinary differential equations with either second-order accurate finite differences or three-point, fourth-order accurate, compact operator approximations, and standard second-order accurate finite difference formulae. The linear stability of the explicit methods has been analyzed, especially for the standard explicit one where it is shown that the Courant number must be smaller than the square root of twice the Fourier number, whereas the latter must be smaller than half the sum of unity and the dispersion number. These restrictions are shown to result in the well-known Courant and Fourier number restrictions for the one-dimensional constant coefficients linear advection–diffusion equation. A rather extensive study of the accuracy of the methods presented here has been performed as a function of the time step and grid spacing for both the EW and RLW equation. From this assessment, it may be concluded that the most accurate technique developed in this paper is the one based on compact operator methods for both the first- and second-order spatial derivatives. The second most accurate technique is the one corresponding to an exponential method that employs compact operators for the first- and second-order spatial derivatives. It is also shown that exponential techniques which use second-order accurate approximations to the spatial derivatives are more accurate than the standard explicit method. Acknowledgments The research reported in this paper was supported by Project FIS2001-03191 from the Ministerio de Eduacio´n y Ciencia of Spain and fondos FEDER (30112005). References [1] T.B. Benjamin, J.L. Bona, J.J. Mahony, Model equations for long waves in nonlinear dispersive systems, Phil. Trans. Roy. Soc. (London), Series A 272 (1972) 47–78. [2] D.H. Peregrine, Calculations of the development of an undular bore, J. Fluid Mech. 25 (1966) 321–330. [3] D.H. Peregrine, Long waves on a beach, J. Fluid Mech. 27 (1967) 815–827. [4] R.K. Dodd, J.C. Eilbeck, J.D. Gibbon, H.C. Morris, Solitons and Nonlinear Wave Equations, Academic Press, New York, 1982. [5] J.C. Lewis, J.A. Tjon, Resonant production of solitons in the RLW equation, Phys. Lett. A 73 (1979) 275–279. [6] K.R. Raslan, Collocation method using quartic B-spline for the equal width (EW) equation, Appl. Math. Comput. 168 (2005) 795– 805. [7] L.R.T. Gardner, G.A. Gardner, A. Dogan, A least squares finite element scheme for the RLW equation, Commun. Numer. Meth. Eng. 12 (1996) 795–804. [8] L.R.T. Gardner, G.A. Gardner, Solitary waves of the regularised long-wave equation, J. Comput. Phys. 91 (1990) 441–459. [9] D.M. Sloan, Fourier pseudospectral solution of the regularised long wave equation, J. Comput. Appl. Math. 36 (1991) 159–179. [10] A. Arau´jo, A. Dura´n, Error propagation in the numerical integration of solitary waves. The regularized long wave equation, Appl. Numer. Math. 36 (2001) 197–217. [11] A. Dura´n, M.A. Lo´pez-Marcos, Conservative numerical methods for solitary wave interactions, J. Phys. A 36 (2003) 7761–7770. [12] J.L. Bona, A. Soyeur, On the stability of solitary wave solutions of model equations for long waves, J. Nonlinear Sci. 4 (1994) 449– 470. [13] Z. Luo, R. Liu, Mixed finite element analysis and numerical solution for the RLW equation, SIAM J. Numer. Anal. 36 (1998) 89–104. [14] P.C. Jain, R. Shankar, T.V. Singh, Numerical solution of the regularized long-wave equation, Commun. Numer. Meth. Eng. 9 (1993) 579–586. [15] S.I. Zaki, Solitary waves of the splitted RLW equation, Comput. Phys. Commun. 138 (2001) 80–91. [16] I. Dag˘, B. Saka, D. Irk, Galerkin method for the numerical solution of the RLW equation using quintic B-splines, J. Comput. Appl. Math., in press. [17] T.S. El-Danaf, M.A. Ramadan, F.E.I. Abd Alaal, The use of Adomian decomposition method for solving the regularized long-wave equation, Chaos Solitons Fract. 26 (2005) 747–757. [18] D. Kaya, S.M. El-Sayed, An application of the decomposition method for the generalized KdV and RLW equations, Chaos Solitons Fract. 17 (2003) 869–877. [19] B.K. Shivamoggi, D.K. Rollins, Evolution of solitary-wave solution of the perturbed regularized long-wave equation, Chaos Solitons Fract. 13 (2002) 1129–1136.

638

J.I. Ramos / Applied Mathematics and Computation 179 (2006) 622–638

[20] J.I. Ramos, Implicit, compact, linearized h-methods with factorization for multidimensional reaction–diffusion equations, Appl. Math. Comput. 94 (1998) 17–43. [21] J.I. Ramos, On diffusive methods and exponentially-fitted techniques, Appl. Math. Comput. 103 (1999) 69–96. [22] J.I. Ramos, A smooth locally-analytical technique for singularly perturbed two-point boundary value problems, Appl. Math. Comput. 163 (2005) 1123–1142.