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

M. S. Ismail, Fakhirah Alotaibi

Department of Math, King Abdulaziz University, Jeddah, Saudi Arabia.

DOI: 10.4236/am.2020.1111077                

Abstract

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.

Keywords

MKdV Equation, Pade Approximation, Nonlinear Numerical Schemes, Linearly Implicit Schemes, Fixed Point Method, Interaction of Solitons

1. Introduction

This article is concerned with, the numerical solution of non-linear MKdV equation [1]

$u_t+\epsilon u^2{u}_x+\mu u_{xxx}=0,\left(x,t\right)\in D=\left[x_l,x_r\right]\times \left(0\le T\right)$ (1)

with the initial condition

$u\left(x\mathrm{,0}\right)=u_0\left(x\right)\mathrm{,}x\in \left[x_l\mathrm{,}{x}_r\right]$ (2)

and boundary conditions

$u\left(x,t\right)=0\text{at}x=x_l,x_r$ (3)

where t and x denote the spatial and temporal variables, respectively, and $u\left(x,t\right)$ is the unknown function, $\epsilon$ and $\mu$ 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.

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 $\left\{\left(x\mathrm{,}t\right)\mathrm{<mo>\; |\; </mo>}{x}_l\le x\le x_r\mathrm{,0}\le t\le T\right\}$. Let $h=\frac{x_r-x_l}{M}$ and $k=\frac{T}{N}$ be uniform step size in the spatial and temporal directions respectively, also we denote $x_m=x_l+mh,t_n=nk$ for $m=0,1,2,\cdots ,M;n=0,1,2,\cdots ,N$. Let $u_m^n$ and $U_m^n$ denote respectively, the exact and the numerical solution at the grid point $\left(x_m\mathrm{,}{t}_n\right)$. 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 [10] [11]

$u_x\approx \frac{B\left(E\right)}{A\left(E\right)}+\frac{1}{480}{h}^6\frac{{\partial }^{7}u}{\partial x^7}\mathrm{,}{u}_{xxx}\approx \frac{C\left(E\right)}{A\left(E\right)}+\frac{1}{240}{h}^4\frac{{\partial }^{7}u}{\partial x^7}\mathrm{,}$ (4)

where

$A\left(E\right)=\frac{1}{120}\left[E^2+26E^1+66+26E^{-1}+E^{-2}\right]\mathrm{,}$

$B\left(E\right)=\frac{1}{24h}\left[E^2+10E^1-10E^{-1}-E^{-2}\right]\mathrm{,}$

$C\left(E\right)=\frac{1}{2h^3}\left[E^2-2E^1+2E^{-1}-E^{-2}\right]\mathrm{,}$

and the shift operator E defined by

$E^j{u}_m=u_{m+j},j=-2,-1,0,1,2$

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

$A\left(E\right){\stackrel{\dot{}}{U}}_m+\frac{\epsilon }{3}B\left(E\right)U_m^3+\mu C\left(E\right)U_m=0.$ (5)

where ${\stackrel{\dot{}}{U}}_m$ denotes the time derivative of U at $x_m$ ; and the boundary conditions

$U_{-1}=U_0=0,U_M=U_{M+1}=0$

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

$\begin{array}{l}A\left(E\right)\left[U_m^{n+1}-U_m^n\right]+p_{1}B\left(E\right)\left[{\left(U^3\right)}_m^{n+1}+{\left(U^3\right)}_m^n\right]\\ +p_{2}C\left(E\right)\left[{\left(U\right)}_m^{n+1}+{\left(U\right)}_m^n\right]=\mathrm{0,}\end{array}$ (6)

where

$p_1=\frac{k\epsilon }{6},p_2=\frac{k\mu }{2},m=0,1,2,\cdots ,M-1$ (7)

The derived scheme in (6) is a nonlinear implicit scheme (denoted by Scheme 1). In order to find the numerical solution ${\left[U_m^{n+1}\right]}_{m=1}^{M-1}$, 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 [13] - [18].

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

