Explicit exact solutions and bifurcation analysis for the mZK equation with truncated M-fractional derivatives utilizing two reliable methods

Pim Malingam; Paiwan Wongsasinchai; Sekson Sirisubtawee

doi:10.1515/phys-2024-0117

Article Open Access

Explicit exact solutions and bifurcation analysis for the mZK equation with truncated M-fractional derivatives utilizing two reliable methods

Pim Malingam , Paiwan Wongsasinchai and Sekson Sirisubtawee

Published/Copyright: February 11, 2025

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 23 Issue 1

Abstract

The ( 2 + 1 ) -dimensional modified Zakharov–Kuznetsov (mZK) partial differential equation is of importance as a model for phenomena in various physical fields such as discrete electrical lattices, electrical waves in cold plasmas, nonlinear optical waves, deep ocean-waves, and the propagation of solitary gravity waves. In this study, the main objective is to give a detailed analysis of exact traveling wave solutions of the mZK equation with truncated M-fractional spatial–temporal partial derivatives. Using an appropriate traveling wave transformation and the homogeneous balance rule, the mZK equation is converted into a corresponding ordinary differential equation (ODE). The ( G ′ ∕ G , 1 ∕ G ) -expansion and Sardar subequation methods are then used to derive exact solutions of the ODE in the form of functions such as hyperbolic, trigonometric, and special generalized hyperbolic and trigonometric functions. The two methods give some novel solutions of the proposed model and are presented here for the first time. The fractional-order effects are studied through numerical simulations, including three-dimensional (3D), two-dimensional (2D), and contour plots. These numerical simulations clearly show physical interpretations of selected solutions. In particular, the generalized hyperbolic and trigonometric function solutions that have been derived by the Sardar subequation method are important as they provide examples of exact traveling wave solutions of various physical types. Furthermore, the results include examples of bifurcations and chaotic behaviors of the model through 2D and 3D plots when parameter values are varied. Finally, the methods of solution described in this study are reliable, powerful, and efficient and can be recommended to obtain traveling wave solutions of other nonlinear partial differential equations with truncated M-fractional derivatives.

Keywords: truncated M-fractional derivative; space–time fractional mZK equation; (G′/G, 1/G)-expansion method; Sardar subequation method; bifurcation analysis

1 Preliminary background

In recent years, nonlinear partial differential equations (NPDEs) have been employed as models for a wide variety of natural phenomena arising in the physical sciences, biological sciences, and engineering. Because of their importance, they have been investigated in many aspects, for instance, modeling NPDEs from natural rules and facts, constructing and developing efficient analytical and numerical methods for solving NPDEs, and interpreting solutions of NPDEs to explain dynamics of real-world phenomena. Some examples of their applications include the following: the quantum Hall effect [1], wave propagation on the surface of water [2], propagation of solitary waves in optical fibers [3], and nanoscale thermal transportation of ternary magnetized Carreau nanofluid using neural network architecture [4]. Exact traveling waves, including solitary waves and solitons, occur frequently in the real world and play an important role in describing many applications in nonlinear dynamics. A solitary wave is a traveling wave comprising a single peak or trough propagating in isolation whose size, shape, or speed does not change [5]. A solitary water wave is usually found in shallow water.

Over the last few decades, a wide range of NPDE models have been extensively studied in terms of their exact traveling wave solutions using many different approaches including the improved generalized Riccati equation method [6], the new extended direct algebraic method [7], the generalized Kudryashov scheme [8], the modified exponential function method [9], the ( G ′ ∕ G , 1 ∕ G ) -expansion method [10], and the Sardar subequation method [11].

The two-dimensional (2D) Zakharov–Kuznetsov (ZK) equation [12] was initially proposed by two Russian mathematicians, namely, Zakharov and Kuznetsov in 1974. It is the canonical 2D extension of the well-known Korteweg-de Vries (K-dV) equation [13]. The ZK equation is a model for weakly nonlinear ion-acoustic waves in intensive geophysical flows and magnetized undamaged plasma in two dimensions. The ZK equation is also used as a model for vulnerable nonlinear ion-acoustic waves arising in a mixture of hot isothermal electrons and cool ions in a uniform magnetic field [14,15]. In addition, this equation has broad applications in mathematical physics and engineering in areas, including ionic temperature, density gradient of ions, existence of dust, compound phenomena in 2D discrete electrical lattices, and inclined propagation. Further details for the ZK equation can be found in previous studies [16–18] and references therein. In recent years, fractional-order ZK models have also been studied with fractional derivatives, including the modified Riemann–Liouville fractional derivative and conformable fractional derivatives [16,19].

The modified Zakharov–Kuznetsov (mZK) equation was first proposed in 1999 [20,21]. The mZK equation is as follows:

(1) ∂ u ∂ t + u 2 ∂ u ∂ x + ∂ 3 u ∂ x 3 + ∂ 3 u ∂ y 2 ∂ x = 0 ,

where u = u ( x , y , t ) has distinct meanings depending upon the physical application. One motivation for its development was in the study of electron plasma. Since the plasma cannot satisfy the Boltzmann distribution, the ZK model was modified. The mZK equation has been of considerable interest for many years because it is an anisotropic 2D generalization of the KdV equation. One example of its application is the analysis of magnetized plasma in the minuscule magnitude Alfven wave in which a threshold angle is considered. In addition, the evolution of diversified solitary waves in magnetized plasma with isothermal multi-components can be effectively studied using the mZK equation.

Recently, there have been many applications of the truncated M-fractional partial derivatives in well-known equations arising in mathematical biology, physics, and engineering research fields. For example, Zhang et al. [22] investigated the dynamical behaviors of nonlinear traveling waves for the generalized reaction Duffing model and density-dependent diffusion–reaction equation in the sense of the M-fractional derivative with respect to temporal and spatial variables. The authors employed the newly extended direct algebraic technique to obtain new exact solutions expressed in terms of rational, hyperbolic, and trigonometric functions. In the study of Siddique et al. [23], the ( 1 ∕ G ′ ) , modified ( G ′ ∕ G 2 ) -expansion methods and new extended direct algebraic methods were used to explore novel exact traveling wave solutions of the time M-truncated fractional modified equal-width (MEW) equation. Moreover, the bifurcation behavior of the proposed equation was analyzed by considering the influence of the system’s parameters on the solution. Further recent literature related to the use of fractional derivatives including the M-truncated fractional derivative and the beta-derivative can be found in previous studies [24–26] and references therein.

In this study, the mZK Eq. (1) is modified by changing its partial derivatives to the truncated M-fractional partial derivatives for both time and space dimensions. After applying the truncated M-fractional partial derivatives to (1), we obtain the space–time fractional mZK equation given by

(2) ∂ M , t α , γ m u + u 2 ∂ M , x β , γ m u + ∂ M , x β , γ m ( ∂ M , x β , γ m ( ∂ M , x β , γ m u ) ) + ∂ M , y β , γ m ( ∂ M , y β , γ m ( ∂ M , x β , γ m u ) ) = 0 ,

where ∂ M , t α , γ m ( ⋅ ) , ∂ M , x β , γ m ( ⋅ ) , and ∂ M , y β , γ m ( ⋅ ) are the truncated M-fractional partial derivative operators of orders 0 < α , β ≤ 1 , respectively, with regard to t > 0 , x > 0 , and y > 0 . A literature review of recent studies on the fractional mZK equation in space–time dimensions is as follows. The fractional mZK equation with a modified derivatives of Jumarie for both space and time was solved to obtain exact analytical solutions employing the fractional sub-equation approach [16]. The solutions of the equation obtained by Ray and Sahoo [16] consisted of generalized hyperbolic and trigonometric functions. In the study of Guo et al. [27], a new time-fractional mZK equation in the sense of the Riesz fractional derivative obtained by a multiscale perturbation technique was introduced. Its solutions were obtained by employing the sech–tanh technique. In the study of Bibi et al. [28], the ( m + ( G ′ ∕ G ) ) -expansion approach was applied to the beta-fractional mZK equation. The results were found in many forms such as kink, dark, W-type soliton, and singularly periodic solutions. The influence of distinct fractional orders on the wave solutions was also examined. Further reviews of the fractional mZK equation can be studied in previous studies [29,30].

In this article, our main aim is to use two reliable methods, namely, the ( G ′ ∕ G , 1 ∕ G ) -expansion and Sardar subequation methods, to derive exact analytical traveling wave solutions of (2). The organization of this article is as follows. In Section 2, a definition of the truncated M-fractional derivative and its interesting properties are given. Also, the important procedures of the ( G ′ ∕ G , 1 ∕ G ) -expansion and Sardar subequation methods are described in this section. In Section 3, the derivation of the exact solutions of problem (2) by the two methods is given. Section 4 shows the bifurcation analysis and chaotic behaviors of the model. Section 5 provides the graphical forms of the chosen solutions. Section 6 provides some vital discussion and conclusions.

2 Methodology

In this section, the truncated one-parameter M-fractional derivative is defined and its main characteristics are described. The details of the ( G ′ ∕ G , 1 ∕ G ) -expansion and Sardar subequation techniques are also given. The Mittag–Leffler function (MLF) for one-parameter γ > 0 [31,32] is defined as E γ ( z ) = ∑ n = 0 ∞ z n Γ ( γ n + 1 ) , where z ∈ C and Γ ( ⋅ ) is the gamma function.

2.1 Truncated M-fractional derivative and its characteristics

Here, definitions of the truncated one-parameter MLF and its derivative are provided. Also, the important properties of the derivative are given.

Definition 2.1

The truncated one-parameter MLF is defined as follows [33–36]:

(3) E γ m ( z ) = ∑ n = 0 m z n Γ ( γ n + 1 ) ,

where z is a complex number and the parameter γ > 0 .

Definition 2.2

Let p be a function from [ 0 , ∞ ) into R . Then, the truncated M-fractional derivative of p with order 0 < β ≤ 1 is defined by [33–36]

(4) D M , t β , γ m p ( t ) = lim τ → 0 p ( t E γ m ( τ t − β ) ) − p ( t ) τ ,

where E γ m is the truncated one-parameter MLF defined in (3).

Some useful characteristics of the derivative (4) are as follows [33,36–42]. Let p ( t ) , q ( t ) be β -truncated M-fractional differentiable functions for all t > 0 , β ∈ ( 0 , 1 ] , and γ > 0 . Then, we have

D M , t β , γ m ( λ ) = 0 , ∀ λ ∈ R .
D M , t β , γ m ( a p ( t ) + b q ( t ) ) = a D M , t β , γ m p ( t ) + b D M , t β , γ m q ( t ) , ∀ a , b ∈ R .
D M , t β , γ m ( p ( t ) q ( t ) ) = p ( t ) D M , t β , γ m q ( t ) + q ( t ) D M , t β , γ m p ( t ) .
D M , t β , γ m p ( t ) q ( t ) = q ( t ) D M , t β , γ m p ( t ) − p ( t ) D M , t β , γ m q ( t ) ( q ( t ) ) 2 , where q ( t ) ≠ 0 .
D M , t β , γ m ( p o q ) ( t ) = p ′ ( q ( t ) ) m D M , t β , γ q ( t ) when p can be differentiated at q ( t ) .
D M , t β , γ m ( p ( t ) ) = t 1 − β Γ ( γ + 1 ) d p ( t ) d t when p is differentiable.

As discussed in the study of Sousa and de Oliveira [33], many important theorems of the classical calculus also apply to the truncated M-fractional derivative (4) including Rolle’s theorem, the mean value theorem, and the inverse property of the corresponding M-fractional integral. The advantage of the operator (4) is that it can be transformed to the conformable fractional derivative [43], the alternative fractional derivative [44], and the generalized fractional derivative [44] by selecting appropriate parameters [33].

