A quintic B-spline finite elements scheme for the KdVB equation

Comput. Methods Appl. Mech. Engrg. 188 (2000) 121±134

www.elsevier.com/locate/cma

A quintic B-spline ®nite elements scheme for the KdVB equation S.I. Zaki Mathematics Department, Science College, Suez Canal University, Ismailia, Egypt Received 4 December 1998; received in revised form 7 March 1999

Abstract An algorithm for solving the Korteweg de±vries Burgers' (KdVB) equation, based on the collocation method with quintic B-spline ®nite elements is set up to simulate the solutions of the KdV, Burgers' and the KdVB equations. The migration of solitary waves and temporal evolution of a Maxwellian initial pulse are studied for the KdV equation. Burgers' equation is solved for di€erent values of Reynolds number. The time evaluation of the solutions of the KdVB equation with di€erent values for the di€usion and dispersion coecients is studied. Invariants and error norms are studied wherever possible to determine the conservation properties of the algorithm. A linear stability analysis shows the scheme to be unconditionally stable. Ó 2000 Elsevier Science S.A. All rights reserved. Keywords: Korteweg de±Vries Burgers' equation; Finite element method; Quintic B-splines; Collocation

1. Introduction The Korteweg de±Vries Burgers' equation (KdVB), derived by Su and Gardner [1], is a model equation for a wide class of nonlinear systems in the weak nonlinearity and long wavelength approximations because it contains both damping and dispersion. This equation has been obtained when including electron inertia e€ects in the description of weak nonlinear plasma waves. The steady-state solution of the KdVB equation has been shown to model weak plasma shocks propagating perpendicular to a magnetic ®eld [2]. When di€usion dominates dispersion, the steady-state solutions of the KdVB equation are monotonic shocks [3], and when dispersion dominates, the shocks are oscillatory. The KdVB equation has also been used in the study of wave propagation through a liquid-®lled elastic tube [4] and for a description of shallow water waves on a viscous ¯uid [5]. Little numerical work has been done to solve the KdVB equation. Brief details of a numerical solution of the KdVB equation using the accurate space derivative method have been given by Canosa and Gazdag [6], who discussed the evolution of nonanalytic initial data into a monotonic shock. Also, Ali et al. [7] have produced a B-spline FE scheme based on Galerkin's method with quadratic B-spline interpolation functions over ®nite elements. In this paper, a B-spline ®nite element algorithm based on the collocation method with trial functions taken as quintic B-spline functions over the elements will be constructed. The present algorithm will be used ®rst to model both the KdV and Burgers' equations and then it will be solved for problems presenting the KdVB equation. The results obtained will be compared where possible with earlier numerical work.

0045-7825/00/$ - see front matter Ó 2000 Elsevier Science S.A. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 9 9 ) 0 0 1 4 2 - 5

122

S.I. Zaki / Comput. Methods Appl. Mech. Engrg. 188 (2000) 121±134

2. The governing equation The KdVB equation derived by Su and Gardner [1] has the form Ut ‡ eUUx ÿ tUxx ‡ lUxxx ˆ 0;

…1†

where e, t and l are positive parameters and the subscripts t and x denote differentiation. When di€usion dominates dispersion (m  l), the numerical solutions of Eq. (1) tend to behave like Burgers' equation [8,9] solutions. Whereas when dispersion dominates (l  m), KdV behavior is observed. On applying a Galerkin approach, we obtain a weak form of Eq. (1) as Z b Wm fUt ‡ eUUx ÿ mUxx ‡ lUxxx g dx ˆ 0; …2† a

where Wm is a weight function. If now the weight functions are identi®ed with the Diarc delta functions Wm ˆ d…xm ÿ x†; then, the resulting set of equations takes the form Ut ‡ eUUx ÿ tUxx ‡ lUxxx xˆxm ˆ 0; m ˆ 0; 1; 2; . . . ; N :

…3†

…4†

These equations may also be identi®ed as point collocation conditions. Boundary conditions will be chosen from the physical boundary conditions associated with the KdVB equation U …a; t† ˆ U …b; t† ˆ 0;

…5a†

and the collocation boundary conditions required for ensuring a unique quintic B-Spline solution Ux …a; t† ˆ Ux …b; t† ˆ 0;

…5b†

Uxx …a; t† ˆ Uxx …b; t† ˆ 0:

