Kawahara solitons in Boussinesq equations using a robust Christov–Galerkin spectral method

Applied Mathematics and Computation 243 (2014) 245–257

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Kawahara solitons in Boussinesq equations using a robust Christov–Galerkin spectral method M.A. Christou ⇑, N.C. Papanicolaou Department of Mathematics, University of Nicosia, 46 Makedonitissas Ave., 1700 Nicosia, Cyprus

a r t i c l e

i n f o

Keywords: Spectral methods Soliton interaction Kawahara solitons Boussinesq equations

a b s t r a c t We develop a robust Christov–Galerkin spectral technique for computing interacting localized wave solutions of and fourth and sixth-order generalized wave equations. To this end, a special complete orthonormal system of functions in L2 ð1; 1Þ is used whose rate of convergence is shown to be exponential for the cases under consideration. For the time-stepping, an implicit algorithm is chosen which makes use of the banded structure of the matrices representing the different spatial derivatives. As featuring examples, the head-on collision of solitary waves is investigated for a sixthorder generalized Boussinesq equation and a fourth-order Boussinesq type equation with a linear term. Its solutions comprise monotone shapes (sech-es) and damped oscillatory shapes (Kawahara solitons). The numerical results are validated against published data in the literature using the method of variational imbedding. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction In recent years a number of physical problems have led to boundary value problems in infinite domains. For these problems, boundary conditions are not specified at given points, but rather the square of the solution (or some other energy norm) is required to be integrable over an infinite domain. Such solutions belong to the space of square-integrable functions L2 ð1; 1Þ. A typical example is furnished by the soliton solutions of nonlinear evolution equations or generalized wave equations. Here, we consider the sixth-order Boussinesq equation and a fourth-order Boussinesq type equation with a linear term, which have been studied using different numerical techniques (see, e.g. [1–3]). There are significant disadvantages in applying finite difference or finite-element methods to problems in L2 ð1; 1Þ. The source of these issues is the inevitable truncation of the infinite interval to a finite one, which introduces artificial eigenvalue problems that are irrelevant to the original infinite domain. Often, the finite-domain problem has a solution only for some enumerable set of intervals of specific length. It may even be the case that each of the finite-domain approximations only admits a trivial solution, whilst the original problem possesses a nontrivial one, or vice versa. The difficulties connected with the reduction of the infinite interval can be overcome if a spectral method is used with a basis system of localized functions which intrinsically satisfies the requirement that the solution belongs to L2 ð1; 1Þ space. Here we make use of the complete orthonormal (CON) system of functions proposed in [4]. These functions are orthogonal without weight and there is an expression for converting the product of two members of the system into series with respect to the system. This property is crucial because it allows the use of a Galerkin type expansion which is much ⇑ Corresponding author. E-mail addresses: [email protected] (M.A. Christou), [email protected] (N.C. Papanicolaou). http://dx.doi.org/10.1016/j.amc.2014.05.076 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.

246

M.A. Christou, N.C. Papanicolaou / Applied Mathematics and Computation 243 (2014) 245–257

simpler and faster in implementation than pseudo-spectral algorithms. In addition, the Galerkin technique only has truncation error, while the pseudospectral technique is also subject to discretization error and aliasing in addition to the truncation error (see [5]). When the physical problem under investigation is very sensitive to error, then Galerkin techniques are preferable due to their higher accuracy and dependability. Until now, the Galerkin technique based on the said CON system has mainly been used in computations involving solitons with monotonically decaying tails. It was also used to calculate the solutions of the Benjamin–Ono equation and the generalized Benjamin–Ono equations with cubic and higher nonlinearity (see [6,7]). Essentially, these are the result of integrable model equations. However, the question here is not whether fully integrable systems admit solutions with solitonic behavior. What spoils integrability is the presence of higher-order dispersive terms which brings into existence more complicated shapes, e.g. solitons with oscillatory non-monotone tails. These were named after Kawahara who discovered them in 1978 [8] whilst investigating the fifth-order KdV equation

ut þ au2 ux þ buxxx  cuxxxxx ¼ 0; where a; b; c are constants: A lot of work has been done on KdV and KdV-type equations. For instance, Wazwaz[9,10] studied modified and higherorder KdV equations with time-dependent coefficients, Biswas[11] the above modified Kawahara-KdV equation and BonaSaut[12] examined the blowup of solutions of the generalized KdV equation. In [13] a ninth-order KdV equation and a sixth-order Boussinesq equation were derived by generalizing their bilinear forms. The main purpose of this paper is to extent and validate the proposed Galerkin technique and utilize it in the numerical investigation of the aforementioned sixth- and fourth- order Boussinesq equations. The first equation considered here models the propagation of shallow water waves. It is known as the sixth-order generalized Boussinesq equation or 6GBE (see [1]) because it contains a sixth-order derivative. The second equation examined here, governs the flexural deformation of an elastic rod on elastic foundation. It is a fourth-order Boussinesq type equation and was named EREFE (see [2]). Both equations reveal very interesting results, in particular, soliton-type solutions with oscillatory shapes. In the case of 6GBE, this is due to the sixth-order derivative, whilst in the case of EREFE, this is due to the presence of a linear term, which acts as a linear restoring force. 2. The Christov–Galerkin method for Lð‘; ‘Þ 2.1. The complete orthonormal set of functions As alluded to earlier, the spectral technique utilized here, is of Fourier–Galerkin type (for other spectral techniques, see, e.g., [5,14]). The Galerkin approach has the advantage of simplicity in implementation compared to the spectral collocation method or the tau-method [5]. The trade-off is that for nonlinear problems, the Galerkin technique requires explicit formulas expressing the products of members of the CON system into series with respect to the system. For instance, Hermite and Laguerre polynomials do not possess this kind of explicit relation. The first system in L2 ð1; 1Þ for which such a product formula exists was proposed in [4]. A Galerkin technique based on the said system was developed in [15] and applied to Korteweg–de Vries (KdV) and Kuramoto–Sivashinsky (KS) equations with quadratic nonlinearity. It has been recently applied to 2D problems [16] and to problems with cubic nonlinearities [17]. In a sequence of papers, Boyd showed a general way of constructing CON systems in L2 ð1; 1Þ based on coordinate transformations to finite intervals and the utilization of Chebyshev polynomials (see the comprehensive monograph of the same author [5] for references). The CON system we use in the present work was constructed in [4], as the real and imaginary parts of the Wiener functions [18]

