Applied Mathematics, 2011, 2, 608-618
doi:10.4236/am.2011.25081 Published Online May 2011 (http://www.SciRP.org/journal/am)
Meshless Method of Lines for Numerical Solution
of Kawahara Type Equations
Nagina Bibi, Syed Ikram Abbas Tirmizi*, Sirajul Haq
Faculty of Engineering Sciences, Ghulam Ishaq Khan Institute, Topi, Pakistan
Received Feburary 22, 2011; revised March 26, 2011; accepted March 30, 2011
Abstract
In this work, an algorithm based on method of lines coupled with radial basis functions namely meshless
method of lines (MMOL) is presented for the numerical solution of Kawahara, modified Kawahara and KdV
Kawahara equations. The motion of a single solitary wave, interaction of two and three solitons and the
phenomena of wave generation is discussed. The results are compared with the exact solution and with the
results in the relevant literature to show the efficiency of the method.
Keywords: Kawahara Type Equations, Meshless Method of Lines (MMOL), Radial Basis Functions (RBFs)
1. Introduction
In this paper we study numerical solution of Kawahara
equation, Modified Kawahara equation and KdV-Kawa-
hara equation respectively given as:
0,, 0
tx xxxxxxxx
uuu uuaxbt  (1)
20,, 0
tx xxxxxxxx
uuu uuaxbt  (2)
0,, 0
t xxxxxxxxxx
uu uuuuaxbt  (3)
A variety of physical phenomena like, magneto acous-
tic waves in plasma [1], shallow water waves with sur-
face tension [2] and capillary gravity water waves [3],
are described by Kawahara Equation (1) and modified
Kawahara Equation (2). KdV-Kawahara Equation (3) is
used to describe the one dimensional evolution of small
but finite amplitude long waves in various problems in
fluid dynamics. This equation is a specific form of Ben-
ney-Lin equation [4,5].
Different analytic and numerical methods including
the tanh-function method [6], Adomian decomposition
method [7], sine-cosine method [8], variational iteration
method, homotopy perturbation method [9], Crank-Ni-
colson Differential quadrature algorithms [10], Predic-
tor corrector methods [11], Dual-Petrov Galerkin method
[12] and RBF collocation method [13] have been pro-
posed for solving the Kawahara type equations.
We shall use method of lines approach [14,15] using
RBFs for numerical solution of the above problems using
radial basis functions. The method of lines [16] is pow-
erful and comprehensive approach for solving time de-
pendent partial differential equations (PDEs). This me-
thod involves two steps, first the approximation of spatial
derivatives converting the PDEs to a system of ordinary
differential equations (ODEs) and in the second step the
resulting system of ODEs is integrated in time. One of
the salient features of the MOL is the use of existing and
well established numerical methods for ODEs. Coupling
method of lines with radial basis functions makes it very
simple by getting rid of mesh generation as compared to
conventional mesh based methods. In 1990 Kansa [17,
18] for the first time used radial basis functions for nu-
merical solution of PDEs. Other researchers like Forn-
berg et al. [19], Hon and Wu [20], Franke and Schaback
[21], Wong et al. [22], Chen and Tanaka [23] also con-
tributed a lot in this area. In order to solve the resultant
collocation method Fasshauer modified the Kansa’s me-
thod to a Hermite type collocation method [24]. For non-
linear time dependent PDEs see also [25,26]. In 1971
Hardy [27], developed MQ to approximate geogra- phi-
cal surfaces. In Franke’s [28] review paper, the MQ was
ranked the best in some thirty scattered data interpolation
methods. For solving PDEs the convergence proofs in
applying the RBF’s is given by Wu [29]. The accuracy of
RBFs methods depends on the choice of a parameter
called shape parameter involved in infinitely smooth
RBFs like, Multiquadric (MQ), Gaussian (GA), Inverse
multiquadric (IMQ) radial basis functions. Tarwater [30]
found that the Root Mean Square of error decreases up to
certain limit and then increases rapidly when
.c

N. BIBI ET AL.609
Golberg, Chen, and Karur [31] and Hickernell and Hon
[32] applied the technique of cross validation to obtain
an optimal value of the shape parameter.
However the methods based on globally supported ra-
dial basis functions (GSRBFs) approach face the prob-
lem of ill-conditioning of the dense interpolation matrix.
In order to overcome these difficulties several alterna-
tives including domain decomposition [33], precondi-
tioning [34] and use of compactly supported RBF [35]
have been introduced. Another useful approach in this
regard is locally supported RBFs instead of globally
supported RBFs. The readers are recommended to visit
Wu and Liu [36], C. K. Lie et al. [37], R. Vertnik and
Bozidar Sarler [38]. It is worth mentioning that global,
infinitely differentiable RBFs typically interpolate sm-
ooth data with spectral accuracy [39-42] and the shape
parameter can be adjusted with the number of centers in
order to produce a interpolation matrix which is well
conditioned enough to be inverted in finite precision ari-
thmetic [43].The globally supported RBFs were used
early on and for the problems whose size does not ex-
ceed 400-500 data points. These methods should still be
the method of choice [44].
This paper is organized in three sections. In Section 2
the numerical scheme is explained and Section 3 con-
tains the numerical examples for the justification of the
method and we conclude in Section 4.
2. Numerical Scheme
For implementation of numerical method, we consider
the Kawahara Equation (1) with the following initial and
boundary conditions
 
,0 ,,uxf xaxb (4)
 


1
,,,uatg tubtgt

2
.
(5)
To apply meshless method of lines, we first use radial ba-
sis functions to approximate space derivatives. The problem
domain [a,b] is divided into nodes ,1,2,,
i
x
iN Out
of these points1
i are interior points
d and
,2,,xi N
,1an
i
x
iN are the boundary points.

,uxtThe approximate solution of is given by
   
T
1
,,
N
nn
jj
j
uxtu xttxx 

 
λ, (6)
where
 
T
T
12
,,,
N
x
xx x   


We have used the following RBFs
 
22,
jj
x
xxc MQ  





2
exp ,
jj
x
cx xGA  
 
22
1,
j
j
x
IMQ
xx c
 

where 1, 2,,jN

and c being shape parameter.
The Equation (6) in the matrix form is

uAλ (7)
where
  
TT
12 12
,,,, ,,,
NN
ututu t    

uλ
and







 
 
11 211
T
112 222
T
2
T
12
...
....
.... ...
....
N
N
N
NNNN
xx x
xxx x
x
x
xx x
 
 
 






























A



Using Equation (7) in Equation (6) it follows that


T1
,() ()uxtxx

Au Su, (8)
where




T1
xx


S
A
Applying Equation (8) to Equation (1), and collocating
at each nodei
x
, we get system of first order ODEs
     
d0,
d
1, 2,,
i
ix ixxx ixxxxxi
uux xx
t
iN



SuS uSu
(9)
where


ii
ut u









 
12
,,, ,1,2,,
,,1,2, ,
xixi xiNxi
jx iji
x
SxSxS xiN
xSxijN
x








S
S
In the similar fashion






12
,,,
1, 2,,
xxx ixxx ixxx iNxxx i
xSxSx Sx
iN
,

 





S
 
3
3,,1,2, ,
jxxx ij i
Sx SxijN
x











12
,,,
1, 2,,
xxxxxixxxxxixxxxxiNxxxxx i
xS xSxSx
iN
,

 





S
 
5
5,,1,2,,
jxxxxx ij i
Sx Sxij
x


N
In order to write the above system of equations in
terms of column vectors, let
T
123 1
... ,
NN
uu uuu


 


U

,
xjxi
N
N
Sx 
  S

N. BIBI ET AL. 610


,
,
xxxjxxxi NN
xxxxxjxxxxx i
N
N
Sx
Sx


 
 
S
S
Equation (9) then can be written as follows
 
d
dx xxx xxxxx
t 0
UUSU SU SU (10)
Rewriting Equation (10) as

d
dG t
UU (11)
where
  
x xxx xxxxx
G UUSUSUSU


and the sy-
mbol (*) denotes component by component multiplica-
tion of two vectors.
The initial condition is
   
00 0
012
,,,
N
tuxuxux

 U (12)
From the boundary conditions described in (5) we get
 


11 2
and N ut gtut gt

 (13)
Now we use classical fourth order Runge-Kutta sch-
eme to solve Equations (11)-(13), namely



1234
1
12 1
324
2( )
6
,2
,
2
nn
nn
nn
tKK KK
t
KG KGK
t
3
K
GKKGt

 











  


UU
UU
UUK

N
(14)
The RK4 scheme does not face problem of stability as
long as the time step Δt is chosen sufficiently small (see
Collatz [45]) which is so in our case. The rule of thumb
in selection of the time step Δt for stability is as follows
[46]: “The method of lines is stable if the eigenvalues of
the (linearized) spatial discretization operator, scaled by
time step Δt, lie in the stability region of the time-
discretization operator”. So the method is stable if all the
eigenvalues of the Jacobian matrix in
(11) lie inside the stability region R (i.e.

, 1,2,,
jj
 
2.78 0
j
t


)
of RK4 scheme. For further details regarding stability of
Runge-Kutta fourth order method, see Lambert [47] and
Jain [48]. As far as the selection of the shape parameter c
is concerned in this work, first we have to find an inter-
val for c in which matrix A of radial basis functions is
invertible and then select a value from that interval
which gives us the most accurate results.
3. Numerical Application
In this section the numerical results for Kawahara, modi-
fied Kawahara and KdV Kawahara equations are pre-
sented. The root mean square norm and maximum error
norms are calculated by the formulas
  
2
2
2
1
N
ex numexnum
j
Luuxuju j

 
 (15)
 
max
ex numexnum
j
Luuuju j

  (16)
For Kawahara equation the lowest three conserved
quantities, defined in [49] are also calculated using
2
12
2232
3
1
d, d,
2
131051 d
8962
xx
IuxIux
x
I
uuu 

 








x
Example 3.1
Consider the Kawahara equation
0
txxxxxxxxx
uuu uu

 
with the following initial and boundary conditions



4
105
, 0sech
169 o
uxkx x


 (17)


,0;,uatubt 0

 (18)
The exact solution given in [10] is,

4
105 36
, sech
169169o
uxtk xtx



 


 (19)
where 11
213
k
For numerical computation we take [a,b] = [–20,30],
2,1, 51
o
xxN

. The simulation is carried out up
to time t = 25. L2 and L∞ norms are calculated at t = 0, 5,
15 and 25, using MQ and GA radial basis functions with
the value of shape parameter found to be in neighbor-
hood of 4 as shown in Figure 1 and similarly for GA the
value is found to be in neighborhood of 0.3. We have
searched the optimal value of the shape parameter by
plotting maximum error verses shape parameter with step
0.01.The three conserved quantities are also shown in the
Table 1.The amplitudes and peak position of the solitary
waves are also calculated. The results with present me-
thod using MQ are better than the polynomial based dif-
ferential quadrature (PDQ) method [10] and are very
close to cosine expansion based differential quadrature
(CDQ) method [10]. While the results obtained by GA
are better than both methods mentioned in [10]. In Fig-
ure 2 the forward motion of the solitary wave in com-
parison with exact solution (19) at different time levels is
shown.
The Point wise rate of convergence in space is calcu-
lated using the following formula:

N. BIBI ET AL. 611
Table 1. Results for Kawahara equation in comparison with [10].
Method Time 3
210L
3
10L 1
I 2 I 3
I Height Peak
Position
CPU
time(s)
MQ(c = 4.3)
GA(c = 0.27)
PDQ[10] (Δt = 0.1)
CDQ[10]
0
5
15
25
0
5
15
25
0
5
15
25
0
5
15
25
0
0.09468
0.15362
0.16818
0
0.10075
0.10113
0.13160
0
1.986
2.543
2.851
0
0.151
0.156
0.159
0
0.04669
0.05939
0.04660
0
0.034297
0.03830
0.03990
0
0.921
1.045
0.863
0
0.043
0.049
0.076
5.97359
5.97348
5.97343
5.97355
5.973599
5.973662
5.973675
5.973532
5.97357
5.97060
5.97014
5.97353
5.97357
5.97372
5.97364
5.97350
1.27250
1.27250
1.27250
1.27250
1.272502
1.272502
1.272502
1.272502
1.27250
1.27250
1.27250
1.27250
1.27250
1.27250
1.27250
1.27250
–0.16458
–0.16458
–0.16458
–0.16458
–0.16458
–0.16458
–0.16458
–0.16458
–0.16458
–0.16458
–0.16458
–0.16457
–0.16458
–0.16458
–0.16458
–0.16458
0.62130
0.62119
0.62038
0.61880
0.621301
0.621201
0.620382
0.618802
0.62130
0.62102
0.62047
0.61872
0.62130
0.62122
0.62037
0.61877
2
3
5
7
2
3
5
7
2
3
5
7
2
3
5
7
0.094
0.188
0.328
0.171
0.172
0.281
1 2 345 67
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Sha pe parame t er for MQ
error
Figure 1. Error verses shape parameter c for Example 3.1.
Figure 2. Travelling wave solution of Kawahara equation
(solid lines showing numerical solution and dot (.) showing
exact solution.



1
10
10 1
log
log
ii
hh
ii
uU uU
hh



where and i
h
U represents the exact solution and the
numerical solution respectively and hi is spatial step size.
We calculate spatial rate of convergence by keeping time
step
u
t0.001

 fixed and varying the number of collo-
cation points (F = 20, 40, 80). From the Table 2 we can
see that the order of convergence decreases with the
smaller spatial step size. In all numerical examples we
have used MQ and GA in order to calculate order of
convergence.
Example 3.2
Consider the modified Kawahara equation
20
txxxxxxxxx
uuuuu

 
with the following initial and boundary conditions

 
2
, 0sechux Dkx
 
(20)



  
2
2
, sech;
, sech
uatDka Bt
ubtDkb Bt



(21)
The above conditions are extracted from the exact so-
lution given in [50],



2
,sechuxtDkx Bt


 (22)
where 311 ,,
25 25
10
DkB 4
neighborhood of 0.5 and 0.0001. The L2 and L∞ norms cal-
culated at t = 0, 5, 15, 25 are shown in Table 3.
The calculation is carried out by taking [a,b] = [–30,30]
with x = 1.We use MQ, GA and IMQ radial basis with
shape parameter found to be in the neighborhood of 3 for
MQ as shown in Figure 3 and for GA and IMQ it is in the

N. BIBI ET AL. 612
en
N Order
Table 2. Spatial rate of convergce at for Example 3.1, t = 25.
L∞ order L2
MQ
2.8100 10–2
9.1143
8.2164
1.0180 10–1
9.2509
8.6358
20
40
80
GA
20
40
80
1×
5.07036×10–5
–5
3.56316×10
1.4909 ×10–2 4
5.01268×10–5
4.77379×10–5
0.5089
0.0704
8×
1.67092×10–4
–4
1.28499×10
5.3877 –2 5×10
1.35443×10–4
1.25773×10–4
0.3788
0.1068
ults for ed Kawation.
method time(s)
Table 3. ResModifihara equa
time L∞ L2 I1 I
2 CPU
M
GA(c 0.43)
IMQ(c = 0.0001)
6.1995×10–5
5.3996×10–5
1.806 10–1
1.7896×10–4
1.0804×10–4
5.2570×10–1
–8.4 2.68
2.
0.18
0.
0.
Q(c = 2.9)
=
0
5
15
25
0
5
15
25
0
5
15
25
0
1.0717×10–4
1.2130×10–4
0
8.6124×10–5
7.8371×10–5
0
59×
5.18281×10–1
7.63711×10–1
0
2.7337×10–4
3.4855×10–4
0
1.9575×10–4
2.5257×10–4
0
5.2819×10–1
2.2681×100
8525
–8.48524
–8.48487
–8.48464
–8.48525
–8.48510
–8.48472
–8.48442
–8.48525
–8.48525
–8.48134
–8.47434
328
2.683176
2.68296
2.68275
2.68328
2.68317
2.68296
2.68275
683281
2.683281
2.683307
2.683351
7
0.313
0.453
171
0.313
0.543
203
0.343
0.469
11.5 22.5 33.5 44.5 55.5 6
0
0. 0 5
0.1
0. 1 5
0.2
0. 2 5
0.3
0. 3 5
0.4
0. 4 5
0.5
Shape paramet er for MQ
error
Figure 3. Error verses shape parameter c for Example 3.2
A and MQ are showing better accuracy than IMQ. The
dV-Kawahara equation
with the following initial and boundary conditions
.
G
order of convergence in space decreases with increasing
N as shown in Table 4. The solitary wave profile at dif-
ferent time levels in comparison with the exact solution
is shown in Figure 4.
Example 3.3
Consider the K
0
txxxxx xxxxx
uuuu uu 
 
4
105 1
,0 sechuxx x


169213 o


(23)

4
1051205
, sech
169 169
213 o
uxaat x



 



4
1051205
, sech
169 169
213 o
uxbbt x



 


 (25)
initial condition and boundary conditions are extracted
from the exact solution given in [51].

4
1051205
, sech
169 169
213 o
uxtxt x




(26)
The calculation is carried out by tak


 
ing [a,b] = [0,200]
with Δx = 1. The discrete root mean square error
and maximum error norm L∞ are calculated us
GA a
norm L2
ing MQ,
nd IMQ for time t = 1 up to 5. From the results
shown in Table 5 we can see that the both MQ and GA
are showing very good agreement with the exact solution.
Shape parameter verses error plot for MQ is shown in
Figure 5. The spatial rate of convergence is shown in
Table 6. The order of convergence decreases by in-
creasing collocation points for a fixed time step Δt =
0.001. The forward movement of the solitary wave at
different time levels in comparison with the exact solu-
tion (26) is shown in Figure 6, same as in [52].
Example 3.4
Considering Equation (1) for interaction of two posi-
tive solitary waves with the following initial condition
 
224
1
, 0sech4
ii
i
ux Ax x





.


(24)
i
A 

We solve the problem by using MQ and GA RBFs taking
[a,b] = [–50,100] with N = 201. The calculation is carried

N. BIBI ET AL. 613
Table 4. Spatial rate of converge
N L∞ order
nce at for Example 3.2, t = 5.
L2 order
MQ
24
48
–1
2.422 10–5
–2
–1
5.670 –5
–2
7.6180
96
GA
24
48
96
2.56860 × 10
1.49043 × 10–3 7.4291
5.20872×10
2.65141×10–3
85 ×
3.47477 × 10
–4
5.00876 × 10
2.17409 × 10–5
5.9428
6.1163
4.5259
16×10
9.95424×10
–4
9.19460×10
5.46571×10–5
5.5472
6.7583
4.0723
Tab Resle 5.ults for Kdwahara eq
thod L2 I1 plitudeCPU (s)
V Kauation at t = 5.
me L∞ I2 Amtime
MQ(c = 2.6)
GA(c = 0.2)
IMQ(c = 0.0001)
3.7697×1 9 141
3.390
3.641
1.0977 × 10–4
4.7924 × 10–6
5.0111 × 10–1
1.8465×10–-5
1.5920×100
5.97368
5.95790
1.27250
1.27250
0.621200
0.619041
0–4 5.9735.272500.62123.266
Tab ratevergt forple le 6.atial Sp of cone ance Exam 3.3, t = 5.
N L∞ order L2order
MQ
40
80
1.568 –4 3.731 –4
4.8553
160
GA
40
80
160
–2
3.12550×10
1.62584×10–3 4.2648
1.33333×10
4.60613×10–3
26×10
2.6983×10
–4
–3
1.4710×10
1.2344×10–5
2.1471
4.1971
3.5749
–1
04×10
1.1620×10
–3
–2
3.2768×10
2.4624×10–4
2.7242
4.7639
3.7341
Figure 4. Solitary wave solution of modified Kawahara
equation in comparison with exact solution (solid lines show
numerical solution and dot (.) showing exact solution).
out up to time t = 50 by taking time step Δt = 0.001. For
our numerical calculations the values of the parameters
involved in above equation are chosen as:
121 2
4108
0, 20,,,
105
xxA A 


 
The two solitary waves propagate towards right as the
time progresses. The process of interaction is shn in
Figure 7. During this process the larger wave catches up
th
ow
e smaller one and then the both waves separate from
each other maintaining their original shape. From Table
7, we can see that the invariants of motion remain almost
conserved as time increases. The variation in the three
1234567
0
0.2
0.4
1.4
0.6
0.8
1
1.2
error
S hape p aram eter for MQ
Figure 5. Error verses shape parameteric for Example 3.3.
conserved quantities is found to be in the range:
onsider
Equation (1) with the initial condition
1
23
1
MQ 40.5092540.48389,
45.8361445.85093, 32.3721932.15991,
GA 40.5092540.41284,
I
II
I

 

23
45.8361445.84364, 32.3721832.14082.II 
Example 3.5
For interaction of three solitary waves we c
 
324
1
, 0sech4
i
ii
i
A
uxAx x



 




N. BIBI ET AL. 614
We use multiquadric and Gaussian radial basis func-
tions in our numerical simulations to solve this problem.
The spatial domain is selected as [a,b] = [–30,120] with
Δx = 0.75. The calculation is carried out up to time t =
50 taking time step Δt = 0.001. The values of the pa-
rameters used in above equation are selected as:
123 123
41086
20, 0, 20,,,,
105
xxx AAA 



 
The three solitary waves propagate towards right. The
process of interaction is shown in Figure 8. During this
rocess the taller wave moves faster and catches up the
m
p
saller waves and then the three waves separate from
each other. The shape of the three solitary waves after
collision is maintained. From Table 8 we can see that the
invariants of motion remain almost conserved as time
increases. The variation in the three conserved quantities
Figure 6. Solitary wave motion of KdV Kawahara equation
(solid lines show numerical solution and dot (.) showing
exact solution).
-50 050 100
0
0. 5
1
1. 5
2
2. 5
3
3. 5
4
x
u(x,t)
t=0
GA
MQ
-50 050 100
0
0. 5
1
1. 5
2
2. 5
3
3. 5
4
4. 5
x
u(x,t)
t=25
GA
MQ
t
= 0
t
= 25
u (x,t)
u (x,t)
Figure 7. Interaction of two solitons for Kawahara equation Example 3.4.
Table 7. Invariants for interaction of two solitons for Example 3.4.
MQ GA
–50
–50
time I1 I2 I3 I1 I
2 I 3
0
10
20
30
40
50
40.509259
40.507987
40.499950
40.5207361
40.550
45.836141
45.839240
45.8427481
45.8359623
45.8443065
–32.37219
–32.12856
–32.73918
–32.58956
–32.1
40.509259
40.492695
40.466653
571742
40.579624
45.836141
45.837775
45.839353
45.840373
45.841964
–32.37218
–32.12820
–32.73746
–32.59560
–32.12278
2.14082
5240
40.4838916 45.8509384 –32.1599140.412842 45.843648 –3
0125
40.
is found to be in the range:
Example 3.6
he phenomena of wave
generation for Equatio). We consider the following
initial conion
1
II
2
3
12
3
MQ 51.4726351.56859,51.10049.15838,
33.63806 33.06973,
GA 51.4726351.47389,51.1004951.10452,
33.63806 33.39361.
I
II
I
 

 

51
In this example we will show t
n (1
dit


41
, 0sech213 o uxx x 




Computational domain [–40,130] with h = 0.625 is
considered. The scheme is run up to time t = 18. We take
x
x

N. BIBI ET AL.615
10

. The sing solilitee s
tary waves. The serves, hosi
of three waves are calcariovel
shown in Table W oflead
itsoci from the other
oo re n in three con-
ies i be wing range:
Table 8. Invariants for interaction of three solitons for Example 3.5.
MQ GA
tary waves. The serves, hosi
of three waves are calcariovel
shown in Table W oflead
itsoci from the other
oo re n in three con-
ies i be wing range:
Table 8. Invariants for interaction of three solitons for Example 3.5.
MQ GA
letary wave sps in to throli-wavoli- wav
conconed quantitied quantitieight and peight and ption twtion tw
ulated at vulated at vus time leus time les as servs as serv
9. 9. ith passageith passage time the time the ing ing
e owing toe owing to faster vel faster velty gets farty gets far
waves as sh waves as shwn in Figu wn in Figu 9. Variatio9. Variatio
ed quantited quantits found tos found toin the folloin the follo
time I1 I2 I3 I1 I
2 I 3
0
10
20
30
40
50
51.472631
51.631982
51.688762
51.664151
51.601572
51.568591
51.100493
51.094986
51.103862
51.127187
51.14098
51.15838
–33.63806
–33.40011
–34.23846
–35.7
–33.3
–33.06973
51.472631
51.835962
51.4
51.100493
51.098240
–33.63806
–33.39015
1904
8561
51.848746
51.690220
51.493752
73898
51.099504
51.101129
51.102653
51.104529
–34.21111
–35.15871
–33.83892
–33.39361
Table 9. Invariants for interaction o
S
f three solitons for Example 3.6.
econd wave Third wave Leading wave
time Height Position Height Position Height Position
MQ
0
2
4
8
12
18
10
13.990228
13.645838
14.024270
1
13.
0
12.5
22.5
43.125
-
7.98375
8.471020
8.347924.375 3.427390 7.5
12.5
20
3.93545
90723
63.75
94.375
8.48316
8.34579
36.25
55
3.49294
3.50220
3
-
5.625
11.875
-
2.964111
3.247767
-
–0.625
1.875
5
-40 -200
4
20 4060 80100120
0
0.5
1
1.5
2
2.5
3
3.5
x
u(x,t)
t=0t=0 GA
MQ
-40
4.5
-20 20 4006080100 120
0
0.5
1
1.5
2
2.5
3
3.5
4
GA
x
u
t=2
(x,t)
5t=25
MQ
-40 -20 02040 60 80 100 120
0
0. 5
1
1. 5
2
2. 5
3
3. 5
4
4. 5
x
u(x,t)
t=40t=40
GA
MQ
-40 -20 020406080 100 120
0
0. 5
1
1. 5
2
2. 5
3
3. 5
4
4. 5
x
u(x,t)
t=60t=60
GA
MQ
ree s
t
= 0
Figure 8. Interaction of tholitons for Example 3.5.
u (x,t)
u (x,t)
u (x,t)
u (x,t)
x
x
x
x
–40 –20 –40 –20
t
= 25
t
= 40
t
= 60
–40 –20
–40 –20

N. BIBI ET AL. 616
-20 020 40 60 80 100
0
1
2
3
4
5
6
7
8
9
10
MQ
GA
t=0
-20 02040 60 80 100
-2
0
2
4
6
8
10
12
14
MQ
GA
t=1
t
= 1
t
= 0
–20
–2
–20
-20 020 40 6080100
-2
0
2
4
6
8
10
12
14
MQ
GA
t=2
14
-20 020 406080100
-2
0
2
4
6
8
10
MQ
GA
12 t=18
Figure 9. Generation of waves for Example 3.6.
t
= 18
t
= 2
12
3
12
3
MQ 96.143654365,329.65039
329.65039, 875.41604701.42681,
GA96.14365 96.14365,329.65039
329.65039, 875.41612706.85789.
96.1
I
I
I
I
I
I
 

 

4. Conclusions
In this paper, we have used method of lines coupled with
radial basis functions for numerical solution of Kawahara
type equations. The numerical results describing motion
of single solitary wave, interaction of two and three soli-
tary wavend phenomena of wave generation have
been discussed. The accuracy of the solution depends
upon the choice of the shape parame
lected experimentally. The numerical results using MQ
and GA for Kawahara equation are better than the Crank
Nicolson differential quadrature algorithm [10]. The in-
variants of motion remained conserved during the proc-
ess of computation for all cawo main advantages of
this method are mesh less property and use of ODE
solvers of high quality and their codes to approach the
solution of PDEs. Also the presented method is simple,
easy to implement because no mesh is required in the
problem domain. Only radial distance between the nodes
is used to approximate the solution.
5. Acknowledgements
The first author is thankful to HEC Pakistan for financial
help through Grant No. 063-112367-Ps3-096. We are
also thankful to the reviewers for their constructive
ts.
6. References
[1] T. Kawahara, “Oscillatory Solitary Waves in Dispersive
Media,” Journal of the Physical Society of Japan, Vol. 33,
s a
ter, which has been
se
ses. T
commen
–20
–2 –2
–2 0

N. BIBI ET AL.617
No. 1, 1972, pp. 260-264. doi:10.1143/JPSJ.33.260
[2] T. Bridges and G. Derks, “Linear Instability of Solitary
Wave Solutions of the Kawahara Equation and itsen-
eralizations,” Society for Industrial and Applied Mathe-
matics: Journal on Mathematical Analysis, Vol. 33, No.
6, 2002, pp. 1356-1378.
036141099361494
d J. Scheurle, “Existence of Perturbed
Solitary Wave Solutions to a Model Equation for Water
and Modified Kawa-
hara Equations by Sine-Cosine Method,” Chaos, So
4, 2008, pp. 1193-1197.
riational Iteration Method and
kmaz and I. Dag, “Crank-Nicolson-Differential
Quadrature Algorithms for the Kawahara Equat
s and Fractals, Vol. 42, No. 1, 2009, pp
1016/j.chaos.2008.10.033
1, 2008, pp. 48-63
10, pp.
991.
-K
. 3-5, 2002, pp.
20)48:8<1187::AI
, Vol. 93, No. 1,
pendent
e,” Computers
dary-Only RBF
cs with Applica-
rsity Press, Proceeding of Chamonix, Nash-
uation,” En-
a Class of Nonlinear
neering Analysis with
G
doi:10.1137/S0
[3] J. K. Hunter an
240
Waves,” Physica D, Vol. 32, No. 2, 1988, pp. 253-268.
doi:10.1016/0167-2789(88)90054-1
[4] D. J. Benny, “Long Waves on Liquid Film,” Journal of
Mathematical Physics, Vol. 45, 1966, pp. 150-155.
[5] S .P. Lin, “Finite Amplitude Side-Band Stability of a
Viscous Film,” Journal of Fluid Mechanics, Vol. 63, No.
3, 1974, pp. 417-429. doi:10.1017/S0022112074001704
[6] E. Yusufoglu and A. Bekir, “Symbolic Computation and
New Families of Exact Travelling Solutions for the Ka-
wahara and Modified Kawahara Equations,” Computers
and Mathematics with Applications, Vol. 55, No. 6, 2008,
pp. 1113-1121. doi:10.1016/j.camwa.2007.06.018
[7] D. Kaya, “An Explicit and Numerical Solutions of Some
Fifth-Order KdV Equation by Decomposition Method,”
Applied Mathematics and Computation, Vol. 144, No.
2-3, 2003, pp. 353-363.
doi:10.1016/S0096-3003(02)00412-5
[8] E. Yusufoglu, A. Bekir and M. Alp, “Periodic and Soli-
tary Wave Solutions of Kawahara
[19]
litons
D-NME942>3.0.CO;2-K
and Fractals, Vol. 37, No.
9] L. Jin, “Application of Va[
Homotopy Perturbation Method to the Modified Kawa-
hara Equation,” Mathematical and Computer Modeling,
Vol. 49, No. 3-4, 2009, pp. 573-578.
doi:10.1016/j.mcm.2008.06.017
[10] A. Kor
ion,”
. and Mathematics with Applications, Vol. 37, No. 8, 1999,
pp. 23-43.
Chaos Soliton
65-73. doi:10.
[11] K. Djidjeli, W. G. Price, E. H. Twizell and Y. Wang,
“Numerical Methods for the Solution of the Third and
Fifth-Order Dispersive Korteweg de Vries Equations,”
Journal of Computational and Applied Mathematics, Vol,
58, No. 3, 1995, pp. 307-336.
doi:10.1016/0377-0427(94)00005-L
[12] J. M. Yuan, J. Shen, W. Jiahong, “A Dual Petrov-Galer-
kin Method for the Kawahara-Type Equations,” Journal
of Scientific Computing, Vol. 34, No. . vil
doi:10.1007/s10915-007-9158-4
[13] S. Haq and M. Uddin, “RBFs Approximation Method for
Kawahara Equation,” Engineering Analysis with Bound-
ary Elements, Vol. 35, No. 3, 2011, pp.575-580.
doi:10.1016/j.enganabound.2010.07.009
[14] Q. Shen, “A Meshless Method of Lines for the Numerical
Solution of KdV Equation Using Radial Basis Functions,”
Engineering Analysis with Boundary Elements, Vol. 33,
No. 10, 2009, pp. 1171-1180.
doi:10.1016/j.enganabound.2009.04.008
[15] S. Haq, N. Bibi, S. I. A. Tirmizi and M. Usman, “Mesh-
less Method of Lines for the Numerical Solution of Gen-
eralized Kuramoto-Sivashinsky Equation,” Applied
Mathematics and Computation, Vol. 217, No. 6, 20
4-2413.
[16] W. E. Schiesser, “The Numerical Method of Lines: Inte-
gration of Partial Differential Equations,” 1st Edition,
Academic Press; San Diego California, 1
[17] E. J. Kansa, “Multiquadrics a Scattered Data Approxima-
tion Scheme with Application to Computational Fluid Dy-
namics-I, ” Computers and Mathematics with Applications,
Vol. 19, No. 8-9, 1990, pp. 127-145.
doi:10.1016/0898-1221(90)90270-T
[18] E. J. Kansa, “Multiquadrics a Scattered Data Approxima-
tion Scheme with Application to Computational Fluid
Dynamics-II,” Computers and MaFFthematics with Ap-
plications, Vol. 19, No. 8-9, 1990, pp.147-161.
doi:10.1016/0898-1221(90)90271
B. Fornberg, T. A. Driscoll, G. Wright and R. Charles,
“Observations on the Behavior of Radial Basis Function
Approximations near Boundaries,” Computers and Ma-
thematics with Applications, Vol. 43, No
473-490.
[20] Y. C. Hon and Z. Wu, “A Quasi-Interpolation Method for
Solving Stiff Ordinary Differential Equations,” Interna-
tional Journal of Numerical Methods for Engineering,
Vol. 48, No. 8, 2000, pp. 1187-97.
doi:10.1002/(SICI)1097-0207(200007
[21] C. Franke and R. Schaback, “Solving Partial Differential
Equations by Collocation Using Radial Basis Functions,”
Applied Mathematics and Computation
1998, pp. 73-82.doi:10.1016/S0096-3003(97)10104-7
[22] A. S. M. Wong, Y. C. Hon, T. S. Li, S. L. Chug and E. J.
Kansa, “Multizone Decomposition of Time De
Problems Using the Multiquadric Schem
doi:10.1016/S0898-1221(99)00098-X
[23] W. Chen and M. Tanaka, “A Meshless Exponential Con-
vergence, Integration-Free and Boun
Technique,” Computers and Mathemati
tions, Vol. 43, No. 3-5, 2002, pp. 379-391.
doi:10.1016/S0898-1221(01)00293-0
[24] G. E. Fasshauer, “Solving Partial Differential Equations by
Collocation with Radial Basis Functions,” Nashville Van-
derbilt Unive
le, 1996.
[25] S. Ul. Islam, A. J. Khattak and I. A. Tirmizi, “A Meshfree
Method for Numerical Solution of KdV Eq
gineering Analysis with Boundary Elements, Vol. 32, No.
10, 2008, pp. 849-855.
[26] A. J. Khattak, S. I. A. Tirmizi and S. Ul. Islam, “Application
of Meshfree Collocation Method to
Partial Differential Equations,” Engi
Boundary Elements, Vol. 33, No. 5, 2009, pp. 661-667.
doi:10.1016/j.enganabound.2008.10.001

N. BIBI ET AL. 618
raph
f Geophysical Re-
. 181-200.
re and Applied Ma-
tudy of Hardy’s Multi-
ved
unction
hug and E. J.
f Time Dependen
E. J. Kansa, “A Least-Square Preconditioner
ise Polynomial, Positive Definite
n of
od (LRPIM) for
ultiquadric
ndary Value Problem
, “Meshless Local Radial Basis
500018411
tions,” Cambridge
er and J.
990, pp.
9-94.
ch-
erential Equa-
-Type Equation,” Bo-
tters A, Vol.
s Solutions for Several
[27] R. L. Hardy, “Multiquadric Equations of Topog
other Irregular Surfaces,” Journal o
y and [39] M. Buhmann and N. Dyn, “Spectral Convergence of
Multiquadric Interpolation,” Proceedings of the Edin-
burgh Mathematical Society, Vol. 36, No. 2, 1993, pp.
319-333.
search, Vol. 76, No. 8, 1971, pp. 1905-1915.
doi:10.1029/JB076i008p01905
[28] C. Franke, “Scattered Data Interpolation: Tests of Some
Methods,” Mathematics of Computation, Vol. 38, No.
157, 1982, pp
[29] Z. M. Wu, “Solving PDE with Radial Basis Function and
the Error Estimation,” In: Z. Chen, Y. Li, C. A. Micchelli,
Y. Xu and M. Dekker, Ed., Advances in Computational
Mathematics, Lecture Notes on Pu
thematics, Volume 202, GuangZhou, 1998.
[30] A. E. Tarwater, “A Parameter S
quadric Method for Scattered Data Interpolation,” Tech-
nical Report, Lawrence Livermore National Laboratory,
Livermore, 1985.
[31] M. A. Golberg, C. S. Chen and S. R. Karur, “Impro
211-
Multiquadric Approximation for Partial Differential Equ-
ations,” Engineering Analysis with Boundary Elements,
Vol. 18, No. 1, 1996, pp. 9-17.
doi:10.1016/S0955-7997(96)00033-1
[32] J. Hickernell, and Y. C. Hon, “Radial Basis F
Approximation of the Surface Wind Field from Scattered
Data,” International Journal of Applied Science and
Computations, Vol. 4, No. 3, 1998, pp. 221-247.
[33] A. S. M. Wong, Y. C. Hon, T. S. Li, S. L. C
Kansa, “Multizone Decomposition ot [46
Problems Using the Multiquadric Scheme,” Computers
and Mathematics with Applications, Vol. 37, No. 8, 1999,
pp. 23-43. doi:10.1016/S0898-1221(99)00098-X
[34] L. Ling and
tion
for Radial Basis Functions Collocation Methods,” Ad-
vances in Computational Mathematics, Vol 23, No. 1-2,
2005, pp. 31-54. doi:10.1007/s10444-004-1809-5
[35] H. Wendland, “Piecew
and Compactly Supported Radial Functions of Minimal
Degree,” Advances in Computational Mathematics, Vol.
4, No. 1, 1995, pp. 389-96. doi:10.1007/BF02123482
[36] Y. L. Wu and G. R. Liu, “A Meshfree Formulatio
Local Radial Point Interpolation Meth
goli
Incompressible Flow Simulation,” Computational Me-
chanics, Vol. 30, No. 5-6, 2003, pp. 355-365.
[37] C. K. Lee, X. Liu and S. C. Fan, “Local M
Approximation for Solving Bou
360,
s,”
F
Computational Mechanics, Vol. 30, No. 5-6, 2003, pp.
396-409. doi:10.1007/s00466-003-0416-5
[38] R. Vertnik and B. Sarler
Function Collocation Method for Convective-Diffusive
Solid-Liquid Phase Change Problems,” International
Journal of Numerical Methods for Heat & Fluid flow,
Vol. 16, No. 5, 2006, pp. 617-640.
doi:10.1108/09615530610669148
doi:10.1017/S0013091
[40] M. D. Buhmann, “Radial Basis Func
University Press, Cambridge, 2003.
doi:10.1017/CBO9780511543241
[41] W. R. Madych and S. A. Nelson, “Error Bounds for Mul-
tiquadric Interpolation,” In: C. Chui, L. Schumak
Ward Eds., Approximation Theory VI, Academic Press,
New York, 1989, pp. 413-416.
[42] W. R. Madych and S. A. Nelson, “Multivariate Interpola-
tion and Conditionally Positive Definite Functions, II,”
Mathematics of Computation, Vol. 54, No. 189, 1
230.doi:10.1090/S0025-5718-1990-0993931-7
[43] S. A. Sarra, “Adaptive Radial Basis Function Method for
Time Dependent Partial Differential Equations,” Applied
Numerical Mathematics, Vol. 54, No. 1, 2005, pp. 7
doi:10.1016/j.apnum.2004.07.004
[44] G. E. Fasshauer, “On the Numerical Solution of Differen-
tial Equations with Radial Basis Functions,” Boundary
Element Technology XIII, Waterford Institute of Te
nology Press, Waterford, 1999, pp. 291-300.
[45] L. Collatz, “The Numerical Treatment of Differential
Equations,” 3rd Edition, Springer, Verlag, 1966.
] J. D. Lambert, “Computational Methods in Ordinary Dif-
ferential Equations”, John Wiley and Sons, J. W. Arrow-
smith Ltd, Bristol, 1983.
[47] M. K. Jain, “Numerical Solution of Diff
s,” New Age International (P) Ltd., Publishers, New
Delhi, 1984.
[48] L. N Trefeten, “Spectral Methods in MATLAB,” Society
for Industrial and Applied Mathematics, 2000.
doi:10.1137/1.9780898719598
[49] R. P. Malik, “On Fifth Order KdV
ubov laboratory of theoretical physics, Joint Institute
for Nuclear Research, Dubna, 1997.
[50] A. M. Wazwaz, “New Solitary Wave Solutions to the
Modified Kawahara Equation,” Physics Le
No. 4-5, 2007, pp. 588-592.
doi:10.1016/j.physleta.2006.08.068
[51] A. M. Wazwaz, “Abundant Soliton
orms of the Fifth-Order KdV Equation by Using the
Tanh Method,” Applied Mathematics and Computation,
Vol. 182, No. 1, 2006, pp. 283-300.
doi:10.1016/j.amc.2006.02.047
[52] J. C. Ceballos, M. Sepulveda and O. P. V. Villagran,
“The Korteweg-de Vries-Kawahara Equation in a
Bounded Domain and Some Numerical Results,” Applied
Mathematics and Computation, Vol. 190, No. 1, 2007, pp.
912-936. doi:10.1016/j.amc.2007.01.107