Discrete Singular Convolution Method for Numerical Solutions of Fifth Order Korteweg-De Vries Equations

A new computational method for solving the fifth order Korteweg-de Vries (fKdV) equation is proposed. The nonlinear partial differential equation is discretized in space using the discrete singular convolution (DSC) scheme and an exponential time integration scheme combined with the best rational approximations based on the Carathéodory-Fejér procedure for time discretization. We check several numerical results of our approach against available analytical solutions. In addition, we computed the conservation laws of the fKdV equation. We find that the DSC approach is a very accurate, efficient and reliable method for solving nonlinear partial differential equations.


Introduction
The study of travelling wave solutions of nonlinear partial differential equations (PDEs) is the major subject in many fields of physical and nonlinear sciences.Concepts like solitons, peakons, kinks, breathers, cusps and compactons have entered into various branches of natural sciences such as chemistry, biology, mathematics, communication and particularly in almost all branches of physics like the fluid dynamics, plasma physics, field theory, nonlinear optics and condensed matter physics.Among these nonlinear PDEs there exists an important class of the fifth order Korteweg-de Vries equations 2 3 2 5 0, where us a better understanding of the erratic and often unpredictable nature of natural phenomena, and soliton theory helps explain natural phenomena that are surprisingly predictable and regular even under conditions that would normally destroy such properties.A soliton is a solitary wave which preserves its shape and velocity after nonlinearly interacting with other solitary waves or (arbitrary) localized disturbances.
In general, Equation (1) does not admit exact solutions, therefore one has to resort to numerical methods.Due to the fifth-order terms in these equations, it is very difficult to compute the solutions of these equations accurately and efficiently.Recently, Shen [5] proposed a new dual-Petrov-Galerkin method for the third and higher oddorder equations.His approach was proven to be very effective for the KdV type equations in bounded domains [5] and in semi-infinite intervals [6].In [7], a numerical scheme based on the dual-Petrov-Galerkin method was proposed and implemented for the Kawahara and modified Kawahara equations.This class includes the well-known Lax [1], Sawada-Kotera (SK) [2], Kaup-Kupershmidt (KK) [3] and Ito [4] equations.The knowledge of close form solutions of nonlinear PDEs facilitates the verification of numerical solvers, aids physicists to better understand the mechanism that governs the physic models, provides knowledge to the physical problem, provides possible applications and helps mathematicians in the stability analysis of solutions.While strange attractors and chaos theory give In this paper, we propose a discrete singular convolution method to solve fifth order Korteweg-de Vries equations.Discrete singular convolution (DSC) methods belong to the family of local spectral (LS) methods.They were proposed by Wei [8] as a potential approach for numerical realization of the Hilbert transform, Abel transform, Random transform and Delta transform.The DSC algorithm has been essential to many practical applications, such as nonlinear equations [9], structural analysis [10,11], compressible and incompressible fluid flows [12,13], electromagnetic wave propagation, scattering [14,15] and image analysis [16].
Recently, Pindza and Maré [17] utilized a combined fourth order exponential time differencing of Adams type and the DSC method to solve the generalized Kortewegde Vries.Their approach revealed exponential convergence.The advantage of the DSC methods is that they exhibit exponential convergence of spectral methods [18] while having the flexibility of local methods for complex boundary conditions [10,19].
The discretization of the generalized Korteweg-de Vries equations in space with the DSC method yields a system of ordinary differential equations (ODE) that needs to be solved by time integration methods.We use the fourth order exponential time differencing Runge Kutta (ETDRK4) [20] for the solution of the resulting semidiscrete equations.The matrix exponential required by the scheme is efficiently computed using best rational approximations based on the Carathéodory-Fejér (CF) procedure [21].
The layout of this paper is as follows.We describe the formulation of the DSC method in Section 2. In Section 3, we discuss the exponential time integration methods for solving the semi-discrete system resulting from the spatial discretization of the nonlinear PDEs.Numerical results illustrating the merits of the new scheme are given in Section 4 and we present our conclusions in Section 5.