1

ðix  1Þ

n

qn ¼ pffiffiffiffi ; n ¼ 0; 1; 2; . . . p ðix þ 1Þnþ1

ð1Þ

which were introduced by Wiener as the Fourier transforms of the Laguerre functions (functions of parabolic cylinder). Higgins [19] defined them for negative indices and proved the completeness and orthogonality of the system. These functions are orthogonal without weight and possess expressions for the product of two members of the system into series with respect to the system

qn qk ¼

qnþk  qnk pffiffiffiffi : 2 p

ð2Þ

Two real-valued subsequences of odd fSn g and even fC n g functions were introduced [4], namely

Sn ¼

qn þ qn1 pffiffiffi i 2

;

Cn ¼

qn  qn1 pffiffiffi 2

;

n ¼ 0; 1; 2; . . .

ð3Þ

Naturally, the product property (2) is inherited by the real sequences fC n g and fSn g, namely,

1 C n C k ¼ pffiffiffiffiffiffiffi ½C nþkþ1  C nþk  C nk þ C nk1 ; 2 2p

ð4Þ

M.A. Christou, N.C. Papanicolaou / Applied Mathematics and Computation 243 (2014) 245–257

247

1 Sn Sk ¼ pffiffiffiffiffiffiffi ½C nþkþ1  C nþk þ C nk  C nk1 ; 2 2p

ð5Þ

1 Sn C k ¼ pffiffiffiffiffiffiffi ½Snþkþ1 þ Snþk þ Snk  Snk1 : 2 2p

ð6Þ

At the time of its discovery in 1982, the above set fC n ; Sn g was the only set of square-integrable functions over ð1; 1Þ which had the product property (2), making it indispensable for creating Galerkin spectral algorithms. Originally, the system (3) was successfully applied for finding localized (solitary-wave) solutions of KdV and Kuramoto–Sivashinsky equations in [15] and to the fifth-order KdV in [20]. The most important result was that a good approximation to the solution could be obtained even when taking fewer than 10 terms in the Galerkin expansion. This made the CON system (3) a very useful tool for low-cost accurate computations of localized solutions. When retaining only 10–20 terms in the expansion, the superiority of the Galerkin technique over the pseudo-spectral approach is even more evident. The following representation for members of (3) is known [15]

C n ðxÞ ¼ ð1Þn

cosðn þ 1Þh þ cos nh pffiffiffi ; 2

Sn ðxÞ ¼ ð1Þnþ1

sinðn þ 1Þh þ sin nh pffiffiffi ; 2

ð7Þ

where x ¼ tanðh=2Þ (or h ¼ 2 arctan x) is a transformation of the independent variable. The importance of this connection to the periodic functions will become apparent in Section 2.3. It is also discussed in [21]. 2.2. Properties of the CON system For the sake of brevity of the notation, we introduce the following formal expressions for the products Sn Sm ; C n C m , and Sn C m , from (4)–(6):

CnCk ¼

1 X bnk;m C m ;

Sn Sk ¼

m¼1

1 X

ank;m C m ; Sn C k ¼

m¼1

1 X

cnk;m Sm ;

ð8Þ

m¼1

1

ank;m ¼ pffiffiffiffiffiffiffi fdm;nþkþ1 þ dm;jnkj  dm;nþk  sgn½jn  kj  0:5dm;½jnkj0:5 g; 2 2p 1 bnk;m ¼ pffiffiffiffiffiffiffi fdm;nþk þ dm;jnkj  dm;nþkþ1  sgn½jn  kj  0:5dm;½jnkj0:5 g; 2 2p 1

cnk;m ¼ pffiffiffiffiffiffiffi fdm;nþk þ sgnðn  kÞdm;jnkj  dm;nþkþ1  sgnðn  kÞdm;jnkj1 g: 2 2p The expansion formulas for the second and fourth derivatives were presented in [4] and are given by

C IIn ¼

1 X

1 X

m¼0

m¼0

vm;n C m ; SIIn ¼

vm;n Sm ;

1 4

ð9Þ 1 4

1 4

vm;n ¼  nðn  1Þdm;n2 þ n2 dm;n1  ðn þ 1Þðn þ 2Þdm;nþ2  n2 þ ð2n þ 1Þ2 þ ðn þ 1Þ2 dm;n þ ðn þ 1Þ2 dm;nþ1 ; and

C IV n ¼

1 X

1 X

m¼0

m¼0

xm;n C m ; SIV n ¼

xm;n ¼

xm;n Sm ;

ð10Þ

1 1 1 nðn  1Þðn  2Þðn  3Þdm;n4  nðn  1Þ2 ðn  2Þdm;n3 þ nðn  1Þð7n2  7n þ 4Þdm;n2 16 2 4 1 2 1 1 2 4 3 2  n ð7n þ 5Þdm;n1 þ ½35n þ 70n þ 85n þ 50n þ 12dm;n  ðn þ 1Þ2 ½7ðn þ 1Þ2 þ 5dm;nþ1 2 8 2 1 1 2 þ ðn þ 1Þðn þ 2Þ½7ðn þ 1Þ þ 7ðn þ 1Þ þ 4dm;nþ2  ðn þ 1Þðn þ 2Þ2 ðn þ 3Þdm;nþ3 4 2 1 þ ðn þ 1Þðn þ 2Þðn þ 3Þðn þ 4Þdm;nþ4 ; 16

respectively. The formulas for the expansions of the sixth derivatives of C n and Sn into members of the series were not available in the literature and had to be derived. They read

