Dynamics and Soliton Propagation in a Modified Oskolkov Equation: Phase Plot Insights

Abstract

This study explores the modified Oskolkov equation, which depicts the behavior of the incompressible viscoelastic Kelvin–Voigt fluid. The primary focus of this research lies in several key areas. Firstly, the Lie symmetries of the considered equation are identified. These symmetries are utilized to transform the discussed model into an ordinary differential equation. Analytical solutions are subsequently derived using the new auxiliary equation technique. Next, a comprehensive analysis of the equation’s dynamic nature is undertaken from multiple aspects. Bifurcation is carried out at fixed points within the system, and chaotic behavior is unveiled by introducing an external force to the dynamic system. Various tools, including 3D and 2D phase plots, time series, Poincaré maps, and multistability analysis, are employed to identify the chaotic nature of the system. Furthermore, the sensitivity of the model is explored across diverse initial conditions. In general, comprehending the dynamic characteristics of systems holds immense significance in forecasting outcomes and innovating new technologies.

1. Introduction

In recent times, there has been a significant surge in interest surrounding nonlinear partial differential equations (NLPDEs). The intricate nonlinearity of natural phenomena is captivating, and numerous experts assert that nonlinear science presents the optimal avenue for comprehending the fundamental principles governing physical laws. A multitude of intricate physical phenomena and scientific domains, such as engineering [1], climate science [2], applied mathematics [3], biological science [4], and chemical reactions [5], find their description through the utilization of NLPDEs. In order to grasp the intricacies of the evaluated system, solving NLPDEs and deriving their numerical as well as analytical solutions assumes paramount importance. The analytical solutions stemming from these equations consistently prove invaluable in deciphering the natural world. A diverse array of efficient methodologies has been adopted to obtain analytical and numerical solutions for NLPDEs. Noteworthy among these are the bilinear residual network technique [6], trial equation approach [7], Lie symmetry [8], Bäcklund transformation [9], Riccati equation approach [10], q-Homotopy technique [11], and generalized auxiliary equation [12].

A wide array of solutions, such as rogue wave solutions, soliton solutions, multisoliton solutions, dromion wave solutions, lumps, and breather solutions, has been discovered by many scholars on different nonlinear models. Several to mention are the Korteweg–de Vries model [13], extended Kadomtsev–Petviashvili model [14], and Fokas–Lenells equation [15]. A particularly fascinating nonlinear occurrence that arises from mathematical modeling is the nonlinear Oskolkov equation [16], which develops a representation of the flow of incompressible viscoelastic fluids. The substantial role of viscoelastic models in applied fields and engineering, particularly in areas such as civil engineering and solid physics, has greatly inspired the impetus behind this research. The (1+1)-dimensional Oskolkov equation can be expressed as follows:

$${\Phi }_t-{\omega }_1{\Phi }_{xxt}-{\omega }_2{\Phi }_{xx}+\Phi {\Phi }_x=0.$$

(1)

The generalized form of Equation (1), also known as the Benjamin–Bona–Mahony–Burgers (BBMB) equation, is given by:

$${\Phi }_t-{\omega }_1{\Phi }_{xxt}-{\omega }_2{\Phi }_{xx}+\Phi {\Phi }_x+{\omega }_3{\Phi }_x=0.$$

(2)

This equation describes surface waves that propagate along the horizontal axis $Ox$. Equations (1) and (2) have been extensively examined in the literature over time to find different types of solutions. Various techniques have been employed in the works of different researchers. Roshid and Roshid [17] examined Equation (1) by adopting a simple equation technique and obtained various soliton patterns, Karakoc et al. [18] investigated Equation (2) by employing the finite element approach, Ghanbari [19] obtained traveling wave profiles for Equations (1) and (2) by using modified methodology, Kaplan et al. [20] implemented the exponential rational function approach and discussed the sensitivity of Equations (1) and (2), and Uddin et al. [21] analyzed the dynamic behavior and discovered traveling waves and rouge wave solutions of Equation (1) by adopting the unified technique.

In this study, we will rigorously direct our attention to the modified form of the Oskolkov model as

$${\Phi }_t-{\omega }_1{\Phi }_{xxt}-{\omega }_2{\Phi }_{xx}+{\Phi }^2{\Phi }_x=0,$$

(3)

