Exploring the fourth-order Boussinesq water wave equation: soliton analysis, modulation instability, sensitivity behavior, and chaotic analysis

In this article, we reveal the novel types of exact solitons to the fourth-order nonlinear (1 + 1)-dimensional Boussinesq water wave equation. This model is obtained under the consideration of the smaller water depth and larger wavelength of the waves. The Boussinesq water wave equation is useful in understanding water wave behavior, harbor design, coastal dynamics, wave propagation in shallow seas, ocean wave models, marine environments, etc. For our aim, we used the Sardar sub-equation technique. As a result, new types of exact wave solitons involving trigonometry, hyperbolic trigonometry, and rational functions are gained. Some gained solutions are represented through 2D, 3D, contour, and density plots. In bifurcation analysis, a new planar dynamical system of the governing model is obtained by applying the Galilean transformation, and all possible phase portraits are discussed. Modulation instability is used to obtain the steady-state solutions of the concerned model. Furthermore, the chaotic behavior of the governing model is analyzed. Sensitivity analysis is utilized to determine the sensitivity behavior of the model. The achieved solutions are fruitful in distinct areas of mathematical physics and engineering fields. At the end, the technique is a useful and reliable approach to solving other important nonlinear partial differential equations. This study applies the Sardar sub-equation method to derive new analytical solutions of the fourth-order nonlinear (1 + 1)-dimensional Boussinesq water wave equation. The method demonstrates greater flexibility than traditional approaches in handling nonlinear terms. However, the results depend on specific parameter conditions, and experimental or numerical validation is left for future investigation.

1 Introduction

Nonlinear partial differential equations (PDEs) are the mathematical form of naturally occurring phenomena. In different fields of science and engineering, there are various PDEs, including the Akbota–Gudekli–Kairat–Zhaidary equation [1], the Kodama equation [2], the extended Kairat-II equation [3], the complex-coupled Kuralay system [4], and the Chaffee–Infante equation [5]. Different techniques have been developed to obtain the various kinds of exact soliton solutions of the nonlinear PDEs, such as the Kumar–Malik technique [6], the modified sub-equation technique [7], and the multivariate generalized exponential rational integral function technique [8]. Water wave equations are utilized to explain the various types of water waves, including sinusoidal waves, nonlinear wave interaction, and shallow water waves. Water wave equations have many applications in different fields, including fluid dynamics, oceanography, and engineering.

The water wave interpolated Boussinesq equation was introduced in 1871 [9] and is given as

$g_{tt}-a{\left(g^2\right)}_{xx}-b{g}_{\mathit{xxxx}}-g_{xx}=0.(1)$

Equation 1 is a standard Boussinesq equation that explains the shallow water wave interaction process solution. This equation includes the various waves and shallow water effects in fluid dynamics, like shoaling, refraction, and weak nonlinearity.

Consider the fourth-order nonlinear Boussinesq water wave equation given in [10]

$g_{tt}-a{\left(g^2\right)}_{xx}-b{g}_{\mathit{xxxx}}+c{g}_{xt}-g_{xx}=0,(2)$

where $g=g(x,t)$ is a wave function, and a, b, and c are the constants. Constant a controls the nonlinearity strength, b is a dispersion coefficient (the term provides high-order dispersion that stabilizes wave steepening), and c represents the damping or mixed effects. The balance between nonlinearity a and dispersion b gives rise to periodic waves. These constants have physical meanings related to buoyancy, pressure, or roughness. Their values can be assumed to obtain specific solutions, such as solitons or other wave solutions. If $b=0$, Equation 1 becomes nonlinear and possibly unstable. If $a=0$, Equation 1 becomes linear, supporting only dispersive linear waves. If $c=0$, Equation 1 reduces to a classical fourth-order Boussinesq equation.

The fourth-order nonlinear Boussinesq water wave equation is of great importance. This equation models the behavior of water waves in shallow water, making it relevant for coastal engineering, oceanography, and tsunami research. The equation’s nonlinearity captures complex wave interactions, leading to fascinating phenomena like wave breaking, soliton formation, and chaotic behavior. It has applications in various fields, including fluid dynamics, coastal engineering, and plasma physics. The equation’s fourth-order nature and nonlinearity make it a rich source of mathematical challenges and opportunities for developing new analytical and numerical methods. Understanding the behavior of water waves is crucial for predicting coastal erosion, flooding, and damage to offshore structures.

This model is obtained under the consideration of the smaller water depth and larger wavelength of the waves. The Boussinesq water wave equation is useful in wave behavior, harbor design, wave propagation in shallow seas, etc. Equation 2 was developed by Wazwaz and Kaur in [18]. In the literature, different solutions of Equation 2 are obtained by using the distinct schemes, including the F-expansion scheme [10], the auxiliary equation scheme [11], the $\exp (-\varphi (\eta ))$-expansion scheme [12], the Jacobi elliptic function expansion scheme [13], the $(G^{\prime }/G)$-expansion scheme [14], the exponential expansion scheme [15], the generalized Arnous method [16], and the physics-informed neural networks technique [17].

We used a simple and useful technique, the Sardar sub-equation (SSE) technique. This technique is applied to achieve various types of exact wave results using the Sawada–Kotera equation [19]; exact solitons of the Fokas–Lenells equation are achieved [20]; some exact wave solitons, including dark, bright, periodic-singular, singular, and dark-bright soliton solutions, are gained for the Zakharov equation [21]; different kinds of optical wave solitons, having periodic wave, dark, bright, and singular solitons, are achieved for the stochastic Schrödinger wave model [22]; singular, kink, and periodic solitons are obtained for the Boiti–Leon–Manna–Pempinelli model [23]; optical solitons, having dark, bright, periodic, and kink, are obtained for the Biswas–Milovic model [24]; and kink, bright, dark, and periodic solitons for the coupled Drinfel’d–Sokolov–Wilson equation are achieved [25].

The fundamental purpose of our work is to explore the distinct exact wave solutions of a (1 + 1)-dimensional Boussinesq water wave equation by utilizing the Sardar sub-equation method. Different analyses, including the modulation instability, bifurcation analysis, chaotic behavior, sensitivity nature, and the Lyapunov exponent of the concerned equation, are performed.

The motivation of this work is to investigate the novel kinds of exact solitons for the fourth-order nonlinear Boussinesq water wave equation by using the Sardar sub-equation technique. For the fourth-order Boussinesq water wave equation, the Sardar sub-equation technique is used for the first time. The obtained solutions do not currently exist in the literature. Some of the dynamical analyses, including modulation instability, bifurcation analysis, chaotic behavior, sensitivity analysis, and Lyapunov exponent analysis, are discussed for the fourth-order nonlinear Boussinesq water wave equation for the first time in the literature.

The article consists of the following sections: The technique is explained in Section 2; the mathematical analysis and exact wave results are mentioned in Section 3; a graphical interpretation is given in Section 4; bifurcation analysis is done in Section 5; chaotic behavior is demonstrated in Section 6; Lyapunov exponent analysis is performed in Section 7; sensitivity nature is discussed in Section 8; modulation instability analysis is performed in Section 9; results and discussion are given in Section 10; and conclusion is provided in Section 11.

2 Methodology

Now, we will represent the Sardar sub-equation method [26] by assuming the nonlinear PDE:

$J\left(g,g_x,g_t,g_{xx},g_{xt},g{g}_{xt},g_{\mathit{xxt}},\dots \right)=0.(3)$

Here, $g=g(x,t)$ represents the function. Putting the given transformations, we get

$g=G\left(\mathrm{\Omega }\right),\mathrm{\Omega }=\lambda x+\mu t.(4)$

The results are given in the form of a nonlinear ordinary differential equation (NLODE):

$Y\left(G,\lambda G^{\prime },\mu G^{\prime },{\lambda }^2{G}^{″},\lambda \mu G^{″},\dots \right)=0.(5)$

