Analytic studies and numerical simulations of the generalized Boussinesq equation

Applied Mathematics and Computation 191 (2007) 320–333 www.elsevier.com/locate/amc

Analytic studies and numerical simulations of the generalized Boussinesq equation Mohamed Ali Hajji a

a,* ,

Kamel Al-Khaled

b

Department of Mathematical Sciences, United Arab Emirates University, P.O. Box 17551, Al-Ain, United Arab Emirates b Department of Mathematics, Jordan University of Science and Technology, IRBID 22110, Jordan

Abstract The modified Adomian decomposition method is used to solve the generalized Boussinesq equation. The equation commonly describes the propagation of small amplitude long waves in several physical contents. The analytic solution of the equation is obtained in the form of a convergent series with easily computable components. For comparison purposes, a numerical algorithm, based on Chebyshev polynomials, is developed and simulated. Numerical results show that the modified Adomian decomposition method proves to be more accurate and computationally more efficient than the Galerkin–Chebyshev method. Ó 2007 Published by Elsevier Inc. Keywords: Singularly perturbed Boussinesq equation; Decomposition method; Solitary wave solutions; Chebyshev polynomials; Numerical simulations; Boundary-value problems

1. Introduction The importance of soliton producing nonlinear wave equations is well understood among theoretical physicists and applied mathematicians. An equation admitting soliton solutions which has received comparatively little attention in the literature is utt ¼ uxx þ ðu2 Þxx þ uxxxx :

ð1:1Þ

It is referred to as the ‘‘bad’’ Boussinesq equation, or the nonlinear beam equation. It describes the motion of long waves in shallow water under gravity in a one-dimensional nonlinear lattices. Eq. (1.1) admits the solitary wave solution pffiffiffiffiffiffiffiffi uðx; tÞ ¼ A sech2 ð A=6ðx  ctÞÞ; ð1:2Þ

*

Corresponding author. E-mail addresses: [email protected] (M.A. Hajji), [email protected] (K. Al-Khaled).

0096-3003/$ - see front matter Ó 2007 Published by Elsevier Inc. doi:10.1016/j.amc.2007.02.090

M.A. Hajji, K. Al-Khaled / Applied Mathematics and Computation 191 (2007) 320–333

321

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where A and c ¼  1 þ 2A=3 are the amplitude and the speed of the solitary wave, respectively. These features of Eq. (1.1) are quite reminiscent of the properties of the Korteweg-de Vries (KdV) equation ut þ uux þ uxxx ¼ 0

ð1:3Þ

in that they both possess solitary wave solutions, except that the KdV equation allows only one-directional wave propagation and the Boussinesq equation describes bi-directional wave propagation. In recent years, a great deal of research has been conducted in the study of Eq. (1.1) from various points of view (see, for example, [6,8,10,12] and the references therein). For example, in [10] an exact formula is given for the interaction of solitary waves and in [8], Hirota has deduced conservation laws and has examined Nsoliton interaction. The representation of periodic waves as sums of solitons has been given by Whitham in [12]. Wazwaz in [11] used a modified algorithm of the Adomian decomposition method (shortly, ADM) to construct soliton solutions of Eq. (1.1) subject to the initial conditions uðx; 0Þ ¼ f ðxÞ;

ut ðx; 0Þ ¼ gðxÞ:

ð1:4Þ

El-Sayed and Kaya [6] studied the solitary-wave solutions, using the ADM, of the (2 + 1)-dimensional Boussinesq equation utt ¼ uxx þ uyy þ ðu2 Þxx þ uxxxx : In this paper, we consider a regularized version of Eq. (1.1) via the singularly perturbed (sixth-order) Boussinesq equation utt ¼ uxx þ ðpðuÞÞxx þ auxxxx þ buxxxxxx ;

ð1:5Þ

where a and b are real numbers (b is small). This equation was originally introduced by Daripa and Hua [4]. The sixth order derivative term provides dispersive regularization. The physical relevance of Eq. (1.5) in the context of water waves was recently addressed by Dash and Daripa [5]. It was shown that Eq. (1.5) actually describes the bi-directional propagation of small amplitude and long capillary-gravity waves on the surface of shallow water. So, it is closely related to the singularly perturbed (fifth-order) KdV equation ut þ uux þ uxxx þ 2 uxxxx ; which can be derived from Eq. (1.5) by using suitable transformations [5]. The fifth-order KdV equation has been studied by Kaya [9] where soliton solutions were found using the ADM. Since Eq. (1.1) has solitary wave solutions, the natural question arises whether Eq. (1.5) also admits solitary wave solutions for small values of b. As Eq. (1.5) can serve as a better model than the classical fourth-order Boussinesq equation (1.1)to describe bi-directional wave propagation on the surface of shallow water, in this paper we consider the generalized Boussinesq equation m X utt ¼ bi uð2iþ2Þx þ ½QðuÞxx ; ð1:6Þ i¼0

where QðuÞ ¼ u þ b0 ur ; r and bi ði ¼ 1; 2; . . . ; mÞ are all real constants and uð2iþ2Þx denotes the (2i + 2)nd derivative of u with respect to x. The modified ADM will be applied to seek soliton solutions of (1.6). Note that the choices m ¼ 1; b0 ¼ 1; b1 ¼ 1 and r = 2 yield Eq. (1.1) and for the choices m ¼ 2; b0 ¼ 1; b1 ¼ a; b2 ¼ b; m ¼ 2, and pðuÞ ¼ ur , Eq. (1.6) becomes the singularly perturbed sixth-order Boussinesq Eq. (1.5). The motivation of this paper is to approach the singularly perturbed Boussinesq equation (1.5) by using the modified ADM [11], the solutions will be calculated in the form of a convergent infinite series with easily computable components. Numerical results will illustrate the rapid convergence of the infinite series. In order to compare the modified ADM with other existing methods, a Galerkin method based on Chebyshev polynomials is proposed to numerically solve (1.6). The proposed Chebyshev–Galerkin method expresses the solution as a linear combination of Chebyshev polynomials with time dependent coefficients. Using the orthogonality of the Chebyshev polynomials, the partial differential equation is reduced to a coupled system of nonlinear second order ordinary differential equations for the time-dependent expansion coefficients. The second order system of ODEs is further written as a larger system of first-order ODEs which is solved numerically using Range–Kutta method of order 4.

