The water wave equation (WWE) was introduced by Boussinesq in 1871 [1]. U t t − U x x − σ U 2 x x − μ U xxxx = 0 . (1.1)

This classic Boussinesq equation (BE) defines the shallow-water wave (SWW) solution interaction process. This equation incorporates various waves and shallow water phenomena in fluid dynamics, including shoaling, diffraction, refraction, and weak non-linearity. In addition to fluid dynamics, it is essential in many disciplines of physics, like ions found in waves in plasma, vibrations in non-linear strings, one-dimensional non-linear lattice waves, and the propagation of long waves in shallow water [2]. This study demonstrates the intricate process of how rogue waves are formed and spread in higher dimensions. In addition, we have created a new BE that can be integrated and has varied dimensions [3]. These equations provide a wide range of soliton solutions, contributing to our understanding of wave processes in many physical environments [4]. The focus will be on the fourth-order nonlinear BE. U t t − σ U 2 x x − μ U xxxx + ν U x t − U x x = 0 . (1.2)

Here, U ( x , t ) represents the surface tension of the water wave, σ denotes the nonlinearity coefficient, μ , and ν are the dispersion coefficients. We see Equation 1.2, originally proposed by Wazwaz and Kaur [4], as completely solvable. Several researchers have investigated different outcomes for nonlinear WWE. For instance, Wang et al. [5] developed advanced Boussinesq-type equations that accurately represent wave dynamics in porous media and apply them to wave propagation in deep water. Fan et al. [6] conducted a study on the use of the widely used G ′ G -expansion method to analyze the unique solutions of non-linear evolution problems, including the sine-Gordon, Klein–Gordon equation, and BE. Numerous travel solutions were introduced by Kumari cite kumari2020abundant. Jun et al. [7] derived the Backlund transformation and Painleve expansion of Equation 2.1 to employ various solutions. Kumar et al. [8] derived the Lie point symmetry along with the lump and breather solution of Equation 2.1. Understanding wave propagation in shallow seas requires an understanding of the fourth-order nonlinear Boussinesq water wave equation model, which incorporates higher-order nonlinear and dispersive features. This increase enhances the accuracy of wave predictions, making the model particularly helpful for evaluating waves with larger amplitudes and longer wavelengths [9]. This field’s applications include coastal engineering and environmental science [10], giving essential insights for coastal development [11], navigation, and disaster prevention [12], including tsunamis and storm surges. The model’s complex features allow more realistic modeling of complicated wave interactions, assisting both theoretical research and practical coastal management.

The nonlinear Integrable Boussinesq Water Wave Equation (IBWWE) has emerged as a significant model in various physical phenomena due to its ability to incorporate both nonlinear and dispersive effects with high order accuracy. Its applications span multiple disciplines, including shallow water wave theory [13], nonlinear lattice wave theory [14], coastal engineering [15], and ion-acoustic wave dynamics in plasmas [16]. Due to its high-order structure, the IBWWE provides a refined representation of wave behavior in optical solitons in fiber media and related photonic systems. Given its broad applicability, continued investigation of the IBWWE’s soliton structures and analytical properties remains a subject of substantial interest.

Solitary waves [17], or solitons, are self-reinforcing waves that retain their shape while moving at a constant speed. These waves occur in certain nonlinear systems and are solutions to specific nonlinear partial differential equations. A key characteristic of solitary waves is that they do not dissipate or spread out as they travel, unlike typical wave packets that tend to disperse and lose their form over time. Recent developments in mathematical modeling reflect a growing emphasis on accurately characterizing the complex behaviors observed in nature and physical systems. In this context, nonlinear partial differential equations (NLPDEs) offer a robust framework for describing diverse dynamical systems. The development of advanced analytical and computational techniques has facilitated the derivation of exact solutions, enabling the deeper understanding of soliton dynamics, nonlinear wave propagation, and pattern formation. Solitary waves play a key role in several scientific and technical sectors owing to their unusual ability to keep their form and speed across vast distances and via interactions. In physics [18] and engineering [19], solitons are used to simulate stable wave phenomena in nonlinear optics [20], fluid dynamics [21], and plasma physics [22, 23], such as optical pulses in fiber-optic [24] and ion-acoustic waves in space plasmas. They are also significant in biological [25] systems for understanding nerve signal transmission and pattern generation, and in chemistry [26] for characterizing reaction-diffusion processes. In mathematics [27], solitons give insight into nonlinear dynamics [28], chaos theory [29, 30], and integrable systems.

