Numerical Simulation of Modified Kortweg-de Vries Equation by Linearized Implicit Schemes

In this paper, we are going to derive four numerical methods for solving the Modified Kortweg-de Vries (MKdV) equation using fourth Pade approximation for space direction and Crank Nicolson in the time direction. Two nonlinear schemes and two linearized schemes are presented. All resulting schemes will be analyzed for accuracy and stability. The exact solution and the conserved quantities are used to highlight the efficiency and the robustness of the proposed schemes. Interaction of two and three solitons will be also conducted. The numerical results show that the interaction behavior is elastic and the conserved quantities are conserved exactly, and this is a good indication of the reliability of the schemes which we derived. A comparison with some existing is presented as well.


Introduction
This article is concerned with, the numerical solution of non-linear MKdV equation [1] where t and x denote the spatial and temporal variables, respectively, and ( ) , u x t is the unknown function, ε and µ are arbitrary constants. The MKdV equation admitted soliton solution which can be defined as a solitary wave solution, highly stable and retains its identity (shape and speed), upon interaction. MKdV equation has a limited number of numerical studies in the literature. Kaya [2], was used the Adomian decomposition method to obtain the higher order modified Korteweg de-Vries equation with initial condition. MKdV equation has been solved by using Galerkin's method with quadratic B-spline finite elements by Biswas et al. [3]. Raslan and Baghdady [4] [5], derived finite difference method and finite element method using collocation method with septic spline for solving the MKdV equation, they show that both schemes are unconditionally stable in the linearized sense. A new variety of (3 + 1)-dimensional MKdV equations and multiple soliton solutions for each new equation were established by Wazwaz [6] [7]. A lumped Galerkin and Petrov-Galerkin methods were applied to solve the MKdV equation by AK et al. [8] [9].
In this paper, we will derive four numerical schemes for solving the MKdV equation are presented; based on the Padé approximation of fourth order accuracy in space, together with Crank-Nicolson scheme of second order accuracy in the time direction. We obtain two nonlinear implicit schemes and two linearly implicit. The resulting system produced is a nonlinear penta-diagonal system in case of the nonlinear method, where Newton's method and Fixed point methods are used to solve these systems. To overcome this difficulty, the solution of the solution of the nonlinear systems, we introduced two linearized implicit schemes, which can be solved directly by using Crout's method. The implicit schemes are unconditionally stable according to the von-Nueumann stability analysis. We have studied the motion of a single solitary wave, interaction of two and three solitary waves to show the performance and efficiency of the suggested methods.
We will present a comparison between our methods and other research.
The rest of the paper is organized as follows. In Section 2, we present four different schemes using fourth order Pade approximation in space direction and second order in time direction using Crank Nicolson. Also, we use the linearized methods to avoid nonlinear term. In Section 3, we present an explanation of the algorithm of the fixed point method. In Section 4, the accuracy of the proposed schemes is given. In Section 5, the stability analysis of schemes is derived using von-Neumann stability analysis. In Section 6, we present various numerical tests which validate the accuracy and the efficiency of the proposed schemes, and the dynamics of the breather solution. In Section 7, birth of solitons is presented by using two tests. Finally, concluding remarks are given Section 8.

Numerical Methods
In this section, we derive a high order compact finite difference method for the initial boundary value problem (1) For more details see [12]. The Pade approximation is more accurate than the finite difference method and produced a highly compact difference scheme with small compact support. Now by using these approximations in MKdV Equations (1) -(3), we obtain the first order differential system in time where m U  denotes the time derivative of U at m x ; and the boundary conditions The first order differential system in (5) can be solved by Crank Nicolson method. It will lead us to a non-linear penta-diagonal system. This system can be solved by using any iterative methods, such as Newton's method or fixed point method. The system in (5) can also be solved by using linearization techniques.
It will give us a linear penta-diagonal system which can be solved by Crout's method directly. In next subsections, details of the proposed schemes will be discussed.

A Linearly Implicit Scheme (Scheme 2)
In the previous subsection, a nonlinear implicit finite difference method is obtained, and at each time step, we need to solve a nonlinear penta-diagonal system. To avoid this difficulty we will use linearization technique, Taylor series approximation of degree one for the nonlinear term in (6), to obtain the following formula By making use of (8) into Equation (6), we obtain the linear implicit scheme where 1 2 , , 1, 2, , The system in (9) now, is a linear penta-diagonal system which can be solved directly by Crout's method directly (no need for iterations) to find the numerical

Scheme 3 (Nonlinear Implicit Scheme)
Another nonlinear implicit scheme for solving the first order differential system where 1 2 , , 1, 2, , 1 12 2 This type of approximation is used in numerical solution of the coupled nonlinear Schrodinger equation (see Ismail [19] [20] [21]). The resulting system is a nonlinear penta-diagonal system. A fixed point method is used to solve this system, the details of this method will be discussed in the next section. The proposed scheme is of second order accuracy in time and fourth order in space.

Scheme 4 (Linearly Implicit Scheme)
As we have seen, the previous Scheme (11), leads us to a nonlinear penta-diagonal system, an iterative scheme is needed to get the numerical solution at each time step. This issue can be considered as a handicapped property of this method. To overcome this difficulty, we proposed a linearization technique for the nonlinear term in (11) in the following manner By using this approximation, we obtain the linearly implicit scheme where 3 2 , , 1, 2, , By collecting the similar terms in (14), this will lead us to the linear pen- The system in (16) is a linear penta-diagonal system which can be solved by Crout's method directly. Scheme (16) is a three time level scheme, and due to this we need two initial conditions ( ) ,0 u x and ( ) , u x t k = , this can be easily obtained from the given initial condition at 0 t = , and any two time level scheme of second order accuracy in time or higher and fourth order in space, like Scheme 1, Scheme 3, can be used to find the numerical solution at t k = . and then, we apply the linearized scheme directly for each time step.

