Applications of the (G′/G²)-Expansion Method for Solving Certain Nonlinear Conformable Evolution Equations

Abstract

The core objective of this article is to generate novel exact traveling wave solutions of two nonlinear conformable evolution equations, namely, the $(2+1)$-dimensional conformable time integro-differential Sawada–Kotera (SK) equation and the $(3+1)$-dimensional conformable time modified KdV–Zakharov–Kuznetsov (mKdV–ZK) equation using the $(G^{\prime }/G^2)$-expansion method. These two equations associate with conformable partial derivatives with respect to time which the former equation is firstly proposed in the form of the conformable integro-differential equation. To the best of the authors’ knowledge, the two equations have not been solved by means of the $(G^{\prime }/G^2)$-expansion method for their exact solutions. As a result, some exact solutions of the equations expressed in terms of trigonometric, exponential, and rational function solutions are reported here for the first time. Furthermore, graphical representations of some selected solutions, plotted using some specific sets of the parameter values and the fractional orders, reveal certain physical features such as a singular single-soliton solution and a doubly periodic wave solution. These kinds of the solutions are usually discovered in natural phenomena. In particular, the soliton solution, which is a solitary wave whose amplitude, velocity, and shape are conserved after a collision with another soliton for a nondissipative system, arises ubiquitously in fluid mechanics, fiber optics, atomic physics, water waves, and plasmas. The method, along with the help of symbolic software packages, can be efficiently and simply used to solve the proposed problems for trustworthy and accurate exact solutions. Consequently, the method could be employed to determine some new exact solutions for other nonlinear conformable evolution equations.

1. Introduction

Nonlinear evolution equations (NLEEs), which are interpreted as the differential law of the development in time of a system and typically expressed in terms of nonlinear partial differential equations (NPDEs), can be utilized to describe many interesting and sophisticated phenomena in physics, mathematical physics, engineering, and other various scientific fields such as fluid mechanics [1,2], plasma physics [3], quantum mechanics [4], biology [5], nonlinear wave theory [6], and fiber optics [7]. Some applications of NLEEs for natural events are as follows [8,9,10]: The nonlinear Schrödinger’s equation explains the dynamics of propagation for solitons through optical fibers. The Korteweg de Vries (KdV) equation can be used to model the shallow water wave dynamics near ocean shore and beaches. The dynamics of an incompressible viscoelastic Kelvin–Voigt fluid can be described by the Oskolkov equation. In addition, an epidemic model on a network such as the present epidemic of COVID-19 comprising susceptible–infected–recovered equations at the nodes, coupled by diffusion using a graph Laplacian, can be analyzed using a system of NLEEs.

To further elucidate some important behaviors of the phenomena modeled by NLEEs, extracting their exact analytic solutions, in particular, solitary wave solutions, is of great significance. The exploration of exact solutions for NLEEs has become quite prominent due to the recent, considerable advances in computational methods and symbolic software packages. Numerous kinds of exact solutions such as solitons, positons, complexitons, dromions, cuspon, rational, kink, periodic, and quasiperiodic solutions have been obtained via solving integrable NLEEs. In the past couple of decades, many robust, efficient, and powerful methods exist that have been developed for finding exact solutions of NLEEs, including the $(G{}^{\prime }/G,1/G)$-expansion method [11], the enhanced $(G^{\prime }/G)$-expansion method [12], the exp-function method [13], the Jacobi elliptic equation method [14], the generalized Kudryashov’s method [15], the sine-Gordon expansion method [16], the sub-equation method [17], the improved $tan(\varphi /2)$-expansion method [18], and the extended direct algebraic method [19,20]. More recently, the $(G^{\prime }/G^2)$-expansion method [4,21,22,23,24,25,26,27,28,29] has attracted a remarkable amount of attention of many researchers who employed the method to construct exact solutions of certain NPDEs. In 2018, Arshed and Sadia [23] used the $(G^{\prime }/G^2)$-expansion method to obtain some new traveling wave solutions for the time-fractional Burgers equation, the fractional biological population model, and the space-time fractional Whitham–Broer–Kaup equations. Sirisubtawee and Koonprasert [24] utilized the method to solve the Benny–Luke equation, the equation of nanoionic currents along microtubules, and the generalized Hirota–Satsuma coupled KdV system for their exact solutions including trigonometric, exponential, and rational function solutions. In 2020, the $(G^{\prime }/G^2)$-expansion approach was employed to construct some novel exact traveling wave solutions of the $(2+1)$-dimensional Boiti–Leon–Pempinelli system [28]. In 2021, Bilal et al. [29] proposed new exact solutions, which consist of shock, singular, shock-singular, and singular periodic wave solutions obtained by the $(G^{\prime }/G^2)$-expansion method, to unidirectional Dullin–Gottwald–Holm (DGH) system describing the prorogation of waves in shallow water.

In this article, we will demonstrate the use of the $\left(G^{\prime }/G^2\right)$-expansion method to construct explicit exact solutions for the following two interesting problems in mathematical physics:

1. The $(2+1)$-dimensional conformable time integro-differential Sawada–Kotera (SK) equation can be expressed as

$$\begin{array}{c}\hfill {\partial }_t^{\alpha }u={\left(u_{xxxx}+5u{u}_{xx}+\frac{5}{3}{u}^3+5u_{xy}\right)}_x-5{\partial }_x^{-1}\left(u_{yy}\right)+5u{u}_y+5u_x{\partial }_x^{-1}\left(u_y\right),\end{array}$$

(1)

where ${\partial }_t^{\alpha }\left(·\right)=\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}\left(·\right)$ is the conformable partial derivative with respect to t of order $\alpha$ with $0<\alpha \le 1$ and ${\partial }_x^{-1}\left(·\right)={\int }_{-\infty }^x\left(·\right)dx$. The dependent variable u in the equation is a multivariable function consisting of three independent variables $x,y$, and $t,$ i.e., $u=u(x,y,t)$. If $\alpha =1$, then Equation (1) reduces into the $(2+1)$-dimensional integro-differential SK equation [30,31,32,33,34,35,36], which was initially established by Konopelcheno and Dubrovsky [37], using the inverse scattering transform method. The $(2+1)$-dimensional SK equation has been investigated extensively and intensively in a number of studies in the literature because of its significance and various applications in two-dimensional quantum gravity field theory, conformal field theory, and conserved current of Liouville equation [38,39,40].

2. The $(3+1)$-dimensional conformable time modified KdV–Zakharov–Kuznetsov (mKdV–ZK) equation reads

$${\partial }_t^{\alpha }u+{\delta }_1{u}^2{u}_x+{\delta }_2{u}_{xxx}+{\delta }_3{(u_{yy}+u_{zz})}_x=0,$$

(2)

where ${\partial }_t^{\alpha }\left(·\right)=\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}\left(·\right)$ represents the conformable partial derivative with respect to t of order $\alpha$ with $0<\alpha \le 1$, and where u is a function of independent variables $x,y,z$, and t. The parameters ${\delta }_1,{\delta }_2,{\delta }_3$ in the equation are real constants. If $\alpha =1$ is inserted into Equation (2), then the equation becomes the (3 + 1)-dimensional mKdV–ZK equation [41,42,43]. The mKdV–ZK equation plays a significant role in explaining dynamics of many branches of physics such as plasma physics, nonlinear optics, fluid dynamics, shallow water waves in oceanography, quantum mechanics and mathematical physics so that fundamental properties of nonlinear propagation for such various physical phenomena are analyzed [41,43,44,45]. Particularly, the (3 + 1)-dimensional mKdV–ZK equation is utilized to control the behavior of weakly nonlinear ion-acoustic waves in magnetized electron–positron plasma including the same hot and cold components of each species [42].

A recent literature review for constructing explicit exact solutions of the integro-differential SK equation using various methods such as the Hirota bilinear method, the $(G{}^{\prime }/G,1/G)$-expansion method, and the generalized Kudryashov method (GKM) can be found in [32,33,34,35,36,46,47]. Furthermore, some scientists have devoted substantial efforts to finding exact solutions of the (3 + 1)-dimensional mKdV–ZK equation in the sense of the classical partial, conformable, and Jumarie’s modified Riemann–Liouville derivatives using different and reliable approaches. In the past few years, the investigation of exact traveling wave solutions for such mKdV–ZK equations has been discussed in [41,42,43,48,49,50,51]. To the best of the authors’ knowledge, there are no research scholars who have found explicit exact solutions for Equations (1) and (2) using the $\left(G^{\prime }/G^2\right)$-expansion method.

The organization of this paper is as follows: We provide a brief description of the conformable derivative and its crucial characteristics in Section 2. Section 3 is devoted to compactly describing the key steps of the $\left(G^{\prime }/G^2\right)$-expansion method. In Section 4, the extraction of exact solutions of Equations (1) and (2) using the proposed technique is illustrated. Some graphical representations of the chosen solutions and their physical explanations are presented in Section 5. The last section summarizes the results of the current study.

2. Conformable Derivative and Its Properties

In this section, a definition of the conformable derivative, which was initially introduced by Khalil et al. [52], and its essential characteristics are briefly presented. They will be utilized for the remaining parts of the present article.

Definition 1.

Let f be a function such that $f:[0,\infty )\to \mathbb{R}$. Then, the conformable derivative of f of order $\alpha ,$ where $0<\alpha \le 1,$ is defined as [52,53,54,55,56,57,58]

$$\begin{array}{c}\hfill D_t^{\alpha }f\left(t\right)=\underset{\epsilon \to 0}{lim}{\displaystyle \frac{f(t+\epsilon t^{1-\alpha })-f\left(t\right)}{\epsilon }},\end{array}$$

(3)

