- Email: [email protected]
A Petrov–Galerkin method for equal width equation
Applied Mathematics and Computation 218 (2011) 2730–2739
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A Petrov–Galerkin method for equal width equation Thoudam Roshan Manipur University, Canchipur 795003, Imphal, Manipur, India
a r t i c l e
i n f o
Keywords: Finite element Petrov–Galerkin Product approximation Quadratic B-spline Solitary wave Soliton
a b s t r a c t The equal width equation is solved numerically by Petrov–Galerkin method using linear hat function and quadratic B-spline function as trial and test functions respectively. Product approximation has been used in this method. A linear stability analysis of the scheme is shown to be conditionally stable. Test problems including the single soliton and the interaction of solitons are used to validate the suggested method, which is found to be accurate and efficient. The Maxwellian initial condition pulse and the development of an undular bore are also studied. Ó 2011 Elsevier Inc. All rights reserved.
1. Introduction The equal width (EW) equation is a model nonlinear partial differential equation for the simulation of one-dimensional wave propagation in nonlinear media with dispersion process. Solitary waves are wave packets or pulses which propagates in nonlinear dispersive media. Due to balance between the nonlinear and dispersive effects these waves retain a stable waveform. A soliton is a very special type of solitary wave which also keeps its wave form after collision with other solitons. Peregrin [1] was the first to derive the regularized long wave (RLW) equation to model the development of an undular bore. Later on, Benjamin et al. [2] proposed the use of RLW equation as preferred alternative to the more classical Kortewegde Vries (KdV) equation to model a larger class of physical phenomena. The EW equation, which is less well-known and proposed by Morrison et al. [3], is an alternative description to the more usual KdV and RLW equations. It has been shown to have solitary wave solution and govern a large number of important physical phenomena such as the nonlinear transverse waves in shallow water, ion-acoustic and magnetohydrodynamic waves in plasma and phonon packets in nonlinear crystals. In the recent years, various degrees of the B-splines are used in getting the numerical solution of some partial differential equations. The EW equation has been studied using Galerkin method with both cubic B-spline finite elements [4,5]. It has also been studied by Galerkin method using linear finite element [6] and a Petrov–Galerkin method using quadratic B-spline finite elements [7]. Zaki [8] used least squares technique with space–time linear finite elements to construct a numerical solution for the EW equation. Collocation method has been successfully implemented for the EW equation using quadratic, cubic B-spline functions [9,10]. In this paper the Petrov–Galerkin method with the used of linear hat and quadratic B-spline as trial and test functions respectively is shown to represent accurately the migration of solitary wave. The interaction of solitary waves and other properties of the EW equation are also studied. 2. Governing equation The equal width (EW) equation for long waves propagating in the positive x-direction, has the form
ut þ uux duxxt ¼ 0;
E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi: 10.1016/j.amc.2011.08.013
ð1Þ
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739
2731
where d is a positive constant, which require the boundary conditions u ? 0 as x ? ±1. In this paper, we shall use periodic boundary conditions for a region a 6 x 6 b. The form of the initial pulse will be chosen so that at large distances from the pulse juj is extremely small and essentially attains free space boundary condition u = 0. In fluid problems, u is related to the vertical displacement of water surface, whereas in the plasma application u is the negative of the electrostatic potential. The EW equation (1) has a solitary wave solution of the form 2
uðx; tÞ ¼ 3c sec h
1 pffiffiffi ðx ct x0 Þ : 2 d
ð2Þ
With zero boundary conditions, solution of the EW equation (1) possess three invariant constants of motion [11] given by
I1 ¼
Z
b
udx;
I2 ¼
a
Z
b
a
2 u þ du2x dx;
I3 ¼
Z
b
u3 dx
ð3Þ
a
corresponding to conservation of mass, momentum and energy. 3. Petrov–Galerkin method (PGM) For our conenience EW equation (1) is rewritten in the form
1 ut þ ðu2 Þx duxxt ¼ 0: 2
ð4Þ
The space interval a 6 x 6 b is discretized by uniform (N + 1) grid points xj = a + jh, where j = 0, 1, 2, . . . , N and the grid spacing is given by h = (b a)/N. Let Uj(t) denote the approximation to the exact solution u(xj, t). Following the method used by SanzSerna and Christie [12] to solve the KdV equation by Petrov–Galerkin method we modified their method to solve Eq. (4). Using Petrov–Galerkin method, we assume the solution of Eq. (4) as
uðx; tÞ ¼
N X
U j ðtÞ/j ðxÞ:
ð5Þ
j¼0
The product approximation technique [13] is used for the nonlinear term in the following manner
u2 ðx; tÞ ¼
N X
U 2j ðtÞ/j ðxÞ;
ð6Þ
j¼0
where /j(x), j = 0, 1, 2, . . . , N are the usual piecewise linear ‘‘hat’’ function given by
8 > < 1 þ ðx jhÞ=h; x 2 ½xj1 ; xj ; /j ðxÞ ¼ 1 ðx jhÞ=h; x 2 ½xj ; xjþ1 ; > : 0; otherwise:
The unknown functions Uj(t), j = 0, 1, 2, . . . , N are determined from the system of ordinary differential equations
1 ðut ; wj Þ þ ððu2 Þx ; wj Þ dðuxxt ; wj Þ ¼ 0; 2
ð7Þ
where wj, j = 0, 1, 2, . . . , N are test functions, which are taken to be quadratic B-splines given by
8 ðx xj1 Þ2 ; > > > > 1 < 2h2 ðxjþ1 xÞ2 ðx xj Þ2 ; wj ðxÞ ¼ 2 h > > ðxjþ2 xÞ2 ; > > : 0;
xj1 6 x 6 xj ; xj 6 x 6 xjþ1 ; xjþ1 6 x 6 xjþ2 ; otherwise
and (,) denotes the usual inner product
ðf ; gÞ ¼
Z
b
f ðxÞgðxÞdx:
a
Integrating by parts and using the fact that w(a) = w(b) = 0, Eq. (7) leads to the formation
1 ðut ; wj Þ þ ððu2 Þx ; wj Þ þ lðuxt ; ðwj Þx Þ ¼ 0: 2
ð8Þ
Performing the integrations on (8) will give the following system of ordinary differential equations (ODEs)
aU_ j1 þ bU_ j þ bU_ jþ1 þ aU_ jþ2 ¼
1 ½ðU j1 Þ2 þ 3ðU j Þ2 3ðU jþ1 Þ2 ðU jþ2 Þ2 ; 9h
ð9Þ
2732
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739
1 1 where a ¼ 18 1 12d 11 þ 12d ; b ¼ 18 and ‘’ means the derivative with respect to time t. h2 h2 Now to solve the ODEs, we assume that U nj to be fully discrete approximation to the exact solution u(xj, tn), where tn = nDt nþ1 n1 and Dt is the time step size. Using the central difference scheme for the time derivative, U_ ¼ U U , Eq. (9) reduces to the 2Dt
system of equations nþ1 nþ1 n1 n1 n1 aU nþ1 þ bU nþ1 þ bU n1 j1 þ bU j jþ1 þ aU jþ2 ¼ aU j1 þ bU j jþ1 þ aU jþ2
þ
2 2 2 2Dt n 2 ; U j1 þ 3 U nj 3 U njþ1 U njþ2 9h
ð10Þ
where j = 1, 2, 3, . . . , N 2 and n = 1, 2, 3, . . . The system (10) required two initial time levels and in our practical calculations, the exact values at time equals zero and time equals Dt are used to obtained the required initial conditions. 4. Stability analysis To investigate the stability of the scheme equation (10) we apply the Von Neumann stability analysis. Assuming u in the ~ , the linearized scheme is nonlinear term uux of the EW equation as locally constant u nþ1 nþ1 n1 n1 n1 aU nþ1 þ bU nþ1 þ bU n1 j1 þ bU j jþ1 þ aU jþ2 ¼ aU j1 þ bU j jþ1 þ aU jþ2 þ
i ~h n 4Dt u U j1 þ 3U nj 3U njþ1 U njþ2 : 9h
ð11Þ
For convenience, each element [xj, xj+1] is transformed into the element [xj1/2, xj+1/2], then our linearized scheme equation (11) would take the form nþ1 nþ1 nþ1 n1 n1 n1 n1 aU nþ1 j3=2 þ bU j1=2 þ bU jþ1=2 þ aU jþ3=2 ¼ aU j3=2 þ bU j1=2 þ bU jþ1=2 þ aU jþ3=2
þ
i ~h n 4Dt u U j3=2 þ 3U nj1=2 3U njþ1=2 U njþ3=2 : 9h
ð12Þ
Substituting the Fourier mode
U nj ¼ nn eijh ;
h¼
pffiffiffiffiffiffiffi hk and i ¼ 1; 2
ð13Þ
where k is the mode number and h the space step size into the system of Eqs. (12), we obtain
nnþ1 ¼ gnn ;
ð14Þ
where the growth factor, g is determined by
~ sin 3h þ sin h 2 Dt u g 2 þ 2i g 1 ¼ 0; 9h a cos 3h þ b cos h
ð15Þ
1 1 and b ¼ 18 . Let the roots of Eq. (15) be g1 and g2 so that g1g2 = 1, and where a ¼ 18 1 12d 11 þ 12d h2 h2
g1 þ g2 ¼
~ i sin 3h þ sin h 4 Dt u : 9h a cos 3h þ b cos h
Dtu~ If we prove that 29h gi
sin 3hþsin h
a cos 3hþb cos h gi
form g1 = e , g2 = e
6 1, then Eq. (15) can be written as g2 + 2sin ggi 1 = 0, from which the roots must be of the
, therefore jg1j = jg2j = 1 and our scheme is conditionally stable. Now we prove the following lemma.
~ and h as defined above, then Lemma. For h; Dt; u
2Dtu ~ sin 3h þ sin h 9h a cos 3h þ b cos h 6 1; for all h; 1 1 and b ¼ 18 . where a ¼ 18 1 12d 11 þ 12d h2 h2 Proof
2Dtu ~ sin 3h þ sin h 4Dt u ~ sin h cos h ¼ : 9h a cos 3h þ b cos h h 1 12d cos2 h þ 2 1 þ 6d h2
Now, suppose x = 1 cos h, then
h2
2733
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739
2 32
~ ~ 2 4 D t u sin h cos h 4 Dt u ð2x x2 Þð1 xÞ2 5 ¼ Sup8h 4 Sup06x62 h i2 : h h 1 12d cos2 h þ 2 1 þ h6d2 ð1 xÞ2 þ 2 1 þ 12d 1 12d 2 2 h2 h
ð2xx2 Þð1xÞ2
Let f ðxÞ ¼ h
i2 and suppose, 1
112d ð1xÞ2 þ2 1þ12d 2 2 h
12d h2
h
2 ð1 xÞ þ 2 1 þ h6d2 ¼ 0, then ð1 xÞ2 ¼ 2
1þ 6d2
h
12d1 h2
and so (1 x)2 6 1 )
h
3 6 0, so we have a contradiction. Hence, 1 12d ð1 xÞ2 þ 2 1 þ 6d2 – 0. Therefore f(x) is continuous on [0, 2]. Now differentiating f with respect to x, to 2 h h ffi ffi pffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 4 2 2 4 2 2 2 find the maximum, then f0 (x) = 0 gives us x1 ¼ 1; x2 ¼ 5h þ12dþ 2 2 5h þ42h dþ72d and x3 ¼ 5h þ12d 2 2 5h þ42h dþ72d Now 5h þ12d
f ðx1 Þ ¼ 00
1
2 1þ 2d2
> 0 and f 00 ðx2 Þ ¼ f 00 ðx3 Þ ¼
ð5h2 þ12dÞ3 216ðh3 þ6hdÞ2
5h þ12d
< 0. fmax ¼ f ðx2 Þ ¼ f ðx3 Þ ¼
h2 . 24ðh2 þ6dÞ
h
Thus,
2 32 2 ~ ~ 2 ~2 4 D t u sin h cos h 4Dt u h 2e2 Dt 2 u 5 ¼ ¼ ¼ R ðsayÞ: Sup8h 4 2 2 h h 2 h þ 2 1 þ 6d 24ðh þ 6dÞ 3ðh þ 6dÞ 1 12d cos 2 2 h h ~ represents single speed and is usually small. So, for The space step h is usually small and the time step Dt is also small. Also u any problem of practical significance R 6 1. Hence the lemma. h
5. Numerical applications In this section, we study test problems concerning the migration of single wave, interactions of two and three solitary waves and the birth of solitons using Maxwellian initial condition. The L2 and L1-error norms are used to measure the accuracy of the present scheme and to compare our results with both exact values, Eq. (2) as well as other results in literature whenever available. 5.1. Single solitary wave First we consider the initial condition
uðx; 0Þ ¼ 3c sec h
2
1 pffiffiffi ðx x0 Þ : 2 d
ð16Þ
Analytical values of the invariants, for a solitary wave of amplitude 3c and width 2p1 ffiffid given by Eq. (2) may be evaluated to give [8]
pffiffiffi I1 ¼ 12c d;
I2 ¼
pffiffiffi pffiffiffi 144c2 d 288c3 d and I3 ¼ : 5 5
ð17Þ
Three sets of parameters have been considered for computational work. The first set is chosen to be c = 0.1, d = 1, h = 0.03, Dt = 0.05 and x0 = 10 with range [0, 30] to coincide with Petrov–Galerkin method of Gardner et al.[5], Galerkin method of Dogan [6] and cubic B-spline collocation method of Dag and Saka [10]. Thus the solitary wave has amplitude 0.3 and simulation are done up to t = 80. Using the values of c and d the value of the three invariants are I1 = 1.2, I2 = 0.288 and I3 = 0.0576. Values of the three invariants as well as L2 and L1-error norms have computed and reported in Table 1. The changes of the invariants from their initial values are less than 0.006%, 0.0003% and 0.0001% respectively. Error deviation are changed in the range of 0.51514 104 6 error 6 0.75975 105. The profile of the solitary wave t = 0 and t = 80 are shown in Fig.1. In the second set we choose c = 0.03, h = 0.05, Dt = 0.05, d = 1 and x0 = 10 with range [0, 30], then the amplitude is 0.09. The simulation are done up to t = 80. Table 2 represent values of the three invariants and error norms to the present scheme in this case. The changes of invariants I1, I2 and I3 from their initial values are less than 0.0014%,
Table 1 Invariants and error norms for single solitary wave: amplitude = 0.3, Dt=0.05, h = 0.03, 0 6 x 6 30 (PGM). t
I1
I2
I3
103 L2
103 L1
0 20 40 60 80
1.9995 1.20004 1.20005 1.20005 1.20004
0.288 0.288 0.288 0.288001 0.288
0.0576 0.0576 0.0576 0.0576 0.0576
0.0 0.0328791 0.0373359 0.0377384 0.0388148
0.0 0.0457024 0.0517833 0.0517373 0.0515138
2734
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739
Fig. 1. Single soliton of amplitude 0.3 at time c = 0.1, t = 0, t = 40 and t = 80.
Table 2 Invariants and error norms for single solitary wave: amplitude = 0.09, Dt = 0.05, h = 0.05, 0 6 x 6 30 (PGM). t
I1
I2
I3
103 L2
103 L1
0 20 40 60 80
0.359983 0.359998 0.360006 0.360011 0.360013
0.02592 0.02592 0.02592 0.02592 0.02592
0.0015552 0.0015552 0.0015552 0.0015552 0.0015552
0.0 0.005465 0.007951 0.009424 0.010249
0.0 0.008969 0.010862 0.012972 0.014128
0.0001% and 0.00001% respectively. The error deviations are changed in the ranges of 1.4128 105 6 error 6 1.4827 106 in this case. In the third set we choose c = 1, h = 0.05, Dt = 0.05, d = 1 and x0 = 15 with range [0, 80], then the amplitude in this case is 3. The simulation are done up to t = 40. Table 3 represent values of the three invariants and error norms to the present scheme in this case.The invariants I1, I2 and I3 are not changed in this case. The numerical results shown in Tables1–3 are compared with the results of the previous publications [5,6,8,10] and it is shown in Table 4. The superiority of the present scheme over the previous ones is clear.
Table 3 Invariants and error norms for single solitary wave: amplitude = 3, Dt = 0.05, h = 0.05, 0 6 x 6 80 (PGM). t
I1
I2
I3
103 L2
103 L1
0 10 20 30 40
11.9994 12.0 12.0 12.0 12.0
28.8 28.8 28.8 28.8 28.8
57.6 57.6 57.6 57.6 57.6
0.0 4.78592 7.21193 9.73328 12.2944
0.0 2.93738 4.28776 5.64959 7.01658
Table 4 Comparisons of results for the single solitary waves for different amplitudes: h = 0.03, Dt = 0.05 for amp. = 0.03 and h = Dt = 0.05 for amp. = 0.09 & 3. Ref.
Amp.
Time
I1
Ours [10] [5] [8] (Dt = 0.03) [6]
0.3
80
1.20004 1.19998 1.1910 1.1964 1.23387
0.288 0.28798 0.2855 0.2858 0.29915
Ours [8] (h = 0.01) [6]
0.09
80
0.36001 0.3597 0.36665
0.02592 0.0259 0.02658
Ours [10] [6]
3
40
12.0 11.9999 11.9977
I2
28.8 28.80644 28.8784
L2 103
L1 103
0.0576 0.05759 0.05582 0.0569 0.06097
0.038815 0.056 3.849 7.444 24.697
0.051514 0.053 2.646 4.373 16.425
0.001555 0.00155 0.00161
0.010249 0.22 2.683
0.014128 0.16 1.836
I3
57.6 57.5983 57.5548
12.2944 17.189 37.365
7.01658 9.3597
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739
2735
5.2. Interaction of two EW solitary waves Here we study the interaction of two well separated EW solitary waves having different amplitudes and traveling in the same direction. The initial condition given by
uðx; 0Þ ¼
2 X
2
3ci sec h
1 pffiffiffi ðx xi Þ ; 2 d
i¼1
ð18Þ
where ci and xi, i = 1, 2 are arbitrary constants. The analytical values of the conservation laws in this case are as follows 2 pffiffiffi X I1 ¼ 12 d ci ; i¼1
pffiffiffi 2 144 d X I2 ¼ c2i ; 5 i¼1
pffiffiffi 2 288 d X I3 ¼ c3i : 5 i¼1
ð19Þ
For computational work we chose d = 1, x1 = 10, x2 = 30, c1 = 0.4 and c2 = 0.2 with interval [0, 100]. The amplitudes are in the ratio 2:1. The analytical values of the invariants are I1 = 7.2, I2 = 5.76 and I3 = 4.1472. Numerical values of the three invariants have been computed and recorded in Table 5 for h = 0.2 and Dt = 0.1. The simulation is done up to t = 160. Changes of the invariants I1, I2 and I3 from their initial values are less than 0.047%, 0.022% and 0.005% respectively. Fig. 2(a)–(d) show the computer plot of the interaction of two solitary waves at different time levels.
Table 5 Invariants for two solitons interactions of amplitudes in the ratio 2:1; h = Dt = 0.1. t
I1
I2
I3
0 40 80 120 160
7.19976 7.20021 7.20021 7.20020 7.20021
5.75988 5.75987 5.75987 5.75998 5.76009
4.14711 4.14712 4.14712 4.14711 4.14713
Fig. 2. Interaction of two solitary waves at (a) t = 0, (b) t = 80, (c) t = 100, and (d) t = 160.
2736
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739
5.3. Interaction of three EW solitary waves For interaction of three solitary waves of different amplitudes and traveling in the same direction we take
uðx; 0Þ ¼
3 X
2
3ci sec h ½0:5ðx xi Þ;
ð20Þ
i¼1
where ci and xi, i = 1, 2, 3 are arbitrary constants. The analytical values of the invariants in this case are 3 pffiffiffi X I1 ¼ 12 d ci ; i¼1
pffiffiffi 3 144 d X I2 ¼ c2i ; 5 i¼1
pffiffiffi 3 288 d X I3 ¼ c3i : 5 i¼1
ð21Þ
For numerical computations, we chose d = 1, x1 = 10, x2 = 25, x3 = 40, c1 = 0.4, c2 = 0.2 and c3 = 0.1 with interval [0, 120]. Thus the amplitudes are in the ratio 4:2:1. The analytical values of the invariants are I1 = 8.4, I2 = 6.048 and I3 = 4.2048. Numerical values of the invariants are calculated using our scheme and reported in Table 6 for h = 0.2 and Dt = 0.1. The simulation is done up to t = 200. Changes of I1, I2 and I3 from their initial values are less than 0.045%, 0.063% and 0.05% respectively. Fig. 3(a)–(d) show the details of interaction of three solitary waves corresponding to the above set of parameters at different time levels.
Table 6 Invariants for three solitons interactions of amplitudes in the ratio 4:2:1; h = 0.1, Dt = 0.5. t
I1
I2
I3
0 40 80 120 160 200
8.39977 8.40020 8.40020 8.40020 8.40020 8.40021
6.04792 6.04789 6.04751 6.04791 6.04793 6.04792
4.20474 4.20470 4.20447 4.20473 4.20473 4.20475
Fig. 3. Interaction of three solitary waves at (a) t = 0, (b) t = 80, (c) t = 140, and (d) t = 200.
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739
2737
5.4. Maxwellian initial condition Evolution of a train of solitary waves of the EW equation has been studied using the Maxwellian initial condition
uðx; 0Þ ¼ Expððx 20Þ2 Þ;
ð22Þ
for various values of d. The computations are carried out and values of the invariants are reported in Table 6 for the cases d = 0.1, 0.05, 0.025 and 0.01 at fixed space–time steps h = 0.05 and Dt = 0.025, respectively. It is to be noted here that when d is reduced more, then more and more solitary waves are formed. A single solitary wave is formed when d = 0.1 as shown in Fig. 4(a). For d = 0.05, the Maxwellian initial pulse breaks up into a train of at least two solitary waves as shown in Fig. 4(b). For d = 0.025, three stable solitary waves are generated as shown in Fig. 4(c). Finally when d = 0.01 the Maxwellian initial condition has split up into four stable solitary waves as shown in Fig. 4(d). All these have been computed at t = 12 and the computed values are reported in Table 7. Changes of the invariantsI2 and I3 from their initial values are less than 0.06% and 0.035% respectively for all values d taken above but I1 does not change in this process. Invariants remain almost constant for larger values of d and changed a little more when we used smaller values of d. 5.5. Development of an undular bore Finally we study the development of an undular bore. For this problem we use the initial condition
h x x i 0 uðx; 0Þ ¼ 0:5u0 1 tanh ; d
ð23Þ
where u(x, 0) denotes the elevation of the water above the equilibrium at time t = 0. The boundary conditions in this case are u ? 0 as x ? 1 and u ? u0 as x ? 1. We choose a range 20 6 x 6 50, with space step h = 0.1 and time step Dt = 0.025 for a simulation with d = 0.16666667, x0 = 0 and d = 2. The development of the undular bore when d = 2 is shown in Fig. 5 for times up to t = 400. Table 8 contains the values of the invariants I1, I2, I3 and position and amplitude of the leading undulation for the present method. The quantities I1, I2, I3 are no longer constant but increase linearly throughout the simulation at the following rates workout by Gardner et al. [5]
Fig. 4. The Maxwellian initial condition at t = 12 with (a) d = 0.1, (b) d = 0.05, (c) d = 0.025, and (d) d = 0.01.
2738
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739 Table 7 Invariants of EW equation using the Maxwellian initial condition at t = 12. d
Time
I1
I2
I3
d = 0.1
0 3 6 9 12
1.77245 1.77245 1.77245 1.77245 1.77245
1.37864 1.37923 1.37880 1.37877 1.37885
1.02333 1.02355 1.02338 1.02336 1.02339
d = 0.05
0 3 6 9 12
1.77245 1.77245 1.77245 1.77245 1.77245
1.31598 1.31648 1.31619 1.31617 1.31612
1.02333 1.02356 1.02340 1.02339 1.02339
d = 0.025
0 3 6 9 12
1.77245 1.77245 1.77245 1.77245 1.77245
1.28464 1.28520 1.28492 1.28418 1.28474
1.02333 1.02357 1.02340 1.02337 1.02337
d = 0.01
0 3 6 9 12
1.77245 1.77245 1.77245 1.77245 1.77245
1.26585 1.26632 1.26599 1.26639 1.26567
1.02333 1.02330 1.02294 1.02295 1.02293
Fig. 5. The development of an undular bore at t = 400 for d = 2.
Table 8 Development of an undular bore with parameters d = 2, U0 = 0.1, h = 0.1, k = 0.025, 20 6 x 6 50. Time
I1
I2
I3
x
u
0 50 100 150 200 250 300 350 400
2.00000 2.24990 2.49940 2.74803 2.99568 3.24238 3.48811 3.73289 3.97672
0.190278 0.223597 0.256823 0.289861 0.322702 0.355344 0.387795 0.420050 0.452121
0.0185000 0.0222486 0.0259850 0.0296895 0.0333601 0.0369971 0.0406011 0.0441725 0.0477118
– 1.04339 3.70855 6.46989 9.34631 12.2749 15.2270 18.1945 21.1660
– 0.113954 0.156979 0.170509 0.174771 0.176732 0.177727 0.178227 0.178410
Z 1 d d 1 I1 ¼ udx ¼ u20 ¼ M1 ; dt dt 1 2 Z 1 d d 2 I2 ¼ ðu2 þ dðux Þ2 Þdx ¼ u30 ¼ M2 ; dt dt 1 3 Z 1 d d 3 4 3 I3 ¼ u dx ¼ u0 ¼ M 3 : dt dt 1 4
ð24Þ
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739
2739
According to the above formulae, the predicted theoretical values are M1 = 5.03 103, M2 = 6.66667 104 and M3 = 7.5 105. The calculated values of the rate of change of I1, I2 and I3 from Table 7 are found to be M1 = 4.9418 103, M2 = 6.5461 104 and M3 = 7.30295 105 respectively. These values are very close to the theoretical values. 6. Conclusion The Petrov–Galerkin method using linear hat function and quadratic B-spline function as test and trial functions respectively was applied to study the solitary waves of the EW equation and it is shown that our scheme is conditionally stable and comparatively accurate than other methods available in the literature. We tested our scheme through single solitary wave in which the analytic solution is known and then extended to study the interaction of solitons where no analytic solution is known. It is also observed that conservation laws were reasonably satisfied for the interaction of two and three solitary waves. The Maxwellian initial condition was also used. Lastly the development of an undular bore was studied. The results indicate that the rates of increase in the invariants I1, I2 and I3 are in good agreement with the theoretical values. Acknowledgement I hereby acknowledge Dr. K.S. Bhamra for extending helping hand in this work and also offer sincere thanks to the reviewer for giving valuable comments and suggestions. References [1] D.H. Peregrine, Calculations of the development of an undular bore, J. Fluid Mech. 25 (2) (1966) 321–330. [2] T.B. Benjamin, J.L. Bona, J.J. Mahony, Model equations for long waves in non-linear dispersive systems, Phil. Trans. Roy. Soc. London A 272 (1972) 47– 78. [3] P.J. Morrison, J.D. Meiss, J.R. Carey, Scattering of RLW solitary waves, Physica 11D (1984) 324–336. [4] L.R.T. Gardner, G.A. Gardner, Solitary waves of equal width equation, J. Comput. Phys. 101 (1992) 218–223. [5] L.R.T. Gardner, G.A. Gardner, F.A. Ayoub, N.K. Amein, Simulation of EW undular bore, Commun. Numer. Methods Eng. 13 (1997) 583–592. [6] A. Dogan, Application of Galerkin’s method to equal width wave equation, J. Appl. Math. Comput. 160 (2005) 65–76. [7] S.I. Zaki, Solitary waves induced by the boundary forced EW equation, Comput. Math. Appl. Mech. Eng. 190 (2001) 4881–4887. [8] S.I. Zaki, A least squares finite element scheme for the EW equation, Comput. Math. Appl. Mech. Eng. 189 (2000) 587–594. [9] D.J. Evan, K.R. Raslan, Solitary waves for the generalized equal width (GEW) equation, Int. J. Comput. Math. 82 (4) (2005) 445–455. [10] I. Dag˘, B. Saka, A cubic B-spline collocation method for EW equation, Math. Comput. Appl. 9 (3) (2004) 381–392. [11] P.J. Olver, Euler operators and conservation laws of the BBM equation, Math. Proc. Camb. Phil. Soc. 85 (1979) 143–159. [12] J.M. Sanz-Serna, I. Christie, Petrov–Galerkin methods for nonlinear dispersive waves, J. Comput. Phys. 39 (1993) 94–103. [13] I. Christie, D.F. Griffths, A.R. Mitchell, J.M. Sanz-Serna, Product approximation for nonlinear problems in finite element methods, IMA J. Numer. Anal. 1 (1981) 253–266.
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A Petrov–Galerkin method for equal width equation Thoudam Roshan Manipur University, Canchipur 795003, Imphal, Manipur, India
a r t i c l e
i n f o
Keywords: Finite element Petrov–Galerkin Product approximation Quadratic B-spline Solitary wave Soliton
a b s t r a c t The equal width equation is solved numerically by Petrov–Galerkin method using linear hat function and quadratic B-spline function as trial and test functions respectively. Product approximation has been used in this method. A linear stability analysis of the scheme is shown to be conditionally stable. Test problems including the single soliton and the interaction of solitons are used to validate the suggested method, which is found to be accurate and efficient. The Maxwellian initial condition pulse and the development of an undular bore are also studied. Ó 2011 Elsevier Inc. All rights reserved.
1. Introduction The equal width (EW) equation is a model nonlinear partial differential equation for the simulation of one-dimensional wave propagation in nonlinear media with dispersion process. Solitary waves are wave packets or pulses which propagates in nonlinear dispersive media. Due to balance between the nonlinear and dispersive effects these waves retain a stable waveform. A soliton is a very special type of solitary wave which also keeps its wave form after collision with other solitons. Peregrin [1] was the first to derive the regularized long wave (RLW) equation to model the development of an undular bore. Later on, Benjamin et al. [2] proposed the use of RLW equation as preferred alternative to the more classical Kortewegde Vries (KdV) equation to model a larger class of physical phenomena. The EW equation, which is less well-known and proposed by Morrison et al. [3], is an alternative description to the more usual KdV and RLW equations. It has been shown to have solitary wave solution and govern a large number of important physical phenomena such as the nonlinear transverse waves in shallow water, ion-acoustic and magnetohydrodynamic waves in plasma and phonon packets in nonlinear crystals. In the recent years, various degrees of the B-splines are used in getting the numerical solution of some partial differential equations. The EW equation has been studied using Galerkin method with both cubic B-spline finite elements [4,5]. It has also been studied by Galerkin method using linear finite element [6] and a Petrov–Galerkin method using quadratic B-spline finite elements [7]. Zaki [8] used least squares technique with space–time linear finite elements to construct a numerical solution for the EW equation. Collocation method has been successfully implemented for the EW equation using quadratic, cubic B-spline functions [9,10]. In this paper the Petrov–Galerkin method with the used of linear hat and quadratic B-spline as trial and test functions respectively is shown to represent accurately the migration of solitary wave. The interaction of solitary waves and other properties of the EW equation are also studied. 2. Governing equation The equal width (EW) equation for long waves propagating in the positive x-direction, has the form
ut þ uux duxxt ¼ 0;
E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi: 10.1016/j.amc.2011.08.013
ð1Þ
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739
2731
where d is a positive constant, which require the boundary conditions u ? 0 as x ? ±1. In this paper, we shall use periodic boundary conditions for a region a 6 x 6 b. The form of the initial pulse will be chosen so that at large distances from the pulse juj is extremely small and essentially attains free space boundary condition u = 0. In fluid problems, u is related to the vertical displacement of water surface, whereas in the plasma application u is the negative of the electrostatic potential. The EW equation (1) has a solitary wave solution of the form 2
uðx; tÞ ¼ 3c sec h
1 pffiffiffi ðx ct x0 Þ : 2 d
ð2Þ
With zero boundary conditions, solution of the EW equation (1) possess three invariant constants of motion [11] given by
I1 ¼
Z
b
udx;
I2 ¼
a
Z
b
a
2 u þ du2x dx;
I3 ¼
Z
b
u3 dx
ð3Þ
a
corresponding to conservation of mass, momentum and energy. 3. Petrov–Galerkin method (PGM) For our conenience EW equation (1) is rewritten in the form
1 ut þ ðu2 Þx duxxt ¼ 0: 2
ð4Þ
The space interval a 6 x 6 b is discretized by uniform (N + 1) grid points xj = a + jh, where j = 0, 1, 2, . . . , N and the grid spacing is given by h = (b a)/N. Let Uj(t) denote the approximation to the exact solution u(xj, t). Following the method used by SanzSerna and Christie [12] to solve the KdV equation by Petrov–Galerkin method we modified their method to solve Eq. (4). Using Petrov–Galerkin method, we assume the solution of Eq. (4) as
uðx; tÞ ¼
N X
U j ðtÞ/j ðxÞ:
ð5Þ
j¼0
The product approximation technique [13] is used for the nonlinear term in the following manner
u2 ðx; tÞ ¼
N X
U 2j ðtÞ/j ðxÞ;
ð6Þ
j¼0
where /j(x), j = 0, 1, 2, . . . , N are the usual piecewise linear ‘‘hat’’ function given by
8 > < 1 þ ðx jhÞ=h; x 2 ½xj1 ; xj ; /j ðxÞ ¼ 1 ðx jhÞ=h; x 2 ½xj ; xjþ1 ; > : 0; otherwise:
The unknown functions Uj(t), j = 0, 1, 2, . . . , N are determined from the system of ordinary differential equations
1 ðut ; wj Þ þ ððu2 Þx ; wj Þ dðuxxt ; wj Þ ¼ 0; 2
ð7Þ
where wj, j = 0, 1, 2, . . . , N are test functions, which are taken to be quadratic B-splines given by
8 ðx xj1 Þ2 ; > > > > 1 < 2h2 ðxjþ1 xÞ2 ðx xj Þ2 ; wj ðxÞ ¼ 2 h > > ðxjþ2 xÞ2 ; > > : 0;
xj1 6 x 6 xj ; xj 6 x 6 xjþ1 ; xjþ1 6 x 6 xjþ2 ; otherwise
and (,) denotes the usual inner product
ðf ; gÞ ¼
Z
b
f ðxÞgðxÞdx:
a
Integrating by parts and using the fact that w(a) = w(b) = 0, Eq. (7) leads to the formation
1 ðut ; wj Þ þ ððu2 Þx ; wj Þ þ lðuxt ; ðwj Þx Þ ¼ 0: 2
ð8Þ
Performing the integrations on (8) will give the following system of ordinary differential equations (ODEs)
aU_ j1 þ bU_ j þ bU_ jþ1 þ aU_ jþ2 ¼
1 ½ðU j1 Þ2 þ 3ðU j Þ2 3ðU jþ1 Þ2 ðU jþ2 Þ2 ; 9h
ð9Þ
2732
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739
1 1 where a ¼ 18 1 12d 11 þ 12d ; b ¼ 18 and ‘’ means the derivative with respect to time t. h2 h2 Now to solve the ODEs, we assume that U nj to be fully discrete approximation to the exact solution u(xj, tn), where tn = nDt nþ1 n1 and Dt is the time step size. Using the central difference scheme for the time derivative, U_ ¼ U U , Eq. (9) reduces to the 2Dt
system of equations nþ1 nþ1 n1 n1 n1 aU nþ1 þ bU nþ1 þ bU n1 j1 þ bU j jþ1 þ aU jþ2 ¼ aU j1 þ bU j jþ1 þ aU jþ2
þ
2 2 2 2Dt n 2 ; U j1 þ 3 U nj 3 U njþ1 U njþ2 9h
ð10Þ
where j = 1, 2, 3, . . . , N 2 and n = 1, 2, 3, . . . The system (10) required two initial time levels and in our practical calculations, the exact values at time equals zero and time equals Dt are used to obtained the required initial conditions. 4. Stability analysis To investigate the stability of the scheme equation (10) we apply the Von Neumann stability analysis. Assuming u in the ~ , the linearized scheme is nonlinear term uux of the EW equation as locally constant u nþ1 nþ1 n1 n1 n1 aU nþ1 þ bU nþ1 þ bU n1 j1 þ bU j jþ1 þ aU jþ2 ¼ aU j1 þ bU j jþ1 þ aU jþ2 þ
i ~h n 4Dt u U j1 þ 3U nj 3U njþ1 U njþ2 : 9h
ð11Þ
For convenience, each element [xj, xj+1] is transformed into the element [xj1/2, xj+1/2], then our linearized scheme equation (11) would take the form nþ1 nþ1 nþ1 n1 n1 n1 n1 aU nþ1 j3=2 þ bU j1=2 þ bU jþ1=2 þ aU jþ3=2 ¼ aU j3=2 þ bU j1=2 þ bU jþ1=2 þ aU jþ3=2
þ
i ~h n 4Dt u U j3=2 þ 3U nj1=2 3U njþ1=2 U njþ3=2 : 9h
ð12Þ
Substituting the Fourier mode
U nj ¼ nn eijh ;
h¼
pffiffiffiffiffiffiffi hk and i ¼ 1; 2
ð13Þ
where k is the mode number and h the space step size into the system of Eqs. (12), we obtain
nnþ1 ¼ gnn ;
ð14Þ
where the growth factor, g is determined by
~ sin 3h þ sin h 2 Dt u g 2 þ 2i g 1 ¼ 0; 9h a cos 3h þ b cos h
ð15Þ
1 1 and b ¼ 18 . Let the roots of Eq. (15) be g1 and g2 so that g1g2 = 1, and where a ¼ 18 1 12d 11 þ 12d h2 h2
g1 þ g2 ¼
~ i sin 3h þ sin h 4 Dt u : 9h a cos 3h þ b cos h
Dtu~ If we prove that 29h gi
sin 3hþsin h
a cos 3hþb cos h gi
form g1 = e , g2 = e
6 1, then Eq. (15) can be written as g2 + 2sin ggi 1 = 0, from which the roots must be of the
, therefore jg1j = jg2j = 1 and our scheme is conditionally stable. Now we prove the following lemma.
~ and h as defined above, then Lemma. For h; Dt; u
2Dtu ~ sin 3h þ sin h 9h a cos 3h þ b cos h 6 1; for all h; 1 1 and b ¼ 18 . where a ¼ 18 1 12d 11 þ 12d h2 h2 Proof
2Dtu ~ sin 3h þ sin h 4Dt u ~ sin h cos h ¼ : 9h a cos 3h þ b cos h h 1 12d cos2 h þ 2 1 þ 6d h2
Now, suppose x = 1 cos h, then
h2
2733
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739
2 32
~ ~ 2 4 D t u sin h cos h 4 Dt u ð2x x2 Þð1 xÞ2 5 ¼ Sup8h 4 Sup06x62 h i2 : h h 1 12d cos2 h þ 2 1 þ h6d2 ð1 xÞ2 þ 2 1 þ 12d 1 12d 2 2 h2 h
ð2xx2 Þð1xÞ2
Let f ðxÞ ¼ h
i2 and suppose, 1
112d ð1xÞ2 þ2 1þ12d 2 2 h
12d h2
h
2 ð1 xÞ þ 2 1 þ h6d2 ¼ 0, then ð1 xÞ2 ¼ 2
1þ 6d2
h
12d1 h2
and so (1 x)2 6 1 )
h
3 6 0, so we have a contradiction. Hence, 1 12d ð1 xÞ2 þ 2 1 þ 6d2 – 0. Therefore f(x) is continuous on [0, 2]. Now differentiating f with respect to x, to 2 h h ffi ffi pffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 4 2 2 4 2 2 2 find the maximum, then f0 (x) = 0 gives us x1 ¼ 1; x2 ¼ 5h þ12dþ 2 2 5h þ42h dþ72d and x3 ¼ 5h þ12d 2 2 5h þ42h dþ72d Now 5h þ12d
f ðx1 Þ ¼ 00
1
2 1þ 2d2
> 0 and f 00 ðx2 Þ ¼ f 00 ðx3 Þ ¼
ð5h2 þ12dÞ3 216ðh3 þ6hdÞ2
5h þ12d
< 0. fmax ¼ f ðx2 Þ ¼ f ðx3 Þ ¼
h2 . 24ðh2 þ6dÞ
h
Thus,
2 32 2 ~ ~ 2 ~2 4 D t u sin h cos h 4Dt u h 2e2 Dt 2 u 5 ¼ ¼ ¼ R ðsayÞ: Sup8h 4 2 2 h h 2 h þ 2 1 þ 6d 24ðh þ 6dÞ 3ðh þ 6dÞ 1 12d cos 2 2 h h ~ represents single speed and is usually small. So, for The space step h is usually small and the time step Dt is also small. Also u any problem of practical significance R 6 1. Hence the lemma. h
5. Numerical applications In this section, we study test problems concerning the migration of single wave, interactions of two and three solitary waves and the birth of solitons using Maxwellian initial condition. The L2 and L1-error norms are used to measure the accuracy of the present scheme and to compare our results with both exact values, Eq. (2) as well as other results in literature whenever available. 5.1. Single solitary wave First we consider the initial condition
uðx; 0Þ ¼ 3c sec h
2
1 pffiffiffi ðx x0 Þ : 2 d
ð16Þ
Analytical values of the invariants, for a solitary wave of amplitude 3c and width 2p1 ffiffid given by Eq. (2) may be evaluated to give [8]
pffiffiffi I1 ¼ 12c d;
I2 ¼
pffiffiffi pffiffiffi 144c2 d 288c3 d and I3 ¼ : 5 5
ð17Þ
Three sets of parameters have been considered for computational work. The first set is chosen to be c = 0.1, d = 1, h = 0.03, Dt = 0.05 and x0 = 10 with range [0, 30] to coincide with Petrov–Galerkin method of Gardner et al.[5], Galerkin method of Dogan [6] and cubic B-spline collocation method of Dag and Saka [10]. Thus the solitary wave has amplitude 0.3 and simulation are done up to t = 80. Using the values of c and d the value of the three invariants are I1 = 1.2, I2 = 0.288 and I3 = 0.0576. Values of the three invariants as well as L2 and L1-error norms have computed and reported in Table 1. The changes of the invariants from their initial values are less than 0.006%, 0.0003% and 0.0001% respectively. Error deviation are changed in the range of 0.51514 104 6 error 6 0.75975 105. The profile of the solitary wave t = 0 and t = 80 are shown in Fig.1. In the second set we choose c = 0.03, h = 0.05, Dt = 0.05, d = 1 and x0 = 10 with range [0, 30], then the amplitude is 0.09. The simulation are done up to t = 80. Table 2 represent values of the three invariants and error norms to the present scheme in this case. The changes of invariants I1, I2 and I3 from their initial values are less than 0.0014%,
Table 1 Invariants and error norms for single solitary wave: amplitude = 0.3, Dt=0.05, h = 0.03, 0 6 x 6 30 (PGM). t
I1
I2
I3
103 L2
103 L1
0 20 40 60 80
1.9995 1.20004 1.20005 1.20005 1.20004
0.288 0.288 0.288 0.288001 0.288
0.0576 0.0576 0.0576 0.0576 0.0576
0.0 0.0328791 0.0373359 0.0377384 0.0388148
0.0 0.0457024 0.0517833 0.0517373 0.0515138
2734
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739
Fig. 1. Single soliton of amplitude 0.3 at time c = 0.1, t = 0, t = 40 and t = 80.
Table 2 Invariants and error norms for single solitary wave: amplitude = 0.09, Dt = 0.05, h = 0.05, 0 6 x 6 30 (PGM). t
I1
I2
I3
103 L2
103 L1
0 20 40 60 80
0.359983 0.359998 0.360006 0.360011 0.360013
0.02592 0.02592 0.02592 0.02592 0.02592
0.0015552 0.0015552 0.0015552 0.0015552 0.0015552
0.0 0.005465 0.007951 0.009424 0.010249
0.0 0.008969 0.010862 0.012972 0.014128
0.0001% and 0.00001% respectively. The error deviations are changed in the ranges of 1.4128 105 6 error 6 1.4827 106 in this case. In the third set we choose c = 1, h = 0.05, Dt = 0.05, d = 1 and x0 = 15 with range [0, 80], then the amplitude in this case is 3. The simulation are done up to t = 40. Table 3 represent values of the three invariants and error norms to the present scheme in this case.The invariants I1, I2 and I3 are not changed in this case. The numerical results shown in Tables1–3 are compared with the results of the previous publications [5,6,8,10] and it is shown in Table 4. The superiority of the present scheme over the previous ones is clear.
Table 3 Invariants and error norms for single solitary wave: amplitude = 3, Dt = 0.05, h = 0.05, 0 6 x 6 80 (PGM). t
I1
I2
I3
103 L2
103 L1
0 10 20 30 40
11.9994 12.0 12.0 12.0 12.0
28.8 28.8 28.8 28.8 28.8
57.6 57.6 57.6 57.6 57.6
0.0 4.78592 7.21193 9.73328 12.2944
0.0 2.93738 4.28776 5.64959 7.01658
Table 4 Comparisons of results for the single solitary waves for different amplitudes: h = 0.03, Dt = 0.05 for amp. = 0.03 and h = Dt = 0.05 for amp. = 0.09 & 3. Ref.
Amp.
Time
I1
Ours [10] [5] [8] (Dt = 0.03) [6]
0.3
80
1.20004 1.19998 1.1910 1.1964 1.23387
0.288 0.28798 0.2855 0.2858 0.29915
Ours [8] (h = 0.01) [6]
0.09
80
0.36001 0.3597 0.36665
0.02592 0.0259 0.02658
Ours [10] [6]
3
40
12.0 11.9999 11.9977
I2
28.8 28.80644 28.8784
L2 103
L1 103
0.0576 0.05759 0.05582 0.0569 0.06097
0.038815 0.056 3.849 7.444 24.697
0.051514 0.053 2.646 4.373 16.425
0.001555 0.00155 0.00161
0.010249 0.22 2.683
0.014128 0.16 1.836
I3
57.6 57.5983 57.5548
12.2944 17.189 37.365
7.01658 9.3597
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739
2735
5.2. Interaction of two EW solitary waves Here we study the interaction of two well separated EW solitary waves having different amplitudes and traveling in the same direction. The initial condition given by
uðx; 0Þ ¼
2 X
2
3ci sec h
1 pffiffiffi ðx xi Þ ; 2 d
i¼1
ð18Þ
where ci and xi, i = 1, 2 are arbitrary constants. The analytical values of the conservation laws in this case are as follows 2 pffiffiffi X I1 ¼ 12 d ci ; i¼1
pffiffiffi 2 144 d X I2 ¼ c2i ; 5 i¼1
pffiffiffi 2 288 d X I3 ¼ c3i : 5 i¼1
ð19Þ
For computational work we chose d = 1, x1 = 10, x2 = 30, c1 = 0.4 and c2 = 0.2 with interval [0, 100]. The amplitudes are in the ratio 2:1. The analytical values of the invariants are I1 = 7.2, I2 = 5.76 and I3 = 4.1472. Numerical values of the three invariants have been computed and recorded in Table 5 for h = 0.2 and Dt = 0.1. The simulation is done up to t = 160. Changes of the invariants I1, I2 and I3 from their initial values are less than 0.047%, 0.022% and 0.005% respectively. Fig. 2(a)–(d) show the computer plot of the interaction of two solitary waves at different time levels.
Table 5 Invariants for two solitons interactions of amplitudes in the ratio 2:1; h = Dt = 0.1. t
I1
I2
I3
0 40 80 120 160
7.19976 7.20021 7.20021 7.20020 7.20021
5.75988 5.75987 5.75987 5.75998 5.76009
4.14711 4.14712 4.14712 4.14711 4.14713
Fig. 2. Interaction of two solitary waves at (a) t = 0, (b) t = 80, (c) t = 100, and (d) t = 160.
2736
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739
5.3. Interaction of three EW solitary waves For interaction of three solitary waves of different amplitudes and traveling in the same direction we take
uðx; 0Þ ¼
3 X
2
3ci sec h ½0:5ðx xi Þ;
ð20Þ
i¼1
where ci and xi, i = 1, 2, 3 are arbitrary constants. The analytical values of the invariants in this case are 3 pffiffiffi X I1 ¼ 12 d ci ; i¼1
pffiffiffi 3 144 d X I2 ¼ c2i ; 5 i¼1
pffiffiffi 3 288 d X I3 ¼ c3i : 5 i¼1
ð21Þ
For numerical computations, we chose d = 1, x1 = 10, x2 = 25, x3 = 40, c1 = 0.4, c2 = 0.2 and c3 = 0.1 with interval [0, 120]. Thus the amplitudes are in the ratio 4:2:1. The analytical values of the invariants are I1 = 8.4, I2 = 6.048 and I3 = 4.2048. Numerical values of the invariants are calculated using our scheme and reported in Table 6 for h = 0.2 and Dt = 0.1. The simulation is done up to t = 200. Changes of I1, I2 and I3 from their initial values are less than 0.045%, 0.063% and 0.05% respectively. Fig. 3(a)–(d) show the details of interaction of three solitary waves corresponding to the above set of parameters at different time levels.
Table 6 Invariants for three solitons interactions of amplitudes in the ratio 4:2:1; h = 0.1, Dt = 0.5. t
I1
I2
I3
0 40 80 120 160 200
8.39977 8.40020 8.40020 8.40020 8.40020 8.40021
6.04792 6.04789 6.04751 6.04791 6.04793 6.04792
4.20474 4.20470 4.20447 4.20473 4.20473 4.20475
Fig. 3. Interaction of three solitary waves at (a) t = 0, (b) t = 80, (c) t = 140, and (d) t = 200.
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739
2737
5.4. Maxwellian initial condition Evolution of a train of solitary waves of the EW equation has been studied using the Maxwellian initial condition
uðx; 0Þ ¼ Expððx 20Þ2 Þ;
ð22Þ
for various values of d. The computations are carried out and values of the invariants are reported in Table 6 for the cases d = 0.1, 0.05, 0.025 and 0.01 at fixed space–time steps h = 0.05 and Dt = 0.025, respectively. It is to be noted here that when d is reduced more, then more and more solitary waves are formed. A single solitary wave is formed when d = 0.1 as shown in Fig. 4(a). For d = 0.05, the Maxwellian initial pulse breaks up into a train of at least two solitary waves as shown in Fig. 4(b). For d = 0.025, three stable solitary waves are generated as shown in Fig. 4(c). Finally when d = 0.01 the Maxwellian initial condition has split up into four stable solitary waves as shown in Fig. 4(d). All these have been computed at t = 12 and the computed values are reported in Table 7. Changes of the invariantsI2 and I3 from their initial values are less than 0.06% and 0.035% respectively for all values d taken above but I1 does not change in this process. Invariants remain almost constant for larger values of d and changed a little more when we used smaller values of d. 5.5. Development of an undular bore Finally we study the development of an undular bore. For this problem we use the initial condition
h x x i 0 uðx; 0Þ ¼ 0:5u0 1 tanh ; d
ð23Þ
where u(x, 0) denotes the elevation of the water above the equilibrium at time t = 0. The boundary conditions in this case are u ? 0 as x ? 1 and u ? u0 as x ? 1. We choose a range 20 6 x 6 50, with space step h = 0.1 and time step Dt = 0.025 for a simulation with d = 0.16666667, x0 = 0 and d = 2. The development of the undular bore when d = 2 is shown in Fig. 5 for times up to t = 400. Table 8 contains the values of the invariants I1, I2, I3 and position and amplitude of the leading undulation for the present method. The quantities I1, I2, I3 are no longer constant but increase linearly throughout the simulation at the following rates workout by Gardner et al. [5]
Fig. 4. The Maxwellian initial condition at t = 12 with (a) d = 0.1, (b) d = 0.05, (c) d = 0.025, and (d) d = 0.01.
2738
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739 Table 7 Invariants of EW equation using the Maxwellian initial condition at t = 12. d
Time
I1
I2
I3
d = 0.1
0 3 6 9 12
1.77245 1.77245 1.77245 1.77245 1.77245
1.37864 1.37923 1.37880 1.37877 1.37885
1.02333 1.02355 1.02338 1.02336 1.02339
d = 0.05
0 3 6 9 12
1.77245 1.77245 1.77245 1.77245 1.77245
1.31598 1.31648 1.31619 1.31617 1.31612
1.02333 1.02356 1.02340 1.02339 1.02339
d = 0.025
0 3 6 9 12
1.77245 1.77245 1.77245 1.77245 1.77245
1.28464 1.28520 1.28492 1.28418 1.28474
1.02333 1.02357 1.02340 1.02337 1.02337
d = 0.01
0 3 6 9 12
1.77245 1.77245 1.77245 1.77245 1.77245
1.26585 1.26632 1.26599 1.26639 1.26567
1.02333 1.02330 1.02294 1.02295 1.02293
Fig. 5. The development of an undular bore at t = 400 for d = 2.
Table 8 Development of an undular bore with parameters d = 2, U0 = 0.1, h = 0.1, k = 0.025, 20 6 x 6 50. Time
I1
I2
I3
x
u
0 50 100 150 200 250 300 350 400
2.00000 2.24990 2.49940 2.74803 2.99568 3.24238 3.48811 3.73289 3.97672
0.190278 0.223597 0.256823 0.289861 0.322702 0.355344 0.387795 0.420050 0.452121
0.0185000 0.0222486 0.0259850 0.0296895 0.0333601 0.0369971 0.0406011 0.0441725 0.0477118
– 1.04339 3.70855 6.46989 9.34631 12.2749 15.2270 18.1945 21.1660
– 0.113954 0.156979 0.170509 0.174771 0.176732 0.177727 0.178227 0.178410
Z 1 d d 1 I1 ¼ udx ¼ u20 ¼ M1 ; dt dt 1 2 Z 1 d d 2 I2 ¼ ðu2 þ dðux Þ2 Þdx ¼ u30 ¼ M2 ; dt dt 1 3 Z 1 d d 3 4 3 I3 ¼ u dx ¼ u0 ¼ M 3 : dt dt 1 4
ð24Þ
T. Roshan / Applied Mathematics and Computation 218 (2011) 2730–2739
2739
According to the above formulae, the predicted theoretical values are M1 = 5.03 103, M2 = 6.66667 104 and M3 = 7.5 105. The calculated values of the rate of change of I1, I2 and I3 from Table 7 are found to be M1 = 4.9418 103, M2 = 6.5461 104 and M3 = 7.30295 105 respectively. These values are very close to the theoretical values. 6. Conclusion The Petrov–Galerkin method using linear hat function and quadratic B-spline function as test and trial functions respectively was applied to study the solitary waves of the EW equation and it is shown that our scheme is conditionally stable and comparatively accurate than other methods available in the literature. We tested our scheme through single solitary wave in which the analytic solution is known and then extended to study the interaction of solitons where no analytic solution is known. It is also observed that conservation laws were reasonably satisfied for the interaction of two and three solitary waves. The Maxwellian initial condition was also used. Lastly the development of an undular bore was studied. The results indicate that the rates of increase in the invariants I1, I2 and I3 are in good agreement with the theoretical values. Acknowledgement I hereby acknowledge Dr. K.S. Bhamra for extending helping hand in this work and also offer sincere thanks to the reviewer for giving valuable comments and suggestions. References [1] D.H. Peregrine, Calculations of the development of an undular bore, J. Fluid Mech. 25 (2) (1966) 321–330. [2] T.B. Benjamin, J.L. Bona, J.J. Mahony, Model equations for long waves in non-linear dispersive systems, Phil. Trans. Roy. Soc. London A 272 (1972) 47– 78. [3] P.J. Morrison, J.D. Meiss, J.R. Carey, Scattering of RLW solitary waves, Physica 11D (1984) 324–336. [4] L.R.T. Gardner, G.A. Gardner, Solitary waves of equal width equation, J. Comput. Phys. 101 (1992) 218–223. [5] L.R.T. Gardner, G.A. Gardner, F.A. Ayoub, N.K. Amein, Simulation of EW undular bore, Commun. Numer. Methods Eng. 13 (1997) 583–592. [6] A. Dogan, Application of Galerkin’s method to equal width wave equation, J. Appl. Math. Comput. 160 (2005) 65–76. [7] S.I. Zaki, Solitary waves induced by the boundary forced EW equation, Comput. Math. Appl. Mech. Eng. 190 (2001) 4881–4887. [8] S.I. Zaki, A least squares finite element scheme for the EW equation, Comput. Math. Appl. Mech. Eng. 189 (2000) 587–594. [9] D.J. Evan, K.R. Raslan, Solitary waves for the generalized equal width (GEW) equation, Int. J. Comput. Math. 82 (4) (2005) 445–455. [10] I. Dag˘, B. Saka, A cubic B-spline collocation method for EW equation, Math. Comput. Appl. 9 (3) (2004) 381–392. [11] P.J. Olver, Euler operators and conservation laws of the BBM equation, Math. Proc. Camb. Phil. Soc. 85 (1979) 143–159. [12] J.M. Sanz-Serna, I. Christie, Petrov–Galerkin methods for nonlinear dispersive waves, J. Comput. Phys. 39 (1993) 94–103. [13] I. Christie, D.F. Griffths, A.R. Mitchell, J.M. Sanz-Serna, Product approximation for nonlinear problems in finite element methods, IMA J. Numer. Anal. 1 (1981) 253–266.