where ${\omega }_1$ and ${\omega }_2$ are real parameters, representing the viscosity and the diffusive coefficient. Numerous studies have been carried out on Equation (3) by many scholars: Roshid et al. [22] discovered soliton solutions of Equation (3) by implementing the simple equation approach, and Akinfe [23] examined Equation (3) by using the hyperbolic tangent function technique. In this work, the Lie symmetry technique [24,25] is adopted to determine the symmetries of Equation (3). Then, by utilizing these symmetries, the discussed equation is translated into an ordinary differential equation (ODE). This manuscript introduces the application of the new auxiliary equation (NAE) technique [26] on a resulting ODE for the analysis of soliton patterns. It is noteworthy that the NAE method has not been previously employed to tackle Equation (3), as demonstrated in earlier studies. Through the utilization of this innovative technique, the outcomes are acquired in explicit form, expressed using trigonometric, hyperbolic, and rational functions. Additionally, the results are visually illustrated, thereby augmenting the clarity of the findings.

In recent years, there has been an increasing focus on the study of the behavior of differential equations (DEs) within the realm of chaos and bifurcation analysis. These fields of inquiry have become invaluable tools for comprehending complex phenomena governed by DEs. More recently, the utilization of bifurcation theory [27,28] in the examination of DEs has been an intriguing field of study. A bifurcation is a qualitative shift in the characteristics of a dynamical system brought on by altering parameters. Concurrently, chaos theory [29,30] provides insights into the chaotic characteristics of a given model, determining whether it exhibits chaotic behavior. This involves scrutinizing the stability or irregularity of solutions when nonlinear phenomena occurring within a medium are subjected to external forces. Numerous approaches for identifying chaos are at our disposal, namely:

• Phase portraits;
• Time series;
• Poincaré maps;
• Multistability.

As far as we know, there have been no prior investigations into the examination of chaos and bifurcation concerning the model being assessed. Furthermore, we have also examined the model’s sensitivity [31]. A system displays sensitive behavior even when a minor modification to the input results in a considerably distinct output.

In this study, we explore the modified Oskolkov model using valuable and effective methodologies. The Lie symmetry approach is employed to ascertain symmetry reductions; then, by utilizing these symmetries, the model is translated into an ODE. The resulting ODE was solved with the aid of the NAE technique for the analysis of soliton patterns. Subsequently, we have delved into the dynamics of the examined equation through the utilization of bifurcation and chaos theories. Consequently, we have observed phase plots depicting bifurcation characteristics, both quasi-periodic and chaos. Furthermore, we conducted a sensitivity analysis of the discussed model under different initial conditions.

In the subsequent sections of the manuscript, Section 2 delves into the Lie symmetry and symmetry reductions of the model being studied. Section 3 and Section 4 deal with analytical solutions along with their graphical representation. Section 5 presents an exploration of the dynamics of the examined equation, encompassing phase plots of bifurcation and chaotic behavior. This investigation is facilitated through both unperturbed and perturbed dynamical systems. In Section 6, we delve into the sensitive nature of the investigated equation. The Section 7 concludes with a comprehensive summary and discussion of the obtained results.

2. Lie Group Analysis and Symmetry Reductions

Lie symmetry analysis [32,33] involves the identification of transformations or symmetries that leave the given DE unchanged. These symmetries typically arise from infinitesimal operators, referred to as Lie operators, constituting Lie algebra. Employing these operators on the DE facilitates the identification of symmetries that preserve its structure. Identifying symmetries and transforming the PDE into lower dimensions or an ODE involves the following key steps.

Step 1:

Identify classical Lie point symmetries, specifically translational symmetries in our case, for the suggested model.

Step 2:

Construct an algebra (in our case, an Abelian algebra) based on the identified symmetries.

Step 3:

Determine the similarity variables corresponding to each symmetry.

Step 4:

Utilize the identified symmetries to reduce the PDE to an ODE.

Step 5:

Derive solutions in the form of traveling waves from the resulting ODE.

2.1. Lie Symmetries

Let us consider the one-parameter Lie group of infinitesimal transformations in $(x,t,\Phi )$ given by:

$$\begin{array}{c}\hfill x*=x+\epsilon {\sigma }^1(x,t,\Phi )+O\left({\epsilon }^2\right),\\ \hfill t*=t+\epsilon {\sigma }^2(x,t,\Phi )+O\left({\epsilon }^2\right),\\ \hfill \Phi *=\Phi +\epsilon \zeta (x,t,\Phi )+O\left({\epsilon }^2\right),\end{array}$$

