Analytical and Numerical Computations of Multi-Solitons in the Korteweg-de Vries (KdV) Equation

Hycienth O. Orapine¹, Emem Ayankop-Andi¹, Godwin J. Ibeh²

¹Department of Mathematical Sciences, Nigerian Defence Academy, Kaduna, Nigeria.
²Department of Physics, Nigerian Defence Academy, Kaduna, Nigeria.

DOI: 10.4236/am.2020.117037                

Abstract

In this paper, an analytical and numerical computation of multi-solitons in Korteweg-de Vries (KdV) equation is presented. The KdV equation, which is classic of all model equations of nonlinear waves in the soliton phenomena, is described. In the analytical computation, the multi-solitons in KdV equation are computed symbolically using computer symbolic manipulator—Wolfram Mathematica via Hirota method because of the lengthy algebraic computation in the method. For the numerical computation, Crank-Nicolson implicit scheme is used to obtain numerical algorithm for the KdV equation. The simulations of solitons in MATLAB as well as results concerning collision or interactions between solitons are presented. Comparing the analytical and numerical solutions, it is observed that the results are identically equal with little ripples in solitons after a collision in the numerical simulations; however there is no significant effect to cause a change in their properties. This supports the existence of solitons solutions and the theoretical assertion that solitons indeed collide with one another and come out without change of properties or identities.

Keywords

Korteweg-de Vries Equation, Solitons, Hirota Method, Crank-Nicolson Method

1. Introduction

The study of nonlinear evolution models which describes a large variety of physical phenomena is found to have two fascinating manifestations of opposite nature: chaos that are the apparent randomness in the behavior of perfectly deterministic systems and special kind of solitary waves called solitons. Soliton theory gives us various significant instances of nonlinear systems behaving in a persistent, quasi-linear pattern. Solitons are therefore, a consequence of a dynamic balance between dispersion and nonlinear effects in any nonlinear evolution models. They are waves of permanent form that preserve their shape while traveling over long distances. The permanent speed and form of a soliton is however not the only special property; but it is said of its special characteristic that, solitons maintain their shape and speed after collisions with other solitons, summing up in the following two basic properties; One, propagating without change of its characteristics (shape, size, velocity etc.), Two, localized waves (stable against mutual collisions and retaining their identities:). The first is a solitary wave condition acknowledged in hydrodynamics since the nineteenth century. The second implies that the wave has the property of a particle [1] [2] [3].

[4] and [5] revealed that solitons in modeling physical phenomena arise in a wide range of areas such as shallow and deep water waves, optics communications, Bose-Einstein condensates and Biological models. They are universal in nature, and can be found in different classifications: from water waves, sound waves and charge-density waves to matter waves and electromagnetic waves [6]. Solitons can be found in a variety of materials including Plasmas, Josephson junction, Polyacetylene (Proteins and DNA) molecules among others and are used for many different purposes. They may carry electrical charges in some substances and when charged, solitons travel through certain polymer chains, which tend to curve. [7] predicted that this property may one day be used in applications such as artificial muscles and also showed that solitons can be used to transfer large amounts of information over large distances with little or no errors in the signal.

Interactions or collision between solitons is perhaps the most captivating features of soliton phenomena. A critical observation was however made in [8] that solitons preservation of identities after collisions is basically suggested by numerical simulations. But in detailed analysis of the results of such numerical simulation, some ripples can be observed after a collision and it therefore seems that the original identity is not completely recovered; it is most imperative to find exact solutions of nonlinear evolution models that admit soliton solutions for proper scrutiny of collisions. [9] in the same vein, studied in contrast, “overtaking interactions” and “head-on collisions” between two electrostatic solitons in Korteweg-de Vries (KdV) equation using an approximate method, and revealed that though numerical method offers valuable insights it still limits the weight of validity of the study of soliton collisions.

In this paper we focus on the existence of more than one soliton solution called multi-solitons, since this enables us to study solitons collision or interactions especially for their said preserved behavior. We look at the famous Korteweg-de Vries (KdV) equation which has soliton solutions, and we describe the context in which this model arises, and solve it both analytically and numerically and discuss properties of the multi-soliton solutions. Korteweg-de Vries (KdV) equation has been extensively solved for single soliton solution using analytical, semi-analytical methods and, multi-solitons for numerical methods. Our motivation here is to apply dissimilar computational techniques to produce its multi-soliton solutions for scrutiny.

2. Method: Computation of Multi-Solitons in Korteweg-de Vries Equation

One of the most famous time evolution models and perhaps the simplest nonlinear system is the third order Korteweg-de Vries (KdV) equation, hereafter abbreviated as KdV equation, given in a general form as

$w_t\pm \beta w{w}_x+\alpha w_{xxx}=0,\alpha ,\beta \in ℝ$ (2.0)

where $w=w\left(x,t\right)$ is a function of two variables which represent the amplitude of the wave at position x, and time t, and α, β are arbitrary. The Equation (2.0) is nonlinear because of the product shown in the second term and of third order for the reason that the third derivative is highest. The KdV Equation (2.0) arises in many physical situations. [10] noted that it very well may be used to portray the investigation of shallow-water waves, gas dynamics, Anharmonic nonlinear grids, Hydro-attractive and Icon-acoustic waves in cool plasma, for instance. The physical derivation of this model is given in [11].