NLPDEs develop as especially significant assets in this scientific quest. Many academics have devoted their efforts to examining distinct NLPDEs to increase their comprehension of the demonstrated behavior in the researched natural phenomena. Recent assessments have involved inquiries into the nonlinear Helmholtz equation [31], complex cubic Nonlinear Schrodinger equation [32], Klein-Fock-Gordon equation [33], Kaup–Newell Model [34], Caudrey–Dodd–Gibbon equation [35]. Studying the single-wave solutions of NLPDEs is crucial for generating improved insights and knowledge of the underlying mechanism and its valuable usage. Therefore, various academics have established novel approaches to investigate these NLPDE replies. Plenty of strong techniques such as EHF technique [36], Darboux transformation [37], exp-function method [38], generalized Kudryashov method [39], extended trial equation method [40], Hirota bilinear method [41], extended Jacobian method [42], extended direct algebraic method [43], NAE method [44], improved extended fan-sub equation method [45], multivariate generalized exponential rational integral function method [46].

Although significant advancements have been made in the computational and symbolic treatment of NLPDEs, analytical exploration of the fourth-order IBWWE, especially in its general form involving dispersive and mixed derivative terms, remains limited. Many of the available methods are limited in scope, only producing restricted forms of solutions. There is a need for comprehensive methods that can produce a border class of exact soliton solutions, including dark, bright, periodic, and compound solitons, while also analyzing the qualitative behavior of stability.

This paper proposes an integrated application of the three advanced solution methods: the Generalized Arnous method [47], the Modified Sub-Equation method [48], and the New Kudryashov method [49], to derive the border spectrum of soliton solutions to the fourth-order IBWWE. Furthermore, a linear stability conducted to assess the robustness of the obtained wave structures. To our knowledge, the combined effects of these three techniques on Equation 1.2, along with he detailed graphical and stability analysis, have not been comprehensively reported in the exciting literature. The Generalized Arnous method is an effective technique for obtaining rational-logarithmic solutions characterized by intricate nonlinear behaviors. The Modified Sub-Equation method utilizes the Riccati type transformations, is particularly suited for the construction of periodic and singular waveforms. The Kudryashov method, recognized for its symbolic strength, is for formulating exact solutions in polynomial-exponential form. Collectively, these methods offer a comprehensive analytical framework to yield a richer and more diverse set of analytical solutions, including mixed and compound solutions.

The article is summarized as follows: Section 2 outlines the mathematical analysis required to transform the nonlinear partial differential problem into an ordinary differential equation. Section 3 examines the Generalized Arnous method, its applications, and includes graphical representations. Section 4 highlights the application of the Modified Sub-Equation method. Section 5 delves into the mathematical framework and applications of the New Kudryashov method. Section 6 focuses on the stability analysis. Section 7 discusses the graphical representation of solutions, and finally, Section 8 concludes the study.

Consider a general NLPDE has the following form [27, 50]: Y U , U t , U x , U x t , U x x , … = 0 . (2.1)

Its NODE will be P Ξ , Ξ ′ , Ξ ′ ′ , … = 0 . (2.2)

Consider the traveling wave ansatz solution to simplify the NLPDEs into NLODEs [51, 52]. U x , t = Ξ ξ ; ξ = ω x − η t . (2.3)

Here, ω represents the wave number, and η denotes the wave speed.

The fourth-order nonlinear differential Boussinesq water wave equation is given as [53]: Ξ t t − σ Ξ 2 x x − μ Ξ xxxx + ν Ξ x t − Ξ x x = 0 (2.4)

Now by using Equation 2.3 in Equation 2.4 we get: μ ω 4 Ξ 4 − Ξ ′ ′ − ν ω η − ω 2 + η 2 + 2 σ ω 2 Ξ Ξ ′ ′ + 2 σ ω 2 Ξ ′ 2 = 0 . (2.5)