with the group parameter $\epsilon \ll 1$, and ${\sigma }^1(x,t,\Phi )$, ${\sigma }^2(x,t,\Phi )$, and $\zeta (x,t,\Phi )$ are coefficient functions. The vector field related to the Lie group mentioned above may be expressed as:

$$\mathrm{R}={\sigma }^1(x,t,\Phi )\frac{\partial }{\partial x}+{\sigma }^2(x,t,\Phi )\frac{\partial }{\partial t}+\zeta (x,t,\Phi )\frac{\partial }{\partial \Phi },$$

where ${\sigma }^1(x,t,\Phi )$, ${\sigma }^2(x,t,\Phi )$, and $\zeta (x,t,\Phi )$ are to be found. The invariance condition for Equation (3) with $\mathrm{R}$ becomes:

$$P{r}^{\left(3\right)}\mathrm{R}\left({\Phi }_t-{\omega }_1{\Phi }_{xxt}-{\omega }_2{\Phi }_{xx}+{\Phi }^2{\Phi }_x\right){\mid }_{Equation3=0}=0.$$

(4)

Here, $P{r}^{\left(3\right)}\mathrm{R}$ is the third prolongation of $\mathrm{R}$, defined as:

$$P{r}^{\left(3\right)}\mathrm{R}=\mathrm{R}+{\zeta }^x\frac{\partial }{\partial {\Phi }_x}+{\zeta }^t\frac{\partial }{\partial {\Phi }_t}+{\zeta }^{xt}\frac{\partial }{\partial {\Phi }_{xt}}+{\zeta }^{xx}\frac{\partial }{\partial {\Phi }_{xx}}+{\zeta }^{xxt}\frac{\partial }{\partial {\Phi }_{xxt}},$$

(5)

where