2.1. Solution of the KdV Equation by Hirota Method

We solve the KdV Equation (2.0) in particular form on the infinite line, stage by stage using Hirota method.

$w_t\pm 6w{w}_x+w_{xxx}=0,-\infty <x<\infty ,0<t<\infty$ (2.1)

with initial condition $w\left(x,0\right)=g\left(x\right).$

The solution of the KdV Equation (2.1) is a travelling wave which has permanent form, occurs due to a balance of its dispersive term ( $w_{xxx}$ ) and its nonlinear term ( $6w{w}_x$ ). This idea of balancing dispersion and nonlinearity is typical of any nonlinear evolution equation that admits soliton solution. The constant coefficients behind each of the terms are unimportant and we express the equation this way for historical reasons, particularly, the factor 6 is just a scaling factor to make the solutions (solitons or solitary waves) easier to be described via Hirota Method [8].

Bilinearization: First the KdV Equation (2.1) is rewritten as

$w_t+3{\left(w^2\right)}_x+w_{xxx}=0$ (2.2)

The supplementary function (the Cole-Hopf transformation) given as

$w\left(x,t\right)=2{\left(\ln f\right)}_{xx}=2\left(\frac{f{f}_{xx}-f_x^2}{f^2}\right)$ (2.3)

is used to transform (2.1) into a bilinear KdV equation. Thus, taking time derivative of $w\left(x,t\right)$ in (2.3), we integrate once with respect to x and substitute in (2.2), after some algebraic manipulations, we obtain the bilinear form of the KdV equation (2.1) as,

$2\left(f{f}_{xt}-f_x{f}_t+3f_{xx}^2-4f_x{f}_{xxx}+f{f}_{xxxx}\right)=0$ (2.4)

The Hirota bilinear form: To put the bilinear KdV Equation (2.4) into Hirota bilinear form in terms of the D-derivative operator as defined as follows, we assume the Hirota differential operator $D_x$, which is a binary operator defined on ordered pairs of functions $f\left(x\right)$ and $g\left(x\right)$, as follows

$D_x^{n}f\cdot g={{\left(\frac{\partial }{\partial x_1}-\frac{\partial }{\partial x_2}\right)}^{n}f\left(x_1\right)g\left(x_2\right)|}_{x_2=x_1=x}$ (2.5)

then

$D_t^m{D}_x^n\left(f\cdot g\right)={{\left(\frac{\partial }{\partial t_1}-\frac{\partial }{\partial t_2}\right)}^m{\left(\frac{\partial }{\partial x_1}-\frac{\partial }{\partial x_2}\right)}^{n}f\left(t_1,x_1\right)g\left(t_2,x_2\right)|}_{x=x_2=x_1,t=t_1=t_2}$ (2.6)

where $m,n$ are nonnegative integers and

${\left(\frac{\partial }{\partial x_1}-\frac{\partial }{\partial x_2}\right)}^n=\underset{r=0}{\overset{n}{{\displaystyle \sum }}}{\left(-1\right)}^r\left(\begin{array}{c}n\\ r\end{array}\right)\frac{{\partial }^{n-r}}{\partial x_1}\cdot \frac{{\partial }^r}{\partial x_2}$ (2.7)

where

$\left(\begin{array}{c}n\\ r\end{array}\right)=\frac{n!}{r!\left(n-r\right)!},0\le r\le n$ (2.8)

More generally we denote some sort of combination of the Hirota D-operator as a polynomial of D-operator $P\left(D\right)$ i.e.

$\begin{array}{c}P\left(D\right)\left\{f\cdot g\right\}=D_x^n\left\{f\cdot g\right\}=\underset{r=0}{\overset{n}{{\displaystyle \sum }}}{\left(-1\right)}^r\left(\begin{array}{c}n\\ r\end{array}\right)f_{\left(n-r\right)x}\cdot g_{rx}\\ =f_{nx}g-n{f}_{\left(n-1\right)x}{g}_x+\frac{n\left(n-1\right)}{2}{f}_{\left(n-2\right)x}{g}_{2x}-\cdots +{\left(-1\right)}^{n}f{g}_{nx}\end{array}$ (2.9)

where the subscripts of the functions f and g define the order of the partial derivatives with respect to x. Thus, D operates on a product of two functions like the Leibniz rule, except for the crucial sign difference.

Thus, from Equation (2.6) we have:

$D_x{D}_t\left(f\cdot g\right)=\left[f_{xt}g-f_t{g}_x\right]-\left[f_x{g}_t-f{g}_{xt}\right]$ (2.10)

Now replacing g with f to have the same function in (2.4) we have that:

$D_x{D}_t\left(f\cdot f\right)=2\left(f{f}_{xt}-f_x{f}_t\right)$ (2.11)

Similarly, from (2.9) we have that

$D_x^4\left(f\cdot f\right)=2\left(f{f}_{xxxx}-4f_x{f}_{xxx}+3f_{xx}^2\right)$ (2.12)

Hence Equation (2.4) becomes