…5c†

The presence of the second and third spatial derivatives in Eq. (4) requires that the interpolation functions and their ®rst and second derivatives must be continuous throughout the region. When using the method of collocation, with collocation points identi®ed with the element nodes, the quintic B-spline interpolation functions can be used with partial di€erential equations containing derivatives up to order four.

3. Quintic B-spline solution A numerical solution to the KdVB Eq. (1) is thought over the ®nite region ‰a; bŠ with boundary conditions chosen from Eqs. (5a)±(5c). A partition of the region ‰a; bŠ by uniformly spaced points xj is set up such that a < x1 < x2 <


< b and /m are those quintic B-splines with knots at the points xm . The set of functions f/ÿ2 ; /ÿ1 ; . . . ; /N ‡2 g forms a basis for functions de®ned over the ®nite region ‰a; bŠ. A global approximation UN …x; t† to the solution U …x; t† is given by UN …x; t† ˆ

N ‡2 X

dm …t†/m …x†;

…6†

mˆÿ2

where the dm are time-dependent parameters to be determined from the boundary conditions (5) and the set of Eq. (4). The B-spline /m …x† and its three principle derivatives varnish outside the region ‰xmÿ3 ; xm‡3 Š. The intervals ‰xm ; xm‡1 Š are identi®ed with the ®nite elements, which are each covered by six B-splines. Over the element ‰xm ; xm‡1 Š the variation of the function U …x; t† is formed from

S.I. Zaki / Comput. Methods Appl. Mech. Engrg. 188 (2000) 121±134

U …x; t† ˆ

m‡3 X

dj …t†/j …x†:

123

…7†

jˆmÿ2

In terms of a local coordinate system f given by h f ˆ x ÿ xm ; where h ˆ …xm‡1 ÿ xm † and 0 6 f 6 1, expressions for the element splines are [10] /mÿ2 ˆ 1 ÿ 5f ‡ 10f2 ÿ 10f3 ‡ 5f4 ÿ f5 ; /mÿ1 ˆ 26 ÿ 50f ‡ 20f2 ‡ 20f3 ÿ 20f4 ‡ 5f5 ; /m ˆ 66 ÿ 60f2 ‡ 30f4 ÿ 10f5 ; /m‡1 ˆ 26 ‡ 50f ‡ 20f2 ÿ 20f3 ÿ 20f4 ‡ 10f5 ;

…8†

/m‡2 ˆ 1 ‡ 5f ‡ 10f2 ‡ 10f3 ‡ 5f4 ÿ 5f5 ; /m‡3 ˆ f5 : The nodal values Um ; Um0 ; Um00 and Um000 at the knot xm are given in terms of dm by Um ˆ dmÿ2 ‡ 26dmÿ1 ‡ 66dm ‡ 26dm‡1 ‡ dm‡2 ;

…9†

hUm0 ˆ 5…dm‡2 ‡ 10dm‡1 ÿ 10dmÿ1 ÿ dmÿ2 †;

…10†

h2 Um00 ˆ 20…dm‡2 ‡ 2dm‡1 ÿ 6dm‡2 ‡ 2dmÿ1 ‡ dmÿ2 †;

…11†

h3 Um000 ˆ 60…dm‡2 ÿ 2dm‡1 ‡ 2dmÿ1 ÿ dmÿ2 †;

…12†

where the dashes denote di€erentiation with respect to x. To obtain a recurrence relationship for the numerical solution of Eq. (4), time center on …n ‡ 1=2†Dt, where Dt is the time step, and use a Crank±Nicholson approach with n‡1=2

…Ut †m

ˆ

1 …U n‡1 ÿ Umn †; Dt m

…13†

1 Umn‡1=2 ˆ …Umn‡1 ‡ Umn †; 2 …UUx †mn‡1=2 ˆ

…14†

 1 …UUx †mn‡1 ‡ …UUx †nm ; 2

…15†

where the subscript m is the node label while the superscripts n and n + 1 are time step labels. Eq. (15) is a quasi-linearisation obtained through an arithmetic mean in time. These results are all second order accurate approximations to the values at time …n ‡ 1=2†Dt. Use Eqs. (9)±(15) to evaluate Um and its space derivatives to dnm and we have for each knot an equation relating parameters at adjacent time levels, dn‡1 m n‡1 n‡1 …1 ÿ L ‡ Q ÿ P †dmÿ2 ‡ …26 ÿ 10L ‡ 2Q ‡ 2P †dmÿ1 ‡ …66 ÿ 6Q†dn‡1 m n‡1 ‡ …26 ‡ 10L ‡ 2Q ÿ 2P †dn‡1 m‡1 ‡ …1 ‡ L ‡ Q ‡ P †dm‡2