Fixed Point Method
The nonlinear schemes (Scheme 1 and Scheme 3) we derived lead us to a nonlinear system, To get the solution, we need to build up an iterative scheme. To accomplish this job the fixed point method is used in the following manner for Scheme 3. We apply the two iterative Schemes (17) and (18) independently till the following condition n s n s The initial guess vector for Scheme 1 and Scheme 3 is given by We note (17) and (18) can be written in a matrix vector form as ( ) where H is a matrix of penta-diagonal structure of constant coefficients. To solve this system, we apply Crout's method by factoring the matrix H as a product of Lower and upper triangular matrices, at the beginning of the calculations, and then at each iteration, we need to solve lower and upper triangular systems which is very cheap and easy.

Accuracy of the Proposed Schemes
To study the accuracy of Scheme 1, we replace the numerical solution Taylor's expansion for all terms in Equation (21) about the grid point ( ) the following expressions are obtained By using these expressions into (21), and collecting similar terms, we will get The first four brackets (25) are zero by the MKdV equation, and the LTE will be reduced into So, Scheme 1 is of second order accuracy in time and fourth order in space; ( ) 2 4 O k h + . Similar analysis can be also applied for the other schemes. We conclude that all the derived schemes are of second order accuracy in time and fourth order accuracy in space. All schemes derived in this paper are consistent with MKdV equation, because

Stability of the Scheme
In this section, we want to study the stability of the proposed schemes [16]. Our u u is locally constant. Using this, we will get the linear version of (6) and this can be given as By substituting (27) into (28), we get after some manipulation the following It is easy to see that, e 1 k α ≤ , thus we can say that Schemes 1, 2, 3 and 4, are unconditionally stable according to Von-Neumann stability analysis. We conclude that all schemes are convergent because of their consistency and unconditional stability. The rate of convergence is calculated and agrees with the order of convergence of all methods, fourth order in space and second order in time.

Numerical Results
To study the behavior of the derived schemes, rewrite Equation (1) as with the homogenous boundary 0 u → as x → ±∞ and the initial condition The exact solution of the MKdV (33) is where A is amplitude and β is the width of the single solitary wave, they are defined as where , , ε µ λ and 0 x are arbitrary constants. The MKdV equation satisfies the conserved quantities: The exact values of the conserved quantities (37) -(39) using the exact solution (35) are

. Single Soliton
To test our numerical methods, we choose the initial condition subject to the homogenous boundary conditions. In order to generate the numerical solutions, the following parameters are used The exact values of the conserved quantities using (40) are In Table 1, we display the errors for Schemes 1 using Newton's method and The results show that two methods have the same results but Newton's method takes more time than Fixed point method. In Table 2, we display the errors for Schemes 1, 2, 3 and 4. The results show that the derived methods produced a highly accurate results.
The results show that the method conserves the conserved quantities exactly.
The simulation of the single soliton is given in Figures 1-3. The cpu time required to produce Tables 3-6 is 0.4, 0.2, 0.7 and 0.3 seconds, respectively. We Table 1. Error norms calculated of Scheme 1.
To calculate the rate convergence of Scheme 1 and Scheme 2 in space using (46), we choose 0.01 k = , for different values of h. we calculate the L ∞ error norm and the (RTC), we displayed the results in Table 7, the fourth order convergence rate in space is observed. To calculate the rate convergence in time for Scheme 1 and Scheme 2, we choose , for different values of k, we calculate the L ∞ error norm and hence the (RTC), we displayed the results in Table   8, second order convergence rate in time is observed. The same procedure can be applied for Schemes 3 and 4, to calculate the rate of convergence for space and time.

Two Solitons Interaction
To study the interaction of two solitons, we choose the initial condition  0 80 x ≤ ≤ . From the numerical results, we have found that the two solitons recover their shapes after the interaction and the computed conserved quantities are in a very good agreement with the exact ones. See Tables 9-12 for Schemes 1 Tables 13-16, we display the conserved quantities for Schemes 1 -4 respectively. All methods showed the conservation of the conserved quantities. In Figures 7-9, we display the interaction scenario of three solitons. We have noticed that the three solitons recover their original shapes after the interaction.
The following parameters are used The exact values of the conserved quantities in this case are In Table 17, we display the errors for 1 t = , 10 and 20. It is very easy to see that our methods are more accurate than [1]. The simulation of the single soliton for all proposed methods are given in Figures 10-12.       By using (56), we have noticed that the conserved quantities are almost conserved after the creation of five solitons and the results are presented in Table 19.
The simulation of this test is given in Figure 15.

Birth of Solitons Test 2
In this second test, we use the initial condition where w is a free parameter.
In order to study the behavior of the solution using the initial condition (58) , The results of the conserved quantities for 0.25 w = and 0.5 w = are given in Table 20 and Table 21. The simulation of this test is given in Figure 16 and     In this test we have noticed, a birth of two solitons in case of 0.5 w = and four solitons in case of 0.25 w = . The conserved quantities are almost conserved.
A very interesting thing we can anticipate the number of created solitons by using the formula ( )

Concluding Remarks
In this work, we have derived four numerical methods for solving the MKdV