$$\left\{\begin{array}{c}{\zeta }^x=D_x\left(\zeta \right)-{\Phi }_x{D}_x\left({\sigma }^1\right)-{\Phi }_t{D}_x\left({\sigma }^2\right),\hfill \\ {\zeta }^t=D_t\left(\zeta \right)-{\Phi }_x{D}_t\left({\sigma }^1\right)-{\Phi }_t{D}_t\left({\sigma }^2\right),\hfill \\ {\zeta }^{xx}=D_x\left({\zeta }^x\right)-{\Phi }_{xx}{D}_x\left({\sigma }^1\right)-{\Phi }_{xt}{D}_x\left({\sigma }^2\right),\hfill \\ {\zeta }^{xt}=D_x\left({\zeta }^t\right)-{\Phi }_{xx}{D}_x\left({\sigma }^1\right)-{\Phi }_{xt}{D}_x\left({\sigma }^2\right),\hfill \\ {\zeta }^{xxt}=D_t\left({\zeta }^{x}x\right)-{\Phi }_{xxx}{D}_t\left({\sigma }^1\right)-{\Phi }_{xxt}{D}_t\left({\sigma }^2\right).\hfill \end{array}\right.$$

(6)

Here, $D_x$ and $D_t$ are the total derivatives with respect to x and t, written as

$$D_j=\frac{\partial }{\partial x^j}+{\Phi }_j\frac{\partial }{\partial \Phi }+{\Phi }_{jk}\frac{\partial }{\partial {\Phi }_k}+\dots ,j=1,2.$$

By plugging Equations (5) and (6) into Equation (4), we yield the following system of determining equations:

$${\sigma }_x^2=0,{\sigma }_x^1=0,{\sigma }_t^2=0,{\sigma }_x^1=0,{\sigma }_{\Phi }^2=0,{\sigma }_{\Phi }^1=0,{\zeta }_{\Phi }=0,$$

(7)

where ${\sigma }_x=\frac{\partial \sigma }{\partial x}$, ${\sigma }_t=\frac{\partial \sigma }{\partial t}$, ${\zeta }_{\Phi }=\frac{\partial \zeta }{\partial \Phi }$.

By solving system (7), we obtain the symmetries for Equation (3) as follows:

$${\sigma }^1=B_{1}t+B_2,{\sigma }^2=B_{1}x+B_3,\zeta =0.$$

(8)

Here, $B_i$, $i=1,2,3$ are constants. Therefore, the infinitesimals of Equation (3) may be represented as:

$${\mathrm{R}}_1=\frac{\partial }{\partial t},{\mathrm{R}}_2=\frac{\partial }{\partial x}.$$

2.2. Symmetry Reductions

In the current section, we will reduce Equation (3) into ODEs by utilizing vector fields.

• Reduction using the symmetry ${\mathrm{R}}_1=\frac{\partial }{\partial t}$

The characteristic equation is:

$$\frac{dx}{0}=\frac{dt}{1}=\frac{d\Phi }{0}.$$

By integrating the aforementioned equation, we obtain

$$\Phi (x,t)=\Theta \left(\rho \right),\mathrm{where}\rho =x.$$

(9)

By plugging Equation (8) into Equation (3), we obtain

$$-{\omega }_2{\Theta }^{″}+{\Theta }^2{\Theta }^{\prime }=0.$$

(10)

Here, ${}^{\prime }$ denotes the derivative with respect to x, and it remains unaffected by changes in time.

• Reduction using the symmetry ${\mathrm{R}}_2=\frac{\partial }{\partial x}$

Under this case, the characteristic equation is:

$$\frac{dx}{1}=\frac{dt}{0}=\frac{d\Phi }{0}.$$

By integrating the above equation, we achieve

$$\Phi (x,t)=\Theta \left(\rho \right),\mathrm{where}\rho =t.$$

(11)

By plugging Equation (11) into Equation (3), we obtain

$${\Theta }^{\prime }=0.$$

(12)

Here, ${}^{\prime }$ denotes the derivative with respect to t, and it remains unaffected by the changes in position coordinate x. Equation (3) has the following solution:

$$\Phi (x,t)=C_1$$

where $C_1$ is the integration constant.

• Reduction using the symmetry ${\mathrm{R}}_1$ +$\mu {\mathrm{R}}_2 = \frac{\partial }{\partial t} + \mu \frac{\partial }{\partial x}$

The characteristic equation for this case is as follows:

$$\frac{dx}{\mu }=\frac{dt}{1}=\frac{d\Phi }{0}.$$

By integrating the above equation, we obtain

$$\Phi (x,t)=\Theta \left(\rho \right),\mathrm{where}\rho =x-\mu t.$$

(13)

By plugging Equation (13) into Equation (3), we obtain

$$\mu {\omega }_1{\Theta }^{″}-{\omega }_2{\Theta }^{\prime }-\mu \Theta +\frac{{\Theta }^3}{3}=0.$$

(14)

The next target is to retrieve the traveling wave structures for Equation (14) by adopting the new auxiliary equation method [30].

3. Analytical Solutions

Under the present section, we will obtain the soliton solutions of Equation (3) using a new auxiliary equation method. This technique has the following initial solution:

$$\Theta \left(\rho \right)=\sum _{i=0}^j{m}_i{\Gamma }^{i\psi \left(\rho \right)},$$

(15)

where $m_i$ is the arbitrary parameter, such that $m_j\ne 0$. The value of the positive integer j can be obtained using the principle of homogeneous balance, which involves balancing the highest order derivative with the nonlinear term in the ODEs. By equating the terms ${\Theta }^{″}$ and ${\Theta }^3$ in Equation (14) as $3j=j+2$, this results in $j=1$. Based on this, for $j=1$, we can take the starting solution as follows:

$$\Theta \left(\rho \right)=m_0+m_1{\Gamma }^{\psi \left(\rho \right)}.$$

(16)

Moreover, $\psi =\psi \left(\rho \right)$, satisfying the following ODE:

$${\psi }^{\prime }\left(\rho \right)=\frac{1}{ln(\Gamma )}\left(\lambda {\Gamma }^{-\psi \left(\rho \right)}+\varrho +\alpha {\Gamma }^{\psi \left(\rho \right)}\right).$$

(17)

Here, $\lambda$, $\varrho$, and $\alpha$ are real constants. By inserting Equations (16) and (17) into Equation (14) and reducing the coefficients of ${\Gamma }^{\psi \left(\rho \right)}$ to zero, we obtain a system of equations. After solving that system, for $m_0$, $m_1$, ${\omega }_1$, ${\omega }_2$, and $\mu$, by using software like Maple-13, the following solutions can be reported.

$$\begin{array}{c}\hfill m_0=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda },m_1=\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda },{\omega }_1=-\frac{2}{{\varrho }^2-4\alpha \lambda },{\omega }_2=0.\end{array}$$

By plugging the above values into Equation (16) and utilizing the transformation presented in Equation (13), the solutions of Equation (14) are as follows:

• Type 1: For $\alpha \ne 0$ and ${\varrho }^2-4\alpha \lambda <0$, we have