322

M.A. Hajji, K. Al-Khaled / Applied Mathematics and Computation 191 (2007) 320–333

The paper is organized as follows. In Section 2, we describe the modified ADM as applied to the generalized Boussinesq equation (1.6). In Section 3, we implement the method on two examples of (1.6) for which an exact solution is known, thereby enabling us to evaluate its accuracy. In Section 4, we consider two linear boundaryvalue problems for which the ADM fails in certain situations. In Section 5, for comparison purposes, we propose a Chebyshev–Galerkin method to numerically solve the examples considered in Section 3. Discussions on the obtained results are presented in Section 6. Finally, the conclusion is presented in Section 7. 2. Analysis of the modified ADM In this section, we consider a general case of Eq. (1.5), namely, the generalized Boussinesq equation m X utt ¼ bi uð2iþ2Þx þ ½QðuÞxx ð2:1Þ i¼0

with the initial conditions (1.4), and present the analysis of the modified ADM. Following [2,3,6,11], we rewrite Eq. (2.1) in operator form as: LðuÞ ¼ RðuÞ þ N ðuÞ;

ð2:2Þ

where L and R are linear operators and N is a nonlinear operator. These operators are defined by m X bi uð2iþ2Þx ; N ðuÞ ¼ ½QðuÞxx : LðuÞ ¼ utt ; RðuÞ ¼ i¼0

We suppose that the linear operator L has an inverse linear operator L1 ðÞ ¼ both sides of Eq. (2.2) and using the initial conditions (1.4) we find

Rt Rt 0

0

ðÞ dt dt. Applying L1 to

uðx; tÞ ¼ f ðxÞ þ tgðxÞ þ L1 ½RðuÞ þ N ðuÞ:

ð2:3Þ

The Adomian decomposition method [2,3] assumes a series solution for the unknown function uðx; tÞ given by 1 X un ðx; tÞ ð2:4Þ uðx; tÞ ¼ n¼0

and the nonlinear operator N ðuÞ ¼ ½QðuÞxx can be decomposed into an infinite series of polynomials given by 1 X An ðu0 ; u1 ; . . . ; un Þ; ð2:5Þ N ðuÞ ¼ n¼0

where the components un ðx; tÞ in (2.4) will be determined recursively from the so-called Adomian polynomials An. The Adomian polynomials, An ; n P 0, can be constructed using the general formula [3], " !# n X 1 dn i An ¼ N k ui ; nP0 ð2:6Þ n! dkn i¼0 k¼0

and can be constructed for all classes of nonlinearity according to algorithms set by Adomian [3]. Substitution of (2.4) and (2.5) into (2.3) yields 2 3 ! Z tZ t X 1 m 1 1 X X X 4 un ¼ f ðxÞ þ tgðxÞ þ bi un þ An 5 dt dt: 0

n¼0

0

i¼0

n¼0

ð2iþ2Þx

ð2:7Þ

n¼0

Following the modified decomposition method, setting the zeroth and first components u0 ðx; tÞ ¼ f ðxÞ; u1 ðx; tÞ ¼ tgðxÞ þ

Z 0

t

Z t "X m 0

i¼0

# ðu0 Þð2iþ2Þx þ A0 dt dt;

ð2:8Þ ð2:9Þ

M.A. Hajji, K. Al-Khaled / Applied Mathematics and Computation 191 (2007) 320–333

323

we obtain the recurrence relation defining the solution components, un ðx; tÞ; n P 2, # Z t Z t "X m unþ1 ðx; tÞ ¼ ðun Þð2iþ2Þx þ An dt dt; n P 1: 0

0

i¼0

Finally and N-term approximate solution is given by /N ðx; tÞ ¼

N 1 X

un ðx; tÞ

n¼0

and the exact solution is uðx; tÞ ¼ lim /N ðx; tÞ: N !1

To show the effectiveness of the proposed ADM and to get a clear overview of methodology, some examples of the generalized Boussinesq equation (1.6) will be discussed in the following section. The convergence of the Adomian decomposition series has been investigated by several authors (see, for example, [1] and the references therein). In the present work, we demonstrate the fast convergence of the method numerically by applying it to examples with known exact solutions and show that the N-term approximate solution is very close to the exact solution. 3. Implementation of the method To illustrate the accuracy and efficiency of the proposed modified ADM in studying (1.6), we consider the singularly perturbed Boussinesq model equation (1.5). In particular, we consider Eq. (1.5) for sufficiently small values of b and apply the technique presented in the previous section to study its soliton solutions where pðuÞ ¼ cu2 ; c constant. Writing Eq. (1.5) in an operator form as in (2.2), we identify the operators L; R and N to be LðuÞ ¼ utt ;

N ðuÞ ¼ cðu2 Þxx :

RðuÞ ¼ uxx þ buxxxx þ auxxxxxx ;

Therefore, Eq. (2.7) corresponding to Eq. (1.5) takes the form ! ! Z t Z t" X 1 1 1 X X un ¼ f ðxÞ þ tgðxÞ þ un þa un 0

n¼0

0

n¼0

xx

n¼0

xxxx

þb

1 X

! þc

un

n¼0

xxxxxx

1 X

# An dt dt:

n¼0

Then, following the modified ADM, the recurrence relation is determined to be: u0 ðx; tÞ ¼ f ðxÞ;

Z tZ t u1 ðx; tÞ ¼ tgðxÞ þ ½ðu0 Þxx þ aðu0 Þxxxx þ bðu0 Þxxxxxx þ cA0  dt dt; 0 0 Z tZ t unþ1 ðx; tÞ ¼ ½ðun Þxx þ aðun Þxxxx þ bðun Þxxxxxx þ cAn  dt dt; n P 1; 0

ð3:1Þ

0

where in this case the An are the Adomian polynomials representing the nonlinear term ðu2 Þxx . The first few An, using (2.6), are calculated to be A0 ¼ ðu20 Þxx ;