248

M.A. Christou, N.C. Papanicolaou / Applied Mathematics and Computation 243 (2014) 245–257

C VI n ¼

1 X wm;n C m ; m¼0

SVI n ¼

1 X wm;n Sm ;

ð11Þ

m¼0

1 3 nðn  1Þðn  2Þðn  3Þðn  4Þðn  5Þdm;n6 þ nðn  1Þðn  2Þ2 ðn  3Þðn  4Þdm;n5 64 16 3 5  nðn  1Þðn  2Þðn  3Þð11n2  33n þ 31Þdm;n4 þ nðn  1Þ2 ðn  2Þð11n2  22n þ 27Þdm;n3 32 16 45 9  nðn  1Þð11n4  22n3 þ 43n2  32n þ 12Þdm;n2 þ n2 ð11n4 þ 35n2 þ 14Þdm;n1 64 8 1  ð231n6 þ 693n5 þ 1680n4 þ 2205n3 þ 1869n2 þ 882n þ 180Þdm;n 16 9 45 þ ðn þ 1Þ2 ð11n4 þ 44n3 þ 101n2 þ 114n þ 60Þdm;nþ1  ðn þ 1Þðn þ 2Þð11n4 þ 66n3 þ 175n2 þ 228n þ 120Þdm;nþ2 8 64 5 3 2 2 þ ðn þ 1Þðn þ 2Þ ðn þ 3Þð11n þ 44n þ 60Þdm;nþ3  ðn þ 1Þðn þ 2Þðn þ 3Þðn þ 4Þð11n2 þ 55n þ 75Þdm;nþ4 16 32 3 1 2 þ ðn þ 1Þðn þ 2Þðn þ 3Þ ðn þ 4Þdm;nþ5  ðn þ 1Þðn þ 2Þðn þ 3Þðn þ 4Þðn þ 5Þðn þ 6Þdm;nþ6 : 16 64

wm;n ¼ 

It is important to point out that the derivatives of the basis functions are reduced to banded matrix operators, which adds to the computational efficiency of the here technique. The positive definiteness of matrix vnm and bounds on the eigenvalues were shown in [21]. These properties are essential for constructing the implicit scheme under consideration. Note that with the exception of the last couple of lines, the matrix x is the square of v and the matrix w is the square of x. 2.3. Convergence analysis An important issue for a spectral method is its rate of convergence. There are different ways to show the rate of convergence for the set fC n ; Sn g. We find the argument employed by Boyd [5, Section 2.11] to show the exponential rate of convergence of Chebyshev series, to be the most straightforward. Quoting directly from [5, Section 2.9]: Theorem 4 (Integration-by-parts Coefficient Bound). If  f ðpÞ ¼ f ðpÞ; f 0 ðpÞ ¼ f 0 ðpÞ; . . . ; f ðk2Þ ðpÞ ¼ f ðk2Þ ðpÞ,  f ðkÞ ðxÞ is integrable, then the coefficients of the Fourier series

f ðxÞ ¼ a0 þ

1 1 X X an cosðnxÞ þ bn sinðnxÞ; n¼1

n¼1

have the upper bounds

jan j 6 F=nk ;

jbn j 6 F=nk ;

for some sufficiently large constant F, which is independent of n. An equivalent way of stating the theorem is that, if the two conditions above are satisfied, then the algebraic index of convergence is at least as large as k. The proof consists of k integrations by parts and implementation of the first condition. What is important is the discussion of the theorem’s consequences (again, see, [5, Section 2.9]) the most important points of which can be quoted as follows:  If the function f ðxÞ is periodic and differentiable to all orders, for large n, the Fourier coefficients decrease faster than any finite power of n, which is known as exponential convergence.  Since a Chebyshev series in x is a Fourier series in y following the transformation x ¼ cosðyÞ a Chebyshev series always has the property of exponential convergence (in x, the Chebyshev argument) if all its derivatives are bounded on the expansion interval, x 2 ½1; 1. To apply the above argument to the set of Christov functions, we use the relationship (7), between the Fourier periodic functions and our system. This implies that any function f ðxÞ expanded in our CON is a periodic function of h with period 2p under the transformation of coordinates (7). In particular, any localized function f ðxÞ can be developed into a series with respect to fC n ; Sn g,

f ðxÞ ¼

1 X an C n ðxÞ þ bn Sn ðxÞ; i¼0

which in turn may be rewritten as a Fourier series for the periodic function

M.A. Christou, N.C. Papanicolaou / Applied Mathematics and Computation 243 (2014) 245–257

249

  X 1 h cosðn þ 1Þh þ cos nh sinðn þ 1Þh þ sin nh pffiffiffi pffiffiffi þ bn ð1Þnþ1 : ¼ f tan an ð1Þn 2 2 2 i¼0 Fourier series expansions of periodic functions are known to have exponential convergence, therefore, the exponential convergence of the fC n ; Sn g series follows. 3. Governing equations Consider the sixth-order Boussinesq equation (6GBE, [1]):

utt ¼ c21 uxx þ ða1 u2 Þxx þ b1 uxxxx þ d1 uxxxxxx ;

ð12Þ

where the sixth-order dispersion coefficient is always positive, d1 > 0. In [1] it was argued that sixth-order dispersion appears naturally both for lattice dynamics and shallow water flows. Note that for a1 ¼ b1 ¼ 1 and d1 ¼ 0, the equation is reduced to the proper Boussinesq equation [22]. Eq. (12) is non-integrable and admits different types of localized solutions for different values of the governing parameters. Therefore, a numerical approach is needed for both stationary propagating waves and for unsteady interactions of localized waves. At first, we solve the problem in the moving frame (stationary problem), a frame of reference which moves with the observer along a trajectory, using the Christov–Galerkin method. Then, the obtained solutions are used as initial conditions to investigate the time-dependent problem. One may easily see that (12) can be written as the following system if the auxiliary function q is introduced

ut ¼ qxx ;

ð13Þ