$\left(D_x{D}_t+D_x^4\right)\left\{f\cdot f\right\}=P\left(D\right)\left\{f\cdot f\right\}=0$ (2.13)

which is the Hirota bilinear form of KdV Equation (2.1).

Application of the Hirota Perturbation: The supplementary function $f=f\left(x,t\right)$, in (2.13) is expressed as:

$f\left(x,t\right)=1+{ϵ}^1{f}_1+{ϵ}^2{f}_2+{ϵ}^3{f}_3+{ϵ}^4{f}_4+\cdots +{ϵ}^n{f}^n=1+{ϵ}^n{\displaystyle \underset{i=1}{\overset{n}{\sum }}{f}_i}$ (2.14)

where $f_1,f_2,\cdots ,f_n$, represent simple exponential functions.

We now insert (2.14) into Equation (2.13) so that

$\begin{array}{l}P\left(D\right)\left\{f\cdot f\right\}\\ =P\left(D\right)\left(1\cdot 1\right)+ϵP\left(D\right)\left(f_1\cdot 1+1\cdot f_1\right)+{ϵ}^{2}P\left(D\right)\left(f_2\cdot 1+f_1\cdot f_1+1\cdot f_2\right)\\ +{ϵ}^{3}P\left(D\right)\left(f_3\cdot 1+f_2\cdot f_1+f_1\cdot f_2+1\cdot f_3\right)+\cdots +{ϵ}^{n}P\left(D\right)\left({\displaystyle \underset{i=0}{\overset{n}{\sum }}{f}_i\cdot f_{n-i}}\right)\\ =0\end{array}$ (2.15)

The coefficient of like powers of $ϵ$ in (2.15) can be equated to zero to obtain the following sets of equations,

$\text{Ο}\left({ϵ}^n\right):P\left(D\right)\left({\displaystyle \underset{i=0}{\overset{n}{\sum }}{f}_i\cdot f_{n-i}}\right)=0,\left(f_0=1\right)$ (2.16)

We make use of the scheme (2.16) to obtain appropriate dispersion relations and coupling coefficients in the KdV Equation (2.13). Thus, if $f\left(x,t\right)$ in (2.14) is a solution of Equation (2.13), then $w\left(x,t\right)$ in Equation (2.3) is a soliton solution to the KdV Equation (2.1).

2.2. Constructing Multi-Soliton Solutions of KdV Equation

2.2.1. The Vacuum Soliton Solution

What we need to find now is a truncated supplementary function $f\left(x,t\right)$ in (2.14) that satisfies (2.13) which, when inserted into Equation (2.3), will yield $w\left(x,t\right)$, which is the solution of KdV Equation (2.1). From (2.14), if we try $f\left(x,t\right)=1$ and substitute this into Equation (2.3), we get the trivial solution $w\left(x,t\right)=0$. This solution is called the vacuum or zero soliton upon which multi-solitons solutions are obtained or propagated. This shows that a soliton can travel in a vacuum.

2.2.2. The Single Soliton Solution

The single soliton solution [2] [12] is given as,

$w\left(x,t\right)=\frac{1}{2}{k}^2{\text{sech}}^2\left[\frac{1}{2}k\left(x-ct\right)+\frac{1}{2}q\right]$ (2.17)

where $k$ is the wave number, c is wave velocity, and q is the initial point of propagation, $c=k^2$ is the dispersion relation. A wolfram Mathematica program is used to compute the single-soliton solution and its simulation. The profile of the single soliton at $q=0$, $t=0$, and $k=2$ is shown in Figure 1 and Figure 2.

Figure 1. 2D Profile of single soliton solution for KdV equation.

Figure 2. 3D profile of single soliton for KdV equation.

2.2.3. The Two-Soliton Solution

To obtain the two-soliton solution, the supplementary function f in (2.14) is truncated after the third term

$f\left(x,t\right)=1+ϵf_1+{ϵ}^2{f}_2,f_i=0,i\ge 3$ (2.18)

To find $f_1$ and $f_2$, use is made of a two-term form of $f_1={\text{e}}^v$ that is usually used for construction of the single-soliton case, i.e. $f_1={\text{e}}^{v_1}+{\text{e}}^{v_2}$ where $v_1=k_1\left(x-c_{1}t\right)+q_1$, $v_2=k_2\left(x-c_{2}t\right)+q_2$. Since two-soliton solution is built from single soliton, and one principle is that for integrable systems one must be able to combine any pair of single-soliton built on top of the same vacuum.

Now, using the perturbation scheme (2.16), we have:

$\begin{array}{l}P\left(D\right)\left(1\cdot f_1\right)+P\left(D\right)\left(f_1\cdot 1\right)\\ =P(\partial )f_1+P(\partial )f_1=2P(\partial )f_1=0\end{array}$ (2.19)

Thus, it follows that,

$\begin{array}{c}2P(\partial )f_1=2\left({\partial }_x{\partial }_t+{\partial }_x^4\right)\left\{{\text{e}}^{v_1}+{\text{e}}^{v_2}\right\}\\ =2\left({\partial }_x{\partial }_t+{\partial }_x^4\right)\left\{{\text{e}}^{v_1}\right\}+2\left({\partial }_x{\partial }_t+{\partial }_x^4\right)\left\{{\text{e}}^{v_2}\right\}\\ =0\end{array}$ (2.20)

