# Differential Quadrature Method to Examine the Dynamical Behavior of Soliton Solutions to the Korteweg-de Vries Equation

Authors: Shubham Mishra; Geeta Arora; Homan Emadifar; Soubhagya Kumar Sahoo; Afshin Ghanizadeh
Journal: Advances in Mathematical Physics
 (2022)

---

## Abstract

Nonlinear evolution equations are crucial for understanding the phenomena in science and technology. One such equation with periodic solutions that has applications in various fields of physics is the Korteweg-de Vries (KdV) equation. In the present work, we are concerned with the implementation of a newly defined quintic B-spline basis function in the differential quadrature method for solving the Korteweg-de Vries (KdV) equation. The results are presented using four experiments involving a single soliton and the interaction of solitons. The accuracy and efficiency of the method are presented by computing theL2 and L∞ norms along with the conservational quantities in the forms of tables. The results show that the proposed scheme not only gives acceptable results but also consumes less time, as shown by the CPU for the elapsed time in two examples. The graphical representations of the obtained numerical solutions are compared with the exact solution to discuss the nature of solitons and their interactions for more than one soliton.


---

## Body

## 1. Introduction

While performing studies to identify the most effective design for canal boats on the Edinburgh-Glasgow canal in 1844, John Scott Russell noticed a phenomenon. He noticed that after one or two miles, the height of water in the canal steadily decreases as it travels along the watercourse. He invented the term “wave of translation” to describe this unique and wonderful phenomenon [1]. This gives rise to the soliton defined as a wave with a defined shape traveling at a constant speed through a given medium. The first wave to exhibit characteristics similar to a soliton was observed by Yuliawati et al. [2]. This was the beginning of an absolutely specific field of research to which scientists and mathematicians have contributed a lot over time. Nowadays, it is known that many equations have soliton solutions. Some of the equations having soliton solution are the KdV equation, Fisher equation, NLS equation, etc.The Korteweg-de Vries (KdV) equation is a nonlinear partial differential equation developed by Gardner and Morikawa in 1895 with respect to plasma waves [3] and then again by Washimi and Taniuti [4] to study acoustic waves in a cold plasma. The KdV equation is used to examine the propagation of low-amplitude water waves in shallow water bodies. The solution to this equation produces solitary waves [5].The KdV equation is given by(1)∂U∂tx,t+εUx,t∂U∂xx,t+μ∂3U∂x3x,t=0,a≤x≤b,t>0,where ε and μ are positive parameters and a,b represents the range under consideration. The KdV equation is a third-order nonlinear evolution equation that characterizes long waves and is widely used in physical and engineering disciplines. For example, it is used in modeling ionic-acoustic solitons in plasma physics [6], in the study of a long wave in subsurface oceans, and shallow sols in geophysical fluid dynamics [7, 8]. It also describes the phenomenon in cluster physics and superdeformed nuclei [9, 10], quantum field theory, and classical general relativity [11]. The solution of the KdV equation has opened enormous possibilities for mathematical concepts.The solutions of nonlinear equations are always of interest to researchers as they are studied using various approaches [12, 13]. In most cases, an analytical solution is not accessible, so numerical aspects are always necessary [14]. Gardner et al. [15] proved both the existence and uniqueness of solutions to the KdV equation. Liu [16] provided an elliptic Jacobi function solution for the KdV equation. In the same research paper, Hufford and Xing [17] reported a numerical solution for the linearized version of the problem as well as superconvergence for the approach used. Trogdon and Deconinck [18] presented a finite-genus solution to the equation. Grava and Klein [19] solved the KdV equation numerically and asymptotically for a small dispersion limit. Leach [20] gives the large-time evolution of the generalized Korteweg-de Vries equation. The wavelet Galerkin approach is used by Kumar and Mehra [21] to find a time-accurate solution of this equation. To solve this equation, Bahadir [22] uses an exponential finite difference technique. Aksan and Özdeş [23] use the Galerkin finite element approach with B-spline functions. This equation was solved numerically and analytically by Özer and Kutluay [24]. Ascher and McLachlan [25] provided a multisymplectic box technique for the KdV equation. Small time solutions of the equation were given by Kutluay et al. [26]. Idrees et al. [27] use the optimal homotopic asymptotic technique to solve this equation. To solve the KdV equation numerically, Gücüyenen and Tanoğlu [28] used the iterative splitting approach. Sarma [29] provided a solitary wave solution for this equation. Van de Fliert and Groesen [30] used a variational methodology, which was further investigated by Yuliawati et al. using the steepest descent approach, to study the solution of the KDV equation in the Hamiltonian condition. In addition, there have been several other successful numerical approaches to the KdV equation, including the spectral method [31], the pseudospectral method, and the collocation method [32].This paper is divided into the following sections. In Section2, the numerical scheme with the weight coefficient calculation procedure is discussed. Section 3 discusses numerical experiments and results, and Section 4 presents the final conclusion.


