Numerical study of soliton behavior of generalised Kuramoto-Sivashinsky type equations with Hermite splines

1.   Introduction

Nonlinear higher-order partial differential equations play an essential role in the physical processes occurring in various sciences and engineering disciplines. The fluid dynamical waves in oceanography and bio-fluid dynamics are usually framed by nonlinear higher-order perturbed forms of partial differential equations. In general, the partial differential problems can be of the form:

The nonlinear higher-order differential problems have been studied so far with various analytical and numerical techniques. Yet, there is no straightforward technique to thoroughly analyze any family of differential problems. The sub-equation method and the generalized Kudryashov method ^[1] discussed the soliton-type solutions of evolution equations. The Homotopy analysis J-transformation method is discussed in ^[2] to analyze the turbulent solutions of the generalized Kuramoto-Sivashinsky equation with the fractional operator. The compact solution is the advantage of any analytic technique but is restricted to specific conditions to study nonlinear differential problems. Various numerical techniques have been designed, such as Legendre wavelets methods ^[3], a formulation of the fifth order Korteweg-de Vries (KdV) equation in the form of Evans function ^[4], a rational spectral method ^[5], mesh-free collocation method ^[6], etc.

In this study, a prominent position occupying fourth-order nonlinear problems is elaborated on with respect to the perturbation parameter.

subjected to the initial and boundary data as follows:

where $\kappa, \delta$ and $\varepsilon$ are the known constant coefficients and $\alpha(u)$ represents a function of $u(\zeta, t)$, describing the physical processes. The functions $\phi(\zeta)$ and $\psi_i (t) : i = 1, 2, 3, 4$ are the known continuous functions restricting $u(\zeta, t)$ at initial time and the spatial boundary. For $\delta = 1$ and $\alpha(u) = u$, the equation(1.2) was introduced by Benny in the study of long waves on thin fluid films ^[7] and is referred to as a special case of the Benny equation. Nowadays, the problem (1.2) has also been studied as the generalized Kuramoto-Shivashinsky equation by the Lattice Boltzmann method ^[8], collocation method with B-splines functions^[9], methods of lines accompanied by radial basis functions ^[10], moving least square meshless method ^[11], and the Chebyshev spectral collocation method ^[12]. The proposed problem has also been considered as the KdV-Burgers-Kuramoto equation in ^[13] considering the traveling fronts for small dissipation.

The proposed problem exhibits rich solutions ranging from stable short waves to turbulent shock waves, and also has applications in various disciplines like the flow of turbulent air over the laminar fluid^[14] and solidification of alloys ^[15]. A numerical study of flame dynamics and instabilities represented by the Kuramoto-Sivashinsky type equation has been conducted by ^[16]. A numerical study accounting for the thin film growth and its roughness modeled by the conserved Kuramoto-Sivashinsky equation has been presented in ^[17]. Earlier, the exact form of solution has been obtained for the aforementioned problem using the Weiss—Tabor—Carnevale method ^[18]. Various other aspects of the problem, like the existence of the solution, its similar solutions and the structure of the traveling wave solution, and the Gevrey regularity and the stability of the solution, have been discussed by ^[19,20,21]. The evolution of one or more traveling pulses with the dominating dispersion in Eq (1.2) has been described by ^[22]. Solitary wave solutions due to dissipation as well as instability has been derived by ^[23]. In higher dimensions, ^[24], the existence and nonlinear stability of solution for Eq (1.2) have been elaborated. A study using equivalence transformation leading to an algebraic system of equations has also been conducted in ^[25] for the proposed problem (1.2).

Earlier, various forms of the traveling waves have been studied by researchers regarding the Kuramoto-Sivashinsky (KS) equation^[26,27], KdV equation^[28,29], KdV Burgers equation^[30], and Benjamin–Bona–Mahony–Burgers' equation ^[31] with conservative compact scheme.

The aforementioned problem incorporates the effect of dispersion as well as the dissipation process due to the Kuramoto-Shivashinsky equation. The addition of dispersion($\delta > 0$) changes the time evolution of the wave significantly.

The balancing effect of both the third and fourth-order derivative terms corresponding to the different initial conditions, especially the $\text{sech}$ function, is thoroughly elaborated concerning the perturbed parameters. It is observed that as the perturbation parameter $\varepsilon \rightarrow 0.5, 0.1$, humps with higher amplitude, yet bounded, and starts to appear due to suppressed dissipative effect only on the negative domain, then for the nonnegative domain, the wave saves its pattern irrespective to the values of $\varepsilon$. Also, the problem was studied with dominant dissipation, where permanent short waves were observed. A detailed study of the dissipation system is found in ^[32,33].