ˆ …1 ‡ L ÿ Q ‡ p†dnmÿ2 ‡ …26 ‡ 10L ÿ 2Q ÿ 2P †dnmÿ2 ‡ …66 ‡ 6Q†dnm ‡ …26 ÿ 10L ÿ 2Q ‡ 2P †dnm‡1 ‡ …1 ÿ L ÿ Q ÿ P †dnm‡2 ;

…16†

where L ˆ 5eDtEmn =2h;

Q ˆ ÿ10mDt=h2 ;

P ˆ 30lDt=h3

and Emn ˆ dnmÿ2 ‡ 26dnmÿ1 ‡ 66dnm ‡ 26dnm‡1 ‡ dnm‡2 :

…17†

124

S.I. Zaki / Comput. Methods Appl. Mech. Engrg. 188 (2000) 121±134

The system Eq. (16) consists of N ‡ 1 quasi-linear equations in N ‡ 5 unknowns (dÿ2 ; dÿ1 ; d0 ; d1 ; . . . ; dN ; dN ‡1 ; dN ‡2 †T . To obtain a unique solution to this system four additional constrains on the derivatives at the end points are obtained from the collocation boundary conditions (5b) and (5c), so that the set Eq. (16) becomes a matrix equation of the form Adn‡1 ˆ Bdn ;

…18†

where A and B are defective pentadiagonal …N ‡ 1†  …N ‡ 1† matrices. The time evolution of the approximate solution UN …x; t† is determined from that of the vector dn which is found by repeatedly solving the matrix Eq. (18) once the initial vector d0 has been computed from the initial conditions. Since the matrices A and B are defective pentadiagonally, a direct numerical scheme for the solution of the system Eq. (18) exists. The numerical scheme Eq. (16) is second order accurate in the time step Dt and fourth order accurate in the space step h [10]. 4. The initial state Initial vector d0 is determined by requiring the initial numerical approximation UN …x; 0† to satisfy the following conditions: (a) It shall agree with the exact initial condition U …x; 0† at the knots; using Eq. (4) this produces N ‡ 1 equations, and (b) The approximate initial condition UN …x; t† shall obey the collocation boundary conditions U0 ˆ U00 ˆ 0 at both ends of the range which produces four further equations. The initial vector d0 is then determined as the solution of the matrix equation, Md0 ˆ b; where

2

20 6 ÿ5 6 6 1 6 6 6 6 6 6 M ˆ6 6 6 6 6 6 6 6 4

…19†

40 ÿ50 26

ÿ120 0 66

40 20 50 5 26 1

1

26

66 26

3

.. 1 . .. .. .. . . . .. .. .. . . . .. .. .. . . .

d0 ˆ …dÿ2 ; dÿ1 ; d0 ; d1 ; . . . ; dN ; dN ‡1 ; dN ‡2 †

1 ÿ5 20

26 ÿ50 40

66 0 ÿ120

26 50 40

7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 1 7 7 5 5 20

…20†

T

and T

b ˆ ‰U 00 …x0 †; U 0 …x0 †; U …x0 †; . . . ; U …xN †; U 0 …xN †; U 00 …xN †Š :

5. Stability analysis The stability analysis is based on the Neumann theory in which the growth factor of the error in a typical n mode of amplitude d^ is n dnj ˆ d^ eijkh ;

…21†

S.I. Zaki / Comput. Methods Appl. Mech. Engrg. 188 (2000) 121±134

125

where K is the mode number and h the element size is determined from a linearisation of the numerical scheme. The di€erential equation (1) is linearised by assuming that the quantity U in the nonlinear term UUx is locally constant. This is equivalent to assuming that in Eq. (17) all the dnm are equal to a local constant d, so that Emn ˆ 120d. Substituting the Fourier mode Eq. (21) into Eq. (16) shows that the growth factor for mod K is gˆ

a ÿ ib ; a ‡ ib

…22†

where a ˆ …1 ‡ Q† cos2kh ‡ …26 ‡ 2Q† cos kh ‡ …33 ‡ 3Q† and b ˆ …L ‡ P † sin 2kh ‡ …10L ÿ 2P † sin kh: Thus jgj is one and the linearised scheme is unconditionally stable. 6. KdV simulations A numerical algorithm for the solution of the KdV equation should possess the following properties [11,12]: (i) The migration of solutions should be adequately described. (ii) the numerical scheme should exhibit the same conservation laws as the di€erential equation. The numerical algorithm developed in Section 3 will be validated by studying test problems concerned with the migration of solitons and the break up of an initial perturbation in the form of Maxwellian initial condition into a train of solitons. We use the L2 and L00 error norms to measure the difference between the numerical and analytic solutions and hence to show how well the scheme predicts the position and amplitude of the solution as the simulation proceeds. The L2 and L00 error norms of the solution are de®ned by Zabusky [12] "

L2 ˆ k U

exact

2 N X ÿ U k2 ˆ h Ujexact ÿ Ujn

#1=2

n

;

…23†

jˆ1

L00 ˆ kU exact ÿ U n k00 ˆ max U exact ÿ Ujn :

…24†

j

Also, the conservation properties of the numerical scheme will be examined by calculating the lowest three invariants corresponding to conservation of mass, momentum and energy, viz., Z I1 ˆ

a

b

Z U dx;

I2 ˆ

a

b

U 2 dx;

Z I3 ˆ

a

b

  3l 2 U 3 ÿ …U 0 † dx: e

…25†

We will now compute solutions to Eq. (1), t ˆ 0, for the following problems. (a) The initial condition: U …x; 0† ˆ 3C sech2 …Ax ‡ D†;

…26†

where A; D and C are constants, together with the boundary conditions U …0; t† ˆ U …2; t† ˆ 0 for all times. In this case the KdV equation has an analytic solution of the form [13] U …x; t† ˆ 3Csech2 …Ax ÿ Bt ‡ D†;

…27†

126

S.I. Zaki / Comput. Methods Appl. Mech. Engrg. 188 (2000) 121±134

provided 1 1=2 A ˆ …ec=l† 2

and

1 1=2 B ˆ ec…ec=l† ; 2

…28†

so that (26) becomes a possible initial condition if A ˆ 12 …e=l†1=2 and in fact represents a single soliton moving to the right with velocity eC. To allow for comparison with an earlier work, we take e ˆ 1; C ˆ 0:3; D ˆ ÿ6, l ˆ 4:84  10ÿ4 , Dt ˆ 0:005; h ˆ 0:001 and t ˆ 0 in Eq. (1). For the present case, the solution will move to the right with speed eC. When the numerical solution and the exact solution are plotted on the same diagram the curves are indistinguishable. The agreement is excellent. To make this observation quantitative, we have recorded the error norms L2 and L00 as well as the ®rst three invariants I1 , I2 and I3 , Table 1, for times up to t ˆ 3.0. As it can be seen in Table 1, all three invariants are satisfactorily constant changing by less than 0.05% during the simulation. The L2 norm is less than 9  10ÿ5 while the L00 norm is less than 4  10ÿ4 and so are satisfactorily small. Our results compare favorably with the work of Gardner et al. [13]. (b) The second test problem we shall study will be the birth of solitons from the initial condition. U …x; 0† ˆ exp …ÿx2 †;

…29†

with boundary conditions U …ÿ15; t† ˆ U …15; t† ˆ 0;

t > 0:

…30†

With the Maxwellian initial condition (29), the behavior of the solution vector depends on some critical value for l, say lc [14]. For l  lc the initial Maxwellian breaks up into a number of solitons according to the theoretical formula N ˆ ‰…1=13l1=2 ] [17]. For l  lc the Maxwellian does not break up into solitons but exhibits rapidly oscillating wave packets. When l  lc a mixed type of solutions is found which consists of a leading soliton and an oscillating tail. For Eq. (29), the critical value of l has been given as lc ˆ 0.0625 [15]. We have considered values for l from lc downwards. For comparisons with earlier work [13,18] we consider for l the set of values 0.04, 0.01, 0.001, 0.0005 and run the simulations up to time t ˆ 12 with e ˆ 1.0, Dt ˆ 0:03 and h ˆ 0:02. For l ˆ 0.04 we observed one solitary wave plus an oscillating tail, Fig. 1(A). The actual velocity Vn has been measured and also computed from the measured amplitude a using the formula Va ˆ ae=3. The calculated value for Va is Va ˆ 1.993 1/3 ˆ 0.3998 and that for Vn is Vn ˆ 7.94/12 ˆ 0.4116, so we have an agreement Vn ˆ Va and con®rmation that the solitary waves are indeed solitons. When l ˆ 0:01 three solitons are formed as shown in Fig. 1(B). We have measured the velocity of the largest solitary wave as Vn ˆ 0.528 and calculated the expected velocity from the observed amplitude 1.5498 as Va ˆ 1.5498 1/3 ˆ 0.5166. These values are again in close agreement. For l ˆ 0.001 nine solitons form, Fig. 1(C), while for l ˆ 0:0005 twelve solitons are formed, Fig. 1(D), in agreement with the theoretical formula. The recorded values of the invariants I1 , I2 and I3 for the last two cases of l ˆ 0.001, 0.0005 are recorded in Table 2. For the case l ˆ 0:001, I1 changes by less than 0.05% while I2 and I3 change by less than 1:6  10ÿ3 % and 0.35%, respectively, while for l ˆ 0:0005 these three invariants change by less than 5:6  10ÿ4 %, 0.013% and 0.68%, respectively. The conservation properties are all very good. The agreement is excellent with an earlier work [13]. Table 1 Invariants and error norms for single soliton Time