for all $t>0$. If the limit in Equation (3) exists, then we can state that f is α-conformable differentiable at a point $t>0$. In addition, if f is α-conformable differentiable in some $(0,a),a>0$ and ${lim}_{t\to 0^{+}}{D}_t^{\alpha }f\left(t\right)$ exists, then we define $D_t^{\alpha }f\left(0\right)={lim}_{t\to 0^{+}}{D}_t^{\alpha }f\left(t\right).$

Let $\alpha \in (0,1]$, and $f\left(t\right),g\left(t\right)$ be $\alpha$-conformable differentiable functions at a point $t>0$. Then, the important properties of the conformable derivative are as follows [52,53,55,56,57,59]:

(1)

$D_t^{\alpha }\left(\lambda \right)=0,where\lambda =constant.$

(2)

$D_t^{\alpha }\left(t^{\mu }\right)=\mu t^{\mu -\alpha },$ for all $\mu \in \mathbb{R}.$

(3)

$D_t^{\alpha }(af\left(t\right)+bg\left(t\right))=a{D}_t^{\alpha }f\left(t\right)+b{D}_t^{\alpha }g\left(t\right),$ for all $a,b\in \mathbb{R}.$

(4)

$D_t^{\alpha }\left(f\left(t\right)g\left(t\right)\right)=f\left(t\right)D_t^{\alpha }g\left(t\right)+g\left(t\right)D_t^{\alpha }f\left(t\right).$

(5)

$D_t^{\alpha }\left({\displaystyle \frac{f\left(t\right)}{g\left(t\right)}}\right)={\displaystyle \frac{g\left(t\right)D_t^{\alpha }f\left(t\right)-f\left(t\right)D_t^{\alpha }g\left(t\right)}{g{\left(t\right)}^2}}.$

(6)

If, in addition, f is differentiable, then $D_t^{\alpha }\left(f\left(t\right)\right)=t^{1-\alpha }{\displaystyle \frac{df\left(t\right)}{dt}}$.

Remark 1.

Using the definition in (3) and the above properties, the conformable derivatives of certain interesting functions are defined as follows [52,53,55,56,57]:

(1) 

$D_t^{\alpha }\left(e^{at}\right)=a{t}^{1-\alpha }{e}^{at},a\in \mathbb{R}$.

(2) 

$D_t^{\alpha }(sinbt)=b{t}^{1-\alpha }cosbt,b\in \mathbb{R}$.

(3) 

$D_t^{\alpha }(cosbt)=-b{t}^{1-\alpha }sinbt,b\in \mathbb{R}$.

(4) 

$D_t^{\alpha }\left(\frac{1}{\alpha }{t}^{\alpha }\right)=1$.

Theorem 1.

[56,57,59,60,61] Suppose $f,g:(0,\infty )\to \mathbb{R}$ are differentiable and also α-conformable differentiable. Further, assume that g is a function defined in the range of f. Then, we have

$$D_t^{\alpha }(f\circ g)\left(t\right)=t^{1-\alpha }{f}^{\prime }\left(g\left(t\right)\right)g^{\prime }\left(t\right),$$

where the prime symbol ${(}^{\prime })$ denotes the classical derivative.

The definition of the conformable derivative and its relevant properties when a fractional order $\alpha \in (n,n+1]$ for some positive integer n are described as follows:

Definition 2.

Let $\alpha \in (n,n+1]$, where n is a positive integer. Further, assume that f is n-times differentiable at $t>0$. Then, the conformable derivative of f of order $\alpha >1$ can be defined as [52,62]

$$\begin{array}{c}\hfill D_t^{\alpha }f\left(t\right)=\underset{\epsilon \to 0}{lim}{\displaystyle \frac{f^{(\lceil \alpha \rceil -1)}(t+\epsilon t^{\lceil \alpha \rceil -\alpha })-f^{(\lceil \alpha \rceil -1)}\left(t\right)}{\epsilon }},\end{array}$$

(4)

where $\lceil \alpha \rceil$ is the smallest integer greater than or equal to α.

Remark 2.

Using the definition in (4) and assuming that f is $(n+1)$-times differentiable at $t>0$, we consequently have [52]

$$D_t^{\alpha }f\left(t\right)=t^{\lceil \alpha \rceil -\alpha }{f}^{\left(\lceil \alpha \rceil \right)}\left(t\right),$$

(5)

where $\alpha \in (n,n+1]$ for some positive integer n.

Remark 3.

Suppose that f is a twice differentiable function at $t>0$.

(1) 

If $\alpha \in (0,1]$, then $D_t^{\alpha }\left(D_t^{\alpha }f\left(t\right)\right)=D_t^{\alpha }\left(t^{1-\alpha }{f}^{\prime }\left(t\right)\right)=t^{2-2\alpha }{f}^{″}\left(t\right)+(1-\alpha )t^{1-2\alpha }{f}^{\prime }\left(t\right)$.

(2) 

If $\alpha \in (0,\frac{1}{2}]$, then $D_t^{2\alpha }f\left(t\right)=t^{1-2\alpha }{f}^{\prime }\left(t\right)$.

(3) 

If $\alpha \in (\frac{1}{2},1]$, then $D_t^{2\alpha }f\left(t\right)=t^{2-2\alpha }{f}^{″}\left(t\right)$.

(4) 

If $\alpha \in (0,1]$, then $D_t^{2\alpha }f\left(t\right)\ne D_t^{\alpha }\left(D_t^{\alpha }f\left(t\right)\right)$.

(5) 

For some positive integer n, further assume that f is $(n+1)$-times differentiable at $t>0$. In general, if $\alpha \in (0,1]$, then $\underset{ntimes}{\underbrace{D_t^{\alpha }\left(D_t^{\alpha }\left(...\left(D_t^{\alpha }f\left(t\right)\right)\right)\right)}}\ne D_t^{n\alpha }f\left(t\right)$.

Using the definition in (3), we can define, for example, the conformable partial derivative of a function $u=u(x,t)$ with respect to t of order $\alpha \in (0,1]$ as

$$\begin{array}{c}\hfill {\partial }_t^{\alpha }u(x,t)=\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}u(x,t)=\underset{\epsilon \to 0}{lim}\frac{u(x,t+\epsilon t^{1-\alpha })-u(x,t)}{\epsilon },t>0.\end{array}$$

(6)

Analogously, using the definition in (4) for $\alpha \in (n,n+1]$, where n is a positive integer, if $u=u(x,t)$ is assumed to be $(n+1)$-times partial differentiable with respect to t, then we obtain for $t>0$

$${\partial }_t^{\alpha }u(x,t)=t^{\lceil \alpha \rceil -\alpha }\frac{{\partial }^{\lceil \alpha \rceil }}{\partial t^{\lceil \alpha \rceil }}u(x,t).$$

(7)

Remark 4.

The reason we do not replace some higher-order classical partial derivatives in Equations (1) and (2) with their corresponding conformable partial derivatives, for instance, replacing the term $u_{xx}$ with ${\partial }_x^{2\alpha }u(x,t),\alpha \in (0,1]$, is because the conformable derivative does not have the sequential derivative property as specified in property (5) of Remark 3. Thus, a conversion of a conformable partial differential equation using a traveling wave transformation to an ordinary differential equation of a new variable may still have some independent variables of the original equation left. This is not what we desire to obtain in the process of seeking exact traveling wave solutions of the conformable NPDEs.

For example, suppose that $u=u(x,t)$ and we use the traveling wave transformation $\xi =\frac{c{x}^{\alpha }}{\alpha }+\frac{k{t}^{\alpha }}{\alpha },$ where $\alpha \in (\frac{1}{2},1]$ and $c,k$ are constants. Introducing a new dependent variable U such that $u(x,t)=U\left(\xi \right)$, if the sequential derivative property for the conformable partial derivative holds, then we would have ${\partial }_t^{2\alpha }u(x,t)={\partial }_t^{\alpha }\left({\partial }_t^{\alpha }u(x,t)\right)=k^2{U}^{″}\left(\xi \right)$. However, the actual conversion of the term ${\partial }_t^{2\alpha }u(x,t)$ is

$${\partial }_t^{2\alpha }u(x,t)=t^{2-2\alpha }{u}_{tt}(x,t)=k^2{U}^{″}\left(\xi \right)+(\alpha -1)k{t}^{-\alpha }{U}^{\prime }\left(\xi \right).$$

(8)

3. Algorithm of the $\left({\mathit{G}}^{\prime }/{\mathit{G}}^{\mathbf{2}}\right)$-Expansion Method

In this section, we briefly describe the $\left(G^{\prime }/G^2\right)$-expansion method, which is discussed in [21,22,23,24,25,26,27,29,63]. Consider the nonlinear conformable partial differential equation of the unknown function $u=u(x_1,x_2,...,x_n,t)$ consisting of the independent variables $x_1,x_2,...,x_n$, and t as follows:

$$\begin{array}{c}\hfill P\left(u,{\partial }_t^{\alpha }u,{\partial }_{x_1}^{{\beta }_1}u,...,{\partial }_{x_n}^{{\beta }_n}u,u_{tt},u_{x_1{x}_1},...,u_{x_n{x}_n},{\partial }_t^{\alpha }\left({\partial }_{x_1}^{{\beta }_1}u\right),...\right)=0,0<\alpha ,{\beta }_1,{\beta }_2,...,{\beta }_n\le 1,\end{array}$$

(9)

where ${\partial }_v^{\gamma }u=\frac{{\partial }^{\gamma }}{\partial v^{\gamma }}u$ is a generic term for the conformable partial derivative of the dependent variable u with respect to the independent variable v of order $\gamma \in (0,1]$, and where the subscript symbols denote the classical partial derivatives, for instance, $u_{tt}=\frac{{\partial }^2}{\partial t^2}u$. The function P in (9) is a polynomial of u and its various partial derivatives. The main steps of the $\left(G^{\prime }/G^2\right)$-expansion method for constructing exact solutions for Equation (9) can be given as follows:

