Meshless Method of Lines for Numerical Solution of Kawahara Type Equations

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.


Introduction
In this paper we study numerical solution of Kawahara equation, Modified Kawahara equation and KdV-Kawahara equation respectively given as: 0, , 0 0, , 0 A variety of physical phenomena like, magneto acoustic waves in plasma [1], shallow water waves with surface 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 Benney-Lin equation [4,5].
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 dependent partial differential equations (PDEs).This method 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 numerical solution of PDEs.Other researchers like Fornberg et al. [19], Hon and Wu [20], Franke and Schaback [21], Wong et al. [22], Chen and Tanaka [23] also contributed a lot in this area.In order to solve the resultant collocation method Fasshauer modified the Kansa's method to a Hermite type collocation method [24].For nonlinear time dependent PDEs see also [25,26].In 1971 Hardy [27], developed MQ to approximate geogra-phical 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   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 radial basis functions (GSRBFs) approach face the problem of ill-conditioning of the dense interpolation matrix.In order to overcome these difficulties several alternatives including domain decomposition [33], preconditioning [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 smooth 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 arithmetic [43].The globally supported RBFs were used early on and for the problems whose size does not exceed 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 contains the numerical examples for the justification of the method and we conclude in Section 4.

Numerical Scheme
For implementation of numerical method, we consider the Kawahara Equation (1) with the following initial and boundary conditions ( To apply meshless method of lines, we first use radial basis functions to approximate space derivatives.The problem domain [a,b] where We have used the following RBFs where 1, 2, , j N   and c being shape parameter.The Equation (6) in the matrix form is  u Aλ (7) where Using Equation (7) in Equation ( 6) it follows that where Applying Equation (8) to Equation (1), and collocating at each node i x , we get system of first order ODEs where , In the similar fashion , , , , 1, 2, , In order to write the above system of equations in terms of column vectors, let Equation ( 9) then can be written as follows Rewriting Equation (10) as where  and the symbol (*) denotes component by component multiplication of two vectors.The initial condition is From the boundary conditions described in (5) we get Now we use classical fourth order Runge-Kutta scheme to solve Equations ( 11)-( 13), namely , 2 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 timediscretization operator".So the method is stable if all the eigenvalues of the Jacobian matrix in (11) lie inside the stability region R (i.e.  ) 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 interval 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.

Numerical Application
In this section the numerical results for Kawahara, modified Kawahara and KdV Kawahara equations are pre-sented.The root mean square norm and maximum error norms are calculated by the formulas For Kawahara equation the lowest three conserved quantities, defined in [49] are also calculated using with the following initial and boundary conditions The exact solution given in [10] is, where 1 1 2 13 k  For numerical computation we take [a,b] = [-20,30], 2, 1, 51 The simulation is carried out up to time t = 25.L 2 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 neighborhood 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 method using MQ are better than the polynomial based differential 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 The Point wise rate of convergence in space is calculated using the following formula:   with the following initial and boundary conditions , sech The above conditions are extracted from the exact solution given in [50],  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    4. The solitary wave profile at different time levels in comparison with the exact solution is shown in Figure 4.
Example 3.3 Consider the K 0 initial condition and boundary conditions are extracted from the exact solution given in [51].
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 L 2 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 increasing 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 solution ( 26) is shown in Figure 6, same as in [52].
Example 3.4 Considering Equation (1) for interaction of two positive solitary waves with the following initial condition We solve the problem by using MQ and GA RBFs taking    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: The two solitary waves propagate towards right as the time progresses.The process of interaction is sh n 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 We use multiquadric and Gaussian radial basis functions 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 parameters used in above equation are selected as:  8 we can see that the invariants of motion remain almost conserved as time increases.The variation in the three conserved quantities In this example we will show t n (1 dit Computational domain [-40,130] with h = 0.625 is considered.The scheme is run up to time t = 18.We take

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 solitary wave nd 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 invariants of motion remained conserved during the proc-ess of computation for all ca wo 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.

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.

Fig- ure 2
the forward motion of the solitary wave in comparison with exact solution (19) at different time levels is shown.

Figure 1 .
Figure 1.Error verses shape parameter c for Example 3.1.

Figure 2 .U
Figure 2. Travelling wave solution of Kawahara equation (solid lines showing numerical solution and dot (.) showing exact solution.

Figure 3 . 2 A
Figure 3. Error verses shape parameter c for Example 3.2 A and MQ are showing better accuracy than IMQ.The [a,b] = [-50,100] with N = 201.The calculation is carried

Figure 4 .
Figure 4. Solitary wave solution of modified Kawahara equation in comparison with exact solution (solid lines show numerical solution and dot (.) showing exact solution).
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 s aller waves and then the three waves separate from each other.The shape of the three solitary waves after collision is maintained.From Table