Numerical Solution of Kortweg-de Vries Equation

Fakhirah Alotaibi, M. S. Ismail

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

DOI: 10.4236/am.2020.114025                

Abstract

In this paper, we are going to derive numerical methods for solving the KdV equation using Pade approximation for space direction, trapezoidal and implicit mid-point rule in the time direction. The schemes will be analyzed for accuracy and stability. The exact solution and the conserved quantities will be used to display the efficiency and the robustness of the proposed schemes. Interaction of two and three solitons will be conducted. The numerical results showed, interaction behavior is elastic and the conserved quantities are conserved which is a good indication of the reliability of the schemes under consideration.

Keywords

KdV, Soliton, Pade Approximation, Numerical Schemes, Interaction of Solitons

1. Introduction

In this work, we consider the well known Kortweg-de Vries (KdV) equation [1]

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

with 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 homogenous boundary conditions

$u\left(x_l,t\right)=u\left(x_l,t\right)=0,u_x\left(x_l,t\right)=u_x\left(x_r,t\right)=0$ (3)

where x and t denote the spatial and temporal variables, respectively, and $u\left(x\mathrm{,}t\right)$ is the unknown function. The KdV equation is one of the mostly analyzed nonlinear partial differential equation. It is used to model the propagation of water waves and also arises in other areas such as hydrodynamics waves in cold plasma, and acoustic waves in harmonic crystals. The importance of this equation lies not only in their applications but also in some special properties which are not expected from a nonlinear partial differential equation. In fact, the KdV equation has certain special solution known as solitary wave solution [2] [3] which retains their shape in time and moves to the right. It can also have multi-soliton solutions, where the individual solutions interact nonlinearly at close distance and then move away from each other without changing in their shapes.

Once these interesting aspects had been discovered, several authors investigated the solution of the KdV equation from both the theoretical and numerical point view. Among the methods available in the literature, Hirota and Adomian methods, these methods are able to provide explicit solution which is useful to check the accuracy of the numerical methods aiming at approximating the solution of the KdV equation. Many researchers have derived efficient numerical methods for investigating the numerical solution of the KdV equation. For instance, in [4] a class of three level finite difference scheme for solving KdV equation is proposed firstly by Zabusky and Kruskal who discovered the concept of soliton, which can be defined as a localized waves with special interaction properties, while studying the results of a numerical computation on the KdV equation. The proposed three-level scheme which requires a starting procedure satisfies momentum and energy conservation laws. In [5], a two-level hopscotch algorithm is derived and the stability and dispersion is analyzed. The Crank-Nicolson extrapolation scheme for the KdV equation is proposed in [6]. A hybrid numerical method for KdV equation is proposed by a combination of finite difference and sinc collocation method is proposed in [7]. In [8], a quadratic B-spline Galerkin finite element method is proposed. The Galerkin method for KdV equation using a new basis of smooth piecewise cubic polynomials is developed for simulating the motions and interaction of solitary waves presented in [9]. Petrov Galerkin method [10] [11] [12] [13] [14] is used to solve KdV and KdV like equations to obtain highly accurate results. Numerical solution of KdV equation by Galerkin B-spline finite element method is presented in [15]. In [16], The conservation and convergence of two finite difference schemes for KdV equations with initial and boundary conditions derived, the first scheme is a two-level nonlinear implicit finite difference scheme which is proved to be unconditional convergent, and the second one is a three-level linearized finite difference scheme and proved to be conditionally convergent.

In this paper, we will derive numerical schemes for the KdV equation; based on the Padé approximation, fourth approximation of the space derivatives is used. The second order implicit schemes using, trapezoidal rule, implicit midpoint rule (both are implicit and unconditionally stable), and the explicit Runge-Kutta method of order four are used to approximate the time derivative. A nonlinear penta-diagonal system is obtained. Newton’s method is used to solve this system in case of the implicit schemes. The implicit schemes are unconditionally stable according to the von-Nueumann stability analysis. Several numerical examples are carefully designed to validate the efficiency and robustness of the proposed schemes, single soliton, interaction of two and three solitons.