A1 ¼ 2u0 u1xx þ 4u0x u1x þ 2u0xx u1 ;

A2 ¼ 2u0 u2xx þ 4u0x u2x þ 2u0xx u2 þ 2u1 u1xx þ 2ðu1x Þ2 ; A3 ¼ ð2u0 u3 þ 2u1 u2 Þxx : In the remainder of this section, we consider two examples of (1.5) for which exact solutions are known. Example 3.1. In this example we consider (1.5) with a ¼ 1; b ¼ 0 and pðuÞ ¼ 3u2 . The model equation with the initial conditions are given as

324

M.A. Hajji, K. Al-Khaled / Applied Mathematics and Computation 191 (2007) 320–333

utt ¼ uxx þ 3ðu2 Þxx þ uxxxx ; 2 kx

2ak e

uðx; 0Þ ¼

ð1 þ aekx Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffi 2ak 1 þ k 2 ð1  aekx Þekx

ð3:2Þ

3

2

ut ðx; 0Þ ¼

;

ð1 þ aekx Þ

3

ð3:3Þ

;

where a and k are arbitrary constants. To find the solution of Eq. (3.2) subject to the conditions in (3.3), we apply the outlined ADM in the previous section, in the same manner as in (3.1), we obtain u0 ðx; tÞ ¼ uðx; 0Þ;

Z

t

Z

t

u1 ðx; tÞ ¼ tut ðx; 0Þ þ ½ðu0 Þxx þ ðu0 Þxxxx þ 3A0  dt dt; 0 0 Z tZ t unþ1 ðx; tÞ ¼ ½ðun Þxx þ ðun Þxxxx þ 3An  dt dt; n P 1: 0

0

Using the software package Mathematica, with a = k = 1, we calculate the first few solution components to be u0 ðx; tÞ ¼

2ex

; ð1 þ ex Þ2 pffiffiffi 2 2ex ð1 þ ex Þ

2ex ð1  4ex þ e2x Þ

t2 ; 3 4 ð1 þ ex Þ ð1 þ ex Þ pffiffiffi 2 2ex ð1 þ ex Þð1  10ex þ e2x Þ 3 u2 ðx; tÞ ¼ t 5 3ð1 þ ex Þ ex ð1  4ex þ e2x Þð1  44ex þ 78e2x  44e3x þ e4x Þ 4 þ t 8 3ð1 þ ex Þ

u1 ðx; tÞ ¼

tþ

and u3 ðx; tÞ ¼

8e2x ð1  10ex þ 20e2x  10e3x þ e4x Þ ð1 þ ex Þ

þe þ

8

4

t þ

pffiffiffi  2ex ð1 þ ex Þð1  56ex þ 246e2x  56e3x þ e4x Þ 15ð1 þ ex Þ7

t5

" x 2x 3x 4x 5x x ð1  452e þ 19149e  207936e þ 807378e  1256568e Þ 45ð1 þ ex Þ12

ð807378e6x  207936e7x þ 19149e8x  452e9x þ e10x Þ 45ð1 þ ex Þ

12

# t6 :

More solution components, un, can be easily generated with the help of Mathematica. Using, the first four computed components, the approximate solution ua ðx; tÞ, in a series form, can be written as uðx; tÞ  u0 ðx; tÞ þ u1 ðx; tÞ þ u2 ðx; tÞ þ u3 ðx; tÞ: The exact solution ue ðx; tÞ of Eq. (3.2) subject to (3.3) is given by [11] pffiffiffiffiffiffiffiffiffiffiffiffiffi ak 2 expðkx þ k 1 þ k 2 tÞ pffiffiffiffiffiffiffiffiffiffiffiffiffi : ue ðx; tÞ ¼ 2 2 ð1 þ a expðkx þ k 1 þ k 2 tÞÞ

ð3:4Þ

The absolute error between the 4-term approximate solution and the exact solution is reported in Table 1. Example 3.2. In this example, we consider (1.5) with a ¼ 1; b ¼ 12, and pðuÞ ¼ u2 . The model equation is 1 utt ¼ uxx þ ðu2 Þxx  uxxxx þ uxxxxxx 2 subject to the initial conditions

ð3:5Þ

M.A. Hajji, K. Al-Khaled / Applied Mathematics and Computation 191 (2007) 320–333

325

Table 1 The error kue  ua k in the approximate solution for Example 3.1 using Modified ADM xi

tj

1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1

0.01

0.02

0.04

0.1

0.2

0.5

2.80886E14 6.27276E14 6.08402E14 1.16573E14 5.53446E14 8.63198E14 5.56222E14 1.14353E14 6.06182E14 6.23945E14 2.79776E14

1.79667 E12 4.01362E12 3.90188E12 7.41129E13 3.53395E12 5.53357E12 3.55044E12 7.14928E13 3.87551E12 3.99519E12 1.78946E12

1.15235E10 2.57471E10 2.25663E10 4.82726E11 2.25663E10 3.54174E10 2.27779E10 4.49107E11 2.47218E10 2.55127E10 1.14307E10

2.83355E8 6.33178E8 6.18024E8 1.23843E8 5.47485E8 8.65197E8 5.60362E8 1.03370E8 5.97562E8 6.18881E8 2.77684E8

1.83899E6 4.10454E6 4.02299E6 8.53800E7 3.47264E6 5.54893E6 3.63600E6 5.93842E7 3.76275E6 3.92220E6 1.76607E6

4.74681E4 1.04489E3 1.03093E3 2.46302E4 8.35783E4 1.37353E3 9.29612E4 9.61260E5 8.79002E4 9.36404E4 4.28986E4

uðx; 0Þ ¼ 

105 sech4 169



 x pffiffiffiffiffi ; 26

210 ut ðx; 0Þ ¼

qffiffiffiffiffi     4 pxffiffiffiffi 194 pxffiffiffiffi sech tanh 13 26 26 2197

:

ð3:6Þ

It is known that Eq. (3.5) subject to (3.6) has an exact solution (see [7]) of the form 105 ue ðx; tÞ ¼  sech4 169

"rffiffiffiffiffi rffiffiffiffiffiffiffiffi # 1 97 ðx  tÞ : 26 169