$\begin{array}{c}\frac{{\left(U^3\right)}_m^{n+1}+{\left(U^3\right)}_m^n}{2}=\frac{1}{2}\left[\left[{\left(U^3\right)}_m^n+3{\left(U^2\right)}_m^n\left(U_m^{n+1}-U_m^n\right)\right]+{\left(U^3\right)}_m^n\right]\\ =\frac{1}{2}\left[3{\left(U^2\right)}_m^n{U}_m^{n+1}-{\left(U^3\right)}_m^n\right].\end{array}$ (8)

By making use of (8) into Equation (6), we obtain the linear implicit scheme

$\begin{array}{l}A\left(E\right)\left[U_m^{n+1}-U_m^n\right]+p_{1}B\left(E\right)\left[3{\left(U^2\right)}_m^n{U}_m^{n+1}-{\left(U^3\right)}_m^n\right]\\ +p_{2}C\left(E\right)\left[{\left(U\right)}_m^{n+1}+{\left(U\right)}_m^n\right]=0\end{array}$ (9)

where

$p_1=\frac{kϵ}{6},p_2=\frac{k\mu }{2},m=1,2,\cdots ,M-1$ (10)

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 ${\left[U_m^{n+1}\right]}_{m=1}^{M-1}$.

2.3. Scheme 3 (Nonlinear Implicit Scheme)

Another nonlinear implicit scheme for solving the first order differential system (5) is given by

$\begin{array}{l}A\left(E\right)\left[U_m^{n+1}-U_m^n\right]+p_{1}B\left(E\right)\left(\left[{\left(U^2\right)}_m^{n+1}+{\left(U^2\right)}_m^n\right]\left[U_m^{n+1}+U_m^n\right]\right)\\ +p_{2}C\left(E\right)\left[{\left(U\right)}_m^{n+1}+{\left(U\right)}_m^n\right]=\mathrm{0,}\end{array}$ (11)

where

$p_1=\frac{kϵ}{12},p_2=\frac{k\mu }{2},m=1,2,\cdots ,M-1$ (12)

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.

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

${\left(U^2\right)}_m^{n+1}+{\left(U^2\right)}_m^n={\left(2U_m^n-U_m^{n-1}\right)}^2+{\left(U^2\right)}_m^n$ (13)

By using this approximation, we obtain the linearly implicit scheme

$\begin{array}{l}A\left(E\right)\left[U_m^{n+1}-U_m^n\right]+p_{1}B\left(E\right)\left(\left[{\left(2U_m^n-U_m^{n-1}\right)}^2+{\left(U^2\right)}_m^n\right]\left[U_m^{n+1}+U_m^n\right]\right)\\ +p_{2}C\left(E\right)\left[{\left(U\right)}_m^{n+1}+{\left(U\right)}_m^n\right]=0\end{array}$ (14)

where

$p_3=\frac{kϵ}{12},p_2=\frac{k\mu }{2},m=1,2,\cdots ,M-1$ (15)

By collecting the similar terms in (14), this will lead us to the linear penta-diagonal system

$\begin{array}{l}A\left(E\right)U_m^{n+1}+p_{1}B\left(E\right)\left[{\left(2U_m^n-U_m^{n-1}\right)}^2+{\left(U^2\right)}_m^n\right]U_m^{n+1}+p_{2}C\left(E\right)U_m^{n+1}\\ =A\left(E\right)\left[U_m^n\right]-p_{1}B\left(E\right)\left(\left[{\left(2U_m^n-U_m^{n-1}\right)}^2+{\left(U^2\right)}_m^n\right]\left[U_m^n\right]\right)-p_{2}C\left(E\right)U_m^n.\end{array}$ (16)

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 $u\left(x\mathrm{,0}\right)$ and $u\left(x\mathrm{,}t=k\right)$, this can be easily obtained from the given initial condition at $t=0$, 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.

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

$\begin{array}{l}\left[A\left(E\right)+p_{2}C\left(E\right)\right]U_m^{n+\mathrm{1,}s+1}\\ =A\left(E\right)U_m^n-p_{1}B\left(E\right)\left[{\left(U^3\right)}_m^{n+\mathrm{1,}s}+{\left(U^3\right)}_m^n\right]-p_{2}C\left(E\right)U_m^n\end{array}$ (17)

for Scheme 1, and