The rest of the paper is organized as follow. In Section 2, we present three different schemes using fourth order Pade approximation in space direction and second order in time direction using trapezoidal and implicit mid- point rule. Also we implement the explicit fourth Runge-Kutta method to solve the first order differential system in time. In Section 3, the accuracy of the proposed schemes is given. In Section 4, the stability analysis of scheme 1 and scheme 2 are derived using von-Neumann stability analysis. In Section 5, we present various numerical tests which validate the accuracy and the efficiency of the proposed schemes. Finally, we conclude with a brief discussion in Section 6.

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 sizes in the spatial and temporal directions respectively, 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 Pade approximation, the following approximations for space derivatives are used [17] [18] [19]

$u_x\approx \frac{B\left(E\right)}{A\left(E\right)}+O\left(h^4\right)\mathrm{,}{u}_{xxx}\approx \frac{C\left(E\right)}{A\left(E\right)}+O\left(h^4\right)$ (4)

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

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

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

and the shift operator E defined by

$E^j{u}_m=u_{m+j}$

Now by using these approximations in KdV Equations (1)-(3), we obtain the first order differential system in time

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

By using (4), the system in (5) can be written as

$\begin{array}{l}\frac{1}{120}\left({\stackrel{\dot{}}{U}}_{m+2}+26{\stackrel{\dot{}}{U}}_{m+1}+66{\stackrel{\dot{}}{U}}_m+26{\stackrel{\dot{}}{U}}_{m-1}+{\stackrel{\dot{}}{U}}_{m-2}\right)\\ +\frac{1}{8h}\left(U_{m+2}^2+10U_{m+1}^2-10U_{m-1}^2-U_{m-2}^2\right)\\ +\frac{1}{2h^3}\left(U_{m+2}-2U_{m+1}+2U_{m-1}-U_{m-2}\right)=0,m=1,2,\cdots ,M-1\end{array}$ (6)

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

$U_{-1}=U_0=0,U_N=U_{N+1}=0$

The system in (6) can be written in a matrix vector form as

$M\stackrel{\dot{}}{U}=F\left(U\right)\mathrm{,}$ (7)

where M is the penta diagonal matrix

$M=\frac{1}{120}\left(\begin{array}{cccccccc}66& 26& 1& 0& \cdots & \cdots & \cdots & 0\\ 26& 66& 26& 1& 0& \cdots & \cdots & 0\\ 1& 26& 66& 26& 1& 0& \cdots & 0\\ 0& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & 1\\ ⋮& \ddots & \ddots & \ddots & 1& 26& 66& 26\\ 0& \cdots & \cdots & \cdots & 0& 1& 26& 66\end{array}\right)\mathrm{,}$

and the elements of $F\left(U\right)$ are

$\begin{array}{l}{F}_m\left(U\right)=-\frac{1}{8h}\left(U_{m+2}^2+10U_{m+1}^2-10U_{m-1}^2-U_{m-2}^2\right)\\ -\frac{1}{2h^3}\left(U_{m+2}-2U_{m+1}+2U_{m-1}-U_{m-2}\right)\mathrm{,}m=\mathrm{1,2,}\cdots \mathrm{,}N-1\end{array}$ (8)

$U={\left[U\mathrm{,}U\mathrm{,}\cdots \mathrm{,}{U}_{N-1}\right]}^t\mathrm{.}$

The previous differential system can be solved by

1) Trapezoidal rule

$M\left[\frac{U^{n+1}-U^n}{k}\right]=\frac{1}{2}\left[F\left(U^{n+1}\right)+F\left(U^n\right)\right]\mathrm{,}$ (9)

Trapezoidal method is of second order, implicit scheme, and A-stable method

2) Implicit mid-point rule

$M\left[\frac{U^{n+1}-U^n}{k}\right]=F\left[\frac{U^{n+1}+U^n}{2}\right]$ (10)

Implicit middle point rule is of second order, implicit scheme, and A-stable method.

The previous methods (Crank-Nicolson like schemes), are usually produce unconditionally stable method, nonlinear penta diagonal system to be solved at each time step. Explicit Runge-Kutta method of fourth order will be used as well, solution of four linear penta diagonal systems are needed for each time step using Crouts method, and the method is conditionally stable. The details of the proposed schemes will be discussed in the following subsections.

2.1. Scheme 1 (Trapezoidal Rule)

Trapezoidal rule method (9) which is of second order accuracy and A-stable is used to solve the ordinary differential system in (7), we assume $U_m^n$ to be the fully discrete approximation of the exact solution $u_m^n$ this will lead to the nonlinear implicit scheme