ð3:7Þ

Following the discussions presented above, and using Eq. (3.1), we find:   105 x 4 p ffiffiffiffiffi sech u0 ðx; tÞ ¼  ; 169 26 rffiffiffiffiffiffiffiffi rffiffiffiffiffi !   105 194 x 2 6 sech pffiffiffiffiffi sinh x t u1 ðx; tÞ ¼ 2197 13 13 26 " rffiffiffiffiffi !#   105 2 x 6 291 þ 194 cosh x sech pffiffiffiffiffi t2 ;  371293 13 26     7 pxffiffiffiffi  pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi 3395sech ð 26Þ x 3x 10816 2522 sinh pffiffiffiffiffi  1664 2522 sinh pffiffiffiffiffi t3 u2 ðx; tÞ ¼ 52206766144 26 26 rffiffiffiffiffi !     x 2 x þ 334200sech5 pffiffiffiffiffi þ 354247 cosh x sech5 pffiffiffiffiffi 13 26 26 rffiffiffiffiffi ! rffiffiffiffiffi !     2 x 2 x  47164 cosh 2 x sech5 pffiffiffiffiffi þ 3201 cosh 3 x sech5 pffiffiffiffiffi 13 13 26 26 rffiffiffiffiffi !  ! # 2 x 388 cosh 4 x sech5 pffiffiffiffiffi t4 ; 13 26 and, similarly, we can find more components un. Table 2, reports the absolute error between the 4-term approximate solution and the exact solution. 4. Two examples where the ADM may fail In this section, we consider two interesting linear boundary-value problems where the Adomian decomposition method may fail in certain situations. The two ODE examples we consider in this section were

326

M.A. Hajji, K. Al-Khaled / Applied Mathematics and Computation 191 (2007) 320–333

Table 2 The error kue  ua k in the approximate solution for Example 3.2 using Modified ADM xi

tj

1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1

0.01

0.02

0.04

0.1

0.2

0.5

7.77156E16 1.11022E16 2.22045E16 1.11022E16 6.66134E16 4.44089E16 5.55112E16 3.33067E16 3.33067E16 3.33067E16 7.77156E16

1.36557E14 1.99840E15 1.09912E14 2.32037E14 3.23075E14 3.49720E14 3.19744E14 2.32037E14 1.12133E14 1.99840E15 1.38778E14

8.57869E13 1.12688E13 7.28861E13 1.50302E12 2.04747E12 2.24365E12 2.04714E12 1.50324E12 7.28528E13 1.13132E13 8.58313E13

2.09264E10 2.73880E11 1.78030E10 3.67002E10 4.99918E10 5.47741E10 4.99820E10 3.66815E10 1.77772E10 2.76944E11 2.09593E10

1.33823E8 1.74288E9 1.14025E8 2.34944E8 3.19983E9 3.50559E8 3.19858E8 2.34706E8 1.13695E8 1.78208E9 1.34244E8

3.25944E6 4.14094E7 2.79028E6 5.74091E6 7.81509E6 8.55935E6 7.80749E6 5.72641E6 2.77022E6 4.41936E7 3.28504E6

considered and suggested by Scott [13]. In fact in [13,14], Scott et al. propose a computer code to solve accurately many two-point boundary-value problems to which other numerical methods may fail. In our study, we shall investigate the reasons behind the failure of the Adomian decomposition method. Example 4.3. The first example we consider is the fourth-order boundary-value problem 1 y 000  ð1 þ cÞy 00 þ cy ¼ cx2  1; 2 with the boundary conditions yð0Þ ¼ 1;

y 0 ð0Þ ¼ 1;

yð1Þ ¼

ð4:1Þ

0 6 x 6 1;

3 þ sinhð1Þ; 2

y 0 ð1Þ ¼ 1 þ coshð1Þ:

ð4:2Þ

It can be easily verified that the exact solution of (4.1) subject to (4.2) is 1 yðxÞ ¼ 1 þ x2 þ sinhðxÞ: 2

ð4:3Þ

It is interesting to note that the exact solution in (4.3) is independent of the parameter c. However, any approximate solution of (4.1) will not be independent of c. This means that the values of c affect greatly the approximate solution. First, we shall apply the ADM and show that the series solution obtained converges to the solution (4.3). Integrating (4.1) we obtain yðxÞ ¼

6 X

ai xi þ ð1 þ cÞI 2 ðyðxÞÞ  cI 4 ðyðxÞÞ;

i¼0

where 1 a2 ¼ ½y 00 ð0Þ  ð1 þ cÞyð0Þ; 2 1 000 1 1 a3 ¼ ½y ð0Þ  ð1 þ cÞy 0 ð0Þ; a4 ¼  ; a5 ¼ 0; a6 ¼ c; 3! 4! 6! Rx Rs and Ik is a k-fold integral, for example I 2 ðyðxÞÞ ¼ 0 0 yðtÞ dt ds. Note that the unknowns are y 00 ð0Þ and y 000 ð0Þ which define a2 and a3. These will be determined later by satisfying the conditions at x =P 1. Pboundary 1 6 Following the ADM procedure, assuming a series expansion yðxÞ ¼ n¼0 y n ðxÞ and letting y 0 ðxÞ ¼ i¼0 ai xi , we find the recursion relation a0 ¼ yð0Þ ¼ 1;

a1 ¼ y 0 ð0Þ ¼ 1;

y n ¼ ½ð1 þ cÞI 2  I 4 ðy n1 Þ;

n P 1:

M.A. Hajji, K. Al-Khaled / Applied Mathematics and Computation 191 (2007) 320–333

Solving for yn in terms of y0, we find y n ðxÞ ¼ ½ð1 þ cÞI 2  cI 4 n ðy 0 ðxÞÞ ¼

327

n   X n ð1Þk ð1 þ cÞnk ck I 2nþ2k ðy 0 ðxÞÞ: k k¼0

We have I 2nþ2k ðy 0 ðxÞÞ ¼

6 X

ai i!

i¼0

xiþ2nþ2k : ði þ 2n þ 2kÞ!

It follows that y n ðxÞ ¼