$\begin{array}{l}\left[A\left(E\right)+p_{2}C\left(E\right)\right]U_m^{n+\mathrm{1,}s+1}\\ =A\left(E\right)U_m^n-p_{1}B\left(E\right)\left(\left[{\left(U^2\right)}_m^{n+\mathrm{1,}s}+{\left(U^2\right)}_m^n\right]\left[U_m^{n+\mathrm{1,}s}+U_m^n\right]\right)-p_{2}C\left(E\right)U_m^n\end{array}$ (18)

for Scheme 3. We apply the two iterative Schemes (17) and (18) independently till the following condition

$\Vert U^{n+\mathrm{1,}s+1}-U^{n+\mathrm{1,}s}\Vert \le {10}^{-10}\mathrm{.}$

The initial guess vector for Scheme 1 and Scheme 3 is given by

$U^{n+\mathrm{1,}\left(0\right)}=U^n$

We note (17) and (18) can be written in a matrix vector form as

$H{U}^{n+\mathrm{1,}s+1}=G\left(U_m^n\mathrm{,}{U}_m^{n+\mathrm{1,}s}\right)$ (19)

where

$H=A\left(E\right)+p_{2}C\left(E\right)$ (20)

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 $U_m^n$ by the exact solution $u_m^n$ in (6) to get the following equation

$A\left(E\right)\left[\frac{u_m^{n+1}-u_m^n}{k}\right]+\frac{\epsilon }{3}B\left(E\right)\left[\frac{{\left(u^3\right)}_m^{n+1}+{\left(u^3\right)}_m^n}{2}\right]+C\left(E\right)\left[\frac{u_m^{n+1}+u_m^n}{2}\right]=0.$ (21)

Taylor’s expansion for all terms in Equation (21) about the grid point $\left(x_m\mathrm{,}{t}_n\right)$, the following expressions are obtained

$\frac{u^{n+1}-u^n}{k}=\frac{\partial u}{\partial t}+\frac{k}{2}\frac{{\partial }^{2}u}{\partial t^2}+\frac{h^2}{4}\frac{{\partial }^{3}u}{\partial t\partial x^2}+\frac{k^2}{6}\frac{{\partial }^{3}u}{\partial t^3}+\frac{k{h}^2}{8}\frac{{\partial }^{4}u}{\partial t^2\partial x^2}+\cdots$ (22)

${\left[\frac{{\left(u^3\right)}^{n+1}+{\left(u^3\right)}^n}{2}\right]}_x=\frac{\partial u^3}{\partial x}+\frac{k}{2}\frac{{\partial }^2{u}^3}{\partial t\partial x}+\frac{h^2}{4}\frac{{\partial }^3{u}^3}{\partial x^3}+\frac{k{h}^2}{8}\frac{{\partial }^4{u}^3}{\partial t\partial x^3}+\frac{7h^4}{120}\frac{{\partial }^5{u}^3}{\partial x^5}+\cdots$ (23)

${\left(\frac{u^{n+1}+u^n}{2}\right)}_{xxx}=\frac{{\partial }^{3}u}{\partial x^3}+\frac{k}{2}\frac{{\partial }^{4}u}{\partial t\partial x^3}+\frac{k{h}^2}{8}\frac{{\partial }^{6}u}{\partial t\partial x^5}+\frac{h^2}{4}\frac{{\partial }^{5}u}{\partial x^5}+\cdots$ (24)

By using these expressions into (21), and collecting similar terms, we will get the local truncation error (LTE)