$\begin{array}{l}\left(U_{m+2}^{n+1}+26U_{m+1}^{n+1}+66U_m^{n+1}+26U_{m-1}^{n+1}+U_{m-2}^{n+1}\right)\\ -\left(U_{m+2}^n+26U_{m+1}^n+66U_m^n+26U_{m-1}^n+U_{m-2}^n\right)\\ +p_1[{\left(U^2\right)}_{m+2}^{n+1}+10{\left(U^2\right)}_{m+1}^{n+1}-10{\left(U^2\right)}_{m-1}^{n+1}-{\left(U^2\right)}_{m-2}^{n+1}\\ +{\left(U^2\right)}_{m+2}^n+10{\left(U^2\right)}_{m+1}^n-10{\left(U^2\right)}_{m-1}^n-{\left(U^2\right)}_{m-2}^n]\\ +p_2\left(U_{m+2}^{n+1}-2U_{m+1}^{n+1}+2U_{m-1}^{n+1}-U_{m-2}^{n+1}+U_{m+2}^n-2U_{m+1}^n+2U_{m-1}^n-U_{m-2}^n\right)=0\end{array}$ (11)

where

$p_1=\frac{15k}{2h},p_2=\frac{30k}{h^3},m=1,2,\cdots ,M,U_{-1}=U_0=0,U_N=U_{N+1}=0$ (12)

The system in (11) is a nonlinear penta-diagonal system in the unknown vector $U^{n+1}$. The solution of this system can be obtained by many methods, like: Newton’s method; the Predictor-Corrector method and linearization techniques. Newton’s method will be adopted in this work.

2.2. Scheme 2 (Implicit Mid-Point Rule)

A second finite difference scheme obtained by using (10) the mid-point rule which is of second order and A-stable, using this will lead us to the nonlinear scheme of the form

The system in (10) can be given by the finite difference equation.

$\begin{array}{l}\left(U_{m+2}^{n+1}+26U_{m+1}^{n+1}+66U_m^{n+1}+26U_{m-1}^{n+1}+U_{m-2}^{n+1}\right)\\ -\left(U_{m+2}^n+26U_{m+1}^n+66U_m^n+26U_{m-1}^n+U_{m-2}^n\right)\\ +p_1\left(\left[V_{m+2}^2+10V_{m+1}^2-10V_{m-1}^2-V_{m-2}^2\right]\right)\\ +p_2(\left(U_{m+2}^{n+1}-2U_{m+1}^{n+1}+2U_{m-1}^{n+1}-U_{m-2}^{n+1}\right)\\ +\left(U_{m+2}^n-2U_{m+1}^n+2U_{m-1}^n-U_{m-2}^n\right))=0,m=1,2,\cdots ,M-1\end{array}$ (13)

where

$V_i=\frac{U_i^{n+1}+U_i^n}{2},p_1=\frac{15k}{4h},p_2=\frac{30k}{h^3},U_{-1}=U_0=0,U_N=U_{N+1}=0$ (14)

The system in (13) is a nonlinear penta-diagonal system which can be solved by many iterative methods, Newton’s method is used to do this job in this work.

2.3. Scheme 3 (Runge-Kutta of Order 4)

We can also solve the first order differential system (7) by using the explicit Runge -Kutta method of fourth order

$U^{n+1}=U^n+\frac{1}{6}\left[K_1+2K_2+2K_3+K_4\right]\mathrm{,}$ (15)

$K_1=k{M}^{-1}F(Un)$

$K_2=k{M}^{-1}F\left(U^n+\frac{1}{2}{K}_1\right)$

$K_3=k{M}^{-1}F\left(U^n+\frac{1}{2}{K}_2\right)$

$K_4=k{M}^{-1}F\left(U^n+K_3\right)$

where

$\begin{array}{l}{F}_m\left(U\right)=-\frac{1}{8h}\left(U_{m+2}^2+10U_{m+1}^2-10U_{m-1}^2-U_{m-2}^2\right)\\ -\frac{1}{2h^3}\left(U_{m+2}-2U_{m+1}+2U_{m-1}-U_{m-2}\right)\mathrm{,}m=\mathrm{1,2,}\cdots \mathrm{,}N-1\end{array}$ (16)

The resulting system is of fourth order accuracy in both directions time and space, it is conditionally stable. In each time step, we need to solve four linear penta diagonal systems to find $K_1\mathrm{,}{K}_2\mathrm{,}{K}_3$ and $K_4$. The solution of the linear penta-diagonal system can be obtained by using Crout’s method, only one LU factorization is needed at the beginning of the calculations which can be done easily.