6 X i¼0

ai i!xi

n   X n k¼0

k

k

ð1Þ ð1 þ cÞ

nk k

c

x2nþ2k : ði þ 2n þ 2kÞ!

Finally, we get the solution yðxÞ formally as ! 6 1 X n   X X n x2nþ2k k nk k i yðxÞ ¼ ai i!x ð1Þ ð1 þ cÞ c : ði þ 2n þ 2kÞ! k i¼0 n¼0 k¼0 The two inner summations can be written as 0 n 1  b 2c  1 X X nk x2n @ : ð1Þk ð1 þ cÞn2k ck A ði þ 2nÞ! k k¼0 n¼0

ð4:4Þ

The inner finite sum in (4.4) simplifies to  bn2c  n X X nk 1  cnþ1 k n2k k : c ¼ ck ¼ ð1Þ ð1 þ cÞ 1c k k¼0 k¼0 The solution yðxÞ is then given by yðxÞ ¼

6 X i¼0

! 1 X 1  cnþ1 x2nþi ai i! ; 1  c ð2n þ iÞ! n¼0

ð4:5Þ

which is clearly a convergent power series. Next, we express the above power series solution in terms of cos h and sin h functions as follows: ! " pffiffiffi 2nþi # 6 1 6 1 1 X X X X 1  cnþ1 x2nþi 1 X x2nþi ðx cÞ 1i=2 c yðxÞ ¼ ai i! ai i! ¼ 1  c i¼0 1  c ð2n þ iÞ! ð2n þ iÞ! ð2n þ iÞ! i¼0 n¼0 n¼0 n¼0 pffiffiffi pffiffiffi 1 ½a0 ðcoshðxÞ  c coshð cxÞÞ þ a1 ðsinhðxÞ  c1=2 sinhð cxÞÞ þ 2!a2 ðcoshðxÞ  1 1c pffiffiffi pffiffiffi pffiffiffi  c0 ðcoshð cxÞ  1ÞÞ þ 3!a3 ðsinhðxÞ  x  c1=2 ðsinhð cxÞ  cxÞÞ þ 4!a4 ðcoshðxÞ  1  x2 =2 pffiffiffi pffiffiffi 2 pffiffiffi pffiffiffi  c1 ðcoshð cxÞ  1  ð cxÞ =2ÞÞ þ 5!a5 ðsinhðxÞ  x  x3 =3!  c3=2 ðsinhð cxÞ  cx pffiffiffi pffiffiffi 2 pffiffiffi 4 pffiffiffi 3  ð cxÞ =3!ÞÞ þ 6!a6 ðcoshðxÞ  1  x2 =2  x4 =4!  c2 ðcoshð cxÞ  1  ð cxÞ =2  ð cxÞ =4!ÞÞ pffiffiffi 1 ½ða0 þ 2!a2 þ 4!a4 þ 6!a6 Þ coshðxÞ  ða0 c þ 2!a2 c0 þ 4!a4 c1 þ 6!a6 c2 Þ coshð cxÞ ¼ 1c pffiffiffi þ ða1 þ 3!a3 þ 5!a5 ÞsinhðxÞ  ða1 c1=2 þ 3!a3 c1=2 þ 5!a5 c3=2 Þ sinhð cxÞ þ ð1  cÞð1 þ x2 =2Þ:

¼

Using the values of ai ; 0 6 i 6 6, the above simplifies to pffiffiffi pffiffiffi 1  00 ðy 000 ð0Þ  1Þ pffiffiffi ðy ð0Þ  1Þ coshðxÞ  ðy 00 ð0Þ  1Þ coshð cxÞðy 000 ð0Þ  cÞsinhðxÞ  yðxÞ ¼ sinhð cxÞ 1c c þ ð1 þ x2 =2Þ: ð4:6Þ

328

M.A. Hajji, K. Al-Khaled / Applied Mathematics and Computation 191 (2007) 320–333

Now applying the two boundary conditions yð1Þ ¼ 3=2 þ sinhð1Þ, y 0 ð1Þ ¼ 1 þ coshð1Þ, we obtain the following system for the unknown y 00 ð0Þ and y 000 ð0Þ: pffiffiffi pffiffiffi pffiffiffi

a11 y 00 ð0Þ þ a12 y 000 ð0Þ ¼ sinhð1Þ þ coshð1Þ  coshð cÞ  sinhð cÞ= c pffiffiffi pffiffiffi pffiffiffi a21 y 00 ð0Þ þ a22 y 000 ð0Þ ¼ sinhð1Þ þ coshð1Þ  coshð cÞ  c sinhð cÞ; where pffiffiffi a11 ¼ coshð1Þ  coshð cÞ; pffiffiffi pffiffiffi a21 ¼ sinhð1Þ  c sinhð cÞ;

pffiffiffi pffiffiffi a12 ¼ sinhð1Þ  sinhð cÞ= c; pffiffiffi a12 ¼ coshð1Þ  coshð cÞ:

The above system has the unique solution y00ð0Þ ¼ y000ð0Þ ¼ 1, for any value of c, which is when substituted in (4.6) give the solution 1 yðxÞ ¼ 1 þ x2 þ sinhðxÞ; 2 which is the exact solution as given in (4.3). It is important to note that even though the ADM procedure was successful in obtaining the exact solution for the above example, the procedure will fail to produce accurate approximate solution for certain values of c. Actually the ADM does not always lead to a closed form exact solution, but one usually constructs an approximate solution by adding up a finite number of solution components y n ðxÞ; 0 6 n 6 N , where N is the truncation level of the series. Our simulations indicate that the ADM does not lead to accurate approximate solution for certain value of c. The reason for this the following. When applying the boundary conditions at x = 1 we obtain a system of two equations for the unknowns y 00 ð0Þ and y 000 ð0Þ which is ill-conditioned for large values of the parameter c. In fact, the condition number of the system increases with increasing c. We have used the software package Mathematica to obtain an approximate solution using the values of c ¼ 1; 10; 100; 1000; 3600; 10 000; 106 ; 108 . We found that for c ¼ 1000; 3600; 10 000; 106 ; and 108, the system to be solved for y 00 ð0Þ and y 000 ð0Þ is ill-conditioned, for any truncation level. However, for c = 1 and c = 10, the system was reasonably well-conditioned and we were able to find an accurate approximate solution to the order of 1010 with a truncation level N = 10. For c = 100 a truncation level N = 20 was needed to obtain and approximate solution accurate to 1010. Example 4.4. In this example, we consider the eight-order boundary-value problem y ð8Þ þ c1 y ð6Þ þ c2 y ð4Þ þ c3 y ð2Þ þ c4 y ¼ 0;