Assuming the result of Equation 5 is given as

$G\left(\mathrm{\Omega }\right)={\displaystyle \sum _{i=0}^m}{b}_i{\psi }^i\left(\mathrm{\Omega }\right).(6)$

Here, $\psi (\mathrm{\Omega })$ fulfills the following equation:

${\psi }^{\prime }\left(\mathrm{\Omega }\right)=\sqrt{\sigma +\kappa {\psi }^2\left(\mathrm{\Omega }\right)+{\psi }^4\left(\mathrm{\Omega }\right)}.(7)$

Here, $\sigma$ and $\kappa$ are parameters.

Putting Equations 6, 7 into Equation 5 and summing up the coefficients of every ${\psi }^i$ term, taking each equal to zero, to gain a set. Simplifying a set, we gain values of undetermined. The solutions of Equation 7 for the different conditions of $\sigma$ and $\kappa$ are given in [27].

Motivation of the method:

This method can effectively handle the nonlinearity of equations, providing solutions in terms of generalized trigonometric and hyperbolic functions. This method can generate different kinds of solutions, including dark, bright, singular, periodic-singular, combined dark-bright, and dark-singular. This technique is considered simple and reliable for solving nonlinear evaluation equations. This method can be applied to various physical systems, including optical fibers, fluid dynamics, and plasma physics, making it a valuable tool for understanding complex phenomena.

2.1 Limitations

The Sardar sub-equation method relies on specific parameter conditions to obtain exact solutions, which might not apply to all cases. The method’s effectiveness is often demonstrated through mathematical derivations and numerical simulations, but experimental validation is necessary to confirm the accuracy of the results. The method might not be applicable to all types of nonlinear partial differential equations (NLPDEs) or systems with complex nonlinearities. The Sardar sub-equation method can become computationally intensive or even intractable for high-dimensional problems. The method might not guarantee finding all possible solutions to the NLPDE, and other methods might be needed to find additional solutions.

3 Mathematical analysis

Consider the given wave transformation:

$g\left(x,t\right)=G\left(\mathrm{\Omega }\right),\mathrm{\Omega }=\left(\mu x-\omega t\right).(8)$

By using Equation 8 in Equation 2, we obtain

$b{\mu }^4{G}^{\left(4\right)}+G^{″}\left(c\mu \omega +{\mu }^2-{\omega }^2\right)+2a{\mu }^{2}G{G}^{″}+2{\mu }^2{\left(G^{\prime }\right)}^2=0.(9)$

By integrating twice and assuming integration constants equal to zero, we get

$a{\mu }^2{G}^2+b{\mu }^4{G}^{″}+\left(c\mu \omega +{\mu }^2-{\omega }^2\right)G=0.(10)$

By using the homogenous balance technique and balancing the terms $G^{″}$ and $G^2$, we achieve $m=2$. Now, we will find the exact wave solutions using the Sardar sub-equation method.

3.1 Exact solitons

In our case, Equation 6 changes into

$G\left(\mathrm{\Omega }\right)=b_0+b_1\psi \left(\mathrm{\Omega }\right)+b_2{\psi }^2\left(\mathrm{\Omega }\right).(11)$

By using Equation 11 in Equation 10 along with Equation 7, we gain solution sets:

Solution set 1:

$\left\{b_0=-\frac{2b{\mu }^2\left(\sqrt{{\kappa }^2-3\sigma }+\kappa \right)}{a},b_1=0,b_2=-\frac{6b{\mu }^2}{a},\omega =\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-16b{\mu }^2\sqrt{{\kappa }^2-3\sigma }}\right)\right\}.(12)$

By using Equations 8, 11, 12, and solutions mentioned in [27], we get the following solutions:

$g\left(x,t\right)=-\frac{2b{\mu }^2}{a}\left(2\kappa +3{\left(\sqrt{\kappa rs}{\text{sech}}_{rs}\left(\sqrt{\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-16b{\mu }^2\kappa }\right)t\right)\right)\right)}^2\right),(13)$

$g\left(x,t\right)=-\frac{2b{\mu }^2}{a}\left(2\kappa +3{\left(\sqrt{\kappa rs}{\text{csch}}_{rs}\left(\sqrt{\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-16b{\mu }^2\kappa }\right)t\right)\right)\right)}^2\right),(14)$

$g\left(x,t\right)=-\frac{2b{\mu }^2}{a}\left(2\kappa +3{\left(\sqrt{-\kappa rs}{\sec }_{rs}\left(\sqrt{-\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-16b{\mu }^2\kappa }\right)t\right)\right)\right)}^2\right),(15)$

$g\left(x,t\right)=-\frac{2b{\mu }^2}{a}\left(2\kappa +3{\left(\sqrt{-\kappa rs}{\csc }_{rs}\left(\sqrt{-\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-16b{\mu }^2\kappa }\right)t\right)\right)\right)}^2\right),(16)$

$g\left(x,t\right)=-\frac{3b{\mu }^2}{a}\left(\kappa +2{\left(\sqrt{-\frac{\kappa }{2}}{\tanh }_{rs}\left(\sqrt{-\frac{\kappa }{2}}\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-8b{\mu }^2\kappa }\right)t\right)\right)\right)}^2\right),(17)$

$g\left(x,t\right)=-\frac{3b{\mu }^2}{a}\left(\kappa +2{\left(\sqrt{-\frac{\kappa }{2}}{\coth }_{rs}\left(\sqrt{-\frac{\kappa }{2}}\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-8b{\mu }^2\kappa }\right)t\right)\right)\right)}^2\right).(18)$

$\begin{array}{ll}\hfill g\left(x,t\right)& =-\frac{3b{\mu }^2}{a}\left(\kappa +2\left(\sqrt{-\frac{\kappa }{2}}\left({\tanh }_{rs}\left(\sqrt{-2\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-8b{\mu }^2\kappa }\right)t\right)\right)\right.\right.\right.\hfill \\ \hfill & \left.{\left.\left.\pm \iota \sqrt{rs}sec{h}_{rs}\left(\sqrt{-2\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-8b{\mu }^2\kappa }\right)t\right)\right)\right)\right)}^2\right),\hfill \end{array}(19)$

$\begin{array}{ll}\hfill g\left(x,t\right)& =-\frac{3b{\mu }^2}{a}\left(\kappa +2\left(\sqrt{-\frac{\kappa }{2}}\left({\coth }_{rs}\left(\sqrt{-2\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-8b{\mu }^2\kappa }\right)t\right)\right)\right.\right.\right.\hfill \\ \hfill & \left.{\left.\left.\pm \sqrt{rs}csc{h}_{rs}\left(\sqrt{-2\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-8b{\mu }^2\kappa }\right)t\right)\right)\right)\right)}^2\right),\hfill \end{array}(20)$

$\begin{array}{ll}\hfill g\left(x,t\right)& =-\frac{3b{\mu }^2}{a}\left(\kappa +2\left(\sqrt{-\frac{\kappa }{8}}\left({\tanh }_{rs}\left(\sqrt{-\frac{\kappa }{8}}\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-8b{\mu }^2\kappa }\right)t\right)\right)\right.\right.\right.\hfill \\ \hfill & \left.{\left.\left.+{\coth }_{rs}\left(\sqrt{-\frac{\kappa }{8}}\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-8b{\mu }^2\kappa }\right)t\right)\right)\right)\right)}^2\right),\hfill \end{array}(21)$

$g\left(x,t\right)=-\frac{3b{\mu }^2}{a}\left(\kappa +2{\left(\sqrt{\frac{\kappa }{2}}{\tan }_{rs}\left(\sqrt{\frac{\kappa }{2}}\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-8b{\mu }^2\kappa }\right)t\right)\right)\right)}^2\right),(22)$