3. Accuracy of the Proposed Schemes

To study the accuracy of Scheme 1, we replace the numerical solution $U_i^n$ by the the exact solution $u_i^n$ in (8) to get

$\begin{array}{l}\left(u_{m+2}^{n+1}+26u_{m+1}^{n+1}+66u_m^{n+1}+26u_{m-1}^{n+1}+u_{m-2}^{n+1}\right)\\ -\left(u_{m+2}^n+26u_{m+1}^n+66u_m^n+26u_{m-1}^n+u_{m-2}^n\right)\\ +p_1[{\left(u^2\right)}_{m+2}^{n+1}+10{\left(u^2\right)}_{m+1}^{n+1}-10{\left(u^2\right)}_{m-1}^{n+1}-{\left(u^2\right)}_{m-2}^{n+1}\\ +{\left(u^2\right)}_{m+2}^n+10{\left(u^2\right)}_{m+1}^n-10{\left(u^2\right)}_{m-1}^n-{\left(u^2\right)}_{m-2}^n]\\ p_2\left(\left(u_{m+2}^{n+1}-2u_{m+1}^{n+1}+2u_{m-1}^{n+1}-u_{m-2}^{n+1}\right)+\left(u_{m+2}^n-2u_{m+1}^n+2u_{m-1}^n-u_{m-2}^n\right)\right)=0\end{array}$ (17)

Taylor’s expansion for all terms in (17) about 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$ (18)

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

${\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$ (20)

By substituting these expressions into (17), and by collecting similar terms, we will get the local truncation error (LTE) as

$\begin{array}{c}LTE=\left[\frac{\partial u}{\partial t}+3\frac{\partial u^2}{\partial x}+\frac{{\partial }^{3}u}{\partial x^3}\right]+\frac{k}{2}\frac{\partial }{\partial t}\left[\frac{\partial u}{\partial t}+3\frac{\partial u^2}{\partial x}+\frac{{\partial }^{3}u}{\partial x^3}\right]\\ +\frac{k{h}^2}{8}\frac{{\partial }^{3}u}{\partial t\partial x^2}\left[\frac{\partial u}{\partial t}+3\frac{\partial u^2}{\partial x}+\frac{{\partial }^{3}u}{\partial x^3}\right]+\frac{h^2}{4}\frac{{\partial }^2}{\partial x^2}\left[\frac{\partial u}{\partial t}+3\frac{\partial u^2}{\partial x}+\frac{{\partial }^{3}u}{\partial x^3}\right]\\ +\left[\frac{k^2}{6}\frac{{\partial }^{3}u}{\partial t^3}+\frac{7h^4}{240}\frac{{\partial }^5{u}^2}{\partial x^5}+\cdots \right]\end{array}$ (21)

The first four brackets are zero by the KdV equation, LTE will be reduced to

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

So, the scheme 1 is of second order in time and fourth order in space $O\left(k^2+h^4\right)$.

Similar analysis can be done for Scheme 2. Scheme 2 is also of second order in time and fourth order in space $O\left(k^2+h^4\right)$. Regarding Scheme 3, the scheme is fourth order in both directions, space and time $O\left(k^4+h^4\right)$.

4. Stability Analysis

In this section, we want study the stability of the proposed schemes [20]. 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}$ (23)

where $i=\sqrt{-1}$, $\beta$ is real number and $\alpha =\alpha \left(\beta \right)$ in general complex. To implement the Fourier stability analysis, the KdV equation needs to be linearized by assuming that $u\left(x\mathrm{,}t\right)$ in the nonlinear term u u_x is locally constant. The linear version of (11) can be given as

$\begin{array}{l}\left(U_{m+2}^{n+1}+26U_{m+1}^{n+1}+66U_m^{n+1}+26U_{m-1}^{n+1}+U_{m-2}^{n+1}\right)\\ -\left(U_{m+2}^n+26U_{m+1}^n+66U_m^n+26U_{m-1}^n+U_{m-2}^n\right)\\ +p_1\omega \left[U_{m+2}^{n+1}+10U_{m+1}^{n+1}-10U_{m-1}^{n+1}-U_{m-2}^{n+1}+U_{m+2}^n+10U_{m+1}^n-10U_{m-1}^n-U_{m-2}^n\right]\\ +p_2\left(U_{m+2}^{n+1}-2U_{m+1}^{n+1}+2U_{m-1}^{n+1}-U_{m-2}^{n+1}+U_{m+2}^n-2U_{m+1}^n+2U_{m-1}^n-U_{m-2}^n\right)=0\end{array}$ (24)