qt ¼ c21 u þ a1 u2  b1 uxx þ d1 uxxxx :

To start our investigation, we first consider the case of stationary waves in the moving frame n ¼ x  ct. After integrating twice, one arrives at the following nonlinear ODE

k1 u þ a1 u2 þ b1 unn þ d1 unnnn ¼ 0;

for k1 ¼ c21  c2 :

ð14Þ

The last equation is subject to the asymptotic boundary conditions

uðt; xÞ ! 0;

for n ! 1:

ð15Þ

The linearized version of (20) admits harmonic solutions of the type

j4 þ bj2 þ k1 ¼ 0;

ejn .

The corresponding dispersion relations read

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i:e; m ¼ b  b2  4k1 =2 where m ¼ j2 :

ð16Þ

Eq. (16) shows that the classification of the stationary solutions should be based on criteria which define the asymptotic pffiffiffiffiffi behavior of the tails. These criteria depend on the values of k1 and b. Now, if b < 2 k1 for k1 > 0, then b2  4k1 > 0 which results in monotone solutions. When this is the case, (20) possesses an exact solution, namely,

uan ðnÞ ¼ ð105b2 =338aÞ sech

4



pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n=2 b=13 ;

jcj ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c2  36b2 =169:

ð17Þ

pffiffiffiffiffi pffiffiffiffiffi If 2 k1 6 b 6 2 k1 then damped oscillatory shapes appear. For b > 0 (e.g. b ¼ 1) and k1 P b2 =4, the wave number k from (16) is complex, say k ¼ kr þ iki . Then the asymptotic behavior of the solution at infinity is

u ¼ sinðki x þ dÞekr x :   This means that the localized solutions are damped oscillations because for b2  4jk1 j  1 (or 4jk1 j ¼ b2 þ e2 ) the wave support L is large

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi 2j ¼  b  ie;

 pffiffiffiffiffiffi kr ’ e= 2 2b ;

)

)

L ’ e1  1:

The second equation under investigation is the equation modeling the flexural deformation of an elastic rod on an elastic foundation (EREFE), namely,



v tt ¼ rv  a2 v 2  c2 v xx

xx

 d2 v ;

d2 > 0:

ð18Þ

EREFE, (18), is only of fourth order, but due to the presence of the linear restoring force, exhibits some of the phenomenology of 6GBE, especially for the case of non-monotonic localized solutions. However, the linear term representing the restoring force, qualitatively changes the equation for stationary waves. The problem in the moving frame (stationary case), was treated successfully with the aid the of method of variational imbedding in [3]. Here, using the aforementioned Christov–Galerkin spectral method, we rediscover the results obtained in [3] and extend our analysis to the time-dependent problem where we study the propagation and interaction of the solitary waves.

250

M.A. Christou, N.C. Papanicolaou / Applied Mathematics and Computation 243 (2014) 245–257

As in the case of the 6GBE, at first we seek traveling wave solutions. We therefore introduce the new variable g ¼ x  e ct into (18), obtaining

k2 v gg  a2 ðv 2 Þgg þ k3 v gggg  d2 v ¼ 0;

for k2 ¼ r  e c ; k3 ¼ c3 :

ð19Þ

Note that re-scaling the spatial variable x does not change the nature of the asymptotic boundary value problem. Introducing x ¼ f1 ^ x into 6GBE and x ¼ f2 ^ x into EREFE we obtain the scaled form of (14) and (19), namely,

k1 u þ a1 u2 þ f21 b1 unn þ f41 d1 unnnn ¼ 0;

for k1 ¼ c21  c2 ;

ð20Þ

and

f22 k2 v gg  f22 a2 ðv 2 Þgg þ f42 k3 v gggg  d2 v ¼ 0;

for k2 ¼ r  e c ; k3 ¼ c3 :

ð21Þ

respectively. The significance of scaling parameters f1 ; f2 is that they can be used to optimize the method in the sense that their introduction allows one to match the typical length scales of the employed system of functions, with that of the support of the sought localized solution. Naturally, this improves the convergence rate of the Galerkin series and accelerates the overall convergence of our method.

4. Application of the Galerkin spectral method 4.1. Stationary problems The problems under investigation, (20) and (21), admit even functions as solutions. Hence, we develop the sought solution u into series with respect to the subsequence of even functions C n only [4], namely

u¼

1 X an C n ðnÞ;

u¼

n¼0

1 X

f an C n ðgÞ:

ð22Þ

n¼0

Then, the various terms that appear in (20) and (21) have the following expansions

uII ¼

1 X 1 X

am vm;n C n ðnÞ;

m¼0 n¼0

u2 ¼

1 X 1 X am xm;n C n ðnÞ; m¼0 n¼0

1 X 1 X 1 X am1 am2 bm1 m2 ;n C n ðnÞ; m1 ¼0m2 ¼0 n¼0

II

uIV ¼

ðu2 Þ ¼

uVI ¼

1 X 1 X am wm;n C n ðnÞ;

ð23Þ

m¼0 n¼0

1 X 1 X 1 X 1 X am1 am2 bm1 m2 ;m vm;n C n ðnÞ; m1 ¼0 m2 ¼0 m¼0 n¼0

bnk;m

 ¼ dm;nþkþ1 þ dm;nþk þ dm;jnkj  sgnðjn  kj  0:5Þdm;½jnkj0:5 =a:

Introducing (22), (23) into (20), (21) and taking successive inner products with C i ; i ¼ 0; 1; 2; 3; . . ., we obtain the followan : ing infinite nonlinear algebraic systems for the unknown coefficients an and f 1 h X

1 1 i X X f41 d1 xm;n þ f21 b1 vm;n þ k1 dm;n am þ a1 am1 am2 bm1 m2 n ¼ 0;

m¼0 1 h X

ð24Þ

m1 ¼0 m2 ¼0 1 1 1 i X X X 2 g ag f42 k3 xm;n þ f22 k2 vm;n  d2 dm;n af m  f2 a2 m1 a m2 bm1 m2 m3 vm3 ;n ¼ 0;

