Exact and numerical solutions of the generalized breaking soliton system: Insights into non-linear wave dynamics

1.   Introduction

The investigation of solitonic solutions for non-linear evolution equations (NLEEs) is essential to advancing several scientific and technological fields, including non-linear movement ^[1], ocean engineering ^[2], fibre optics ^[3,4], high-frequency waves ^[5,6], plasma physics ^[7], nanoscale components ^[8], and fluid interaction ^[9,10]. Such singular wave phenomena, called solitons, maintain their shape and energy over large distances, which is crucial for understanding the activities of complex non-linear systems. Solitons are essential for the study of non-linear systems due to their durability and ability to capture essential knowledge on system behavior ^[11]. The designation and scrutiny of solitonic solutions not only enhance our abstract understanding but also encourage pragmatic improvements in innovations based on non-linear phenomena ^[12]. Beyond the span, absorbing phenomena such as rogue waves and transient and localized waves that appear unexpectedly and disperse without trace have been detected in deep probes in non-linear situations ^[13]. Rogue surges have been detected in a variety of contexts ^{[14,15,16,17]}. Given the significant importance of solitonic solutions in understanding and advanced non-linear frameworks throughout multiple objectives, considerable efforts have been made to develop robust approaches to their compilation and verification. Over the past few decades, scientists have developed various methods to study and construct solitonic solutions for NLEEs. These methods include adjusting and developing classical Bäcklund transformations ^[18], Darboux transformations ^[19], as well as the Breather limit approach ^[20], Hirota bilinear method ^[21], generalized exponential rational function method ^{[22,23,24,25]}, exponential rational function approach, sinh cosh approach ^[26], direct algebraic approach ^[27,28,29], improved modified extended tanh function approach ^[30,31], fundamental (G'/G)-expansion method ^[32], and many others.

The main objective of the present research is to examine the (2+1)-dimensional generalized breaking soliton system (GBSS). This system is a robust theoretical framework to model the dispersion, non-linearity, and cross-wave coupling in media such as shallow water waves, optical fibre systems and magnetized plasma. While earlier studies focused on particular cases of GBSS, such as ^{[18,25,33,34,35]}, the present work aims to derive the equation in a more complete and generalized form. By broadening the scope, this study seeks to uncover dynamic behaviors and solutions that may have been indifferent in previous studies of particular cases. The generalized form then ideally offers a unified perspective, which allows one to appreciate different non-linear phenomena that the GBSS encompasses and understand them theoretically better. This work not only improves upon earlier approaches but also contributes to the growing body of literature on more compact and efficient ways of studying non-linear evolution equations. The foundation of this work is the GBSS, which is formally expressed by the following equations:

where the parameters $\alpha, \beta, \gamma, \delta$, and $\lambda$ are real parameters, and the functions $\varPsi = \varPsi(x, y, t)$ and $\varPhi = \varPhi(x, y, t)$ describe the values propagated in the direction of the two structures, respectively. $\varPsi_t$ is what governs the way the wave changes with time. The third derivative with respect to $x$ gives $\alpha \varPsi_{x x x}$, which denotes $x$-axial dispersion. Dispersion is a process that spreads out the wave so that the steepening of the wave due to the non-linear processes is reduced. The mixed derivative $\beta \varPsi_{x x y}$ describes the spread in $x$ and $y$ directions. Therefore, it shows the complicated interaction of waves in both directions. A non-linear term $\gamma \varPsi \varPsi_x$ represents the self-interaction of the wave. It implies that it is increasingly becoming steeper as the wave gets fatter in amplitude. This explains why waves with higher amplitude would travel faster, and the wave would become steep. In addition, the terms $\delta \varPsi \varPhi_x$ and $\lambda \varPsi_x \varPhi$ are interactions in both directions of the two wave fields that integrate the motion of the waves with respect to both axes. They show how waves travelling along one direction affect the motions of the waves travelling along a perpendicular direction. The second equation on the system 1.1, $\varPsi_y = \varPhi_x$, relates the two wave fields $\varPsi$ and $\varPhi$, implying that if there is a wave propagating in the $y$ direction, there must also be a wave propagating in the $x$ direction.

In a fluid dynamic context, the GBSS pertains to the interaction of surface waves that ripple in the water's surface with internal waves that develop in the water below the surface as a result of mechanical force in strata of different densities. Surface waves, which are caused by wind and gravity, typically have a short period and propagate quickly, while internal waves have larger wavelengths, a slower speed, and larger amplitudes. The GBSS captures these interactions by balancing non-linear effects from wave steepening (terms at the $\gamma \varPsi \varPsi_x$ level) with dispersive effects that drive wave dispersion (the $\alpha \varPsi_{x x x}$ term). The cross-interaction terms $\delta \varPsi \varPhi_x$ and $\lambda \varPsi_x \varPhi$ express the energy flow between surface and interior waves while in contact. The situation might include internal waves being magnified by surface waves leading to a similar soliton non-linear interaction where the wave structure changes due to this interaction, such as breaking or waterfalls ^[36]. In the context of non-linear optics as well, this (2+1)-dimensional GBSS describes the interaction of light pulses, or solitons, travelling along optical fibers, where the non-linear refractive index and dispersion are balanced so that solitons can be generated. When these light pulses propagate through the fibre, they experience the influence of a non-linear medium, which induces pulse steepening and the effect of dispersion, which leads to pulse broadening. The GBSS integrates such dynamics through concepts such as $\alpha \varPsi_{x x x}$ for dispersion and $\gamma \varPsi \varPsi_x$, which represents non-linear effects. The cross-interaction terms $\delta \varPsi \varPhi_x$ and $\lambda \varPsi_x \varPhi$ represent the coupling between two light pulses, which can propagate in different polarization modes or neighboring fibers, allowing for energy transfer and complex interactions. One such interaction is the steepening or a breakup into sub-pulses (the breaking soliton phenomena) of one pulse. Such modeling is critical for optical communications: A profound understanding of soliton behavior and interaction can transfer data with better quality and over longer distances ^[37]. Similarly, other applications of GBSS include alfvén waves in magnetized plasma, earthquake dynamics and seismic waves, and the interaction of gravity waves.

In the sections that follow, we analyze the exact and numerical solutions of the system (1.1). In Section 2, we adopt the F-exponential method to obtain the exact solutions to (1.1), which demonstrates that the technique is effective for constructing the analytical solutions for the system. Section 3 deals with obtaining numerical solutions using the finite difference method that supplements the analytical methodology, which simulates the system's behavior under certain conditions.

2.   Exact solutions of the GBSS

The development equation related to the physical fields $\varPsi(x, y, t)$ and $\varPhi(x, y, t)$ in the variables $x$, $y,$ and $t$ is expressed as follows:

First, we derive the travelling-wave solutions of system (1.1), which are formulated as follows:

where $k$, $c$ are arbitrary constants and $w$ is the wave speed. Subsequently, the non-linear formation (2.1) is simplified to the following ordinary differential equation:

where $P_2$ is a polynomial in $\varphi(\zeta), \; \rho(\zeta),$ and their total derivatives with respect to $\zeta$.

In accordance with the modified F-expansion method, the solutions of (2.3) are given by the form

that $\Pi(\zeta)$ is a solution of the following differential equation:

where $\nu_0$, $\nu_1$, $\nu_2$ are specified in Table 1 ^[38], and $\mu_k$, $\sigma_k$ are to be identified after that.

The modified F-expansion method gives a successful approach to obtain solutions using the function $\Pi(\zeta)$. Transformations (2.2) will be used in (1.1) to simplify it into a system of ODEs given by

Integrating the second equation from the system (2.6), $c\, \rho_\zeta = k\, \varphi_\zeta$, with respect to $\zeta$, gives us

By substituting Eq (2.7) into the first equation of system (2.6) and performing one integration with regard to $\zeta$, we obtain

Balancing $\rho_{\zeta\zeta}$ with $\rho^2$ in Eq (2.8) gives the value of $N = 2$. According to the modified F-expansion method with $N = 2$, the general solutions of (2.8) are

Here $\Pi(\zeta)$ is a solution of the following differential equation:

where $\nu_0$, $\nu_1$, $\nu_2$ are specified in Table 1. To investigate the analytic solutions to (2.8), the following steps should be undertaken: First, substitute (2.9) and (2.10) into (2.8), group every coefficient of $\Pi(\zeta)^j$, $j$ = $-4, ..., 4,$ and set their values to zero. The result will give a system of equations for $\mu_0, \; \mu_i, \; \sigma_i$, $i = 1, 2$. Next, solve the resulting system using Wolfram Language. Finally, select the values of $\nu_0, \; \nu_1$, and $\nu_2$, as well as the function $\Pi(\zeta)$ from Table 1, and replace these along with $\mu_0, \; \mu_i, \; \sigma_i$ for $i = 1, 2$ into (2.9) to yield a collection of trigonometric functions and rational solutions to (2.8). With the methodology now established, we are able to proceed with adopting it. This enables us to reach the following conclusions:

Note that we will represent the mathematical statement $k (k \alpha + c \beta)$ with the symbol $\eta$ for the reason of simplifying the notation in the solutions.

All of these solutions refer to types of wave structures predicted by the GBSS. They mainly refer to hyperbolic and trigonometric functions, including $\text{sech}^2$, $\text{csch}^2$, $\text{sec}^2$, and $\text{csc}^2$, each one characterizing localized solitons or periodic waveforms. However, the hyperbolic functions usually represent localized solitons that decay exponentially or quickly. While the trigonometric functions show periodic wave structures with periodic peaks or troughs, the constants and parameters control these solitons' speed, shape, and direction, traveling in the $x$- and $y$-directions over time. While $\varPsi_1$ and $\varPsi_2$ demonstrate basic solitons with exponential decay in space (while also being strongly affected by boundary conditions), $\varPsi_3$ and $\varPsi_4$ show the general complexity of solitons, where they can already be multiple-peaked or even multi-dipped. However, solutions of the form $\varPsi_7$, $\varPsi_8$, and $\varPsi_{11, 12}$, which incorporate trigonometric functions, have a more oscillatory or periodic nature. These solutions are composed of waveforms with sharply formed and highly arcane peaks or oscillations across space; they lead to periodic activity rather than simply exponential decay. Due to combinations of constants with secant or cosecant functions, some solutions like $\varPsi_{9}$ and $\varPsi_{10}$ produce solitons where noticeable peaks or dips are present in their center. These solutions can describe different wave behaviors, including localized solitons and periodic wave trains, and they consist of real physical representations due to the types of functions used (trigonometric/hyperbolic) in many of them.

For the next sections, the analysis will concern the properties and types of soliton solutions. This will be enhanced by an analysis of their dynamics, evolution, and importance within the framework of the non-linear system being investigated (1.1).

Figure 1 illustrates the exact solutions $\varPsi_1(x, t)$ and $\varPsi_3(x, t)$ for a set of values of $n$ taken with step 0.5. $\varPsi_1$ is symmetric and localized around the center at $x = 0$, whose amplitude increases drastically with $n$. This behavior reveals a soliton-like structure that is more localized for larger values of $n$. However, the right figure of Figure 1 shows $\varPsi_3$, which has more localized peaks along with flat regions asymptotically. The solution will have a sign change, creating a dip (negative parts) and a peak (positive parts) around $x = 0$ and, for higher numbers, $n$, sharp changes. This indicates that parameter $n$ controls the generalized shape and general amplitude of the solutions, which can be helpful in finding solutions to non-linear wave propagation that resembles waves such as solitons or localized energy regions.

Figure 2 shows the impact of varying parameters $\alpha$, $\beta$, $\gamma$, $\delta, \lambda$ on the behavior and dynamics of soliton at $t = 2$ for $\varPsi_1$ and $\varPsi_3$. The top and bottom left figures show the effects of $\alpha$ and $\beta$ on $t = 2$, saying that when we increase both parameters for a fixed time moment, the soliton is more enhanced and local with a significant positive change in amplitude. The soliton for negative values becomes wider and smoother, with a lower peak, indicating that the energy is spread over a larger area. The positive values of $\alpha$ and $\beta$ yield peaks higher and sharper, thus indicating energy more localized. The figures on the top and bottom right demonstrate the effects of $\gamma$, $\delta$, and $\lambda$ at $t = 2$, showing an identical behavior. Higher values of these parameters result in sharper solitons with increased amplitude, while lower or negative values lead to wider and less pronounced waves. In both cases, the soliton behavior shows that increasing the parameter values presses into a more compact high energy field configuration and decreasing them causes wave flattening and diluting.

Figure 3 shows the evolution of the soliton solutions $\varPsi_1(x, y, t)$ and $\varPsi_3(x, y, t)$ over increasing time, with five time-steps shown: $t = 0, 4, 8, 12, 16$, while the parameters used in these figures are as follows: $\alpha = 2.0, \, \beta = 0.5, \, \gamma = 0.6, \, \delta = -1.0, \, \lambda = 1.2, \, k = c = 1.0,$ and $n = 1.0$. In figure $(a)$, $\varPsi_1$ moves in the positive $x$-direction in time while keeping the same shape and localization. In the figure $(b)$, $\varPsi_3$ instead evolves in the negative $x$-direction over time. This behavior illustrates the time and space evolution of the soliton solutions dynamic and directionality of $\varPsi_1$ and $\varPsi_3$ as they continue to evolve in opposite temporal and spatial directions.

3.   Numerical solutions of GBSS

This section discusses numerical solutions to the resulting ODE (2.8) using various numerical approaches. This approach aims to validate the accuracy of the analytical solutions. As a representative model, we selected one of the above analytic solutions, (2.11). The non-linear shooting and boundary value problem (BVP) methods, implemented at $t = 0$, utilize the value of $\rho$ at the right endpoint of the domain $\zeta = 0$ while making an initial guess for $\rho_\zeta$. The goal is to get the second boundary condition of $\rho$ at the other endpoint. After getting the numerical results, we compare them with the analytical solution. The MATLAB solvers ODE15s and FSOLVE ^[39] are utilized to explore the numerical solution. To implement the BVP method, we discretize the resultant ODE (2.8) as follows:

To apply the non-linear shooting method, the ODE (2.8) is discretized as follows:

Figure 4 compares the respective analytical solution to the numerical solution obtained by the shooting and the BVP methods. The numerical solutions are evaluated at $t = 0, \, y = 1.5,$ and the spatial range $x\in[-55, 15]$. The analytical solution roughly matches both numerical solutions, including the narrow symmetric peak after $x = -1.55$. There is a comparison of methods studied and numerical problems solution efficiency. Results from the shooting or the BVP approaches provide a numerical solution, which will be used as the basis for the following numerical scheme. More specifically, these results will serve as initial conditions for the numerical scheme presented in the next section. This guarantees that the numerical scheme always begins with the correct and cohesive data, aiding in convergence for the controlling equations around the problem parameters and boundary conditions.

3.1.   Numerical scheme

This section employs the finite-difference method to derive the numerical results of system (1.1) within the domain $[x_0, \, x_M]\times [y_0, \, y_N]$. In this context, $x_0$, $y_0$, $x_M$, and $y_N$ define the bounds of the rectangular domain along the $x$- and $y$-axes, with $T_f$ signifying a specific time. The domain is divided into $(M+1)\times (N+1)$ mesh points:

where $h_x$ and $h_y$ denote the step sizes for the $x$ and $y$ domains, respectively. The system (1.1) is transformed into an ODE system by discretising the spatial derivatives while maintaining the continuity of the time derivative. Performing this results in

subject to the boundary conditions:

where $\alpha, \beta, \gamma, \delta, \lambda, k$, and $c$ are user-defined parameters. Equation (3.4) allows for the use of fictitious points in the calculation of spatial derivatives at the boundaries of the domain. In the numerical results provided in this section, the parameter values are established as follows: $\alpha = 0.6, \, \beta = 1.1, \, \gamma = -0.6, \, \delta = 1.3, \, \lambda = 0.3, \, x_0 = -20, \, y = 0 \to 5, \, x = -55 \to 15,$ and $t = 0 \to 20$.

We used MATLAB version R2023b on a PC equipped with an Intel® Core^TM i7 processor operating at 3.40 GHz, 16 GB of RAM, and running Windows 11 to carry out our simulations. To solve implicit ODEs (3.3), the solver ODE15i uses the implicit, multi-step numerical method called the backward differentiation formulas (BDF). The technique has proven very effective in stiff systems and equations with complicated implicit conditions as it approximates derivatives with the previous time steps of the solution. Algorithms like Newton's iteration are used to solve non-linear systems resulting from the implicit formulation. The ODE15i solver is well-known for its accuracy and efficiency, especially with systems that are challenging for explicit methods. This is demonstrated by Figure 5.

Figure 6 illustrates the comparison of the exact solution and the numerical solution of $\varPsi_3(x, y = 1.5, t = 20)$ as the number of points varies. As can be seen, Figure 6 provides an appropriate agreement between numerical and exact solutions, which is evidence of the high reliability of numerical models. It also shows improved numerical precision as the spatial resolution is now finer.

Figure 7 offers a graphical illustration of how closely the numerical computations approximated the third wave function $\varPsi_3$. The comparison enables us to appreciate the degree of effectiveness and correctness of the numerical method and highlights the inherent problems and quantifies the extent of stability of the solution within the specified domain. The following values were taken for the calculations: $\alpha = 0.6$, $\beta = 1.1$, $\gamma = -0.6$, $\delta = 1.3$, and $\lambda = 0.3$. Additionally, the spatio-temporal limits are $x \in [-55, 15], y \in [0, 5]$, and $t \in [0, 20]$, respectively. Finally, in the respect to the resolution numerical parameters, a total of $M = 2300$ and $N = 150$ are set for the $x$- and $y$-directions, respectively, enabling for high precision to be maintained throughout the range set earlier on.

3.2.   Stability of the numerical scheme

Von Neumann analysis is applied to determine the stability of the scheme (3.3). Von Neumann analysis, sometimes referred to as Fourier analysis, is applied exclusively in the context of linear systems. Therefore, we propose that the linear version is represented by

where $r_0 = \gamma \varPsi$, $r_1 = \delta \varPsi$, and $r_2 = \lambda \varPhi$ are constant quantities defined as

We assume a solution form:

As a result, we are able to write

Substituting (3.6) into (3.5) and doing some operations yields

Hence,

where

Thus,

Equation (3.9) indicates the absence of growth, thus the stability condition of the von Neumann analysis is satisfied. As a result, the scheme is unconditionally stable.

3.3.   Error analysis

To effectively evaluate the accuracy of numerical scheme (3.3), we perform a comprehensive investigation of the truncation error, $T(x, y, t)$, using Taylor series expansions. This method allows a systematic analysis of the error's behavior and provides a basis for evaluating the scheme's accuracy under various circumstances. Suppose that the error is

where $\varPsi(x_m, y_n, t_{k+1})$ denotes the analytical solution, whereas $\varPsi_{m, n}^{k+1}$ denotes the approximate numerical solution. Substituting (3.10) into (3.3) gives

We are given that

Thus

Consequently, the truncation error of the numerical scheme is

Table 2 and Figure 8 demonstrate the relationship between numerical accuracy and computational effort across varying ($h_x$) at $t = 20$. Table 2 illustrates the relative error in the $L_2$ norm with corresponding CPU time, which indicates that fine ($h_x$) will lead to increased accuracy, although with more significant computational needs. This trade-off points out the computational expense that is linked to achieving reduced errors. Figure 8 displays the convergence of the numerical method, which exhibits a slope of roughly 1.9, which indicates its second-order accuracy. The analysis emphasises a specific variable ($y = 0.5$) and computational domain ($x \in [-55, 15]$), which provides an accurate measurement of the method's performance. This finding highlights the necessary balance between accuracy and efficiency in computational simulations.

4.   Conclusions

We addressed numerical and analytical investigations for the (2+1)-dimensional GBSS model. In this paper, we applied more improved analytical and computational approaches, such as the modified F-expansion, and were able to extract exact solutions. These solutions exhibit the remarkable property of solitons, which is the ability to keep their shape and remain stable over particularly large distances for broad periods of time, which plays a key role in cases of non-linear wave dynamics. Theoretical models were validated by the use of robust numerical algorithms such as finite difference and boundary-value problem approaches. The results that were predicted during the numerical simulations corresponded to predictions made by the theories that were posed, and our analytical approaches were also proved correct. Moreover, the results suggest potential practical applications of these solutions as solitons. Consider the case where solitons can travel through fibers, neither changing shape nor dissipating energy, in optical communication systems. If so, they could increase channel capacity and lower bit-error rates, enabling data to travel more smoothly. Similarly, the insights gained may be helpful in fluid mechanics and plasma physics, where managing non-linear wave propagation is crucial. The contribution of analytical and numerical solutions not only contributes to a deeper understanding of the complex processes describing the GBSS but also shows the efficiency of our approaches in modeling complicated non-linear phenomena. In general, this study serves as a theoretical foundation for further scientific advancements and offers practical means for technological advancements in spheres that require extensive insights into non-linear wave behaviors.

Author contributions

Amer Ahmed: Conceptualization, methodology, formal analysis, validation, and writing original draft; Abdulghani R. Alharbi (arharbi@taibahu.edu.sa): Conceptualisation, investigation, data curation, and visualization; Ishak Hashim (ishak_h@ukm.edu): Supervision and writing review and editing. All authors have read and approved the final version of the manuscript for publication.

Use of Generative-AI tools declaration

The authors declare that they have not used artificial intelligence (AI) tools to create this article.

Acknowledgments

The author, Amer Ahmed, is grateful to Taibah University for its generous support in the form of a fully paid scholarship, which significantly helped in successfully completing this work.

Conflict of interest

The authors declare no conflicts of interest.