This implies

$\left(2k_1^4-2c_1{k}_1^2\right){\text{e}}^{v_1}=0$ and $\left(2k_2^4-2c_2{k}_2^2\right){\text{e}}^{v_2}=0$ (2.21)

i.e. the dispersion relation is obtained as:

$c_i=k_i^2,i=1,2$ (2.22)

Now, we obtain $f_2$ from the perturbation scheme (2.16), as

$\begin{array}{l}P\left(D\right)\left(1\cdot f_2+f_2\cdot 1\right)+P\left(D\right)\left(f_1\cdot f_1\right)\\ =2P(\partial )f_2+P\left(D\right)\left(f_1\cdot f_1\right)=0\end{array}$ (2.23)

Thus,

$2P(\partial )f_2+P\left(D\right)\left\{\left({\text{e}}^{v_1}+{\text{e}}^{v_2}\right)\cdot \left({\text{e}}^{v_1}+{\text{e}}^{v_2}\right)\right\}=0$ (2.24)

or,

$2P(\partial )f_2-6k_1{k}_2{\left(k_1-k_2\right)}^2{\text{e}}^{v_1+v_2}=0$ (2.25)

Clearly, Equation (2.25) holds true, if and only if $f_2$ is of the form

$f_2=c_{12}{\text{e}}^{v_1+v_2}$ (2.26)

where $c_{12}$ is a coupling constant yet to be determined.

Now substituting (2.26) back into (2.25) we have that

$2c_{12}\left\{-\left(k_1+k_2\right)\left(k_1{c}_1+k_2{c}_2\right)+{\left(k_1+k_2\right)}^4\right\}{\text{e}}^{v_1+v_2}-6k_1{k}_2{\left(k_1-k_2\right)}^2{\text{e}}^{v_1+v_2}=0$ (2.27)

Substituting (2.22) in (2.27) and simplifying we have the coupling constant $c_{12}$ as:

$c_{12}=\frac{{\left(k_1-k_2\right)}^2}{{\left(k_1+k_2\right)}^2}$ (2.28)

Thus, $f_1$ and $f_2$ are determined as simple exponential functions in supplementary function in (2.18), without loss of generality, we set $ϵ=1$, and use (2.18) into (2.3) in constructing the two-soliton solution as

$\begin{array}{c}w\left(x,t\right)=\frac{2\left(k_1^2{\text{e}}^{v_1}+k_2^2{\text{e}}^{v_2}+c_{12}{\left(k_1+k_2\right)}^2{\text{e}}^{v_1+v_2}\right)}{\left(1+{\text{e}}^{v_1}+{\text{e}}^{v_2}+c_{12}{\text{e}}^{v_1+v_2}\right)}\\ -\frac{2{\left(k_1{\text{e}}^{v_1}+k_2{\text{e}}^{v_2}+c_{12}\left(k_1+k_2\right){\text{e}}^{v_1+v_2}\right)}^2}{{\left(1+{\text{e}}^{v_1}+{\text{e}}^{v_2}+c_{12}{\text{e}}^{v_1+v_2}\right)}^2}\end{array}$ (2.29)

where

$v_i=k_i\left(x-c_{i}t\right)+q_i,i=1,2$ and $c_{12}=\frac{{\left(k_1-k_2\right)}^2}{{\left(k_1+k_2\right)}^2}$ (2.30)

We now have a fully defined two-soliton solution for KdV Equation (2.2), A Wolfram Mathematica program was again used for computation and simulation of the two-soliton solution obtained in Equation (2.29). A study of this solution is given in Figures 3-8 where $k_1=1$, $k_2=1.5$ in Equation (2.29). It is clearly seen that the Solitons move left to right with the taller faster overtaking the shorter-slower soliton.

2.2.4. The Three-Soliton Solution

This process is similar to the two-soliton problem except that here we need to find $f_1,f_2$ and $f_3$ to form the supplementary function truncated from (2.14)

$f\left(x,t\right)=1+ϵf_1+{ϵ}^2{f}_2+{ϵ}^3{f}_3,f_i=0,i\ge 4$ (2.31)

And since we are seeking for a three-soliton solution, we use three-term forms of $f_1$ and $f_2$ that were successful for the two-soliton case, i.e.

$f_1={\text{e}}^{v_1}+{\text{e}}^{v_2}+{\text{e}}^{v_3}$ (2.32)

and

$f_2=c_{12}{\text{e}}^{v_1+v_2}+c_{13}{\text{e}}^{v_1+v_3}+c_{23}{\text{e}}^{v_2+v_3}$ (2.33)

where $c_{12},c_{13}$ and $c_{23}$ are coupling constants yet to be determined in terms of $k_1,k_2$ and $k_3$.

And we shall deal with finding $f_3$ in the computation process, similar to the $f_2$ in the two-soliton case.

Again, using the perturbation scheme (2.16), we have:

$P\left(D\right)\left(1\cdot f_1+f_1\cdot 1\right)=P\left(D\right)\left(1\cdot f_1\right)+P\left(D\right)\left(f_1\cdot 1\right)=2P(\partial )f_1=0$ (2.34)

Thus,

$\left(k_1^4-k_1^2{c}_1\right){\text{e}}^{v_1}+\left(k_2^4-k_2^2{c}_2\right){\text{e}}^{v_2}+\left(k_3^4-k_3^2{c}_3\right){\text{e}}^{v_3}=0$ (2.35)

and from (2.35) we obtain the dispersive relation:

$c_i=k_i^2,i=1,2,3$ (2.36)

Also

$P\left(D\right)\left(1\cdot f_2+f_2\cdot 1\right)+P\left(D\right)\left(f_1\cdot f_1\right)=2P(\partial )f_2+P\left(D\right)\left(f_1\cdot f_1\right)=0$ (2.37)

which implies,

Figure 3. 3D profile of 2-soliton solution.

Figure 4. 3D profile of 2-solitons from the top view.

Figure 5. Profile of 2 solitons before interactions.

Figure 6. Profile of 2 solitons as the taller soliton fast approaches for interactions.

Figure 7. Profile of 2 solitons during interactions.

Figure 8. Profile of 2 solitons just after interactions.

$\begin{array}{c}2P(\partial )f_2=\left(12c_{12}{k}_1^2{k}_2^2+6c_{12}{k}_1^3{k}_2+6c_{12}{k}_1{k}_2^3\right){\text{e}}^{v_1+v_2}\\ +\left(12c_{13}{k}_1^2{k}_3^2+6c_{13}{k}_1{k}_3^3+6c_{13}{k}_1^3{k}_3\right){\text{e}}^{v_1+v_3}\\ +\left(12c_{23}{k}_2^2{k}_3^2+6c_{23}{k}_2{k}_3^3+6c_{23}{k}_2^3{k}_3\right){\text{e}}^{v_2+v_3}\end{array}$ (2.38)

and,

$\begin{array}{c}P\left(D\right)\left(f_1\cdot f_1\right)=\left(12k_2^2{k}_1^2-6k_2{k}_1^3-6k_2^3{k}_1\right){\text{e}}^{v_1+v_2}\\ +\left(12k_1^2{k}_3^2-6k_3^3{k}_1-6k_3{k}_1^3\right){\text{e}}^{v_1+v_3}\\ +\left(12k_2^2{k}_3^2-6k_2{k}_3^3-6k_2^3{k}_3\right){\text{e}}^{v_2+v_3}\end{array}$ (2.39)

Adding Equations (2.38) and (2.39), and simplifying we have

$c_{12}=\frac{{\left(k_1-k_2\right)}^2}{{\left(k_1+k_2\right)}^2},c_{13}=\frac{{\left(k_1-k_3\right)}^2}{{\left(k_1+k_3\right)}^2},c_{23}=\frac{{\left(k_2-k_3\right)}^2}{{\left(k_2+k_3\right)}^2}$ (2.40)

We see here that all the coupling constants for $f_2$ are determined, we now have to find $f_3$.

Again

$P\left(D\right)\left(1\cdot f_3+f_3\cdot 1\right)+P\left(D\right)\left(f_1\cdot f_2+f_2\cdot f_1\right)=0$ (2.41)

Thus,

$P\left(D\right)\left(1\cdot f_3+f_3\cdot 1\right)-\frac{6{\left(k_1-k_2\right)}^2{\left(k_1-k_3\right)}^2{\left(k_2-k_3\right)}^2\left(k_1+k_2+k_3\right)}{(k_1+k_2)\left(k_1+k_3\right)\left(k_2+k_3\right)}{\text{e}}^{v_1+v_2+v_3}=0$ (2.42)

Clearly, Equation (2.42) holds if and only if, $f_3$ is of the form:

$f_3=B\cdot {\text{e}}^{v_1+v_2+v_3}$ (2.43)

Substituting (2.43) in the first term of Equation (2.42) we have

$\begin{array}{l}P\left(D\right)\left(1\cdot f_3+f_3\cdot 1\right)\\ =6B\left(k_1+k_2\right)\left(k_1+k_3\right)\left(k_2+k_3\right)\left(k_1+k_2+k_3\right){\text{e}}^{v_1+v_2+v_3}\end{array}$ (2.44)

Substituting (2.46) back into (2.44) and simplifying, the constant B is determined as:

$B=\frac{{\left(k_1-k_2\right)}^2{\left(k_1-k_3\right)}^2{\left(k_2-k_3\right)}^2}{{\left(k_1+k_2\right)}^2{\left(k_1+k_3\right)}^2{\left(k_2+k_3\right)}^2}=c_{12}\cdot c_{13}\cdot c_{23}$ (2.45)

Thus, we have successfully obtained $f_1,f_2$ and $f_3$ so that the supplementary function (2.31) by setting $ϵ=1$ becomes:

$f\left(x,t\right)=1+{\text{e}}^{v_1}+{\text{e}}^{v_2}+{\text{e}}^{v_3}+c_{12}{\text{e}}^{v_1+v_2}+c_{13}{\text{e}}^{v_1+v_3}+c_{23}{\text{e}}^{v_2+v_3}+B{\text{e}}^{v_1+v_2+v_2}$ (2.46)