After integrating Equation 2.5 twice w. r.t ξ we get [54, 55]: σ ω 2 Ξ 2 + μ ω 4 Ξ ′ ′ − Ξ − ν ω η − ω 2 + η 2 = 0 . (2.6)

The basic steps of the generalized Arnous (GA) method are as follows [47].

Step 1: The (GA) method provides the solution of Equation 2.3 as follows:

Here ρ , A, and γ are real constants.

Step 2: By balancing the non-linear term with the highest order derivative in Equation 2.6, the positive integer N is determined for Equation 3.1.

The basic steps of the Modified sub-equation (MSE) method are as follows [56].

Step 1: The (MSE) method provides the solution of Equation 2.3 as follows:

c 0 , c j (for j = 1, 2, … , N) are non zero constants, with the condition c N ≠ 0 , and the function g ( ξ ) in Equation 4.1 satisfied the relation: g ′ ξ = λ 2 g 4 ξ + λ 1 g 2 ξ + λ 0 . (4.2)

Here λ 0 , λ 1 , and λ 2 ≠ 0 are real constants. The answer to Equation 4.2 as follows.

Case 1: When λ 0 = 0 , λ 1 > 0 , and λ 2 ≠ 0 then,

Here are some important steps of the new Kudryashov method (NK).

Step 1: The NK method provides the solution of Equation 2.6 as:

where the coefficients l i for i = 0, 1, 2,…, N are constants to be determined such that l N ≠ 0 , and g ( ξ ) = 1 a B δ ξ + b B − δ ξ is the solution of the following non-linear ODE: g ′ ξ 2 = δ ⁡ ln B g ξ 2 1 − 4 a b g 2 ξ . (5.2) g ′ ′ ξ = δ 2 ⁡ ln B 2 g ξ 1 − 8 a b g 2 ξ . (5.3)

here constants a, b, δ , and B are all non-zero, with B > 0 and B ≠ 1 .

Step 2: Using the homogeneous balance principle, we may get the positive integer N by balancing the highest-order derivative and nonlinear variables in Equation 2.3.

In this section, we will discuss the stability of Equation 1.2. Consider a perturbed solution of Equation 1.2 has the form [57, 58]. Ξ x , t = P + λ U x , t . (6.1)

For any constant value of P, it is obvious that Equation 1.2 possesses a stable solution. U is the function of x,t, and λ is a real constant. By Inserting Equation 6.1 in Equation 1.2, we obtain the following result λ U t t − σ λ 2 U x x 2 − μ λ U xxxx − ν λ U x t − λ U x x = 0 . (6.2)

Linearized Equation 6.2. λ U t t − μ λ U xxxx − ν λ U x t − λ U x x = 0 . (6.3)

Suppose that Equation 6.3 has the solution of the from U x , t = e ι m x − s t . (6.4)

Here, m represents the normalized wave numbers, and s represents the dispersion relation. By inserting Equation 6.4 into Equation 6.3, the following result is obtained s m = μ m 4 − m 2 . s m = − m θ − μ m 4 + m 2 . (6.5)

Now we’ll look at the dispersed characteristics shown in Equation 6.5. The dispersion is stable if the real component of Equation 6.5 is negative for all m values. If it is positive, the dispersion is unstable. If it is zero, the dispersion is minimal.

