- Email: [email protected]
Numerical investigation for the solitary waves interaction of the “good” Boussinesq equation
Applied Numerical Mathematics 45 (2003) 161–173 www.elsevier.com/locate/apnum
Numerical investigation for the solitary waves interaction of the “good” Boussinesq equation H. El-Zoheiry Department of Engineering Mathematics and Physics, Faculty of Engineering, Cairo University, Giza, Egypt
Abstract The “good” Boussinesq equation is studied numerically using an iterative implicit finite-difference scheme. The stability and accuracy of the proposed method are discussed. Soliton solutions are shown to exist for a finite range of amplitude size. The features of break-up and blow-up of solution are also witnessed. The reported results are in conformity with the available results. 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: Nonlinear; One-dimensional; Finite difference; Soliton; Collision
1. Introduction The importance of soliton-producing nonlinear wave equations is well understood among theoretical physicists and applied mathematicians. An example of a soliton producing equation is the Boussinesq equation (1) vt t = vxx + vxxxx + v 2 xx , which describes shallow water waves propagating in both directions. Eq. (1) which is known as the “bad” Boussinesq equation, has been studied by Hirota [5]. The related equation (2) vt t = vxx − vxxxx + v 2 xx , is the nonlinear string equation [13] or the “good” Boussinesq equation [9]. Closely related to Eq. (1) is the Korteweg–de Vries equation (KdV) (3) vt + vxxx + v 2 xx = 0. Eq. (2) is similar to the KdV equation in that it gives rise to solitons, however, it differs from this wellknown equation in a number of features, such as the finite range of velocities for the solitary waves and E-mail address: [email protected] (H. El-Zoheiry). 0168-9274/02/$30.00 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. doi:10.1016/S0168-9274(02)00187-3
162
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
the possibility of interactions leading to singular solutions. Manoranjan et al. [7] obtained a closed form solution for the two soliton interactions of Eq. (2) and carried out numerical experiments based on the Petrov–Galerkin method with linear “hat” functions to demonstrate the possibility of the break-up of an initial pulse into two solitons. In [8] it has been shown that the “good” Boussinesq equation possesses a highly complicated mechanism for the solitary waves interaction. The nonlinear stability and the convergence of some simple finite-difference schemes for the numerical solution of problems involving the “good” Boussinesq equation are presented in [10]. Iskandar and Jain [6] solved some equations analogous to Boussinesq equations. Recent studies related to Boussinesq equations are found in [12,1,4]. In this paper, the combined approach of linearization and finite-difference is used for the numerical integration of the “good” Boussinesq equation. The organization of the article is as follows: In Section 2, we present an unconditionally stable implicit scheme. In Section 3, the numerical procedure of the proposed method is discussed. Section 4 is devoted to some numerical experiments. 2. The finite-difference scheme Putting v = u − 1/2, it can easily be shown that the linear term vxx is removed from Eq. (2) resulting in the equation (4) ut t = −uxxxx + u2 xx . For the purpose of numerical analysis, it is simpler to use (4) instead of (2). Let u(n) denotes the value of u at the nth iteration and u(0) be the initial guess. Then the Picard’s approximation [3] for Eq. (4) is 2 (n) = −u(n+1) n = 0, 1, 2, . . . (5) u(n+1) tt xxxx + u xx , in the region R = [xmin x xmax ]x[t 0], subject to the initial conditions u(x, 0) = f (x),
ut (x, 0) = g(x),
(6)
and the boundary conditions u(0, t) = β1 (t),
u(Nh, t) = β2 (t),
ux (0, t) = γ1 (t),
ux (Nh, t) = γ2 (t).
(7)
The functions in the sequence {u(n) } satisfy the boundary conditions specified for u. Eq. (5) along with the given initial and boundary conditions is solved using the finite-difference method. A finite interval for the variable x is considered, the domain in the x − t plane is discretized by a grid with step length h and time step k. The exact value of u at the grid point (x, t) = (rh, sk) is denoted by ur,s (r = 0, 1, 2, . . . , N and s = 0, 1, 2, . . .), while the numerical value is designated by Ur,s , we take 1 2 (8) ut t |r,s = 2 δur,s − 2 ur,s + O k 2 , k k 2 1 1 (9) u xx |r,s = 2 δu2r+1,s + δu2r−1,s − 2 δu2r,s + O h2 + k 2 , 2h h 1 2 3 (10) uxxxx |r,s = 4 (δur+2,s + δur−2,s ) − 4 (δur+1,s + δur−1,s ) + 4 δur,s + O h2 + k 2 , 2h h h where δur,s = ur,s+1 + ur,s−1 .
(11)
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
163
The derivatives in the boundary conditions (7) are replaced by 1 ∂u = (Ur+1,s − Ur−1,s ), r = 0, N and s = 1, 2, . . . . ∂x r,s 2h
(12)
Substituting Eqs. (8)–(10) in Eq. (5) and neglecting the truncation error, the totality of equations can be expressed in the matrix form Aw = −Av + P By + P Bz + 2u + e
(13)
where
W1 W2 w = .. , . WN−1
V1 V2 v = .. , . VN−1
U1 U2 u = .. , . UN−1
Y12
Y2 2 y = . , .. 2 YN−1
V12
V2 2 z = . . .. 2 VN−1 (14)
The vector e has (N − 1) components involving the boundary conditions and A = I + QD where D and B are matrices of order (N − 1) given by 7 −4 1 −4 6 −4 1 −2 1 1 −4 6 −4 1 1 −2 1 ... ... ... ... , (15) D= , B = . . . . . . ... ... 1 −4 6 −4 1 1 −2 1 1 −4 6 −4 1 −2 1 −4 7 and I is the unit matrix of order (N − 1). The algebraic system of equations of (13) takes the form (1 + 6Q)Wr + Q(Wr+2 − 4Wr+1 − 4Wr−1 + Wr−2 ) 2 2 − 2Yr2 + Yr−1 = 2Ur − (1 + 6Q)Vr − Q(Vr+2 − 4Vr+1 − 4Vr−1 + Vr−2 ) + P Yr+1 2 2 − 2Vr2 + Vr−1 ), + P (Vr+1
r = 1, 2, . . . , N, s = 1, 2, . . . , n = 0, 1, 2, . . .
(16)
where (n+1) = Wr , Ur,s+1
(n) Ur,s+1 = Yr ,
(ns ) Ur,s = Ur ,
(n
)
s−1 Ur,s−1 = Vr ,
P=
k2 , 2h2
Q=
k2 , 4h4 (17)
and the superscript ns denotes the final number of iterations required to obtain an acceptable approximation to the value Ur,s at the grid points on the line t = sk subject to the criterion (18) max U (n+1) − U (n) 10−6 , 1 r N. r
r,s
r,s
164
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
A similar definition is valid for the superscript ns−1 . It should be noted that, scheme (16) represents a linear set of equations in the unknown values Wr . To study the stability of scheme (16), we apply the Fourier stability method [2]. We drop the iteration indices and linearize (16) by replacing the nonlinear terms Yr2 and Vr2 in the right-hand side of (16) by CYr and CVr with C as a constant such that
(19) max |Ur,s+1 |, |Ur,s |, |Ur,s−1| < C. r
According to Fourier’s stability method, let a solution of the linearized form of (16) be √ −1 mrh
Ur,s = ψ s e
,
(20)
where ψ = eαk , m is any real number and α = α(m), is in general, complex. Setting (20) into the linearized form of (16) results in k2 k2 −1 2 (21) ψ + ψ = 2 1 + C 2 (1 − cos mh) + 4 (1 − cos mh) , h h which can be expressed as the quadratic ψ 2 − 2µψ + 1 = 0,
(22)
in ψ, where k2 k2 1 = 1 + C 2 (1 − cos mh) + 4 (1 − cos mh)2 . (23) µ h h If we are to avoid an increasing exponential solution as s → ∞, it is necessary that |ψ| 1 for all values of m (the Von-Neumann’s necessary condition for stability [11]). From (23) it is clear that |µ| 1 for all values of m and for arbitrary k/ h. This means that the discriminant of (22) is non-positive. Hence the roots of (22) have modulus unity, thus the linearized form of the finite-difference scheme (16) is unconditionally stable. Expanding all terms of scheme (16) about (r, s) using Taylor’s method, the principal part of the truncation error was found to be k 2 ∂ 4 u h2 ∂ 6 u h2 ∂ 4 (u2 ) − + , (24) 12 ∂t 4 6 ∂x 6 12 ∂x 4 which is of O(h2 + k 2 ). To define (24), the exact solution ur,s of Eq. (5) subject to the given initial and boundary conditions is assumed to exist and to have six continuous derivatives with respect to x and four continuous derivatives with respect to t; that is, u ∈ C 6,4 . 3. Numerical procedure The difference scheme (16) can be rewritten as QWr+2 − 4QWr+1 + (1 + 6Q)Wr − 4QWr−1 + QWr−2 = F, r = 1, 2, . . . , N − 1, s = 1, 2, . . . , n = 0, 1, 2, . . .
(25)
where F = 2Ur − (1 + 6Q)Vr + Q(−Vr+2 + 4Vr+1 + 4Vr−1 − Vr−2 ) 2 2 2 2 − 2Yr2 + Yr−1 − 2Vr2 + Vr−1 + P Vr+1 . + P Yr+1
(26)
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
165
We use the second condition of (6) to find values on the line t = −k by employing a false boundary and the second-order central difference formula ur,1 − ur,−1 (27) + O k2 . ut |r,0 = 2k Writing g(rh) = gr , we have the approximation Ur,−1 = Ur,1 − 2kgr .
(28)
From (28), we replace Ur,−1 (≡ Vr ) by (Wr − 2kgr ) in the linear terms of (25) and by (Yr − 2kgr ) in the nonlinear terms of (25) to obtain QWr+2 − 4QWr+1 + (1 + 6Q)Wr − 4QWr−1 + QWr−2 = G, r = 1, 2, . . . , N − 1, s = 1, n = 0, 1, 2, . . .
(29)
where
2 2 − 2Yr2 + Yr−1 G = Ur + kQ(gr+2 − 4gr+1 − 4gr−1 + gr−2 ) + kgr (1 + 6Q) + P Yr+1 2 2 − 2gr2 + gr−1 − 2kP (Yr+1 gr+1 − 2Yr gr + Yr−1 gr−1 ). + 2k 2 P gr+1
(30)
Let the solution of (29) takes the recurrence relation Wr = Ar Wr+1 + Br Wr+2 + Cr ,
1 r N − 1.
(31)
We use (31) to find Wr−1 and Wr−2 and substitute them in (29) and compare the coefficients of the resulting equation with (31) we get E = 1 + 6Q − 4QAr−1 + QBr−2 + QAr−2 Ar−1 , Br = −Q/E, Ar = Q(4 + 4Br−1 − Ar−2 Br−1 )/E, Cr = (G + 4QCr−1 − QCr−2 − QAr−2 Cr−1 )/E, r = 2, 3, . . . , N − 1.
(32)
In view of the first boundary condition of (7), we have C0 = W0 ,
A0 = 0,
B0 = 0.
(33)
The third boundary condition can be discretized as W1 − W−1 = 2hγ1 (k),
(34)
with r = 1 and using (34), Eq. (31) will be identical with Eq. (29) if 4Q −Q (Gr=1 ) + 4QW0 + 2hγ1 (k) , B1 = , C1 = . (35) A1 = 1 + 7Q 1 + 7Q 1 + 7Q The scalars Ar , Br and Cr are computed from (32) in a forward sweep manner subject to the starting (n+1) is obtained values in (33) and (35). Using the stored values of Ar , Br and Cr , the solution Wr = Ur,s+1 by a backward sweep manner subject to the initial values WN (known) and WN−1 which can be evaluated using (31) with r = N − 1 and the discretization of the fourth boundary condition in (7). Scheme (16) together with the given initial and boundary conditions, is solved using the following iterative procedure: (I) For the first time level (s = 1) (a) Take Ur = f (rh), Yr = Ur for all r. (b) Solve (29) for Wr using the recurrence relation (31).
166
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
(c) Set Yr = Wr for all r. (d) Calculate improved values of Wr from (29). (e) Repeat steps (b)–(d) till condition (18) is satisfied. (II) For the time levels {s +1 | s = 1, 2, . . .} we repeat steps (b)–(d) except that in step (b), we solve (25) for Wr using (31) subject to the known initial values Ur and Vr for all r and the known boundary values U0,s and UN,s for all s. Scheme (16) is a penta-diagonal system of linear equations having constant matrix of coefficients. This means that only the right-hand of this scheme changes at each time level and this saves computational time and minimizes machine storage.
4. Numerical examples and results 4.1. Example 1 The analytic expression of solutions of Eq. (4) is [8] σ (x ± αt) − x0 + 0.5, u(x, t) = −ω sec h2 2
(36)
where x0 and σ > 0 are real parameters and the amplitude ω and velocity α of the wave are related to σ through the relation 3 α = 1 − σ 2. (37) ω = σ 2, 2 It should be noted that due to the square root in (37), the solitary waves (36) only exist for a finite range of velocities. Using the numerical schemes presented in the previous section, we now solve Eq. (4) together with the initial data 2 σ x + 0.5, g(x) = 0, (38) f (x) = −ω sec h 2 and zero boundary conditions. This is carried out for different values of ω and variety of mesh lengths h and k. As time evolves the initial pulse either splits into two pulses moving in opposite directions or blows-up according to the value of ω (see [7] for more details). In Fig. 1, an initial stationary pulse of amplitude 1.48 breaking-up into two pulses moving in opposite directions is shown. Fig. 2 displays the “blow-up” of an initial stationary pulse of amplitude 1.5. As time evolves, the amplitude of the pulse increases while its breadth decreases. 4.2. Example 2 In this example the interaction of two solitary waves travelling towards each other and colliding headon is illustrated. Consider Eq. (4) with the following initial conditions f (x) = −ω1 sec h2 Ω1 − ω2 sec h2 Ω2 , g(x) = −ω1 σ1 sec h2 Ω1 tan hΩ1 − ω2 σ2 sec h2 Ω2 tan hΩ2 ,
(39)
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
167
Fig. 1. Break-up of the pulse, h = 0.5, k = 0.05, −75 x 75, 0 t 60, x0 = 0, ω = 1.48.
Fig. 2. Blow-up of the pulse. h = 0.5, k = 0.1, −75 x 75, 0 t 10, x0 = 0, ω = 1.5.
where
3 2 α1 αi = 1 − σi2 , Ω1 = x + x0 , ωi = σi , 2 2 the boundary conditions taken at x = xmin are β1 (x) = −ω1 sec h2 ξ1 − ω2 sec h2 ξ2 + 1, γ1 (x) = −ω1 σ1 sec h2 ξ1 tan hξ1 − ω2 σ2 sec h2 ξ2 tan hξ2 ,
Ω2 =
α2 x − x0 , 2
i = 1, 2,
(40)
(41)
168
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
Fig. 3. Collision of two pulses. h = 0.5, k = 0.1, −75 x 75, 0 t 60, x0 = 5, ω1 = 0.25, ω2 = 0.25.
Fig. 4. Collision of two pulses. h = 0.5, k = 0.1, −75 x 75, 0 t 60, x0 = 4, ω1 = 0.35, ω2 = 0.35.
where σ1 σ2 (xmin − α1 t) + x0 , ξ2 = (xmin − α2 t) − x0 . (42) 2 2 On replacing xmin by xmax in (42) we get the boundary conditions β2 (t) and γ2 (t) taken at x = xmax . For the computational work suitable values of σ1 , σ2 , x0 and the number of points N (x-interval) were chosen so that the accuracy within reasonable run time is maintained. With the help of the computer plots, the ξ1 =
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
169
Fig. 5. Collision of two pulses. h = 0.5, k = 0.1, −75 x 75, 0 t 60, x0 = 5, ω1 = 0.35, ω2 = 0.15.
Fig. 6. Collision of two pulses. h = 0.5, k = 0.1, −75 x 75, 0 t 60, x0 = 10, ω1 = 0.38, ω2 = 0.38.
nature of the “good” Boussinesq solitary waves interaction is investigated for different amplitudes of these waves. Figs. 3–6 show typical interactions of the “good” Boussinesq two solitary waves plotted at successive time steps for different amplitude values. In Figs. 3–5, the pulses emerge with their original shape and speed after interaction as if the collision had not taken place. In Fig. 6, “blow-up” of the solution is obtained. Figs. 7 and 8 display the head on collision at selected time steps. It has been shown by Manoranjan et al. [7,8] that the nature of the interaction between any pair of solitary waves is dependent
170
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
Fig. 7. Head-on collision at selected time levels. h = 0.5, k = 0.1, −75 x 75, x0 = 5, ω1 = 0.35, ω2 = 0.35.
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
Fig. 8. Head-on collision at selected time levels. h = 0.5, k = 0.05, −75 x 75, x0 = 10, ω1 = 0.38, ω2 = 0.38.
171
172
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
Fig. 9. Joint amplitude (ωj ) for ω = 0.05(0.0025)0.35, h = 0.5, k = 0.05, x0 = 4, −75 x 75.
on their amplitudes ω1 and ω2 . They reported that, for real solitons interaction the amplitude values ω1 , and ω2 must satisfy certain conditions and when the initial amplitudes of the waves are equal (ω, say) “blow-up” of the solution is witnessed for ω greater than some critical value. It should be noted that when the two waves overlap (Figs. 3 and 4) then the joint amplitude (ωj ) is greater than twice the amplitude of an individual solitary wave (ω). Fig. 9 displays the relation between ωj and ω. The numerical computations presented in this paper were carried out on a personal computer and three to six iterations at each time level were required to satisfy condition (18).
Conclusions On the basis of linearization and finite-difference techniques, an iterative implicit scheme for the “good” Boussinesq equation is considered. Computational time is saved by using this method since a constant coefficient penta-diagonal system of equations is to be solved at each time level. Based on the numerical results and computer plots, the “good” Boussinesq equation has been found to possess an interesting soliton-interaction mechanism. For moderate amplitude values, the incoming waves emerge of the interaction without changing shape or velocity (clean interaction). However, if the amplitudes are large, the interaction leads to the creation of “blow-up” solutions. The numerical method developed in this paper provides results which are in conformity with the general ideas reported in [7] and [8].
References [1] Y. Agnon, P.A. Madsen, H.A. Schoffer, A new approach to high-order Boussinesq models, J. Fluid Mech. 399 (1999) 319–333. [2] W.F. Ames, Numerical Methods for Partial Differential Equations, 2nd Edition, Academic Press, New York, 1965.
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
173
[3] R.E. Bellman, R.E. Kalaba, Quasilinearization and Nonlinear Boundary Value Problems, Elsevier, New York, 1965. [4] A.G. Bratsos, The solution of the Boussinesq equation using the method of lines, Comput. Methods Appl. Mech. Engrg. 157 (1998) 33–44. [5] R. Hirota, Exact N -soliton solutions of the wave equation of long waves in shallow-water and in nonlinear lattices, J. Math. Phys. 14 (1973) 810–814. [6] L. Iskandar, P.C. Jain, Numerical solution of equations having inelastic soliton wave interaction, Comput. Math. Appl. 6 (1980) 373–383. [7] V.S. Manoranjan, A.R. Mitchell, J.L. Morris, Numerical solutions of the “good” Boussinesq equation, SIAM J. Sci. Statist. Comput. 5 (1984) 946–957. [8] V.S. Manoranjan, T. Ortega, J.M. Sang-Serna, Soliton and antisoliton interactions in the “good” Boussinesq equation, J. Math. Phys. 29 (1988) 1964–1968. [9] H.P. Mckean, Boussinesq’s equation on the circle, Comm. Pure Appl. Math. 34 (1981) 599–691. [10] T. Ortega, J.M. Sanz-Sema, Nonlinear stability and convergence of finite-difference methods for the “good” Boussinesq equation, Numer. Math. 58 (1990) 215–229. [11] R.D. Richtmyer, K.W. Morton, Difference Methods for Initial Value Problems, 2nd Edition, Wiley–Interscience, New York, 1967. [12] G. Schneider, The long wave limit for a Boussinesq equation, SIAM J. Appl. Math. 58 (1992) 1237–1245. [13] V.B. Zakharov, Randomization problem for one-dimensional chains of nonlinear operators, Zh. Eksper. Theoret. Fiz. 65 (1973) 219–228 (in Russian).
Numerical investigation for the solitary waves interaction of the “good” Boussinesq equation H. El-Zoheiry Department of Engineering Mathematics and Physics, Faculty of Engineering, Cairo University, Giza, Egypt
Abstract The “good” Boussinesq equation is studied numerically using an iterative implicit finite-difference scheme. The stability and accuracy of the proposed method are discussed. Soliton solutions are shown to exist for a finite range of amplitude size. The features of break-up and blow-up of solution are also witnessed. The reported results are in conformity with the available results. 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: Nonlinear; One-dimensional; Finite difference; Soliton; Collision
1. Introduction The importance of soliton-producing nonlinear wave equations is well understood among theoretical physicists and applied mathematicians. An example of a soliton producing equation is the Boussinesq equation (1) vt t = vxx + vxxxx + v 2 xx , which describes shallow water waves propagating in both directions. Eq. (1) which is known as the “bad” Boussinesq equation, has been studied by Hirota [5]. The related equation (2) vt t = vxx − vxxxx + v 2 xx , is the nonlinear string equation [13] or the “good” Boussinesq equation [9]. Closely related to Eq. (1) is the Korteweg–de Vries equation (KdV) (3) vt + vxxx + v 2 xx = 0. Eq. (2) is similar to the KdV equation in that it gives rise to solitons, however, it differs from this wellknown equation in a number of features, such as the finite range of velocities for the solitary waves and E-mail address: [email protected] (H. El-Zoheiry). 0168-9274/02/$30.00 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. doi:10.1016/S0168-9274(02)00187-3
162
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
the possibility of interactions leading to singular solutions. Manoranjan et al. [7] obtained a closed form solution for the two soliton interactions of Eq. (2) and carried out numerical experiments based on the Petrov–Galerkin method with linear “hat” functions to demonstrate the possibility of the break-up of an initial pulse into two solitons. In [8] it has been shown that the “good” Boussinesq equation possesses a highly complicated mechanism for the solitary waves interaction. The nonlinear stability and the convergence of some simple finite-difference schemes for the numerical solution of problems involving the “good” Boussinesq equation are presented in [10]. Iskandar and Jain [6] solved some equations analogous to Boussinesq equations. Recent studies related to Boussinesq equations are found in [12,1,4]. In this paper, the combined approach of linearization and finite-difference is used for the numerical integration of the “good” Boussinesq equation. The organization of the article is as follows: In Section 2, we present an unconditionally stable implicit scheme. In Section 3, the numerical procedure of the proposed method is discussed. Section 4 is devoted to some numerical experiments. 2. The finite-difference scheme Putting v = u − 1/2, it can easily be shown that the linear term vxx is removed from Eq. (2) resulting in the equation (4) ut t = −uxxxx + u2 xx . For the purpose of numerical analysis, it is simpler to use (4) instead of (2). Let u(n) denotes the value of u at the nth iteration and u(0) be the initial guess. Then the Picard’s approximation [3] for Eq. (4) is 2 (n) = −u(n+1) n = 0, 1, 2, . . . (5) u(n+1) tt xxxx + u xx , in the region R = [xmin x xmax ]x[t 0], subject to the initial conditions u(x, 0) = f (x),
ut (x, 0) = g(x),
(6)
and the boundary conditions u(0, t) = β1 (t),
u(Nh, t) = β2 (t),
ux (0, t) = γ1 (t),
ux (Nh, t) = γ2 (t).
(7)
The functions in the sequence {u(n) } satisfy the boundary conditions specified for u. Eq. (5) along with the given initial and boundary conditions is solved using the finite-difference method. A finite interval for the variable x is considered, the domain in the x − t plane is discretized by a grid with step length h and time step k. The exact value of u at the grid point (x, t) = (rh, sk) is denoted by ur,s (r = 0, 1, 2, . . . , N and s = 0, 1, 2, . . .), while the numerical value is designated by Ur,s , we take 1 2 (8) ut t |r,s = 2 δur,s − 2 ur,s + O k 2 , k k 2 1 1 (9) u xx |r,s = 2 δu2r+1,s + δu2r−1,s − 2 δu2r,s + O h2 + k 2 , 2h h 1 2 3 (10) uxxxx |r,s = 4 (δur+2,s + δur−2,s ) − 4 (δur+1,s + δur−1,s ) + 4 δur,s + O h2 + k 2 , 2h h h where δur,s = ur,s+1 + ur,s−1 .
(11)
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
163
The derivatives in the boundary conditions (7) are replaced by 1 ∂u = (Ur+1,s − Ur−1,s ), r = 0, N and s = 1, 2, . . . . ∂x r,s 2h
(12)
Substituting Eqs. (8)–(10) in Eq. (5) and neglecting the truncation error, the totality of equations can be expressed in the matrix form Aw = −Av + P By + P Bz + 2u + e
(13)
where
W1 W2 w = .. , . WN−1
V1 V2 v = .. , . VN−1
U1 U2 u = .. , . UN−1
Y12
Y2 2 y = . , .. 2 YN−1
V12
V2 2 z = . . .. 2 VN−1 (14)
The vector e has (N − 1) components involving the boundary conditions and A = I + QD where D and B are matrices of order (N − 1) given by 7 −4 1 −4 6 −4 1 −2 1 1 −4 6 −4 1 1 −2 1 ... ... ... ... , (15) D= , B = . . . . . . ... ... 1 −4 6 −4 1 1 −2 1 1 −4 6 −4 1 −2 1 −4 7 and I is the unit matrix of order (N − 1). The algebraic system of equations of (13) takes the form (1 + 6Q)Wr + Q(Wr+2 − 4Wr+1 − 4Wr−1 + Wr−2 ) 2 2 − 2Yr2 + Yr−1 = 2Ur − (1 + 6Q)Vr − Q(Vr+2 − 4Vr+1 − 4Vr−1 + Vr−2 ) + P Yr+1 2 2 − 2Vr2 + Vr−1 ), + P (Vr+1
r = 1, 2, . . . , N, s = 1, 2, . . . , n = 0, 1, 2, . . .
(16)
where (n+1) = Wr , Ur,s+1
(n) Ur,s+1 = Yr ,
(ns ) Ur,s = Ur ,
(n
)
s−1 Ur,s−1 = Vr ,
P=
k2 , 2h2
Q=
k2 , 4h4 (17)
and the superscript ns denotes the final number of iterations required to obtain an acceptable approximation to the value Ur,s at the grid points on the line t = sk subject to the criterion (18) max U (n+1) − U (n) 10−6 , 1 r N. r
r,s
r,s
164
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
A similar definition is valid for the superscript ns−1 . It should be noted that, scheme (16) represents a linear set of equations in the unknown values Wr . To study the stability of scheme (16), we apply the Fourier stability method [2]. We drop the iteration indices and linearize (16) by replacing the nonlinear terms Yr2 and Vr2 in the right-hand side of (16) by CYr and CVr with C as a constant such that
(19) max |Ur,s+1 |, |Ur,s |, |Ur,s−1| < C. r
According to Fourier’s stability method, let a solution of the linearized form of (16) be √ −1 mrh
Ur,s = ψ s e
,
(20)
where ψ = eαk , m is any real number and α = α(m), is in general, complex. Setting (20) into the linearized form of (16) results in k2 k2 −1 2 (21) ψ + ψ = 2 1 + C 2 (1 − cos mh) + 4 (1 − cos mh) , h h which can be expressed as the quadratic ψ 2 − 2µψ + 1 = 0,
(22)
in ψ, where k2 k2 1 = 1 + C 2 (1 − cos mh) + 4 (1 − cos mh)2 . (23) µ h h If we are to avoid an increasing exponential solution as s → ∞, it is necessary that |ψ| 1 for all values of m (the Von-Neumann’s necessary condition for stability [11]). From (23) it is clear that |µ| 1 for all values of m and for arbitrary k/ h. This means that the discriminant of (22) is non-positive. Hence the roots of (22) have modulus unity, thus the linearized form of the finite-difference scheme (16) is unconditionally stable. Expanding all terms of scheme (16) about (r, s) using Taylor’s method, the principal part of the truncation error was found to be k 2 ∂ 4 u h2 ∂ 6 u h2 ∂ 4 (u2 ) − + , (24) 12 ∂t 4 6 ∂x 6 12 ∂x 4 which is of O(h2 + k 2 ). To define (24), the exact solution ur,s of Eq. (5) subject to the given initial and boundary conditions is assumed to exist and to have six continuous derivatives with respect to x and four continuous derivatives with respect to t; that is, u ∈ C 6,4 . 3. Numerical procedure The difference scheme (16) can be rewritten as QWr+2 − 4QWr+1 + (1 + 6Q)Wr − 4QWr−1 + QWr−2 = F, r = 1, 2, . . . , N − 1, s = 1, 2, . . . , n = 0, 1, 2, . . .
(25)
where F = 2Ur − (1 + 6Q)Vr + Q(−Vr+2 + 4Vr+1 + 4Vr−1 − Vr−2 ) 2 2 2 2 − 2Yr2 + Yr−1 − 2Vr2 + Vr−1 + P Vr+1 . + P Yr+1
(26)
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
165
We use the second condition of (6) to find values on the line t = −k by employing a false boundary and the second-order central difference formula ur,1 − ur,−1 (27) + O k2 . ut |r,0 = 2k Writing g(rh) = gr , we have the approximation Ur,−1 = Ur,1 − 2kgr .
(28)
From (28), we replace Ur,−1 (≡ Vr ) by (Wr − 2kgr ) in the linear terms of (25) and by (Yr − 2kgr ) in the nonlinear terms of (25) to obtain QWr+2 − 4QWr+1 + (1 + 6Q)Wr − 4QWr−1 + QWr−2 = G, r = 1, 2, . . . , N − 1, s = 1, n = 0, 1, 2, . . .
(29)
where
2 2 − 2Yr2 + Yr−1 G = Ur + kQ(gr+2 − 4gr+1 − 4gr−1 + gr−2 ) + kgr (1 + 6Q) + P Yr+1 2 2 − 2gr2 + gr−1 − 2kP (Yr+1 gr+1 − 2Yr gr + Yr−1 gr−1 ). + 2k 2 P gr+1
(30)
Let the solution of (29) takes the recurrence relation Wr = Ar Wr+1 + Br Wr+2 + Cr ,
1 r N − 1.
(31)
We use (31) to find Wr−1 and Wr−2 and substitute them in (29) and compare the coefficients of the resulting equation with (31) we get E = 1 + 6Q − 4QAr−1 + QBr−2 + QAr−2 Ar−1 , Br = −Q/E, Ar = Q(4 + 4Br−1 − Ar−2 Br−1 )/E, Cr = (G + 4QCr−1 − QCr−2 − QAr−2 Cr−1 )/E, r = 2, 3, . . . , N − 1.
(32)
In view of the first boundary condition of (7), we have C0 = W0 ,
A0 = 0,
B0 = 0.
(33)
The third boundary condition can be discretized as W1 − W−1 = 2hγ1 (k),
(34)
with r = 1 and using (34), Eq. (31) will be identical with Eq. (29) if 4Q −Q (Gr=1 ) + 4QW0 + 2hγ1 (k) , B1 = , C1 = . (35) A1 = 1 + 7Q 1 + 7Q 1 + 7Q The scalars Ar , Br and Cr are computed from (32) in a forward sweep manner subject to the starting (n+1) is obtained values in (33) and (35). Using the stored values of Ar , Br and Cr , the solution Wr = Ur,s+1 by a backward sweep manner subject to the initial values WN (known) and WN−1 which can be evaluated using (31) with r = N − 1 and the discretization of the fourth boundary condition in (7). Scheme (16) together with the given initial and boundary conditions, is solved using the following iterative procedure: (I) For the first time level (s = 1) (a) Take Ur = f (rh), Yr = Ur for all r. (b) Solve (29) for Wr using the recurrence relation (31).
166
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
(c) Set Yr = Wr for all r. (d) Calculate improved values of Wr from (29). (e) Repeat steps (b)–(d) till condition (18) is satisfied. (II) For the time levels {s +1 | s = 1, 2, . . .} we repeat steps (b)–(d) except that in step (b), we solve (25) for Wr using (31) subject to the known initial values Ur and Vr for all r and the known boundary values U0,s and UN,s for all s. Scheme (16) is a penta-diagonal system of linear equations having constant matrix of coefficients. This means that only the right-hand of this scheme changes at each time level and this saves computational time and minimizes machine storage.
4. Numerical examples and results 4.1. Example 1 The analytic expression of solutions of Eq. (4) is [8] σ (x ± αt) − x0 + 0.5, u(x, t) = −ω sec h2 2
(36)
where x0 and σ > 0 are real parameters and the amplitude ω and velocity α of the wave are related to σ through the relation 3 α = 1 − σ 2. (37) ω = σ 2, 2 It should be noted that due to the square root in (37), the solitary waves (36) only exist for a finite range of velocities. Using the numerical schemes presented in the previous section, we now solve Eq. (4) together with the initial data 2 σ x + 0.5, g(x) = 0, (38) f (x) = −ω sec h 2 and zero boundary conditions. This is carried out for different values of ω and variety of mesh lengths h and k. As time evolves the initial pulse either splits into two pulses moving in opposite directions or blows-up according to the value of ω (see [7] for more details). In Fig. 1, an initial stationary pulse of amplitude 1.48 breaking-up into two pulses moving in opposite directions is shown. Fig. 2 displays the “blow-up” of an initial stationary pulse of amplitude 1.5. As time evolves, the amplitude of the pulse increases while its breadth decreases. 4.2. Example 2 In this example the interaction of two solitary waves travelling towards each other and colliding headon is illustrated. Consider Eq. (4) with the following initial conditions f (x) = −ω1 sec h2 Ω1 − ω2 sec h2 Ω2 , g(x) = −ω1 σ1 sec h2 Ω1 tan hΩ1 − ω2 σ2 sec h2 Ω2 tan hΩ2 ,
(39)
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
167
Fig. 1. Break-up of the pulse, h = 0.5, k = 0.05, −75 x 75, 0 t 60, x0 = 0, ω = 1.48.
Fig. 2. Blow-up of the pulse. h = 0.5, k = 0.1, −75 x 75, 0 t 10, x0 = 0, ω = 1.5.
where
3 2 α1 αi = 1 − σi2 , Ω1 = x + x0 , ωi = σi , 2 2 the boundary conditions taken at x = xmin are β1 (x) = −ω1 sec h2 ξ1 − ω2 sec h2 ξ2 + 1, γ1 (x) = −ω1 σ1 sec h2 ξ1 tan hξ1 − ω2 σ2 sec h2 ξ2 tan hξ2 ,
Ω2 =
α2 x − x0 , 2
i = 1, 2,
(40)
(41)
168
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
Fig. 3. Collision of two pulses. h = 0.5, k = 0.1, −75 x 75, 0 t 60, x0 = 5, ω1 = 0.25, ω2 = 0.25.
Fig. 4. Collision of two pulses. h = 0.5, k = 0.1, −75 x 75, 0 t 60, x0 = 4, ω1 = 0.35, ω2 = 0.35.
where σ1 σ2 (xmin − α1 t) + x0 , ξ2 = (xmin − α2 t) − x0 . (42) 2 2 On replacing xmin by xmax in (42) we get the boundary conditions β2 (t) and γ2 (t) taken at x = xmax . For the computational work suitable values of σ1 , σ2 , x0 and the number of points N (x-interval) were chosen so that the accuracy within reasonable run time is maintained. With the help of the computer plots, the ξ1 =
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
169
Fig. 5. Collision of two pulses. h = 0.5, k = 0.1, −75 x 75, 0 t 60, x0 = 5, ω1 = 0.35, ω2 = 0.15.
Fig. 6. Collision of two pulses. h = 0.5, k = 0.1, −75 x 75, 0 t 60, x0 = 10, ω1 = 0.38, ω2 = 0.38.
nature of the “good” Boussinesq solitary waves interaction is investigated for different amplitudes of these waves. Figs. 3–6 show typical interactions of the “good” Boussinesq two solitary waves plotted at successive time steps for different amplitude values. In Figs. 3–5, the pulses emerge with their original shape and speed after interaction as if the collision had not taken place. In Fig. 6, “blow-up” of the solution is obtained. Figs. 7 and 8 display the head on collision at selected time steps. It has been shown by Manoranjan et al. [7,8] that the nature of the interaction between any pair of solitary waves is dependent
170
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
Fig. 7. Head-on collision at selected time levels. h = 0.5, k = 0.1, −75 x 75, x0 = 5, ω1 = 0.35, ω2 = 0.35.
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
Fig. 8. Head-on collision at selected time levels. h = 0.5, k = 0.05, −75 x 75, x0 = 10, ω1 = 0.38, ω2 = 0.38.
171
172
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
Fig. 9. Joint amplitude (ωj ) for ω = 0.05(0.0025)0.35, h = 0.5, k = 0.05, x0 = 4, −75 x 75.
on their amplitudes ω1 and ω2 . They reported that, for real solitons interaction the amplitude values ω1 , and ω2 must satisfy certain conditions and when the initial amplitudes of the waves are equal (ω, say) “blow-up” of the solution is witnessed for ω greater than some critical value. It should be noted that when the two waves overlap (Figs. 3 and 4) then the joint amplitude (ωj ) is greater than twice the amplitude of an individual solitary wave (ω). Fig. 9 displays the relation between ωj and ω. The numerical computations presented in this paper were carried out on a personal computer and three to six iterations at each time level were required to satisfy condition (18).
Conclusions On the basis of linearization and finite-difference techniques, an iterative implicit scheme for the “good” Boussinesq equation is considered. Computational time is saved by using this method since a constant coefficient penta-diagonal system of equations is to be solved at each time level. Based on the numerical results and computer plots, the “good” Boussinesq equation has been found to possess an interesting soliton-interaction mechanism. For moderate amplitude values, the incoming waves emerge of the interaction without changing shape or velocity (clean interaction). However, if the amplitudes are large, the interaction leads to the creation of “blow-up” solutions. The numerical method developed in this paper provides results which are in conformity with the general ideas reported in [7] and [8].
References [1] Y. Agnon, P.A. Madsen, H.A. Schoffer, A new approach to high-order Boussinesq models, J. Fluid Mech. 399 (1999) 319–333. [2] W.F. Ames, Numerical Methods for Partial Differential Equations, 2nd Edition, Academic Press, New York, 1965.
H. El-Zoheiry / Applied Numerical Mathematics 45 (2003) 161–173
173
[3] R.E. Bellman, R.E. Kalaba, Quasilinearization and Nonlinear Boundary Value Problems, Elsevier, New York, 1965. [4] A.G. Bratsos, The solution of the Boussinesq equation using the method of lines, Comput. Methods Appl. Mech. Engrg. 157 (1998) 33–44. [5] R. Hirota, Exact N -soliton solutions of the wave equation of long waves in shallow-water and in nonlinear lattices, J. Math. Phys. 14 (1973) 810–814. [6] L. Iskandar, P.C. Jain, Numerical solution of equations having inelastic soliton wave interaction, Comput. Math. Appl. 6 (1980) 373–383. [7] V.S. Manoranjan, A.R. Mitchell, J.L. Morris, Numerical solutions of the “good” Boussinesq equation, SIAM J. Sci. Statist. Comput. 5 (1984) 946–957. [8] V.S. Manoranjan, T. Ortega, J.M. Sang-Serna, Soliton and antisoliton interactions in the “good” Boussinesq equation, J. Math. Phys. 29 (1988) 1964–1968. [9] H.P. Mckean, Boussinesq’s equation on the circle, Comm. Pure Appl. Math. 34 (1981) 599–691. [10] T. Ortega, J.M. Sanz-Sema, Nonlinear stability and convergence of finite-difference methods for the “good” Boussinesq equation, Numer. Math. 58 (1990) 215–229. [11] R.D. Richtmyer, K.W. Morton, Difference Methods for Initial Value Problems, 2nd Edition, Wiley–Interscience, New York, 1967. [12] G. Schneider, The long wave limit for a Boussinesq equation, SIAM J. Appl. Math. 58 (1992) 1237–1245. [13] V.B. Zakharov, Randomization problem for one-dimensional chains of nonlinear operators, Zh. Eksper. Theoret. Fiz. 65 (1973) 219–228 (in Russian).