$$\left\{\begin{array}{c}{\Theta }_{1,1}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }\left(\frac{-\varrho }{2\alpha }+\frac{\sqrt{4\lambda \alpha -{\varrho }^2}}{2\alpha }tan\left(\frac{\sqrt{4\lambda \alpha -{\varrho }^2}}{2}\rho \right)\right),\hfill \\ \\ {\Theta }_{1,2}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }\left(\frac{-\varrho }{2\alpha }-\frac{\sqrt{4\lambda \alpha -{\varrho }^2}}{2\alpha }cot\left(\frac{\sqrt{4\lambda \alpha -{\varrho }^2}}{2}\rho \right)\right).\hfill \end{array}\right.$$

(18)

• Type 2: For $\alpha \ne 0$ and ${\varrho }^2-4\alpha \lambda >0$, we have

$$\left\{\begin{array}{c}{\Theta }_{2,1}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }\left(\frac{-\varrho }{2\alpha }-\frac{\sqrt{{\varrho }^2-4\lambda \alpha }}{2\alpha }tanh\left(\frac{\sqrt{{\varrho }^2-4\lambda \alpha }}{2}\rho \right)\right),\hfill \\ \\ {\Theta }_{2,2}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }\left(\frac{-\varrho }{2\alpha }-\frac{\sqrt{{\varrho }^2-4\lambda \alpha }}{2\alpha }coth\left(\frac{\sqrt{{\varrho }^2-4\lambda \alpha }}{2}\rho \right)\right).\hfill \end{array}\right.$$

(19)

• Type 3: For $\alpha \ne 0$, $\alpha =-\lambda$ and ${\varrho }^2+4{\lambda }^2<0$, we have

$$\left\{\begin{array}{c}{\Theta }_{3,1}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }{\left(\frac{\varrho }{2\lambda }-\frac{\sqrt{-{\varrho }^2-4{\lambda }^2}}{2\lambda }tan\left(\frac{\sqrt{-{\varrho }^2-4{\lambda }^2}}{2}\rho \right)\right)}^2,\hfill \\ \\ {\Theta }_{3,2}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }{\left(\frac{\varrho }{2\lambda }+\frac{\sqrt{-{\varrho }^2-4{\lambda }^2}}{2\lambda }cot\left(\frac{\sqrt{-{\varrho }^2-4{\lambda }^2}}{2}\rho \right)\right)}^2.\hfill \end{array}\right.$$

(20)

• Type 4: For $\alpha \ne 0$, $\alpha =-\lambda$ and ${\varrho }^2+4{\lambda }^2>0$, we have

$$\left\{\begin{array}{c}{\Theta }_{4,1}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }{\left(\frac{\varrho }{2\lambda }+\frac{\sqrt{{\varrho }^2+4{\lambda }^2}}{2\lambda }tanh\left(\frac{\sqrt{{\varrho }^2+4{\lambda }^2}}{2}\rho \right)\right)}^2,\hfill \\ \\ {\Theta }_{4,2}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }{\left(\frac{\varrho }{2\lambda }+\frac{\sqrt{{\varrho }^2+4{\lambda }^2}}{2\lambda }coth\left(\frac{\sqrt{{\varrho }^2+4{\lambda }^2}}{2}\rho \right)\right)}^2.\hfill \end{array}\right.$$

(21)

• Type 5: For $\alpha =\lambda$ and ${\varrho }^2-4{\lambda }^2<0$, we have

$$\left\{\begin{array}{c}{\Theta }_{5,1}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }\left(\frac{-\varrho }{2\lambda }+\frac{\sqrt{4{\lambda }^2-{\varrho }^2}}{2\lambda }tan\left(\frac{\sqrt{4{\lambda }^2-{\varrho }^2}}{2}\rho \right)\right),\hfill \\ \\ {\Theta }_{5,2}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }\left(\frac{-\varrho }{2\lambda }-\frac{\sqrt{4{\lambda }^2-{\varrho }^2}}{2\lambda }cot\left(\frac{\sqrt{4{\lambda }^2-{\varrho }^2}}{2}\rho \right)\right).\hfill \end{array}\right.$$

(22)

• Type 6: For $\alpha$=$\lambda$ and ${\varrho }^2-4{\lambda }^2>0$, we have