In this manuscript, the traveling waves under the influences of dissipation and dispersion are investigated. The rest of the study is manipulated as follows: In Section(2), a weighted finite difference scheme has been studied. Section(3) discusses the fully discretization of the proposed technique. Section(4) comprises the stability analysis using the von Neumann method. The numerical illustrations are presented in Section(5). The conclusion of the entire study is discussed in Section(6).

2.   Space splitting and semi discretization technique

In order to meet the requirement that the approximate solution should satisfy the boundary conditions imposed on the solution, the fourth order equation is reduced to a coupled system of second order by assuming a differentiable function $w = -u_{\zeta\zeta}$^[34]. The reduced system reads:

To numerically approximate the solution, the coupled system of equations, after reduction of the order of the proposed problem, is treated with the Crank-Nicholson scheme for time variable coupled with the collocation method for spatial approximation. Earlier, the traveling wave behavior of the Kuramoto-Shivashinsky equation ^[35], Burgers Huxley and Burgers Fisher ^[36] and Benjamin-Bona-Mahony-Burgers equations ^[37], singularly perturbed problems^[38], etc. have been studied by the weighted finite difference scheme with quintic Hermite spline collocation method. Here, the simultaneous impact of the dispersive and dissipative mechanisms on the traveling waves and the corresponding perturbed problems have been considered.

In the present study, the Crank-Nicolson (CN) scheme as proposed by ^[39] operates on the mean of function over the interval $[t_{j}\; \; t_{j+1}]$ with uniform distribution of points ^[35,38].

Define the partition $\pi_{t}:0 = t_{0} < t_{1} < \ldots < t_{M} = T$ with $\tau = t_{j+1}-t_{j}$. After applying the CN scheme on the coupled system of equations given by Eq (2.1), following system of equations appears:

and the nonlinear function $f(u, u_\zeta)^{j+1} = (\alpha(u)u_\zeta)^{j+1}$ is quasi-linearized with truncation error of order $\mathcal{O}(\tau^2)$ by

The non-linearity is treated with different formulations using implicit and compact finite difference schemes in ^[40,41]. Also, a special case of quasi-linearization of the nonlinear term, i.e., $\alpha(u) = u$, is studied by ^[42]. After rearranging the terms, the following system is obtained:

For convenience, write $u(\zeta, t_{j+1}) = u^{j+1}$.

3.   Fully discretization technique

Orthogonal splines are the piecewise orthogonal polynomials that interpolate the function at node points. Hermite splines are orthogonal splines and are usually followed to interpolate the function. Earlier, the essence of Hermite splines has been discussed in many ways. Hermite splines are considered to be the extension of the Lagrangian interpolating polynomials. Hermite interpolating polynomials of order $'2k+1'$ interpolate the function and its $k^{th}$ order derivative at node points. This feature of Hermite interpolating polynomials makes it superior to Lagrangian interpolating polynomials. In the present study, quintic Hermite splines, which are of order 5, are followed. Piecewise, the basis includes six splines and interpolates the function as well as its first and second-order derivatives at intermediate node points. The detailed study of quintic Hermite splines is given hereunder.

Consider an interval $I = [a, b]$, let $\Pi$ be the partition of $I$ such that:

$\Pi: a = \zeta_0 < \zeta_1 < \zeta_2 < ... < \zeta_{m-1} < \zeta_m = b$.

Let $I_\gamma = [\zeta_{\gamma-1}, \zeta_\gamma]$ such that $\gamma = 1(1)m$ is the $\gamma^{\text{th}}$ subinterval of the partitioned domain. Let $h = \zeta_{\gamma}-\zeta_{\gamma-1}: 1\leq \gamma \leq m$ represent the constant step size of the partition $\Pi$. The piecewise interpolating Hermite splines are represented as follows:

In order to implement the collocation, the Hermite splines within each sub-domain $I_\gamma$ are mapped onto the interval $[0, 1]$ using a linear mapping $\xi = \frac{\zeta-\zeta_{\gamma-1}}{h}$. This linear map reduces the structure of splines mentioned in (3.1) to the tabular form given in Table (1) for $0 \leq \xi\leq 1$.