2.2 Methods

In this subsection, the algorithms of the ( G ′ ∕ G , 1 ∕ G ) -expansion and Sardar subequation approaches are summarized. First, the common procedures between these two methods are as follows. We consider the following truncated M-fractional NPDE:

(5) F ( v , ∂ M , t α , γ m v , ∂ M , x β , γ m v , ∂ M , y δ , γ m v , ∂ M , t α , γ m ( ∂ M , x β , γ m v ) , ∂ M , t α , γ m ( ∂ M , y δ , γ m v ) , ∂ M , x β , γ m ( ∂ M , y δ , γ m v ) , … ) = 0 ,

where 0 < α , β , δ ≤ 1 , γ > 0 , and ∂ M , t α , γ m v , ∂ M , x β , γ m v , and ∂ M , y δ , γ m v are the truncated M-fractional derivatives of the response variable v = v ( x , y , t ) in the explanatory variables t , x , and y . The arguments of the polynomial function F are v and its partial derivatives. The next step is to apply the following traveling wave transformation to (5):

(6) V ( ϱ ) = v ( x , y , t ) , ϱ = Γ ( γ + 1 ) k x β β + l y δ δ + c t α α + d ,

where the constants k , l , c , and d are to be discovered later. Eq. (5) turns out to be an ODE in V = V ( ϱ ) as follows:

(7) H ( V , V ′ , V ″ , … ) = 0 ,

where the polynomial H has the arguments V = V ( ϱ ) and its derivatives V ′ = d V d ϱ , V ″ = d 2 V d ϱ 2 , … . Next, we introduce a brief description of the two methods for solving the ODE (7).

2.2.1 ( G ′ ∕ G , 1 ∕ G ) -expansion technique

Before giving an algorithm of the ( G ′ ∕ G , 1 ∕ G ) -expansion approach, we need to use the following auxiliary equation and its solutions [10,42,45–48]. Consider the following auxiliary equation defined as the linear ODE:

(8) G ″ ( ϱ ) + λ G ( ϱ ) = μ ,

where λ and μ are constants. We then define Φ and Ψ as

(9) Φ ( ϱ ) = G ′ ( ϱ ) G ( ϱ ) and Ψ ( ϱ ) = 1 G ( ϱ ) .

Then, using Eq. (9), we can convert Eq. (8) into the system of two nonlinear first-order ODEs

(10) Φ ′ = − Φ 2 + μ Ψ − λ , Ψ ′ = − Φ Ψ .

Explicit solutions of Eq. (8) depend on values of λ and can be classified as the following three cases.

Case 1: When λ < 0 , the solution of (8) is expressed as

(11) G ( ϱ ) = A 1 sinh ( ϱ − λ ) + A 2 cosh ( ϱ − λ ) + μ λ ,

with the following relationship:

(12) Ψ 2 = − λ λ 2 σ 1 + μ 2 ( Φ 2 − 2 μ Ψ + λ ) ,

where σ 1 = A 1 2 − A 2 2 when A 1 and A 2 are arbitrary real constants.

Case 2: When λ > 0 , the general solution of (8) is obtained as

(13) G ( ϱ ) = A 1 sin ( ϱ λ ) + A 2 cos ( ϱ λ ) + μ λ ,

and the associated relation for this case is

(14) Ψ 2 = λ λ 2 σ 2 − μ 2 ( Φ 2 − 2 μ Ψ + λ ) ,

where σ 2 = A 1 2 + A 2 2 when A 1 and A 2 are arbitrary real constants.

Case 3: When λ = 0 , we have the solution of (8) as follows:

(15) G ( ϱ ) = μ 2 ϱ 2 + A 1 ϱ + A 2 ,

with the following associated relation:

(16) Ψ 2 = 1 A 1 2 − 2 μ A 2 ( Φ 2 − 2 μ Ψ ) ,

where A 1 and A 2 are arbitrary real constants.

The main steps of the ( G ′ ∕ G , 1 ∕ G ) -expansion technique [10,42,45–48] are as follows.

Step 1: A solution of Eq. (7) is assumed to be a polynomial of Φ and Ψ as

(17) V ( ϱ ) = a 0 + ∑ m = 1 N a m Φ m + ∑ m = 1 N b m Φ m − 1 Ψ ,

where a 0 , a m , and b m for m = 1 , 2 , … , N are the coefficients to be found with a N 2 + b N 2 ≠ 0 and where Φ = Φ ( ϱ ) and Ψ = Ψ ( ϱ ) are implicitly involved with Eq. (8) through the relationship (9).

The description of Steps 2–5 of the method is the same as Steps 2–5 of our previous publication [42]. We will omit these steps here to avoid repetition.

2.2.2 Sardar subequation approach

In this subsection, the main steps of the Sardar subequation method [11,49–52] are briefly described as follows.

Step 1: After transforming Eq. (5) into Eq. (7) via the transformation (6), we set up, by the Sardar method, a solution of Eq. (7) as

(18) V ( ϱ ) = ∑ m = 0 N ϖ m Q m ( ϱ ) ,

where ϖ m , ( m = 0 , 1 , 2 , … N ) are the coefficients of Q m with ϖ N ≠ 0 . In addition, Q ( ϱ ) is assumed to satisfy the following auxiliary equation:

(19) ( Q ′ ( ϱ ) ) 2 = μ + a Q 2 ( ϱ ) + Q 4 ( ϱ ) ,

where constants μ and a are real. The general solutions of Eq. (19) can be written as

Case 1: If a > 0 and μ = 0 , then

(20) Q 1 ± ( ϱ ) = ± − p q a sech p q ( a ϱ ) , Q 2 ± ( ϱ ) = ± p q a csch p q ( a ϱ ) ,

where sech p q ( ϱ ) = 2 p e ϱ + q e − ϱ and csch p q ( ϱ ) = 2 p e ϱ − q e − ϱ .

Case 2: If a < 0 and μ = 0 , then

(21) Q 3 ± ( ϱ ) = ± − p q a sec p q ( − a ϱ ) , Q 4 ± ( ϱ ) = ± − p q a csc p q ( − a ϱ ) ,

where sec p q ( ϱ ) = 2 p e i ϱ + q e − i ϱ and csc p q ( ϱ ) = 2 i p e i ϱ − q e − i ϱ .

Case 3: If a < 0 and μ = a 2 4 , then

(22) Q 5 ± ( ϱ ) = ± − a 2 tanh p q − a 2 ϱ , Q 6 ± ( ϱ ) = ± − a 2 coth p q − a 2 ϱ , Q 7 ± ( ϱ ) = ± − a 2 ( tanh p q ( − 2 a ϱ ) ± i p q sech p q ( − 2 a ϱ ) ) , Q 8 ± ( ϱ ) = ± − a 2 ( coth p q ( − 2 a ϱ ) ± p q csch p q ( − 2 a ϱ ) ) , Q 9 ± ( ϱ ) = ± − a 8 tanh p q − a 8 ϱ + coth p q − a 8 ϱ ,

where tanh p q ( ϱ ) = p e ϱ − q e − ϱ p e ϱ + q e − ϱ and coth p q ( ϱ ) = p e ϱ + q e − ϱ p e ϱ − q e − ϱ .

Case 4: If a > 0 and μ = a 2 4 , then

(23) Q 10 ± ( ϱ ) = ± a 2 tan p q a 2 ϱ , Q 11 ± ( ϱ ) = ± a 2 cot p q a 2 ϱ , Q 12 ± ( ϱ ) = ± a 2 ( tan p q ( 2 a ϱ ) ± p q sec p q ( 2 a ϱ ) ) , Q 13 ± ( ϱ ) = ± a 2 ( cot p q ( 2 a ϱ ) ± p q csc p q ( 2 a ϱ ) ) , Q 14 ± ( ϱ ) = ± a 8 tan p q a 8 ϱ − cot p q a 8 ϱ ,

where tan p q ( ϱ ) = − i p e i ϱ − q e − i ϱ p e i ϱ + q e − i ϱ and cot p q ( ϱ ) = i p e i ϱ + q e − i ϱ p e i ϱ − q e − i ϱ .

Step 2: Compute the value of N in the solution form (18) using the homogeneous balance principle as stated in Step 2 of the study of Malingam et al. [42].

Step 3: Replacing the value of N into solution (18) and substituting the resulting solution and its necessary derivatives in Eq. (7) with a use of Eq. (19), a polynomial equation in Q m ( ϱ ) is obtained. Equating the coefficients of Q m ( ϱ ) of the obtained polynomial to zero, a system of nonlinear equations in ϖ m (for m = 0 , 1 , 2 , … N ), k , l , c and d is derived. This system of nonlinear equations can then be solved for their symbolic values using Maple 17.

Step 4: The exact solutions of Eq. (5) can then be derived by replacing the wave transformation (6), the solutions Q ( ϱ ) of (19), and the results obtained by Step 3 into solution form (18). These explicit solutions can be obtained using Maple 17.

These two reliable methods are simple to implement and can provide abundant solutions to the NPDEs, including hyperbolic, trigonometric, and exponential rational functions.

3 Implementation of the schemes

This section investigates the use of the two analytical approaches, namely, the ( G ′ ∕ G , 1 ∕ G ) -expansion and Sardar subequation methods, as described earlier, to derive explicit exact solutions of Eq. (2). Before applying the methods to Eq. (2), its exact solution is assumed to be in the form

(24) u ( x , y , t ) = U ( ϱ ) ,

where the traveling wave transformation ϱ is

(25) ϱ = Γ ( γ + 1 ) x β β + y β β − ρ t α α ,

where 0 < α , β ≤ 1 , the parameter γ > 0 , and ρ ≠ 0 is a constant, which will be found later.

Substituting solutions (24) and (25) into the proposed problem (2) and using the properties of the derivative of U as mentioned earlier, we have the ODE

(26) − ρ U ′ ( ϱ ) + U 2 ( ϱ ) U ′ ( ϱ ) + 2 U ‴ ( ϱ ) = 0 .

Integrating Eq. (26) with regard to ϱ once, we obtain

(27) − ρ U ( ϱ ) + 1 3 U 3 ( ϱ ) + 2 U ″ ( ϱ ) = 0 .

Using the balance principle and the degree formulas expressed in the study of Malingam et al. [42], we obtain the balance number N = 1 .

3.1 Application of the ( G ′ ∕ G , 1 ∕ G ) -expansion scheme

Following the ( G ′ ∕ G , 1 ∕ G ) -expansion scheme and using N = 1 , we set up the solution form of Eq. (27) as

(28) U ( ϱ ) = a 0 + a 1 Φ ( ϱ ) + b 1 Ψ ( ϱ ) ,

where a 0 , a 1 , and b 1 are the constant coefficients such that a 1 2 + b 1 2 ≠ 0 , which will be determined later. As described in Section 2.2.1, there are three cases for the functions Φ ( ϱ ) and Ψ ( ϱ ) in (28) that are associated with equations (11), (13), and (15) according to the sign of λ .

Case 1 ( λ < 0 ): For this case, we have hyperbolic function solutions, which we can obtain as follows. We insert Eq. (28) along with Eqs (10) and (12) into Eq. (27). Then, the resulting equation of Eq. (27) has the left-hand expression in Φ ( ϱ ) and Ψ ( ϱ ) . Equating all coefficients of the obtained polynomial to zero, a system of nonlinear equations is obtained as follows:

Φ 3 : λ 4 A 1 4 a 1 3 − 2 λ 4 A 1 2 A 2 2 a 1 3 + λ 4 A 2 4 a 1 3 + 12 λ 4 A 1 4 a 1 − 24 λ 4 A 1 2 A 2 2 a 1 + 12 λ 4 A 2 4 a 1 + 2 λ 2 μ 2 A 1 2 a 1 3 − 2 λ 2 μ 2 A 2 2 a 1 3 − 3 λ 3 A 1 2 a 1 b 1 2 + 3 λ 3 A 2 2 a 1 b 1 2 + 24 λ 2 μ 2 A 1 2 a 1 − 24 λ 2 μ 2 A 2 2 a 1 + μ 4 a 1 3 − 3 λ μ 2 a 1 b 1 2 + 12 μ 4 a 1 = 0 , Φ 2 : 3 λ 4 A 1 4 a 0 a 1 2 − 6 λ 4 A 1 2 A 2 2 a 0 a 1 2 + 3 λ 4 A 2 4 a 0 a 1 2 + 6 λ 2 μ 2 A 1 2 a 0 a 1 2 − 6 λ 2 μ 2 A 2 2 a 0 a 1 2 − 3 λ 3 A 1 2 a 0 b 1 2 + 3 λ 3 A 2 2 a 0 b 1 2 + 6 λ 3 μ A 1 2 b 1 − 6 λ 3 μ A 2 2 b 1 + 3 μ 4 a 0 a 1 2 − 2 λ 2 μ b 1 3 − 3 λ μ 2 a 0 b 1 2 + 6 λ μ 3 b 1 = 0 , Φ 2 Ψ : 3 λ 4 A 1 4 a 1 2 b 1 − 6 λ 4 A 1 2 A 2 2 a 1 2 b 1 + 3 λ 4 A 2 4 a 1 2 b 1 + 12 λ 4 A 1 4 b 1 − 24 λ 4 A 1 2 A 2 2 b 1 + 12 λ 4 A 2 4 b 1 + 6 λ 2 μ 2 A 1 2 a 1 2 b 1 − 6 λ 2 μ 2 A 2 2 a 1 2 b 1 − λ 3 A 1 2 b 1 3 + λ 3 A 2 2 b 1 3 + 24 λ 2 μ 2 A 1 2 b 1 − 24 λ 2 μ 2 A 2 2 b 1 + 3 μ 4 a 1 2 b 1 − λ μ 2 b 1 3 + 12 μ 4 b 1 = 0 ,

Φ : 3 λ 4 A 1 4 a 0 2 a 1 − 6 λ 4 A 1 2 A 2 2 a 0 2 a 1 + 3 λ 4 A 2 4 a 0 2 a 1 + 12 λ 5 A 1 4 a 1 − 24 λ 5 A 1 2 A 2 2 a 1 + 12 λ 5 A 2 4 a 1 − 3 λ 4 ρ A 1 4 a 1 + 6 λ 4 ρ A 1 2 A 2 2 a 1 − 3 λ 4 ρ A 2 4 a 1 − 3 λ 4 A 1 2 a 1 b 1 2 + 3 λ 4 A 2 2 a 1 b 1 2 + 6 λ 2 μ 2 A 1 2 a 0 2 a 1 − 6 λ 2 μ 2 A 2 2 a 0 2 a 1 + 24 λ 3 μ 2 A 1 2 a 1 − 24 λ 3 μ 2 A 2 2 a 1 − 6 λ 2 μ 2 ρ A 1 2 a 1 + 6 λ 2 μ 2 ρ A 2 2 a 1 − 3 λ 2 μ 2 a 1 b 1 2 + 3 μ 4 a 0 2 a 1 + 12 λ μ 4 a 1 − 3 μ 4 ρ a 1 = 0 ,

Φ Ψ : 6 λ 4 A 1 4 a 0 a 1 b 1 − 12 λ 4 A 1 2 A 2 2 a 0 a 1 b 1 + 6 λ 4 A 2 4 a 0 a 1 b 1 − 18 λ 4 μ A 1 4 a 1 + 36 λ 4 μ A 1 2 A 2 2 a 1 − 18 λ 4 μ A 2 4 a 1 + 6 λ 3 μ A 1 2 a 1 b 1 2 − 6 λ 3 μ A 2 2 a 1 b 1 2 + 12 λ 2 μ 2 A 1 2 a 0 a 1 b 1 − 12 λ 2 μ 2 A 2 2 a 0 a 1 b 1 − 36 λ 2 μ 3 A 1 2 a 1 + 36 λ 2 μ 3 A 2 2 a 1 + 6 λ μ 3 a 1 b 1 2 + 6 μ 4 a 0 a 1 b 1 − 18 μ 5 a 1 = 0 , Ψ : 3 λ 4 A 1 4 a 0 2 b 1 − 6 λ 4 A 1 2 A 2 2 a 0 2 b 1 + 3 λ 4 A 2 4 a 0 2 b 1 + 6 λ 5 A 1 4 b 1 − 12 λ 5 A 1 2 A 2 2 b 1 + 6 λ 5 A 2 4 b 1 − 3 λ 4 ρ A 1 4 b 1 + 6 λ 4 ρ A 1 2 A 2 2 b 1 − 3 λ 4 ρ A 2 4 b 1 − λ 4 A 1 2 b 1 3 + λ 4 A 2 2 b 1 3 + 6 λ 3 μ A 1 2 a 0 b 1 2 − 6 λ 3 μ A 2 2 a 0 b 1 2 + 6 λ 2 μ 2 A 1 2 a 0 2 b 1 − 6 λ 2 μ 2 A 2 2 a 0 2 b 1 − 6 λ 2 μ 2 ρ A 1 2 b 1 + 6 λ 2 μ 2 ρ A 2 2 b 1 + 3 λ 2 μ 2 b 1 3 + 6 λ μ 3 a 0 b 1 2 + 3 μ 4 a 0 2 b 1 − 6 λ μ 4 b 1 − 3 μ 4 ρ b 1 = 0 ,

Φ 0 : λ 4 A 1 4 a 0 3 − 2 λ 4 A 1 2 A 2 2 a 0 3 + λ 4 A 2 4 a 0 3 − 3 λ 4 ρ A 1 4 a 0 + 6 λ 4 ρ A 1 2 A 2 2 a 0 − 3 λ 4 ρ A 2 4 a 0 − 3 λ 4 A 1 2 a 0 b 1 2 + 3 λ 4 A 2 2 a 0 b 1 2 + 2 λ 2 μ 2 A 1 2 a 0 3 − 2 λ 2 μ 2 A 2 2 a 0 3 + 6 λ 4 μ A 1 2 b 1 − 6 λ 4 μ A 2 2 b 1 − 6 λ 2 μ 2 ρ A 1 2 a 0 + 6 λ 2 μ 2 ρ A 2 2 a 0 − 2 λ 3 μ b 1 3 − 3 λ 2 μ 2 a 0 b 1 2 + μ 4 a 0 3 + 6 λ 2 μ 3 b 1 − 3 μ 4 ρ a 0 = 0 .

Using Maple 17 to solve the aforementioned system, we then obtain the outputs as follows:

Result 1

(29) λ = λ , μ = μ , ρ = λ , a 0 = 0 , a 1 = ± 3 i , b 1 = ± 3 λ σ 1 + 3 μ 2 λ ,

where σ 1 = A 1 2 − A 2 2 and λ ( < 0 ) , μ , A 1 , and A 2 are arbitrary constants. Utilizing Eqs (9), (11), (24), (28), and (29), the exact solution of Eq. (2) is

(30) u ± ( x , y , t ) = ± 3 λ ( A 1 cosh ( ϱ − λ ) + A 2 sinh ( ϱ − λ ) ) A 1 sinh ( ϱ − λ ) + A 2 cosh ( ϱ − λ ) + μ λ ± 3 λ σ 1 + 3 μ 2 λ A 1 sinh ( ϱ − λ ) + A 2 cosh ( ϱ − λ ) + μ λ ,

where ϱ = Γ ( γ + 1 ) x β β + y β β − ρ t α α .

Result 2

(31) λ = λ , μ = 0 , ρ = − 2 λ , a 0 = 0 , a 1 = 0 , b 1 = ± 2 3 λ σ 1 ,

where σ 1 = A 1 2 − A 2 2 and λ ( < 0 ) , A 1 , and A 2 are arbitrary constants. Using Eqs (9), (11), (24), (28), and (31), the exact solution of Eq. (2) can be expressed as

(32) u ± ( x , y , t ) = ± 2 3 λ σ 1 A 1 sinh ( ϱ − λ ) + A 2 cosh ( ϱ − λ ) ,

where ϱ = Γ ( γ + 1 ) x β β + y β β − ρ t α α .

Result 3

(33) λ = λ , μ = 0 , ρ = 4 λ , a 0 = 0 , a 1 = ± 2 3 i , b 1 = 0 ,

where λ ( < 0 ) is an arbitrary constant. Using Eqs (9), (11), (24), (28), and (33), the explicit exact solution for Eq. (2) is

(34) u ± ( x , y , t ) = ± 2 3 λ ( A 1 cosh ( ϱ − λ ) + A 2 sinh ( ϱ − λ ) ) A 1 sinh ( ϱ − λ ) + A 2 cosh ( ϱ − λ ) ,

where A 1 and A 2 are arbitrary constants and

ϱ = Γ ( γ + 1 ) x β β + y β β − ρ t α α .

Case 2 ( λ > 0 ): For this case, we obtain trigonometric function solutions by doing some steps as follows. Inserting Eq. (28) along with Eqs (10) and (14) into Eq. (27), the left-hand expression of Eq. (27) turns out to be a polynomial of Φ ( ϱ ) and Ψ ( ϱ ) . Setting all coefficients of the obtained result to zero, the nonlinear equations in λ , μ , ρ , a 0 , a 1 , and b 1 are obtained as

Φ 3 : λ 4 A 1 4 a 1 3 + 2 λ 4 A 1 2 A 2 2 a 1 3 + λ 4 A 2 4 a 1 3 + 12 λ 4 A 1 4 a 1 + 24 λ 4 A 1 2 A 2 2 a 1 + 12 λ 4 A 2 4 a 1 − 2 λ 2 μ 2 A 1 2 a 1 3 − 2 λ 2 μ 2 A 2 2 a 1 3 + 3 λ 3 A 1 2 a 1 b 1 2 + 3 λ 3 A 2 2 a 1 b 1 2 − 24 λ 2 μ 2 A 1 2 a 1 − 24 λ 2 μ 2 A 2 2 a 1 + μ 4 a 1 3 − 3 λ μ 2 a 1 b 1 2 + 12 μ 4 a 1 = 0 , Φ 2 : 3 λ 4 A 1 4 a 0 a 1 2 + 6 λ 4 A 1 2 A 2 2 a 0 a 1 2 + 3 λ 4 A 2 4 a 0 a 1 2 − 6 λ 2 μ 2 A 1 2 a 0 a 1 2 − 6 λ 2 μ 2 A 2 2 a 0 a 1 2 + 3 λ 3 A 1 2 a 0 b 1 2 + 3 λ 3 A 2 2 a 0 b 1 2 − 6 λ 3 μ A 1 2 b 1 − 6 λ 3 μ A 2 2 b 1 + 3 μ 4 a 0 a 1 2 − 2 λ 2 μ b 1 3 − 3 λ μ 2 a 0 b 1 2 + 6 λ μ 3 b 1 = 0 , Φ 2 Ψ : 3 λ 4 A 1 4 a 1 2 b 1 + 6 λ 4 A 1 2 A 2 2 a 1 2 b 1 + 3 λ 4 A 2 4 a 1 2 b 1 + 12 λ 4 A 1 4 b 1 + 24 λ 4 A 1 2 A 2 2 b 1 + 12 λ 4 A 2 4 b 1 − 6 λ 2 μ 2 A 1 2 a 1 2 b 1 − 6 λ 2 μ 2 A 2 2 a 1 2 b 1 + λ 3 A 1 2 b 1 3 + λ 3 A 2 2 b 1 3 − 24 λ 2 μ 2 A 1 2 b 1 − 24 λ 2 μ 2 A 2 2 b 1 + 3 μ 4 a 1 2 b 1 − λ μ 2 b 1 3 + 12 μ 4 b 1 = 0 ,

