Multi-scale dynamics of predator-prey systems with Holling-IV functional response

1.   Introduction

In many physical fields, nonlinear phenomena can be described by nonlinear Schrödinger equation (NLSE). For example, in plasma physics, the coupled Schrödinger-Korteweg-de Vries equation can effectively describe various processes in dusty plasma, such as Langmuir waves, dust-acoustic waves, and electromagnetic waves ^[1]; in fluid mechanics, NLSE can describe nonlinear wave phenomena ^[2]. In addition, NLSE is also an important model for describing Bose-Einstein condensation ^[3,4].

In particular, in the field of optical fiber communication, the NLSE can describe the evolution of optical solitons in optical fibers ^[5]. Optical solitons are the result of the combination of group velocity dispersion (GVD) and self-phase modulation (SPM) and are the exact balance of GVD and SPM effects in the anomalously dispersive region ^[6]. Due to the nonlinear effect of optical fibers, strong self-interaction and self-modulation effects occur when the optical pulse propagates in the medium, resulting in changes in the shape, frequency, and phase of the optical pulse. These changes can be described and analyzed by the NLSE. Therefore, the NLSE is an important model to describe the propagation of optical solitons ^[7,8,9].

In fact, the propagation modes of optical solitons are different in different nonlinear mediums. For example, many scholars have studied the propagation of optical solitons based on Kerr's law ^{[10,11,12,13]}. With the development of research, scholars began to be interested in optical solitons under non-Kerr's law, such as Kudryashov's law ^[14,15], parabolic law ^[16,17], power law ^[18,19], etc. Therefore, the standard NLSE has been extended to coupled, high-dimensional, high-order forms or other special forms that satisfy various specific conditions. In more detail, for example, the Manakov model, which is an extended form of Schrödinger equation, is a coupled partial differential equation for controlling dual-mode transmission in fiber optic communications and deep-sea transport with double-layer constraints ^[20]; the cascaded model is a model combining NLSE, Lakshmanan-Porsezian-Daniel and Sasa-Satsuma equation, which is used to describe the propagation of optical solitons over long-distance ^[21]; the anti-cubic nonlinear form of NLSE can describe the propagation of pulses in optical metamaterials with anti-cubic nonlinearity ^[22]; the nonlinear Schrödinger-Hirota equation is derived from the NLSE by the Lie transform, which can accurately describe the propagation of dispersive optical solitons when GVD is small ^[23].

From a mathematical point of view, optical solitons can be regarded as completely integrable solutions of the NLSE. By controlling the parameters of soliton solutions, the modulation and transmission of optical signals can be realized. Therefore, it is very important to study the exact solutions of NLSE. There are many effective methods to solve NLSE, such as the improved extended tanh function method ^[24,25], the Jacobi elliptic function expansion method ^[26,27], and enhanced Kudryashov's approach ^[28], etc.

In this paper, we study the extended (3+1)-dimensional NLSE proposed by Elsherbeny et al. ^[29], which is built on the SPM form of Kudryashov. It extends the NLSE to (3+1)-dimensions by means of a comprehensive dispersion model and introduces a new SPM effect. The equation is as follows:

where the coefficients $a_{i} (i = 1, 2, 3, 4, 5, 6)$ and $b_{j} (j = 1, 2)$ are constants, which are parameters that characterize the propagation pulse. $n$ is an arbitrary parameter related to the pulse propagation characteristics in the fiber. $i = \sqrt{-1}$, $t$ represents the temporal variable, while $x, y, z$ represent the spatial variables. Lastly, $a_{2}, a_{3}, a_{4}, a_{5}$ and $a_{6}$ come from the cross-spatial dispersion. When $a_{2} = a_{3} = a_{4} = a_{5} = a_{6} = 0$, Eq (1.1) obtains Kudryashov's equation. If $b_{1} = b_{2} = 0$, Eq (1.1) demonstrates the dual-power law nonlinearity, and then $n = 2$, Eq (1.1) exhibits the parabolic law nonlinearity.

The Eq (1.1) describes the propagation of optical fiber under the condition of interaction between SPM and spatial dispersion effect, which is of great help to understanding the propagation of optical in complex media and plays an important role in the field of optical fiber communication. In particular, the study of the optical solitons problem of this equation is of great significance to the development of nonlinear optics and to the solution of practical problems such as optical signal transmission. Therefore, it is an important work to obtain the exact solutions of Eq (1.1).