Step 1: Convert nonlinear conformable partial differential Equation (9) into an ordinary differential equation (ODE) via the fractional complex traveling wave transformation in a variable $\xi$,

$$\begin{array}{c}\hfill u(x_1,x_2,...,x_n,t)=U\left(\xi \right),\xi =\frac{c_1{x}_1^{{\beta }_1}}{{\beta }_1}+\frac{c_2{x}_2^{{\beta }_2}}{{\beta }_2}+...+\frac{c_n{x}_n^{{\beta }_n}}{{\beta }_n}+\frac{k{t}^{\alpha }}{\alpha },\end{array}$$

(10)

where $c_1,c_2,...,c_n,k$ are nonzero constants that will be determined at a later step. Applying transformation (10) to (9) and then integrating the resulting equation with respect to $\xi$ as many as possible, we obtain the following ODE in $U=U\left(\xi \right)$:

$$\begin{array}{c}\hfill Q(U,U^{\prime },U^{″},U^{‴},...)=0,\end{array}$$

(11)

where Q is a polynomial function of $U\left(\xi \right)$, and its various integer-order derivatives. The prime notation (${}^{\prime }$) denotes the ordinary derivative with respect to $\xi$.

Step 2: Suppose that the general solution of the above ODE can be expressed in terms of $\left(G^{\prime }/G^2\right)$ as

$$\begin{array}{c}\hfill U\left(\xi \right)=a_0+\sum _{j=1}^N\left[a_j{\left({\displaystyle \frac{G^{\prime }}{G^2}}\right)}^j+b_j{\left({\displaystyle \frac{G^{\prime }}{G^2}}\right)}^{-j}\right],\end{array}$$

(12)

where $G=G\left(\xi \right)$ satisfies the simple Riccati equation:

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{G^{\prime }}{G^2}}\right)}^{\prime }=\mu +\lambda {\left({\displaystyle \frac{G^{\prime }}{G^2}}\right)}^2,\end{array}$$

(13)

in which $\mu \ne 1$ and $\lambda \ne 0$ are arbitrary integers. The unknown constants $a_N$ or $b_N$ may be zero, but both of them cannot be zero simultaneously. The coefficients $a_0,a_j,b_j$$(j=1,2,...,N)$ are unknown constants to be determined in Step 3. The value of the positive integer N can be computed using the homogeneous balance principle, in other words, by balancing between the highest order derivatives and the nonlinear terms appearing in Equation (11). More precisely, if the degree of $U\left(\xi \right)$ is $\mathrm{Deg}\left[U\right(\xi \left)\right]=N$, then the degree of the following terms can be calculated using the following formulas [57]:

$$\begin{array}{c}\hfill \mathrm{Deg}\left[{\displaystyle \frac{d^{q}U\left(\xi \right)}{d{\xi }^q}}\right]=N+q,\mathrm{Deg}\left[{\left(U\left(\xi \right)\right)}^p{\left({\displaystyle \frac{d^{q}U\left(\xi \right)}{d{\xi }^q}}\right)}^s\right]=Np+s(N+q).\end{array}$$

(14)

Step 3: Substituting Equation (12), along with Equation (13), into Equation (11), we have a polynomial in $\left(G^{\prime }/G^2\right)$. Collecting all coefficients of the same power of ${\left(G^{\prime }/G^2\right)}^i$ where $i=0,\pm 1,\pm 2,...,\pm M$ where M is some positive integer and then setting all of the obtained coefficients to zero, we obtain a system of nonlinear algebraic equations for the unknown constants $a_0,a_j,b_j(j=1,2,...,N),c_1,c_2,...,c_n$ and k. Assume that the resulting algebraic system can be possibly solved for the unknown constants using symbolic software packages such as Maple.

Step 4: The general solutions of Equation (13) can be separated into the following three cases depending on the values of $\mu$ and $\lambda$:

If $\mu \lambda >0$, then (13) has the trigonometric function solution as

$$\begin{array}{c}\hfill {\displaystyle \frac{G^{\prime }}{G^2}}=\sqrt{{\displaystyle \frac{\mu }{\lambda }}}\left({\displaystyle \frac{Ccos\left(\sqrt{\mu \lambda }\xi \right)+Dsin\left(\sqrt{\mu \lambda }\xi \right)}{Dcos\left(\sqrt{\mu \lambda }\xi \right)-Csin\left(\sqrt{\mu \lambda }\xi \right)}}\right),\end{array}$$

(15)

where $C,D$ are arbitrary nonzero constants.

If $\mu \lambda <0$, then (13) has the exponential function solution as

$$\begin{array}{c}\hfill {\displaystyle \frac{G^{\prime }}{G^2}}=\frac{1}{2\lambda }\left(2\sqrt{\left|\mu \lambda \right|}-{\displaystyle \frac{4C\sqrt{\left|\mu \lambda \right|}{e}^{2\xi \sqrt{\left|\mu \lambda \right|}}}{C{e}^{2\xi \sqrt{\left|\mu \lambda \right|}}-D}}\right),\end{array}$$

(16)

which is equivalent to the hyperbolic function solution as

$$\begin{array}{c}\hfill {\displaystyle \frac{G^{\prime }}{G^2}}=-{\displaystyle \frac{\sqrt{\left|\mu \lambda \right|}}{\lambda }}\left({\displaystyle \frac{Csinh\left(2\sqrt{\left|\mu \lambda \right|}\xi \right)+Ccosh\left(2\sqrt{\left|\mu \lambda \right|}\xi \right)+D}{Csinh\left(2\sqrt{\left|\mu \lambda \right|}\xi \right)+Ccosh\left(2\sqrt{\left|\mu \lambda \right|}\xi \right)-D}}\right),\end{array}$$

where $C,D$ are arbitrary nonzero constants.

If $\mu =0$ and $\lambda \ne 0$, then (13) has the rational function solution as

$$\begin{array}{c}\hfill {\displaystyle \frac{G^{\prime }}{G^2}}=-{\displaystyle \frac{C}{\lambda \left(C\xi +D\right)}},\end{array}$$

(17)

where $C,D$ are arbitrary nonzero constants.

The exact traveling wave solutions of Equation (9) can be obtained by inserting the obtained values of $a_0,a_j,b_j(j=1,2,...,N),c_1,c_2,...,c_n,k$ and the solutions (15)–(17) into Equation (12) with the transformation (10).

4. Applications of the $\left({\mathit{G}}^{\prime }/{\mathit{G}}^{\mathbf{2}}\right)$-Expansion Method

In this section, we will implement the $\left(G^{\prime }/G^2\right)$-expansion method to solve Equations (1) and (2) for their exact solutions.

4.1. The $(2+1)$-Dimensional Conformable Time Integro-Differential Sawada–Kotera Equation

In this subsection, the $\left(G^{\prime }/G^2\right)$-expansion method will be utilized to extract exact traveling wave solutions of the $(2+1)$-dimensional conformable time integro-differential SK Equation (1). Using the transformation $u(x,y,t)=v_x(x,y,t)$ to convert (1) to a new nonlinear PDE, we obtain

$$\begin{array}{c}\hfill {\partial }_t^{\alpha }\left(v_x\right)=v_{xxxxxx}+5{\left(v_x{v}_{xxx}\right)}_x+\frac{5}{3}{\left(v_x^3\right)}_x+5v_{xxxy}-5v_{yy}+5v_x{v}_{xy}+5v_{xx}{v}_y.\end{array}$$

(18)

Applying the traveling wave transformation

$$\begin{array}{c}\hfill v(x,y,t)=V\left(\xi \right),\xi =x+y-\frac{k{t}^{\alpha }}{\alpha },\end{array}$$

(19)

where k is a constant to the resulting PDE with the use of Theorem 1, Equation (18) is transformed to the following ODE:

$$\begin{array}{c}\hfill -k{V}^{″}=V^{\left(6\right)}+5{\left(V^{\prime }{V}^{‴}\right)}^{\prime }+\frac{5}{3}{\left({\left(V^{\prime }\right)}^3\right)}^{\prime }+5V^{\left(4\right)}-5V^{″}+10V^{\prime }{V}^{″},\end{array}$$

(20)

where the prime notation ${(}^{\prime })$ denotes the ordinary derivative with respect to $\xi$. Integrating (20) with respect to $\xi$, we obtain

$$\begin{array}{c}\hfill V^{\left(5\right)}+5V^{\prime }{V}^{‴}+\frac{5}{3}{\left(V^{\prime }\right)}^3+5V^{‴}+(k-5)V^{\prime }+5{\left(V^{\prime }\right)}^2+c=0,\end{array}$$

(21)

where c is a constant of integration. Taking the transformation $W=V^{\prime }$, we obtain

$$\begin{array}{c}\hfill W^{\left(4\right)}+5W{W}^{″}+\frac{5}{3}{W}^3+5W^{″}+(k-5)W+5W^2+c=0.\end{array}$$

(22)

On the basis of Equation (12), we assume that the general solution of (22) takes the form

$$\begin{array}{c}\hfill W\left(\xi \right)=a_0+\sum _{j=1}^N\left[a_j{\left({\displaystyle \frac{G^{\prime }}{G^2}}\right)}^j+b_j{\left({\displaystyle \frac{G^{\prime }}{G^2}}\right)}^{-j}\right],\end{array}$$

(23)

for which $\mathrm{Deg}\left[W\right(\xi \left)\right]=N$, and the function G satisfies (13). Balancing the highest order derivative $W^{\left(4\right)}$ in (22) with the nonlinear term $W^3$ via the formulas (14), we obtain the balancing number $N=2$. In consequence, the solution form of ODE (22) can be written as