m¼0

ð25Þ

m1 ¼0 m2 ¼0 m3 ¼0 1

f where fan g1 n¼0 and f an gn¼0 are the unknown coefficients for (20) and (21) respectively. In the numerical calculations, a truncated version of the above system is used, where terms are retained up to a specified finite index N. To investigate (13) and (18), we need initial conditions. If exact solutions exist, a single solution or the superposition of two or more such solutions may be used as initial conditions. However, in general, such solutions are not available. Therefore, to solve the time-dependent problem we use as initial condition the numerical solutions obtained in the stationary problem. In our case, this means the computational solutions of (24) and (25). The investigation of an overcoming or head-on interaction of solitons requires the evaluation of the two solitons in two different positions in order to obtain a composite initial shape. Then, the initial spectral coefficients are obtained with the aid of numerical quadratures. A similar approach is used when studying single soliton propagation. The initial condition for function q is obtained by solving the first equation of system (13) numerically.

251

M.A. Christou, N.C. Papanicolaou / Applied Mathematics and Computation 243 (2014) 245–257

4.2. Non-stationary problems To solve (13) and (18) we use the following semi-implicit time-stepping schemes of Crank–Nicholson type for the linear terms,



 lþ1 ulþ1  ul =s ¼ f21 qxx 2 ;

ð26Þ

 1 





lþ1  2 1 4 l1 l1 qlþ2  ql2 =s ¼ c21 ulþ1 þ ul1 =2 þ a1 ðul Þ  f21 b1 ulþ1 xx þ uxx =2 þ f1 uxxxx þ uxxxx =2; and







l1 v lþ1  2v l þ v l1 =s2 ¼ f22 r v lþ1 xx þ v xx



2

=2  f22 a2 ðv l Þxx  f42 c2



l1 v lþ1 xxxx þ v xxxx



=2  d2



v lþ1 þ v l1

 =2:

ð27Þ

respectively. We develop, for (26), the sought solutions u and q into series with respect to the subsequences C n and Sn , namely,

ul ðzÞ ¼

1 X l aln C n ðzÞ þ bn Sn ðzÞ;

1

qlþ2 ðzÞ ¼

n¼0

v l ðzÞ ¼

ð28Þ

m¼0

In like fashion, we develop the solution 1 X

1 X lþ1 lþ1 dm 2 C m ðzÞ þ em 2 Sm ðzÞ:

v of (27) in the form

e b ln Sn ðzÞ: a ln C n ðzÞ þ e

ð29Þ

n¼0

Then, we substitute the spectral expansions (28), (29) into systems (26), (27) and make use of the orthogonality of the system of functions, to derive algebraic systems for the coefficients. The systems for the even and odd function coefficients for (26), respectively, read:



1 X  lþ1 2 l alþ1 dk 2 vm;k ; m  am =s ¼ f1

ð30Þ

k¼0 1 1  1  X X

lþ1 





 lþ l1 l1 dm 2  dm 2 =s ¼ f21 ðb1 =2Þ ak þ al1 vk;m þ f41 ðd1 =2Þ alþ1 þ al1 xk;m þ c21 =2 alþ1 k k m þ am k k¼0

k¼0

1 X 1   X l l þ a1 alk1 alk2 bk1 k2 ;m þ bk1 bk2 ak1 k2 ;m ; k1 ¼0 k2 ¼0 1   X lþ1 lþ1 l bm  bm =s ¼ f21 ek 2 vm;k ;

ð31Þ

k¼0 1  1   lþ1      X X l1 lþ1 l1 lþ1 l1 lþ1 l1 em 2  em 2 =s ¼ f21 ðb1 =2Þ bk þ bk vk;m þ f41 ðd1 =2Þ bk þ bk xk;m þ ðc21 =2Þ bm þ bm k¼0

þ a1

1 X 1  X

k¼0

 l l alk1 bk2 ck1 k2 ;m þ bk1 alk2 ck1 k2 ;m :

k1 ¼0 k2 ¼0

Similarly, the algebraic systems for the even and odd function coefficients for (27) is given by



1 1 X X 

lþ1 

lþ1 

 e el e l1 =s2 ¼ f22 ðr=2Þ e e a lþ1 an þ e a l1 an þ e a l1 vm;n  f42 ðc2 =2Þ xm;n  ðd2 =2Þ ea lþ1 þ ea l1 m  2am þ am n n n¼0

 f22 a2

1 X 1 X 1  X

n¼0

e b lk1 e b lk2 ak1 k2 ;m a lk1 e a lk2 bk1 k2 ;m þ e



vm;n ;

ð32Þ

n¼0 k1 ¼0 k2 ¼0

and



1  1       X X e el e l1 =s2 ¼ f2 ðr=2Þ e e l1 v  f4 ðc =2Þ e e l1 xm;n  ðd2 =2Þ e b lþ1 b lþ1 b lþ1 b lþ1 þ e b l1 2 2 2 m  2bm þ bm n þ bn n þ bn m;n n¼0

 f22 a2

1 X 1 X 1  X n¼0 k1 ¼0 k2 ¼0

respectively.

n¼0

e b lk2 bk1 k2 ;m þ e b lk1 e a lk1 e a lk2 ak1 k2 ;m



vm;n :

ð33Þ

252

M.A. Christou, N.C. Papanicolaou / Applied Mathematics and Computation 243 (2014) 245–257 lþ12

lþ1

1

; elþ2 are inextricably coupled in the above systems. In fact, for the purpose n o n o lþ1 lþ12 lþ1 and bm ; em 2 are assembled. Then (30) and (31) of Gaussian elimination for each system, two composite vectors alþ1 m ; dm Note here that the unknowns alþ1 ; d

and b