Therefore substituting (2.46) in (2.3) we obtain the three-soliton solution as

$w=\frac{2\left(f{f}_{xx}-f_x^2\right)}{f^2}=\frac{2N\left(x,t\right)}{D\left(x,t\right)}$ (2.47)

where,

$\begin{array}{l}N\left(x,t\right)\\ ={\text{e}}^{v_1+v_2}\left[2{\left(k_1-k_2\right)}^2+2{\left(k_1-k_2\right)}^2{c}_{13}{c}_{23}{\text{e}}^{2v_1}+c_{12}{k}_1^2{\text{e}}^{v_2}+c_{12}{k}_2^2{\text{e}}^{v_1}\right]\\ +{\text{e}}^{v_1+v_3}\left[2{\left(k_1-k_3\right)}^2+2{\left(k_1-k_3\right)}^2{c}_{12}{c}_{23}{\text{e}}^{2v_2}+c_{13}{k}_1^2{\text{e}}^{v_3}+c_{13}{k}_3^2{\text{e}}^{v_1}\right]\\ +{\text{e}}^{v_2+v_3}\left[2{\left(k_2-k_3\right)}^2+2{\left(k_2-k_3\right)}^2{c}_{12}{c}_{23}{\text{e}}^{2v_1}+c_{23}{k}_2^2{\text{e}}^{v_3}+c_{23}{k}_3^2{\text{e}}^{v_2}\right]\\ +k_1^2{\text{e}}^{v_1}+k_2^2{\text{e}}^{v_2}+k_3^2{\text{e}}^{v_3}+B{\text{e}}^{v_1+v_2+v_3}\left[c_{12}{k}_3^2{\text{e}}^{v_1+v_2}+c_{13}{k}_2^2{\text{e}}^{v_1+v_3}+c_{23}{k}_1^2{\text{e}}^{v_2+v_3}\right]\\ +{\text{e}}^{v_1+v_2+v_3}[c_{12}\left(k_1^2+k_2^2+k_3^2+2k_1{k}_2-2k_1{k}_3-2k_2{k}_3\right)\\ +c_{13}\left(k_1^2+k_2^2+k_3^2+2k_1{k}_3-2k_1{k}_2-2k_2{k}_3\right)\\ +c_{23}\left(k_1^2+k_2^2+k_3^2+2k_2{k}_3-2k_1{k}_2-2k_1{k}_3\right)\\ +B\left(k_1^2+k_2^2+k_3^2+2k_2{k}_3-2k_1{k}_2-2k_1{k}_3\right)]\end{array}$ (2.48)

and

$D\left(x,t\right)={\left(1+{\text{e}}^{v_1}+{\text{e}}^{v_2}+{\text{e}}^{v_3}+c_{12}{\text{e}}^{v_1+v_2}+c_{13}{\text{e}}^{v_1+v_3}+c_{23}{\text{e}}^{v_2+v_3}+B{\text{e}}^{v_1+v_2+v_2}\right)}^2$ (2.49)

and $v_i=k_1\left(x-c_{i}t\right)+q_1,=1,2,3$ ; $c_{ij}=\frac{{\left(k_i-k_j\right)}^2}{{\left(k_i+k_j\right)}^2},i,j=1,2,3,i<j$ and $B=c_{12}\cdot c_{13}\cdot c_{23}$

Now we have a three-soliton solution for KdV Equation (2.1), a Wolfram Mathematica program was again used for computation and simulation of the three-solitons solution obtained in Equation (2.49). A study of this solution is given in Figures 9-13 respectively, with $k_1=1$, $k_2=2$, $k_3=3$ in equation (2.49). The solitons move left to right with the taller faster overtaking the shorter slower soliton.

Figure 9. 3D profile of 3 solitons before interactions.

Figure 10. Profile of 3 solitons before interactions.

Figure 11. Profile of 3 solitons during interactions.

Figure 12. Profile of 3 solitons just after interactions.

Figure 13. Profile of 2 solitons after interactions.

2.3. Numerical Solution of the KdV Equation Using Crank-Nicolson Method

Based on the Crank Nicolson scheme [13], we proposed to evaluate each term in the Equation (2.0) not only at time level $\left(n+\frac{1}{2}\right)$, but also at spatial location $\left(i+\frac{1}{2}\right)$ (i.e. at the midpoint of every subinterval) so as to ensure the first order

time derivative and third order space derivative (dispersive term) terms are both suitably centered. This ensures that we will be employing a cell-edge grid, but that the spatial finite differences in the scheme will be cell centered [12].

$w_t=\frac{w_i^{n+1}-w_i^n}{2\tau }+\frac{w_{i+1}^{n+1}-w_{i+1}^n}{2\tau }$ (2.50)

$w_{xxx}=\frac{w_{i+2}^n-3w_{i+1}^n+3w_i^n-w_{i-1}^n}{2h^3}+\frac{w_{i+2}^{n+1}-3w_{i+1}^{n+1}+3w_i^{n+1}-w_{i-1}^{n+1}}{2h^3}$ (2.51)