I1

I2

I3

L2  103

L00  103

0.0 1.0 2.0 3.0

0.14460 0.14461 0.14461 0.14460

0.08676 0.08677 0.08677 0.08676

0.04685 0.04686 0.04686 0.04685

0.0001 0.04963 0.03967 0.02984

0.0005 0.18851 0.11571 0.07525

S.I. Zaki / Comput. Methods Appl. Mech. Engrg. 188 (2000) 121±134

127

Fig. 1. KdV type solutions, Maxwellian initial condition, state at time t ˆ 12.5, different values for l. (A) One soliton plus an oscillating tail, l ˆ 0:04. (B) Three solitons, l ˆ 0:01. (C) Nine solitons, l ˆ 0:001. (D) Twelve solitons, l ˆ 0:0005.

Table 2 Invariants for Maxwellian initial condition Time

0.0 3.0 6.0 9.0 12.0

l ˆ 510ÿ4

l ˆ 0.001 I1

I2

I3

I1

I2

I3

1.77245 1.77245 1.77245 1.77245 1.77335

1.25331 1.25331 1.25331 1.25331 1.25333

1.01957 1.02286 1.02320 1.02321 1.02316

1.77245 1.77245 1.77242 1.77239 1.77244

1.25331 1.25331 1.25326 1.25320 1.25315