The Hermite spline collocation method is applied to both the functions $u^n, w^n$ simultaneously. Thus, assuming the following piecewise approximations:

where $a^{\gamma_i, n}$, $b^{\gamma_k, n}$ are the controlling constant coefficients to be determined. To implement the orthogonal collocation, the roots of orthogonal polynomials are taken as collocation points within each sub-domain. Traditionally, the zeros of Jacobi polynomials are taken as collocation points. The zeros of Chebyshev or Legendre polynomials of degree six, which are special cases of Jacobi polynomials, are usually chosen as collocation points. In the present study, the zeros of shifted Legendre polynomials of degree six have been taken as collocation points ^[43]. Using the piecewise approximations defined in Eqs (3.2) and (3.3) in the semi-discretized coupled Eqs (2.5) and (2.6), the following system of equations is obtained:

The matrix formulation of the algebraic system of Eqs (3.4) and (3.5) can be represented in the following form:

The structure of each block of the matrices $\mathcal{P}$ and $\mathcal{Q}$ is given below:

For $r = 1, m$ the order of block $P_r$ and $Q_r$ is $8\times10$ due the homogeneous boundary constraints. The above components of the matrix blocks are defined as below:

The matrix formulation defined in Eq (3.6) is solved using mathematical software such as MATLAB. With the increment in the order of the algebraic system and decrement in the step size in the time direction, the computational cost increases. In order to balance the computational cost, the optimal choice of partitions has been chosen, and good accuracy is achieved.

4.   Stability analysis

In this section, the stability of the proposed technique is discussed using the Von-Neumann criterion of stability. For quasi-linearization, let us assume $\alpha(u)$ is locally bounded by $\eta$ over the given domain. Therefore, the linearized coupled system reads:

In the form of operator, the linearized problem in Eqs (4.1) and (4.2) can be written as

It is clear from Eq (4.1) that $(w+u_{\zeta\zeta})^{j+1}+(w+u_{\zeta\zeta})^j = 0$; also, initially $(w+u_{\zeta\zeta})^0 = 0$ this implies $(w+u_{\zeta\zeta})^j = 0 \forall j$. Thus, the boundedness of both $w$ and $u_{\zeta\zeta}$ throughout the given domain is ensured.

Application of the piecewise Hermite spline collocation technique as discussed in Eqs (3.2) and (3.3), in Eq (4.1), one gets:

Similarly, considering the relation between $w$ and $u_{\zeta\zeta}$, it can be derived that:

Again, considering the second equation of the coupled system stated in Eq (4.2), it is simple to say that:

Now, using the relations given in Eq (4.4), the following equation is obtained:

further rearranging the terms, Eq (4.5) can be written as:

Now, assume $a_j^{\gamma, n} = \rho^n\exp{(i\mu j)}$, where $\mu$ represents the mode number, $\rho$ is the amplification factor, and $i = \sqrt{-1}$, and substituting in Eq (4.6), one obtains:

where

It is clear from the structure of $M_p(\xi)$ and $N_p(\xi)$ that $\mid N_p(\xi)\mid \leq \mid M_p(\xi)\mid$. Therefore,

This implies that the given technique is unconditionally stable.

5.   Numerical illustrations

The efficiency of any technique is incomplete if it does not study the error norms. Define the $L_2$-norm and $L_{\infty}$-norm as:

$\Vert u\Vert_{\infty} = \max_{\zeta \epsilon \Omega}.\mid u(\zeta, t_j)\mid$,

$\Vert u\Vert_{2} = \sqrt{h\sum_{k = 0}^{m} u(\zeta_k, t_j)^2}$.

$\Vert e\Vert_{\infty} = \Vert u_{num} - u_{exact}\Vert_{\infty}$,

$\Vert e\Vert_{2} = \Vert u_{num}-u_{exact} \Vert_{2}$.

$GRE(t_j) = \frac{\sum_{k = 0}^m \mid u_{exact}(\zeta_k, t_j)- {u}_{num}(\zeta_k, t_j) \mid }{\sum_{k = 1}^m \mid u_{exact}(\zeta_k, t_j) \mid }$

where $e$ represents the difference between exact and approximate solution, $u_{exact}$ is the exact solution, and $u_{num}$ is the numerical approximation of the proposed problem with quintic Hermite splines. GRE represents the global relative error.

Example 1. In this test problem, Eq (1.2) having the compact form of solution as:

is considered, where $c_1 = \frac{1}{2}; k_1 = 6; \bar{\zeta} = -10$.