Discrete Singular Convolution Methods
Discrete singular convolution (DSC) methods are relatively new numerical techniques in the field of nonlinear equations.They are defined as follow.Consider a distribution, and an element of the space of test functions.A singular convolution can be defined by where is a singular kernel.For many science and engineering problems, an appropriate choice of has to be done.For instance, in the field of interpolation of surfaces and curves the singular kernel of delta type is very important.For numerical solutions of partial differential equations, the kernel is essential, where the subscript n denotes the -th order derivative of the distribution with respect to parameter , 0,1, n   n , x .While using the DSC method, numerical approximations of a function and its derivatives can be treated as convolutions with some kernels.According to the DSC method, the -th derivative of a function where is the grid spacing, k h x is the set of discrete grid points which are centered around x , and 2 1 M  is the effective kernel, or computational bandwidth; and is usually smaller than the whole computational domain.
In the present paper, we focus our attention on the regularized Shannon kernel (RSK) to provide discrete approximations to the singular convolution kernels of the delta type (3).The required derivatives of the DSC kernels can be easily obtained using ( [12]) The error estimation of the regularized Shannon kernel (RSK) delivers very small truncation errors when it uses the above convolution algorithm ( [23]) Theorem 2.1 (Qian [23]).Let Here is the number of grid points.The N L  error given by (6) decays exponentially with respect to the increase of the DSC band width .M The proof of the above theorem is beyond the scope of this paper.We refer the reader to [23] for a detailed discussion on the Shannon's sampling theorem.
Using (4) and ( 5), the entries of the first, second, third and fifth differentiation matrices , , and are given explicitly by Note that the differentiation matrix in ( 5) is in general banded.This gives rise to great advantage in large scale computations.Extension to higher dimensions can be realized by tensor products.
The choice of M ,  and was suggested by Qian and Wei [23].For instance, if the norm error is set to the following relations must be satisfied where r h   and is the frequency bound of the underlying function B f .To illustrate the procedure of discretization of PDEs by the DSC method, we consider the computation of fifth order KdV equations given by 2 3 2 5 0, where  .This equation was previously considered in [24] where its properties were studied and its analytical soliton solutions were revealed.In the present paper we mainly focus on numerical solutions of Equation ( 12) via the use of DSC method.
The semi-discretized version, at the th row, of the equation in consideration is obtained by substituting the relations (3) and ( 5) into (12), yielding 13) can be expressed in the following matrix form where represents the linear part of the system and represents the nonlinear part.
The main difficulty when dealing with systems of the type ( 14) is that the use of explicit time integrators is inefficient because the system typically suffers from instability due to the higher order derivative.This was emphasized by Pindza [25].Consequently, the time step size must be significantly reduced in order to fulfill the drastic stability condition present in explicit time integrators.In this paper we use the fourth order exponential time differencing Runge-Kutta method.

Exponential Time Differencing
Exponential time differencing (ETD) schemes are known for a long time in computational electrodynamics; see [26] for a comprehensive review of ETD methods and their history.In this section, we describe the exponential time differencing fourth-order Runge-Kutta (ETD4RK) method which was proposed by Cox-Matthews [27].
The main idea of the ETD methods is to multiply both sides of a differential equation by some integrating factor, then we make a change of variable that allows us to solve the linear part exactly and, finally, we use a numerical method of our choice to solve the transformed nonlinear part.