$\begin{array}{c}LTE=\left[A\left(E\right)\frac{\partial u}{\partial t}+\frac{\epsilon }{3}B\left(E\right)\frac{\partial u^3}{\partial x}+C\left(E\right)\frac{{\partial }^{3}u}{\partial x^3}\right]\\ +\frac{k}{2}\frac{\partial }{\partial t}\left[A\left(E\right)\frac{\partial u}{\partial t}+\frac{\epsilon }{3}B\left(E\right)\frac{\partial u^3}{\partial x}+C\left(E\right)\frac{{\partial }^{3}u}{\partial x^3}\right]\\ +\frac{k{h}^2}{8}\frac{{\partial }^{3}u}{\partial t\partial x^2}\left[A\left(E\right)\frac{\partial u}{\partial t}+\frac{\epsilon }{3}B\left(E\right)\frac{\partial u^3}{\partial x}+C\left(E\right)\frac{{\partial }^{3}u}{\partial x^3}\right]\\ +\frac{h^2}{4}\frac{{\partial }^2}{\partial x^2}\left[A\left(E\right)\frac{\partial u}{\partial t}+\frac{\epsilon }{3}B\left(E\right)\frac{\partial u^3}{\partial x}+C\left(E\right)\frac{{\partial }^{3}u}{\partial x^3}\right]\\ +\left[\frac{k^2}{6}\frac{{\partial }^{3}u}{\partial t^3}+\frac{7h^4}{120}\frac{{\partial }^5{u}^3}{\partial x^5}+\cdots \right]\end{array}$ (25)

The first four brackets (25) are zero by the MKdV equation, and the LTE will be reduced into

$LTE=\left[\frac{k^2}{6}\frac{{\partial }^{3}u}{\partial t^3}+\frac{7h^4}{120}\frac{{\partial }^5{u}^3}{\partial x^5}+\cdots \right]$ (26)

So, Scheme 1 is of second order accuracy in time and fourth order in space; $O\left(k^2+h^4\right)$. 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

$LTE=\left[\frac{k^2}{6}\frac{{\partial }^{3}u}{\partial t^3}+\frac{7h^4}{120}\frac{{\partial }^5{u}^3}{\partial x^5}+\cdots \right]\to 0\text{as}h\text{and}k\to 0$

5. Stability of the Scheme

In this section, we want to study the stability of the proposed schemes [16]. Our stability analysis is based on the Von Neumann theory in which the growth factor of atypical Fourier mode defined as

$U_i^n={\text{e}}^{\alpha nk}{\text{e}}^{i\beta mh}$ (27)

where $i=\sqrt{-1}$, $\beta$ is real number and $\alpha =\alpha \left(\beta \right)$ in general complex. To implement the Fourier stability analysis, the MKdV equation needs to be linearized by assuming the nonlinear term $u^2$ in $u^2{u}_x$ is locally constant. Using this, we will get the linear version of (6) and this can be given as

$A\left(E\right)\left[U_m^{n+1}-U_m^n\right]+p_1\omega B\left(E\right)\left[U_m^{n+1}+U_m^n\right]+p_{2}C\left(E\right)\left[U_m^{n+1}+U_m^n\right]=0.$ (28)

where

$p_1=\frac{k\epsilon }{6},p_2=\frac{k\mu }{2},\omega =\max \left|{\left(U^2\right)}_m^n\right|$ (29)

By substituting (27) into (28), we get after some manipulation the following equation

$\left[\tilde{A}+i{p}_1\omega \tilde{B}+i{p}_2\tilde{C}\right]{\text{e}}^{\alpha k}=\left[\tilde{A}-i{p}_1\omega \tilde{B}-i{p}_2\tilde{C}\right]$ (30)

where

$\tilde{A}=\frac{1}{120}\left[2cos\left(2\theta \right)+52cos\left(\theta \right)+66\right]\mathrm{,}$

$\tilde{B}=\frac{1}{24h}\left[2isin\left(2\theta \right)+20isin\left(\theta \right)\right]\mathrm{,}$ (31)

$\tilde{C}=\frac{1}{2h^3}\left[2isin\left(2\theta \right)-4isin\left(\theta \right)\right]\mathrm{,}$

By solving Equation (30) for ${\text{e}}^{\alpha k}$, we can get

${\text{e}}^{\alpha k}=\frac{\left[\tilde{A}-i\left(p_1\omega \tilde{B}+p_2\tilde{C}\right)\right]}{\left[\tilde{A}+i\left(p_1\omega \tilde{B}+p_2\tilde{C}\right)\right]}$ (32)

It is easy to see that, $\left|{\text{e}}^{\alpha k}\right|\le 1$, 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

$u_t+\frac{\epsilon }{3}{\left(u^3\right)}_x+\mu u_{xxx}=0,\left(x,t\right)\in D=\left[x_l,x_r\right]\times \left(0\le T\right]$ (33)