ð4:7Þ

0 6 x 6 L;

with the boundary conditions y ðiÞ ð0Þ ¼ 0; 0 6 i 6 3;

y ðiÞ ðLÞ ¼ sðiÞ ðLÞ; 0 6 i 6 3;

ð4:8Þ

where c1 ¼ 914; c2 ¼ 12649; c3 ¼ 44136; c4 ¼ 32400, and sðxÞ ¼ ex  2e2x þ e3x is the exact solution of (4.7) subject to (4.8). Following the ADM procedure as in the previous example, we find yðxÞ ¼

7 X

ai xi  ðc1 I 2 þ c2 I 4 þ c3 I 6 þ c4 I 8 ÞyðxÞ;

i¼0

where, as before, Ik is a k-fold integral, and y ð4Þ ð0Þ y ð5Þ ð0Þ y ð6Þ ð0Þ þ c1 y ð4Þ ð0Þ y ð7Þ ð0Þ þ c1 y ð5Þ ð0Þ ; a5 ¼ ; a6 ¼ ; a7 ¼ : ai ¼ 0; 0 6 i 6 3; a4 ¼ 4! 5! 6! 7! P1 P7 Letting yðxÞ ¼ n¼0 y n ðxÞ and the zeroth component y 0 ðxÞ ¼ i¼0 ai xi , we find that the nth component y n ðxÞ is given by n

n

y n ðxÞ ¼ ð1Þ ðc1 I 2 þ c2 I 4 þ c3 I 6 þ c4 I 8 Þ y 0 ðxÞ   n X k X nk   X n k nk n ¼ ð1Þ cl3 c4nkl ; c c cm1 ckm 2 k m l k¼0 m¼0 l¼0

I 8n4k2m2l ðy 0 ðxÞÞ;

M.A. Hajji, K. Al-Khaled / Applied Mathematics and Computation 191 (2007) 320–333

329

where I 8n4k2m2l ðy 0 ðxÞÞ ¼

7 X

i!ai

i¼0

xiþ8n4k2m2l : ði þ 8n  4k  2m  2lÞ!

Then the solution yðxÞ is given formally by the power series yðxÞ ¼

7 X

i!ai x

i

i¼0

!     n k n  k cm1 ckm cl3 c4nkl x8n4k2m2l 2 ð1Þ ; ði þ 8n  4k  2m  2lÞ! k m l l¼0

1 X n X k X nk X n¼0 k¼0 m¼0

n

which can be shown to be convergent. An approximate solution, with a truncation level N, is y approx ðxÞ ¼

N X

y n ðxÞ:

n¼0

Our numerical simulations on this example reveal that for any truncation level N, we were not able to get an accurate approximate solution. For moderate N, the system of linear equations to be solved for the unknowns y ðiÞ ð0Þ; 4 6 i 6 7, was reasonably well-conditioned but this led to inaccurate approximate solution. For high truncation level N, the system became ill-conditioned. The above two examples suggest that the Adomian decomposition method is not always suitable for certain boundary value problems. In Example 4.3, the series solution converges to the exact solution, but an accurate approximate solution was not possible for certain values of the parameter c. In Example 4.4, the series solution is convergent, but again an accurate approximate solution was not possible. The reason for this drawback in the ADM is due to the ill-conditioning of the system defining the unknowns. 5. A Chebyshev–Galerkin solution In this section, we formulate a Galerkin approach based on Chebyshev polynomials that may be applied to solve numerically Eq. (1.5). We illustrate the method by solving the Boussinesq equation utt ¼ uxx þ uxxxx þ 3ðu2 Þxx ;

t P 0; 1 6 x 6 1;

ð5:1Þ

subject to the initial conditions 2ex

; 1 6 x 6 1; 2 ð1 þ ex Þ pffiffiffi 2 2ex ð1  ex Þ ut ðx; 0Þ ¼ gðxÞ ¼  ; 1 6 x 6 1: 3 ð1 þ ex Þ

uðx; 0Þ ¼ f ðxÞ ¼

Eq. (5.1) subject to (5.2) and (5.3) admits an exact travelling wave solution pffiffi 2exþ 2t pffiffi ue ðx; tÞ ¼ : 2 ð1 þ exþ 2t Þ

ð5:2Þ ð5:3Þ

ð5:4Þ

The Chebyshev–Galerkin method is based on approximating the solution uðx; tÞ by a linear combination of Chebyshev polynomials: uðx; tÞ ¼

N X

y i ðtÞT i ðxÞ

ð5:5Þ

i¼0

with time-dependent coefficients y i ðtÞ, where Tis are the monic Chebyshev polynomials defined, recursively, by

330

M.A. Hajji, K. Al-Khaled / Applied Mathematics and Computation 191 (2007) 320–333

1 T 0 ðxÞ ¼ 1; T 1 ðxÞ ¼ x; T 2 ðxÞ ¼ x2  ; and 2 1 T iþ1 ðxÞ ¼ xT i ðxÞ  T i1 ðxÞ for each i P 2: 4 It is known that these polynomials are orthogonal, with respect to the weight function 1 wðxÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi : 1  x2 Substituting (5.5) into (5.1), we obtain N X

y 00i ðtÞT i ðxÞ ¼

i¼0

N X

0000

y i ðtÞðT 00i ðxÞ þ T i ðxÞÞ þ 3

i¼0

N X

ðy i ðtÞy j ðtÞÞðT i ðxÞT j ðxÞÞ

00

i;j¼0

and multiplying by a test function T m ðxÞ and integrating from 1 to 1, exploiting the orthogonality of Ti, we arrive at the nonlinear system of second order ODEs for y m ðtÞ, y 00m ðtÞ ¼