where

$p_1=\frac{15k}{h},p_2=\frac{30k}{h^3},\omega =\max \left|U_i^n\right|$ (25)

Substituting (23) into (24), we get after some manipulation the following equation

$\left[{\gamma }_1+\left(p_1\omega {\gamma }_2+p_2{\gamma }_3\right)i\right]e^{\alpha k}=\left[{\gamma }_1-\left(p_1\omega {\gamma }_2+p_2{\gamma }_3\right)i\right]$ (26)

where

$\begin{array}{l}{\gamma }_1=66+2\cos \left(2\beta h\right)+52\cos \left(\beta h\right),\\ {\gamma }_2=2\sin \left(2\beta h\right)+20\sin \left(\beta h\right),\\ {\gamma }_3=2\sin \left(2\beta h\right)+4\sin \left(\beta h\right)\end{array}$ (27)

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

${\text{e}}^{\alpha k}=\frac{{\gamma }_1-\left(p_1\omega {\gamma }_2+p_2{\gamma }_3\right)i}{{\gamma }_1+\left(p_1\omega {\gamma }_2+p_2{\gamma }_3\right)i}$ (28)

It is easy to see from (28) that, $\left|{\text{e}}^{\alpha k}\right|=1$, thus we can say that scheme 1 and scheme 2 are unconditionally stable according von- Neumann stability analysis. For scheme 3 the numerical results showed that the scheme is conditionally stable (since it is explicit).

5. Numerical Results

In this part, several numerical examples are carefully designed to validate the efficiency of the proposed methods. In the physical opinion, the motion and interactions of solitons will be considered. In addition, some conservation laws that KdV equation satisfies will be examined by numerically calculating (using trapezoidal rule) the following three invariants corresponding to conservation of mass, momentum and energy as defined [7]

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

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

$I_3={\displaystyle {\int }_{-\infty }^{\infty }}\left[2u^3\left(x\mathrm{,}t\right)-{\left(u_x\left(x\mathrm{,}t\right)\right)}^2\right]\text{d}x=\text{Val}3$ (31)

To demonstrate the efficiency and accuracy of the presented method for the KdV equation, we use the maximum error norm given by

$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{,}$ (32)

5.1. Single Soliton

To test our numerical methods, we choose the initial condition

$u\left(x\mathrm{,0}\right)=2{\lambda }^2{\text{sech}}^2\left(\lambda x-x_0\right)\mathrm{,}0\le x\le 100$ (33)

subject to the homogenous boundary conditions. The exact solution for this test is

$u\left(x\mathrm{,}t\right)=2{\lambda }^2{\text{sech}}^2\left(\lambda \left(x-4{\lambda }^{2}t-x_0\right)\right)\mathrm{,}$ (34)

In order to generate the numerical solutions, the following parameters are used

$h=0.1,k=0.001,x_0=30,\lambda =0.5,TOL={10}^{-10}.$ (35)

The exact values of the conserved quantities (29)-(31) using (34) are

$I_1=4\lambda ,I_2=\frac{16}{3}{\lambda }^3,I_3=\frac{64}{5}{\lambda }^5$ (36)

In Table 1 and Table 2, we display the conserved quantities and the error for Scheme 1 and Scheme 2, respectively. The results show that the methods conserve the conserved quantities, with high accuracy. The simulation of the single soliton is given in Figure 1.

In Table 3, we present the results for a single solution using Scheme 3. These results are obtained by using the same previous parameters except $k=0.0001,x_0=20$. From the numerical results, we have noticed that the method is conserves the conserved quantities exactly with high accuracy. The motion of single soliton displayed in Figure 2.

Table 1. Scheme 1 (Single soliton).

Table 2. Scheme 2 (Single soliton).

Table 3. Scheme 3 (Single soliton).

Figure 1. The evolution of the numerical solution of scheme 1 with $h=0.1,k=0.001,\lambda =0.5$.

Figure 2. The evolution of the numerical solution of scheme 3 with $h=0.1,k=0.0001,\lambda =0.5$.

