This article is concerned with, the numerical solution of non-linear MKdV equation 
with the initial condition
and boundary conditions
where t and x denote the spatial and temporal variables, respectively, and 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 , 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. . Raslan and Baghdady  , 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  . A lumped Galerkin and Petrov-Galerkin methods were applied to solve the MKdV equation by AK et al.  .
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.
2. Numerical Methods
In this section, we derive a high order compact finite difference method for the initial boundary value problem (1) - (3). We first describe our solution domain and its grids. The solution domain is defined to be . Let and be uniform step size in the spatial and temporal directions respectively, also we denote for . Let and denote respectively, the exact and the numerical solution at the grid point . Using the generalized Pade approximation where the numerator and denominator of approximant are extended from polynomial functions to a series of any kind of functions or operators. The following Pade approximations of space derivatives are used  
and the shift operator E defined by
For more details see . 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 denotes the time derivative of U at ; 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.
2.1. Scheme 1 (Nonlinear Scheme)
The Crank Nicolson scheme for solving the first order differential system (5) can be displayed as
The derived scheme in (6) is a nonlinear implicit scheme (denoted by Scheme 1). In order to find the numerical solution , two iterative methods can be used, Newton’s method and fixed point method to solve the nonlinear penta-diagonal system obtained. We want to point out that Scheme 1, can be also obtained using Petrov-Galerkin with cubic spline as the test functions and linear splines as trial function  - .
2.2. 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
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 solution .
2.3. Scheme 3 (Nonlinear Implicit Scheme)
Another nonlinear implicit scheme for solving the first order differential system (5) is given by
This type of approximation is used in numerical solution of the coupled nonlinear Schrodinger equation (see Ismail    ). 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.
2.4. 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
By collecting the similar terms in (14), this will lead us to the linear penta-diagonal system
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 and , this can be easily obtained from the given initial condition at , 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 . and then, we apply the linearized scheme directly for each time step.
3. 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 1, and
for Scheme 3. We apply the two iterative Schemes (17) and (18) independently till the following condition
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.
4. Accuracy of the Proposed Schemes
To study the accuracy of Scheme 1, we replace the numerical solution by the exact solution in (6) to get the following equation
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 local truncation error (LTE)
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; . 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
5. Stability of the Scheme
In this section, we want to study the stability of the proposed schemes . Our stability analysis is based on the Von Neumann theory in which the growth factor of atypical Fourier mode defined as
where , is real number and in general complex. To implement the Fourier stability analysis, the MKdV equation needs to be linearized by assuming the nonlinear term in 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 equation
By solving Equation (30) for , we can get
It is easy to see that, , 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.
6. Numerical Results
To study the behavior of the derived schemes, rewrite Equation (1) as
with the homogenous boundary as 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 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
To demonstrate the efficiency and accuracy of the presented methods for solving the MKdV equation, we use the error norm
and the error norm
Two values of are used. The first case: and , in this case we test the efficiency of the proposed methods by conducting different experiments. The motion of single soliton and the interactions of two solitons will be considered. Trapezoidal rule is used to calculate the conserved quantities. In the second case, we choose and , and in this case, a comparison between our methods with .
6.1. The First Case: ( )
6.1.1. 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 Fixed point method and the cpu time for each method is 0.4 and 0.3, respectively. 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.
In Tables 3-6, we display the conserved quantities for Schemes 1, 2, 3 and 4. 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.
Figure 1. The evolution of the numerical solution of Scheme 1 with ε = 6, μ = 1, h = 0.1, k = 0.01, λ = 0.25.
Figure 2. The evolution of the numerical solution of Scheme 3 with ε = 6, μ = 1, h = 0.1, k = 0.01, λ = 0.25.
Figure 3. The evolution of the numerical solution of Schemes (2, 4) with ε = 6, μ = 1, h = 0.1, k = 0.01, λ = 0.25.
note that, Scheme 3 takes more time than the other method. In Scheme 1, three iterations are needed for convergence with tolerance 10−10, while Scheme 3 needs seven iterations for convergence.
Table 2. Error norms calculated of Schemes 1, 2, 3 and 4.
Table 3. Scheme 1 (Single soliton with ε = 6, μ = 1, h = 0.1, k = 0.01, λ = 0.25).
Table 4. Scheme 2 (Single soliton with ε = 6, μ = 1, h = 0.1, k = 0.01, λ = 0.25).
Table 5. Scheme 3 (Single soliton with ε = 6, μ = 1, h = 0.1, k = 0.01, λ = 0.25).
Table 6. Scheme 4 (Single soliton with ε = 6, μ = 1, h = 0.1, k = 0.01, λ = 0.25).
6.1.2. Rate of Convergence
To calculate the order of the proposed numerical schemes. We calculate the rate of convergence (RTC) using the formula
To calculate the rate convergence of Scheme 1 and Scheme 2 in space using (46), we choose , for different values of h. we calculate the 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 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.
6.1.3. Two Solitons Interaction
To study the interaction of two solitons, we choose the initial condition
and are arbitrary constants. To ensure an interaction of two solitary waves for Schemes 1, 2, 3 and 4, we choose the set of parameters , , , , , , and , over the interval
Table 7. Rate of convergence in space with k = 0.01, λ = 0.25.
Table 8. Rate of convergence in space with h = 0.1, λ = 0.25.
. 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, 2, 3 and 4. The interaction scenario is displayed in Figures 4-6.
Figure 4. Two solitons interaction for Scheme 1 with ε = 6, μ = 1, h = 0.1, k = 0.01, λ1 = 2.0, λ2 = 1.0.
Figure 5. Two solitons interaction for Scheme 3 with ε = 6, μ = 1, h = 0.1, k = 0.01, λ1 = 2.0, λ2 = 1.0.
Figure 6. Two solitons interaction for Schemes 2 and 4 with parameters ε = 6, μ = 1, h = 0.1, k = 0.01, λ1 = 2.0, λ2 = 1.0.
Table 9. Scheme 1: Interaction of two solitons (Conserved quantities).
Table 10. Scheme 2: Interaction of two solitons (Conserved quantities).
Table 11. Scheme 3: Interaction of two solitons (Conserved quantities).
Table 12. Scheme 4: Interaction of two solitons (Conserved quantities).
6.1.4. Three Solitons Interaction
In this test, we want to study the interaction scenario of three solitons. Accordingly, we choose the initial condition
subject to the homogenous boundary conditions. We have considered the Schemes 1, 2, 3 and 4 with the set of parameters , , , , , , , , and over the interval . In 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.
Figure 7. Three solitons interaction for Scheme 1 with ε = 6, μ = 1, h = 0.1, k = 0.01, λ1 = 2.0, λ2 = 1.0, λ3 = 0.5.
Figure 8. Three solitons interaction for Scheme 3 with ε = 6, μ = 1, h = 0.1, k = 0.01, λ1 = 2.0, λ2 = 1.0, λ3 = 0.5.
Figure 9. Three solitons interaction for Schemes (2, 4) with ε = 6, μ = 1, h = 0.1, k = 0.01, λ1 = 2.0, λ2 = 1.0, λ3 = 0.5.
Table 13. Scheme 1: Three solitons interaction (Conserved quantities).
Table 14. Scheme 2: Three solitons interaction (Conserved quantities).
Table 15. Scheme 3: Three solitons interaction (Conserved quantities).
Table 16. Scheme 4: Three solitons interaction (Conserved quantities).
6.2. The Second Case ( and )
In order to compare our results with , we choose the initial condition
The following parameters are used
The exact values of the conserved quantities in this case are
In Table 17, we display the errors for , 10 and 20. It is very easy to see that our methods are more accurate than . The simulation of the single soliton for all proposed methods are given in Figures 10-12.
Figure 10. The evolution of the numerical solution of Scheme 1 with ε = 3, μ = 1, h = 0.1, k = 0.01, λ = 0.845.
Figure 11. The evolution of the numerical solution of Scheme 3 with ε = 3, μ = 1, h = 0.1, k = 0.01, λ = 0.845.
Figure 12. The evolution of the numerical solution of Schemes 2 and 4 with ε = 3, μ = 1, h = 0.1, k = 0.01, λ = 0.845.
Table 17. Comparison of error norms for single soliton with ε = 3, μ = 1, h = 0.1, k = 0.01, λ = 0.845.
6.3. Breather Dynamics for MKdV Equation
In this section, we want to study the dynamics of the breather solution for the MKdV equation which conserve their energy during propagation. Sometimes breathers are called pulsating solitons because they propagate as isolated perturbations without losses. The exact breather solution for Equation (1) has the following form 
and p and q are arbitrary constants. A breather has nontrivial time periodic behavior. To study the behavior of the breather we choose our initial condition as
with the set of parameters
The exact values of the conserved quantities I1 and I2 are 0.0 and 8q, respectively. The conserved quantities are given in Table 18.
7. Birth of Solitons
7.1. Birth of Solitons Test 1
To study the birth of solitons, we choose the initial condition
together with set of parameters
Figure 13. The simulation of the breather solution.
Figure 14. Contours of the breather solution of the MKdV equation.
Table 18. Conserved quantities (breather solution with h = 0.0, k = 0.0001, p = 0.2, q = 1).
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.
7.2. 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) , we choose the set of parameters
Table 19. Brith of solitons test 1 with h = 0.1, k = 0.001.
Table 20. Brith of solitons test 2 with w = 0.25.
Table 21. Brith of solitons test 2 with w = 0.5.
Figure 15. Birth of ve. solitons using test 1.
Figure 16. Birth of solitons using test 2 with w = 0.25.
Figure 17. Birth of solitons using test 2 with w = 0.5.
In this test we have noticed, a birth of two solitons in case of and four solitons in case of . The conserved quantities are almost conserved. A very interesting thing we can anticipate the number of created solitons by using the formula
8. Concluding Remarks
In this work, we have derived four numerical methods for solving the MKdV equation. All methods are of fourth order accuracy in space and second order in time, the schemes are unconditionally stable. The exact solution and the conserved quantities are used to check the efficiency and the robustness of the proposed methods. The dynamics of the interaction of two and three solitons are discussed. Also we study the simulation of the breather solution of the MKdV, finally we study the birth of solitons using two different tests. We conclude that the schemes we have derived are highly accurate and can be easily generalized to solve any KdV like equation as well as the coupled KdV equation.
 Kaya, D. (2005) An Application for the Higher Order Modified KdV Equation by Decomposition Method. Communications in Nonlinear Science and Numerical Simulation, 10, 693-702.
 Wazwaz, A.-M. (2015) New (3+1)-Dimensional Nonlinear Evolution Equations with mKdV Equation Constituting Its Main Part: Multiple Soliton Solutions. Chaos, Solitons and Fractals, 76, 93-97.
 Ak, T., Karakoc, S.B.G. and Biswas, A. (2017) A New Approach for Numerical Solution of Modified Korteweg-de Vries Equation. Iranian Journal of Science and Technology, Transactions A: Science, 41, 1109-1121.
 Ak, T., Karakoc, S.B.G. and Biswas, A. (2017) Application of Petrov-Galerkin Finite Element Method to Shallow Water Waves Model: Modified Korteweg-de Vries equation, Scientia Iranica B, 24, 1148-1159.
 Mihaila, B., Cardenas, A., Cooper, F. and Saxena, A. (2010) Stability and Dynamical Properties of Rosenau-Hyman Compactons Using Padé Approximants. Physical Review E, 81, 056708.
 Cooper, F., Hyman, J.M. and Khare, A. (2001) Compacton Solutions in a Class of Generalized Fifth-Order Korteweg-de Vries Equations. Physical Review E, 64, 026608.
 Ismail, M.S. (2008) Numerical Solution of Complex Modified Korteweg-de Vries Equation by Petrov-Galerkin Method. Applied Mathematics and Computation, 202, 520-531.
 Ismail, M.S. (2009) Numerical Solution of a Coupled Korteweg-de Vries Equations by Collocation Method. Numerical Methods for Partial Differential Equations, 25, 275-291.
 Ismail, M.S. (2009) Numerical Solution of Complex Modified Korteweg-de Vries Equation by Collocation Method. Communication in Nonlinear Science and Numerical Simulation, 14, 749-759.
 Mashat, D.S., Wazzan, L.A. and Ismail, M.S. (2006) Petrov-Galerkin Method and K(2,2) Equation. International Journal of Computer Mathematics, 83, 331-343.
 Ismail, M.S. and Alamri, S.Z. (2004) Highly Accurate Finite Difference Method for Coupled Nonlinear Schrodinger Equation. International Journal of Computer Mathematics, 81, 333-351.
 Ismail, M.S. and Taha, T.R. (2001) Numerical Simulation of Coupled Nonlinear Schrodinger Equation. Mathematics and Computer Simulations, 56, 547-562.