$g\left(x,t\right)=-\frac{3b{\mu }^2}{a}\left(\kappa +2{\left(\sqrt{\frac{\kappa }{2}}{\cot }_{rs}\left(\sqrt{\frac{\kappa }{2}}\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-8b{\mu }^2\kappa }\right)t\right)\right)\right)}^2\right),(23)$

$\begin{array}{ll}\hfill g\left(x,t\right)& =-\frac{3b{\mu }^2}{a}\left(\kappa +2\left(\sqrt{\frac{\kappa }{2}}\left({\tan }_{rs}\left(\sqrt{2\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-8b{\mu }^2\kappa }\right)t\right)\right)\right.\right.\right.\hfill \\ \hfill & \left.{\left.\left.\pm \sqrt{rs}{\sec }_{rs}\left(\sqrt{2\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-8b{\mu }^2\kappa }\right)t\right)\right)\right)\right)}^2\right),\hfill \end{array}(24)$

$\begin{array}{ll}\hfill g\left(x,t\right)& =-\frac{3b{\mu }^2}{a}\left(\kappa +2\left(\sqrt{\frac{\kappa }{2}}\left({\cot }_{rs}\left(\sqrt{2\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-8b{\mu }^2\kappa }\right)t\right)\right)\right.\right.\right.\hfill \\ \hfill & \left.{\left.\left.\pm \sqrt{rs}{\csc }_{rs}\left(\sqrt{2\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-8b{\mu }^2\kappa }\right)t\right)\right)\right)\right)}^2\right),\hfill \end{array}(25)$

$\begin{array}{ll}\hfill g\left(x,t\right)& =-\frac{3b{\mu }^2}{a}\left(\kappa +2\left(\sqrt{\frac{\kappa }{8}}\left({\tan }_{rs}\left(\sqrt{\frac{\kappa }{8}}\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-8b{\mu }^2\kappa }\right)t\right)\right)\right.\right.\right.\hfill \\ \hfill & \left.{\left.\left.+{\cot }_{rs}\left(\sqrt{\frac{\kappa }{8}}\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{c^2+4-8b{\mu }^2\kappa }\right)t\right)\right)\right)\right)}^2\right).\hfill \end{array}(26)$

Solution set 2:

$\left\{b_0=-\frac{2b{\mu }^2\left(\kappa -\sqrt{{\kappa }^2-3\sigma }\right)}{a},b_1=0,b_2=-\frac{6b{\mu }^2}{a},\omega =\frac{\mu }{2}\left(c\pm \sqrt{16b{\mu }^2\sqrt{{\kappa }^2-3\sigma }+c^2+4}\right)\right\}.(27)$

By using Equations 8, 11, 27 and the solutions mentioned in [27], we get the following solutions:

$g\left(x,t\right)=-\frac{6b{\mu }^2}{a}{\left(\sqrt{\kappa rs}{\text{sech}}_{rs}\left(\sqrt{\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{16b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right)}^2,(28)$

$g\left(x,t\right)=-\frac{6b{\mu }^2}{a}{\left(\sqrt{\kappa rs}{\text{csch}}_{rs}\left(\sqrt{\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{16b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right)}^2,(29)$

$g\left(x,t\right)=-\frac{6b{\mu }^2}{a}{\left(\sqrt{-\kappa rs}{\sec }_{rs}\left(\sqrt{-\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{16b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right)}^2,(30)$

$g\left(x,t\right)=-\frac{6b{\mu }^2}{a}{\left(\sqrt{-\kappa rs}{\csc }_{rs}\left(\sqrt{-\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{16b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right)}^2,(31)$

$g\left(x,t\right)=-\frac{b{\mu }^2\kappa }{a}\left(1+6{\left(\sqrt{-\frac{\kappa }{2}}{\tanh }_{rs}\left(\sqrt{-\frac{\kappa }{2}}\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{8b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right)}^2\right),(32)$

$g\left(x,t\right)=-\frac{b{\mu }^2\kappa }{a}\left(1+6{\left(\sqrt{-\frac{\kappa }{2}}{\coth }_{rs}\left(\sqrt{-\frac{\kappa }{2}}\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{8b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right)}^2\right),(33)$

$\begin{array}{ll}\hfill g\left(x,t\right)& =-\frac{b{\mu }^2\kappa }{a}\left(1+6\left(\sqrt{-\frac{\kappa }{2}}\left({\tanh }_{rs}\left(\sqrt{-2\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{8b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right.\right.\right.\hfill \\ \hfill & \left.{\left.\left.\pm \iota \sqrt{rs}sec{h}_{rs}\left(\sqrt{-2\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{8b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right)\right)}^2\right),\hfill \end{array}(34)$

$\begin{array}{ll}\hfill g\left(x,t\right)& =-\frac{b{\mu }^2\kappa }{a}\left(1+6\left(\sqrt{-\frac{\kappa }{2}}\left({\coth }_{rs}\left(\sqrt{-2\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{8b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right.\right.\right.\hfill \\ \hfill & \left.{\left.\left.\pm \sqrt{rs}csc{h}_{rs}\left(\sqrt{-2\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{8b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right)\right)}^2\right),\hfill \end{array}(35)$

$\begin{array}{ll}\hfill g\left(x,t\right)& =-\frac{b{\mu }^2\kappa }{a}\left(1+6\left(\sqrt{-\frac{\kappa }{8}}\left({\tanh }_{rs}\left(\sqrt{-\frac{\kappa }{8}}\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{8b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right.\right.\right.\hfill \\ \hfill & \left.{\left.\left.+{\coth }_{rs}\left(\sqrt{-\frac{\kappa }{8}}\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{8b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right)\right)}^2\right),\hfill \end{array}(36)$

$g\left(x,t\right)=-\frac{b{\mu }^2\kappa }{a}\left(1+6{\left(\sqrt{\frac{\kappa }{2}}{\tan }_{rs}\left(\sqrt{\frac{\kappa }{2}}\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{8b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right)}^2\right),(37)$

$g\left(x,t\right)=-\frac{b{\mu }^2\kappa }{a}\left(1+6{\left(\sqrt{\frac{\kappa }{2}}{\cot }_{rs}\left(\sqrt{\frac{\kappa }{2}}\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{8b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right)}^2\right),(38)$

$\begin{array}{ll}\hfill g\left(x,t\right)& =-\frac{b{\mu }^2\kappa }{a}\left(1+6\left(\sqrt{\frac{\kappa }{2}}\left({\tan }_{rs}\left(\sqrt{2\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{8b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right.\right.\right.\hfill \\ \hfill & \left.{\left.\left.\pm \sqrt{rs}{\sec }_{rs}\left(\sqrt{2\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{8b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right)\right)}^2\right),\hfill \end{array}(39)$

$\begin{array}{ll}\hfill g\left(x,t\right)& =-\frac{b{\mu }^2\kappa }{a}\left(1+6\left(\sqrt{\frac{\kappa }{2}}\left({\cot }_{rs}\left(\sqrt{2\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{8b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right.\right.\right.\hfill \\ \hfill & \left.{\left.\left.\pm \sqrt{rs}{\csc }_{rs}\left(\sqrt{2\kappa }\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{8b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right)\right)}^2\right),\hfill \end{array}(40)$

$\begin{array}{ll}\hfill g\left(x,t\right)& =-\frac{b{\mu }^2\kappa }{a}\left(1+6\left(\sqrt{\frac{\kappa }{8}}\left({\tan }_{rs}\left(\sqrt{\frac{\kappa }{8}}\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{8b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right.\right.\right.\hfill \\ \hfill & \left.{\left.\left.+{\cot }_{rs}\left(\sqrt{\frac{\kappa }{8}}\left(\mu x-\frac{\mu }{2}\left(c\pm \sqrt{8b{\mu }^2\kappa +c^2+4}\right)t\right)\right)\right)\right)}^2\right).\hfill \end{array}(41)$

4 Graphical interpretation

In this section, the results gained are demonstrated through 2- and 3-dimensional, contour, and density figures with the use of Mathematica software.

Figure 1a demonstrates the 2D graph of a bright soliton when $-10<x<10$ at different values of t. The purple color represents the graph at $t=0$, the yellow color represents the graph at $t=1$, and the pink color represents the graph at $t=2$. We can observe that the wave solution is time-dependent because the phase of the wave solution has shifted with time. Figure 1b demonstrates the 3D plot when $-10<x<10$ and $0<t<10$. This shows that the wave solution has a symmetric property because the wave solution has a sharp transition in both the x direction and the t direction. Figure 1c indicates a contour graph when $-10<x<10$ and $0<t<10$. Figure 1d indicates a density graph when $-10<x<10$ and $0<t<10$.

Figure 1

Figure 1. (Bright soliton) Graphical representation of $|g(x,t)|$ appearance in Equation 13 for $a=1,\hspace{0.17em}\hspace{0.17em}b=1,\hspace{0.17em}\hspace{0.17em}c=1,\hspace{0.17em}\hspace{0.17em}\mu =1,\hspace{0.17em}\hspace{0.17em}r=1,\hspace{0.17em}\hspace{0.17em}s=1$, and $\kappa =0.1$. (a) 2D plot. (b) 3D surface plot. (c) Contour plot. (d) Density plot.

Figure 2a demonstrates the 2-D graph of a singular soliton when $-5<x<5$ at different values of t. The purple color represents the graph at $t=0$, the yellow color represents the graph at $t=1$, and the pink color represents the graph at $t=2$. We can observe that the wave solution is time-dependent because the phase of the wave solution has shifted with time. Figure 2b demonstrates the 3D plot when $-5<x<5$ and $0<t<10$. This shows that the wave solution has a symmetric property because the wave solution has a sharp transition in both the x direction and the t direction. Figure 2c indicates a contour graph when $-5<x<5$ and $0<t<10$. Figure 2d indicates a density graph when $-5<x<5$ and $0<t<10$.

Figure 2

Figure 2. (Singular soliton) Graphical representation of $|g(x,t)|$ appearance in Equation 14 for $a=-1,\hspace{0.17em}\hspace{0.17em}b=1,\hspace{0.17em}\hspace{0.17em}c=1,\hspace{0.17em}\hspace{0.17em}\mu =1,\hspace{0.17em}\hspace{0.17em}r=1,\hspace{0.17em}\hspace{0.17em}s=1$, and $\kappa =0.3$. (a) 2D plot. (b) 3D surface plot. (c) Contour plot. (d) Density plot.

Figure 3a demonstrates the 2-D graph of a periodic soliton when $-12<x<12$ at different values of t. The purple color represents the graph at $t=0$, the yellow color represents the graph at $t=1$, and the pink color represents the graph at $t=2$. We can observe that the wave solution is time-dependent because the phase of the wave solution has shifted with time. Figure 3b demonstrates the 3D plot when $-12<x<12$ and $0<t<10$. This shows that the wave solution has a symmetric property because the wave solution has pole-like waves in both the x direction and the t direction. Figure 3c indicates a contour graph when $-10<x<10$ and $0<t<10$. Figure 3d indicates a density graph when $-10<x<10$ and $0<t<10$.

Figure 3

Figure 3. (Periodic wave solution) Graphical representation of $|g(x,t)|$ appearance in Equation 15 for $a=-1,\hspace{0.17em}\hspace{0.17em}b=1,\hspace{0.17em}\hspace{0.17em}c=1,\hspace{0.17em}\hspace{0.17em}\mu =1,\hspace{0.17em}\hspace{0.17em}r=1,\hspace{0.17em}\hspace{0.17em}s=1$, and $\kappa =-0.1$. (a) 2D plot. (b) 3D surface plot. (c) Contour plot. (d) Density plot.

Figure 4a demonstrates the 2D graph of a dark soliton when $-12<x<12$ at different values of t. The purple color represents the graph at $t=0$, the yellow color represents the graph at $t=1$, and the pink color represents the graph at $t=2$. We can observe that the wave solution is time-dependent because the phase of the wave solution has shifted with time. Figure 4b demonstrates the 3D plot when $-10<x<10$ and $0<t<10$. This shows that the wave solution has a symmetric property because the wave solution has pole-like waves in both the x direction and the t direction. Figure 4c indicates a contour graph when $-10<x<10$ and $0<t<10$. Figure 4d indicates a density graph when $-10<x<10$ and $0<t<10$.

Figure 4

Figure 4. (Cuspon soliton) Graphical representation of $|g(x,t)|$ appearance in Equation 16 for $a=1,\hspace{0.17em}\hspace{0.17em}b=1,\hspace{0.17em}\hspace{0.17em}c=1,\hspace{0.17em}\hspace{0.17em}\mu =1,\hspace{0.17em}\hspace{0.17em}r=1,\hspace{0.17em}\hspace{0.17em}s=1$, and $\kappa =-1$. (a) 2-D plot. (b) 3-D surface plot. (c) Contour plot. (d) Density plot.

Figure 5a demonstrates the 2-D graph of a kink-like soliton when $-12<x<12$ at different values of t. The purple color represents the graph at $t=0$, the yellow color represents the graph at $t=1$, and the pink color represents the graph at $t=2$. We can observe that the wave solution is time-dependent because the phase of the wave solution has shifted with time. Figure 5b demonstrates the 3D plot when $-10<x<10$ and $0<t<10$. This shows that the wave solution has a symmetric property because the wave solution has a sharp transition in both the x direction and the t direction. Figure 5c indicates a contour graph when $-10<x<10$ and $0<t<10$. Figure 5d indicates a density graph when $-10<x<10$ and $0<t<10$.

Figure 5

Figure 5. (Dark soliton) Graphical representation of $|g(x,t)|$ appearance in Equation 17 for $a=1,\hspace{0.17em}\hspace{0.17em}b=1,\hspace{0.17em}\hspace{0.17em}c=1,\hspace{0.17em}\hspace{0.17em}\mu =1,\hspace{0.17em}\hspace{0.17em}r=1,\hspace{0.17em}\hspace{0.17em}s=1$, and $\kappa =-0.08$. (a) 2D plot. (b) 3D surface plot. (c) Contour plot. (d) Density plot.

Figure 6a demonstrates the 2D graph of a complex soliton when $-12<x<12$ at different values of t. The purple color represents the graph at $t=0$, the yellow color represents the graph at $t=1$, and the pink color represents the graph at $t=2$. We can observe that the wave solution is time-dependent because the phase of the wave solution has shifted with time. Figure 6b demonstrates the 3D plot when $-10<x<10$ and $0<t<10$. This shows that the wave solution has a symmetric property because the wave solution has a sharp transition in both the x direction and the t direction. Figure 6c indicates a contour graph when $-10<x<10$ and $0<t<10$. Figure 6d indicates a density graph when $-10<x<10$ and $0<t<10$.

Figure 6

Figure 6. (Dark-bright soliton) Graphical representation of $g(x,t)$ appearance in Equation 19 for $a=1,\hspace{0.17em}\hspace{0.17em}b=1,\hspace{0.17em}\hspace{0.17em}c=1,\hspace{0.17em}\hspace{0.17em}\mu =1,\hspace{0.17em}\hspace{0.17em}r=1,\hspace{0.17em}\hspace{0.17em}s=1$, and $\kappa =-0.08$. (a) 2D plot. (b) 3D surface plot. (c) Contour plot. (d) Density plot.

Figure 7a demonstrates the 2-D graph of a bright soliton when $-12<x<12$ at different values of t. The purple color represents the graph at $t=0$, the yellow color represents the graph at $t=1$, and the pink color represents the graph at $t=2$. We can observe that the wave solution is time-dependent because the phase of the wave solution has shifted with time. Figure 7b demonstrates the 3D plot when $-10<x<10$ and $0<t<10$. This shows that the wave solution has a symmetric property because the wave solution has a sharp transition in both the x direction and the t direction. Figure 7c indicates a contour graph when $-10<x<10$ and $0<t<10$. Figure 7d indicates a density graph when $-10<x<10$ and $0<t<10$.

Figure 7

Figure 7. (Bright soliton) Graphical representation of $g(x,t)$ appearance in Equation 28 for $a=-2,\hspace{0.17em}\hspace{0.17em}b=1,\hspace{0.17em}\hspace{0.17em}c=1,\hspace{0.17em}\hspace{0.17em}\mu =1,\hspace{0.17em}\hspace{0.17em}r=1,\hspace{0.17em}\hspace{0.17em}s=1$, and $\kappa =0.1$. (a) 2D plot. (b) 3D surface plot. (c) Contour plot. (d) Density plot.

5 Bifurcation analysis

The idea of bifurcation denotes the mathematical changes in a system, as well as the quality of the results gained by a system of differential equations. This analysis is common in research into mathematical models of dynamical systems. Bifurcation phenomena take place when a small change in parametric values leads to a sudden change in behavior. This concept may be used for a problem containing a split quality. This analysis delves into standard models like stability and into the composition of dividing solutions briefly.

Here, we will give a new planar dynamical system obtained by Equation 10 by using a Galilean transformation. By utilizing a Galilean transformation in Equation 10, we get

$\begin{array}{ll}\hfill \frac{dG\left(\mathrm{\Omega }\right)}{d\mathrm{\Omega }}& =F\left(\mathrm{\Omega }\right),\hfill \\ \hfill \frac{dF\left(\mathrm{\Omega }\right)}{d\mathrm{\Omega }}& =-A_1{G}^2\left(\mathrm{\Omega }\right)+A_{2}G\left(\mathrm{\Omega }\right).\hfill \end{array}(42)$

Here, $A_1=\frac{a}{b{\mu }^2}$ and $A_2=-\frac{c\mu \omega +{\mu }^2-{\omega }^2}{b{\mu }^4}$.

Bifurcation analysis includes the phase portraits of the governing system shown in Equation 42. First, one obtains a Hamiltonian function for the governing system in Equation 42, which is given as follows:

$H\left(G,F\right)=\frac{F^2}{2}+\frac{A_1{G}^3}{3}-\frac{A_2{G}^2}{2}=h.(43)$

Here, h represents the Hamiltonian constant.

For the purpose of obtaining the equilibrium points, we assume a new system given as

$\begin{array}{l}\hfill 0=F\left(\mathrm{\Omega }\right),\\ \hfill 0=-A_1{G}^2\left(\mathrm{\Omega }\right)+A_{2}G\left(\mathrm{\Omega }\right).\end{array}(44)$

By solving the above system, we get the equilibrium points given as

$E_1=(\mathrm{0,0})$ and $E_2=(\frac{A_2}{A_1},0)$.

The determinant of the Jacobian matrix of the system given by Equation 44 is

$D(G,Y)=\left|\begin{array}{cc}\hfill 0\hfill & \hfill 1\hfill \\ \hfill -A_1{G}^2+A_{2}G\hfill & \hfill 0\hfill \end{array}\right|=A_1{G}^2-A_{2}G$.

According to [30], we get the following classification conditions for the equilibrium points:

1-When $D(G,0)<0$, the (G, 0) point is called a saddle.

2-When $D(G,0)>0$, the (G, 0) point is called a center.

3-When $D(G,0)=0$, the (G, 0) point is called cuspidal.

4-When $D(G,0)>0$ and ${(\tau (D(G,Y)))}^2-4D(G,Y)>0$, the point (G, Y) is called a node, where $\tau$ denotes the trace of the system as shown by Equation 44.

Case 1: if $A_1>0$ and $A_2>0$.

By using the parametric values $a=1,b=1,c=-1,\mu =1$, and $\omega =1$, we achieved the two equilibrium points: $(\mathrm{0,0})$ and $(1.007,0.027)$, as represented in Figure 7a In this figure, point $(\mathrm{0,0})$ represents the saddle point, while the point $(1.007,0.027)$ represents the center point.

Case 2: $A_1>0$ and $A_2<0$.

By using the parametric values $a=0.6,b=1,c=3,\mu =1$, and $\omega =1$, we gained the only non-complex equilibrium point as (0,0), as shown in Figure 6b. In Figure 6b, (0,0) shows the center point.

Case 3: $A_1<0$ and $A_2>0$.

By using the parametric values $a=-0.6,b=1,c=1,\mu =1$, and $\omega =1$, we gained the two equilibrium points: (0,0) and $(1.667,0.463)$, as shown in Figure 6c. In Figure 6c, (0,0) shows the center point, while the point $(1.667,0.463)$ represents the saddle point.

Case 4: $A_1<0$ and $A_2<0$.

By using the parametric values $a=-0.3,b=1,c=3,\mu =1$, and $\omega =1$, we gained the only one equilibrium point: (0,0), as shown in Figure 6d. In Figure 6d, point (0,0) shows the saddle point.

In bifurcation analysis and phase portraits, different parameter choices can significantly impact the phase portrait topology.

1. Varying parameters can create or destroy equilibrium points or change their stability properties (e.g., from stable to unstable or vice versa).

2. Parameters can be tuned to critical values, leading to bifurcations, which are sudden changes in the qualitative behavior of the system.

3. Parameters can influence the topology of the phase portrait, such as:

3.1. Creating or destroying limit cycles (closed orbits).

3.2. Changing the stability of limit cycles.

3.3. Creating or destroying homoclinic or heteroclinic orbits.

4. Different parameter regimes can lead to distinct qualitative behaviors, such as:

4.1. Oscillatory vs. non-oscillatory behavior.

4.2. Stable vs. unstable behavior.

6 Chaotic behaviors

Here, we will discuss the chaotic behaviors of the governing model. Chaotic behavior describes the complex, seemingly random, and unpredictable patterns found in systems that follow deterministic rules. We can observe that small changes in the ICs can lead to vastly different outcomes, making long-term predictions difficult.

By introducing the perturbation term $\nu \cos (\varphi t)$ in the dynamical system defined by Equation 42, we get the following perturbed dynamical system according to [31]:

$\begin{array}{ll}\hfill \frac{dG\left(\mathrm{\Omega }\right)}{d\mathrm{\Omega }}& =F\left(\mathrm{\Omega }\right),\hfill \\ \hfill \frac{dF\left(\mathrm{\Omega }\right)}{d\mathrm{\Omega }}& =-A_1{G}^2\left(\mathrm{\Omega }\right)+A_{2}G\left(\mathrm{\Omega }\right)+\nu \cos \left(\varphi t\right).\hfill \end{array}(45)$

Here, $\nu$ and $\varphi$ are the intensity and frequency of the external perturbation term.

The perturbation term represents an external forcing or disturbance that affects the system’s behavior. This term can be interpreted in the context of water wave dynamics. The perturbation term can model the effect of wind on the water surface, where $\nu$ represents the wind stress, and $\varphi$ is the frequency of the wind forcing. The term can also represent the effect of surface tension on the water surface, where $\nu$ is related to the surface tension coefficient. The perturbation term can also model external disturbances, such as waves generated by a paddle or a ship. The perturbation term reflects realistic applications in water wave dynamics. The system can model the behavior of ocean waves under the influence of wind, currents, or other external factors. Understanding the effects of external forcing on water waves is crucial for designing coastal structures, such as seawalls or breakwaters. The system can be used to study the behavior of waves in wave energy harvesting systems, where the perturbation term represents the external forcing that drives the energy conversion. The perturbation term can significantly impact the system’s behavior. The system can exhibit nonlinear resonance, where the external forcing amplifies the system’s response. The perturbation term can lead to chaotic behavior, where the system’s response becomes unpredictable and sensitive to initial conditions. The system can exhibit pattern formation, where the external forcing leads to the emergence of complex spatial or temporal patterns.

We use 2D phase portrait, 3D phase portrait, time series, and Poincaré section to obtain the chaotic and quasi-periodic structures. A perturbation term is taken in the dynamical model defined by Equation 45, which is not taken in the dynamical system defined by Equation 42. This analysis will explain how the frequency term affects the concerned equation. We will investigate the effects of force and frequency of the perturbations while taking the other physical attributes of the overall evaluation as constants.

7 Lyapunov exponent

Here, we aim to explore the Lyapunov exponent of the concerned model. The Lyapunov characteristic exponent (LCE), or Lyapunov exponent, is a tool through which we can determine whether the nearby trajectories in a model converge or diverge. The Russian mathematician Aleksandr Lyapunov, who created the theory of stability of dynamical systems in the late 19th century, is credited with naming the Lyapunov exponent.

In the phase space of the dynamical system, the average distance rate of neighboring trajectories is represented by a real number called the Lyapunov exponent. Numerous applications of the Lyapunov exponent exist in various fields, including biology, engineering, physics, fluid flow, weather patterns, and financial markets. This analysis is used for many models of different fields, including a Konno–Onno model [32], a Schrödinger equation with cubic nonlinearity [29], and a Wazwaz Kaur Boussinesq model [33].

We observed the link between the Lyapunov exponent results and the observed phase portraits. Positive Lyapunov exponents correspond to chaotic regions in phase portraits, characterized by complex, aperiodic trajectories. Negative Lyapunov exponents correspond to stable regions, featuring periodic or quasi-periodic trajectories. Changes in Lyapunov exponents can signal bifurcations, where the system’s behavior changes qualitatively. Lyapunov exponents can help understand the topology of phase portraits, including the existence of attractors, repellers, or saddle points.

8 Sensitivity nature

Here, we discuss the sensitivity of the dynamical model described by Equation 42. The specific values of parameters $a=-2,b=0.6,d=1,p=1,q=1,e=1$, and $\lambda =1$ are selected for this purpose. Moreover, we suppose the following different initial conditions (ICs) of the dynamical system.

(i) $(0.1,0)$; (ii) $(0.01,0.1)$; (iii)$(0.02,0.01)$; (iv) $(1.0,0.1)$; (v) $(1.03,0.5)$; (vi) (0.01,0.001).

The results are explained in Figure 7 according to the abovementioned ICs. In the figure, the red graph denotes G, and the blue graph represents F. It is observed in Figure 7 that small changes in ICs result in a large effect on the concerned model.

9 Modulation instability

Assuming a solution of a (1 + 1)-dimensional Boussinesq water wave model is represented in [28, 29].

$g\left(x,t\right)=\left(\sqrt{\tau }+G\left(x,t\right)\right){\mathrm{e}}^{\iota \tau t}.(46)$

Here, $\tau$ denotes an arbitrary real constant and G(x,t) is a complex-valued function of x and t. Putting Equation 46 into Equation 2, we obtain an equation for G by linearity, given as

$-b{G}_{\text{xxxx}}+c\iota \tau G_x+c{G}_{\text{xt}}-G{\tau }^2+2\iota \tau G_t+G_{\text{tt}}-G_{\text{xx}}-{\tau }^{5/2}=0.(47)$

Now consider a new transformation given as

$G\left(x,t\right)=G_1{\mathrm{e}}^{\iota \left(\rho x-t\lambda \right)}+G_2{\mathrm{e}}^{-\iota \left(\rho x-t\lambda \right)}.(48)$

Here, $\rho$ is a real disturbance wave number, $\lambda$ represents a frequency, while $G_1$ and $G_2$ are the coefficients of linear combination. By using Equation 48 in Equation 47, we get homogeneous equations given as

$\begin{array}{ll}\hfill G_1\left(-b{\rho }^4+c\lambda \rho -c\rho \tau -{\lambda }^2+2\lambda \tau +{\rho }^2-{\tau }^2\right)& =0,\hfill \\ \hfill G_2\left(-b{\rho }^4+c\lambda \rho +c\rho \tau -{\lambda }^2-2\lambda \tau +{\rho }^2-{\tau }^2\right)& =0.\hfill \end{array}(49)$

When the determinant of the system of Equation 49 is set equal to 0, we get the following relation:

$\left(-b{\rho }^4+c\lambda \rho -c\rho \tau -{\lambda }^2+2\lambda \tau +{\rho }^2-{\tau }^2\right)\left(-b{\rho }^4+c\lambda \rho +c\rho \tau -{\lambda }^2-2\lambda \tau +{\rho }^2-{\tau }^2\right)=0.(50)$

Assuming Equation 50, we can discuss types of modulation instability (MI) of Equation 2 given as

$\lambda =\frac{1}{2}\left(c\rho \pm 2\tau \pm \sqrt{-4b{\rho }^4+c^2{\rho }^2+4{\rho }^2}\right).(51)$

A steady-state stable solution is found by Equation 51.

If $\lambda$ has an imaginary part, then the steady-state solution is not stable because the perturbation increases exponentially.

If $\lambda$ is not an imaginary part, then the steady-state solution is stable because the perturbation is small.

Hence, MI of Equation 51 can occur if

$-4b{\rho }^4+c^2{\rho }^2+4{\rho }^2<0.(52)$

Therefore, we obtain the MI gain spectrum given as

$G\left(\rho \right)=2\mathit{I}\mathit{m}\left(\lambda \right)=2\mathit{I}\mathit{m}\left(\frac{1}{2}\left(c\rho \pm 2\tau \pm \sqrt{-4b{\rho }^4+c^2{\rho }^2+4{\rho }^2}\right)\right).(53)$

Physically, MI in shallow water waves can be interpreted as follows:

1. When a wave train propagates in shallow water, it can become unstable due to the interplay between nonlinearity and dispersion. MI can cause the wave train to break down into smaller-scale structures.

2. MI can contribute to the formation of freak waves or rogue waves, which are unusually high and short-lived waves that can pose a significant threat to coastal structures and marine vessels.

3. MI can lead to a redistribution of energy within the wave spectrum, potentially influencing coastal erosion, sediment transport, and wave-induced forces on structures.

4. MI is a manifestation of nonlinear wave interactions, which play a crucial role in shaping the evolution of wave fields in shallow water.

10 Results and discussion

Here, we will represent the obtained results and discussion by comparing them with the existing results. In [10], dark, bright, dark-periodic, and singular-periodic soliton solutions are obtained by using the modified $(G^{\prime }/G^2)$-expansion and F-expansion techniques. In [11], symmetrical, non-symmetrical kink solutions, solitary wave solutions, and Jacobi and Weierstrass elliptic function solutions are gained by applying the extended auxiliary equation scheme. In [12], solitary wave solutions are achieved by utilizing the $\exp (-\varphi (\eta ))$-expansion scheme. In [13], periodic shock wave solutions are obtained with the use of the Jacobi elliptic function expansion scheme. In [16], kink, bright, and dark soliton solutions are achieved by using the generalized Arnous method. In our research, we investigate the dark-bright, dark, bright, periodic, periodic-kink, singular, dark-singular, and other exact soliton solutions by using the Sardar sub-equation technique as shown in Figures 1–7. We performed the modulation instability to obtain the steady-state solutions as shown in Figure 8. We assess sensitivity using the sensitivity analysis as shown in Figure 9. By using bifurcation analysis as shown in Figures 10–13, chaotic analysis as shown in Figures 14–18, and Lyapunov exponent analysis as shown in Figure 19, we discussed the different behaviors of the governing model. These various analyses have not been performed on the fourth-order Boussinesq water wave equation in the literature. The results obtained have applications in different fields of science and engineering.

Figure 8

Figure 8. Gain spectrum of modulation instability for $\rho =\mathrm{2,4,6,8}$ and $b=0.1,\hspace{0.17em}\hspace{0.17em}c=0.02$ in Equation 53. (a) 2D plot. (b) 3D surface plot. (c) Contour plot. (d) Density plot.

Figure 9

Figure 9. Graph of sensitivity demonstration of the concerned model, considering the values of constants along with ICs: (a) $(0.1,0)$, (b) $(0.01,0.1)$, (c) $(0.02,0.01)$, (d) $(1.0,0.1)$, (e) $(1.03,0.5)$, and (f) (0.01,0.001).

Figure 10

Figure 10. Phase portraits of the system shown by Equation 42. (a) 2D streamline plot. (b) Contour plot.

Figure 11

Figure 11. Phase portraits of the system shown by Equation 42. (a) 2D streamline plot. (b) Contour plot.

Figure 12

Figure 12. Phase portraits of the system shown by Equation 42. (a) 2D streamline plot. (b) Contour plot.

Figure 13

Figure 13. Phase portraits of the system shown by Equation 42. (a) 2D streamline plot. (b) Contour plot.

Figure 14

Figure 14. Graph of chaotic behavior of concerned dynamical system given in Equation 45 upon assuming values of parameters $a=-0.6,b=1,c=1,\mu =1,\omega =1,\nu =0.1$, and $\varphi =0.5$ along with the initial condition, $(0.01,0.1)$. (a) 2D Phase portrait. (b) 3D Phase portrait. (c) Time series. (d) Poincaré section.

Figure 15

Figure 15. Graph of chaotic behavior of the concerned dynamical system given in Equation 45 upon assuming values of parameters $a=-0.6,b=1,c=1,\mu =1,\omega =1,\nu =0.01$, and $\varphi =0.2$ along with the initial condition, $(0.01,0)$. (a) 2D phase portrait. (b) 3D phase portrait. (c) Time series. (d) Poincaré section.

Figure 16

Figure 16. Graph of the chaotic behavior of the concerned dynamical system given in Equation 45 upon assuming values of parameters $a=-0.6,b=1,c=1,\mu =1,\omega =1,\nu =0.2$, and $\varphi =2$ along with the initial condition, $(0.01,0)$. (a) 2D phase portrait. (b) 3D phase portrait. (c) Time series. (d) Poincaré section.

Figure 17

Figure 17. Graph of the chaotic behavior of the concerned dynamical system given in Equation 45 upon assuming values of parameters $a=-0.6,b=1,c=1,\mu =1,\omega =1,\nu =0.2$, and $\varphi =2.5$ along with the initial condition, $(0.01,0.03)$. (a) 2D phase portrait. (b) 3D phase portrait. (c) Time series. (d) Poincaré section.

Figure 18

Figure 18. Graph of the chaotic behavior of the concerned dynamical system given in Equation 45 upon assuming values of parameters $a=-0.6,b=1,c=1,\mu =1,\omega =1,\nu =0.04$, and $\varphi =1.5$ along with the initial condition, $(0.03,0.01)$. (a) 2D phase portrait. (b) 3D phase portrait. (c) Time series. (d) Poincaré section.

Figure 19

Figure 19. Graph of the Lyapunov exponent of the concerned dynamical system given in Equation 45 upon assuming values of parameters $a=-0.6,17b=1,c=1,\mu =1,\omega =1,\nu =0.1$, and $\varphi =0.5$. In this figure, the orange graph is for the initial condition (1,0), and the blue graph is for the initial condition (0,1).

11 Conclusion

It is concluded that the Sardar sub-equation scheme was utilized for the concerning model in obtaining distinct kinds of exact solitons to the Boussinesq water wave equation. The results gained are demonstrated with the use of 2D, 3D, contour, and density plots. The results gained have not been studied earlier.

Modulation instability is used to obtain the steady-state solutions for the concerned equation. By using bifurcation analysis, all the phase portraits are discussed. Chaotic behavior is discussed. Sensitivity analysis is used to discuss the sensitivity behavior of the model. The solutions obtained are useful in different fields, including coastal engineering, harbor design, and waves in shallow waters. The Boussinesq water wave equation is useful in the study of water wave behavior, harbor design, coastal dynamics, wave propagation in shallow seas, ocean wave models, and marine environments.

In the future, we can compare the obtained exact solutions with the numerical solutions. We can also obtain the results experimentally in the laboratory. We can study the Sardar sub-equation method by conducting experiments to validate the results obtained. We can perform numerical simulations to verify the accuracy and stability of the solutions, develop modifications or extensions to the Sardar sub-equation method to handle high-dimensional problems, and compare the results obtained using the Sardar sub-equation method with other analytical or numerical methods.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

MR: Conceptualization, Methodology, Writing – original draft. AZ: Conceptualization, Formal Analysis, Methodology, Supervision, Writing – original draft. AA: Funding acquisition, Writing – review and editing. NA: Methodology, Project administration, Writing – review and editing.

Funding

The authors declare that financial support was received for the research and/or publication of this article. This work was supported by the Deanship of Scientific Research, Vice Presidency for Graduate Studies and Scientific Research, King Faisal University, Saudi Arabia (KFU252591).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The authors declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Ibrahim RA, Bekir A, Zahran EHM. On the exact traveling wave solutions to the Akbota-Gudekli-Kairat-Zhaidary equation and its numerical solutions. Eur Phys J Plus (2025) 140(5):354–13. doi:10.1140/epjp/s13360-025-06289-x

CrossRef Full Text | Google Scholar

2. Farooq K, Tedjani AH, Li Z, Hussain E. Soliton dynamics of the nonlinear kodama equation with M-Truncated derivative via two innovative schemes: the generalized arnous method and the Kudryashov method. Fractal and Fract. (2025) 9(7):436. doi:10.3390/fractalfract9070436

CrossRef Full Text | Google Scholar

3. Muhammad J, Rehman SU, Nasreen NYU, Bilal M. Exploring the fractional effect to the optical wave propagation for the extended Kairat-II equation. Nonlinear Dyn (2025) 113:1501–12. doi:10.1007/s11071-024-10139-3

CrossRef Full Text | Google Scholar

4. Kopçasız B. Unveiling new exact solutions of the complex-coupled kuralay system using the generalized Riccati equation mapping method. J Math Sci Model (2024) 7(3):146–5. doi:10.33187/jmsm.1475211

CrossRef Full Text | Google Scholar

5. Zhu C, Al-Dossari M, Rezapour S, Gunay B. On the exact soliton solutions and different wave structures to the (2+ 1) dimensional Chaffee–Infante equation. Results Phys (2024) 57(1):107431. doi:10.1016/j.rinp.2024.107431

CrossRef Full Text | Google Scholar

6. Farooq K, Hussain E, Younas U, Mukalazi H, Khalaf TM, Mutlib A, et al. Exploring the wave’s structures to the nonlinear coupled system arising in surface geometry. Sci Rep (2025) 15:11624. doi:10.1038/s41598-024-84657-w

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Farooq K, Hussain E, Abujabal HA, Alshammari FS. Propagation of nonlinear dispersive waves in shallow water and acoustic media in the framework of integrable Schwarz–Korteweg–de Vries equation. AIMS Math (2025) 10:17543–66. doi:10.3934/math.2025784

CrossRef Full Text | Google Scholar

8. Hussain E, Tedjani AH, Farooq K, Beenish A. Modeling and exploration of localized wave phenomena in optical fibers using the generalized kundu–eckhaus equation for Femtosecond pulse transmission. Axioms (2025) 14:513. doi:10.3390/axioms14070513

CrossRef Full Text | Google Scholar

9. Boussinesq J. Théorie de l’intumescence liquide appelée onde solitaire ou de translation se propageant dans un canal rectangulaire. C R Acad Sci Paris (1871) 72:755–9.

Google Scholar

10. Razzaq W, Zafar A, Nazir A, Junjua M, Awwad FA, Ismail EA. Optical soliton solutions of nonlinear differential Boussinesq water wave equation via two analytical techniques. Results Phys (2024) 64:107898. doi:10.1016/j.rinp.2024.107898

CrossRef Full Text | Google Scholar

11. Helal M, Atef H, Seadawy AR, Zekry MH. Stability analysis of solitary wave solutions for the fourth-order nonlinear Boussinesq water wave equation. Appl Math Comput (2014) 232:1094–103. doi:10.1016/j.amc.2014.01.066

CrossRef Full Text | Google Scholar

12. Akbar MA, Ali NHM. Solitary wave solutions of the fourth order Boussinesq equation through the exp(−ϕ(η))-expansion method. SpringerPlus (2014) 3:344–6. doi:10.1186/2193-1801-3-344

CrossRef Full Text | Google Scholar

13. Haider JA, Muhammad N, Nadeem S, Asghar S. Analytical analysis of the fourth-order Boussinesq equation by traveling wave solutions. Int J Mod Phys B (2023) 37:22350170. doi:10.1142/S0217979223501709

CrossRef Full Text | Google Scholar

14. Naher H, Abdullah FA. The basic (G′/G)-expansion method for the fourth order Boussinesq equation. Appl Math (2012) 3:1144–52.

Google Scholar

15. Aktar A, Rahhman MM, Roy KC. Solitary and periodic wave solutions of the fourth order boussinesq equation through the novel exponential expansion method. Am J Appl Math (2019) 7:49–57. doi:10.11648/j.ajam.20190702.12

CrossRef Full Text | Google Scholar

16. Farooq K, Alshammari FS, Li Z, Hussain E. Soliton dynamics and stability in the Boussinesq equation for shallow water applications. Front Phys (2025) 13:1637491. doi:10.3389/fphy.2025.1637491

CrossRef Full Text | Google Scholar

17. Cheng Y, Dong C, Zheng S, Hu W. Solving nonlinear Boussinesq equation of second-order time derivatives with physics-informed neural networks. Commun Theor Phys (2025) 77:105001. doi:10.1088/1572-9494/adcc8e

CrossRef Full Text | Google Scholar

18. Wazwaz AM, Kaur L. New integrable Boussinesq equations of distinct dimensions with diverse variety of soliton solutions. Nonlinear Dyn (2019) 97:83–94. doi:10.1007/s11071-019-04955-1

CrossRef Full Text | Google Scholar

19. Kong D, Rezazadeh H, Ullah N, Vahidi J, Tariq KU, Akinyemi L. New soliton wave solutions of a (2+1)-dimensional Sawada-Kotera equation. J Ocean Eng Sci (2022). doi:10.1016/j.joes.2022.03.007

CrossRef Full Text | Google Scholar

20. Cinar M, Secer A, Ozisik M, Bayram M. Derivation of optical solitons of dimensionless Fokas-Lenells equation with perturbation term using Sardar sub-equation method. Opt Quan Electron (2022) 54:402. doi:10.1007/s11082-022-03819-0

CrossRef Full Text | Google Scholar

21. Rehman HU, Habib A, Abro KA, Awan AU. Study of langmuir waves for Zakharov equation using Sardar sub-equation method. Int J Nonlinear Anal Appl (2023) 14:9–18. doi:10.22075/ijnaa.2023.27106.3500

CrossRef Full Text | Google Scholar

22. Rehman HU, Akber R, Wazwaz AM, Alshehri HM, Osman MS. Analysis of Brownian motion in stochastic Schrödinger wave equation using sardar sub-equation method. Optik (2023) 289:171305. doi:10.1016/j.ijleo.2023.171305

CrossRef Full Text | Google Scholar

23. Pleumpreedaporn C, Pleumpreedaporn S, Moore EJ, Sirisubtawee S, Sungnul S. Novel exact traveling wave solutions for the (2+1)-Dimensional boiti-leon-manna-pempinelli equation with Atangana’s space and time beta-derivatives via the sardar subequation method. Thai J Math (2024) 22:1–18.

Google Scholar

24. Hossain MN, El-Rashidy K, Alsharif F, Kanan M, Ma WX, Miah MM. New optical soliton solutions to the Biswas–Milovic equations with power law and parabolic law nonlinearity using the Sardar-subequation method. Opt Quan Electron (2024) 56:1163. doi:10.1007/s11082-024-07073-4

CrossRef Full Text | Google Scholar

25. Hossain MN, Miah MM, Alosaimi M, Alsharif F, Kanan M. Exploring novel soliton solutions to the time-fractional coupled Drinfel’d–Sokolov–Wilson equation in industrial engineering using two efficient techniques. Fractal Fract (2024) 8:352. doi:10.3390/fractalfract8060352

CrossRef Full Text | Google Scholar

26. Rezazadeh H, Inc M, Baleanu D. New solitary wave solutions for variants of (3+1)-dimensional wazwaz-benjamin-bona-mahony equations. Front Phys (2020) 8:332. doi:10.3389/fphy.2020.00332

CrossRef Full Text | Google Scholar

27. Alsharidi AK, Bekir A. Discovery of new exact wave solutions to the M-fractional complex three coupled Maccari’s system by sardar sub-equation scheme. Symmetry (2023) 15:1567. doi:10.3390/sym15081567

CrossRef Full Text | Google Scholar

28. Alomair A, Al Naim AS, Bekir A. Exploration of soliton solutions to the special Korteweg–De vries equation with a stability analysis and modulation instability. Mathematics (2024) 13:54. doi:10.3390/math13010054

CrossRef Full Text | Google Scholar

29. Rahaman MS, Islam MN, Ullah MS. Bifurcation, chaos, modulation instability, and soliton analysis of the schrödinger equation with cubic nonlinearity. Sci Rep (2025) 15:11689. doi:10.1038/s41598-025-96327-6

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Alqahtani A, Akram S, Alosaimi M. Study of bifurcations, chaotic structures with sensitivity analysis and novel soliton solutions of non-linear dynamical model. J Taibah Univ Sci (2024) 18:2399870. doi:10.1080/16583655.2024.2399870

CrossRef Full Text | Google Scholar

31. Jhangeer A, Almusawa H, Hussain Z. Bifurcation study and pattern formation analysis of a nonlinear dynamical system for chaotic behavior in traveling wave solution. Results Phys (2022) 37:105492. doi:10.1016/j.rinp.2022.105492

CrossRef Full Text | Google Scholar

32. Chahlaoui Y, Ali A, Ahmad J, Javed S. Dynamical behavior of chaos, bifurcation analysis and soliton solutions to a Konno-Onno model. PLoS One (2023) 18:e0291197. doi:10.1371/journal.pone.0291197

CrossRef Full Text | Google Scholar

33. Zhao S. Chaos analysis and traveling wave solutions for fractional (3+1)-dimensional Wazwaz Kaur Boussinesq equation with beta derivative. Sci Rep (2024) 14:23034. doi:10.1038/s41598-024-74606-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Nomenclature

2D Two-dimensional

3D Three-dimensional

g(x,t) Wave function

ICs Initial conditions

LCE Lyapunov characteristic exponent

MI Modulation instability

NLODE Nonlinear ordinary differential equation

$\mathit{\nu }$ Intensity of the external perturbation term

$\mathit{\omega }$ Velocity of the wave.

PDE Partial differential equation

$\mathit{\varphi }$ Frequency of the external perturbation term

SSE Sardar sub-equation

t Temporal coordinate

x Spatial coordinate