In recent years, many scholars have studied the complex, high-dimensional NLSE. It should be mentioned that Elsherbeny et al. used the enhanced Kudryashov's approach for the Eq (1.1) to retrieve the optical bullets and domain walls in the presence of cross-spatial dispersion effects ^[29]. Biswas et al. used the variational principle to obtain the parametric dynamics of the multidimensional solitons of a generalized (3+1)-dimensional NLSE ^[30]. Hosseini et al. used the new expansion method based on the Jacobi elliptic equation to study the (3+1)-dimensional resonance NLSE, the which obtained the exact solutions of several Jacobi elliptic functions and exponential functions ^[31]. Kumar et al. studied the (3+1)-dimensional NLSE with parabolic law, which used the extended generalized Riccati equation mapping method to show the dynamic behavior of isolated waves ^[32]. Rabie et al. used the modified extended mapping approach to study NLSE, for which rational solutions, exponential solutions, hyperbolic wave solutions, singular periodic solutions, and Jacobi elliptic function solutions are given ^[33].

In the research process of the above methods, the assumptions of the form of solutions are often needed, and the solving process has great limitations on parameters. However, the trial equation method and the complete discrimination for polynomial method by Liu are systematic and efficient, which can directly derive the solutions of nonlinear differential equations (NDE) without such assumptions ^{[34,35,36,37]}. In particular, compared with other existing methods, the method proposed by Liu can not only obtain richer solution forms but can also be directly written out when given parameters. In addition, this method can also be used to analyze the properties of dynamic systems ^[38,39]. Now many scholars have used these methods to solve many important NDEs ^{[40,41,42,43]}.

In this paper, we apply Liu's method to solve Eq (1.1) and obtain multiple types of solutions, such as solitary wave solutions, trigonometric function solutions, singular solutions, etc. In particular, the elliptic function double periodic solutions of this equation are also obtained, which is a new discovery. In addition, we also apply this method to study the dynamic characteristics of Eq (1.1) and obtain the existence conditions of solitons and periodic solutions of the equation. In particular, we analyze the chaotic behaviors of Eq (3.4) and find that given appropriate perturbations, such as Gaussian perturbation or triangular perturbation, the system has chaotic behaviors, which is a meaningful topic ^[44,45].

The framework of this paper is as follows: In Section 2, we transform Eq (1.1) into an ordinary differential equation (ODE) by the traveling wave transformation and the trial equation method. In Section 3, we write the dynamic system of Eq (1.1) and analyze the dynamic properties by the complete discrimination system for the polynomial method, which obtains the existence conditions of soliton solutions and periodic solutions. In Section 4, we obtain the exact solutions by the complete discrimination system for the polynomial method and plot 3-dimensional phase diagrams of the solutions to show the physical characteristics. In Section 5, for the perturbed system given appropriate perturbation terms, we verify the chaotic behaviors of the system through the largest Lyapunov exponents (LLEs) and the phase diagrams. Finally, in Section 6, we make a summary.

2.   Mathematical analysis

To obtain the solutions to Eq (1.1), we use the following transformation ^[29],

where

Along the $x, y$, and $z$ directions, the wave numbers are represented by $k_{j}$ for $j = 1, 2, 3$, respectively. Additionally, $\theta$ represents the phase constant, and $\omega$ represents the frequency. Substituting Eqs (2.1) and (2.2) into Eq (1.1) and splitting imaginary and real parts, we obtain

and

where

Applying the restriction $\Gamma(\xi) = {V(\xi)}^{\frac{1}{n}}$ to Eq (2.4), then we obtain the ODE as follows:

For Eq (2.7), it cannot be solved directly by integral, so we choose the trial equation method to solve it.

First, we set the trial equation ^[34],

We take the derivative of Eq (2.8) and obtain

Bringing Eqs (2.8) and (2.9) into Eq (2.7), then balancing the order of $n^{2}b^{4}V^{4}$ and $-nAk^{2}VV^{\prime\prime}$, we obtain $n_{1} = 4$, then

Bringing Eqs (2.10) and (2.11) into Eq (2.7), let each coefficient of $V, V^{2}, V^{3}$, and $V^{4}$ be zero to form an algebraic equation system. We can obtain the values of the coefficient $z_{i}, (i = 0, 1, 2, 3, 4)$ of Eq (2.10) as follows:

provided $n\neq 1, n\neq 2.$

3.   Qualitative analysis

In this section, the dynamic system of the equation is written, and then the Hamiltonian is obtained. According to the complete discrimination system for the polynomial method, the existence conditions of soliton solutions and periodic solutions of Eq (1.1) are obtained.

Substituting the following transformation into Eq (2.10)

we have

where

The dynamic system of Eq (3.2) as follows:

where $\lambda$ represents the generalized momentum, while $\eta$ represents the generalized coordinate. Then the Hamiltonian is

This is an autonomous conservation system, which satisfies

Therefore, the trajectory of the dynamical system is the contour of the Hamiltonian (3.5). In addition, we can get the following potential energy

Since the dynamic characteristics of Eq (3.2) can be analyzed through the equilibrium point, it is necessary to know the specific equilibrium point. The equilibrium point is the calculated root of the following equation:

We discuss the case of roots according to the complete discrimination system for the cubic polynomial method. The complete discrimination system of cubic polynomials is as follows ^[34]:

we will discuss it in the following four cases.

Case 1. $\triangle = 0, p < 0$, then,

The dynamical system (3.2) has two equilibrium points $(\alpha, 0), (\beta, 0)$. The discussion is divided into two scenarios. When $z_{4} < 0$, $(\beta, 0)$ is a center, and $(\alpha, 0)$ is a cusp. For example, when $z_{4} = -1$, $p = -9, q = 6\sqrt{3}$, we obtain $\alpha = \sqrt{3}, \beta = -2\sqrt{3}$. The phase diagram of this example is shown in Figure 1. From trajectories Ⅰ, we can see that $(-2\sqrt{3}, 0)$ is the center, and from trajectories Ⅱ, $(\sqrt{3}, 0)$ is the cusp. At the same time, it can be seen that trajectories Ⅰ, Ⅱ, and Ⅲ are closed orbits, so the Eq (3.2) has periodic solutions. Analogously, while $z_{4} > 0$, Figure 2 is obtained, which can be seen from the red and blue trajectories we obtain: $(\beta, 0)$ is a cusp and $(\alpha, 0)$ is a saddle point.

Case 2. $\triangle = 0, p = 0$, then,

The dynamical system (3.2) has one equilibrium point $(\alpha, 0)$. When $z_{4} < 0$, $(\alpha, 0)$ is a center. For example, when $z_{4} = -1, p = 0, q = 0$, we get $\alpha = 0$. The phase diagram of this example is shown in Figure 3, from which it can be seen that trajectories Ⅰ, Ⅱ, and Ⅲ are closed orbits centered on $(0, 0)$, so the Eq (3.2) has periodic solutions. Analogously, while $z_{4} > 0$, we can see from Figure 4 that $(0, 0)$ is the center.

Case 3. $\triangle > 0, p < 0$, then

The dynamical system (3.2) has three equilibrium points $(\alpha, 0), (\beta, 0), (\gamma, 0)$. When $z_{4} < 0$, $(\alpha, 0)$ and $(\gamma, 0)$ are centers, and $(\beta, 0)$ is a saddle point. For example, when $z_{4} = -1, p = -8, q = 0$, we get $\alpha = 2\sqrt{2}, \beta = 0, \gamma = -2\sqrt{2}$. At this point, that is, $\beta = 0$, we find that the phase diagram is symmetric. It can be seen from Figure 5 that there are two homologous orbitals, namely trajectories Ⅲ and Ⅳ, which present the "Figure-eight loop". This shows that Eq (3.2) has bell-shaped solitons (bright and dark solitons) solutions. While $z_{4} = 1, p = 8, q = 0$, then $\alpha = 2\sqrt{2}, \beta = 0, \gamma = -2\sqrt{2}$, the phase diagram for this example is shown in Figure 6. It can be seen that the trajectories Ⅰ and Ⅱ are centered on $(-2\sqrt{2}, 0), (2\sqrt{2}, 0)$, and trajectories Ⅱ and Ⅲ are two heterologous trajectories, which shows that equations Eq (3.2) have kink and anti-kink solitaire wave solutions.

Case 4. $\triangle < 0$, then

The dynamical symstem (3.2) has one equilibrium point $(\alpha, 0)$. When $z_{4} < 0$, $(\alpha, 0)$ is a center. This case is similar to Case $2$ in this section. Therefore, it is omitted here.