Overview of the Method
In order to elaborate on this approach, let us consider the following semi-linear partial differential equation where and are the linear and nonlinear operators, respectively.The semi-linear partial differential equation is discretized in space with the discrete singular convolution method.Therefore, we obtain a system of ordinary differential equations (ODEs) The exponential time differencing (ETD) methods can be obtained by integrating Equation ( 16) exactly between the time steps and with respect to t , to obtain There exist various ETD methods for the evaluation of (17).The purpose of this work is not to give a complete classification of ETD methods.We focus specifically on the fourth order exponential time differencing Runge-Kutta (ETDRK4) given by The main computational challenge in the implementation of exponential time differencing (ETD) methods is the need for fast and stable evaluations of exponential and related  -functions i.e., functions of the form   . The computation of these functions depends significantly on the structure and the range of eigenvalues of the linear operator and the dimensionality of the semi-discretized PDE.Unfor-tunately, for DSC methods the linear part have eigenvalues approaching zero, which leads to complications in the computation of the coefficients.Saad [28], and Hochbruck and Lubich [29] introduced Kyrlov methods to compute  -functions.Kassam and Trefethen [20] used Cauchy integral representation on a circle for a stable computation of  -functions.Our evaluation of expo- nential and related  -matrix functions follows the idea of Schmelzer and Trefethen [30].This method is based on computing optimal rational approximations to the matrix functions on the negative real axis using the Carathéodory-Fejér (CF) procedure [21], closely.The  - functions ( 19) can be computed explicitly by a recursive formula Another way to compute the functions j  is to use the Taylor series representation.Therefore, for all complex numbers , we have z However, it is known that the computation of these functions in their explicit or Taylor series form suffers from computational inaccuracy for matrices whose eigenvalues are equal to or approaching zero.This is generally the case when the spatial discretization is based on spectral methods.In order to overcome the numerical difficulties encountered in the computation of ( 20) and ( 21), a different tactic for evaluating the function was proposed in [20].The key idea is to approximate the functions (for matrices or scalars) by means of contour integrals in the complex plane where 2π If the size of the matrix is large, it is more advantageous to compute the product of the functions L   j z  and vectors rather than to compute explicitly.We have where s  and are the poles and the residues, respectively.The sum in (24) is evaluated by solving at c  most shifted linear systems.The poles and the residues are computed efficiently in standard precision by the Carathéodory-Fejér method [21,30].The stability region is four-dimensional, if both and L  are complex.The two-dimensional stability region is obtained if both and L  are purely imaginary or purely real, or if  is complex and is fixed and real.