$$\begin{array}{c}\hfill W\left(\xi \right)=a_0+a_1\left(\frac{G^{\prime }}{G^2}\right)+b_1{\left(\frac{G^{\prime }}{G^2}\right)}^{-1}+a_2{\left(\frac{G^{\prime }}{G^2}\right)}^2+b_2{\left(\frac{G^{\prime }}{G^2}\right)}^{-2},\end{array}$$

(24)

where $a_0,a_1,a_2,b_1,$ and $b_2$ are unknown constants that will be found at a later step. Substituting Equation (24) into Equation (22), along with Equation (13), and then collecting all the coefficients of the same power of ${\left(G^{\prime }/G^2\right)}^i,(i=0,\pm 1,\pm 2,...)$, and finally setting these resulting coefficients to zero, we obtain the following system of algebraic equations in $a_0,a_1,a_2,b_1,b_2,k$:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\left(\frac{G^{\prime }}{G^2}\right)}^{-6}:& \frac{5b_2^3}{3}+30b_2^2{\mu }^2+120b_2{\mu }^4=0,\hfill \\ \hfill {\left(\frac{G^{\prime }}{G^2}\right)}^{-5}:& 24{\mu }^4{b}_1+40{\mu }^2{b}_1{b}_2+5b_1{b}_2^2=0,\hfill \\ \hfill {\left(\frac{G^{\prime }}{G^2}\right)}^{-4}:& 240\lambda {\mu }^3{b}_2+40\lambda \mu b_2^2+30{\mu }^2{a}_0{b}_2+10{\mu }^2{b}_1^2+30{\mu }^2{b}_2+5a_0{b}_2^2+5b_1^2{b}_2+5b_2^2=0,\hfill \\ \hfill {\left(\frac{G^{\prime }}{G^2}\right)}^{-3}:& \frac{5b_1^3}{3}+10b_1{\mu }^2+10b_1{b}_2+5a_1{b}_2^2+50b_1{b}_2\mu \lambda +10a_0{b}_1{\mu }^2+30a_1{b}_2{\mu }^2+10a_0{b}_1{b}_2\hfill \\ & +40b_1{\mu }^3\lambda =0,\hfill \end{array}\end{array}$$

(25)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\left(\frac{G^{\prime }}{G^2}\right)}^{-2}:& 136{\lambda }^2{\mu }^2{b}_2+10{\lambda }^2{b}_2^2+40\lambda \mu a_0{b}_2+10\lambda \mu b_1^2+10{\mu }^2{a}_1{b}_1+40{\mu }^2{a}_2{b}_2+40\lambda \mu b_2\hfill \\ & +5a_0^2{b}_2+5a_0{b}_1^2+10a_1{b}_1{b}_2+5a_2{b}_2^2+k{b}_2+10a_0{b}_2+5b_1^2-5b_2=0,\hfill \\ \hfill {\left(\frac{G^{\prime }}{G^2}\right)}^{-1}:& 16{\lambda }^2{\mu }^2{b}_1+10{\lambda }^2{b}_1{b}_2+10\lambda \mu a_0{b}_1+50\lambda \mu a_1{b}_2+20{\mu }^2{a}_2{b}_1+10\lambda \mu b_1+5a_0^2{b}_1\hfill \\ & +10a_0{a}_1{b}_2+5a_1{b}_1^2+10a_2{b}_1{b}_2+k{b}_1+10a_0{b}_1+10a_1{b}_2-5b_1=0,\hfill \\ \hfill {\left(\frac{G^{\prime }}{G^2}\right)}^0:& 10b_2{\lambda }^2+16a_2{\mu }^3\lambda +10a_1{b}_1+10a_2{b}_2+5a_1^2{b}_2+5b_1^2{a}_2+\frac{5a_0^3}{3}+c+5a_0^2+k{a}_0\hfill \\ & +10a_0{a}_2{b}_2+10a_0{a}_1{b}_1+10a_0{a}_2{\mu }^2+10a_0{b}_2{\lambda }^2+16b_2\mu {\lambda }^3+10a_2{\mu }^2-5a_0\hfill \\ & +80a_2{b}_2\mu \lambda +20a_1{b}_1\mu \lambda =0,\hfill \\ \hfill \left(\frac{G^{\prime }}{G^2}\right):& 16{\lambda }^2{\mu }^2{a}_1+20{\lambda }^2{a}_1{b}_2+10\lambda \mu a_0{a}_1+50\lambda \mu a_2{b}_1+10{\mu }^2{a}_1{a}_2+10\lambda \mu a_1+5a_0^2{a}_1\hfill \\ & +10a_0{a}_2{b}_1+5a_1^2{b}_1+10a_1{a}_2{b}_2+k{a}_1+10a_0{a}_1+10a_2{b}_1-5a_1=0,\hfill \\ \hfill {\left(\frac{G^{\prime }}{G^2}\right)}^2:& 136{\lambda }^2{\mu }^2{a}_2+10{\lambda }^2{a}_1{b}_1+40{\lambda }^2{a}_2{b}_2+40\lambda \mu a_0{a}_2+10\lambda \mu a_1^2+10{\mu }^2{a}_2^2+40\lambda \mu a_2\hfill \\ & +5a_0^2{a}_2+5a_0{a}_1^2+10a_1{a}_2{b}_1+5a_2^2{b}_2+k{a}_2+10a_0{a}_2+5a_1^2-5a_2=0,\hfill \\ \hfill {\left(\frac{G^{\prime }}{G^2}\right)}^3:& \frac{5a_1^3}{3}+5b_1{a}_2^2+10a_1{\lambda }^2+10a_1{a}_2+50a_1{a}_2\mu \lambda +10a_0{a}_1{\lambda }^2+30b_1{a}_2{\lambda }^2+10a_0{a}_1{a}_2\hfill \\ & +40a_1{\lambda }^3\mu =0,\hfill \\ \hfill {\left(\frac{G^{\prime }}{G^2}\right)}^4:& 240{\lambda }^3\mu a_2+30{\lambda }^2{a}_0{a}_2+10{\lambda }^2{a}_1^2+40\lambda \mu a_2^2+30{\lambda }^2{a}_2+5a_0{a}_2^2+5a_1^2{a}_2+5a_2^2=0,\hfill \\ \hfill {\left(\frac{G^{\prime }}{G^2}\right)}^5:& 24{\lambda }^4{a}_1+40{\lambda }^2{a}_1{a}_2+5a_1{a}_2^2=0,\hfill \\ \hfill {\left(\frac{G^{\prime }}{G^2}\right)}^6:& \frac{5a_2^3}{3}+30a_2^2{\lambda }^2+120a_2{\lambda }^4=0.\hfill \end{array}\end{array}$$

Solving the algebraic system in (25) with the assistance of the Maple package program, we have the following three results:

Result 1:

$$\begin{array}{c}\hfill \begin{array}{cc}& a_0=-\frac{1}{2}-4\mu \lambda +\frac{{\left(\omega \left(c,\mu ,\lambda \right)\right)}^{1/3}}{10}+\frac{\left(24{\lambda }^2{\mu }^2+\frac{5}{2}\right)}{{\left(\omega \left(c,\mu ,\lambda \right)\right)}^{1/3}},a_1=0,a_2=0,b_1=0,b_2=-6{\mu }^2,\hfill \\ & k=\frac{25}{4}-20{\lambda }^2{\mu }^2-\frac{{\left(\omega (c,\mu ,\lambda )\right)}^{1/3}}{2}-\frac{{\left(\omega (c,\mu ,\lambda )\right)}^{2/3}}{20}-\frac{\left(\frac{25}{2}+120{\lambda }^2{\mu }^2\right)}{{\left(\omega (c,\mu ,\lambda )\right)}^{1/3}}\hfill \\ & -\frac{\left(2880{\lambda }^4{\mu }^4+600{\lambda }^2{\mu }^2+\frac{125}{4}\right)}{{\left(\omega (c,\mu ,\lambda )\right)}^{2/3}},\hfill \end{array}\end{array}$$

(26)

where

$$\begin{array}{cc}\omega \left(c,\mu ,\lambda \right)\hfill & =3200{\left(\mu \lambda \right)}^3+1200{\left(\mu \lambda \right)}^2+150c-125\hfill \\ & +10\left(-35,840{\left(\mu \lambda \right)}^6+76,800{\left(\mu \lambda \right)}^5-28,800{\left(\mu \lambda \right)}^4+\left(9600c-8000\right){\left(\mu \lambda \right)}^3\right.\hfill \\ & {\left.+\left(3600c-7500\right){\left(\mu \lambda \right)}^2+225c^2-375c\right)}^{\frac{1}{2}},\hfill \end{array}$$

(27)

where c is an arbitrary constant and $\mu ,\lambda$ are arbitrary integers. Substituting (26) into solution form (24), along with the ratio $\left(G^{\prime }/G^2\right)$ in (15)–(17), depending on the values of $\mu$ and $\lambda$, we obtain the solution $W\left(\xi \right)$. Using the relation $v(x,y,t)=V\left(\xi \right),$ $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$, the transformation $u(x,y,t)=v_x(x,y,t)$ and $W=V^{\prime }$, we consequently have $u(x,y,t)=V^{\prime }\left(\xi \right)·{\xi }_x=V^{\prime }\left(\xi \right)=W\left(\xi \right),$ where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$.

If $\mu \lambda >0$, then the trigonometric function solution of (22) can be written as