## 2. Numerical Scheme

Bellman et al. [33] first introduced the differential quadrature method (DQM) for the numerical solution of partial differential equations in 1972. Due to its simplicity, the approach has recently attracted much attention. The concept of the method is to use basis functions whose derivatives at the nodes are known [34]. Numerous researchers have used various test functions to construct different types of DQMs [35–38].The differential quadrature method involves estimating a derivative of a given function using linear summation of its components at different nodes of the problem domain. The domaina,b can be simply partitioned into uniformly distributed finite nodes xi with distance h, such that
(2)a=x0<x1<x2<⋯<xn−1<xn=b.LetBix be the quintic B-splines with knots at points xn, n=0,1,2,⋯,N. The arrangement of splines B−1,B0,B1,⋯,BN,BN+1 forms the basis for any function on a,b. For i=1,2,⋯,N, the solution at each time point of the node xi is Uxi,t. The estimated derivative parameters are calculated as follows:
(3)Ux=∑j=1Npijuxj,t,Uxx=∑j=1Nqijuxj,t,Uxxx=∑j=1Nrijuxj,t,for i=1,2,⋯,N. The derivatives are approximated by pij,qij,andrij. Once the values of pij are fixed as described in the next section, the weighting coefficients qijandrij can be easily calculated. The method for calculating the other coefficients is as follows:
(4)∂2ui∂x2=∂∂x∂u∂x=∑k=1Npik∂u∂xx=xk=∑k=1Npik∑j=1Npkjuxj,t=∑k=1N∑j=1Npikpkjuxj,t=∑j=1Nqijuxj,t,i=1,2,3,⋯,N.Sinceqij is calculated using pij,rij can be calculated in a similar manner.Fori=−2,−1,0,⋯,N+2, Bix, the quintic B-spline basis function, describes a piecewise-defined function with the properties of continuity and division of unity.The following equations can be used to calculate the basis functions.(5)Bix=1h5x−xi−35,x∈xi−3,xi−2,x−xi−35−6x−xi−25,x∈xi−2,xi−1,x−xi−35−6x−xi−25+15x−xi−15,x∈xi−1,xi,xi+3−x5−6xi+2−x5+15xi+1−x5,x∈xi,xi+1,xi+3−x5−6xi+2−x5,x∈xi+1,xi+2,xi+3−x5,x∈xi+2,xi+3,0,otherwise,where B−2,B−1,B0,B1,⋯,0i+1,Bi+2 are the bases formed over the region a≤x≤b. Each quintic B-spline covers six elements, so that a total of six quintic B-splines cover one element. Table 1 summarizes the values of Bix and the first four derivatives.Table 1
Values ofBix and its derivatives at the nodes [34].
xxi−3xi−2xi−1xixi+1xi+2xi+3Bix0126662610Bi′x05/h50/h0−50/h5/h0Bi″x020/h240/h2−120/h240/h220/h20Bi‴x060/h3−120/h30120/h3−60/h30Biivx0120/h4−480/h4720/h4−480/h4120/h40The first-order approximation of the derivative can be estimated using the following relation:(6)Bi′xi=∑j=1NpijBixj,fori=1,2,⋯,N.As a result, a matrix system emerges as follows:(7)Ap⟶i=s⟶i.Here,A is the coefficient matrix given by
(8)66261000·0266626100·0126662610·0········0·0126662610·0012666260·00012666,representing the vector, corresponding to node point xi. The unknown coefficients are p⟶i=pi1,pi2,⋯,piNT, i=1,2,⋯,N, with the right-hand side given as follows:
(9)s⟶1=0,f,g,0,⋯,0T,s⟶2=−f,0,f,g,0,⋯,0T,s⟶3=−g,−f,0,f,g,⋯,0T,··s⟶N−2=0,⋯,−g,−f,0,f,gT,s⟶N−1=0,⋯,0,−g,−f,0,fT,s⟶N=0,⋯,0,−g,−f,0T.Here,f=50/h and g=5/h.The coefficientspi1,pi2,⋯,piN for i=1,2,⋯,N were calculated using MATLAB 2014 to solve the given five-band matrix system. Substituting approximate values for the derived first- and third-order spatial derivatives in equation (1) yields the following system:
(10)ut=−ϵu∑j=1Npijuj−μ∑j=1Nrijuj.The SSP-RK43 scheme [39] is then used to solve this system of ordinary differential equations, which offers numerical solutions at various time levels.