1.02145 1.02813 1.02871 1.02860 1.02848

128

S.I. Zaki / Comput. Methods Appl. Mech. Engrg. 188 (2000) 121±134

7. Burgers' type solutions To solve the KdVB equation (1) as a Burgers' type equation (l ˆ 0), it is very useful to consider as an initial condition the function [8,16] x=t ; …31† U …x; t† ˆ 1 ‡ …t=t0 † exp …x2 =4tt† where t0 ˆ exp(1=8t), evaluated at t ˆ 1, and solve the system of equations for different values of t with boundary conditions U …a; t† ˆ U …b; t† ˆ 0;

8t P 1:

…32†

The initial condition (31) is very useful since the resulting analytic solution can be expressed in closed form so that the L2 and L00 error norms are easily calculated for any value of t. The length of the region of solution will be detected by the spread of the solution. To permit for comparisons with earlier work [16], we consider for t the set of values 0.5, 0.05, 0.005, 5 10ÿ4 . Fig. 2(A) represents the development of the initial

Fig. 2. Burgers' type solutions taken every 50 time steps with Dt ˆ 0.02, top curve at time t ˆ 1 and the bottom at time t ˆ 4. (A) t ˆ 0:5. (B) t ˆ 0:05. (C) t ˆ 0:005. (D) t ˆ 0:0005.

S.I. Zaki / Comput. Methods Appl. Mech. Engrg. 188 (2000) 121±134

129