To carefully deal with the nonlinear term $\beta w{w}_x$, we assume that the leading coefficient $\beta w$ is known so as to escape the nonlinearity in the system of algebraic equations that will be obtained in the scheme. We designate it by $\beta \hat{w}$ and decide how to suitably estimate its value. Reviewing again that we have to

calculate at time level $\left(n+\frac{1}{2}\right)$ and space area $\left(i+\frac{1}{2}\right)$, the nonlinear term of the Equation (2.0) becomes:

$\beta w{w}_x=\frac{\beta {\hat{w}}_{i+1}+\beta {\hat{w}}_i}{4}\left(\frac{w_{i+1}^{n+1}-w_i^{n+1}}{h}+\frac{w_{i+1}^n-w_i^n}{h}\right)$ (2.52)

We note here that the derivatives in (2.52) are centered in time at $\left(n+\frac{1}{2}\right)$

while the $\hat{w}$ term is not. This will be estimated using predictor-corrector technique. Thus the KdV Equation (2.0) becomes a system of algebraic equations:

$\begin{array}{l}\frac{w_i^{n+1}-w_i^n}{2\tau }+\frac{w_{i+1}^{n+1}-w_{i+1}^n}{2\tau }+\frac{\beta {\hat{w}}_{i+1}+\beta {\hat{w}}_i}{4h}\left(w_{i+1}^{n+1}-w_i^{n+1}+w_{i+1}^n-w_i^n\right)\\ +\alpha \left(\frac{w_{i+2}^{n+1}-3w_{i+1}^{n+1}+3w_i^{n+1}-w_{i-1}^{n+1}}{2h^3}+\frac{w_{i+2}^n-3w_{i+1}^n+3w_i^n-w_{i-1}^n}{2h^3}\right)=0\end{array}$ (2.53)

Estimation in the Nonlinear Term using Predictor-Corrector Technique

We use a predictor-corrector technique [14] to properly center Crank-Nicolson in time between $t_n$ and $t_{n+1}$. Thus we approximately obtain $\hat{w}={\hat{w}}^{n+\frac{1}{2}}$ in this

procedure: The Crank Nicolson is applied two times in each time step. In the initial step (the predictor step) we essentially substitute $\hat{w}$ for $w^n$, the present estimated value of w, and call the subsequent new value (after Crank Nicolson is applied) ${\hat{w}}^{n+1}$, the predicted future value. In the second step we join this

predicted value with current value to approximately build ${\hat{w}}^{n+\frac{1}{2}}$ using $w^{n+\frac{1}{2}}=\frac{1}{2}\left(w^n+w^{n+1}\right)$, then reconstruct iteratively Crank Nicolson again.

Thus, the Algorithm for the Crank-Nicolson scheme (2.53) was developed and implemented in MATLAB [12] to study the behavior of solitons in the KdV equation in detail. Figures 14-18 show three solitons computed using this algorithm, computed at $\tau =0.01$, $\alpha =0.134$ and $\beta =1$.

3. Discussion

In Figure 1, Figures 5-8, and Figures 10-13, it is seen that Solitons move from left to right based on the positive sign effect of the nonlinear term of the KdV equation. In addition, Solitons in our simulations are propagated at different point locations so that the tallest one is located at the further left, but eventually move faster and over take the ones with smaller amplitudes.

In the numerical computation, The Crank-Nicolson implicit scheme was used to compute multi-soliton solutions in KdV equation. In Figures 14-18 we implemented the Crank-Nicolson algorithm in MATLAB, and studied the interactions of three solitons with soliton initial profile $6k_n^2\alpha {\text{sech}}^{2}k\left(x-x_{0n}\right),n=1,2,3$. Computed at $\tau =0.01$, $\alpha =0.134$, $\beta =1$, $k_1=1.5$, $k_2=2$, $k_3=2.5$ and

initial points $x_{01}=\frac{3L}{4}$, $x_{02}=\frac{L}{4}$, $x_{03}=\frac{L}{2}$ of propagation at an interval, L, at

four different times. We note here however that $\tau =0.01$ is chosen sufficiently small because of stability issues (to avoid blowup solutions which we have seen for larger time step $\tau$ ) in the approximation of the quantity $\hat{w}$ in the nonlinear term, estimated by predictor corrector technique.

Solitons in the numerical simulations move either left to right or right to left depending on the sign (positive or negative) in nonlinear term in KdV equation

Figure 14. Profile of 3 solitons comparing Analytical (Blue) and Numerical results (Red).

Figure 15. Profile of 3 solitons before interactions.

Figure 16. Profile of 3 solitons as the taller soliton fast approaches for interactions.

Figure 17. Profile of 3 solitons during interactions of 2 solitons.

Figure 18. Profile of 3 solitons after interactions of 2 solitons.