Φ : 3 λ 4 A 1 4 a 0 2 a 1 + 6 λ 4 A 1 2 A 2 2 a 0 2 a 1 + 3 λ 4 A 2 4 a 0 2 a 1 + 12 λ 5 A 1 4 a 1 + 24 λ 5 A 1 2 A 2 2 a 1 + 12 λ 5 A 2 4 a 1 − 3 λ 4 ρ A 1 4 a 1 − 6 λ 4 ρ A 1 2 A 2 2 a 1 − 3 λ 4 ρ A 2 4 a 1 + 3 λ 4 A 1 2 a 1 b 1 2 + 3 λ 4 A 2 2 a 1 b 1 2 − 6 λ 2 μ 2 A 1 2 a 0 2 a 1 − 6 λ 2 μ 2 A 2 2 a 0 2 a 1 − 24 λ 3 μ 2 A 1 2 a 1 − 24 λ 3 μ 2 A 2 2 a 1 + 6 λ 2 μ 2 ρ A 1 2 a 1 + 6 λ 2 μ 2 ρ A 2 2 a 1 − 3 λ 2 μ 2 a 1 b 1 2 + 3 μ 4 a 0 2 a 1 + 12 λ μ 4 a 1 − 3 μ 4 ρ a 1 = 0 ,

Φ Ψ : 6 λ 4 A 1 4 a 0 a 1 b 1 + 12 λ 4 A 1 2 A 2 2 a 0 a 1 b 1 + 6 λ 4 A 2 4 a 0 a 1 b 1 − 18 λ 4 μ A 1 4 a 1 − 36 λ 4 μ A 1 2 A 2 2 a 1 − 18 λ 4 μ A 2 4 a 1 − 6 λ 3 μ A 1 2 a 1 b 1 2 − 6 λ 3 μ A 2 2 a 1 b 1 2 − 12 λ 2 μ 2 A 1 2 a 0 a 1 b 1 − 12 λ 2 μ 2 A 2 2 a 0 a 1 b 1 + 36 λ 2 μ 3 A 1 2 a 1 + 36 λ 2 μ 3 A 2 2 a 1 + 6 λ μ 3 a 1 b 1 2 + 6 μ 4 a 0 a 1 b 1 − 18 μ 5 a 1 = 0 ,

Ψ : 3 λ 4 A 1 4 a 0 2 b 1 + 6 λ 4 A 1 2 A 2 2 a 0 2 b 1 + 3 λ 4 A 2 4 a 0 2 b 1 + 6 λ 5 A 1 4 b 1 + 12 λ 5 A 1 2 A 2 2 b 1 + 6 λ 5 A 2 4 b 1 − 3 λ 4 ρ A 1 4 b 1 − 6 λ 4 ρ A 1 2 A 2 2 b 1 − 3 λ 4 ρ A 2 4 b 1 + λ 4 A 1 2 b 1 3 + λ 4 A 2 2 b 1 3 − 6 λ 3 μ A 1 2 a 0 b 1 2 − 6 λ 3 μ A 2 2 a 0 b 1 2 − 6 λ 2 μ 2 A 1 2 a 0 2 b 1 − 6 λ 2 μ 2 A 2 2 a 0 2 b 1 + 6 λ 2 μ 2 ρ A 1 2 b 1 + 6 λ 2 μ 2 ρ A 2 2 b 1 + 3 λ 2 μ 2 b 1 3 + 6 λ μ 3 a 0 b 1 2 + 3 μ 4 a 0 2 b 1 − 6 λ μ 4 b 1 − 3 μ 4 ρ b 1 = 0 ,

Φ 0 : λ 4 A 1 4 a 0 3 + 2 λ 4 A 1 2 A 2 2 a 0 3 + λ 4 A 2 4 a 0 3 − 3 λ 4 ρ A 1 4 a 0 − 6 λ 4 ρ A 1 2 A 2 2 a 0 − 3 λ 4 ρ A 2 4 a 0 + 3 λ 4 A 1 2 a 0 b 1 2 + 3 λ 4 A 2 2 a 0 b 1 2 − 2 λ 2 μ 2 A 1 2 a 0 3 − 2 λ 2 μ 2 A 2 2 a 0 3 − 6 λ 4 μ A 1 2 b 1 − 6 λ 4 μ A 2 2 b 1 + 6 λ 2 μ 2 ρ A 1 2 a 0 + 6 λ 2 μ 2 ρ A 2 2 a 0 − 2 λ 3 μ b 1 3 − 3 λ 2 μ 2 a 0 b 1 2 + μ 4 a 0 3 + 6 λ 2 μ 3 b 1 − 3 μ 4 ρ a 0 = 0 .

Utilizing Maple 17 to solve the aforementioned system, we have

Result 1

(35) λ = λ , μ = μ , ρ = λ , a 0 = 0 , a 1 = ± 3 i , b 1 = ± − 3 λ σ 2 + 3 μ 2 λ ,

where σ 2 = A 1 2 + A 2 2 and λ ( > 0 ) , μ , A 1 , and A 2 are arbitrary constants. Using Eqs (9), (13), (24), (28), and (35), the explicit exact solution of Eq. (2) is

(36) u ± ( x , y , t ) = ± 3 λ ( A 1 cos ( ϱ λ ) − A 2 sin ( ϱ λ ) ) A 1 sin ( ϱ λ ) + A 2 cos ( ϱ λ ) + μ λ i ± − 3 λ σ 2 + 3 μ 2 λ A 1 sin ( ϱ λ ) + A 2 cos ( ϱ λ ) + μ λ ,

where ϱ = Γ ( γ + 1 ) x β β + y β β − ρ t α α .

Result 2

(37) λ = λ , μ = 0 , ρ = 4 λ , a 0 = 0 , a 1 = ± 2 3 i , b 1 = 0 ,

where λ ( > 0 ) is an arbitrary constant. Employing Eqs (9), (13), (24), (28), and (37), we obtain the exact solution of Eq. (2) for this case as

(38) u ± ( x , y , t ) = ± 2 3 λ ( A 1 cos ( ϱ λ ) − A 2 sin ( ϱ λ ) ) A 1 sin ( ϱ λ ) + A 2 cos ( ϱ λ ) i ,

where A 1 and A 2 are arbitrary constants and

ϱ = Γ ( γ + 1 ) x β β + y β β − ρ t α α .

Case 3 ( λ = 0 ): For this case, we obtain rational function solutions by proceeding as follows. Substituting Eq. (28) along with Eqs (10) and (16) into Eq. (27), the left-hand expression of Eq. (27) is converted into a polynomial function of Φ ( ϱ ) and Ψ ( ϱ ) . Letting all coefficients of this polynomial be zero, a set of nonlinear algebraic equations for λ , μ , ρ , a 0 , a 1 , and b 1 is obtained as

Φ 3 : 4 μ 2 A 2 2 a 1 3 − 4 μ A 1 2 A 2 a 1 3 + A 1 4 a 1 3 + 48 μ 2 A 2 2 a 1 − 48 μ A 1 2 A 2 a 1 − 6 μ A 2 a 1 b 1 2 + 12 A 1 4 a 1 + 3 A 1 2 a 1 b 1 2 = 0 , Φ 2 : 12 μ 2 A 2 2 a 0 a 1 2 − 12 μ A 1 2 A 2 a 0 a 1 2 + 3 A 1 4 a 0 a 1 2 − 6 μ A 2 a 0 b 1 2 + 3 A 1 2 a 0 b 1 2 + 12 μ 2 A 2 b 1 − 6 μ A 1 2 b 1 − 2 μ b 1 3 = 0 , Φ 2 Ψ : 12 μ 2 A 2 2 a 1 2 b 1 − 12 μ A 1 2 A 2 a 1 2 b 1 + 3 A 1 4 a 1 2 b 1 + 48 μ 2 A 2 2 b 1 − 48 μ A 1 2 A 2 b 1 − 2 μ A 2 b 1 3 + 12 A 1 4 b 1 + A 1 2 b 1 3 = 0 , Φ : 12 μ 2 A 2 2 a 0 2 a 1 − 12 μ A 1 2 A 2 a 0 2 a 1 + 3 A 1 4 a 0 2 a 1 − 12 μ 2 ρ A 2 2 a 1 + 12 μ ρ A 1 2 A 2 a 1 − 3 ρ A 1 4 a 1 = 0 , Φ Ψ : 24 μ 2 A 2 2 a 0 a 1 b 1 − 24 μ A 1 2 A 2 a 0 a 1 b 1 + 6 A 1 4 a 0 a 1 b 1 − 72 μ 3 A 2 2 a 1 + 72 μ 2 A 1 2 A 2 a 1 + 12 μ 2 A 2 a 1 b 1 2 − 18 μ A 1 4 a 1 − 6 μ A 1 2 a 1 b 1 2 = 0 , Ψ : 12 μ 2 A 2 2 a 0 2 b 1 − 12 μ A 1 2 A 2 a 0 2 b 1 + 3 A 1 4 a 0 2 b 1 − 12 μ 2 ρ A 2 2 b 1 + 12 μ 2 A 2 a 0 b 1 2 + 12 μ ρ A 1 2 A 2 b 1 − 6 μ A 1 2 a 0 b 1 2 − 3 ρ A 1 4 b 1 − 24 μ 3 A 2 b 1 + 12 μ 2 A 1 2 b 1 + 4 μ 2 b 1 3 = 0 , Φ 0 : 4 μ 2 A 2 2 a 0 3 − 4 μ A 1 2 A 2 a 0 3 + A 1 4 a 0 3 − 12 μ 2 ρ A 2 2 a 0 + 12 μ ρ A 1 2 A 2 a 0 − 3 ρ A 1 4 a 0 = 0 .

Utilizing Maple 17 to solve the aforementioned system, it results in an impractical case, i.e., ρ = 0 . This leads to an absence of the independent variable t , which is not allowed for the proposed equation. Thus, this case does not provide any exact solution for problem (2).

3.2 Application of the Sardar subequation technique

In this section, explicit exact solutions of Eq. (2) are obtained using the Sardar subequation method. Based on the method described in Section 2.2.2, we assume that a solution form of Eq. (27) with N = 1 is

(39) U ( ϱ ) = ϖ 0 + ϖ 1 Q ( ϱ ) ,

where ϖ 0 and ϖ 1 ≠ 0 are the constant coefficients, which will be discovered later, and the function Q ( ϱ ) satisfies Eq. (19) whose analytical solutions are expressed in Eqs (20)–(23). Substituting Eq. (39) and its required derivatives into Eq. (27) along with the aid of Eq. (19) and then letting all coefficients of Q ( ϱ ) be zero, a system of nonlinear equations in ϖ 0 , ϖ 1 , a , and ρ is obtained as follows:

(40) Q 0 : − ρ ϖ 0 + 1 3 ϖ 0 3 = 0 , Q 1 : ϖ 0 2 ϖ 1 + 2 ϖ 1 a − ρ ϖ 1 = 0 , Q 2 : ϖ 0 ϖ 1 2 = 0 , Q 3 : 1 3 ϖ 1 3 + 4 ϖ 1 = 0 .