Stability Analysis
In this section, we investigate the linear stability of the ETDRK4 method for the nonlinear autonomous system of ODEs, In the paper, we follow the analysis employed in [   where u is now the perturbation of and However, one can observe that the computation of 1 , 2 , 3 and suffers from computational inaccuracy for values of respectively, where is the number of interior points, N j u and j u u are the exact and computed values of the solution at point .j In this paper, we consider two case studies depending on the set of parameters of ( 25) that provide multi-soliton solutions.We evaluate the performance the DSC algorithms for different time increment , spatial discretization , the support size of DSC kernels t N M and regularization parameter  .
In our computation, the first set of parameters that we select are given by 5, 5, 5 In this case, the fifth order KdV Equation ( 11) is known as the Sawada-Kotera (SK) [2] equation and is given by 2 3 2 5 5 5 5 The SK (30) admits multi-soliton solutions [31].The derivation of these soliton solutions is beyond the scope of this paper.We only list them here for testing numerical procedures purposes.Single and two soliton solutions are given by where and In our computational work, we use the collocation points The SK equation possesses infinite conservation laws [31].The first three conservation laws are given as follow related to the mass, momentum and energy.The quantities 1 I , 2 I and 3 I are applied to measure the conservation properties of the collocation scheme, calculated by The second set of parameters are chosen as 10, 25, 20 . This is well-known as the Kaup-Kupershmidt (KK) [3] Multi-soliton solutions can be generated by the following nonlinear transformation of the dependent variable, For one soliton solution, the dependent variable function is given by For two soliton solutions, the dependent variable function is with The KK equation possesses infinite conservation laws [31], the first three are given as follows The quantities 1 I , 2 I and 3 I are applied to measure the conservation properties of the collocation scheme, calculated by In next sections, we study the propagation and the interaction of single and two soliton solutions, respectively.

Propagation of Single Solitons
In our numerical experiments, we first model the motion of a single soliton of the SK (30) and KK (38) equations.For the SK equation, the initial condition is taken from the exact solutions (32) and ( 31) at initial profile.Whereas for the KK equation, the initial condition is taken from the exact solutions (40) and (39) at initial profile.The boundary conditions in both cases are chosen so that In the first computation, we would like to investigate the convergence of the DSC method with respect to the number of grid points and the DSC bandwidth N M .The values of the parameters used in our numerical experiments are: .Figure 2 illustrates the convergence of the DSC with respect to the number of the grid points and the DSC bandwidth N M .We observe that numerical soliton solu- tions of the DSC method converge towards the exact soliton solutions as the number of grid points increases.We remark that the convergence of the DSC method also relies on the bandwidth N M .The results in Figure 2 shows that the case gives a better convergence, the case gives the worst convergence, whereas when we have an intermediate convergence.In fact when the bandwidth M is large, the DSC me- thod behaves like a global and detains exponential accuracy, whereas for a small value of M , the DSC behaves like a local method such as finite difference methods.This result is stated by Theorem 2.1.Figure 3 represents numerical propagation of one soliton solutions of the SK (a) and the KK (b) equations.These propagations occur for a long period of time with no spurious oscillations.
In the next experiment, we compute the error norms L  , 2 and conservation quantities 1 L I , 2 I and 3 I .The results are shown in Table 1 for one soliton solution of the SK equation and in Table 2 for one soliton solution of the KK equation.
From the numerical results given in Tables 1 and 2 it   I at a given time are equal to those of the initial value.Our scheme conserve, the mass, momentum and energy.t

Interaction of Two Solitons
This computational work is related to the interaction of two soliton solutions of SK (30) and KK (38) equations having different amplitudes and travelling in the same direction.For the SK equation, the initial condition is taken from the exact solutions (33) and ( 31) at initial profile; whereas for the KK equation, the initial condition is taken from the exact solutions (41) and (39) at initial profile.The boundary conditions in both cases are chosen so that To allow the interaction to occur, the experiment was run from 0 t  to in the region 400   100,100  We also investigate the convergence of the DSC method with respect to the number of the grid points and the DSC bandwidth N M as we did in the case of one soliton solutions.
All the results are shown in Figure 5.We observe that numerical soliton solutions obtained by means of the DSC method converge to the exact soliton solutions as the number of grid points increases.We also observe that the convergence of the DSC method relies on the bandwidth N M .The results on Figure 5 show that the case 64 M  gives a better convergence, the case 16 M  gives the worst convergence, whereas when

M 
we have an intermediate convergence.In fact when the bandwidth M is large, the DSC method be- haves like a global and detains exponential accuracy, whereas for a small value of M , the DSC behaves like a local method such as finite difference methods.
In addition, we compute the error norms , 2 and conservation quantities 1 L  L I , 2 I and 2 I are computed.The result are shown in Table 3 for two soliton solutions of the SK equation and in Table 4 for two soliton so-lutions of the KK equation.
From the numerical results given in Table 3 it is observed that throughout the simulation, the error norms L  and 2 are of magnitude at a long period of time , whereas the error norms and 2 (Table 4) are of magnitude  three quantities remain constant with respect to time.

Conclusion
We studied the application of the combined DSC scheme in space discretization and the ETDRK4 for time discretization to solve the SK and KK equations.We considered the case of the propagation of a single soliton and the interaction of two solitons.Numerical results showed that the DSC method converges exponentially with respect to the number of grid points and the bandwidth N M .Numerical checks on the conservation mass, mo- mentum and energy revealed that the three quantities remain constant with respect to time .The DSC scheme is a robust and reliable numerical method of the fifth order KdV equation.We are currently investigating the utility of the DSC method to solve the GRLW equation.t and  are real numbers.
or approaching zero.Therefore, it is important to make use of their Taylor expansions complex x plane where 1 r  .Hence, the boundary of the stability region is determined by writing

Figure 1 .
Figure 1.Stability regions in the complex x plane.The curves correspond to , from the outer curve to the inner curve respectively.The inner red curve corresponds to y = 0. ; ; ; ; ; 8 7 6 4 1 0.5 y        equation of the SK and KK equations.In each case, the soliton moves to the right across the space interval the DSC bandwidth M and the regular- izer parameter  is done according to the conditions (

Figure 2 .
Figure 2. Convergence the DSC method for the propagation of single soliton solution of the SK (a) and the KK (b) equations at 10 t  with , and 1 0.4 k  1 0   0.001 δt 

Figure 3 .
Figure 3. Propagation of single soliton solution of the SK (a) and the KK (b) equations with , 1 0.4 k 

4
shows the interaction of two soliton solutions of the SK (top) and KK (bottom) equations for , , 0,100 .It can be seen that the faster pulse interacts with and emerges ahead of the lower pulse with the shape and velocity of each soliton retained.

Figure 4 .
Figure 4. Interaction of two soliton solutions of the SK (top) and the KK (bottom) equations with , , 1 0.4 k 

8 10Figure 5 . 3 .
Figure 5. Convergence the DSC method for the interaction of two soliton solutions of the SK (a) and the KK (b) equations with , , 1 0.4 k 