N X

bmi y i ðtÞ þ

i¼0

N X

cmij y i ðtÞy j ðtÞ;

0 6 m 6 N;

ð5:6Þ

i;j¼0

where bmi ¼

1 am

Z

1

1

0000

wðxÞT m ðxÞðT 00i ðxÞ þ T i ðxÞÞ dx;

Z 1 3 00 wðxÞT m ðxÞðT i ðxÞT j ðxÞÞ dx; cmij ¼ am 1 Z 1 wðxÞT 2m ðxÞ dx: am ¼ 1

To solve system (5.6), we first rewrite it as coupled system of first-order ODEs. Let y 1;m ¼ y m and y 2;m ¼ y 0m . Then y 00m ¼ y 02;m and system (5.6) becomes 8 0 y ¼ y 2;m ; > > < 1;m N N > y 0 ¼ P b y ðtÞ þ P c y ðtÞy ðtÞ; > mi 1;i mij 1;i 1;j : 2;m i¼0

0 6 m 6 N;

ð5:7Þ

i;j¼0

which is now a nonlinear system of 2ðN þ 1Þ first order ODEs. We propose to use Range–Kutta method of ðsÞ ðsÞ order 4 to solve (5.7). Let z1;m  y 1;m ðshÞ; z2;m  y 2;m ðshÞ, where h is the time step. We then have the explicit four-order difference equations 1 ð1Þ ðsþ1Þ ðsÞ ð2Þ ð3Þ ð4Þ z1;m ¼ z1;m þ ðK 1;m þ 2K 1;m þ 2K 1;m þ K 1;m Þ; 6 ðsþ1Þ z2;m

¼

ðsÞ z2;m

1 ð1Þ ð2Þ ð3Þ ð4Þ þ ðK 2;m þ 2K 2;m þ 2K 2;m þ K 2;m Þ; 6

where the Ks are given by

ð5:8Þ

M.A. Hajji, K. Al-Khaled / Applied Mathematics and Computation 191 (2007) 320–333

331

Table 3 The error kue  ua k in the approximate solution for Example 3.1 using Chebyshev–Galerkin xi

ti

1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1

0.01

0.02

0.04

0.1

0.2

0.5

2.4017E6 1.99012E7 1.16573E7 8.82439E8 2.58332E8 9.55971E8 6.5813E8 6.61518E8 1.40176E7 2.57542E7 2.61255E 6

9.22269E6 6.94226E7 2.03366E7 2.41799E7 2.16049E8 2.11864E7 1.27462E7 1.91194E7 2.91472E7 9.80217E7 1.0699E5

3.61065E5 4.2212E6 7.26982E7 1.28865E6 1.04947E7 9.31864E7 3.54842E7 1.21154E6 3.13147E7 6.2021E6 4.74423E5

2.65583E4 1.00858E4 5.66502E5 2.43507E5 7.61264E6 1.72816E5 2.27783E6 3.27314E5 6.12082E5 1.36886E4 4.35211E4

1.86075E3 1.27366E3 7.36563E4 2.0634E4 1.64814E4 2.30693E4 2.4065E5 4.85432E4 1.03433E3 1.74915E3 3.11884E3

1.16237E2 8.41167E3 1.54798E3 4.85949E3 8.10358E3 6.82062E3 8.61072E4 8.92505E3 2.09809E2 3.3402E2 4.43529E2

ð1Þ

ðsÞ

K 1;m ¼ hz2;m ; ð1Þ K 2;m

¼h

N X

ðsÞ bmi z1;i

þ

i¼0 ð2Þ

ðsÞ

N X

! ðsÞ ðsÞ cmij z1;i z1;j

;

i;j¼0 ð1Þ

K 1;m ¼ hðz2;m þ K 2;m =2Þ; ð2Þ K 2;m

¼h

N X

ðsÞ bmi ½z1;i

þ

ð1Þ K 1;i =2

þ

i¼0 ð3Þ

ðsÞ

N X

! ðsÞ cmij ½z1;i

þ

ð1Þ ðsÞ K 1;i =2½z1;j

þ

ð2Þ ðsÞ K 1;i =2½z1;j

þ

ð1Þ K 1;j =2

þ

ð2Þ K 1;j =2

ð2Þ

K 1;m ¼ hðz2;m þ K 2;m =2Þ; ð3Þ K 2;m

¼h

N X

ðsÞ bmi ½z1;i

þ

ð2Þ K 1;i =2

i¼0 ð4Þ

ðsÞ

þ

N X

! ðsÞ cmij ½z1;i

¼h

N X

;

i;j¼0 ð3Þ

K 1;m ¼ hðz2;m þ K 2;m Þ; ð4Þ K 2;m

;

i;j¼0

ðsÞ bmi ½z1;i

þ

ð3Þ K 1;i 

þ

i¼0

N X

! ðsÞ cmij ½z1;i

þ

ð3Þ ðsÞ K 1;i ½z1;j

þ

ð3Þ K 1;j 

:

i;j¼0

The initial conditions for system (5.8) are ð0Þ

z1;m ¼

hf ðxÞ; Lm ðxÞiw ; hLm ðxÞ; Lm ðxÞiw

ð0Þ

z2;m ¼

hgðxÞ; Lm ðxÞiw ; hLm ðxÞ; Lm ðxÞiw

where the inner product h; iw is given by Z 1 hf ; giw ¼ f ðxÞgðxÞwðxÞ dx: 1

Numerical simulation of the above outlined Galerkin method has been carried out to solve numerically (5.1) subject to (5.2) and (5.3) using the parameter values a ¼ k ¼ 1; N ¼ 8; h ¼ 0:01. The accuracy of the method has been evaluated by calculating the absolute error between the approximate and the exact solution at various time steps and various spacial nodes xi. The results are reported in Table 3. 6. Numerical evaluation and discussions In Section 3 the modified ADM has been applied to approximate the solution of the fourth order Boussinesq equation (3.2) and of the singularly perturbed sixth order Boussinesq equation (3.5). The approximate P4 solution was obtained by adding up only four Adomian solution components; i.e., ua ðx; tÞ ¼ i¼0 ui ðx; tÞ.