Solving system (40) using Maple 17, we obtain the following result:

(41) μ = μ , ρ = 2 a , ϖ 0 = 0 , ϖ 1 = ± 2 3 i ,

where μ and a are arbitrary constants. Using Eqs (20)–(23) depending on the conditions of a and μ in (19), the solution form (39), and the result (41), we obtain the explicit exact solutions of Eq. (2) as shown in the following.

Case 1: when a > 0 and μ = 0 , the solutions of Eq. (2) are

(42) u 1 ± ( x , y , t ) = ± 2 3 p q a sech p q ( a ϱ ) , u 2 ± ( x , y , t ) = ± 2 i 3 p q a csch p q ( a ϱ ) ,

where ϱ = Γ ( γ + 1 ) x β β + y β β − ρ t α α .

Case 2: When a < 0 and μ = 0 , the explicit exact solutions of Eq. (2) are

(43) u 3 ± ( x , y , t ) = ± 2 3 p q a sec p q ( − a ϱ ) , u 4 ± ( x , y , t ) = ± 2 3 p q a csc p q ( − a ϱ ) ,

where ϱ = Γ ( γ + 1 ) x β β + y β β − ρ t α α .

Case 3: When a < 0 and μ = a 2 4 , the solutions of Eq. (2) are

(44) u 5 ± ( x , y , t ) = ± 6 a tanh p q − a 2 ϱ , u 6 ± ( x , y , t ) = ± 6 a coth p q − a 2 ϱ , u 7 ± ( x , y , t ) = ± 6 a ( tanh p q ( − 2 a ϱ ) ± i p q sech p q ( − 2 a ϱ ) ) , u 8 ± ( x , y , t ) = ± 6 a ( coth p q ( − 2 a ϱ ) ± p q csch p q ( − 2 a ϱ ) ) , u 9 ± ( x , y , t ) = ± 3 2 a tanh p q − a 8 ϱ + coth p q − a 8 ϱ ,

where ϱ = Γ ( γ + 1 ) x β β + y β β − ρ t α α .

Case 4: When a > 0 and μ = a 2 4 , the explicit exact solutions of Eq. (2) are

(45) u 10 ± ( x , y , t ) = ± i 6 a tan p q a 2 ϱ , u 11 ± ( x , y , t ) = ± i 6 a cot p q a 2 ϱ , u 12 ± ( x , y , t ) = ± i 6 a ( tan p q ( 2 a ϱ ) ± p q sec p q ( 2 a ϱ ) ) , u 13 ± ( x , y , t ) = ± i 6 a ( cot p q ( 2 a ϱ ) ± p q csc p q ( 2 a ϱ ) ) , u 14 ± ( x , y , t ) = ± i 3 2 a tan p q a 8 ϱ − cot p q a 8 ϱ ,

where ϱ = Γ ( γ + 1 ) x β β + y β β − ρ t α α .

4 Bifurcation analysis and chaotic behavior

This section is devoted to analyzing the bifurcations [24,53,54] of Eq. (27). To achieve this, we employ the Galilean transformation to transform Eq. (27) to an equivalent dynamical system. Later, we will analyze the phase portraits of the corresponding dynamical system whose behaviors can be studied using MATLAB. Now, defining χ ( ϱ ) = U ( ϱ ) and Λ ( ϱ ) = U ′ ( ϱ ) , we obtain

(46) χ ′ ( ϱ ) = Λ ( ϱ ) , Λ ′ ( ϱ ) = ρ 2 χ ( ϱ ) − 1 6 χ 3 ( ϱ ) .

The phase portrait analysis depends on the dynamical behavior of system (46). Moreover, on the basis of the phase plane, we can plot different trajectories to show the potential state of the aforementioned system. We proceed to derive the auxiliary Hamiltonian function, which is considered as the first integral of system (46) given by

(47) H ( χ , Λ ) = 1 2 Λ 2 − ρ 4 χ 2 + 1 24 χ 4 = w ,

where w is the Hamiltonian constant of motion. The function (47) is regarded as the energy of a particle with an effective mass m ( χ ) = 1 moving with a potential U ˆ ( χ ) = − ρ 4 χ 2 + 1 24 χ 4 . It can be observed that the stationary behavior of system (46) will depend on the value of the parameter ρ initially in the wave variable determining the time of the wave propagation. Now, we are going to analyze the phase portraits of system (46) based on its equilibrium points. The system admits the following equilibrium points: Ξ 1 = ( 0 , 0 ) , Ξ 2 = ( 3 ρ , 0 ) , and Ξ 3 = ( − 3 ρ , 0 ) . Therefore, the Jacobian of system (46) at equilibrium point Ξ ( χ , 0 ) is

(48) J ( χ , Λ ) = 0 1 ρ 2 − 1 2 χ 2 ( ϱ ) 0 ,

whose determinant is

(49) ∣ J ( χ , Λ ) ∣ = − ρ 2 + 1 2 χ 2 ( ϱ ) .

It is known from the theory of planar dynamical systems that we can analyze the equilibrium point Ξ ( χ , 0 ) as follows:

If ∣ J ( χ , Λ ) ∣ > 0 , then ( χ , Λ ) is a center point.
If ∣ J ( χ , Λ ) ∣ = 0 , then ( χ , Λ ) is a degenerate point.
If ∣ J ( χ , Λ ) ∣ < 0 , then ( χ , Λ ) is a saddle point.

Clearly, based on the conditions outlined earlier and the presence of several equilibrium points, one can now analyze the bifurcation of the phase portrait. The system orbit determines which solitary waves propagate in the governing equation. When two equilibrium points are formed from a heteroclinic orbit, the nonlinear system may generate kink or anti-kink solitary waves. In contrast, shallow (gray) solitons can arise in the governing equation if the phase portrait depicts a homoclinic orbit from a point other than the origin. The behavior of systems is solely dependent on the value of the parameter ρ . If we set ρ = 1 , then three equilibrium points emerge: ( 0 , 0 ) , ( 3 , 0 ) , and ( − 3 , 0 ) as shown in Figure 1(a). The heteroclinic orbits are identified in the system to confirm the propagation of kink or anti-kink solitary waves. Furthermore, the scenario in Figure 1(b) reveals that the point ( 0 , 0 ) represents a saddle point.

$Figure 1 (a) and (b) Bifurcation of phase portrait for system (46) by varying the parameter ρ \rho .$

Figure 1

(a) and (b) Bifurcation of phase portrait for system (46) by varying the parameter ρ .

Next, we aim to investigate the chaotic behavior of system (46). To achieve that, we have to add a periodic external force to system (46) to become

(50) χ ′ ( ϱ ) = Λ ( ϱ ) , Λ ′ ( ϱ ) = ρ 2 χ ( ϱ ) − 1 6 χ 3 ( ϱ ) + ε cos ( υ t ) ,

where υ is the frequency of the system and ε is the amplitude of the external acting force. It is obvious that for ε = 0 , the dynamical system (50) corresponds to the unperturbed system (46). In Figures 2 and 3, we utilized the same values as in Figure 1, and we selected the value of the external force amplitude ε = 0.75 to display time-domain waveforms. The waveforms exhibit alternate splitting and resting for υ = 0.1 , followed by modest oscillation at 0.2. This indicates that the magnitude of the externally used perturbation causes chaotic behavior in the nonlinear system. Adjusting the external force amplitude to ε = 0.25 , Figure 4(a)--(c) shows a similar perturbed nonlinear system, indicating chaotic and oscillatory dynamics for the system with frequency υ = 0.5 and ρ = − 5 . In Figure 5, we vary the external force amplitude to 0.95. The provided figures show chaotic behaviors, confirming that the nonlinear system can have solitary and periodic wave solutions for various homoclinic and heteroclinic orbits. The system is very sensitive to changes in the parameters ε and υ , offering substantial insight into the impact of the perturbation term ε cos ( υ t ) on the entire system behavior.

$Figure 2 (a) 2D phase portrait, (b) 3D phase portrait, and (c) time-dependent waveforms of (50) with external force exhibit chaotic behavior for ρ = 1 \rho =1 , ε = 0.75 \varepsilon =0.75 , and υ = 0.1 \upsilon =0.1 .$

Figure 2

(a) 2D phase portrait, (b) 3D phase portrait, and (c) time-dependent waveforms of (50) with external force exhibit chaotic behavior for ρ = 1 , ε = 0.75 , and υ = 0.1 .

$Figure 3 (a) 2D phase portrait, (b) 3D phase portrait, and (c) time-dependent waveforms of (50) with external force exhibit chaotic behavior for ρ = 1 \rho =1 , ε = 0.75 \varepsilon =0.75 , and υ = 0.2 \upsilon =0.2 .$

Figure 3

(a) 2D phase portrait, (b) 3D phase portrait, and (c) time-dependent waveforms of (50) with external force exhibit chaotic behavior for ρ = 1 , ε = 0.75 , and υ = 0.2 .

$Figure 4 (a) 2D phase portrait, (b) 3D phase portrait, and (c) time-dependent waveforms of (50) with external force exhibit chaotic behavior for ρ = − 5 \rho =-5 , ε = 0.25 \varepsilon =0.25 , and υ = 0.5 \upsilon =0.5 .$

Figure 4

(a) 2D phase portrait, (b) 3D phase portrait, and (c) time-dependent waveforms of (50) with external force exhibit chaotic behavior for ρ = − 5 , ε = 0.25 , and υ = 0.5 .

$Figure 5 (a) 2D phase portrait, (b) 3D phase portrait, and (c) time-dependent waveforms of (50) with external force exhibit chaotic behavior for ρ = − 5 \rho =-5 , ε = 0.95 \varepsilon =0.95 , and υ = 0.5 \upsilon =0.5 .$

Figure 5

(a) 2D phase portrait, (b) 3D phase portrait, and (c) time-dependent waveforms of (50) with external force exhibit chaotic behavior for ρ = − 5 , ε = 0.95 , and υ = 0.5 .

5 Graphical results

In this section, some interesting graphical plots of selected explicit exact solutions of (2), derived via the ( G ′ ∕ G , 1 ∕ G ) -expansion and Sardar subequation methods, are given. All graphs were plotted using the Maple 17 software package. The fractional orders in the graphs are as follows: α = 0.7 , 0.8 , 0.9 and β = 0.9 .

5.1 Solution graphs generated using the ( G ′ ∕ G , 1 ∕ G ) -expansion technique

In this section, we show plots obtained from the ( G ′ ∕ G , 1 ∕ G ) -expansion method, and we demonstrate plots of 3D, 2D, and contour of solutions (32) and (36) of (2). In Figure 6, plots of u ( x , y , t ) in (32) are shown for D 1 = { x ∈ [ 0 , 5 ] , y = 1 , t ∈ [ 0 , 5 ] } and D 2 = { x ∈ [ 0 , 5 ] , y = 1 , t = 1 } for the 3D and 2D plots, respectively. Contour plots for solution (32), representing a 3D surface by plotting ( x , t ) contours for some fixed values of u , are plotted for y = 1 . The following parameter values: λ = − 1 , μ = 0 , a 0 = 0 , a 1 = 0 , A 1 = 3 , A 2 = 5 , and γ = 1.5 are used for plotting in the figure. Particularly, Figure 6(a)–(c), (d)–(f), and (g)–(i) displays the graphs of 3D, 2D, and contour of solution (32) evaluated at { α = β = 0.9 } , { α = 0.8 , β = 0.9 } , and { α = 0.7 , β = 0.9 } , respectively. As can be observed from the 3D plots of Figure 6, the physical behavior of solution (32) is a bell-shaped solution or a bright-soliton wave. Figure 7 shows the effect of the time-fractional derivative on the solution (32) when the time-fractional order α is varied ranging from 0.7 to 0.9.

