A Quartic Subdomain Finite Element Method for the Modified KdV Equation

In this article, we have obtained numerical solutions of the modified Korteweg-de Vries (MKdV) equation by a numerical technique attributed on subdomain finite element method using quartic B-splines. The proposed numerical algorithm is controlled by applying three test problems including single solitary wave, interaction of two and three solitary waves. To inspect the performance of the newly applied method, the error norms, L2 and L∞, as well as the four lowest invariants, I1, I2, I3 and I4 have been computed. Linear stability analysis of the algorithm is also examined.


Introduction
The modified Korteweg de-Vries (MKdV) equation which will be studied in this article is related to the following Korteweg de-Vries (KdV) equation The terms U U x and U xxx in the Eq.( 1) represent the nonlinear convection and dispersion, respectively.KdV equation is one of the main mathematical model for describing the theory of water waves in shallow channels.Some important physical phenomena for example propagation of long waves in shallow water waves, bubbleliquid mixtures, ion acoustic plasma waves and wave phenomena in enharmonic crystals can be described by the KdV equation which was first suggested by Korteweg and de Vries [1].The exact solutions of the equation obtained by [2,3].At first, Zabusky and Kruskal solved the KdV equation numerically using the finite difference method [4].Gardner et al. [5] showed the existence and uniqueness of solutions of the KdV equation.Many researches have used various numerical methods including finite difference method [6,7], finite element method [8][9][10][11][12][13][14][15][16][17][18], pseudospectral method [3] and heat balance integral method [19] to solve the equation.MKdV equation is a special case of the generalized Korteweg de-Vries (GKdV) equation having the form where p is a positive integer.In this study, we will consider the MKdV equation, a special form of (2) with the choice p = 2, ε = 3 and µ = 1,

610
A QUARTIC SUBDOMAIN FINITE ELEMENT METHOD with the homogeneous boundary conditions and an initial condition where t is time, x is the space coordinate and f (x) is a detected function.MKdV equation have a limited number of numerical studies in the literature.Kaya [20], used the Adomian decomposition method to obtain the higher order modified Korteweg de-Vries equation with initial condition.MKdV equation have been solved by using Galerkins' method with quadratic B-spline finite elements by Biswas et al. [21].Raslan and Baghdady [22,23], showed the accuracy and stability of the difference solution of the MKdV equation and they obtained the numerical aspects of the dynamics of shallow water waves along lakes' shores and beaches modeled by the MKdV equation.
A new variety of (3 + 1)-dimensional MKdV equation and multiple soliton solutions for each new equation were established by Wazwaz [24,25].A lumped Galerkin and Petrov Galerkin methods were applied to the MKdV equation by Ak et al. [26,27].
In the present study, we have developed Subdomain finite element method for the MKdV equation using quartic B-spline functions.Motion of a single solitary wave, interaction of two and three solitary waves are examined to show the performance and efficiency of the suggested method.We show that the suggested method is unconditionally stable applying the von-Neumann stability analysis.

Quartic B-splines
Quartic B-splines ϕ m (x), (m = −2, −1, ..., N, N + 1) at the knots x m are designated over the interval [a, b] by Prenter [28] Each quartic B-spline covers five elements, thus each element [x m , x m+1 ] is covered by five B-splines.A typical finite interval [x m , x m+1 ] is mapped to the interval [0, 1] by a local coordinate transformation defined by hξ = x − x m , 0 ≤ ξ ≤ 1.So quartic B-splines (6) in terms of ξ over [0, 1] can be given as follows: Stat., Optim.Inf.Comput.Vol. 6, December 2018 For the problem, the finite elements are identified with the interval [x m , x m+1 ].Using Eq.( 6) and Eq.( 7), the nodal values of U m , U ′ m , U ′′ m and U ′′′ m are given in terms of the element parameters δ m by and the variation of U over the element [x m , x m+1 ] is given by

Subdomain Method For The MKdV Equation
For the numerical algorithm, solution area of the problem is limited over an interval a ≤ x ≤ b.Let the partition of the space interval [a, b] into equally sized finite elements of length h at the points where δ m (t) are time dependent parameters and will be defined from the both boundary and weighted residual conditions.Applying the numerical approach to Eq.( 3) with the weight function we get the following equation When we use the hξ = x − x m transformation into the weak form of Eq.( 12) and integrating it term by term with some regulation by parts, guides to h 5 where the dot denotes differentiation with respect to t and If time parameters δ i and its time derivatives δi in Eq.( 13) are discretized by the Crank-Nicolson formula and usual forward difference approach respectively, we derive a repetition relationship between two time levels n and n + 1 relating two unknown parameters where The system (16) involves of N linear equations containing (N + 4) unknown coefficients (δ −2 , δ −1 , . . . ,δ N , δ N +1 ) T .We need four additional restraints to obtain a unique solution for this system.These are obtained from the boundary conditions (4) and can be used to remove δ −2 , δ −1 , δ N and δ N +1 from the systems (16) which occures a matrix equation for the N unknowns d n = (δ 0 , δ 1 , . . ., δ N −1 ) T of the form This matrix equation have been solved by using a variant of the Thomas algorithm.However, two or three inner iterations are implemented to the term at each time step to overcome the non-linearity caused by Z m .Before the solution process begins iteratively, the initial vector d 0 must be established by using the initial condition and following derivatives at the boundaries; So we have the following matrix form for the initial vector d 0 ;