$$\begin{array}{cc}\hfill W_1^1\left(\xi \right)=& -\frac{1}{2}-4\mu \lambda +\frac{{\left(\omega \left(c,\mu ,\lambda \right)\right)}^{1/3}}{10}+\frac{\left(24{\lambda }^2{\mu }^2+\frac{5}{2}\right)}{{\left(\omega \left(c,\mu ,\lambda \right)\right)}^{1/3}}\hfill \\ \hfill & -\frac{6\mu \lambda {\left(Dcos\left(\sqrt{\mu \lambda }\xi \right)-Csin\left(\sqrt{\mu \lambda }\xi \right)\right)}^2}{{\left(Ccos\left(\sqrt{\mu \lambda }\xi \right)+Dsin\left(\sqrt{\mu \lambda }\xi \right)\right)}^2}.\hfill \end{array}$$

(28)

Using the fact that $u(x,y,t)=W\left(\xi \right),$ where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$; thus, the exact solution of (1), written in terms of trigonometric functions, is

$$\begin{array}{c}\hfill u_1^1(x,y,t)=W_1^1\left(\xi \right),\end{array}$$

(29)

where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$ and k is defined in (26).

If $\mu \lambda <0$, then the exponential function solution of (22) can be expressed as

$$\begin{array}{c}\hfill W_2^1\left(\xi \right)=-\frac{1}{2}-4\mu \lambda +\frac{{\left(\omega \left(c,\mu ,\lambda \right)\right)}^{1/3}}{10}+\frac{\left(24{\lambda }^2{\mu }^2+\frac{5}{2}\right)}{{\left(\omega \left(c,\mu ,\lambda \right)\right)}^{1/3}}+-\frac{{6\lambda }^2{\mu }^2{\left(C{\mathrm{e}}^{2\xi \sqrt{\left|\mu \lambda \right|}}-D\right)}^2}{\left|\mu \lambda \right|{\left(C{\mathrm{e}}^{2\xi \sqrt{\left|\mu \lambda \right|}}+D\right)}^2}.\end{array}$$

(30)

Utilizing the fact that $u(x,y,t)=W\left(\xi \right),$ where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$; hence, the exact solution of (1), written in terms of exponential functions, is

$$\begin{array}{c}\hfill u_2^1(x,y,t)=W_2^1\left(\xi \right),\end{array}$$

(31)

where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$ and k is given in (26).

If $\mu =0,\lambda \ne 0$, then the rational function solution of (22) can be exhibited as

$$\begin{array}{c}\hfill W_3^1\left(\xi \right)=\frac{1}{2{\left(\tau \left(c\right)\right)}^{1/3}}\left({\frac{\left(\tau \left(c\right)\right)}{5}}^{2/3}-{\left(\tau \left(c\right)\right)}^{1/3}+5\right),\end{array}$$

(32)

where

$$\begin{array}{c}\hfill \tau \left(c\right)=\omega \left(c,0,\lambda \right)=-125+150c+50\sqrt{9c^2-15c}.\end{array}$$

(33)

Using the fact that $u(x,y,t)=W\left(\xi \right),$ where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$; thus, the exact solution of (1) for this case is

$$\begin{array}{c}\hfill u_3^1(x,y,t)=W_3^1\left(\xi \right),\end{array}$$

(34)

which is a constant, as shown in (32).

Result 2:

$$\begin{array}{c}\hfill \begin{array}{cc}& a_0=-\frac{1}{2}-4\mu \lambda +\frac{{\left(\omega \left(c,\mu ,\lambda \right)\right)}^{1/3}}{10}+\frac{\left(24{\lambda }^2{\mu }^2+\frac{5}{2}\right)}{{\left(\omega \left(c,\mu ,\lambda \right)\right)}^{1/3}},a_1=0,a_2=-6{\lambda }^2,b_1=0,b_2=0,\hfill \\ & k=\frac{25}{4}-20{\lambda }^2{\mu }^2-\frac{{\left(\omega \left(c,\mu ,\lambda \right)\right)}^{1/3}}{2}-\frac{{\left(\omega \left(c,\mu ,\lambda \right)\right)}^{2/3}}{20}-\frac{\left(\frac{25}{2}+120{\lambda }^2{\mu }^2\right)}{{\left(\omega \left(c,\mu ,\lambda \right)\right)}^{1/3}}\hfill \\ & -\frac{\left(2880{\lambda }^4{\mu }^4+600{\lambda }^2{\mu }^2+\frac{125}{4}\right)}{{\left(\omega \left(c,\mu ,\lambda \right)\right)}^{2/3}},\hfill \end{array}\end{array}$$

(35)

where $\omega \left(c,\mu ,\lambda \right)$ is defined in (27), c is an arbitrary constant and $\mu ,\lambda$ are arbitrary integers. Substituting (35) into solution (24), along with the ratio $\left(G^{\prime }/G^2\right)$, in (15)–(17), depending on the values of $\mu$ and $\lambda$, we obtain the solution $W\left(\xi \right)$.

If $\mu \lambda >0$, then the trigonometric function solution of (22) can be written as

$$\begin{array}{cc}\hfill W_1^2\left(\xi \right)=& -\frac{1}{2}-4\mu \lambda +\frac{{\left(\omega \left(c,\mu ,\lambda \right)\right)}^{1/3}}{10}+\frac{\left(24{\lambda }^2{\mu }^2+\frac{5}{2}\right)}{{\left(\omega \left(c,\mu ,\lambda \right)\right)}^{1/3}}\hfill \\ \hfill & -\frac{6\mu \lambda {\left(Ccos\left(\sqrt{\mu \lambda }\xi \right)+Dsin\left(\sqrt{\mu \lambda }\xi \right)\right)}^2}{{\left(Dcos\left(\sqrt{\mu \lambda }\xi \right)-Csin\left(\sqrt{\mu \lambda }\xi \right)\right)}^2}.\hfill \end{array}$$

(36)

Using the fact that $u(x,y,t)=W\left(\xi \right),$ where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$, so the exact solution of (1), written in terms of trigonometric functions, is

$$\begin{array}{c}\hfill u_1^2(x,y,t)=W_1^2\left(\xi \right),\end{array}$$

(37)

where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$ and k is defined in (35).

If $\mu \lambda <0$, then the exponential function solution of (22) can be expressed as

$$\begin{array}{c}\hfill W_2^2\left(\xi \right)=-\frac{1}{2}-4\mu \lambda +\frac{{\left(\omega \left(c,\mu ,\lambda \right)\right)}^{1/3}}{10}+\frac{\left(24{\lambda }^2{\mu }^2+\frac{5}{2}\right)}{{\left(\omega \left(c,\mu ,\lambda \right)\right)}^{1/3}}-\frac{6\left|\mu \lambda \right|{\left(C{\mathrm{e}}^{2\xi \sqrt{\left|\mu \lambda \right|}}+D\right)}^2}{{\left(C{\mathrm{e}}^{2\xi \sqrt{\left|\mu \lambda \right|}}-D\right)}^2}.\end{array}$$

(38)

Utilizing the fact that $u(x,y,t)=W\left(\xi \right),$ where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$, hence the exact solution of (1), written in terms of exponential functions, takes the form

$$\begin{array}{c}\hfill u_2^2(x,y,t)=W_2^2\left(\xi \right),\end{array}$$

(39)

where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$ and k is given in (35).

If $\mu =0,\lambda \ne 0$, then the rational function solution of (22) can be expressed as

$$\begin{array}{c}\hfill W_3^2\left(\xi \right)=\frac{\left(-5C^2{\xi }^2-10CD\xi -60C^2-5D^2\right){\left(\tau \left(c\right)\right)}^{1/3}+\left({\left(\tau \left(c\right)\right)}^{2/3}+25\right){\left(C\xi +D\right)}^2}{10{\left(\tau \left(c\right)\right)}^{1/3}{\left(C\xi +D\right)}^2},\end{array}$$

(40)

where $\tau \left(c\right)$ is shown in (33). Using the fact that $u(x,y,t)=W\left(\xi \right),$ where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$, thus the rational function solution of (1) is

$$\begin{array}{c}\hfill u_3^2(x,y,t)=W_3^2\left(\xi \right),\end{array}$$

(41)

where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$ in which k is reduced as

$$\begin{array}{c}\hfill k=\frac{-1}{20{\left(\tau \left(c\right)\right)}^{2/3}}\left({\left(\tau \left(c\right)\right)}^{4/3}+10\tau \left(c\right)-125{\left(\tau \left(c\right)\right)}^{2/3}+250{\left(\tau \left(c\right)\right)}^{1/3}+625\right).\end{array}$$

(42)

Result 3:

$$\begin{array}{c}\hfill \begin{array}{cc}& a_0=-\frac{1}{2}-4\mu \lambda +\frac{{\left(\theta \left(c,\mu ,\lambda \right)\right)}^{1/3}}{10}+\frac{\left(384{\lambda }^2{\mu }^2+\frac{5}{2}\right)}{{\left(\theta \left(c,\mu ,\lambda \right)\right)}^{1/3}},a_1=0,a_2=-6{\lambda }^2,b_1=0,b_2=-6{\mu }^2,\hfill \\ & k=\frac{25}{4}-320{\lambda }^2{\mu }^2-\frac{{\left(\theta \left(c,\mu ,\lambda \right)\right)}^{1/3}}{2}-\frac{{\left(\theta \left(c,\mu ,\lambda \right)\right)}^{2/3}}{20}-\frac{\left(1920{\lambda }^2{\mu }^2+\frac{25}{2}\right)}{{\left(\theta \left(c,\mu ,\lambda \right)\right)}^{1/3}}\hfill \\ & -\frac{\left(737,280{\lambda }^4{\mu }^4+9600{\lambda }^2{\mu }^2+\frac{125}{4}\right)}{{\left(\theta \left(c,\mu ,\lambda \right)\right)}^{2/3}},\hfill \end{array}\end{array}$$

(43)

where