condition (31) with time for e ˆ 1, l ˆ 0, t ˆ 0.5, h ˆ 0:01 and Dt ˆ 0.02 for 0 6 x 6 1. The simulation is run up to time t ˆ 5.0 and recorded every 50 time steps. The top curve is at t ˆ 1 while the bottom curve is at t ˆ 4.0. To assess convergence, the error norms are recorded in Table 3. The L2 norm is less than 210ÿ6 while the L00 norm is less than 4.10ÿ6 , very satisfactory results. When t is reduced to 0.05 the propagation front becomes steeper. Fig. 2(B) represents the development of the solution with time up to t ˆ 5.0 recorded every 50 time steps for 0 6 x 6 3, h ˆ 0.01 and Dt ˆ 0:02. The error norms are recorded also in Table 4. L2 is less than 1:0  10ÿ6 and L00 is less than 2:0  10ÿ6 . Fig. 4(C) represents the development of the solution for t ˆ 0.005 with 0 6 x 6 1:2; h ˆ 0:001 and Dt ˆ 0:02. The error norms for this case are recorded in Table 3, L2 is less than 3 10ÿ6 while L1 is less than 2 10ÿ5 . Also for a considerable small viscosity t ˆ 5 10ÿ4 , the maximum of the L2 norm is less than 5 10ÿ6 while that of the L00 norm is less than 1 10ÿ5 . Fig. 4(D) represents the solution vectors for this case. All the results reported for different values of t are in very good agreement with the work of Ali et al. [16]. As a second example, consider the particular solution of Burgers' equation [16] U …x; t† ˆ ‰a ‡ b ‡ …b ÿ a† exp …c†Š=‰1 ‡ exp…c†Š;

…33†

where c ˆ a…x ÿ bt ÿ k). So that comparison can be made with earlier work, the constants in Eq. (33) are chosen to have the values a ˆ 0:4; b ˆ 0:6; k ˆ 0:125 and t ˆ 0:01. The initial condition is obtained from Eq. (33) when t ˆ 0. The boundary conditions are U …0; t† ˆ 1 and U …1; t† ˆ 0:2. Fig. 3, represents snap shots of the solution vector taken every 250 time steps with h ˆ 0:0125, Dt ˆ 0:001 and 0 6 x 6 1. The ®rst curve to the left is at time t ˆ 0, while the last curve to the right is at time t ˆ 1. From Fig. 3, it can be seen that the solution represents a traveling wave, initially situated at x ˆ k, moving to the right with speed b.

Table 3 Error norms for di€erent values of t Time

1 2 3 4 5

m ˆ 0:5

m ˆ 0:05

m ˆ 10ÿ4

m ˆ 0:005

L2  105

L1  105

L2  105

L1  105

L2  105

L1  105

L2  105

L1  105

0.0 0.165 0.138 0.116 0.110

0.0 0.274 0.278 0.122 0.389

0.0 0.130 0.114 0.101 0.092

0.0 0.281 0.180 0.130 0.101

0.0 0.012 0.011 0.013 0.212

0.0 0.058 0.060 0.064 1.423

0.0 0.103 0.037 0.197 0.465

0.0 1.163 0.386 4.401 9.952

Table 4 Invariants for KdV type simulation Time

I1

I2

I3

0 100 200 300 400 500 600 700 800

50.00021 50.00034 50.00058 50.00612 50.00237 49.99435 49.97857 49.96607 49.96331

45.00055 45.00003 44.99962 44.99999 44.99921 44.99850 44.99820 44.99815 44.99803

42.30074 42.30028 42.30098 42.30227 42.30135 42.30030 42.29995 42.29979 42.29974

130

S.I. Zaki / Comput. Methods Appl. Mech. Engrg. 188 (2000) 121±134

Fig. 3. Burgers' type solutions taken every 250 time steps with Dt ˆ 0:02, the ®rst curve to the left is at time t ˆ 0 while the last curve to the right is at time t ˆ 0, m ˆ 0.01.

The numerical results for the present method, L1 ˆ 0:0026, and the results quoted from [16], for alternative approaches including collocation with cubic B-spline, L1 ˆ 0:005, for a standard Galerkin approach, L1 ˆ 0:096, a product approximation Galerkin method, L1 ˆ 0:082, and a compact ®nite di€erence method, L1 ˆ 0:151. We see that the present method produces a stable solution closer to the analytical one, over the whole range, than any of the other mentioned methods. Also the L2 norm for the present case is of magnitude 0.684 10ÿ3 . 8. The KDVB type solutions We now examine more carefully the behavior of the KdVB equation (1) and study the e€ect of using di€erent values of l and t on the solution vector. To allow for such a study, we use as an initial condition [7]   j xj ÿ x0 ; …34† U …x; 0† ˆ 0:5 1 ÿ tanh d and boundary conditions U …ÿ50; t† ˆ U …150; t† ˆ 0;