$Figure 6 Graphs for u + ( x , y , t ) {u}^{+}\left(x,y,t) in (32) by the ( G ′ ∕ G , 1 ∕ G ) \left(G^{\prime} /G,1/G) -expansion technique: (a)–(c) when β = α = 0.9 \beta =\alpha =0.9 ; (d)–(f) when β = 0.9 \beta =0.9 and α = 0.8 \alpha =0.8 ; (g)–(i) when β = 0.9 \beta =0.9 and α = 0.7 \alpha =0.7 .$

Figure 6

Graphs for u + ( x , y , t ) in (32) by the ( G ′ ∕ G , 1 ∕ G ) -expansion technique: (a)–(c) when β = α = 0.9 ; (d)–(f) when β = 0.9 and α = 0.8 ; (g)–(i) when β = 0.9 and α = 0.7 .

$Figure 7 Effect of the time-fractional order α \alpha on the solution u + ( x , y , t ) {u}^{+}\left(x,y,t) in (32) obtained by the ( G ′ ∕ G , 1 ∕ G ) \left(G^{\prime} /G,1/G) -expansion technique when α = 0.9 \alpha =0.9 , α = 0.8 \alpha =0.8 , and α = 0.7 \alpha =0.7 .$

Figure 7

Effect of the time-fractional order α on the solution u + ( x , y , t ) in (32) obtained by the ( G ′ ∕ G , 1 ∕ G ) -expansion technique when α = 0.9 , α = 0.8 , and α = 0.7 .

Figure 8 shows the magnitudes of u ( x , y , t ) in (36), which are drawn in 3D for D 3 = { x ∈ [ − 10 , 10 ] , y = 1 , t ∈ [ − 10 , 10 ] } , and in 2D for D 4 = { x ∈ [ − 10 , 10 ] , y = 1 , t = 1 } . Contours of the magnitude of solution (36) are plotted for y = 1 . The following parameter values: λ = 1 , μ = 1 , a 0 = 0 , a 1 = 3 i , A 1 = 3 , A 2 = 5 , and γ = 1.5 are utilized for plotting in Figure 8. In particular, Figure 8(a)–(c), (d)–(f), and (g)–(i) shows the magnitude graphs of 3D, 2D, and contour for solution (36) computed at { α = β = 0.9 } , { α = 0.8 , β = 0.9 } , and { α = 0.7 , β = 0.9 } , respectively. The magnitude of solution (36) as shown in the three-dimensional (3D) graphs of Figure 8 can be physically interpreted as a multiple-soliton wave. Figure 9 illustrates the influence of the time-fractional order α on the solution (36), with values ranging from 0.7 to 0.9.

$Figure 8 Graphs for the magnitude of u + {u}^{+} in (36) by the ( G ′ ∕ G , 1 ∕ G ) \left(G^{\prime} /G,1/G) -expansion approach: (a)–(c) when β = α = 0.9 \beta =\alpha =0.9 ; (d)–(f) when β = 0.9 \beta =0.9 and α = 0.8 \alpha =0.8 ; (g)–(i) when β = 0.9 \beta =0.9 and α = 0.7 \alpha =0.7 .$

Figure 8

Graphs for the magnitude of u + in (36) by the ( G ′ ∕ G , 1 ∕ G ) -expansion approach: (a)–(c) when β = α = 0.9 ; (d)–(f) when β = 0.9 and α = 0.8 ; (g)–(i) when β = 0.9 and α = 0.7 .

$Figure 9 Effect of the time-fractional order α \alpha on ∣ u + ( x , y , t ) ∣ | {u}^{+}\left(x,y,t)| , where u + ( x , y , t ) {u}^{+}\left(x,y,t) in (36) obtained by the ( G ′ ∕ G , 1 ∕ G ) \left(G^{\prime} /G,1/G) -expansion technique when α = 0.9 \alpha =0.9 , α = 0.8 \alpha =0.8 , and α = 0.7 \alpha =0.7 .$

Figure 9

Effect of the time-fractional order α on ∣ u + ( x , y , t ) ∣ , where u + ( x , y , t ) in (36) obtained by the ( G ′ ∕ G , 1 ∕ G ) -expansion technique when α = 0.9 , α = 0.8 , and α = 0.7 .

5.2 Solution graphs generated using the Sardar subequation scheme

In this section, the results from (2) derived by the Sardar subequation method are plotted for selected solutions, namely, u 1 + in Eq. (42), u 4 + in Eq. (43), u 5 + in Eq. (44), and u 12 + in Eq. (45). Graphs for 3D, 2D, and contours are plotted for each solution. Graphical plots of u 1 + in (42) are plotted in Figure 10 on D 5 = { 0 ⩽ x ⩽ 20 , y = 1 , 0 ⩽ t ⩽ 20 } for 3D graphs and on D 6 = { 0 ⩽ x ⩽ 20 , y = 1 , t = 1 } for 2D plots. Solution contours of (42), representing a 3D surface by plotting ( x , t ) contours for some fixed values of u 1 , are plotted for y = 1 . The following parameter values: a = 1 , μ = 0 , ϖ 0 = 0 , ϖ 1 = 2 3 i , p = 1.2 , q = 1.4 , and γ = 1.5 are used for the graphs in the figure. Particularly, Figure 10(a)–(c) shows the plots for 3D, 2D, and contour of (42) when { α = β = 0.9 } . Similarly, Figure 10(d)–(f) and (g)–(i) shows the assorted graphs of the solution for { α = 0.8 , β = 0.9 } and { α = 0.7 , β = 0.9 } , respectively. As can be observed from the 3D displays of Figure 10, the physical behavior of solution (42) can be interpreted as an anti-bell-shaped soliton solution or a dark-soliton wave. Figure 11 illustrates the influence of the time-fractional order α on solution (42), with values ranging from 0.7 to 0.9 .

$Figure 10 Graphs for u 1 + ( x , y , t ) {u}_{1}^{+}\left(x,y,t) in (42) by the Sardar subequation method: (a)–(c) when β = α = 0.9 \beta =\alpha =0.9 ; (d)–(f) when β = 0.9 \beta =0.9 and α = 0.8 \alpha =0.8 ; (g)–(i) when β = 0.9 \beta =0.9 and α = 0.7 \alpha =0.7 .$

Figure 10

Graphs for u 1 + ( x , y , t ) in (42) by the Sardar subequation method: (a)–(c) when β = α = 0.9 ; (d)–(f) when β = 0.9 and α = 0.8 ; (g)–(i) when β = 0.9 and α = 0.7 .

$Figure 11 Effect of the time-fractional order α \alpha on the solution u 1 + ( x , y , t ) {u}_{1}^{+}\left(x,y,t) in (42) obtained by the Sardar subequation method when α = 0.9 \alpha =0.9 , α = 0.8 \alpha =0.8 , and α = 0.7 \alpha =0.7 .$

Figure 11

Effect of the time-fractional order α on the solution u 1 + ( x , y , t ) in (42) obtained by the Sardar subequation method when α = 0.9 , α = 0.8 , and α = 0.7 .

Magnitudes of u 4 + ( x , y , t ) in (43) for parameter values D 7 = { 0 ⩽ x ⩽ 5 , y = 1 , 0 ⩽ t ⩽ 5 } for the 3D plots and D 8 = { 0 ⩽ x ⩽ 5 , y = 1 , t = 1 } for the 2D plots are shown in Figure 12. The contours of magnitude of u 4 + are portrayed with y = 1 . The following parameter values: a = − 1 , μ = 0 , ϖ 0 = 0 , ϖ 1 = 2 3 i , p = 1.2 , q = 1.4 , and γ = 1.5 are used for plotting in Figure 12. Particularly, Figure 12(a)–(c) shows the magnitude graphs for 3D, 2D, and contour of (43) calculated when { α = β = 0.9 } . Similarly, assorted magnitude graphs of this solution computed at { α = 0.8 , β = 0.9 } and { α = 0.7 , β = 0.9 } are shown in Figure 12(d)–(f) and (g)–(i), respectively. As can be seen from the 3D graphs of Figure 12, the physical performance of u 4 + behaves as a periodical breather-type wave solution. Figure 13 shows the effect of the time-fractional derivative on the magnitude of solution (43) when the time-fractional order α is varied ranging from 0.7 to 0.9.

$Figure 12 Graphs for ∣ u 4 + ( x , y , t ) ∣ | {u}_{4}^{+}\left(x,y,t)| where u 4 + {u}_{4}^{+} is expressed in (43) by the Sardar subequation technique: (a)–(c) when β = α = 0.9 \beta =\alpha =0.9 ; (d)–(f) when β = 0.9 \beta =0.9 and α = 0.8 \alpha =0.8 ; (g)–(i) when β = 0.9 \beta =0.9 and α = 0.7 \alpha =0.7 .$

Figure 12

Graphs for ∣ u 4 + ( x , y , t ) ∣ where u 4 + is expressed in (43) by the Sardar subequation technique: (a)–(c) when β = α = 0.9 ; (d)–(f) when β = 0.9 and α = 0.8 ; (g)–(i) when β = 0.9 and α = 0.7 .

$Figure 13 Effect of the time-fractional order α \alpha on ∣ u 4 + ( x , y , t ) ∣ | {u}_{4}^{+}\left(x,y,t)| , where u 4 + {u}_{4}^{+} in (43) obtained by the Sardar subequation technique when α = 0.9 \alpha =0.9 , α = 0.8 \alpha =0.8 , and α = 0.7 \alpha =0.7 .$

Figure 13

Effect of the time-fractional order α on ∣ u 4 + ( x , y , t ) ∣ , where u 4 + in (43) obtained by the Sardar subequation technique when α = 0.9 , α = 0.8 , and α = 0.7 .

Magnitudes of u 5 + in (44) are plotted for the 3D and 2D graphs on D 9 = { x ∈ [ − 10 , 10 ] , y = 1 , t ∈ [ − 10 , 10 ] } and on D 10 = { x ∈ [ − 10 , 10 ] , y = 1 , t = 1 } , respectively, as shown in Figure 14. Contours of magnitude of u 5 + are portrayed for y = 1 . The following parameter values: a = − 1 , ϖ 0 = 0 , ϖ 1 = 2 3 i , p = 1.2 , q = 1.4 , and γ = 1.5 are used for plotting in Figure 14. Particularly, Figure 14(a)–(c), (d)–(f), and (g)–(i) shows the assorted graphs of u 5 + calculated at { α = β = 0.9 } , { α = 0.8 , β = 0.9 } , and { α = 0.7 , β = 0.9 } , respectively. As can be seen from the 3D graphs of Figure 14 displaying the magnitude of u 5 + in (44), its physical behavior can be characterized as a singular multiple-soliton. Figure 15 illustrates the influence of the time-fractional derivative on the magnitude of solution (44) when the time-fractional order α is varied ranging from 0.7 to 0.9.

$Figure 14 Graphs for ∣ u 5 + ( x , y , t ) ∣ | {u}_{5}^{+}\left(x,y,t)| where u 5 + {u}_{5}^{+} is shown in (44) by the Sardar subequation method: (a)–(c) when β = α = 0.9 \beta =\alpha =0.9 ; (d)–(f) when β = 0.9 \beta =0.9 and α = 0.8 \alpha =0.8 ; (g)–(i) when β = 0.9 \beta =0.9 and α = 0.7 \alpha =0.7 .$

Figure 14

Graphs for ∣ u 5 + ( x , y , t ) ∣ where u 5 + is shown in (44) by the Sardar subequation method: (a)–(c) when β = α = 0.9 ; (d)–(f) when β = 0.9 and α = 0.8 ; (g)–(i) when β = 0.9 and α = 0.7 .