yield two coupled sixteen-diagonal algebraic systems for these vectors. This makes the scheme strongly implicit. As we have already mentioned, our initial conditions are superpositions of numerical solutions of the stationary prob n 0o 1 lems. The initial spectral coefficients a0n ; bn and a1n , fbn g are calculated for t ¼ 0 and t ¼ s by evaluating the inner n 1o n 1o product of these solutions with C n or Sn using quadrature. In turn, the initial conditions for the coefficients d2n , e2n of q are computed via numerical quadrature after the first equation of the system (13) is multiplied by C n or Sn . Thus, having n o 3 prepared the initial conditions we begin the iterative process with solving (30) and (31) for l ¼ 1, i.e. we compute a2n ; d2n n o 3 2 and bn ; e2n . After these systems are solved and one time step is completed, the time index l is reset, and the process is repeated. The same process is carried out for coefficients e bn. a l and e 5. Results and discussion 5.1. Numerical verification of spectral convergence In Section 2, we proved that the convergence rate of a Fourier–Galerkin series using the CON system of Christov functions is exponential. It is important to verify this numerically. We explore the convergence of the spectral series expansions our sought functions by taking a close look of the way a superposition of two sech solitons is expanded into series. As a featuring example we take the case cl ¼ 0:95; cr ¼ 0:9 for 6GBE and cel ¼ 0:1; cer ¼ 0:95 for EREFE. Note that for both cases the solubn are non-trivial. The profile under consideration is a superposition of two an ; f tion is neither odd nor even and an ; bn and f sech-es and it is expanded into the CON system using numerical quadratures with 100,000 grid points. This ensures that the truncation error of the numerical quadrature is negligible. Since in the initial moment the solitons are separated by the largest distance (the toughest case for the spectral expansion), we can use the initial profile to tune the parameters of the method. The effect of scaling parameter f can be examined by computing the coefficients of the expansion of the initial condition for different values of f. Naturally, its optimal value is different for different initial configurations of the system of solitons. When the solitons are situated further from each other, one is faced with a system which is not tightly localized. Then the optimal value of f is smaller. Conversely, if the initial configuration is tight enough, the value of f which brings the scales of the sought function and the CON system in a closest rapport, tends to be larger. In the present work we consider systems of solitons that are well separated (in order not to overlap significantly), but not excessively far from each other (so that we do not loose localization). This is the most import case from a physical point of view. In each particular case the investigation begins with choosing the optimal f. Figs. 1 and 2 show the dependence of the spectral coefficients on the number of terms for the soliton systems described earlier. For the particular case of 6GBE discussed here, the round-off error of the calculations is reached using 180 terms, whereas for 4BTE the round-off error was reached in about 500 terms. Because the results presented in Figs. 1 and 2 are scaled by the maximum values of the respective coefficients (not necessarily the lowest-number coefficients), the best-fit curves have different amplitudes. In Tables 1 and 2 we present the best-fit formulas for the lines drawn in Figs. 1 and 2 respectively.

1 0.01 0.0001

0.0001

1e-06

1e-06

1e-08 1e-10

1e-08 1e-10

1e-12

1e-12

1e-14

1e-14

1e-16

1e-16

1e-18

0

50

100

150 n

200

0.3 0.1 0.06

0.01

scaled |bn|

scaled |an|

1

0.3 0.1 0.06

250

1e-18

0

50

100

150

200

250

n

Fig. 1. Log-linear plot of scaled spectral coefficients for the initial condition of 6GBE for different values of the scale parameter f. Left panel: magnitude of coefficients jan j of of even functions C n scaled by their maximum value maxn jan j. Right panel: coefficients jbn j of odd functions Sn scaled by maxn jbn j. The formulas of the best-fit curves are given in Table 1.

253

1

1

0.01

0.01

0.0001

0.0001 scaled |b~n|

scaled |a~n|

M.A. Christou, N.C. Papanicolaou / Applied Mathematics and Computation 243 (2014) 245–257

1e-06 1e-08 1e-10

1e-14

0

1e-08 1e-10

0.5 0.1 0.03

1e-12

1e-06

0.5 0.1 0.03

1e-12

100

200

300

400

1e-14

500

0

100

200

300

n

400

500

n

Fig. 2. Log-linear plot of scaled spectral coefficients for the initial condition of 4BTE for different values of f. Left panel: scaled coefficients j f an j=maxn j f an j of C n . Right panel: scaled coefficients j f bn j=maxn j f bn j of Sn . The best-fit curve formulas can be found in Table 2.

Table 1 Best-fit formulas for spectral coefficients of initial condition of 6GBE. f

Coefficients of C n (Fig. 1 left)

Coefficients of Sn (Fig. 1 right)

0.3 0.1 0.06

0:2 expð0:045nÞ 8 expð0:18nÞ 7 expð0:27nÞ

0:2 expð0:045nÞ 8 expð0:18nÞ 7 expð0:27nÞ

Table 2 Best fit curve formulas for spectral coefficients of initial condition of 4BTE. f

Coefficients of C n (Fig. 2 left)

Coefficients of Sn (Fig. 2 right)

0.5 0.1 0.03

0:75 expð0:005nÞ 2 expð0:035nÞ 1:5 expð0:06nÞ

0:25 expð0:005nÞ 2 expð0:035nÞ 1:5 expð0:06nÞ

5.2. Method validation The method is validated against published numerical data in the literature and against exact solutions for some cases. An exact solution is known for some cases of the 6BGE, (17), and for these particular values we defined the error function to be the absolute error, i.e. the absolute value of the difference between the exact and the numerical solution. Fig. 3 shows the error function for two different values of the scaling parameter f, i.e. f ¼ 1 and f ¼ 8. The error function was evaluated using 80 members from the infinite series. We can see that the error is of order 104 in the case where f ¼ 1 and of order 107 6e-07

1.5e-04

5e-07 4e-07

5e-05

error function

error function

1e-04

0 -5e-05

2e-07 1e-07

-1e-04 -1.5e-04 -100

3e-07

0

-50

0 ξ

50

100

-1e-07 -100

-50

0

50

100

ξ