…35†

where ÿ50 6 x 6 150, d ˆ 5 and x0 ˆ 25 will be considered in all the simulations. Fig. 4(A), represents the solution vector after a very long run, time t ˆ 800 with h ˆ 0.05, Dt ˆ 0:4; e ˆ 0:2; l ˆ 0:1; m ˆ 0. In this case Eq. (1) is a KdV type equation and a train of 10 solitons is formed. The invariants I1 , I2 and I3 for this case are recorded in Table 4. As it can be seen from Table 4 the invariants change by less than 0.074%, 0.006% and 0.0024%, respectively, of their original values during this very long run and so can be considered constant. Fig. 4(B) represents the solution vector at time t ˆ 800 with the same set of data except that t has been increased to the new value t ˆ 0:0001, very small viscosity. In fact this graph is indistinguishable from that of Fig. 4(A). Also a train of 10 solitons is formed. The leading soliton in the train has a measured amplitude

S.I. Zaki / Comput. Methods Appl. Mech. Engrg. 188 (2000) 121±134

131

Fig. 4. KdVB type solutions, state at time t ˆ 800, with Dt ˆ 0:4, h ˆ 0.05, e ˆ 0:2; l ˆ 0:1, different values for t. (A) t ˆ 0. (B) t ˆ 0:0001. (C) t ˆ 0:005. (D) t ˆ 0:01. (E) t ˆ 0:03. (F) t ˆ 0:05. (G) t ˆ 0:2. (H) t ˆ 0:4.

a ˆ 1.9288 and a measured velocity V ˆ 0.1328 giving V/a ˆ 0.06885 in agreement with the theoretical value of the ratio for the KdV solution V/a ˆ e/3 ˆ 0.06667 [13]. To study the e€ect of increasing the viscosity and hence the dispersion term on the solution vector, we ®x all the data except for t that takes the values 0.005, 0.01, 0.03, 0.05, 0.2, 0.4. Fig. 4(C)±(H) represent the solution pro®les for these cases at time t ˆ 800, respectively. As it can be seen from these graphs, the more we increase t the more the solution vector for the KdVB equation (1) tends to behave like a solution of

132

S.I. Zaki / Comput. Methods Appl. Mech. Engrg. 188 (2000) 121±134

Burgers' equation (l ˆ 0). This can be seen clearly from Fig. 4(G) and (H) where the solution vectors end up behaving like traveling waves for which the amplitudes are damped. Now, reconsider the initial condition (34) with the boundary conditions U …0; t† ˆ 1

…36†

U …150; t† ˆ 0:

…37†

and

Fig. 5. KdVB type solutions, state at time t ˆ 800, with Dt ˆ 0:4, h ˆ 0.05, e ˆ 0:2; l ˆ 0:1, different values for t. (A) t ˆ 5. (B) t ˆ 0:5. (C) t ˆ 0:4. (D) t ˆ 0:3. (E) t ˆ 0:2. (F) t ˆ 0:1.

S.I. Zaki / Comput. Methods Appl. Mech. Engrg. 188 (2000) 121±134

133