5.2. Rate of Convergence

To calculate the order of the proposed numerical schemes. We define the rate of convergence (RTC) as

$RTC=\frac{\ln \left[\frac{E\left(h_1\right)}{E\left(h_2\right)}\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|,$ (37)

To calculate the rate convergence in space using (37), we choose $k=0.001$, with different values of h. we calculate the $L_{\infty }$ error norm and the (RTC) for different values of h, and we displayed these results in Table 4, and fourth order convergence in space is observed.

To calculate the rate convergence in time, we choose $h=0.05$, with different values of k. we calculate the $L_{\infty }$ error norm and the (RTC for different values of k, and we displayed these results in Table 5, and second order convergence in time is observed.

Table 4. Rate of convergence in space with k = 0.001, λ = 0.5.

Table 5. Rate of convergence in time with h = 0.05, λ = 0.5.

5.3. Two Solitons Interaction

As a second problem, we have discussed the behavior of the interaction of two solitary waves for Scheme (1) and (2) for different analytical solutions.

5.3.1. Sum of Two Single Solitons (a)

We choose initial condition a sum of two well separated single solitons:

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

${\lambda }_{{}^j}$ and $x_j$ are arbitrary constants. To ensure an interaction of two solitary waves, we select the set of parameters $h=0.1$, $k=0.001$, $x_1=10$, $x_2=40$, ${\lambda }_1=1$, ${\lambda }_2=0.5$, $0\le x\le 100$. From the numerical result 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 Table 6 and Table 7. For the interaction scenario, see Figure 3.

5.3.2. Exact Two Solitons Solution (b)

In this test we pick our initial condition from the exact two soliton solution [2] [3] [21]

$u\left(x\mathrm{,}t\right)=2{\left(logf\left(x\mathrm{,}t\right)\right)}_{xx}$

$f\left(x,t\right)=1+{\text{e}}^{{\eta }_1}+{\text{e}}^{{\eta }_2}+{\left(\frac{{\alpha }_1-{\alpha }_2}{{\alpha }_1+{\alpha }_2}\right)}^2{\text{e}}^{{\eta }_1+{\eta }_2},{\eta }_i={\alpha }_{i}x-{\alpha }_i^{3}t,i=1,2.$

${\eta }_i={\alpha }_{i}x-{\alpha }_i^{3}t,i=1,2.$ (39)

together with the set of parameters $h=0.1$, $k=0.001$, ${\alpha }_1=0.5$, ${\alpha }_2=0.8$, $-30\le x\le 30$. The error and conserved quantities presented in Table 8 and Table 9. The results are highly accurate and the values of conserved quantities are in good agreement with their analytical values. The simulation of this interaction is given in Figure 4 for $-10\le t\le 20$. We have noticed that, the two waves collide each other and leave the interaction region without any disturbance in their identities.

Table 6. Scheme 1 (Interaction of two solitons (a)).

Table 7. Scheme 2 (Interaction of two solitons (a)).

Table 8. Scheme 1 (Interaction of two solitons (b)).

Table 9. Scheme 2 (Interaction of two solitons (b)).

Figure 3. Two solitons interaction with $h=0.1,k=0.001,{\lambda }_1=1.0,{\lambda }_2=0.5$.

Figure 4. Numerical solution of two solitary waves for scheme 2 with $h=0.1,k=0.001,{\alpha }_1=0.5,{\alpha }_2=0.8$.

5.3.3. Two Solitons Solution (c)

We consider also the interaction of two solitary waves using the initial condition [3]

$u\left(x\mathrm{,0}\right)=12\frac{3+4cosh\left(2x\right)+cosh\left(4x\right)}{{\left[3cosh\left(x\right)+cosh\left(3x\right)\right]}^2}$ (40)

where the exact solution

$u\left(x,t\right)=12\frac{3+4\cosh \left(2x-8t\right)+\cosh \left(4x-64t\right)}{{\left[3\cosh \left(x-28t\right)+\cosh \left(3x-36t\right)\right]}^2}$ (41)

By using (40) and the set of parameters $h=0.05$, $k=0.0001$, $-0.1\le t\le 0.1$, $-10\le x\le 10$. The errors and conserved quantities are displayed in Table 10 and Table 11. The results are highly accurate and the conserved quantities are almost constants. The interaction scenarios are displayed in Figure 5. We have noticed that the two waves interact and emerge after the interaction without any disturbance in their identities.

Using Scheme 3 with the initial condition (38) and the set of parameters

$h=0.1,k=0.0001,x_1=10,x_2=30,{\lambda }_1=0.8,{\lambda }_2=0.4,0\le x\le 70$

To check the interaction of two solitons, in Table 12, we display the conserved quantities which exactly conserved. The interaction scenario displayed in Figure 6.

Table 10. Scheme 1 (Interaction of two solitons (c)).

Table 11. Scheme 2 (Interaction of two solitons (c)).

Table 12. Scheme 3 (Runge-Kutta of order 4).

Figure 5. Numerical solution of two solitary waves for scheme 1 with $h=0.05,k=0.0001$.

Figure 6. Numerical solution of two solitary waves for scheme 3 with $h=0.1,k=0.0001,{\lambda }_1=0.8,{\lambda }_2=0.4$.

5.4. Three Solitons Interaction

In this test, we want to study the interaction of three solitons having different amplitudes and traveling in the same direction. We choose the initial condition as sum of three well separated solitons of the form

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

subject to the homogenous boundary conditions. We have considered the Scheme 1 and Scheme 2, with parameters $h=0.1$, $k=0.001$, ${\lambda }_1=2$, ${\lambda }_2=1$, ${\lambda }_3=0.5$, $x_1=10$, $x_2=40$ and $x_3=50$ over the interval $0\le x\le 100$, while we have considered Scheme 3 with the same previous parameters except $k=0.0001$. In Tables 13-15, we presented the conserved quantities for the proposed schemes. All methods showed the conservation of the conserved quantities. In Figure 7 and Figure 8, we display the interaction scenario of the three solitons. We have noticed that the three solitons recover their shapes after the interaction

5.5. Single Soliton Using Periodic Boundary Conditions

The final test in this work is to study the behavior of single soliton with periodic boundary conditions (for long time simulation) using Scheme 1. We select the set of parameters

$h=0.1,k=0.001,x_0=20,\lambda =0.5,tol={10}^{-10}$ (43)

In Table 16, we presented the error and conserved quantities. We have noticed a high accuracy and exact conservation of the conserved quantities. In Figure 9, we display the motion of the single soliton for $t=\mathrm{0,}\cdots ,50$.

6. Conclusion

The numerical schemes for the KdV equation presented using Pade approximation of fourth order accuracy in space direction and second order accuracy in time using trapezoidal and implicit middle point rule. The resulting system is a nonlinear penta diagonal system, where Newton’s method is used to find the numerical solutions. Scheme 1 and Scheme 2 are unconditionally stable, while the scheme 3 is conditionally stable with fourth order accuracy in both directions.

Table 13. Scheme 1 (Three solitons).

Table 14. Scheme 2 (Three solitons).

Table 15. Scheme 3 (Three solitons).

Table 16. Scheme 1 (Periodic boundary conditions).

Figure 7. Three solitons interaction for scheme (1) and (2) with $h=0.1,k=0.001,{\lambda }_1=2.0,{\lambda }_2=1.0,{\lambda }_3=0.5$.

Figure 8. Three solitons interaction for scheme 3 with $h=0.1,k=0.0001,{\lambda }_1=2.0,{\lambda }_2=1.0,{\lambda }_3=0.5$.

Figure 9. The evolution of the numerical solution of scheme 1 with $h=0.1,k=0.001,\lambda =0.5$.

Different tests are used to check the accuracy and efficiency of the proposed schemes. The exact solutions and the conserved quantities are used to check the accuracy and the efficiency of the proposed schemes. The interaction of three solitons is presented. For long time simulation, periodic boundary conditions equipped with scheme 1 are presented. The numerical results keep the agreement with the physical phenomena with the proposed methods. Based on our numerical tests, we can say that the present methods give accurate numerical solutions and are useful for simulation of the conservation properties, motion of solitons, and the interactions of solitons for the KdV equation.

Conflicts of Interest

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

References

[1]	Korteweg, D.J. and de Vried, G. (1895) On the Change of Form of Long Waves Advancing in a Rectangular Canal and on a New Type of Long Stationary Wave. Philosophical Magazine, 39, 422-443. https://doi.org/10.1080/14786449508620739
[2]	Wazwaz, A.-M. (2009) Partial Differential Equations and Solitary Wave Theory. Springer-Verlag, Berlin Heidelberg. https://doi.org/10.1007/978-3-642-00251-9
[3]	Drazian, P.G. and Jhonson, R.S. (1988) Solitons: An Introduction. Cambridge Texts in Applied Mathematics.
[4]	Zabusky, N.J. and Kruskal, M.D. (1965) Interaction of Solitons in a Collisionless Plasma and the Recurrence of Initial States. Physical Review Letters, 15, 240-243. https://doi.org/10.1103/PhysRevLett.15.240
[5]	Greig, I.S. and Morris, J.I. (1976) A Hopscotch Method for Korteweg-de Vries. Journal of Computational Physics, 20, 64-80. https://doi.org/10.1016/0021-9991(76)90102-9
[6]	Wang, P.F. and Huang, P.Z. (2019) Convergence of the Crank-Nicolson Extrapolation Scheme for the Korteweg-de Vries Equation. Applied Numerical Mathematics, 141, 88-96. https://doi.org/10.1016/j.apnum.2019.04.001
[7]	Kong, D.S., Xu, Y.F. and Zheng, Z.S. (2019) A Hybrid Numerical Method for KdV Equation by Finite Difference and Sinc Collocation Method. Applied Mathematics and Computation, 355, 61-72. https://doi.org/10.1016/j.amc.2019.02.031
[8]	Askan, E.N. and Ozdes, A. (2006) Numerical Solution of Korteweg-de Vries Equation by Galerkin B-Spline Finite Element Method. Applied Mathematics and Computation, 175, 1256-1265. https://doi.org/10.1016/j.amc.2005.08.038
[9]	Hao, S.Y. and Xie, S.S. (2012) The Galerkin Method for the KdV Equation Using a New Basis of Smooth Piecewise Cubic Polynomials. Applied Mathematics and Computation, 218, 8659-8671. https://doi.org/10.1016/j.amc.2012.02.027
[10]	Ma, H.P. and Sun, W.W. (2001) Optimal Error Estimates of the Legender-Petrov-Galerkin Method for the Korteweg-de Vries. SIAM Journal on Numerical Analysis, 39, 1380-1394. https://doi.org/10.1137/S0036142900378327
[11]	Shen, J. (2003) A New Dual-Petrov-Galerkin Method for Third and Higher Order Odd Order Differential Equation: Application to the KdV Equation. SIAM Journal on Numerical Analysis, 41, 1595-1619. https://doi.org/10.1137/S0036142902410271
[12]	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
[13]	Ismail, M.S. (2009) Numerical Solution of Complex Modified Korteweg-de Vries Equation by Collocation Method. Communications in Nonlinear Science and Numerical Simulation, 14, 749-759. https://doi.org/10.1016/j.cnsns.2007.12.005
[14]	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
[15]	Irk, D., Idris, D. and Bulent, S. (2006) A Small Time Solution for the Korteweg-de Vries Equation Using Spline Approximation. Applied Mathematics and Computation, 173, 834-846. https://doi.org/10.1016/j.amc.2005.04.018
[16]	Shen, J., Wang, X.P. and Sun, Z.-Z. (2020) The Conservation and Convergence of Two Finite Difference Schemes for KdV Equations with Initial and Boundary Conditions. Numerical Mathematics: Theory, Methods and Applications, 13, 253-280. https://doi.org/10.4208/nmtma.OA-2019-0038
[17]	Garralon, J., Rus, F.F. and Villatoro, F.R. (2013) Numerical Interactions between Compactons and Kovatons of the Rosenau-Pikovsky K(cos) Equation. Communications in Nonlinear Science and Numerical Simulation, 18, 1576-1588. https://doi.org/10.1016/j.cnsns.2012.10.016
[18]	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, Article ID: 056708. https://doi.org/10.1103/PhysRevE.81.056708
[19]	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, Article ID: 026608. https://doi.org/10.1103/PhysRevE.64.026608
[20]	Mitchell, A.R. and Griffiths, D.F. (1980) The Finite Difference Method in Partial Differential Equations. John Wiley and Sons, Hoboken.
[21]	Taha, T.R. and Ablowttz, M.J. (1984) Analytical and Numerical Aspects of Certain Nonlinear Evolution Equations III. Numerical Korteweg-de Vries Equation. Journal of Computational Physics, 15, 231-253. https://doi.org/10.1016/0021-9991(84)90004-4