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   .
 and  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 . 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.  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  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.  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
where 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.  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 .
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.
with initial condition
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 ( ) and its nonlinear term ( ). 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 .
Bilinearization: First the KdV Equation (2.1) is rewritten as
The supplementary function (the Cole-Hopf transformation) given as
is used to transform (2.1) into a bilinear KdV equation. Thus, taking time derivative of 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,
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 , which is a binary operator defined on ordered pairs of functions and , as follows
where are nonnegative integers and
More generally we denote some sort of combination of the Hirota D-operator as a polynomial of D-operator i.e.
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:
Now replacing g with f to have the same function in (2.4) we have that:
Similarly, from (2.9) we have that
Hence Equation (2.4) becomes
which is the Hirota bilinear form of KdV Equation (2.1).
Application of the Hirota Perturbation: The supplementary function , in (2.13) is expressed as:
where , represent simple exponential functions.
We now insert (2.14) into Equation (2.13) so that
The coefficient of like powers of in (2.15) can be equated to zero to obtain the following sets of equations,
We make use of the scheme (2.16) to obtain appropriate dispersion relations and coupling coefficients in the KdV Equation (2.13). Thus, if in (2.14) is a solution of Equation (2.13), then 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 in (2.14) that satisfies (2.13) which, when inserted into Equation (2.3), will yield , which is the solution of KdV Equation (2.1). From (2.14), if we try and substitute this into Equation (2.3), we get the trivial solution . 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   is given as,
where is the wave number, c is wave velocity, and q is the initial point of propagation, 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 , , and 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
To find and , use is made of a two-term form of that is usually used for construction of the single-soliton case, i.e. where , . 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:
Thus, it follows that,
i.e. the dispersion relation is obtained as:
Now, we obtain from the perturbation scheme (2.16), as
Clearly, Equation (2.25) holds true, if and only if is of the form
where is a coupling constant yet to be determined.
Now substituting (2.26) back into (2.25) we have that
Substituting (2.22) in (2.27) and simplifying we have the coupling constant as:
Thus, and are determined as simple exponential functions in supplementary function in (2.18), without loss of generality, we set , and use (2.18) into (2.3) in constructing the two-soliton solution as
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 , 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 and to form the supplementary function truncated from (2.14)
And since we are seeking for a three-soliton solution, we use three-term forms of and that were successful for the two-soliton case, i.e.
where and are coupling constants yet to be determined in terms of and .
And we shall deal with finding in the computation process, similar to the in the two-soliton case.
Again, using the perturbation scheme (2.16), we have:
and from (2.35) we obtain the dispersive relation:
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.
Adding Equations (2.38) and (2.39), and simplifying we have
We see here that all the coupling constants for are determined, we now have to find .
Clearly, Equation (2.42) holds if and only if, is of the form:
Substituting (2.43) in the first term of Equation (2.42) we have
Substituting (2.46) back into (2.44) and simplifying, the constant B is determined as:
Thus, we have successfully obtained and so that the supplementary function (2.31) by setting becomes:
Therefore substituting (2.46) in (2.3) we obtain the three-soliton solution as
and ; and
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 , , 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 , we proposed to evaluate each term in the Equation (2.0) not only at time level , but also at spatial location (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 .
To carefully deal with the nonlinear term , we assume that the leading coefficient 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 and decide how to suitably estimate its value. Reviewing again that we have to
calculate at time level and space area , the nonlinear term of the Equation (2.0) becomes:
We note here that the derivatives in (2.52) are centered in time at
while the term is not. This will be estimated using predictor-corrector technique. Thus the KdV Equation (2.0) becomes a system of algebraic equations:
Estimation in the Nonlinear Term using Predictor-Corrector Technique
We use a predictor-corrector technique  to properly center Crank-Nicolson in time between and . Thus we approximately obtain in this
procedure: The Crank Nicolson is applied two times in each time step. In the initial step (the predictor step) we essentially substitute for , the present estimated value of w, and call the subsequent new value (after Crank Nicolson is applied) , the predicted future value. In the second step we join this
predicted value with current value to approximately build using , then reconstruct iteratively Crank Nicolson again.
Thus, the Algorithm for the Crank-Nicolson scheme (2.53) was developed and implemented in MATLAB  to study the behavior of solitons in the KdV equation in detail. Figures 14-18 show three solitons computed using this algorithm, computed at , and .
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 . Computed at , , , , , and
initial points , , of propagation at an interval, L, at
four different times. We note here however that is chosen sufficiently small because of stability issues (to avoid blowup solutions which we have seen for larger time step ) in the approximation of the quantity 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.
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.
 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.
 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.
 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.
 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.
 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.