332

M.A. Hajji, K. Al-Khaled / Applied Mathematics and Computation 191 (2007) 320–333

Tables 1 and 2 clearly show that the modified ADM gave excellent results in that the errors at various spacial and temporal nodes are very small. It is worth noting that the error can be made smaller simply by adding up more solution components. It can be seen as well from Table 2, that the errors in the case of (3.5) are of smaller order than in the case of (3.2). This can be attributed to the nature of the solution of (3.5), where we need fewer Adomian solution components to accurately approximate the solution. It is important to note that this method does not require any discretization neither of the spacial nor of the time variable. For comparison purposes, in Section 5, we have applied a Galerkin method based on monic Chebyshev polynomials to solve numerically the (3.2) in the interval ½1; 1. Runge-Kutta method of order 4 was used to solve the obtained system of ODES. As can be seen from the results shown in Table 3, the Galerkin method has limited accuracy; the accuracy depends on the time step size, h. In fact, this limitation is experienced by many numerical schemes. 7. Conclusion In this paper, the modified ADM was used to solve the generalized Boussinesq equation (2.1) subject to initial conditions, where we considered the case r ¼ 2; m ¼ 2. It is worthwhile to point out here that r and m may be any arbitrary integers, and so our method can be applied to some other nonlinear differential equations. The solution was approximated by adding up few Adomian solution components produced by the modified ADM. Very good accuracy was achieved with only 4 terms. The Galerkin method, on the other hand, had limited accuracy, due to time discretization. A clear conclusion can be draw from the numerical results is that the ADM provides highly accurate numerical solutions for certain nonlinear partial differential equations. It is also worth noting that the Adomian decomposition methodology, for ceratin problems, displays a fast convergence of the solution. On the other hand (see Section 4), the ADM may not be suitable for certain boundary-value problems. In particular, in Example 4.4, the ADM fails to produce an accurate approximate solutions. However, in Example 4.3, accurate approximate solution was only possible for limit values of the parameter c. This drawback in the ADM is attributed to the large coefficients in the equation. Precisely, using the ADM one has to solve a system of linear equations and this system becomes ill-conditioned if the coefficients are too large. For these kinds of boundary-value problems one might use other numerical algorithms such as the ones discussed in [13–17]. This present work affirms that the modified ADM is an easy straight forward method to solve nonlinear partial differential equations, but may not be suitable for certain boundary-value problems. Finally, we point out that the ADM can be generalized to two and higher-dimensional problems. For example, one might consider the two-dimensional generalized Boussinesq equation utt ¼ ½QðuÞxx þ ½RðuÞyy þ

m X i¼1

bi uð2iþ2Þx þ

k X

cj uð2jþ2Þy ;

j¼1

where RðuÞ ¼ us ; and s; r; bi ; cj are all real numbers. Acknowledgement The authors are grateful to Dr. M. Scott for his interesting remarks and suggestions. References [1] K. Abbaoui, M.J. Pujol, Y. Cherruault, N. Himoun, P. Grimalt, A new formulation of adomian method: convergence result, Kybernetes 30 (2001) 1183–1191. [2] G. Adomian, Solving Frontier Problems of Physics, The Decomposition Method, Kluwer Academic Publishers, Boston, MA, 1994. [3] G. Adomian, A review of the decomposition method in applied mathematics, J. Math. Anal. Appl. 135 (1988) 501–544. [4] P. Daripa, W. Hua, A numerical method for solving an illposed Boussinesq equation arising in water waves and nonlinear lattices, Appl. Math. Comput. 101 (1999) 159–207. [5] R.K. Dash, P. Daripa, Analytical and numerical studies of a singularly perturbed Boussinesq equation, Appl. Math. Comput. 126 (1) (2002) 1–30, Feb.

M.A. Hajji, K. Al-Khaled / Applied Mathematics and Computation 191 (2007) 320–333

333

[6] S.M. El-Sayad, D. Kaya, The decomposition method for solving (2 + 1)-dimensional Boussinesq equation and (3 + 1)-dimensional KP equation, Appl. Math. Comput. 157 (2) (2004) 523–534. [7] Z. Feng, Traveling solitary wave solutions to the generalized Boussinesq equation, Wave Motion 37 (2003) 17–23. [8] R. Hirota, Exact N-soliton solution of the wave equation of long waves in shallow-water and in nonlinear lattices, J. Math. Phys. 14 (1973) 810–814. [9] D. Kaya, S.M. El-Sayed, An application of the decomposition method for the generalized KdV and RLW equations, Chaos Solitons Fract. 17 (5) (2003) 869–877. [10] V.S. Manoranjan, A.S. Mitchell, J.L. Morris, Numerical solution of the good Boussinesq equation, SIAM J. Sci. Stat. Comput. 5 (1984) 946–957. [11] A.M. Wazwaz, Construction of soliton solutions and periodic solutions of the Boussinesq equation by the modified decomposition method, Chaos Solitons Fract. 12 (2001) 1549–1556. [12] G.B. Whitham, Comments on periodic waves and solitons, IMA J. Appl. Math. 32 (1–3) (1984) 353–366. [13] M.R. Scott, H.A. Watts, Computational solution of linear two-point boundary value problems via orthonormalization, SIAM J. Numer. Anal. 14 (1977) 40–70. [14] M.R. Scott, H.A. Watts, SUPORT-a computer code for two-point boundary-value problems via orthonormalization, SAND75-0198, Sandia Laboratories, Albuquerque, NM, 1975. [15] Layne T. Watson, Melvin R. Scott, Solving spline-collocation approximations to nonlinear two-point boundary-value problems by a homotopy method, Appl. Math. Comput. 24 (1987) 333–357. [16] Layne Watson, Engineering applications of the Chow–Yorke algorithm, Appl. Math. Comput. 9 (1981) 111–133. [17] M.R. Scott, H.A. Watts, A systematized collection of codes for solving two-point boundary-value problems, in: L. Lapidus, W.E. Schiesser (Eds.), Numerical Methods for Differential Systems, Academic, New York, 1976.