$Figure 15 Effect of the time-fractional order α \alpha on ∣ u 5 + ( x , y , t ) ∣ | {u}_{5}^{+}\left(x,y,t)| , where u 5 + {u}_{5}^{+} in (44) obtained by the Sardar subequation technique when α = 0.9 \alpha =0.9 , α = 0.8 \alpha =0.8 , and α = 0.7 \alpha =0.7 .$

Figure 15

Effect of the time-fractional order α on ∣ u 5 + ( x , y , t ) ∣ , where u 5 + in (44) obtained by the Sardar subequation technique when α = 0.9 , α = 0.8 , and α = 0.7 .

In Figure 16, magnitudes of u 12 + in (45) are plotted for the 3D and 2D plots on different domains, i.e., D 11 = { x ∈ [ − 5 , 5 ] , y = 1 , t ∈ [ − 5 , 5 ] } and D 12 = { x ∈ [ − 5 , 5 ] , y = 1 , t = 1 } . The magnitude plots of its contour are shown for y = 1 . The following parameter values: a = 1 , ϖ 0 = 0 , ϖ 1 = 2 3 i , p = 1.2 , q = 1.4 , and γ = 1.5 are used for plotting in Figure 16. Specifically, Figure 16(a)–(c) illustrates 3D, 2D, and contour graphs of ∣ u 12 + ∣ evaluated at { α = β = 0.9 } . Similarly, Figure 16(d)–(f) and (g)–(i) shows the assorted magnitude plots of this solution for { α = 0.8 , β = 0.9 } and { α = 0.7 , β = 0.9 } , respectively. According to 3D plots of this figure, u 12 + in (45) shows the behavior of a singularly periodic wave or a multiple-soliton wave. Figure 17 shows the effect of the time-fractional derivative on the magnitude of solution (45) when the time-fractional order α is varied ranging from 0.7 to 0.9.

$Figure 16 Graphs for ∣ u 12 + ( x , y , t ) ∣ | {u}_{12}^{+}\left(x,y,t)| where u 12 + {u}_{12}^{+} is shown in (45) by the Sardar subequation method: (a)–(c) when β = α = 0.9 \beta =\alpha =0.9 ; (d)–(f) when β = 0.9 \beta =0.9 and α = 0.8 \alpha =0.8 ; (g)–(i) when β = 0.9 \beta =0.9 and α = 0.7 \alpha =0.7 .$

Figure 16

Graphs for ∣ u 12 + ( x , y , t ) ∣ where u 12 + is shown in (45) by the Sardar subequation method: (a)–(c) when β = α = 0.9 ; (d)–(f) when β = 0.9 and α = 0.8 ; (g)–(i) when β = 0.9 and α = 0.7 .

$Figure 17 Effect of the time-fractional order α \alpha on ∣ u 12 + ( x , y , t ) ∣ | {u}_{12}^{+}\left(x,y,t)| , where u 12 + {u}_{12}^{+} in (45) obtained by the Sardar subequation technique when α = 0.9 \alpha =0.9 , α = 0.8 \alpha =0.8 , and α = 0.7 \alpha =0.7 .$

Figure 17

Effect of the time-fractional order α on ∣ u 12 + ( x , y , t ) ∣ , where u 12 + in (45) obtained by the Sardar subequation technique when α = 0.9 , α = 0.8 , and α = 0.7 .

From all graphs shown previously, it can be seen that varying the fractional orders α and β typically affects only their amplitude or translation, but not their overall graphical structure.

6 Discussion and conclusions

The major findings of this article are to obtain a variety of novel exact traveling wave solutions of (2) via the ( G ′ ∕ G , 1 ∕ G ) -expansion and Sardar subequation techniques and to investigate bifurcation and chaotic behaviors for the model problem. Since Eq. (2) is new, all the results derived here have been presented for the first time. Utilizing these two methods along with the assistance of Maple 17, we have successfully obtained exact solutions of the proposed equation, which have been displayed in terms of the wave transformation ϱ in (25) and the fractional orders 0 < α , β ≤ 1 . Specifically, the solutions of the equation derived utilizing the ( G ′ ∕ G , 1 ∕ G ) -expansion method have been expressed as hyperbolic and trigonometric functions, but not as a rational function. Applying this method to (2), we obtain five solutions as mentioned in Section 3.1. The Sardar subequation method has given explicit exact solutions of the problem with regard to the special generalized hyperbolic and trigonometric functions. It is worth noting that those special solutions can be reduced to regular hyperbolic and trigonometric functions by selecting p = q = 1 . As described in Section 3.2, the Sardar subequation method provides 14 exact traveling solutions for (2). Bifurcation analysis has been carried out by converting the governing equation into a dynamical system. The chaotic behaviors of the dynamical system with a perturbed term have been classified via the phase portraits. In addition, the plots of 3D, 2D, and contour graphs of selected solutions for a range of values of fractional orders α and β have been plotted using Maple 17 in order to study the effects on the physical behavior of selected solutions. As a result, the solutions of (2), selected for plotting in this article, provide examples of different physical behaviors, namely, a bell-shaped solution, a singularly periodic wave, and a singular multiple-soliton. It is worth noting that the smoothness of the obtained solutions in this study depends on their singularities. For convenience, when values of the parameters are set in the solutions, their singularities can be determined at which the smoothness of those solutions is absent. Furthermore, the influence of varying values of the fractional orders of the proposed Eq. (2) has been investigated in Figures 7, 9, 11, 13, 15, and 17. For each of these figures, varying the time-fractional order α primarily affects the amplitude or translation of the solution graphs such as the shifting positions of singularities, while the overall graphical structure remains quite unchanged.

With the aid of Maple 17, all of the obtained exact solutions have been checked by substituting them back in the original problem. According to the results reported here, we have shown that both of the methods used here along with the aid of symbolic software such as Maple are fruitful, efficient, and reliable mathematical methods for generating exact traveling wave solutions of NPDEs. It has been recommended in future that the ( 3 + 1 ) -dimensional mZK equation [55] could be solved using these two methods because they can produce several types of explicit exact solutions. In fact, the ( G ′ ∕ G , 1 ∕ G ) -expansion method can be reduced to the ( G ′ ∕ G ) -expansion method so solutions derived by the ( G ′ ∕ G , 1 ∕ G ) -expansion method recover the solutions generated employing the ( G ′ ∕ G ) -expansion method [56]. Moreover, we have seen that applying the methods to the truncated M-fractional also gives good solutions. We therefore believe that the methods discussed in this article could be successfully applied to a range of NPDEs with fractional derivatives of various types.

Funding information: The authors greatly appreciate referees for their many significant comments and valuable suggestions. Sekson Sirisubtawee was supported by the Faculty of Applied Science, King Mongkut s University of Technology North Bangkok (Grant No. 672168).
Author contributions: Conceptualization of this study, P.W., S.S.; methodology, P.M., P.W.; formal analysis, P.W., S.S.; investigation, P.M., P.W., S.S.; software, P.W., S.S.; visualization, P.M., P.W., S.S.; validation, P.W., S.S.; data curation, S.S.; writing–original draft preparation, P.M., P.W., S.S.; writing–review and editing, P.M., P.W., S.S.; supervision project administration, S.S.; funding acquisition, S.S. All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
Conflict of interest: The authors state no conflict of interest.
Data availability statement: All data generated or analysed during this study are included in this published article.

References

[1] Torvattanabun M, Khansai N, Sirisubtawee S, Koonprasert S, Tuan NM. New exact traveling wave solutions of the (3.1)-dimensional Chiral nonlinear Schrödinger equation using two reliable techniques. Thai J Math. 2024;22(1):145–63. Search in Google Scholar

[2] Tuan NM, Kooprasert S, Sirisubtawee S, Meesad P. The bilinear neural network method for solving Benney-Luke equation. Partial Differ Equ Appl Math. 2024;10:100682. 10.1016/j.padiff.2024.100682Search in Google Scholar

[3] Boakye G, Hosseini K, Hincccal E, Sirisubtawee S, Osman M. Some models of solitary wave propagation in optical fibers involving Kerr and parabolic laws. Optical Quantum Electronics. 2024;56(3):345. 10.1007/s11082-023-05903-5Search in Google Scholar

[4] Alqudah M, Shah SZH, Riaz MB, Khalifa HAEW, Akgül A, Ayub A. Neural network architecture to optimize the nanoscale thermal transport of ternary magnetized Carreau nanofluid over 3D wedge. Results Phys. 2024;59:107616. 10.1016/j.rinp.2024.107616Search in Google Scholar

[5] Manukure S, Booker T. A short overview of solitons and applications. Partial Differ Equ Appl Math. 2021;4: 100140. 10.1016/j.padiff.2021.100140Search in Google Scholar

[6] Sabíu J, Sirisubtawee S, Inc M. Optical soliton solutions for the Chavy-Waddy-Kolokolnikov model for bacterial colonies using two improved methods. J Appl Math Comput. 2024;70:1–24. 10.1007/s12190-024-02169-2Search in Google Scholar

[7] Rezazadeh H. New solitons solutions of the complex Ginzburg-Landau equation with Kerr law nonlinearity. Optik. 2018;167:218–27. 10.1016/j.ijleo.2018.04.026Search in Google Scholar

[8] Akbar MA, Akinyemi L, Yao SW, Jhangeer A, Rezazadeh H, Khater MM, et al. Soliton solutions to the Boussinesq equation through sine-Gordon method and Kudryashov method. Results Phys. 2021;25: 104228. 10.1016/j.rinp.2021.104228Search in Google Scholar

[9] León-Ramírez A, González-Gaxiola O, Chacón-Acosta G. Analytical solutions to the Chavy-Waddy-Kolokolnikov model of bacterial aggregates in phototaxis by three integration schemes. Mathematics. 2023;11(10):2352. 10.3390/math11102352Search in Google Scholar

[10] Sirisubtawee S, Koonprasert S, Sungnul S. Some applications of the (G′∕G,1∕G)-expansion method for finding exact traveling wave solutions of nonlinear fractional evolution equations. Symmetry. 2019;11(8):952. 10.3390/sym11080952Search in Google Scholar

[11] Pleumpreedaporn C, Pleumpreedaporn S, Moore EJ, Sirisubtawee S, Sungnul S. Novel exact traveling wave solutions for the (2+1)-dimensional Boiti-Leon-Manna-Pempinelli equation with Atangana’s space and time beta-derivatives via the Sardar subequation method: Annual Meeting in Mathematics 2023. Thai J Math. 2024;22(1):1–18. Search in Google Scholar

[12] Zakharov V, Kuznetsov EA. Three-dimensional solutions. Soviet Phys Uspekhi. 1974 Jan;39:285–6. Search in Google Scholar

[13] Kadomstev BB, Petviashvilli VI. On the stability of solitary waves in weakly dispersing media. Soviet Phys Doklady. 1970;15:539. Search in Google Scholar

[14] Mushtaq A, Shah H. Nonlinear Zakharov-Kuznetsov equation for obliquely propagating two-dimensional ion-acoustic solitary waves in a relativistic, rotating magnetized electron-positron-ion plasma. Phys Plasmas. 2005;12(7):072306. 10.1063/1.1946729Search in Google Scholar

[15] Biswas A, Zerrad E. Solitary wave solution of the Zakharov-Kuznetsov equation in plasmas with power law nonlinearity. Nonl Anal Real World Appl. 2010;11(4):3272–4. 10.1016/j.nonrwa.2009.08.007Search in Google Scholar