with the homogenous boundary $u\to 0$ as $x\to \pm \infty$ and the initial condition

$u\left(x\mathrm{,0}\right)=u_0\left(x\right)\mathrm{,}x\in \left[x_l\mathrm{,}{x}_r\right]$ (34)

The exact solution of the MKdV (33) is

$u\left(x,t\right)=A\text{sech}\left(\beta \left(x-\lambda t-x_0\right)\right),$ (35)

where A is amplitude and $\beta$ is the width of the single solitary wave, they are defined as

$A=\sqrt{\frac{6\lambda }{\epsilon }},\beta =\sqrt{\frac{\lambda }{\mu }}$ (36)

where $\epsilon \mathrm{,}\mu \mathrm{,}\lambda$ and $x_0$ are arbitrary constants. The MKdV equation satisfies the conserved quantities:

$I_1={\displaystyle {\int }_{-\infty }^{\infty }}u\left(x\mathrm{,}t\right)\text{d}x$ (37)

$I_2={\displaystyle {\int }_{-\infty }^{\infty }}{u}^2\left(x\mathrm{,}t\right)\text{d}x$ (38)

$I_3={\displaystyle {\int }_{-\infty }^{\infty }}\left(u^4\left(x\mathrm{,}t\right)-\frac{6\mu }{\epsilon }{u}_x^2\left(x\mathrm{,}t\right)\right)\text{d}x$ (39)

The exact values of the conserved quantities (37) - (39) using the exact solution (35) are

$I_1=\pi \sqrt{\frac{6\mu }{\epsilon }},I_2=\frac{12\sqrt{\mu \lambda }}{\epsilon },I_3=\frac{24{\lambda }^2}{{\epsilon }^2}\sqrt{\frac{\mu }{\lambda }},$ (40)

To demonstrate the efficiency and accuracy of the presented methods for solving the MKdV equation, we use the $L_{\infty }$ error norm

$L_{\infty }=\Vert U_m^n-u_m^n\Vert =\underset{1\le m\le M}{max}\left|U_m^n-u_m^n\right|\mathrm{,}$ (41)

and the $L_2$ error norm

$L_2={\Vert U_m^n-u_m^n\Vert }_2\simeq \sqrt{h{\displaystyle \underset{n=1}{\overset{M}{\sum }}}{\left|U_m^n-u_m^n\right|}^2}\mathrm{.}$ (42)

Two values of $\epsilon$ are used. The first case: $\epsilon =6$ and $\mu =1$, 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 $\epsilon =3$ and $\mu =1$, and in this case, a comparison between our methods with [1].

6.1. The First Case: ( $\epsilon =6,\mu =1$ )

6.1.1. Single Soliton

To test our numerical methods, we choose the initial condition

$u\left(x,0\right)=A\text{sech}\left(\beta \left(x-x_0\right)\right),0\le x\le 80$ (43)

subject to the homogenous boundary conditions. In order to generate the numerical solutions, the following parameters are used

$\epsilon =6,\mu =1,h=0.1,k=0.01,x_0=40,\lambda =0.25,TOL={10}^{-10},$ (44)

The exact values of the conserved quantities using (40) are

$I_1=\pi ,I_2=2\sqrt{\lambda },I_3=\frac{2}{3}\sqrt{{\lambda }^3}.$ (45)

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⁻¹⁰, 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

$RTC=\frac{\ln \left[\frac{E_1}{E_2}\right]}{\ln \left[\frac{h_1}{h_2}\right]},E\left(h\right)=\underset{1\le m\le M}{\max }\left|U_m^n-u_m^n\right|.$ (46)

To calculate the rate convergence of Scheme 1 and Scheme 2 in space using (46), we choose $k=0.01$, for different values of h. we calculate the $L_{\infty }$ 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 $h=0.1$, for different values of k, we calculate the $L_{\infty }$ 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

$u\left(x,0\right)={\displaystyle \underset{j=1}{\overset{2}{\sum }}}{A}_{{}^j}\text{sech}\left({\lambda }_{{}^j}\left(x-x_j\right)\right),$ (47)