## 3. Numerical Experiments

In this section, the accuracy of the proposed method is shown by calculating theL2andL∞ errors defined as follows:
(11)L2=Uex−UN2≃h∑j=1NUjex.−UNj2,L∞=Uex−UN∞≈maxjUjex.−UNj,j=1,2,⋯,N−1.The lowest three invariants related to mass, momentum, and energy conservation are also calculated by the following equations:(12)I1=∫abUdx,I2=∫cbU2dx,I3=∫abU3−3μεU′2dx.


### 3.1. Experimental of Evaluation of a Single Soliton

Consider the KdV equation with the exact solution given as [40] follows:
(13)Ux,t=3Csech2Ax−Bt+D,Here,(14)A=12εCμ1/2,B=12εCεCμ1/2,so that (13) offers a single soliton with amplitude 3C and velocity εC moving towards the right.The equation is solved with the initial state taken from analytic solution (13) as follows:
(15)Ux,0=3Csech2Ax+D,and the boundary conditions U0,t=U2,t=0 for t≥0. ε=1,μ=4.84×10−4,C=0.3,D=−6 is employed in order to create a comparison with other investigations. To demonstrate the evolution of the current technique using a modified quintic B-spline DQM, Tables 2 and 3 show the error norm and invariant values, respectively. Moreover, at different values of ∆t=0.0005and0.001, numerical and exact solutions are represented by Figures 1 and 2, respectively.Table 2
Experimental evaluation of single soliton:∆t=0.0005.
MethodNTL2L∞I1I2I3Present scheme1510.252.3524×10−66.4229×10−60.14460.08680.04690.504.0430×10−61.2235×10−50.14460.08680.04690.756.0076×10−61.9232×10−50.14460.08680.04691.008.1958×10−62.5869×10−50.14460.08680.04692.004.3005×10−59.0073×10−50.14460.08680.04693.008.3086×10−40.0022×10−60.14440.08680.0469MQ_DQM [41]2010.251.01×10−52.66×10−50.14450.08670.04680.501.11×10−52.59×10−50.14450.08670.04680.751.33×10−53.94×10−50.14450.08670.04681.001.43×10−54.08×10−50.14450.08670.04682.002.14×10−56.74×10−50.14450.08670.04683.002.86×10−58.15×10−50.14460.08670.0468Table 3
Experimental evaluation of single soliton:∆t=0.001.
MethodNTL2L∞I1I2I3Present scheme910.256.2478×10−52.1816×10−40.14460.08680.04690.509.6636×10−52.0265×10−40.14460.08680.04690.751.9160×10−43.0993×10−40.14460.08680.04691.004.4543×10−47.0164×10−40.14460.08680.0469MQ_DQM [41]2010.250.0000100.0000270.14450.08670.04680.500.0000100.0000210.14450.08670.04680.750.0000120.0000340.14450.08670.04681.000.0000120.0000320.14450.08670.0468[42]2000.250.00522—0.500.01200—0.1445900.0867590.0468710.750.01220—1.000.02220—0.1445900.0867590.046873Figure 1
Simulations of single solitons:∆t=0.0005.Figure 2
Simulations of single solitons:∆t=0.001.