[16] Ray SS, Sahoo S. New exact solutions of fractional Zakharov-Kuznetsov and modified Zakharov-Kuznetsov equations using fractional sub-equation method. Commun Theoret Phys. 2015;63(1):25. 10.1088/0253-6102/63/1/05Search in Google Scholar

[17] Al-Amin M, Islam MN, Ilhan OA, Akbar MA, Soybaş D, et al. Solitary wave solutions to the modified Zakharov-Kuznetsov and the (2+1)-dimensional Calogero-Bogoyavlenskii-Schiff models in mathematical physics. J Math. 2022;2022(1):5224289. 10.1155/2022/5224289Search in Google Scholar

[18] Afrin FU. Solitary wave solutions and investigation the effects of different wave velocities of the nonlinear modified Zakharov-Kuznetsov model for the wave propagation in nonlinear media. Partial Differ Equ Appl Math. 2023;8:100583. 10.1016/j.padiff.2023.100583Search in Google Scholar

[19] Alam MN, Tunccc C. Constructions of the optical solitons and other solitons to the conformable fractional Zakharov-Kuznetsov equation with power law nonlinearity. J Taibah Univ Sci. 2020;14(1):94–100. 10.1080/16583655.2019.1708542Search in Google Scholar

[20] Munro S, Parkes E. The derivation of a modified Zakharov-Kuznetsov equation and the stability of its solutions. J Plasma Phys. 1999;62(3):305–17. 10.1017/S0022377899007874Search in Google Scholar

[21] Munro S, Parkes E. Stability of solitary-wave solutions to a modified Zakharov-Kuznetsov equation. J Plasma Phys. 2000;64(4):411–26. 10.1017/S0022377800008771Search in Google Scholar

[22] Zhang XZ, Siddique I, Mehdi KB, Elmandouh A, Inc M. Novel exact solutions, bifurcation of nonlinear and supernonlinear traveling waves for M-fractional generalized reaction Duffing model and the density dependent M-fractional diffusion reaction equation. Results Phys. 2022;37:105485. 10.1016/j.rinp.2022.105485Search in Google Scholar

[23] Siddique I, Mehdi KB, Jaradat MM, Zafar A, Elbrolosy ME, Elmandouh AA, et al. Bifurcation of some new traveling wave solutions for the time-space M-fractional MEW equation via three altered methods. Results Phys. 2022;41:105896. 10.1016/j.rinp.2022.105896Search in Google Scholar

[24] Han T, Zhang K, Jiang Y, Rezazadeh H. Chaotic pattern and solitary solutions for the (2+1)-dimensional beta-fractional double-chain DNA system. Fractal Fractional. 2024;8(7):415. 10.3390/fractalfract8070415Search in Google Scholar

[25] Mehdi KB, Mehdi Z, Samreen S, Siddique I, Elmandouh AA, Elbrolosy ME, et al. Novel exact traveling wave solutions of the space-time fractional Sharma Tasso-Olver equation via three reliable methods. Partial Differ Equ Appl Math. 2024;11:100784. 10.1016/j.padiff.2024.100784Search in Google Scholar

[26] Siddique I, Jaradat M, Zafar A, Mehdi KB, Osman M. Exact traveling wave solutions for two prolific conformable M-Fractional differential equations via three diverse approaches. Results Phys. 2021;28:104557. 10.1016/j.rinp.2021.104557Search in Google Scholar

[27] Guo M, Dong H, Liu J, Yang H. The time-fractional mZK equation for gravity solitary waves and solutions using sech-tanh and radial basic function method. Nonl Anal Model Control. 2019;24(1):1–19. 10.15388/NA.2019.1.1Search in Google Scholar

[28] Bibi A, Shakeel M, Khan D, Hussain S, Chou D. Study of solitary and kink waves, stability analysis, and fractional effect in magnetized plasma. Results Phys. 2023;44:106166. 10.1016/j.rinp.2022.106166Search in Google Scholar

[29] Fu L, Chen Y, Yang H. Time-space fractional coupled generalized Zakharov-Kuznetsov equations set for rossby solitary waves in two-layer fluids. Mathematics. 2019;7(1):41. 10.3390/math7010041Search in Google Scholar

[30] Baleanu D, Inç M, Yusuf A, Aliyu AI. Lie symmetry analysis, exact solutions and conservation laws for the time fractional modified Zakharov-Kuznetsov equation. Nonlinear Anal Model Control. 2017;22(6):861–76. 10.15388/NA.2017.6.9Search in Google Scholar

[31] Gorenflo R, Kilbas AA, Mainardi F, Rogosin SV. Mittag-Leffler functions, related topics and applications. Berlin: Springer; 2020. 10.1007/978-3-662-61550-8Search in Google Scholar

[32] Mittag-Leffler GM. Sur la nouvelle fonction Ealpha(x). CR Acad Sci Paris. 1903;137(2):554–8. Search in Google Scholar

[33] Sousa J, de Oliveira EC. A new truncated M-fractional derivative type unifying some fractional derivative types with classical properties. Int J Anal Appl. 2018;16(1):83–96. Search in Google Scholar

[34] Riaz MB, Awrejcewicz J, Jhangeer A. Optical solitons with beta and M-truncated derivatives in nonlinear negative-index materials with Bohm potential. Materials. 2021;14(18):5335. 10.3390/ma14185335Search in Google Scholar PubMed PubMed Central

[35] Yao SW, Manzoor R, Zafar A, Inc M, Abbagari S, Houwe A. Exact soliton solutions to the Cahn-Allen equation and Predator-Prey model with truncated M-fractional derivative. Results Phys. 2022;37:105455. 10.1016/j.rinp.2022.105455Search in Google Scholar

[36] Chauhan RB, Chudasama MH. A study of the left local general truncated M-fractional derivative. Appl Math E-Notes. 2023;23:100–23. Search in Google Scholar

[37] Onder I, Cinar M, Secer A, Bayram M. Analytical solutions of simplified modified Camassa-Holm equation with conformable and M-truncated derivatives: A comparative study. J Ocean Eng Sci. 2024;9(3):240–50. 10.1016/j.joes.2022.06.012Search in Google Scholar

[38] Raheel M, Razzaq W, Alsharidi AK, Zafar A. Exact solitons to M-fractional (2.1)-dimensional CNLSE based on three different methods. Results Phys. 2022;42:105983. 10.1016/j.rinp.2022.105983Search in Google Scholar

[39] Mohammed WW, Cesarano C, Al-Askar FM. Solutions to the (4+1)-dimensional time-fractional Fokas equation with M-truncated derivative. Mathematics. 2022;11(1):194. 10.3390/math11010194Search in Google Scholar

[40] Roshid HO, Roshid MM, Hossain MM, Hasan MS, Munshi MJH, Sajib AH. Dynamical structure of truncated M- fractional Klein-Gordon model via two integral schemes. Results Phys. 2023;46:106272. 10.1016/j.rinp.2023.106272Search in Google Scholar

[41] Ozdemir N. M-truncated soliton solutions of the fractional (4+1)-dimensional Fokas equation. Int J Optimiz Control Theories Appl (IJOCTA). 2023;13(1):123–9. 10.11121/ijocta.2023.1321Search in Google Scholar

[42] Malingam P, Wongsasinchai P, Sirisubtawee S, Koonprasert S. Exact solutions of the paraxial wave dynamical model in Kerr media with truncated M-fractional derivative using the (G′∕G,1∕G)-expansion method. WSEAS Trans Syst Control. 2023;18:498–512. 10.37394/23203.2023.18.53Search in Google Scholar

[43] Khalil R, Al Horani M, Yousef A, Sababheh M. A new definition of fractional derivative. J Comput Appl Math. 2014;264:65–70. 10.1016/j.cam.2014.01.002Search in Google Scholar

[44] Katugampola UN. A new fractional derivative with classical properties. Classical Analysis and ODEs. 2014;1–8. https://arxiv.org/abs/1410.6535. Search in Google Scholar

[45] Demiray S, Ünsal Ö, Bekir A. Exact solutions of nonlinear wave equations using (G′G,1G)-expansion method. J Egypt Math Soc. 2015;23(1):78–84. 10.1016/j.joems.2014.02.011Search in Google Scholar

[46] Miah MM, Ali HS, Akbar MA, Seadawy AR. New applications of the two variable (G′∕G,1∕G)-expansion method for closed form traveling wave solutions of integro-differential equations. J Ocean Eng Sci. 2019;4(2):132–43. 10.1016/j.joes.2019.03.001Search in Google Scholar

[47] Alurrfi KA, Shahoot AM, Elhasadi OI. Exact solutions for the GKdV-mKdV equation with higher-order nonlinear terms using the generalized (G′∕G,1∕G)-expansion method and the generalized Liénardequation. Ricerche di Matematica. 2024;73:887–905. 10.1007/s11587-021-00637-6Search in Google Scholar

[48] Sadaf M, Arshed S, Akram G, Iqra. Exact soliton and solitary wave solutions to the Fokas system using two variables (G′∕G,1∕G)-expansion technique and generalized projective Riccati equation method. Optik. 2022;268:169713. 10.1016/j.ijleo.2022.169713Search in Google Scholar

[49] Ahmad S, Mahmoud EE, Saifullah S, Ullah A, Ahmad S, Akgül A, et al. New waves solutions of a nonlinear Landau-Ginzburg-Higgs equation: The Sardar-subequation and energy balance approaches. Results Phys. 2023;51:106736. 10.1016/j.rinp.2023.106736Search in Google Scholar

[50] Hussain R, Imtiaz A, Rasool T, Rezazadeh H, Inç M. Novel exact and solitary solutions of conformable Klein-Gordon equation via Sardar-subequation method. J Ocean Eng Sci. 2022;7(1):1–7. 10.1016/j.joes.2022.04.036Search in Google Scholar

[51] Asjad MI, Inc M, Iqbal I. Exact solutions for new coupled Konno-Oono equation via Sardar subequation method. Optical Quantum Electron. 2022;54(12):798. 10.1007/s11082-022-04208-3Search in Google Scholar

[52] Rezazadeh H, Inc M, Baleanu D. New solitary wave solutions for variants of (3+1)-dimensional Wazwaz-Benjamin-Bona-Mahony equations. Front Phys. 2020;8:332. 10.3389/fphy.2020.00332Search in Google Scholar

[53] Tang L. Bifurcation analysis and optical soliton solutions for the fractional complex Ginzburg-Landau equation in communication systems. Optik. 2023;276:170639. 10.1016/j.ijleo.2023.170639Search in Google Scholar

[54] Nguetcho AT, Li J, Bilbault JM. Bifurcations of phase portraits of a singular nonlinear equation of the second class. Commun Nonl Sci Numer Simul. 2014;19(8):2590–601. 10.1016/j.cnsns.2013.12.022Search in Google Scholar

[55] Ali KK, Mehanna M. Computational and analytical solutions to modified Zakharov-Kuznetsov model with stability analysis via efficient techniques. Opt Quantum Electron. 2021;53(12):723. 10.1007/s11082-021-03363-3Search in Google Scholar

[56] Zayed E, Alurrfi K. The (G′∕G,1∕G)-expansion method and its applications to two nonlinear Schrödinger equations describing the propagation of femtosecond pulses in nonlinear optical fibers. Optik. 2016;127(4):1581–89. 10.1016/j.ijleo.2015.11.027Search in Google Scholar

Received: 2024-07-04

Revised: 2024-11-16

Accepted: 2025-01-10

Published Online: 2025-02-11

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/phys-2024-0117

Keywords for this article

truncated M-fractional derivative; space–time fractional mZK equation; (G′/G, 1/G)-expansion method; Sardar subequation method; bifurcation analysis

Creative Commons

BY 4.0