${\lambda }_{{}^j}$ and $x_j$ are arbitrary constants. To ensure an interaction of two solitary waves for Schemes 1, 2, 3 and 4, we choose the set of parameters $\epsilon =6$, $\mu =1$, $h=0.1$, $k=0.01$, ${\lambda }_1=2.0$, ${\lambda }_2=1.0$, $x_1=15$ and $x_2=25$, 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.

$0\le x\le 80$. 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, λ₁ = 2.0, λ₂ = 1.0.

Figure 5. Two solitons interaction for Scheme 3 with ε = 6, μ = 1, h = 0.1, k = 0.01, λ₁ = 2.0, λ₂ = 1.0.

Figure 6. Two solitons interaction for Schemes 2 and 4 with parameters ε = 6, μ = 1, h = 0.1, k = 0.01, λ₁ = 2.0, λ₂ = 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

$u\left(x,0\right)={\displaystyle \underset{j=1}{\overset{3}{\sum }}}{A}_{{}^j}\text{sech}\left({\lambda }_{{}^j}\left(x-x_j\right)\right),$ (48)

subject to the homogenous boundary conditions. We have considered the Schemes 1, 2, 3 and 4 with the set of parameters $\epsilon =6$, $\mu =1$, $h=0.1$, $k=0.01$, ${\lambda }_1=2$, ${\lambda }_2=1$, ${\lambda }_3=0.5$, $x_1=10$, $x_2=25$ and $x_3=35$ over the interval $0\le x\le 80$. 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, λ₁ = 2.0, λ₂ = 1.0, λ₃ = 0.5.

Figure 8. Three solitons interaction for Scheme 3 with ε = 6, μ = 1, h = 0.1, k = 0.01, λ₁ = 2.0, λ₂ = 1.0, λ₃ = 0.5.

Figure 9. Three solitons interaction for Schemes (2, 4) with ε = 6, μ = 1, h = 0.1, k = 0.01, λ₁ = 2.0, λ₂ = 1.0, λ₃ = 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 ( $\epsilon =3$ and $\mu =1$ )

Single Soliton

In order to compare our results with [1], we choose the initial condition

$u\left(x,0\right)=A\text{sech}\left(\beta \left(x-x_0\right)\right),0\le x\le 80,$ (49)

The following parameters are used

$\epsilon =3,\mu =1,h=0.1,k=0.01,x_0=20,\lambda =0.845,TOL={10}^{-10}$ (50)

The exact values of the conserved quantities in this case are

$I_1=\pi \sqrt{2},I_2=4\sqrt{\lambda },I_3=\frac{8}{3}\sqrt{{\lambda }^3}$ (51)

In Table 17, we display the errors for $t=1$, 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.

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 [22]

$u\left(x,t\right)=-4q\text{sech}\left(R\right)\left[\frac{\cos \left(f\right)-z\sin \left(f\right)\tanh \left(R\right)}{1+z^2{\sin }^2\left(f\right){\text{sech}}^2\left(R\right)}\right],$ (52)

where

$\begin{array}{l}R=2qx+8q\left(3p^2-q^2\right)t,\\ f=2px+8p\left(p^2-3q^2\right)t,z=\frac{p}{q},\end{array}$ (53)

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

$u\left(x,0\right)=-4q\text{sech}\left(2qx\right)\left[\frac{\cos \left(2px\right)-z\sin \left(2px\right)\tanh \left(2qx\right)}{1+z^2{\sin }^2\left(2px\right){\text{sech}}^2\left(2qx\right)}\right],$ (54)

with the set of parameters

$h=0.0,k=0.0001,p=0.2,q=1,t=0,2,4,\cdots ,10$ (55)

The exact values of the conserved quantities I₁ and I₂ are 0.0 and 8q, respectively. The conserved quantities are given in Table 18.

Which is exactly conserved, and this is a credit for our proposed schemes. The simulation of the breather is given in Figure 13 and Figure 14.

7. Birth of Solitons

7.1. Birth of Solitons Test 1

To study the birth of solitons, we choose the initial condition

$u\left(x\mathrm{,0}\right)=exp\left(-0.01x^2\right)$ (56)

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).