In conclusion, according to the complete discrimination system for the polynomial method, the existence conditions of soliton solutions and periodic solutions are obtained. The next step is to construct a specific traveling wave solution of Eq (3.2).

4.   Quantitative analysis

In this section, we divided the solution cases into nine categories, and thirteen exact solutions of the equation are obtained by the complete discrimination system for the quartic polynomial method, among which the elliptic function double periodic solutions are new solutions, and the existence of the solutions is proved by examples.

4.1.   Exact solutions

Here, since the state analysis processes for $z_{4} > 0$ and $z_{4} < 0$ are similar, we will only consider the state when $z_{4} > 0$. Taking the following transformation into Eq (2.10):

we obtain

where

Simplifying Eq (4.2) into the integral form

where

Now, we can use the complete discrimination system for quartic polynomials to solve Eq (4.4). The discrimination system is as follows ^[34]:

Then we obtain thirteen different modes of optical wave in total.

Case 1. $D_{2} = D_{3} = D_{4} = 0$, then $F(\phi) = \phi^{4}$, we have a singular rational solution as follows:

Case 2. $E_{2} < 0, D_{2} < 0, D_{3} = D_{4} = 0$, then $F(\phi) = [(\phi-\gamma_{1})^{2}+\gamma_{2}^{2}]^{2}$, we have a singular periodic solution as follows:

where $\gamma_{1}$ and $\gamma_{2}$ are real constants.

Case 3. $D_{2} > 0, D_{3} > 0, D_{4} = 0$, then $F(\phi) = (\phi-\gamma_{1})^{2}(\phi-\gamma_{2})(\phi-\gamma_{3})$, when $\frac{\gamma_{1}-\gamma_{3}}{\gamma_{2}-\gamma_{1}} < 0$, and $\gamma_{1}$, $\gamma_{2}$, and $\gamma_{3}$ are real constants, we have two solitary wave solutions as follows:

and

where $B = ((\frac{n{b_{4}}}{(1+n)Ak^{2}})^{\frac{1}{4}}\xi-\xi_{0})\sqrt{(\gamma_{1}-\gamma_{2})(\gamma_{1}-\gamma_{3})}$, $\gamma_{1}$, $\gamma_{2}$ and $\gamma_{3}$ are real constants.

When $\frac{\gamma_{1}-\gamma_{3}}{\gamma_{2}-\gamma_{1}} > 0$, and $\gamma_{1}$, $\gamma_{2}$, and $\gamma_{3}$ are real constants. Then we have a solitary wave solution as follows:

where $B = ((\frac{n{b_{4}}}{(1+n)Ak^{2}})^{\frac{1}{4}}\xi-\xi_{0})\sqrt{(\gamma_{2}-\gamma_{1})(\gamma_{1}-\gamma_{3})}$.

Case 4. $E_{2} > 0, D_{2} > 0, D_{3} = D_{4} = 0$, then $F(\phi) = (\phi-\gamma_{1})^{2}(\phi-\gamma_{2})^{2}$, where $\gamma_{1}$ and $\gamma_{2}$ are real numbers and $\gamma_{1} > \gamma_{2}$.

When $\phi > \gamma_{1}$ or $\phi < \gamma_{2}$, we obtain the solitary wave solution as follows:

when $\gamma_{2} < \phi < \gamma_{1}$, we obtain the solitary wave solution as follows:

Case 5. $D_{3} < 0, D_{4} = 0$, then $F(\phi) = (\phi-\gamma_{1})^{2}[(\phi-\gamma_{2})^{2}+\gamma_{3}^{2}]$, we have a solitary wave solution in exponential form as follows:

where $c = \exp[\pm\sqrt{(\gamma_{1}-\gamma_{2})^{2}+\gamma_{3}^{2}}((\frac{n{b_{4}}}{(1+n)Ak^{2}})^{\frac{1}{4}}\xi-\xi_{0})]$, $l = \frac{\gamma_{1}-2\gamma_{2}}{\sqrt{(\gamma_{1}-\gamma_{2})^{2}+\gamma_{3}^{2}}}$, $\gamma_{1}$, $\gamma_{2}$ and $\gamma_{3}$ are real constants.

Case 6. $E_{2} = 0, D_{2} > 0, D_{3} = D_{4} = 0$, then $F(\phi) = (\phi-\gamma_{1})^{3}(\phi-\gamma_{2})$, we have a rational singular mode as follows:

where $\gamma_{1}$ and $\gamma_{2}$ are real constants.

Case 7. $D_{2} > 0, D_{3} > 0, D_{4} > 0$, then $F(\phi) = (\phi-\gamma_{1})(\phi-\gamma_{2})(\phi-\gamma_{3})(\phi-\gamma_{4})$, we have the elliptic function double periodic solutions as follows:

where $C = \frac{\sqrt{(\gamma_{1}-\gamma_{3})(\gamma_{2}-\gamma_{4})}}{2}((\frac{n{b_{4}}}{(1+n)Ak^{2}})^{\frac{1}{4}}\xi-\xi_{0})$, $M = \sqrt{\frac{(\gamma_{1}-\gamma_{4})(\gamma_{2}-\gamma_{3})}{(\gamma_{1}-\gamma_{3})(\gamma_{2}-\gamma_{4}})}$, $\gamma_{1}$, $\gamma_{2}$, $\gamma_{3}$ and $\gamma_{4}$ are real constants and $\gamma_{1} > \gamma_{2} > \gamma_{3} > \gamma_{4}$.

Case 8. $D_{4} > 0, D_{2}D_{3}\leq0$, then $F(\phi) = [(\phi-\zeta_{1})^{2}+\zeta_{2}^{2}][(\phi-\zeta_{3})^{2}+\zeta_{4}^{2}]$. When $\zeta_{2}\geq\zeta_{4} > 0$, we have the elliptic function double periodic solutions as follows:

where $\zeta_{1}$, $\zeta_{2}$, $\zeta_{3}$, and $\zeta_{4}$ are real constants, and

Case 9. $D_{4} < 0, D_{2}D_{3}\geq0$, then $F(\phi) = (\phi-\zeta_{1})(\phi-\zeta_{2})[(\phi-\zeta_{3})^{2}+\zeta_{4}^{2}]$, we have the elliptic function double periodic solutions as follows:

where $\zeta_{1}$, $\zeta_{2}$, $\zeta_{3}$ and $\zeta_{4}$ are real constants, and

4.2.   Numerical simulation of exact solutions

In this section, we calculate several characteristic solutions under the conditions of given parameter values and plot envelope graphs to obtain a more intuitive conclusion.

Example 1. Triangular function periodic solutions

Taking $n = A = B = a_{3} = 3$, $a_{1} = a_{2} = a_{4} = a_{5} = a_{6} = k = k_{1} = k_{2} = k_{3} = \theta = \gamma_{2} = \xi_{0} = B_{3} = 1$, $B_{1} = B_{2} = \gamma_{1} = 0, w = \frac{53}{24}, v = 2$, $b_{4} = \frac{3}{4}, b_{3} = \frac{5}{6}, b_{2} = -\frac{17}{96}, b_{1} = -\frac{289}{384}$, we get

The 3D phase diagram of $|u_{2}|$ is shown in Figure 7.

Example 2. Rational solutions

Taking $n = A = B = a_{3} = 3, a_{1} = a_{2} = a_{4} = a_{5} = a_{6} = k = B_{3} = k_{1} = k_{2} = k_{3} = \theta = \xi_{0} = 1$, $B_{1} = B_{2} = \gamma_{1} = \gamma_{2} = \gamma_{3} = \gamma_{4} = 0$, $w = \frac{23}{8}, v = 2, b_{4} = \frac{3}{4}, b_{3} = \frac{5}{6}, b_{2} = -\frac{1}{96}, b_{1} = -\frac{1}{384}$, we get

The 3D phase diagram of $|u_{1}|$ is shown in Figure 8.

Example 3. Solitary wave solutions

Taking $n = A = B = a_{3} = 3, \gamma_{1} = a_{1} = a_{2} = a_{4} = a_{5} = a_{6} = k = B_{3} = k_{1} = k_{2} = k_{3} = \theta = \xi_{0} =$$1$, $\gamma_{2} = -1, B_{1} = B_{2} = 0, w = \frac{85}{24}, v = 2, b_{4} = \frac{3}{4}$, $b_{3} = \frac{5}{6}, b_{2} = -\frac{5}{32}, b_{1} = -\frac{75}{128}$, we get

The 3D phase diagram of $|u_{6}|$ is shown in Figure 9.

Example 4. Solitary wave solutions