$$\begin{array}{cc}\theta \left(c,\mu ,\lambda \right)=\hfill & 204,800{\left(\mu \lambda \right)}^3+19,200{\left(\mu \lambda \right)}^2+150c-125\hfill \\ & +10\left(-146,800,640{\left(\mu \lambda \right)}^6+78,643,200{\left(\mu \lambda \right)}^5-7,372,800{\left(\mu \lambda \right)}^4+614,400\left(c-\frac{5}{6}\right){\left(\mu \lambda \right)}^3\right.\hfill \\ & {\left.+57,600\left(c-\frac{25}{12}\right){\left(\mu \lambda \right)}^2+225c^2-375c\right)}^{\frac{1}{2}},\hfill \end{array}$$

(44)

where c is an arbitrary constant and $\mu ,\lambda$ are arbitrary integers. Substituting (43) into solution form (24), along with the ratio $\left(G^{\prime }/G^2\right)$, in (15)–(17), depending on the values of $\mu$ and $\lambda$, we obtain the solution $W\left(\xi \right)$.

If $\mu \lambda >0$, then the trigonometric function solution of (22) is written as

$$\begin{array}{cc}\hfill W_1^3\left(\xi \right)=& -\frac{1}{2}-4\mu \lambda +\frac{{\left(\theta \left(c,\mu ,\lambda \right)\right)}^{1/3}}{10}+\frac{\left(384{\lambda }^2{\mu }^2+\frac{5}{2}\right)}{{\left(\theta \left(c,\mu ,\lambda \right)\right)}^{1/3}}\hfill \\ \hfill & -\frac{6\mu \lambda {\left(Ccos\left(\sqrt{\mu \lambda }\xi \right)+Dsin\left(\sqrt{\mu \lambda }\xi \right)\right)}^2}{{\left(Dcos\left(\sqrt{\mu \lambda }\xi \right)-Csin\left(\sqrt{\mu \lambda }\xi \right)\right)}^2}.\hfill \end{array}$$

(45)

Using the fact that $u(x,y,t)=W\left(\xi \right),$ where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$; thus, the exact solution of (1), expressed in terms of trigonometric functions, is

$$\begin{array}{c}\hfill u_1^3(x,y,t)=W_1^3\left(\xi \right),\end{array}$$

(46)

where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$ and k is defined in (43).

If $\mu \lambda <0$, then the exponential function solution of (22) is

$$\begin{array}{c}\hfill W_2^3\left(\xi \right)=-\frac{1}{2}+2\mu \lambda +\frac{{\left(\theta \left(c,\mu ,\lambda \right)\right)}^{1/3}}{10}+\frac{\left(768{\lambda }^2{\mu }^2+5\right)}{2{\left(\theta \left(c,\mu ,\lambda \right)\right)}^{1/3}}-\frac{6\left|\mu \lambda \right|{\left(C{\mathrm{e}}^{2\xi \sqrt{\left|\mu \lambda \right|}}+D\right)}^2}{{\left(C{\mathrm{e}}^{2\xi \sqrt{\left|\mu \lambda \right|}}-D\right)}^2}.\end{array}$$

(47)

Utilizing the fact that $u(x,y,t)=W\left(\xi \right),$ where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$; hence, the exact solution of (1), written in terms of exponential functions, takes the form

$$\begin{array}{c}\hfill u_2^3(x,y,t)=W_2^3\left(\xi \right),\end{array}$$

(48)

where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$ and k is given in (43).

If $\mu =0,\lambda \ne 0$, then the rational function solution of (22) can be expressed as

$$\begin{array}{cc}\hfill W_3^3\left(\xi \right)=& \frac{1}{2{\left(\tau \left(c\right)\right)}^{1/3}{\left(C\xi +D\right)}^2}\hfill \\ \hfill & \times \left(\frac{{\left(C\xi +D\right)}^2{\left(\tau \left(c\right)\right)}^{2/3}}{5}-\left(\left({\xi }^2+12\right)C^2+2CD\xi +D^2\right){\left(\tau \left(c\right)\right)}^{1/3}+5{\left(C\xi +D\right)}^2\right).\hfill \end{array}$$

(49)

where $\tau \left(c\right)=\theta \left(c,0,\lambda \right)\left(=\omega \left(c,0,\lambda \right)\right)$ is expressed in (33). Utilizing the fact that $u(x,y,t)=W\left(\xi \right),$ where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$; thus, the rational function solution of (1) is

$$\begin{array}{c}\hfill u_3^3(x,y,t)=W_3^3\left(\xi \right),\end{array}$$

(50)

where $\xi =x+y-\frac{k{t}^{\alpha }}{\alpha }$ in which k is simplified, as shown in (42).

4.2. The (3 + 1)-Dimensional Conformable Time Modified KdV–Zakharov–Kuznetsov Equation

Before constructing exact traveling wave solutions of the $(3+1)$-dimensional space-time fractional modified KdV–Zakharov–Kuznetsov Equation (2) by means of the $\left(G^{\prime }/G^2\right)$-expansion method, we must convert the equation to an ordinary differential equation via the fractional complex transformation:

$$u(x,y,z,t)=U\left(\xi \right),\xi =\epsilon x+\beta y+\sigma z-\frac{k{t}^{\alpha }}{\alpha },$$

(51)

where $\epsilon ,\beta ,\sigma ,$ and k are constants. After performing algebraic manipulations, Equation (2) is transformed into the ODE in the variable $U=U\left(\xi \right)$ as

$$-kU{}^{\prime }+\epsilon {\delta }_1{U}^2{U}^{\prime }+\epsilon \left({\epsilon }^2{\delta }_2+\left({\beta }^2+{\sigma }^2\right){\delta }_3\right)U{}^{‴}=0,$$

(52)

where the prime notation (${}^{\prime }$) represents the ordinary derivative with respect to $\xi$. Integrating (52) with respect to $\xi$, we have the following ODE:

$$-kU+\frac{\epsilon {\delta }_1}{3}{U}^3+\epsilon \left({\epsilon }^2{\delta }_2+\left({\beta }^2+{\sigma }^2\right){\delta }_3\right)U{}^{″}+c=0,$$

(53)

where c is a constant of integration. Using the solution form (12) of the technique, the general solution $U\left(\xi \right)$ of (53) has the degree N. After balancing the highest order derivative $U{}^{″}$ in (53) with the nonlinear term of the highest order, i.e., $U^3$, we obtain $N=1$. Consequently, the solution of Equation (53) has the following form:

$$\begin{array}{c}\hfill U\left(\xi \right)=a_0+a_1\left(\frac{G^{\prime }}{G^2}\right)+b_1{\left(\frac{G^{\prime }}{G^2}\right)}^{-1},\end{array}$$

(54)

where $a_0,a_1,$ and $b_1$ are unknown constants that will be determined. Replacing Equation (54) into Equation (53), along with Equation (13), and then collecting all the coefficients of similar power of ${\left(G^{\prime }/G^2\right)}^i,(i=0,\pm 1,\pm 2,...)$, and ultimately setting these resulting coefficients to zero, we have the following system of algebraic equations in $a_0,a_1,b_1,\epsilon ,\beta ,\sigma ,k:$

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\left(\frac{G^{\prime }}{G^2}\right)}^{-3}:& 2{\mu }^2{\epsilon }^3{b}_1+\frac{\delta \epsilon b_1^3}{3}+2{\beta }^2{\mu }^2\epsilon b_1+2{\mu }^2{\sigma }^2\epsilon b_1=0,\hfill \\ \hfill {\left(\frac{G^{\prime }}{G^2}\right)}^{-2}:& \delta \epsilon a_0{b_1}^2=0,\hfill \\ \hfill {\left(\frac{G^{\prime }}{G^2}\right)}^{-1}:& 2{\beta }^2\lambda \mu \epsilon b_1+2\mu \lambda {\sigma }^2\epsilon b_1+2\mu \lambda {\epsilon }^3{b}_1+\delta \epsilon a_0^2{b}_1+\delta \epsilon a_1{b}_1^2-k{b}_1=0,\hfill \\ \hfill {\left(\frac{G^{\prime }}{G^2}\right)}^0:& \frac{\delta \epsilon a_0^3}{3}-k{a}_0+c+2\delta \epsilon a_0{a}_1{b}_1=0,\hfill \\ \hfill \left(\frac{G^{\prime }}{G^2}\right):& 2{\beta }^2\lambda \mu \epsilon a_1+2\lambda \mu {\sigma }^2\epsilon a_1+2\lambda \mu {\epsilon }^3{a}_1+\delta \epsilon a_0^2{a}_1+\delta \epsilon a_1^2{b}_1-k{a}_1=0,\hfill \\ \hfill {\left(\frac{G^{\prime }}{G^2}\right)}^2:& \delta \epsilon a_0{a}_1^2=0,\hfill \\ \hfill {\left(\frac{G^{\prime }}{G^2}\right)}^3:& 2{\lambda }^2{\epsilon }^3{a}_1+\frac{\delta \epsilon a_1^3}{3}+2{\beta }^2{\lambda }^2\epsilon a_1+2{\lambda }^2{\sigma }^2\epsilon a_1=0.\hfill \end{array}\end{array}$$

(55)

Solving the algebraic system in (55) with the help of Maple, we have the following three results for the exact solutions of (2):

Result 1:

$$\begin{array}{cc}\hfill & a_0=0,a_1=0,b_1=\pm \sqrt{-\frac{6\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}}\mu ,k=2\lambda \mu \epsilon \left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\delta }_2{\epsilon }^2\right),c=0,\hfill \end{array}$$

(56)

where ${\delta }_1,{\delta }_2,{\delta }_3,\epsilon ,\beta ,\sigma$ are arbitrary constants and $\mu ,\lambda$ are arbitrary integers. Substituting (56) into solution (54), along with the ratio $\left(G^{\prime }/G^2\right)$, in (15)–(17), the exact solutions of (2) with $\xi =\epsilon x+\beta y+\sigma z-\frac{k{t}^{\alpha }}{\alpha },$ where k is expressed in (56), are described depending on the values of $\mu ,\lambda$ as follows:

If $\mu \lambda >0$, then the trigonometric function solution of (2), which is formulated using the ratio (15), can be written as

$$\begin{array}{c}\hfill u_1^1(x,y,z,t)=\pm \sqrt{-\frac{6\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}}\left(\frac{\sqrt{\mu \lambda }\left(Dcos\left(\sqrt{\mu \lambda }\xi \right)-Csin\left(\sqrt{\mu \lambda }\xi \right)\right)}{Ccos\left(\sqrt{\mu \lambda }\xi \right)+Dsin\left(\sqrt{\mu \lambda }\xi \right)}\right),\end{array}$$

(57)

where $\frac{\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}<0$ and $C,D$ are arbitrary nonzero constants.

If $\mu \lambda <0$, then the exponential function solution of (2), which is formulated using the ratio (16), can be expressed as

$$\begin{array}{c}\hfill u_2^1(x,y,z,t)=\mp \sqrt{-\frac{6\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}}\left(\frac{\lambda \mu \left(C{\mathrm{e}}^{2\xi \sqrt{\left|\lambda \mu \right|}}-D\right)}{\sqrt{\left|\lambda \mu \right|}\left(C{\mathrm{e}}^{2\xi \sqrt{\left|\lambda \mu \right|}}+D\right)}\right),\end{array}$$

(58)

where $\frac{\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}<0$ and $C,D$ are arbitrary nonzero constants.

If $\mu =0,\lambda \ne 0$, then the rational function solution of (2), which is constructed using the ratio (17), can be consequently shown as

$$\begin{array}{c}\hfill u_3^1(x,y,z,t)=0.\end{array}$$

(59)

Result 2:

$$\begin{array}{c}\hfill a_0=0,a_1=\pm \sqrt{-\frac{6\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}}\lambda ,b_1=0,k=2\lambda \mu \epsilon \left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\delta }_2{\epsilon }^2\right),c=0,\end{array}$$

(60)

where ${\delta }_1,{\delta }_2,{\delta }_3,\epsilon ,\beta ,\sigma$ are arbitrary constants and $\mu ,\lambda$ are arbitrary integers. Substituting (60) into solution (54), along with the ratio $\left(G^{\prime }/G^2\right)$, in (15)–(17), the explicit exact solutions of (2) with $\xi =\epsilon x+\beta y+\sigma z-\frac{k{t}^{\alpha }}{\alpha },$ where k is expressed in (60), are exhibited depending on the values of $\mu ,\lambda$ as follows:

If $\mu \lambda >0$, then the trigonometric function solution of (2), which is formulated using the ratio (15), can be expressed as

$$\begin{array}{c}\hfill u_1^2(x,y,z,t)=\pm \sqrt{-\frac{6\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}}\left(\frac{\sqrt{\mu \lambda }Ccos\left(\sqrt{\mu \lambda }\xi \right)+Dsin\left(\sqrt{\mu \lambda }\xi \right)}{Dcos\left(\sqrt{\mu \lambda }\xi \right)-Csin\left(\sqrt{\mu \lambda }\xi \right)}\right),\end{array}$$

(61)

where $\frac{\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}<0$ and $C,D$ are arbitrary nonzero constants.

If $\mu \lambda <0$, then the exponential function solution of (2), which is constructed using the ratio (16), can be demonstrated as

$$\begin{array}{c}\hfill u_2^2(x,y,z,t)=\mp \sqrt{-\frac{6\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}}\left(\frac{\sqrt{\left|\lambda \mu \right|}\left(C{\mathrm{e}}^{2\xi \sqrt{\left|\lambda \mu \right|}}+D\right)}{C{\mathrm{e}}^{2\xi \sqrt{\left|\lambda \mu \right|}}-D}\right),\end{array}$$

(62)

where $\frac{\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}<0$ and $C,D$ are arbitrary nonzero constants.

If $\mu =0,\lambda \ne 0$, then the rational function solution of (2), which is formulated using the ratio (17), can be expressed as

$$\begin{array}{c}\hfill u_3^2(x,y,z,t)=\mp \sqrt{-\frac{6\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}}\left(\frac{C}{C\xi +D}\right),\end{array}$$

(63)

where $\frac{\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}<0,$$\xi =\epsilon x+\beta y+\sigma z$ and $C,D$ are arbitrary nonzero constants.

Result 3: The set of the parameter values for this result is separated as two sub-categories, namely, Result 3.1 and Result 3.2. In order to prevent confusion from selecting the sign in front of each exact solution formulated using the parameter values in this result, then

Table 1. Sign of $u(x,y,z,t)$ selected from ± or ∓ in front of $u(x,y,z,t)$.

Result	Sign of ${\mathit{a}}_1$	Sign of ${\mathit{b}}_1$	Sign Selected from ± or ∓ in Front of $\mathit{u}(\mathit{x},\mathit{y},\mathit{z},\mathit{t})$
			$\mathbf{\mu }\mathbf{\lambda }>\mathbf{0}$	$\mathbf{\mu }\mathbf{\lambda }<\mathbf{0}$	$\mathbf{\mu }=\mathbf{0},\mathbf{\lambda }\ne \mathbf{0}$
3.1	+	+	−	−	−
	−	−	+	+	+
3.2	+	−	+	+	−
	−	+	−	−	+

shows the correct signs of $u(x,y,z,t)$ selected from ± or ∓ in front of $u(x,y,z,t)$ in each case.

Result 3.1:

$$\begin{array}{c}\hfill \begin{array}{cc}& a_0=0,a_1=\pm \sqrt{-\frac{6\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}}\lambda ,b_1=\pm \frac{\mu a_1}{\lambda },\hfill \\ & k=-4\epsilon \lambda \mu \left({\delta }_2{\epsilon }^2+{\delta }_3\left({\beta }^2+{\sigma }^2\right)\right),c=0,\hfill \end{array}\end{array}$$

(64)

where ${\delta }_1,{\delta }_2,{\delta }_3,\epsilon ,\beta ,\sigma$ are arbitrary constants and $\mu ,\lambda$ are arbitrary integers. Substituting (64), along with the ratio $\left(G^{\prime }/G^2\right)$, in (15)–(17), depending upon the values of $\mu ,\lambda$ into the solution form (54), the exact traveling wave solutions of (2) with $\xi =\epsilon x+\beta y+\sigma z-\frac{k{t}^{\alpha }}{\alpha },$ where k is defined in (64), are exhibited as follows:

When $\mu \lambda >0$, the trigonometric function solution of (2), which is formulated using the ratio (15), can be expressed as

$$\begin{array}{c}\hfill u_1^{3,1}(x,y,z,t)=\mp \frac{\sqrt{-\frac{6\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}}\sqrt{\mu \lambda }(C^2+D^2)}{sin\left(\sqrt{\mu \lambda }\xi \right)\left(C^2-D^2\right)cos\left(\sqrt{\mu \lambda }\xi \right)-2CD{\left(cos\left(\sqrt{\mu \lambda }\xi \right)\right)}^2+CD},\end{array}$$

(65)

where $\frac{\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}<0$ and $C,D$ are arbitrary nonzero constants.

When $\mu \lambda <0$, the exponential function solution of (2), which is constructed using the ratio (16), can be demonstrated as

$$\begin{array}{cc}\hfill u_2^{3,1}(x,y,z,t)=& \mp \sqrt{-\frac{6\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}}\hfill \\ \hfill & \times \left(\frac{\mu \lambda {\left(C{\mathrm{e}}^{2\sqrt{\left|\lambda \mu \right|}\xi }-D\right)}^2+{\left(C{\mathrm{e}}^{2\sqrt{\left|\lambda \mu \right|}\xi }+D\right)}^2\left|\lambda \mu \right|}{\sqrt{\left|\lambda \mu \right|}\left(C{\mathrm{e}}^{2\sqrt{\left|\lambda \mu \right|}\xi }-D\right)\left(C{\mathrm{e}}^{2\sqrt{\left|\lambda \mu \right|}\xi }+D\right)}\right),\hfill \end{array}$$

(66)

where $\frac{\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}<0$ and $C,D$ are arbitrary nonzero constants.

When $\mu =0,\lambda \ne 0$, the rational function solution of (2), which is obtained using the ratio (17), can be written as

$$\begin{array}{c}\hfill u_3^{3,1}(x,y,z,t)=\mp \sqrt{-\frac{6\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}}\left(\frac{C}{C\xi +D}\right),\end{array}$$

(67)

where $\frac{\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}<0,$ $\xi =\epsilon x+\beta y+\sigma z$ and $C,D$ are arbitrary nonzero constants.

Result 3.2:

$$\begin{array}{c}\hfill \begin{array}{cc}& a_0=0,a_1=\pm \sqrt{-\frac{6\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}}\lambda ,b_1=\pm \frac{\mu a_1}{\lambda },\hfill \\ & k=8\epsilon \lambda \mu \left({\delta }_2{\epsilon }^2+{\delta }_3\left({\beta }^2+{\sigma }^2\right)\right),c=0,\hfill \end{array}\end{array}$$

(68)

where ${\delta }_1,{\delta }_2,{\delta }_3,\epsilon ,\beta ,\sigma$ are arbitrary constants and $\mu ,\lambda$ are arbitrary integers. Substituting (68), along with the ratio $\left(G^{\prime }/G^2\right)$, in (15)–(17), depending upon the values of $\mu ,\lambda$ into the solution form (54), the exact traveling wave solutions of (2) with $\xi =\epsilon x+\beta y+\sigma z-\frac{k{t}^{\alpha }}{\alpha },$ where k is defined in (68), are described as follows:

When $\mu \lambda >0$, the trigonometric function solution of (2), which is formulated using the ratio (15), can be expressed as