The initial and the boundary conditions to approximate the solution are derived from the equation (5.1). Different values of $\tau, h$ have been considered to analyze the problem. Figure 1 represents the graphical behavior of the soliton solution generated by the proposed technique. Figures 2 and 3 represent the exact solution and the absolute error showing the agreement of the approximate solution. Soliton solution behavior at different levels of time is represented in Figure 4, whereas Figures 5 and 6 depict the computed measures of error. The GRE and the $\Vert e\Vert_{2}$ both seem to decrease after $T = 3$, ensuring the accuracy of the solution even after a large number of iterations. Computationally, the error norms and GRE have been considered. In Table 2, the progress of GRE, $\Vert e\Vert_{2}$, and the $\Vert e\Vert_{\infty}$ with respect to time is studied. The effect of technique parameters like $\tau, m$ simultaneously on the approximation are represented in the form of GRE, $\Vert e\Vert_{2}$ in Table 3. Considering $\tau = 0.0001$, the behavior of GRE, $\Vert e\Vert_{2}$, and $\Vert e\Vert_{\infty}$ is presented in Table 4. A comparative study of the GRE generated from the technique with other methods is presented in Table 5. Clearly, the generated results are found to be better than the existing ones.

Example 2. In this test problem, the KdV-Burgers equation corresponding to $\varepsilon = 0$ and $\alpha(u) = u$ is considered in Eq (1.2), having the compact form of solution as:

where $\gamma = \frac{\kappa}{10\delta}; \rho = \frac{6\kappa^2}{25\delta}$.

For computational purpose, $\kappa = -9\times 10^{-4}, \delta = 2\times 10^{-5}$ are chosen. The initial and the boundary conditions to approximate the solution are derived from the Eq (5.2). The problem has been analyzed for different values of $\tau, h$. The shock wave profile of solution has been observed in Figures 7 and 8 for approximate and exact solution up to $T = 200$ with $\tau = 0.1, m = 50$. The absolute error between the approximate and exact solution is shown in Figure 9, which signifies that good accuracy has been achieved in approximation. The progress of shock wave with respect to different levels of time $T$ has been studied in Figure 10, whereas Figures 11 and 12 represents the growth of GRE and the $\Vert e\Vert_{2}$ with respect to time $T$. In Table 6, the computed behavior of GRE, $\Vert e\Vert_{2}$, and the $\Vert e\Vert_{\infty}$ with respect to time $T$, choosing $\tau = 0.1, m = 50$ is shown. Comparison of the computed GRE with those available in ^[10] and ^[44] has been studied in Table 7. The method proposed by ^[10] has studied the three values of GRE as shown in Table 7 corresponding to the different technique parameters. Clearly, the approximation is better with the proposed technique as compared to those available in ^[10] and ^[44].

Example 3. In this problem, the Gardner equation as a special case of Eq (1.2) with $\alpha(u) = u-5u^2$ and $\delta = 1, \varepsilon = 0$ is considered with the compact form of solution is defined in ^[45] as:

where $l_1 = \frac{\sqrt{30}}{60}$.

The initial and the boundary conditions are derived from the Eq (5.3). The problem has been studied for different values of $\tau$ and $h$. The kink-type behavior of the solution has been recorded. Figure 13 represents the exact solution, and Figure 14 represents the approximate solution computed with the proposed technique. The agreement of both exact and approximate solutions at the fully discrete temporal-spatial domain is plotted in Figure 15. As can be seen, the approximation accuracy is good. GRE up to $T = 15$ has been plotted in Figure 16. Table 8 compares the GRE and $\Vert e\Vert_{\infty}$ at different time levels with the results given in ^[45]. For $\tau = 0.001$ and $m = 50$ over the domain $[-50, 50]$, the behavior of GRE, $\Vert e\Vert_{2}$, and $\Vert e\Vert_{\infty}$ with respect to time is shown in Table 9. The effect of $\tau, m$ at $T = 1$ on the approximation is studied in Table 10 by GRE, $\Vert e\Vert_{2}$, and the $\Vert e\Vert_{\infty}$. It can be observed that the decrement in $\tau$ does not affect the accuracy of approximation very much for this problem, whereas the changes in $m$, i.e., decrement in $h$, improve the accuracy.

Example 4. In this test problem, Eq (1.2) is considered with the following initial function as :