Fig. 3. Absolute error of the 6GBE for the case when c ¼ 0:88712 for two different values of the scaling parameter f. Left panel: f ¼ 1. Right panel: f ¼ 8.

254

M.A. Christou, N.C. Papanicolaou / Applied Mathematics and Computation 243 (2014) 245–257

8

6

our results Christov et al

7 6

4

4

3

3

2

υ

υ

5

2

1

1

0

0

-1

-1 -2

our results Christov et al

5

0

2

4

6

8

-2

10

0

2

4

η

6

8

10

η

Fig. 4. Comparison with data published in [3] using the same values for the parameters. Here the linear term coefficient is d ¼ 1:6. Left panel: k2 ¼ 0:5. Right panel: k2 ¼ 0:5.

when f ¼ 8. The results clearly show the efficiency of the method since only 80 functions were used, but also demonstrate the need importance of the scaling parameter. For the 4BTE there is no exact solution available and the comparison is done with the computational solutions obtained in [3]. In [3] the stationary shapes of the 4BTE (21) are computed using the method of variational imbedding (MVI). The latter was introduced in [23] and applied to finding the homoclinic solution of the Lorentz equations. In this method the original (inverse) problem is replaced by the Euler–Lagrange equations for the minimization of a quadratic functional of the governing equations. In [24], the MVI was successfully applied to an equation governing Benard–Marangoni convection finding the phase speed for which a solitary wave exists as well as its profile. In [2], the MVI was applied to find solitary-wave solutions of the sixth-order generalized Boussinesq equation 6GBE. The comparison of our results with the results of the MVI is depicted in Fig. 4. Here we are comparing our results with results published in [3] for two different cases of the parameter k2 for the same value of d. In both cases our results are found to be in very good agreement with the available exact solutions and with data already published in the literature obtained using other numerical methods. 5.3. Stationary solitons The strategy to solve the problems at hand is to first solve both problems in the moving frame, meaning that our problems are essentially time-independent. We then use the obtained solutions as initial condition to solve the time-dependent problem and examine the head-on collision. In Fig. 5 we present some of our numerical results for 6GBE for various values of the governing parameters. In the left panel we consider b1 ¼ 1; d1 ¼ 1 and change the phase speed c. For this set of parameters, we observe that no solutions exist if the threshold c ¼ 0:866 is exceeded. However, for different values of b and d increasing c beyond the threshold leads

1.2

0.4

c=0.1 c=0.5 c=0.7 c=0.86

1.0

0.3

0.8

0.25 0.2 u

u

0.6 0.4

0.15 0.1 0.05

0.2

0

0 -0.2 -30

β=0.2 β=0.4 β=0.6 β=0.8

0.35

-0.05 -20

-10

0 ξ

10

20

30

-0.1 -30

-20

-10

0 ξ

10

Fig. 5. Kawahara solitons obtained solving the 6GBE. Left panel: using different c. Right panel: using different b.

20

30

255

M.A. Christou, N.C. Papanicolaou / Applied Mathematics and Computation 243 (2014) 245–257

3.5

2.5

0.8

2.0

0.6

1.5

0.4

1.0

0.2

0.5

0

0

-0.2

-0.5

-0.4

-1.0

-0.6

-1.5 -30

-20

-10

0

10

20

δ=0.9 δ=0.94 δ=0.98 δ=1.0

1.0

υ

υ

1.2

c=1.6 c=1.7 c=1.72 c=1.734

3.0

30

-0.8 -30

-20

-10

0

η

10

20

30

η

Fig. 6. Kawahara solitons obtained solving the EREFE. Left panel: using different c. Right panel: using different b.

to solutions which become more oscillatory as c increases. In the right panel of the same figure, we plot the transition of solitons from monotone shapes to solitons with oscillatory tails. This is achieved if we fix the value of c (here c ¼ 1:7) and change b, the coefficient of the fourth derivative. As b increases from negative to positive the oscillatory aspect of the tails becomes more pronounced. In Fig. 6 we depict results obtained from the EREFE equation. We can see how the profile of our Kawahara solitons is affected by varying the values of the governing parameters. The effect of the phase speed c is depicted in the left panel, and the effect of the coefficient of the linear term d2 is shown in the right panel. Our results are in good agreement with the results published in [3]. The shapes of the EREFE and the 6GBE look alike even though they are two totally different equations. In the 6GBE the term responsible for these kind of shapes is the sixth order derivative. This term is one of the reasons for the stability of the linear problem [1]. In the EREFE the term responsible for solitons with oscillatory tails is the linear term. The linear term is acting as a linear restoring force, qualitatively changing the equation for the stationary waves [3]. 5.4. Moving solitons The results presented in the previous section are obtained by solving Eqs. (20) and (21). We are now ready to set up our initial condition and examine Eqs. (12) and (18). The initial condition in this paper is the superposition of two numerical solutions which are situated far enough from each other allowing us to neglect their intersection at t ¼ 0. 0.3

u

0.2

t=0

t=400

t=900

t=1200

t=600

0.1

0

-0.1 0.3

u

0.2

t=1700

0.1 0 -0.1

-120-80 -40 0 40 80 120 x

-120-80 -40 0 40 80 120 x

-120-80 -40 0 40 80 120 x

Fig. 7. Interaction of Kawahara solitons for c1 ¼ 0:81 and c2 ¼ 0:86.

256

M.A. Christou, N.C. Papanicolaou / Applied Mathematics and Computation 243 (2014) 245–257

1.6 1.2

t=0

t=300

t=450

t=600

t=900

t=950

-60 -40 -20 0 20 40 60 x

-60 -40 -20 0 20 40 60 x

υ

0.8 0.4 0 -0.4

1.6 1.2 υ

0.8 0.4 0 -0.4 -60 -40 -20 0 20 40 60 x

Fig. 8. Interaction of Kawahara solitons for c1 ¼ 1:73 and c2 ¼ 1:62.