$h=0.1,k=0.001,t=0,0.4,0.8,\cdots ,20$ (57)

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

$u\left(x,0\right)=\text{sech}\left(\omega x\right),$ (58)

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

$h=0.1,k=0.001,\omega =0.5,0.25,t=0,0.4,\cdots ,20$ (59)

The results of the conserved quantities for $w=0.25$ and $w=0.5$ are given in Table 20 and Table 21. The simulation of this test is given in Figure 16 and Figure 17.

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 $w=0.5$ and four solitons in case of $w=0.25$. The conserved quantities are almost conserved. A very interesting thing we can anticipate the number of created solitons by using the formula

$Ns\left(\omega \right)=\frac{1}{\omega }\mathrm{.}$ (60)

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.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Karakoc, S.B. (2018) Numerical Solutions of the Modified KdV Equation with Collocation Method. Malaya Journal of Mathematik, 6, 835-842.
[2]	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. https://doi.org/10.1016/j.cnsns.2003.12.009
[3]	Biswas, A. and Raslan, K.R. (2011) Numerical Simulation of the Modified Korteweg-de Vries Equation. Physics of Wave Phenomena, 19, 142-147. https://doi.org/10.3103/S1541308X11020105
[4]	Raslan, K.R. and Baghdady, H.A. (2015) A Finite Difference Scheme for the Modified Korteweg-de Vries Equation. General Mathematics Notes, 27, 101-113.
[5]	Raslan, K.R. and Baghdady, H.A. (2014) New Algorithm for Solving the Modified Korteweg-de Vries (mKdV) Equation. International Journal of Research and Reviews in Applies Sciences, 18, 59-64.
[6]	Wazwaz, A.-M. (2014) A Variety of (3+1)-Dimensional mKdV Equations Derived by Using the mKdV Recursion Operator. Computers and Fluids, 93, 41-45. https://doi.org/10.1016/j.compfluid.2014.01.010
[7]	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. https://doi.org/10.1016/j.chaos.2015.03.018
[8]	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. https://doi.org/10.1007/s40995-017-0238-5
[9]	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. https://doi.org/10.24200/sci.2017.4096
[10]	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. https://doi.org/10.1103/PhysRevE.81.056708
[11]	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. https://doi.org/10.1103/PhysRevE.64.026608
[12]	Wazwaz, A.M. (2009) Partial Differential Equations and Solitary Waves Theory. Nonlinear Physical Science, Springer. https://doi.org/10.1007/978-3-642-00251-9
[13]	Ismail, M.S. (2008) Numerical Solution of Complex Modified Korteweg-de Vries Equation by Petrov-Galerkin Method. Applied Mathematics and Computation, 202, 520-531. https://doi.org/10.1016/j.amc.2008.02.033
[14]	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. https://doi.org/10.1002/num.20343
[15]	Ismail, M.S. and Taha, T.R. (1998) A Numerical Study of Compactons. Mathematics and Computers in Simulation, 47, 519-530. https://doi.org/10.1016/S0378-4754(98)00132-3
[16]	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. https://doi.org/10.1016/j.cnsns.2007.12.005
[17]	Ismail, M.S. and Ashi, H.A. (2014) A Numerical for Hirota-Satsuma Coupled KdV Equation. Abstract and Applied Analysis, 2014, Article ID: 819367. https://doi.org/10.1155/2014/819367
[18]	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. https://doi.org/10.1080/00207160600747938
[19]	Ismail, M.S. and Alaseri, S.H. (2016) Computational Methods for Three Coupled Nonlinear Schrodinger Equations. Applied Mathematics, 7, 2110-2131. https://doi.org/10.4236/am.2016.717168
[20]	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. https://doi.org/10.1080/00207160410001661339
[21]	Ismail, M.S. and Taha, T.R. (2001) Numerical Simulation of Coupled Nonlinear Schrodinger Equation. Mathematics and Computer Simulations, 56, 547-562. https://doi.org/10.1016/S0378-4754(01)00324-X
[22]	Didenkulova, E. and Pelinovsky, E. (2020) Breather’s Properties within the Framework of the Modified Koreweg-de Vries Equation. Symmetry, 12, 638. https://doi.org/10.3390/sym12040638