1. Introduction
In the present age of technological advancements, significant progress has been analyzed on the subject of solitons. They play a key role in transferring data over extended distances through fiber optics. Solitons also matter significantly as they are pivotal for the advancement of computer systems’ computing power, while they present a wide variety of applications. Such applications cover image processing, data analysis, neurology, and fluids, among other fields [
1,
2,
3]. As a result, experts have actively explored the evolution of soliton solutions for nonlinear partial differential equations (NLPDEs). Various effective methods have been employed to obtain soliton solutions for NLPDEs, including the Darboux transform [
4], Hirota’s approach [
5], the Painlevé test [
6], the sub-equation method [
7], the unified approach [
8], the bilinear approach [
9], and multiple techniques discussed in [
10].
Among these approaches, Lie analysis stands out as one of the systematic approaches for deriving closed-form solutions to NLPDEs. In recent decades, the application of Lie symmetries has been widespread, extending to numerous instances in both mathematics and engineering [
11]. Upon conducting Lie analysis on nonlinear models, symmetry reductions and group-invariant solutions become accessible. These resultant symmetries frequently prove instrumental in diminishing both the number of independent variables in PDEs and the order of an ordinary differential equation (ODE). The academic community has demonstrated significant interest in the exploration of symmetries, leading to a substantial body of noteworthy research in this field. In their work, Kumar et al. [
11] applied the Lie approach to examine the generalized Kadomtsev Petviashvili (KP) equation, while the authors in reference [
12] scrutinized the Oskolkov model using Lie analysis and extracted several soliton structures.
Substantial investigations have previously explored the behavior of solitons, utilizing various nonlinear frameworks like the Sakovich equation [
13], Schrödinger equation [
14], and Zakharov–Kuznetsov modified equal-width equation [
15]. The integrable equations in (3 + 1) dimensions have been the subject of increased research attention lately. These equations play a key role in elucidating the physics behind numerous scientific and engineering phenomena. Various nonlinear extended equations in (3 + 1) dimensions, including the Korteweg–de Vries (KdV) equation, modified KP equation, and others, have been formulated in response to this growing interest.
Akinyemi [
16] examined the equation in (2 + 1) dimensions:
and its extended version, comprising two additional linear terms:
Here, the constants
, and
in the given context are unrestricted real values. It is worth highlighting that when
, Equation (
2) simplifies to the Boussinesq equation:
Equations (
1) and (
2) have been demonstrated to exhibit Painlevé integrability, and their multiple solitons have also been obtained [
16]. In the ongoing study, our intention is to address an extended version of Equation (
1) [
17]:
where
. It is evident that Equation (
4) is constructed by introducing four additional linear terms, specifically,
, and
, to Equation (
1). Furthermore, the coefficients
, and
m are random real values. In their study [
17], the authors investigated the integrability criteria for the discussed models (
4) through the use of the Painlevé test. The study concluded by analyzing a set of lump solutions for the suggested equation.
This study deals with the use of the new auxiliary equation (NAE) method [
12,
18] in order to analyze as well as to generate soliton structures for the considered model. It is worth mentioning here that the methodology discussed in this work has not been applied in prior research to address Equation (
4). Through the adoption of this novel approach, the obtained results manifest into different forms, such as exponential, trigonometric, hyperbolic, and rational functions. Additionally, the findings are depicted graphically using three-dimensional, two-dimensional, and density plots, thereby improving the clarity of the results. The applications of dynamical systems are diverse, covering fields such as biology, engineering, and economics [
19].
In general, scholars have demonstrated a growing interest in investigating many aspects of dynamical systems, including bifurcation analysis, chaos, and sensitivity. Bifurcation studies how an orbit varies with respect to different parameters. The majority of continuous dynamical systems using differential equations (DEs) have parameters. It is possible for a little change in a parameter to have a big effect on the response. It is a qualitative examination of the planar dynamical system that takes place for distinct values of the physical parameters that are present in the system. Equilibria, periodic orbits, and other invariant sets typically experience changes in their local stability characteristics during a bifurcation.(Please revise this part to reduce duplication.) It has various applications in different fields such as climate science, ecology, economics, and finance. In deterministic nonlinear systems that behave irregularly and erratically, chaos can also arise. Numerous researchers have studied the chaos and bifurcation phenomena within planar systems. For instance, Kazmi et al. [
20] reported the soliton solutions and examined the dynamics of bifurcation and chaos in the q-deformed Sinh–Gordon model. Similarly, Rafiq et al. [
21] investigated the qualitative nature of shallow water equations, uncovering multiwave solitons. Additionally, Hosseini et al. [
22] analyzed the dynamic characteristics of the Schrödinger equation. The chaotic phenomena of the suggested equation is analyzed through the presentation of two-dimensional plots, time plots, and Lyapunov exponents analysis. Furthermore, a sensitivity [
23] of the system has been conducted employing the Runge–Kutta (RK) method.
In this work, we an explored extended integrable wave equation (EIWE) (
4) from three distinct viewpoints:
Initially, the Lie approach is utilized to identify translational symmetries. Subsequently, leveraging these symmetries, the model undergoes a transformation into an ODE. The solution to the resulting ODE is obtained through the application of the NAE approach, facilitating the examination of soliton structures.
Secondly, it is essential to explore the parameter-dependent analysis of the studied system through the lenses of bifurcation and chaos. A detailed explanation of bifurcation for the undisturbed dynamic system is provided, illustrated via phase plots. Additionally, it explores chaotic tendencies in the perturbed dynamical system employing the phase plot, time plot, and Lyapunov exponents tool.
Lastly, we examine the sensitivity of the considered model under various initial conditions. Slight modifications to the initial values can lead to a subtle change in the outcome of the system. Therefore, our observation suggests that the suggested system exhibits sensitivity, though not to an extreme extent.
The structure of the manuscript is as follows: In
Section 2, we identify the Lie symmetries governing the equation and transform the model into an ODE using the extracted symmetries. The analytical solution of the model and results obtained from the analysis are discussed in
Section 3. A comprehensive investigation of the dynamic characteristics of the proposed equation, utilizing phase portraits of bifurcation, is conducted in
Section 4. Chaotic phenomena are explored in
Section 5.
Section 6 is dedicated to a discussion on the multistability analysis of the system. Sensitivity analysis is covered in
Section 7. Finally, the conclusion is provided in the last section.
2. Lie Symmetries
Lie analysis focuses on discerning transformations that maintain the integrity of a given DE. These symmetries are often derived from Lie operators, which collectively form a Lie algebra [
24,
25]. The application of these operators to the DE allows for the identification of symmetries that uphold its inherent structure. Extracting symmetries and reducing the dimensions of the PDE involves the following key steps:
Step 1: Identify classical Lie point symmetries for the suggested model.
Step 2: Construct an algebra (in our case, an abelian algebra) based on the identified symmetries.
Step 3: Determine similarity variables corresponding to each symmetry.
Step 4: Utilize the obtained symmetries to reduce the PDE to an ODE.
Step 5: Derive solutions in the form of traveling waves from the resulting ODE.
Let us examine a group of infinitesimal transformations characterized by a single parameter in Lie theory.
where group parameter
and
,
,
,
, and
are the coefficient functions. An infinitesimal generator associated with the aforementioned Lie group can be expressed as follows:
The invariance condition for Equation (
4) with
becomes
where
is the fourth prolongation of
, defined as
where
Here
,
,
, and
are total derivatives with respect to
x,
y,
z, and
t; it can written as
Combining system (
7) and Equation (
6) into Equation (
5) produces a set of equations for determination. Solving this system results in the derivation of the following symmetries for Equation (
4):
Commutator table:
| | | | |
| 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 |
where
Adjoint table:
| | | | |
| | | | |
| | | | |
| | | | |
| | | | |
where
Traveling Wave Structures
One can notice that the set of
forms an abelian subalgebra. Our objective is to identify the traveling wave structures in the given model, which are associated with translational symmetries of the form:
Through the application of the above linear combination denoted as the
of symmetries, we establish the transformation in the following manner:
Here,
and
represent the characteristics of the traveling wave, specifically referring to its shape and velocity. Additionally,
serve as random parameters. Upon substituting the expression from Equation (
8) into Equation (
4), we obtain the resulting equation:
By integrating Equation (
9) twice with respect to
, and setting the integration constants to zero, we obtain
3. Analytical Solutions
In the following section, we examine the soliton solutions of Equation (
4) by employing the NAE approach [
12,
18]. This approach encompasses an initial solution that can be expressed as follows:
where
is a random real parameter such that
. The positive integer
j serves as the balancing number, determined through the principle of homogeneous balance. This principle entails balancing the highest order derivative with the nonlinear term in the ODEs. By setting the terms
and
equal in Equation (
10) as
, it yields
. Utilizing this approach for
, the initial solution becomes
Here,
,
, and
are arbitrary parameters to be determined. Additionally, an auxiliary function
is introduced, which satisfies the following ODE:
In this context,
,
, and
represent arbitrary constants. Upon substituting the values from Equations (
12) and (
13) into Equation (
10) and equating the coefficients of
to zero, a set of equations is obtained. Solving this system, with the assistance of software such as Mathematica, yields solutions for
,
,
, and
a. The obtained solutions are as follows:
By substituting the aforementioned values into Equation (
12) and applying the transformation outlined in Equation (
8), we derive the solutions for Equation (
4).
Type 1: For
and
, we have
Type 2: For
and
, we have
Type 3: For
,
, and
, we have
Type 4: For
,
, and
, we have
Type 5: For
and
, we have
Type 6: For
=
and
, we have
Type 7: For
, we have
Type 8: For
,
, and
, we have
Type 9: For
,
, we have
Type 10: For
,
, we have
Type 11: For
, we have
Type 12: For
, we have
Type 13: For
, we have
Type 14: For
, we have
Type 15: For
, we have
Type 16: For
, we have
Type 17: For
,
, we have
Type 18: For
, we have
Type 19: For
,
, we have
Type 20: For
, we have
Graphical Representation
In this part, we analyze a diverse set of acquired solutions, evaluating their unique features. By investigating the physical structures, we interpret the dynamic nature of these solutions concerning the EIWE (
4). Employing the NAE method, this study generates a range of soliton structures. To deepen our comprehension of the dynamic wave behavior exhibited by the acquired solutions, we present graphical depictions in 3D, 2D, and density plots. Additionally, by employing suitable parameter values, we identify various soliton structures, including bright, kink, and dark solitons. In
Figure 1a, the solution
is visually presented in three dimensions in the
-plane, accompanied by its two-dimensional and density plots. The parameter values for this display are as follows:
,
,
,
,
,
,
,
,
, and
. The 3D graph spans the ranges
and
. The 2D plot is between dependent variable
and independent variable
x within the interval
at various time values as
. In
Figure 1b, 3D, density, and 2D visualizations of
are presented in the
-plane with the same parameters as in the previous case. The 3D graph is plotted within the ranges of
and
. The 2D plot illustrates the relationship between the dependent variable
and the independent variable
z within the interval
at
. In
Figure 1c, a visualization of the same solution is presented in the
-plane. The 3D graph is plotted within the ranges of
and
. The 2D plot illustrates the relationship between the dependent variable
and the independent variable
y within the interval
at
. In this context, the solutions manifest as bright solitons. In nonlinear optics, bright solitons are primary self-localized patterns found in optical fields within dispersive nonlinear mediums, such as optical fibers. They manifest as concentrated intensity peaks amid a uniform background.
Moving to
Figure 2a, the solution
is depicted in three dimensions in the
-plane, along with its 2D and density graphs. The parameter values for this display are as follows:
,
,
,
,
,
,
,
,
,
, and
. The 3D graph spans the ranges
and
. The 2D plot is between dependent variable
and independent variable
x within the interval
at
. In Figure
2b, the 3D, density, and 2D visualizations of
are presented in the
-plane with the same parameters as in the previous case. The 3D graph is plotted within the ranges of
and
. The 2D plot illustrates the relationship between the dependent variable
and the independent variable
t within the interval
at
. In Figure
2c, a visualization of same solution is presented in the
-plane. The 3D graph is plotted within the ranges of
and
. The 2D plot illustrates the relationship between the dependent variable
and the independent variable
z within the interval
at
. In this case, the solutions are displayed as kink solitons. In nonlinear optics, kink waves are propagating waves that rise or incline from one asymptotic position to another, eventually stabilizing as they extend to infinity.
Furthermore, in
Figure 3a, the solution
is depicted in three dimensions in the
-plane, along with its 2D and density graphs. The parameter values for this display are as follows:
,
,
,
,
,
,
,
,
,
, and
. The 3D graph spans the ranges
and
. The 2D plot is between dependent variable
and independent variable
x within the interval
at
. In
Figure 3b, the 3D, density, and 2D visualizations of
are presented in the
-plane with the same parameters as in the previous case. The 3D graph is plotted within the ranges of
and
. The 2D plot illustrates the relationship between the dependent variable
and the independent variable
z within the interval
at
. In
Figure 3c, a visualization of the same solution is presented in the
-plane. The 3D graph is plotted within the ranges of
and
. The 2D plot illustrates the relationship between the dependent variable
and the independent variable
x within the interval
at
. In this case, the solutions are displayed as dark solitons. Dark solitons attract considerable attention in optics due to their stable transmission. In fiber optics, a stable soliton typically emerges when there is a balance between Kerr nonlinearity and group velocity dispersion (GVD). Within the anomalous dispersion, solitons manifest as bright solitons, while in the normal dispersion, they are transmitted as dark solitons. All solutions exhibit practical applicability in optical communications. Furthermore, these solutions demonstrate the ability to maintain both their shape and velocity during prolonged travel over considerable distances.
4. Bifurcation Phenomena
In the current section, we analyze Equation (
4) within the framework of bifurcation analysis [
20,
21,
22]. Upon employing the Galilean transformation to Equation (
10), we derive the planar system as follows:
where
,
. The Hamiltonian of (
34) is expressed as follows:
which satisfies
In practice, Equation (
34) defines a two-dimensional system, and the trajectories within its phase space are governed by the vector field associated with Equation (
34). Hence, it is crucial to analyze different phase portraits of (
34) by varying its parameters. Equilibria and periodic orbits typically experience changes in their local stability characteristics during a bifurcation. For this reason, we initially pinpoint the critical points
(where j = 1, 2) for system (
35). These points of equilibrium occur when
and
, yielding two critical points determined by specific parameter values.
Furthermore, the Jacobian of the system will be
The determinant and trace of (
36) at the critical point
are indicated by
and
, respectively, and these values are provided as follows:
The fixed point is classified as a saddle when , a central point when and , a cusp when , and a node if and .
In the subsequent discussion, we depict symbols used for classifying distinct orbits within the phase portraits of a planner system (
34) as follows:
A nonlinear periodic trajectory is denoted as .
A nonlinear homoclinic trajectory is denoted as .
Here, and refer to the total fixed points and the overall separation layers encompassed within the trajectory depicted in the phase portrait. This gives rise to various scenarios as follows:
Case 1: Let and .
For
,
, and
, system (
34) has two fixed points,
and
. In this case,
is center and
is saddle. These points are presented in
Figure 4a.
Case 2: Let and .
For
,
, and
, system (
34) has two fixed points,
and
. In this case,
is saddle and
is center. These points are presented in
Figure 4b.
Case 3: Let and .
For
,
, and
, system (
34) has fixed points as
and
. In this case,
is saddle and
is center. These points are presented in
Figure 5a.
Case 4: Let and .
For
,
, and
, system (
34) has fixed points as
and
. In this case,
is center and
is saddle. These points are presented in
Figure 5b.
Case 5: Let and .
For
,
, and
, system (
34) has only one fixed point as
. In this case,
is a cuspidal point, and it is presented in
Figure 6a.
Case 6: Let and .
For
,
, and
, system (
34) has only one fixed point as
. In this case,
is a cuspidal point, and it is presented in
Figure 6b.
5. Chaotic Phenomena
In this section, we introduce an outward periodic force into system (
34) to examine the dynamics of both quasi-periodic and chaotic phenomena [
20,
21,
22]. The modified system can be expressed in the following manner:
The above system is autonomous with
. In the modified system (
37), the perturbing force is characterized by two components referred to as
and
. In this context,
represents the intensity, while
signifies the frequency of the perturbed term introduced into the system (
34). We have systematically examined and illustrated the chaotic dynamics of (
37) through diverse tools, including phase plots, time plots, and Lyapunov exponents. To delve deeper into this issue, we examine the impact of both the intensity
and the frequency
. We keep all other parameters constant at
and
during our analysis.
In
Figure 7, a two-dimensional plot and a time plot are presented, deliberately opting for small values for both intensity and frequency, specifically, assigning
and
. Under these chosen parameter values, periodic behavior is observed in system (
37). Subsequently, a slight increment is made to these parameters, setting them to
and
. During this adjustment, the system (
37) experiences a subtle disturbance and transitions into a quasi-periodic state, as visually depicted in
Figure 8. Continuing the analysis, the intensity is further elevated to
, and the frequency of the perturbing force is increased to
. In this scenario, system (
37) displays increased irregularity, revealing a quasi-periodic and chaotic nature, as illustrated in
Figure 9. Moving forward, the intensity and frequency are set to
and
. With these increased factors, system (
37) demonstrates heightened randomness, unveiling a chaotic nature, as presented in
Figure 10.
The Lyapunov exponent is employed to assess the rate at which nearby trajectories diverge within a dynamical system. It is a scalar metric that signifies the level of chaos within the system. A positive Lyapunov exponent denotes chaotic behavior, while a negative value indicates stability. A higher Lyapunov exponent value corresponds to a more chaotic system, showcasing a faster separation of neighboring trajectories. Subsequently, we visually represent the evolution of these exponents over time to gain insights into the dynamics of the system (
37). In
Figure 11, we illustrate the Lyapunov exponents obtained against the temporal evolution, aiming to identify the chaotic behavior of the modified dynamical system with the parameters
and initial condition (0.5, 0.5, 0.5).
6. Multistability
In the current part, our primary emphasis is on investigating the potential presence of multistability [
20] within (
37). When the system undergoes perturbation, it demonstrates a compelling ability to exhibit multistability. This suggests that it can display the simultaneous coexistence of various dynamic patterns, depending on specific parameter configurations and distinct initial conditions. These dynamic patterns encompass periodicity, quasi-periodicity, and chaos.
To visually depict this phenomenon of multistability, we refer to
Figure 12, which illustrates the diverse responses of the system to varying initial values. In
Figure 12a the red plot, corresponding to initial values (0.4, 0.3, 0), reveals a periodic nature. Meanwhile, the blue and green plots simultaneously demonstrate quasi-periodic behavior for initial values (−0.5, 0.1, 0) and (−0.03, 0.4, 0). In
Figure 12b, the red plot indicates periodic behavior for (−0.5, −0.03, 0), while the blue and green plots exhibit chaotic phenomena at (−0.4, −0.3, 0) and (−0.03, 0.4, 0).
7. Sensitivity Analysis
In this section, we investigate how the suggested equation responds to changes in initial conditions [
22,
23]. To evaluate the sensitivity of the model, we analyze three distinct sets of initial conditions. The first set, represented by the green curve, corresponds to
; the second set, illustrated by the red curve, corresponds to
; and the third set, indicated by the pink curve, involves
. In
Figure 13a, two solutions are observed: one for the initial condition
in green and the other for
in red. Likewise, in
Figure 13b, we present two solutions:
in green and
in pink.
Figure 13c also displays two solutions:
in red and
in pink.
We conducted a comparative analysis using different initial values, specifically (0.01, 0.01), (0.05, 0.03), and (0.4, 0.08), as shown in
Figure 13d. It is apparent that even slight modifications in the initial conditions can result in subtle shifts in the outcomes of a dynamic system. Therefore, we conclude that the proposed system demonstrates little sensitivity.
8. Conclusions
In summary, this study investigated the extended integrable wave equation, a widely utilized concept in physics, engineering, and diverse scientific disciplines. Our primary concern was to comprehensively analyze this equation, exploring various aspects including the extraction of Lie symmetries; solitons; an analysis of the phenomena of bifurcation; chaos; and an analysis of the sensitive nature of the suggested equation. Initially, the identification of the translational symmetries of the model was carried out through the utilization of the Lie approach. Then, the NAE method was applied to elucidate the soliton dynamics, and subsequent visualizations of the solutions were created utilizing Mathematica software, including three-dimensional, two-dimensional, and density plots.
Subsequently, an analysis of the dynamical nature of the model was conducted, encompassing various angles such as bifurcation, chaos, and sensitivity. The bifurcation phenomenon was analyzed at critical points within a dynamical system, accompanied by the application of an outward force, which unveiled the emergence of chaotic phenomena. Two-dimensional plots, time plots, multistability, and Lyapunov exponents were presented to illustrate these chaotic trajectories as depicted in
Figure 7,
Figure 8,
Figure 9,
Figure 10,
Figure 11 and
Figure 12. Furthermore, the sensitivity of the investigated model was executed utilizing the RK4 method. This analysis confirms that the stability of the solution is minimally affected by small changes in initial conditions. The outcomes of our investigation constitute a significant addition to the domain, offering new perspectives on the relevance and applicability of the extended water equation. Going beyond the confines of soliton theory, our findings contribute to a more comprehensive exploration of dynamical systems.