The first interaction of moving solitons examined here is from the 6GBE. In Fig. 7 we considered a case of two solitons with different amplitude. We placed the first soliton with c1 ¼ 0:81 at x ¼ 80 and the second one with c2 ¼ 0:86 at x ¼ 80. We then allow the solitons to propagate. After 600 time steps they start meeting and are fully separated after 1700 time steps. We observe that the separation was achieved without any changes in the initial hight or shape of both solitons. The second example of soliton collision is from the EREFE. In Fig. 8 we examine the interaction of two solitons with different amplitude initially placed at x ¼ 50 and x ¼ 50. The corresponding values of c are c1 ¼ 1:73 and c2 ¼ 1:62 respectively. Full separation is obtained after 950 time steps. 6. Conclusions The unsteady problem governing the interaction of solitary wave solutions of the sixth-order generalized Boussinesq equation and a fourth-order Boussinesq type equation were solved with the aid of a Fourier–Galerkin spectral technique based on a complete orthonormal basis system in L2 ð1; 1Þ. An important application of the technique was presented here, where solitons with non-monotonic behavior were evaluated. The necessary details of the technique were discussed for the time-dependent case. A scaling parameter, was introduced allowing for the fine-tuning of the method and optimization of the convergence rate. Semi-implicit time-stepping schemes exhibiting very small phase error were constructed. The schemes were proved to be stable for relatively large time increments. The special structure of the matrices to be inverted made the phase errors for large s much smaller than the corresponding ones for finite-difference techniques. Numerical experiments with different numbers of terms were conducted demonstrating the exponential convergence rate of the method. Highly accurate results were obtained for the time-dependent problem with as few as N ¼ 80 terms. This shows that the above proposed technique is an efficient, robust and accurate numerical approach, suitable for the solution of time-dependent problems modeling soliton interactions. References [1] C.I. Christov, G.A. Maugin, M.G. Velarde, On the well-posed Boussinesq Paradigm with purely spatial higher-order derivatives, Phys. Rev. E 54 (1996) 3621–3638. [2] T.T. Marinov, C.I. Christov, R.S. Marinova, Novel numerical approach to solitary-wave solutions identification of Boussinesq and Korteweg–de Vries equations, Int. J. Bifurcation Chaos 15 (2005) 557–565. [3] C.I. Christov, T.T. Marinov, R.S. Marinova, Identification of solitary-wave solutions as an inverse problem: application to shapes with oscillatory tails, Math. Comput. Simul. 80 (2009) 56–65. [4] C.I. Christov, A complete orthonormal sequence of functions in L2 ð1; 1Þ space, SIAM J. Appl. Math. 42 (1982) 1337–1344. [5] J.P. Boyd, Chebyshev and Fourier Spectral Methods, Second ed. Revised, Springer-Verlag, New York, 2001. [6] J.P. Boyd, Z. Xu, Comparison of three spectral methods for the Benjamin–Ono equation: fourier pseudospectral, rational christov functions and gaussian radial basis functions, Wave Motion 48 (2011) 702–706. [7] J.P. Boyd, Z. Xu, Numerical and perturbative computations of solitary waves of the Benjamin–Ono equation with higher order nonlinearity using christov rational basis functions, J. Appl. Phys. 231 (2012) 1216–1229.

M.A. Christou, N.C. Papanicolaou / Applied Mathematics and Computation 243 (2014) 245–257

257

[8] T. Kawahara, Oscillatory solitary waves in dispersive media, J. Phys. Soc. Jpn. 33 (1972) 260–264. [9] A.M. Wazwaz, Soliton solutions for the fifth-order KdV equation and the Kawahara equation with time-dependent coefficients, Phys. Scr. 82 (2010) 035009. [10] A.M. Wazwaz, Non-integrable variants of Boussinesq equation with two solitons, Appl. Math. Comp. 217 (2010) 820–825. [11] A. Biswas, Soliton perturbation theory for the modified Kawahara equation, Appl. Appl. Math. 3 (2008) 218–223. [12] J.L. Bona, J.C. Saut, Dispersive blowup of solutions of generalized Korteweg–de Vries equations, J. Differ. Equ. 103 (1993) 3–57. [13] A.M. Wazwaz, Multiple-soliton solutions for the ninth-order KdV equation and the sixth-order Boussinesq equation, Appl. Math. Comp. 203 (2008) 277–283. [14] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamics, Springer-Verlag, New York, 1988. [15] C.I. Christov, K.L. Bekyarov, A Fourier-series method for solving soliton problems, SIAM J. Sci. Stat. Comput. 11 (1990) 631–647. [16] C.I. Christov, Fourier–Galerkin algorithm for 2D localized solutions, Annu. de l’Univ. Sofia Livre 2 – Math. Appl. Inf. 89 (1995) 169–179. [17] M.A. Christou, C.I. Christov, Fourier–Galerkin method for localized solutions of equations with cubic nonlinearity, J. Compu. Anal. Appl. 4 (1) (2002) 63– 77. [18] N. Wiener, Extrapolation, Interpolation and Smoothing of Stationary Time Series, Technology Press MIT and John Wiley, New York, 1949. [19] J.R. Higgins, Completeness and Basis Properties of Sets of Special Functions, Cambridge University Press, London, 1977. [20] K.L. Bekyarov, C.I. Christov, Fourier–Galerkin numerical technique for solitary waves of the fifth-order KdV equation, Chaos, Solitons Fractals 1 (1992) 423–430. [21] J.A.C. Weideman, The eigenvalues of Hermite and rational spectral differentiation matrices, Numer. Math. 61 (1992) 409–431. [22] C.I. Christov, An energy–consistent dispersive shallow-water model, Wave Motion 34 (2001) 161–174. [23] C.I. Christov, A method for identification of homoclinic trajectories, in: Proc. 14th Spring Conf. Union of Bulg. Mathematicians, Sofia, Bulgaria, 1985, pp. 571–577. [24] C.I. Christov, M.G. Velarde, On localized solutions of an equation governing Benard–Marangoni convection, Appl. Math. Model. 17 (1993) 311–320.