$$\left\{\begin{array}{c}{\Theta }_{6,1}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }\left(\frac{-\varrho }{2\lambda }-\frac{\sqrt{{\varrho }^2-4{\lambda }^2}}{2\lambda }tanh\left(\frac{\sqrt{{\varrho }^2-4{\lambda }^2}}{2}\rho \right)\right),\hfill \\ \\ {\Theta }_{6,2}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }\left(\frac{-\varrho }{2\lambda }-\frac{\sqrt{{\varrho }^2-4{\lambda }^2}}{2\lambda }coth\left(\frac{\sqrt{{\varrho }^2-4{\lambda }^2}}{2}\rho \right)\right).\hfill \end{array}\right.$$

(23)

• Type 7: For $\alpha \lambda <0$, $\varrho =0$, and $\alpha \ne 0$, we have

  $${\Theta }_{7,1}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }\left(-\sqrt{\frac{-\lambda }{\alpha }}tanh\left(\sqrt{-\lambda \alpha }\rho \right)\right).$$

  (24)

• Type 8: For $\lambda =-\alpha$, $\varrho =0$, we have

$${\Theta }_{8,1}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }\left(-\frac{1+e^{-2\alpha }}{1-e^{-2\alpha }}\rho \right).$$

(25)

• Type 9: For $\alpha =\varrho =K$, $\lambda =0$, we have

$${\Theta }_{9,1}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }\left(\frac{e^{K\rho }}{1-e^{K\rho }}\right).$$

(26)

• Type 10: For $\varrho =\lambda +\alpha$, we have

$${\Theta }_{10,1}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }\left(-\frac{1-\lambda e^{(\lambda -\alpha )\rho }}{1-e^{(\lambda -\alpha )\rho }}\right).$$

(27)

• Type 11: For $\varrho =-\lambda -\alpha$, we have

$${\Theta }_{11,1}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }\left(\frac{e^{(-\alpha +\lambda )\rho }-\lambda }{e^{(-\alpha +\lambda )\rho }-\alpha }\right).$$

(28)

• Type 12: For $\lambda =\varrho =\alpha \ne 0$, we have

$${\Theta }_{12,1}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }\left(\frac{1}{2}\left[\sqrt{3}tan\left(\frac{\sqrt{3}}{2}\lambda \rho \right)-1\right]\right).$$

(29)

• Type 13: For $\lambda =0$, we have

$${\Theta }_{13,1}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }\left(\frac{\varrho e^{\varrho \rho }}{1-\alpha e^{\varrho \rho }}\right).$$

(30)

• Type 14: For $\lambda =\alpha$ and $\varrho =0$, we have

$${\Theta }_{14,1}(x,t)=-\frac{\varrho \sqrt{3\mu }}{{\varrho }^2-4\alpha \lambda }+\frac{2\alpha \sqrt{3\mu ({\varrho }^2-4\alpha \lambda )}}{-{\varrho }^2+4\alpha \lambda }\left(tan\left(\lambda \rho \right)\right).$$

(31)

4. Graphical Illustration

The outcomes of the suggested model are visually presented in this section. Through the application of the new auxiliary equation technique, diverse soliton structures have been generated by varying the parameters for each result. The acquired data are depicted through three-dimensional, two-dimensional, and contour plots. Periodic, kink, and singular solitons have been documented for ${\Theta }_{1,1}(x,t)$, ${\Theta }_{2,2}(x,t)$, and ${\Theta }_{6,2}(x,t)$, respectively. The outcomes acquired from the experiments are exhibited in Figure 1, Figure 2 and Figure 3.

Figure 1 illustrates the identification of the periodic solitons solution for ${\Theta }_{1,1}(x,t)$, corresponding to the parameters $\varrho =2$, $\alpha =1$, $\lambda =0.1$, and $\mu =0.2$. Kink solitons are examined in Figure 2 for ${\Theta }_{2,2}(x,t)$, exploring various parameter values such as $\varrho =-2$, $\alpha =1.3$, $\lambda =0.25$, and $\mu =0.2$. Additionally, Figure 3 presents the singular solitons solution for ${\Theta }_{6,2}(x,t)$, with the selected parameters $\varrho =3$, $\alpha =-0.05$, $\lambda =-0.05$, and $\mu =0.2$. The kink solitons represent a noteworthy class that is capable of traversing nonlinear media, showcasing the advanced effects arising from dispersion and nonlinearity. Within the realm of fiber optics, a kink soliton refers to a shock front traversing the dispersive medium without undergoing distortion. Furthermore, each of these outcomes has the ability to maintain its form and velocity when traversing long distances and holds substantial importance in applied scientific domains.

5. Dynamics of the Investigating Model