Stability Analysis
To investigate the stability analysis of the suggested algorithm, it is properly use Von-Neumann theory based on Fourier analysis.Supposing the quantity U 2 in the nonlinear term U 2 U x is locally constant.Substituting the Fourier mode δ n m = ξ n e iσmh ,(i = √ −1) into the form of (16) we attain, where σ is mode number, h is the element size, θ = σh If we simplify the Eq.(19), Stat., Optim.Inf.Comput.Vol. 6, December 2018 S.B.G.KARAKOC

613
According to the Fourier stability analysis, for the given scheme to be stable, the condition |ξ| < 1 must be satisfied.Using a symbolic programming software or using simple calculations, since a 2 + b 2 = a 2 + (−b) 2 it becomes evident that the modulus of |ξ| is 1.Therefore the linearized scheme is unconditionally stable.

Numerical Results and Discussion
In this part, to realise the correction of our algorithm, some numerical exercises have been calculated: the motion of single solitary wave whose analytical solution is known and extended the scheme to the study of two and three solitary waves, whose analytical solutions are unknown during the interaction.The initial boundary value problem (3) − (5) possesses following conservative quantities [29,30]; To calculate the difference between exact and numerical solutions at some specified times, L 2 and L ∞ error norms have been used.

The motion of single solitary wave
For this problem, Eq.( 3) is evaluated with the boundary conditions U → 0 as x → ±∞ and the initial condition where µ and A is amplitude, k is the width of the single solitary wave.The exact solution of the MKdV equation is given as where c and x 0 are arbitrary constants.For this problem, the exact values of the invariants can be given as [21] For this test problem, we have taken the parameters ε = 3, µ = 1, h = 0.1, ∆t = 0.01, c = 0.845 and 0.3 through the interval 0 ≤ x ≤ 80, so the solitary waves have amplitudes 1.3 and 0.7746, respectively.The numerical algorithms are actuated to t = 1 to get the error norms L 2 , L ∞ and conserved quantities I 1 , I 2 , I 3 and I 4 .Invariants and error norms provided by the suggested method are given in Table (1) and Table (2).From these tables, it is obviously seen that the computed values of invariants are in good agreement with their analytical values.Solitary wave profiles are demonstrated at different time levels in Fig. (1) in which the soliton moves toward the right having nearly unchanged properties for example speed and amplitude.

Interaction of two solitary waves
As a second problem, we have studied the behavior of the interaction of two solitary waves having different amplitudes and travelling in the same direction.Initial condition of two well-seperated solitary waves of different amplitudes has the following form: (j = 1, 2) c j and x j are arbitrary constants.We have chosen the parameters ε = 3, µ = 1, h = 0.1, ∆t = 0.01, c 1 = 2, c 2 = 1, x 1 = 15 and x 2 = 25 over the interval 0 ≤ x ≤ 80.The run of the algorithm is carried up to time

Interaction of three solitary waves
As a last problem, we have considered the behavior of the interaction of three solitary waves having different amplitudes and traveling in the same direction.We examined the interaction of three solitary waves by using the initial condition together with boundary conditions U → 0 as x → ±∞.We have conceived the problem with parameters ε = 3, µ = 1, h = 0.1, ∆t = 0.01, c 1 = 2, c 2 = 1, c 3 = 0.5, x 1 = 15, x 2 = 25 and x 3 = 35 over the interval 0 ≤ x ≤ 80.
The experiment is run from t = 0 to t = 5 and values of the quantities are listed in Table (4).Table (4) indicates that invariants are nearly constant as the time increases.The behavior of the interaction of three solitary waves denotes at different times in Figure (3).

Conclusion
In this article, we have applied a quartic B-spline subdomain finite element method to the MKdV equation.Three different test problems have been solved.To demonstrate the efficiency of numerical scheme, the error norms L 2 , L ∞ and conserved quantities I 1 , I 2 , I 3 and I 4 have been calculated for the test problems.According to the tables in the paper, one can has easily seen that our error norms are enough small and the obtained invariants are acceptable in good agreement with their exact values.Our numerical scheme is also unconditionally stable.So, we can say our numerical algorithm is a reliable method for getting the numerical solutions of the physically important non-linear partial differential equations.