$$\begin{array}{cc}\hfill u_1^{3,2}(x,y,z,t)=& \pm \sqrt{\mu \lambda }\sqrt{-\frac{6\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}}\hfill \\ \hfill & \times \left(\frac{2\left(C^2-D^2\right){\left(cos\left(\sqrt{\mu \lambda }\xi \right)\right)}^2-4CDcos\left(\sqrt{\mu \lambda }\xi \right)sin\left(\sqrt{\mu \lambda }\xi \right)-C^2+D^2}{sin\left(\sqrt{\mu \lambda }\xi \right)\left(C^2-D^2\right)cos\left(\sqrt{\mu \lambda }\xi \right)+2CD{\left(cos\left(\sqrt{\mu \lambda }\xi \right)\right)}^2-CD}\right),\hfill \end{array}$$

(69)

where $\frac{\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}<0$ and $C,D$ are arbitrary nonzero constants.

When $\mu \lambda <0$, the exponential function solution of (2), which is constructed using the ratio (16), can be shown as

$$\begin{array}{cc}\hfill u_2^{3,2}(x,y,z,t)=& \pm \sqrt{-\frac{6\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}}\hfill \\ \hfill & \times \left(\frac{\mu \lambda {\left(C{\mathrm{e}}^{2\sqrt{\left|\lambda \mu \right|}\xi }-D\right)}^2-{\left(C{\mathrm{e}}^{2\sqrt{\left|\lambda \mu \right|}\xi }+D\right)}^2\left|\lambda \mu \right|}{\sqrt{\left|\lambda \mu \right|}\left(C{\mathrm{e}}^{2\sqrt{\left|\lambda \mu \right|}\xi }-D\right)\left(C{\mathrm{e}}^{2\sqrt{\left|\lambda \mu \right|}\xi }+D\right)}\right),\hfill \end{array}$$

(70)

where $\frac{\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}<0$ and $C,D$ are arbitrary nonzero constants.

When $\mu =0,\lambda \ne 0$, the rational function solution of (2), which is formulated using the ratio (17), can be written as

$$\begin{array}{c}\hfill u_3^{3,2}(x,y,z,t)=\mp \sqrt{-\frac{6\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}}\left(\frac{C}{C\xi +D}\right),\end{array}$$

(71)

where $\frac{\left({\beta }^2{\delta }_3+{\sigma }^2{\delta }_3+{\epsilon }^2{\delta }_2\right)}{{\delta }_1}<0,$$\xi =\epsilon x+\beta y+\sigma z$ and $C,D$ are arbitrary nonzero constants.

5. Graphical Representations of the Selected Solutions

In this portion, we will manifest interesting graphical representations of the selected exact solutions of the $(2+1)$-dimensional conformable time integro-differential Sawada–Kotera Equation (1) and of the (3 + 1)-dimensional conformable time modified KdV–Zakharov–Kuznetsov (mKdV–ZK) Equation (2) obtained using the algorithm of the $(G^{\prime }/G^2)$-expansion method. The time-fractional order $\alpha$ for the equations is changed in order to study graphical behaviors of the exact solutions chosen from the previous section. Particularly, the values of the time-fractional order used for the following simulations are $\alpha =1,0.7$ and $0.3$. Solutions (41) and (46) of Equation (1) and solution (65) of Equation (2) are selected to present in terms of 3D, 2D, and contour plots according to the values of $\alpha$. All of the 3D solution graphs of (1) and (2) are portrayed on the domain $\left\{\right(x,y,t):0\le x,t\le 10,y=1\}$ and $\left\{\right(x,y,z,t):0\le x,t\le 10,y=z=1\}$, respectively. The 2D solution graphs, demonstrating a relation between $u\left(x\right)$ and x, of (1) and (2) are depicted on $\left\{\right(x,y,t):0\le x\le 10,y=t=1\}$ and $\left\{\right(x,y,z,t):0\le x\le 10,y=z=t=1\}$, respectively. Moreover, the contour plots, showing a 3D surface by plotting constant u slices on a 2D plane, are drawn to connect the $(x,t)$ coordinates when the values of u are given and $y=1$ for (1) but $y=z=1$ for (2). In addition, physical descriptions of the displayed graphs will be mentioned in this section.

Figure 1 and Figure 2 show the solution graphs of the exact solutions $u_3^2(x,y,t)$ in (41) and $u_1^3(x,y,t)$ in (46) for problem (1), respectively. They are unfolded in different aspects, i.e., the 3D, 2D, and contour plots. Varying the values of $\alpha =1,0.7,0.3$, the solutions (41) and (46) are evaluated using the parameter sets $\{c=2,\mu =0,\lambda =1,C=0.5,D=-10\}$ and $\{c=1300,\mu =2,\lambda =1,C=1,D=1\}$, respectively, to plot their graphs on the domains. Particularly, Figure 1a–c shows the 3D, 2D, and contour plots for solution (41), respectively, when $\alpha =1$. Figure 1d–f and Figure 1g–i are plotted in the same manner as before except using $\alpha =0.7$ and $\alpha =0.3,$ respectively. By classifying the shapes of the 3D and 2D graphs in Figure 1, it can be identified that solution (41) is a singular single-soliton solution that is a solitary wave with discontinuous derivatives occurring at some domain regions, as observed in the contour plots of Figure 1. In addition, Figure 2a–c shows the 3D, 2D, and contour plots for solution (46), respectively, when $\alpha =1$. Figure 2d–f and Figure 2g–i are drawn in a similar manner to the above plots except using $\alpha =0.7$ and $\alpha =0.3,$ respectively. As noticed in the 3D and 2D graph structures in Figure 2, their physical behavior is considered as a singular periodic wave solution (or, a singular wavetrain), which is spatiotemporal oscillations with discontinuous derivatives. It can be roughly observed from the 3D graphs that the number of oscillations of the singular periodic wave solutions gradually increases as the fractionality $\alpha \in (0,1)$ decreases.

The solution graphs of the exact solution $u_1^{3,1}(x,y,z,t)$ in (65) for problem (2) are presented in Figure 3. To be clear, exact solution (65), in which the top sign of ∓ is chosen, is computed utilizing the parameter set $\{{\delta }_1={\delta }_2=-1,{\delta }_3=1,\epsilon =\beta =\sigma =1,\mu =2,\lambda =1,C=0.5,D=-1\}$ to plot its 3D, 2D, and contour graphs on the domains according to the used values of $\alpha$. Specifically, Figure 3a–c shows the 3D, 2D, and contour graphs of the solution (65) when $\alpha =1$ is used. However, the 3D, 2D, and contour plots of the solution (65) when $\alpha =0.7$ and $\alpha =0.3$ are exhibited in Figure 3d–f and Figure 3g–i, respectively. The physical behavior of these graphs is characterized as a singularly double periodic wave solution. The significant part of the doubly periodic wave solution represents a traveling wave whose envelope of emerging oscillations is bounded by a pattern periodic in both time and space. In addition, the number of oscillations of this solution type is inversely proportional to the value of $\alpha$.

6. Conclusions and Future Work

In our study, we have determined closed-form traveling wave solutions for the two nonlinear conformable evolution equations, which are the $(2+1)$-dimensional conformable time integro-differential SK Equation (1) and the (3 + 1)-dimensional conformable time mKdV–ZK Equation (2) by means of the $(G^{\prime }/G^2)$-expansion method. After eliminating the trivial and disqualified solutions, Equations (1) and (2) have three main results; each result provides the following three types of solutions: trigonometric, exponential (or, equivalently, hyperbolic), and rational function solutions. All of the exact solutions obtained in this paper were substituted back into their corresponding equations with the help of the Maple package program and their satisfactions confirm the validity of the solutions expressed in the current article. After visualizing the graphs of some solutions, they present some physical behaviors such as the singular single-soliton solution, the singular periodic wave solution and the singularly double periodic wave solution. These characteristics of the solutions are favorable for investigating certain nonlinear phenomena arising in physics, applied mathematics, and engineering. In particular, the soliton is a self-reinforcing wave packet maintaining its shape while propagating at a constant velocity. In other words, solitons are unscathed in shape and speed by a collision with other solitons and are often studied in quantum mechanics, nuclear physics, and waves along a weakly anharmonic mass-spring chain. Moreover, periodic traveling waves play a fundamental role in several mathematical physics including self-oscillatory systems, reaction–diffusion–advection systems, and excitable chemical reactions. Specifically, the family of doubly periodic wave solutions is of great importance in several physical phenomena such as modulation instability applied to the classical nonlinear Schrödinger equation (NLSE) and applications both in optics and deep water waves [64]. Since Equation (1), involving the time conformable partial derivative, is first proposed, then its obtained solutions are new and informed here for the first time. Equation (2) is an extension of the $(3+1)$-dimensional conformable time mKdV–ZK equation in Equation (48) of [65] for which the coefficients ${\delta }_2$ and ${\delta }_3$ are added. It is worth comparing our outcomes for (2) and the exact solutions of Equation (48) of [65] as follows: In [65], the first integral method and the functional variable method were used to solve the equation for which the trigonometric and hyperbolic function solutions were established. The mathematical structures of their results agree with those of our solutions for (2) except the $(G^{\prime }/G^2)$-expansion method additionally provides the rational function solutions for the equation. In conclusion, the performance of the method is direct, concise, reliable, and effective, and the method gives some interestingly particular types of solutions. Therefore, we deduce that the proposed methods can be extensively employed to solve many conformable NPDEs arising in the theory of solitons or other physics and engineering fields. Lastly, future studies could fruitfully explore the use of the $(G^{\prime }/G^2)$-expansion method further by applying it to the proposed problems with an extension of the spatiotemporal conformable partial derivatives or to NLEEs involving with sequential conformable partial derivatives.