The current section is used to depict the dynamical nature of the discussed model. For this, we will consider Equation (14) to discuss the behavior of bifurcation [34], chaotic dynamics, and multistability analysis through phase portraits [35].

5.1. Bifurcation Analysis

From Equation (14), we obtain the following planar dynamical system:

$$\left\{\begin{array}{c}\frac{d\Theta }{d\rho }=\mathrm{Y},\hfill \\ \frac{d\mathrm{Y}}{d\rho }={\kappa }_1\mathrm{Y}+{\kappa }_2\Theta -{\kappa }_3{\Theta }^3,\hfill \end{array}\right.$$

(32)

where ${\kappa }_1=\frac{{\omega }_2}{\mu {\omega }_1}$, ${\kappa }_2=\frac{1}{{\omega }_1}$, and ${\kappa }_3=\frac{1}{3\mu {\omega }_1}$. The above system has the integral as follows:

$$L(\Theta ,\mathrm{Y})=\frac{(1-{\kappa }_1)}{2}{\mathrm{Y}}^2-{\kappa }_2\frac{{\Theta }^2}{2}+{\kappa }_3\frac{{\Theta }^4}{4}=l,$$

(33)

where l is the real parameter. System (32) has the following fixed points on the $\Theta$-axis.

$$N_1=(0,0),N_2=(\sqrt{\frac{{\kappa }_2}{{\kappa }_3}},0),N_3=(-\sqrt{\frac{{\kappa }_2}{{\kappa }_3}},0).$$

Furthermore, the Jacobian of (32) is:

$$J(\Theta ,\mathrm{Y})=\left|\begin{array}{cc}0& 1\\ {\kappa }_2-3{\kappa }_3{\Theta }^2& {\kappa }_1\end{array}\right|$$

$$=-{\kappa }_2+3{\kappa }_3{\Theta }^2,$$

Thus $(\Theta ,\mathrm{Y})$ is a saddle point for $J(\Theta ,\mathrm{Y})<0$, a center for $J(\Theta ,\mathrm{Y})>0$, and a cusp if $J(\Theta ,\mathrm{Y})=0$.

In order to examine the nature of system (32) at critical points, the following outcomes are observed.

• Case 1: Let ${\kappa }_2>0$ and ${\kappa }_3>0$.

For ${\kappa }_2=1$, ${\kappa }_3=1$, ${\omega }_1=1$, and $\mu =0.3$, system (32) exhibits three fixed points: $N_1=(0,0)$, $N_2=(1,0)$, and $N_3=(-1,0)$. In this case, $N_1$ is the saddle point and $N_2,N_3$ are central points. Additionally, the impact of ${\kappa }_1\mathrm{Y}$ has been discussed by taking different values of ${\kappa }_1$. The phase plots are displayed in Figure 4.

• Case 2: Let ${\kappa }_2<0$ and ${\kappa }_3<0$.

For ${\kappa }_2=-1$, ${\kappa }_3=-1$, ${\omega }_1=-1$, and $\mu =0.3$, the system (32) exhibits three fixed points: $N_1=(0,0)$, $N_2=(1,0)$, and $N_3=(-1,0)$. In this case, $N_1$ is the center and $N_2,N_3$ are saddle points. Additionally, the impact of ${\kappa }_1\mathrm{Y}$ had been discussed by taking different values of ${\kappa }_1$. The phase plots are displayed in Figure 5.

5.2. Behavior of Chaotic Motion

This section examines the nature of quasi-periodic and chaotic dynamics of considered model (3). For this, we incorporate a periodic force into the system (32). Hence, the modified form of the system (32) is:

$$\left\{\begin{array}{c}\frac{d\Theta }{d\rho }=\mathrm{Y},\hfill \\ \frac{d\mathrm{Y}}{d\rho }={\kappa }_1\mathrm{Y}+{\kappa }_2\Theta -{\kappa }_3{\Theta }^3+{\delta }_{0}cos\left(G\right),\hfill \\ \frac{dG}{d\rho }=\upsilon .\hfill \end{array}\right.$$

(34)

The aforementioned system is autonomous with $G=\upsilon \rho$. In the system (34), the external perturbed force consists of two elements, denoted as ${\delta }_0$ and $\upsilon$. Here, ${\delta }_0$ and $\upsilon$ indicate the amplitude and frequency, respectively, of the perturbed term. An unpredictable and random behavior can be observed by adding an eternal force to the planar system. In our work, we have analyzed this phenomenon within a system (32) that demonstrates a chaotic nature, characterized by trajectories diverging from regular patterns and displaying randomness with time. We have investigated the chaotic dynamics of the system (34) using various tools, such as phase plots, time plots, and Poincaré maps. In order to examine the behavior, we will evaluate the influence of the parameters ${\delta }_0$ and $\upsilon$.

In Figure 6, the three-dimensional and two-dimensional profiles, time plot, and Poincaré map are shown for the parametric values ${\kappa }_1=0$, ${\kappa }_2=0.02$, ${\kappa }_3=2$, ${\delta }_0=0.1$, and $\upsilon =2$. It is noted that system (34) is quasi-periodic. In Figure 7, by increasing the amplitude and frequency as ${\delta }_0=1.1$ and $\upsilon =3$ while the other parameters remain the same as in the last case, system (34) exhibits irregular patterns. It is evident that the disturbed system demonstrates a chaotic nature under these parameter values. Additionally, the Poincaré map shows numerous scattered points, providing further confirmation of the chaotic nature of the system.

5.3. Multistability Analysis

When a perturbation is introduced to the planar system, it can demonstrate multistability, implying the presence of several possible dynamic behaviors. These behaviors manifest for the same parameter values but distinct initial circumstances. This spectrum of behaviors encompasses chaos, quasi-periodicity, and periodicity, all of which may emerge within the system under varying conditions. In Figure 8, we conducted an analysis of multistability concerning the disturbed system (34) under diverse initial circumstances. Our observations indicate a pronounced sensitivity of system (34) to the initial conditions that result in a chaotic nature.

6. Sensitivity Analysis

In this segment, our motivation was to analyze the sensitivity of the system (32). To achieve this, we took into account three distinct initial conditions: $(\Theta ,\mathrm{Y})=(0.4,0)$ in solid green, $(\Theta ,\mathrm{Y})=(0.6,0)$ in red (dashed), and $(\Theta ,\mathrm{Y})=(0.8,0)$ in solid blue. In Figure 9, two solutions are depicted: $(\Theta ,\mathrm{Y})=(0.4,0)$ in green (solid) and $(\Theta ,\mathrm{Y})=(0.6,0)$ in red (dashed). In Figure 10, two solutions are presented: $(\Theta ,\mathrm{Y})=(0.4,0)$ in green (solid) and $(\Theta ,\mathrm{Y})=(0.8,0)$ in blue (solid). In Figure 11, two solutions are displayed: $(\Theta ,\mathrm{Y})=(0.6,0)$ in red (dashed) and $(\Theta ,\mathrm{Y})=(0.8,0)$ in solid blue. A comparison was conducted using various initial values, specifically (0.4, 0), (0.6, 0), and (0.8, 0), as depicted in Figure 12. The outcomes distinctly indicate that minor differences in the initial values significantly influence the final outcome. Hence, it can be deduced that the examined model has a sensitive nature.

7. Concluding Remarks

In this work, we successfully investigated the modified Oskolkov equation, which elucidates the behavior of the incompressible viscoelastic Kelvin–Voigt fluid. Through our analysis, we examined this equation from multiple perspectives, including the Lie symmetry technique, bifurcation analysis, and chaotic behaviors, as well as the sensitive nature of the considered model. First of all, we identified the Lie symmetries of the equation and, utilizing similarity reduction, we changed the partial differential equation into an ordinary differential equation. Solving the obtained ordinary equation with the assistance of the new auxiliary technique, we successfully derived soliton solutions, which are visually depicted in Figure 1, Figure 2 and Figure 3. Furthermore, by performing some calculations, we derived a second-order ordinary equation, yielding a system of first-order DEs (32).

By subjecting the considered model to a bifurcation analysis, we have comprehensively explored its dynamic nature at critical points. This culminated in the depiction of phase portraits, visually represented in Figure 4 and Figure 5. Furthermore, we introduced an external periodic force into the dynamical system, subsequently employing a range of tools to identify chaotic behavior. The utilization of phase plots, Poincaré maps, time plots, and multistability analysis has effectively revealed the presence of quasi-periodic and chaotic dynamics within the disturbed system. These findings are depicted in Figure 6, Figure 7 and Figure 8. In addition, our investigation has highlighted the system’s sensitivity to the initial conditions. Remarkably, even minor adjustments in the initial values lead to considerable variations in the system’s behavior. In conclusion, the outcomes of this study are not only captivating but also underscore the efficacy of the proposed methodologies in assessing soliton dynamics and phase patterns across diverse nonlinear models.