### 3.2. Experimental Evaluation of the Interaction of Two Solitons

Consider this second experiment [43] with the initial condition stated as follows:
(16)U=∑i=123Cisech2Aix+xi,Ai=12εCiμ1/2,i=1,2,with boundary conditions
(17)U0,t=U2,t=0,where ε=1,μ=4.84×10−4,C1=0.3,C2=0.1,x1=x2=−6 is considered in all simulations. The same parameters as in the previous study [43] are used for numerical calculations using MATLAB R2015b (32 bit) in Windows 10 version 21H2 for 64x, with N=91 and ∆t=0.005. Table 4 displays the error norm and invariant value. Moreover, Figure 3 demonstrates the numerical solution at different values of T.Table 4
Experimental evaluation of interaction of two solitons:∆t=0.005.
MethodNTI1I2I3CPU time (sec)Present scheme910.750.22810.10710.05330.2081.500.22790.10710.05330.2442.250.22780.10710.05330.2833.000.22380.10740.05330.316MQ_DQM [41]910.750.22810.10700.05331.500.22800.10700.05332.250.22790.10700.05333.000.22770.10700.05332000.750.22800.10700.05351.500.22800.10700.05343.000.22790.10700.0532Figure 3
Simulations of two solitons:∆t=0.005.


### 3.3. Experimental Evaluation of Interaction of Three Solitons

The numerical solution is calculated for the interaction of three solitons having the initial condition [44] given as follows:
(18)Ux,0=∑i=1312Ci2sech2Cix−xi,with the zero boundary conditions for domain [−100,100] with ε=1.0,μ=1.0,C1=0.3,C2=0.25,C3=0.2,x1=−60,x2=−44,x3=−26. The same parameters as in the previous study [44] at ∆t=0.1 and a much smaller number of grid points N=251 than those in the previous study [44] N=481 were used in the numerical computations. Table 5 displays the error norm and invariant value, and the numerical solution at different values of T is presented in Figure 4.Table 5
Experimental evaluation of interaction of three solitons:∆t=0.1.
MethodNTI1I2I3Present scheme2515618.00049.82745.261511217.99919.82755.263316818.00299.82755.262722418.00989.82765.262428018.14459.89895.2627MQ_DQM [41]4815618.00029.82735.262211217.99949.82735.262116817.99899.82745.262322417.99889.82745.262328018.00069.82745.2623MQ20005618.00189.59365.032811217.99749.51384.965116817.99719.32284.780822417.99859.06974.536228017.99958.83274.3141Figure 4
Simulations of three solitons:∆t=0.1.


### 3.4. Experimental Evaluation of Interaction of Four Solitons

In this example, the interaction of four solitons is presented with initial condition [44] given as follows:
(19)U=∑i=1412Ci2sech2Cix−xi,along with zero boundary conditions for domain [−150,150] with ε=1.0,μ=1.0,C1=0.3,C2=0.25,C3=0.2,C4=0.15,x1=−85,x2=−60,x3=−35,x4=−10.The considered parameters are the same as those in the previous study [44] at ∆t=0.1, and a much smaller number of grid points N=401 than that in the previous study [44] N=451 were used in the numerical computations. Table 6 displays the error norm and invariant value at different time levels with the physical behavior being shown in Figure 5 for different values of T.Table 6
Experimental evaluation of the interaction of four solitons:∆t=0.1.
MethodNTI1I2I3CPU time (sec)Present scheme4018021.600010.38875.26872.22516021.599910.38875.26833.81924021.600210.38875.27025.31632021.599210.38875.26917.26340021.536010.39735.27079.173MQ_DQM [41]4518021.600010.38875.268816021.600010.38865.268724021.600010.38865.268832021.599810.38875.268840021.600010.38875.26881500802.160289.97234.859416021.60499.74484.642624021.60119.70234.603532021.60079.47744.394340021.60749.19224.1368Figure 5
Simulations of four solitons at∆t=0.1.