The boundary conditions, in this case, are homogeneous. The problem has been studied for different values of dispersion and dissipation parameters. KdV-Burgers equation for the case $\varepsilon = 0$ has been analyzed for different values of $\kappa, \delta$. In Figure 17 to Figure 20, the solution behavior has been presented concerning the impact of dispersive and diffusive parameters. Further, the combined effect of dissipative, dispersive, and diffusive parameters has been presented from Figure 21 to Figure 26. It has been observed from Figure 21 to Figure 23 that for $\kappa = \delta = 1$, as the values of $\varepsilon$ decrease from 1 to 0.1, more waves arise with rising amplitudes. The simultaneous decrement in the values of dispersive and dissipative parameters has been presented from Figure 24 to Figure 26. The dispersive parameter increases the peaks of the traveling wave. $\Vert u\Vert_2$-norm has also been studied corresponding to dissipative, dispersive, and diffusive parameters. In Table 11, $\varepsilon = 0$, i.e., dissipation is zero, with decreasing values of dispersive parameter, a slight decrement in the $\Vert u\Vert_2$-norm has been observed, but in the simultaneous decrement in the diffusive and dispersive parameter, the $\Vert u\Vert_2$-norm increases as for $\kappa = 0.1, \delta = 0.01$. In Table 12, the effects of dissipation and dispersion parameters are observed. It is observed that with decreasing $\varepsilon$ as well as $\delta$, the $\Vert u \Vert_2$-norm decreases. The effect of decreasing $\varepsilon$ is greater as compared to decreasing $\delta$ as for $\kappa = 1, \delta = 1, \varepsilon = 0.1$ and $\kappa = 1, \delta = 0.1, \varepsilon = 0.1$.

Example 5. In this test problem, the KdV equation is considered as a special case of Eq (1.2) with the initial function as:

The problem has been analyzed for the long-time behavior. Also, the effect of domain range has been considered. For time direction, the behavior of the solution from a small time range as in Figure 27 to long-time behavior of the traveling wave as in Figure 30 is examined. For spatial direction, the domain has been extended from $(-10, 10)$ to $(-60, 60)$ and the effect is noticed correspondingly. Different values of $\tau$ and $m$ have been considered to analyze the solution for the long-time domain as well as the wider spatial domain. The decreasing wave amplitude has been presented in Figures 28, 29, and 30 for the time and spatial domain. The computed values of $\Vert u\Vert_2$-norm are presented in Table 13 and Table 14. The effect of the domain and $\tau$ on the $\Vert u\Vert_2$-norm at large time levels is studied in Table 13 whereas Table 14 represents the $\Vert u\Vert_2$-norm at small time. At the higher times, the $\Vert u\Vert_2$-norm seems to decrease irrespective of the domain.

6.   Conclusions

In this study, higher-order nonlinear dispersive dissipative partial differential equations have been studied. The computational study has been conducted using a hybrid technique comprising a weighted finite difference scheme followed by a collocation method using quintic Hermite splines. The main focus of the present study is toward the traveling wave solution exhibited by the dispersive-dissipative systems. The simultaneous impact of the dispersive-dissipative process, dominating the dispersive process over the dissipative process and vice-versa, has been studied thoroughly.

It is observed that the dominance of dispersion raises the nonlinear waves, whereas the dominance of the dissipative process raises the amplitude of nonlinear waves without dispersion, as studied in Figure 21 to Figure 26. Also, it has been studied that as the dispersion parameter, i.e., $\delta$ decreases, the nonlinear waves disappear and make the amplitude higher with the sharp peak, as studied from Figure 24 to Figure 26.

Soliton solution, shock wave solution, and Kink-type solution have also been studied with good approximation accuracy. The stability analysis and the convergence analysis of the adopted technique have also been studied. The computed results show that the applied technique gives fruitful results in studying the traveling wave behavior of nonlinear systems. A comparative study of the computed results with those available in the literature presents the accuracy and efficiency of the technique as well. A good accuracy has been achieved for the small number of partitions in the space as well as in the time direction.

Author contributions

Conceptualization, Shelly Arora and Abdul-Majeed Ayebire; Methodology, Priyanka; Software, Priyanka and Saroj Sahani.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Dr. Shelly Arora is thankful to SERB-POWER to provide grant to complete this research via grant number SPG/2022/001269. Mr. Abdul Majeed Ayebire is thankful to ICCR for providing grant via reference number LY7930408566174

Conflict of interest

Authors declare that they have no conflict of interest.