For such boundary conditions, Canosa and Gazdag [6] have shown that the steady-state solution of the KdVB equation exhibits di€erent behaviors depending on the relative values of t and l. For t2 P 4l the solution is a shock wave decreasing monotonically from upstream to downstream, while for t2 < 4l the solution is a shock wave, which becomes oscillatory upstream and monotonic downstream. These observations will be considered in the following simulations. Fig. 5 represents the solution pro®les at time t ˆ 800 with x0 ˆ 75, d ˆ 5, 0 6 x 6 150, h ˆ 0.05, Dt ˆ 0:4; e ˆ 0:2; l ˆ 0:1 for different values of t. For t ˆ 5 the solution pro®le is a shock wave decreasing monotonically from upstream to downstream, Fig. 5(A), con®rming the theoretical predications. For t ˆ 0:5 the shock wave becomes oscillatory upstream and monotonic down stream, Fig. 5(B). As t is decreased further to the values 0.4, 0.3, 0.2, and 0.1 the computed pro®les become more oscillatory, see Fig. 5(C)±(F). All results obtained are consistent with the theoretical predications of Canosa and Gazdag [6]. 9. Conclusion The KdVB equation (1) is a transient nonlinear di€usive dispersive equation so that any numerical scheme that simulates this equation must represent faithfully all the features of the KdV equation (t ˆ 0), Burgers' equation (l ˆ 0) and at the same time cope with the balance between e; l and t of Eq. (1). To ful®ll these requirements, we have constructed a one space dimension B-spline ®nite elements scheme based on the collocation method together with trial functions taken as quintic B-spline functions to cope with the second and third spatial derivatives of Eq.(1). Linear stability analysis proved that the presented scheme is unconditionally stable. All the results obtained represented the physics of the problem, as one would expect. When di€usion dominates the numerical solutions tend to behave like Burgers' equation, whereas when dispersion dominates KdV type behavior is obtained. The numerical scheme presented has been shown to cope with very long runs, t ˆ 800, considered with the simulations of the KdVB equation. This gives us con®dence that the scheme can be used for runs of the KdVB equation which are of long duration. We believe that the scheme presented can be useful for other applications where continuity of derivatives is essential. Acknowledgements The author wishes to thank the referees for several useful suggestions, especially that the ®nal problem studied in Section 8 and also the ®nal problem of Section 9 should be included. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]

C.H. Su, C.S. Gardner, Derivation of the Korteweg de±Vries and Burgers' equation, J. Math. Phys. 10 (1969) 536. H. Grad, P.N. Hu, Uni®ed shock pro®le in a plasma, Phys. Fluids 10 (1967) 2596. J.D. Cole, on a quasilinear parabolic equation occurring in aerodynamics, Quart. Appl. Math. 9 (1951) 225. R.S. Johnson, A nonlinear equation incorporating damping and dispersion, J. Fluid Mech. 42 (1970) 49. R.S. Johnson, Shallow water waves in a viscous ¯uid, The undular bore Phys. Fluids 15 (1972) 693. J. Canosa, J. Gazdag, The Korteweg de±Vries Burgers' equation, J. Comp. Phys. 23 (1977) 393. H.A. Ali, L.R.T. Gardner, G.A. Gardner, Numerical study of the KdVB equation using B-spline ®nite elements, J. Math. Phys. Sci. 27-1 (1993) 37. E. Hopf, The partial di€erential equation Ut ‡ Uux ˆ lUxx , Comm. Pure Appl. Math. 9 (1950) 210. J. Burgers, A mathematical model illustrating the theory of turbulence, Adv. in Appl. Math., Acdemic Press, New York, 1948. P.M. Prenter, Splines, Variational Methods, Wiley, New York, 1975. B. Fornberg, G.B. Whitham, A numerical study of certain nonlinear wave phenomena, Phil. Trans Roy. Soc. 289 (1978) 73. N.J. Zabusky, A synergetic approach to problem of non-linear dispersive wave propagation and interaction, in: Proceedings of the Symposium on Non linear PDES, Academic Press, New York, 1967. L.R.T. Gardner, G.A. Gardner, A.H.A. Ali, Simulations of solitons using quadratic spline ®nite elements, Comp. Meth. Appl. Mech. Eng. 92 (1991) 231. A. Je€rey, T. Kakutani, weak nonlinear dispersive waves: a discussion centered around the KdV equation, SIAM Rev. 14 (1972) 522.

134

S.I. Zaki / Comput. Methods Appl. Mech. Engrg. 188 (2000) 121±134

[15] YU.A. Berezin, V.I. Karpman, Nonlinear evolution of disturbances in plasmas and other dispersive media, Soviet Phys. JETP 24 (1967) 1049. [16] A.H.A. Ali, G.A. Gardner, L.R.T. Gardner, A Collocation solution for Burgers' equation using cubic B-spline ®nite elements, Comput. Meth. Appl. Mech. Eng. 100 (1992) 325. [17] V.I. Karpman, An asymptotic solution to the KdV equation, Soviet Phys. JEPT 25A (1967) 70. [18] K. Goda, On stability of some ®nite di€erence schemes for the KdV equation, J. Phys. Soc. Japan 39 (1975) 229.