## 3.1. Experimental of Evaluation of a Single Soliton

Consider the KdV equation with the exact solution given as [40] follows:
(13)Ux,t=3Csech2Ax−Bt+D,Here,(14)A=12εCμ1/2,B=12εCεCμ1/2,so that (13) offers a single soliton with amplitude 3C and velocity εC moving towards the right.The equation is solved with the initial state taken from analytic solution (13) as follows:
(15)Ux,0=3Csech2Ax+D,and the boundary conditions U0,t=U2,t=0 for t≥0. ε=1,μ=4.84×10−4,C=0.3,D=−6 is employed in order to create a comparison with other investigations. To demonstrate the evolution of the current technique using a modified quintic B-spline DQM, Tables 2 and 3 show the error norm and invariant values, respectively. Moreover, at different values of ∆t=0.0005and0.001, numerical and exact solutions are represented by Figures 1 and 2, respectively.Table 2
Experimental evaluation of single soliton:∆t=0.0005.
MethodNTL2L∞I1I2I3Present scheme1510.252.3524×10−66.4229×10−60.14460.08680.04690.504.0430×10−61.2235×10−50.14460.08680.04690.756.0076×10−61.9232×10−50.14460.08680.04691.008.1958×10−62.5869×10−50.14460.08680.04692.004.3005×10−59.0073×10−50.14460.08680.04693.008.3086×10−40.0022×10−60.14440.08680.0469MQ_DQM [41]2010.251.01×10−52.66×10−50.14450.08670.04680.501.11×10−52.59×10−50.14450.08670.04680.751.33×10−53.94×10−50.14450.08670.04681.001.43×10−54.08×10−50.14450.08670.04682.002.14×10−56.74×10−50.14450.08670.04683.002.86×10−58.15×10−50.14460.08670.0468Table 3
Experimental evaluation of single soliton:∆t=0.001.
MethodNTL2L∞I1I2I3Present scheme910.256.2478×10−52.1816×10−40.14460.08680.04690.509.6636×10−52.0265×10−40.14460.08680.04690.751.9160×10−43.0993×10−40.14460.08680.04691.004.4543×10−47.0164×10−40.14460.08680.0469MQ_DQM [41]2010.250.0000100.0000270.14450.08670.04680.500.0000100.0000210.14450.08670.04680.750.0000120.0000340.14450.08670.04681.000.0000120.0000320.14450.08670.0468[42]2000.250.00522—0.500.01200—0.1445900.0867590.0468710.750.01220—1.000.02220—0.1445900.0867590.046873Figure 1
Simulations of single solitons:∆t=0.0005.Figure 2
Simulations of single solitons:∆t=0.001.


## 3.2. Experimental Evaluation of the Interaction of Two Solitons

Consider this second experiment [43] with the initial condition stated as follows:
(16)U=∑i=123Cisech2Aix+xi,Ai=12εCiμ1/2,i=1,2,with boundary conditions
(17)U0,t=U2,t=0,where ε=1,μ=4.84×10−4,C1=0.3,C2=0.1,x1=x2=−6 is considered in all simulations. The same parameters as in the previous study [43] are used for numerical calculations using MATLAB R2015b (32 bit) in Windows 10 version 21H2 for 64x, with N=91 and ∆t=0.005. Table 4 displays the error norm and invariant value. Moreover, Figure 3 demonstrates the numerical solution at different values of T.Table 4
Experimental evaluation of interaction of two solitons:∆t=0.005.
MethodNTI1I2I3CPU time (sec)Present scheme910.750.22810.10710.05330.2081.500.22790.10710.05330.2442.250.22780.10710.05330.2833.000.22380.10740.05330.316MQ_DQM [41]910.750.22810.10700.05331.500.22800.10700.05332.250.22790.10700.05333.000.22770.10700.05332000.750.22800.10700.05351.500.22800.10700.05343.000.22790.10700.0532Figure 3
Simulations of two solitons:∆t=0.005.