with the tallest one moving faster. It is also observed that as two solitons collide with each other they momentarily form a single soliton pulse as represented in Figure 7, Figure 11 and Figure 17, afterward they separate and continue with their original identities. This formation of a single pulse is however only typical of an interaction between solitons with sufficiently large difference in amplitudes. It is also noticed that during interactions of solitons in all the simulations carried out, the amplitude of the supposed single soliton pulse formed by two or more solitons is observed to be smaller than the amplitude of the larger soliton in the interactions. For instance, in Figure 14, the larger soliton amplitude is 10 units while the amplitude of single-soliton pulse formed is 7.5 units in Figure 16, but regains their heights or amplitudes after interactions. This similarly can also be seen in Figures 10-13 and Figures 5-8.

Again there is a phase shift after interactions, since the smaller soliton that is in front becomes behind and the larger or taller one becomes further head as opposed to linear waves. This confirms that solitons do not obey superposition principle but interact nonlinearly with each other. That notwithstanding we again checked to prove if superposition of two linear waves will form a soliton but this was not possible. Testing initial profile to be superposition of two linear waves, only a blowup solution appeared, justifying the nonlinearity in KdV equation and all models that admit this kind of solution are nonlinear in nature, balanced with dispersive term. The simulations results were observed with little ripples after collision particularly in the numerical computation after critical inspections as can be seen in figures16 and 18 but it has no significant effects to cause change in them.

4. Conclusions

In this paper, we have performed several computations both analytical and numerical. We obtained exact solutions of KdV equation via Hirota Method for one, two and three solitons. In order to ease computation process we wrote some simple computer codes to compute the three solutions. To study interactions or collision of solitons we preformed simulations using computer symbolic programming language—Wolfram Mathematica and MATLAB.

It was observed how two or more of solitons interact or collide with one another. It was fascinating to note in the study of analytical solution and numerical results obtained by Crank-Nicolson scheme that solitons pass through each other and come out unchanged. Thus we affirm that solitons indeed collide with one another and come without change of properties or identities. The analytical and numerical comparison of solutions agreement is very good.

We also uphold that in physical application, the study of Multi-Soliton solution in KdV equation gives us a clue why data receivers’ sets such as Radio sets get access to (signal transmitted easily during tuning) data and information transmission centers say Radio Stations with very high mast antennas faster than those with lower mast antennas since solitons of higher amplitude travel faster than solitons of shorter amplitude. A bit of noise of such data transmission can be perceived when different signal pulses transmitted collide but after, fade away as the tall signal pulse moves faster separating itself from the short signal pulse both of which however do not dissipate energy and continue with their identity and this is typical of all data transmission processes.

Conflicts of Interest

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

References

[1]	Miki, W. (2001) Introduction to Solitons. Pramana—Journal of Physics, 57, 841-847. https://doi.org/10.1007/s12043-001-0002-3
[2]	Rudy, L.H. (2002) A Brief Introduction to Soliton Theory in a Class of Nonlinear PDEs. Department of Mathematics and Computer Science, California State University, Hayward. https://www.researchgate.net/publication/280007998
[3]	Zabusky, N.J. and Kruskal, M.D. (1965) Interactions of “Solitons” in a Collisionless Plasma and the Recurrence of Initial States. Physics Review Letters, 15, 240. https://doi.org/10.1103/PhysRevLett.15.240
[4]	Dodd, R.K., Eilbeck, J.C., Gibbon, J.D. and Morris, H.C. (1982) Solitons and Nonlinear Wave Equations. Academic Press, New York.
[5]	Drazin, P.G. and Johnson, R.S. (1989) Solitons: An Introduction. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9781139172059
[6]	Stegeman, G.I. and Segev, M. (1999) Optical Spatial Solitons and Their Interactions: Universality and Diversity. Science, 286, 1518-1523. https://doi.org/10.1126/science.286.5444.1518
[7]	Pam, F.G. (2010) Solitons Could Power Molecular Electronics, Artificial Muscles. http://www.physorg.com/news71409967.html
[8]	Ryogo, H. (2004) The Direct Method in Soliton Theory. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511543043
[9]	Verheest, F., Hellberg, M.A. and Hereman, W.A. (2012) Head-On Collisions of Electrostatic Solitons in Non-Thermal Plasmas. Physical Review E, 82, Article ID: 057406. https://doi.org/10.1103/PhysRevE.86.036402
[10]	Masitah, S., Maureen, S.F.C. and Hajah, H.M.N. (2015) Applying Explicit Schemes to the Korteweg-deVries Equation. Modern Applied Science, 9, 200. https://doi.org/10.5539/mas.v9n4p200
[11]	Newell, A.C. (1985) Solitons in Mathematics and Physics. The Society for Industrial and Applied Mathematics (SIAM), Philadelphia, 25. https://doi.org/10.1137/1.9781611970227
[12]	Orapine, H.O. (2020) Comparative Study of the Analytical and Numerical Computations of Multi-Solitons in Korteweg-de Vries Model. Unpublished MSc. Thesis, Department of Mathematical Sciences, Nigerian Defence Academy, Kaduna.
[13]	Morton, K.W. and Mayers, D.F. (1994) Numerical Solutions of Partial Differential Equations. Cambridge University Press, Cambridge.
[14]	Conte, S.D. and de Boor, C. (1980) Elementary Numerical Analysis: An Algorithmic Approach. International Series in Pure and Applied Mathematics, the Solution of Differential Equation. McGraw-Hill, New York, 346-401.