The graphical solutions produced using the Modified Sub-Equation, Generalized Arnous, and Kudryashov method to illustrate the presence of a wide range of soliton solutions within the framework of IBWWE. These visualizations also demonstrate how key parameters affect the wave behavior. Specifically, increasing the dispersion coefficient μ leads to sharper and narrower wave fronts, reflecting the intensification of dispersive mechanisms. Similarly, a higher nonlinearity parameter σ yields increased amplitude and steepness, highlighting enhanced nonlinear interactions. The parameter ν plays a critical role in shaping the symmetry and phase behavior of the solutions the solution, sometimes introducing asymmetry or a deformed wave shape. These findings suggest that adjusting the model parameters allows control over wave localization, structural properties, and stability, offering practical insight into physical systems modeled by the IBWWE. The pictorial appearance of the solutions produced is investigated in this section. Specific values are supplied to the unknown constants to construct 3D and 2D graphs of the resulting solutions. The figures depicted in part (a) reflect a 3D plot, while part (b) represents the 2D line graph of the solutions, part (c) displays the Contour graph, and part (d) depicts the Polar plot. In Figure 1, the wave solution is visualized through 3D surfaces, 2D profiles, and contour maps under the parameter configuration Ξ 01 ( x , t ) :   ρ = 1 , θ = 0.2 , μ = 1 , ω = 0.3 , β 2 = 1 , σ = 1 , η = 1 , A = 0.5 , with the phase variable defined as ξ = ω x − η t . The plots are generated for time slices t = − 1,0,1 . Figure 2 presents the wave evolution corresponding to a dark solitary structure under the same parameter configuration Ξ 02 ( x , t ) , where the solution is illustrated through 3D plots, 2D line profiles, and contour maps for t = − 3,0,3 . In Figure 3, the bright solitary wave behavior is captured using the parameters Ξ 01 ( x , t ) :   λ 0 = 0 , λ 1 = 1 , λ 2 = 2 , μ = 2 , σ = 1 , ρ = 0.8 , η = 0.5 , ω = 1 , ν = − 8 , with ξ = ω x − η t , and evaluated over the time domain t = − 5,0,5 . Figure 4 displays a dark-bright solitary wave structure governed by the constants Ξ 04 ( x , t ) :     λ 0 = 0.005 , λ 1 = − 0.1 , λ 2 = 0.5 , μ = 2 , σ = − 2.5 , ρ = 1 , η = 0.5 , ω = 0.5 , ν = 0.56 , with ξ = ω x − η t . The visual representation is provided for t = − 5,0,5 . In Figure 5, the anti-kink solitary wave is visualized using Ξ 06 ( x , t ) : λ 0 = 0.001 , λ 1 = − 0.08 , λ 2 = 1.2 , μ = 2 , σ = − 2.5 , ρ = 2 , η = 2 , ω = 1.4 , ν = 0.72 − 0.65 i , with ξ = ω x − η t . The plots correspond to time levels t = − 3,0,3 . Figure 6 illustrates the bright solitary wave pattern for Ξ 14 ( x , t ) : λ 0 = 0 , λ 1 = 0.5 , λ 2 = − 0.5 , μ = 2 , σ = 1 , ρ = − 1.5 , η = 1 , ω = 0.5 , ν = 1 , with the phase ξ = ω x − η t and plots evaluated for t = − 3,0,3 . In Figure 7, the solution evolves under the parametric structure Ξ 15 ( x , t ) : λ 0 = 0 , λ 1 = 0.9 , λ 2 = − 0.1 , μ = 1 , σ = − 1 , ρ = − 1 , η = 1 , ω = 0.5 , ν = 1.05 , with ξ = ω x − η t . The visualization is provided for time slices t = − 1,0,1 . Figure 8 depicts a bright solitary wave structure governed by Ξ 1 ( x , t ) : λ 0 = 0 , λ 1 = 0.9 , λ 2 = − 0.1 , μ = 1 , σ = − 1 , ρ = − 1 , η = 1 , ω = 0.5 , ν = 1.05 , with ξ = ω x − η t , and evaluated at t = − 3,0,3 . In Figure 9, the dark solitary wave pattern is illustrated for Ξ 2 ( x , t ) :   λ 0 = 0 , λ 1 = 0.9 , λ 2 = − 0.1 , μ = 1 , σ = − 1 , ρ = − 1 , η = 1 , ω = 0.5 , ν = 1.05 , using the same phase ξ = ω x − η t . The solution behavior is shown for time levels t = − 3,0,7 . Finally, Figure 10 provides the 3D surface and 2D projection visualizations reflecting the stability features of the system described by Equation 6.5, evaluated under the parameters μ = 0.5 , m ∈ [ − 10,10 ] , and θ ∈ [ − 20,20 ] .