Taking $n = k_{3} = 3, \gamma_{1} = a_{1} = a_{2} = a_{4} = a_{5} = a_{6} = k = B_{3} = k_{1} = k_{2} = \theta = 1$, $\gamma_{2} = -1, B_{1} = B_{2} = \xi_{0} = 0, w = \frac{1406}{54}, v = 1, b_{4} = -\frac{64}{9}$, $b_{3} = -\frac{5}{9}, b_{2} = \frac{8}{9}, b_{1} = \frac{21}{9}, A = a_{3} = \frac{1}{6}, B = 27$, we get

The 3D phase diagram of $|u_{7}|$ is shown in Figure 10.

5.   Chaotic behaviors

The study of chaotic phenomena in models is of great significance to physics. In the previous sections, we analyze Eq (1.1) qualitatively and quantitatively, and no chaotic phenomenon is found. However, we found that by adding a specific type of perturbation term, the dynamical system can appear to have chaotic behaviors, which is a new result. The specific analysis process is as follows:

First, consider dynamic system (3.2) with perturbation term

where $H(\xi)$ is the perturbed term and $s = 2z_{4}, m = Q_{2}, e = \frac{1}{2}Q_{1}$.

Secondly, the cases where the perturbation terms are trigonometric function and Gaussian function are discussed, respectively.

Case 1. Gaussian function: $H(\xi) = 100\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}(0.025\xi)^{2}}$.

Letting $s = -30, m = -20, e = -5$, we respectively show the LLEs corresponding to the three parameters, as shown in Figures 11–13. It can be seen from the figures that the parameters corresponding to the cubic term have the greatest influence on the chaotic behaviors of the perturbed system. Meanwhile, the corresponding three-dimensional phase diagram is shown in Figure 14. In addition, it can be seen from Figure 14 that the trajectories of the perturbed system (5.1) intersect, indicating that the solution is not unique if the initial value is given at the intersection.

Case 2. Trigonometric function: $H(\xi) = 0.5\sin(0.1\xi)$.

Letting $s = -50, m = -30, e = -5$, the LLEs correspond to the three parameters, as shown in Figures 15–17. Similar to the results of Gaussian perturbation analysis, it can be seen from the figures that the LLEs of the parameter $s$ corresponding to the cubic term are the largest, indicating that the parameter $s$ has a greater influence on the behavior of the perturbed system than $m$ and $e$. Meanwhile, the corresponding three-dimensional phase portrait is shown in Figure 18. We can find that the trajectories of the perturbed system (5.1) intersect and conclude that the solution is not unique.

To sum up, when given a suitable perturbation term, the perturbed system has chaotic behaviors, which are different for different perturbations. In the present study, we cannot determine which type of perturbation term will cause the system to exhibit chaotic behaviors. and we will continue to focus on this problem in the future.

6.   Conclusions

In this paper, the (3+1) dimensional NLSE is studied, which is important in describing the transmission of optical fibers under the complex interaction of self-phase modulation and spatial dispersion effects. The properties of dynamical systems are analyzed by using the trial equation method and the complete discrimination system for the polynomial method, for which the existence of soliton solutions and periodic solutions is proved. According to the complete discrimination system for the quartic polynomial method, thirteen exact solutions are obtained, including solitary wave solutions, triangular function solutions, rational solutions, and the elliptic function double periodic solutions. Especially, elliptic function double periodic solutions are new solutions. Then the exact solutions of the equation are obtained by the complete discrimination system for the quartic polynomial method. In particular, the elliptic function double periodic solutions are obtained, which are new solutions. The physical behaviors of the solutions are shown as 3D phase diagrams. In addition, for dynamical systems with suitable perturbation terms, we get the conclusion that chaotic behaviors exist in the system according to the relationship diagram of the LLEs and the corresponding three-dimensional phase diagrams, and the chaotic behaviors of the system vary with the types of perturbation terms. This is important for physics.

In summary, the methods presented in this paper are powerful tools for solving similar NDEs. We can determine the forms and types of the solutions based on the physical parameters, and the exact solutions obtained are more abundant. The discovery of chaotic behaviors in perturbation systems also provides us with new insight for the study of such equations.

Acknowledgments

This project is supported by the Special Programm for the Ability Promotion of the Basic and Scientific Research (No. 2023JCYJ-01).

Conflict of interest

The author declares no conflict of interest.