## 3.3. Experimental Evaluation of Interaction of Three Solitons

The numerical solution is calculated for the interaction of three solitons having the initial condition [44] given as follows:
(18)Ux,0=∑i=1312Ci2sech2Cix−xi,with the zero boundary conditions for domain [−100,100] with ε=1.0,μ=1.0,C1=0.3,C2=0.25,C3=0.2,x1=−60,x2=−44,x3=−26. The same parameters as in the previous study [44] at ∆t=0.1 and a much smaller number of grid points N=251 than those in the previous study [44] N=481 were used in the numerical computations. Table 5 displays the error norm and invariant value, and the numerical solution at different values of T is presented in Figure 4.Table 5
Experimental evaluation of interaction of three solitons:∆t=0.1.
MethodNTI1I2I3Present scheme2515618.00049.82745.261511217.99919.82755.263316818.00299.82755.262722418.00989.82765.262428018.14459.89895.2627MQ_DQM [41]4815618.00029.82735.262211217.99949.82735.262116817.99899.82745.262322417.99889.82745.262328018.00069.82745.2623MQ20005618.00189.59365.032811217.99749.51384.965116817.99719.32284.780822417.99859.06974.536228017.99958.83274.3141Figure 4
Simulations of three solitons:∆t=0.1.


## 3.4. Experimental Evaluation of Interaction of Four Solitons

In this example, the interaction of four solitons is presented with initial condition [44] given as follows:
(19)U=∑i=1412Ci2sech2Cix−xi,along with zero boundary conditions for domain [−150,150] with ε=1.0,μ=1.0,C1=0.3,C2=0.25,C3=0.2,C4=0.15,x1=−85,x2=−60,x3=−35,x4=−10.The considered parameters are the same as those in the previous study [44] at ∆t=0.1, and a much smaller number of grid points N=401 than that in the previous study [44] N=451 were used in the numerical computations. Table 6 displays the error norm and invariant value at different time levels with the physical behavior being shown in Figure 5 for different values of T.Table 6
Experimental evaluation of the interaction of four solitons:∆t=0.1.
MethodNTI1I2I3CPU time (sec)Present scheme4018021.600010.38875.26872.22516021.599910.38875.26833.81924021.600210.38875.27025.31632021.599210.38875.26917.26340021.536010.39735.27079.173MQ_DQM [41]4518021.600010.38875.268816021.600010.38865.268724021.600010.38865.268832021.599810.38875.268840021.600010.38875.26881500802.160289.97234.859416021.60499.74484.642624021.60119.70234.603532021.60079.47744.394340021.60749.19224.1368Figure 5
Simulations of four solitons at∆t=0.1.


## 4. Conclusion

Due to the numerous applications of the KDV equation in the physical phenomena, in recent years, this equation has become a point of attraction for the researchers who want to find a numerical solution for this equation using various methods. In this paper, the newly defined quintic B-spline basis function is presented to solve the equation using the differential quadrature method. The advantage of this approach is involved in transforming the partial differential equation to an ordinary differential equation which can be solved by any numerical technique for the solution of the ordinary differential equation. In the present work, the SSP-RK43 is implemented to solve the obtained system of the ordinary differential equation which is a combination of the RK method of orders four and five that is a strong-stability-preserving scheme. Numerical results in terms of conservation variables and errors are calculated for the single soliton and extended till interaction of four solitons. The results are compared with numerical solutions from the literature. The obtained results agree well with those obtained earlier. The advantage of the proposed method is its ease of implementation compared to the previous methods. Thus, the present approach can be utilized to solve a variety of nonlinear physical models with extension and application to two-dimensional problems.



---

*Source: 8479433-2022-07-09.xml*