# Applications of Homogenous Balance Principles Combined with Fractional Calculus Approach and Separate Variable Method on Investigating Exact Solutions to Multidimensional Fractional Nonlinear PDEs

Authors: Ruichao Ren; Shunli Zhang; Weiguo Rui
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

We investigate the exact solutions of multidimensional time-fractional nonlinear PDEs (fnPDEs) in this paper. In terms of the fractional calculus properties and the separate variable method, we present a new homogenous balance principle (HBP) on the basis of the (1 + 1)-dimensional time fnPDEs. Taking advantage of the new types of HBP together with fractional calculus formulas that subtly avoid the chain rule, the fnPDEs can be reduced to spatial PDEs, and then we solve these PDEs by the fractional calculus method and the separate variable approach. In this way, some new type exact solutions of the certain time-fractional (2 + 1)-dimensional KP equation, (3 + 1)-dimensional Zakharov–Kuznetsov (ZK) equation, and Jimbo–Miwa (JM) equation are explicitly obtained under both Riemann–Liouville derivatives and Caputo derivatives. The dynamical analysis of solutions is shown by numerical simulations of taking property parameters as well.


---

## Body

## 1. Introduction

The fractional PDEs (fPDEs) have been more and more widely followed with interest up to now since they can be used to accurately describe many nonlinear stochastic phenomena which depend on both time instant and the previous time history in the real-time problem [1–4]. Models set up from the fPDEs play very important roles in a range of scientific fields, such as viscoelastic flow [5, 6], signal processing [7, 8], control systems [9, 10], material diffusion including normal diffusion and anomalous diffusion(superdiffusion, subdiffusion, fast diffusion, and slow diffusion) [11–15], biological mathematics [16, 17], and magnetohydrodynamics (MHD) [18, 19]. Similarly, as some time-fractional PDEs, the space fractional models are also very frequently used in some elastic materials (see [20–24]) in the field of mechanical engineering. Therefore, investigating their solutions also draws much attention from the mathematical and physical points of view, and it can help to concisely characterize and well understand the qualitative features of the concerned phenomena and nonlinear processes in various areas of natural science and engineering, which involved the ubiquitous time memory effects [20, 21] (including some short memories) and space viscoelastic effects (see [22–26]), in particular many complex random walks and material motions in the microscopic space, such as dynamical behaviors of fractional diffusion, particles spread in heat bath, and soft matter interaction with viscoelasticity.In the past few years, there were excessive studies on the (1 + 1)-dimensional fractional nonlinear PDEs(fnPDEs), and a lot of excellent tools were used for solving them, which include Adomian decomposition method [27, 28], homotopy analysis method (HAM) [29, 30], fractional variational method [31, 32], Lie symmetry method [33–40], invariant subspace method [40–44], and homogenous balance principle (HBP) [45, 46]. Recently, a few of the above methods were also applied to solve several (2 + 1)-dimensional fnPDEs [47–50]. Although approximate analytic solutions or exact solutions of some (1 + 1)-dimensional fnPDEs and (2 + 1)-dimensional fnPDEs can be successfully obtained by using the above methods, this is far from enough, and these methods still have many limitations in solving more complex multidimensional fnPDEs. In this paper, we suggest a new technique to solve the following type of fnPDEs with certain time-fractional derivatives:(1)∂tαu=Fx,t,u,∂xu,…,∂xmu,where x=x1,…,xN,∂xku=∂i1+⋯+iNu/∂x1i1…∂xNiN,1≤k≤m,1≤i1+iN≤k, and the index 0<α≤1. The α-th order time derivative α∈0,1 is well defined as an abnormal derivative in the real applications for explaining short memories of evolutional physical systems; there are several kinds of definitions of fractional derivatives, and the two frequently used, classical, and very widely influential definitions are still Riemann–Liouville definition and Caputo definition (see Definitions 1 and 2). Indeed, with the help of Riemann–Liouville derivative and Caputo derivative, the results derived from the time-fractional PDE models are more precise and more general in nature than those of the integer-order ones, while the smaller the α is, the faster effective time memory enjoys [20, 21]. When α=1, expressions (3) and (4) are exactly in accordance with the classical derivatives ∂u/∂t; however, compared with the Riemann–Liouville type of derivative, the Caputo type of derivative possesses weaker singularity for handling some fractional initial problems.Moreover, it is necessary to point out that the two travelling wave transformationsu=Uξξ=∑i=1Nkixi−ωtα/Γ1+α and u=Uξξ=∑i=1Nkixi−ωt, which were often used to reduce the integer-order PDE to ODE, are actually not valid for the fractional-order one since it had been successfully verified that the compound fractional derivatives disagree with the following chain rule (refer to [45, 46, 51, 52]):(2)Dxαfgx=f′gxDxαgx=Dgxαfgxg′xα,thus, we hardly take the fractional derivative of compound function straightly in accordance with the classical chain rule, and the normal chain rule (see [51]) of fractional derivatives was ineffectively applied to solve equation (1) since it contains infinite series. That is to say, the exact solutions of the compound function type of equation (1) were impossible to be obtained by the invalid chain rule (2), and we hardly find the travelling wave solutions (even the exact soliton solutions) of fnPDEs via the two transformations mentioned above or some other complex transformations such as Darboux transformation [53], bilinear method [54], and F-expansion method [55].To the best of our knowledge, avoiding the invalid fractional chain rule (2), there is no more direct and effective method to obtain the exact solutions for equation (1). The exact solutions of higher-dimensional nonlinear PDEs with time-fractional derivatives were not well obtained. The main difficulty is how to construct solutions to reduce the multidimensional fnPDEs to the classical spatial PDEs. Inspired by the previous homogenous balance principle (HBP) [45–47] including fractional calculus method [40–44, 56] and separate variable approach, we improve the way [45–47] and introduce a new type of HBPs for (1) so that the solutions can be assumed as a general separated variable form and (1) can be reduced to spatial PDEs, and then new type exact solutions of some multidimensional fnPDEs will be successfully obtained by solving these reduced PDEs in the variable separate way under both Riemann–Liouville derivatives and Caputo derivatives.The main contents of this paper are organized as follows. The definitions and properties of the fractional calculus and fractional Laplace transformation are briefly described in Section2. In Section 3, based on the fractional derivative formulas and the method of separate variable, the new homogenous balance principle (HBP) is suggested for (1), and in this way, the certain time-fractional (2 + 1)-dimensional KP equation, (3 + 1)-dimension Zakharov–Kuznetsov (ZK) equation, and (3 + 1)-dimension Jimbo–Miwa (JM) equation are reduced to the spatial PDEs and explicitly solved in general separated variable forms; we can see that some of these solutions possess new types including some arbitrary functions which were never attained before by other way, which means more general solutions are obtained. Furthermore, there are real differences between the Riemann–Liouville case and the Caputo case: the singularity occurs in solutions under Riemann–Liouville derivatives but no singularity appears in solutions under Caputo derivatives. The dynamical profiles of these solutions are displayed as can be seen from Figures 1–13 with property parameters, and we also analyze the long time behaviors for some of them as well. The last section is the conclusions of our works.Figure 1
Dynamical profiles of solution (28) on the (x, t, u)-plane when (a)α=1/3 and (b)α=2/3.
(a)
(b)Figure 2
Dynamical profiles of solution (35) on (a) the (x, t, u)-plane and (b) the (y, t, u)-plane.
(a)
(b)Figure 3
Dynamical profiles of solution (36) on (a) the (x, t, u)-plane and (b) the (y, t, u)-plane.
(a)
(b)Figure 4
Dynamical profiles of solution (37) on (a) the (x, t, u)-plane and (b) the (y, t, u)-plane.
(a)
(b)Figure 5
Dynamical profiles of solution (49) on (a) the (x, t, u)-plane, (b) the (y, t, u)-plane, and (c) the (z, t, u)-plane.
(a)
(b)
(c)Figure 6
Dynamical profiles of solution (54) on (a) the (y, t, u)-plane and (b) the (z, t, u)-plane.
(a)
(b)Figure 7
Dynamical profiles of solution (66) on (a) the (x, t, u)-plane, (b) the (y, t, u)-plane, and (c) the (z, t, u)-plane.
(a)
(b)
(c)Figure 8
Dynamical profiles of solution (75) on (a) the (x, t, u)-plane, (b) the (y, t, u)-plane, and (c) the (z, t, u)-plane.
(a)
(b)
(c)Figure 9
Dynamical profiles of solution (76) on the (y, t, u)-plane.Figure 10
Dynamical profiles of solution (83) on the (x, t, u)-plane.Figure 11
Dynamical profiles of solution (89) on (a) the (x, t, u)-plane and (b) the (z, t, u)-plane.
(a)
(b)Figure 12
Dynamical profiles of solution (91) on (a) the (x, t, u)-plane and (b) the (z, t, u)-plane.
(a)
(b)Figure 13
Dynamical profiles of solution (110) on the (x, t, u)-plane when (a)α=1/3 and (b)α=2/3.
(a)
(b)


## 2. Preliminaries

In this section, we recall some useful definitions, properties, and theorem.


### 2.1. Definitions and Properties of Two Type Fractional Derivatives

Definition 1.
The Riemann–Liouville fractional derivative of orderα>0 is defined by the following expression:(3)∂RLtαux,t=1Γn−α∂n∂tn∫0tt−τn−α−1ux,τdτ,n=α+1,∂nux,t∂tn,n=α.Definition 2.
The Caputo fractional derivative of orderα>0 is defined by the following expression:(4)∂Ctαux,t=1Γn−α∫0tt−τn−α−1∂nux,τ∂τndτ,n=α+1,∂nux,t∂tn,n=α.Proposition 1.
Some properties of the fractional derivative:(i)
DRtαc=ct−α/Γ1−α (c is a constant)(ii)
DCtαc=0 (c is a constant)(iii)
Dtα∑k=1nckfkt=∑k=1nckDtαfkt (ck are constants)(iv)
Dtαtμ=Γ1+μ/Γ1+μ−αtμ−α,μ≠α−1 holds for μ>−1 under the Riemann–Liouville case and μ>0 under the Caputo case (5).Proposition 2.
The rules of the partial fractional derivative of a separate variable form:(i)
Forφt≠Const,Dtαfxφt=fxDtαφt.(ii)
Forφt=1,(5)Dtαfx=fxt−αΓ1−α,under Riemann–Liouville derivative,0,under Caputo derivative.


### 2.2. Definitions of Mittag-Leffler Function and Some Properties

Definition 3.
The Mittag-Leffler function is defined by the infinity series expression:(6)Eα,βz=∑k=0∞zkΓαk+β,with α>0andβ>0.Proposition 3.
The properties of the Mittag-Leffler function:(i)
E1,1z=ez,(ii)
Dtγtαk+β−1Eα,βkλtα=tαk+β−γ−1Eα,β−γkλtα,(iii)
DCtαEα,1λtα=λEα,1λtα,(iv)
DRtαEα,1λtα=t−α/Γ1−α+λEα,1λtα,(v)
The Mittag-Leffler functionEα,1λtα is increasing as λ>0 while decreasing as λ<0 when α>0 (Definition 3).
The brief proof of (iii) and (iv) is given as follows.
For the Caputo derivatives, we have(7)DCtαEα,1λtα=DCtα∑k=0∞λtαkΓαk+1=∑k=0∞λkDCtαtkαΓαk+1=λ∑k=1∞λk−1Γαk+1/Γαk+1−αtk−1αΓαk+1=λ∑n=0∞λtαnΓαn+1=λEα,1λtα.
For the Riemann–Liouville derivatives, we obtain(8)DRtαEα,1λtα=DRtα∑k=0∞λtαkΓαk+1=∑k=0∞λkDRtαtkαΓαk+1=t−αΓ1−α+λ∑k=1∞λk−1Γαk+1/Γαk+1−αtk−1αΓαk+1=t−αΓ1−α+λ∑n=0∞λtαnΓαn+1=t−αΓ1−α+λEα,1λtα.
These achieve (iii) and (iv) in Proposition3.


## 2.1. Definitions and Properties of Two Type Fractional Derivatives

Definition 1.
The Riemann–Liouville fractional derivative of orderα>0 is defined by the following expression:(3)∂RLtαux,t=1Γn−α∂n∂tn∫0tt−τn−α−1ux,τdτ,n=α+1,∂nux,t∂tn,n=α.Definition 2.
The Caputo fractional derivative of orderα>0 is defined by the following expression:(4)∂Ctαux,t=1Γn−α∫0tt−τn−α−1∂nux,τ∂τndτ,n=α+1,∂nux,t∂tn,n=α.Proposition 1.
Some properties of the fractional derivative:(i)
DRtαc=ct−α/Γ1−α (c is a constant)(ii)
DCtαc=0 (c is a constant)(iii)
Dtα∑k=1nckfkt=∑k=1nckDtαfkt (ck are constants)(iv)
Dtαtμ=Γ1+μ/Γ1+μ−αtμ−α,μ≠α−1 holds for μ>−1 under the Riemann–Liouville case and μ>0 under the Caputo case (5).Proposition 2.
The rules of the partial fractional derivative of a separate variable form:(i)
Forφt≠Const,Dtαfxφt=fxDtαφt.(ii)
Forφt=1,(5)Dtαfx=fxt−αΓ1−α,under Riemann–Liouville derivative,0,under Caputo derivative.


## 2.2. Definitions of Mittag-Leffler Function and Some Properties

Definition 3.
The Mittag-Leffler function is defined by the infinity series expression:(6)Eα,βz=∑k=0∞zkΓαk+β,with α>0andβ>0.Proposition 3.
The properties of the Mittag-Leffler function:(i)
E1,1z=ez,(ii)
Dtγtαk+β−1Eα,βkλtα=tαk+β−γ−1Eα,β−γkλtα,(iii)
DCtαEα,1λtα=λEα,1λtα,(iv)
DRtαEα,1λtα=t−α/Γ1−α+λEα,1λtα,(v)
The Mittag-Leffler functionEα,1λtα is increasing as λ>0 while decreasing as λ<0 when α>0 (Definition 3).
The brief proof of (iii) and (iv) is given as follows.
For the Caputo derivatives, we have(7)DCtαEα,1λtα=DCtα∑k=0∞λtαkΓαk+1=∑k=0∞λkDCtαtkαΓαk+1=λ∑k=1∞λk−1Γαk+1/Γαk+1−αtk−1αΓαk+1=λ∑n=0∞λtαnΓαn+1=λEα,1λtα.
For the Riemann–Liouville derivatives, we obtain(8)DRtαEα,1λtα=DRtα∑k=0∞λtαkΓαk+1=∑k=0∞λkDRtαtkαΓαk+1=t−αΓ1−α+λ∑k=1∞λk−1Γαk+1/Γαk+1−αtk−1αΓαk+1=t−αΓ1−α+λ∑n=0∞λtαnΓαn+1=t−αΓ1−α+λEα,1λtα.
These achieve (iii) and (iv) in Proposition3.


## 3. HBP for Some Multidimensional Time-Fractional PDEs

In the beginning of this section, we improve the previous way to construct exact solutions to some multidimensional time-fractional nonlinear PDEs (1) according to ideas of homogenous balance and the fractional properties (5).The main procedures are stated as follows:Step 1. As can be seen from the properties of power function and Mittag-Leffler function in Section 2 and motivated by the HBP of (1 + 1)-dimensional time-fractional nonlinear PDE, for equation (1), we suppose the exact solutions as the following two types:(1)
For the Riemann–Liouville case:(9)u=a0tμ1+a1vx1,…,xNtμ2,μ1,μ2≠α−1.(2)
For the Caputo case:(10)u=a0+a1vx1,…,xNEα,1λtα.Hereμ1,μ2, and λ are undetermined parameters and vx1,…,xN is a undetermined function.Step 2. Substituting (9) or (10) into equation (1) by comparing the powers of t or the coefficients of Eα,1λtα, we have parameters μ1,μ2, or λ and the reduced spatial PDE system for vx1,…,xN:(11)vx1,…,xN=Φx,t,v,∂xv,…,∂xmv,where(12)x=x1,…,xN,∂xkv=∂i1+⋯+iNv∂x1i1…∂xNiN,1≤k≤m,1≤i1+iN≤k.Step 3. Solving the reduced PDE system (10) by using complex transformation ξ=c0+∑i=1Ncixi or a separate variable approach (in product or summary form) leads to the exact solutions.Remark 1.
Whenμ1,μ2=α−1, singular solutions will appear under the Riemann–Liouville derivative and while solutions appear, no singularity will be obtained in the Caputo sense. t should involve minus power, thereby the assumption (9) fits for the Riemann–Liouville derivative, and since the fractional derivative of Mittag-Leffler function includes an additional term t−α/Γ1−α under the Riemann–Liouville derivative, this will complicate the calculation of solutions, and the assumption (10) fits for the Caputo derivative.
We will study time-fractional (2 + 1)-dimensional KP equation, (3 + 1)-dimensional ZK equation, and JM equation in the following paragraphs.


### 3.1. HBP for the (2 + 1)-Dimensional Time-Fractional KP Equation

We first begin with the (2 + 1)-dimensional fractional KP equation:(13)∂tαu−6uux+uxxxx+3uyy=0,0<α≤1,which developed from the classical KdV equation, describes long dispersive wave propagating in two dimension in shallow water.


#### 3.1.1. Exact Solutions under Riemann–Liouville Derivatives

According to (9), we suppose the exact solutions as(14)u=a0tμ1+a1vx,ytμ2,where a0,a1,μ1,andμ2 are the undetermined constants.Balancing the coefficients oft-power yields the following two subcases:Subcase I.(15)μ2−α=μ1+μ2=2μ2.Thus,μ1=μ2=−α.Substituting (15) into (14) by comparing the coefficients of t−αandt−2α yields(16a)Ω0vx=6a0+a1vvxx+6a1vx2,Ω0=Γ1−αΓ1−2α,(16b)vxxxx+3vyy=0.Then, integrating (16a) corresponding to variable x once by selecting an integrate constant 0 leads to(17)dvdx=Ω0v6a0+a1v.This gives rise to the implicit relation:(18)6a0lnv+a1v=Ω0x+ϕy,whereϕy is an undetermined function.Calculating the derivatives of implicit function (18) reads that(19a)vxxxx=a0Ω046a12v2−8a0a1v+a02v1296a17v7+7a0a16v6+21a02a15v5+35a03a14v4+35a04a13v3+21a05a12v2+7a06a1v+a07,(19b)vyy=Ω0Ω0a0ϕ′y2+6a0+a1v2ϕ′′yv2v36a13v3+3a0a12v2+3a02a1v+a03.Plugging (19a) and (19b) into (15), we have(20)a0=0,ϕ′′y=0,which leads to the solution of (15)(21)vx,y=Ω06a1x+b1y+b2,with arbitrary constantsa1,b1,andb2  and the exact solution of (11)(22)ut,x,y=a0+Γ1−α6Γ1−2αx+b1y+b2t−α,α≠12,wherec1andc2 are the arbitrary constants.Subcase II. μ2−α=2μ2andμ1+μ2=μ2.Thus,μ1=0andμ2=−α.Similarly as above, substituting (13) into (11) by comparing the coefficients of t−αandt−2α yields(23a)Ω0vx=6a1vvxx+6a1vx2,Ω0=Γ1−αΓ1−2α,(23b)vxxxx+6a0vxx+3vyy=0.Then, integrating (23a) corresponding to variable x by choosing an integrate constant 0 once leads to(24)Ω0v=6a1vvx.Sincev≠0 (only nontrivial solutions are considered), we have(25)dvdx=Ω06a1,which give rises to(26)v=Ω06a1x+ϕy.Plugging (26) into (23b) yields(27)ϕy=c1y+c2.By rewritingb1=a1c1andb2=a2c2, we have the exact solution(28)ut,x,y=a0+Γ1−α6Γ1−2αx+b1y+b2t−α,α≠12,wherea0,b1,andb2 are the arbitrary constants.Remark 2.
We have the same type solutions under the complex transformationξ=c1x+c2y+c0, so this case is omitted here.


#### 3.1.2. Exact Solutions under Caputo Derivatives

According to (10), we assume the exact solutions formed as(29)u=a0+a1vx,yEα,1λtα.Substituting (29) into (11) by comparing the coefficients of Eα,1λtα leads to(30a)λa1vx+a1vxxxx+a0vxx+3a1vyy=0,(30b)vvxx=vx2.We consider the following two cases:Case I. By using the total differential rule, from (30b), we have(31)vvxx−vx2v2=vxvx=0.Integrating (31) leads to(32)vx=Cyv.Thus,(33)v=ϕyeCyx,whereϕyandCy are the undetermined functions.Plugging (33) into (30a) yields(34)ϕ′′y+cλ+c4+a0c23ϕy=0.Then, we have following three type solutions:(i)
δ=cλ+c4+a0c2<0.We obtainϕy=c1e−−δ/3y+c2e−δ/3y and the exact solution(35)ut,x,y=a0+c1ecx−−δ/3y+c2ecx+−δ/3yEα,1λtα,wherea0,c,c1,andc2 are the arbitrary constants.(ii)
δ=cλ+c4+a0c2>0.We obtain the exact solution(36)ut,x,y=a0+c1cosδ3y+c2sinδ3yecxEα,1λtα,wherea0,c,c1,andc2 are the arbitrary constants.(iii)
δ=cλ+c4+a0c2=0.We obtainϕy=c1y+c2 and the exact solution(37)ut,x,y=a0+c1y+c2ecxEα,1λtα,wherea0,c,c1,andc2 are the arbitrary constants.Case II. Using the complex differential transformation ξ=c1x+c2y+c0 on (30a) and (30b), we have(38a)c14v4ξ+a0c12+3c22v′′ξ+λc1v′ξ=0,(38b)vξv′′ξ=v′2ξ.Solving (38b) similarly as (25) yields(39)vξ=recξ.Taking (39) into (38a), we obtain c14c3+a0c12+3c22c+λc1=c0 and the exact solution by rewriting a=a1recc0,b1=cc1,b2=cc2(40)ut,x,y=a0+aeb1x+b2yEα,1λtα,wherea0,a,b1,andb2 are the arbitrary constants.Remark 3.
(33) contains product forms of separate variable v=fxgy, and for omitting the trivial solution, we do not consider the summary form of a separate variable v=fxgy.Remark 4.
The singularity only appears whenα=1/2 in solutions (22) and (28) in the Riemann–Liouville case.


#### 3.1.3. Dynamical Analysis of Exact Solution for Fractional KP Model

Under the Riemann–Liouville case:As can be seen from Figure1, taking parameters α=1/3 and α=2,a0=1,b1=2,b2=1,y=1, and interval x∈−10,10,t∈0,10, it is shown that solution (28) increases with increase in the spatial variable and decreases with increase in time. When x and y are fixed, the solution tends to 0 at an α decay rate as t⟶∞, the larger the α is, the faster the u decay is.Under the Caputo case:(i)
As can be seen from Figure2, choosing c1=1,c2=−2,c=1,λ=−2,α=1/2,δ=−1,y=2,x∈−10,10,t∈0,10, and x=2,y∈−10,10,t∈0,10, it is shown that solution (35) decreases with increase in time-space since slope ux<0 or uy<0; however, it is unbounded as x or y tends to ∞.(ii)
As we see from Figure3, selecting c1=1,c2=−2,c=2,λ=12,α=1/2,δ=−1,y=2,x∈−10,10,t∈0,10, and x=2,y∈−30,10,t∈0,10, solution (36) decreases as time increases, periodic bounded as y, t tends to ∞ if x is fixed while unbounded as x, and t tends to ∞ if y is fixed.(iii)
As we see from Figure4, choosing c1=1,c2=−2,c=1,λ=−1,α=1/2,δ=0,y=3,x∈−10,10,t∈0,10, and x=2,y∈−10,10,t∈0,10, solution (37) decreases as time increases and unbounded as time-space tend to ∞.


### 3.2. HBP for the (3 + 1)-Dimensional Fractional ZK Equation

#### 3.2.1. Exact Solutions under Riemann–Liouville Derivatives

The (3 + 1)-dimensional fractional ZK equation is investigated in this section.(41)∂tαu+puux+quzzz+ruxxz+suyyz=0,0<α≤1,which represented the acoustic dynamics in a magnetized plasma in three-dimensional space with a low pressure.The exact solution can be assumed as(42)u=a0tμ1+a1vx,y,ztμ2.Balancing the coefficients oft-power by plugging (42) into (41) yields the following two cases:Cases I. μ1−α=μ1+μ2=2μ2That is,(43)μ1=μ2=−α.By comparing the coefficients oft−αandt−2α, we have(44a)Ω0+pa1vx=0,Ω0=Γ1−αΓ1−2α,(44b)qvzzz+rvxxz+svyyz=0.By solving (44a) directly, we obtain(45)vx,y,z=−Ω0pa1x+ϕy,z.Substituting (45) into (44b) yields(46)qϕzzz+sϕyyz=0.Subcase I. Assume η=b1y+b2z+b0, then we have(47)qb22+sb12ϕ′′′η=0.Solving (47) leads to the following two solutions:(i)
Whenqb22+sb12=0,sq<0, then(48)ux,y,z,t=a0−Γ1−αpΓ1−2αx+ϕb1y+b2z+b0t−α,α≠12,wherea0andb0 are the arbitrary constants and ϕ is an arbitrary function.(ii)
Whenϕ=c2η2/2+c1η+c0, then(49)ux,y,z,t=a−Γ1−αpΓ1−2αx+c2b1y+b2z+b022+c1b1y+b2zt−α,α≠12,wherea=a0+c0+c1b0,c1,andc2 are the arbitrary constants.Subcase II. Note that (46) is a linear equation, and if sq<0, by the linear PDE method, we have(50)ϕy,z=f1y−2−sqsf1y2−sz2−sq+f2−sqsy+z,which gives rise to exact solution whensq<0:(51)ux,y,z,t=a0−Γ1−αpΓ1−2αx+a1f1y−2−sqsf1y2−sz2−sq+f2−sqsy+zt−α,α≠12,wherea0anda1 are the arbitrary constants and f1andf2 are the two arbitrary functions.Subcase III. According to the linear PDE method, we assume the solution of (46) as a product form of a separate variable ϕy,z=fygz, thus(52)qg′′′z=c0g′z,sf′′y+c0fy=0,wherec0 is an arbitrary constant.By solving them directly, we obtain the following four solutions (67)–(70) where a0,a1,c1,c2,c3,c4,andc5 are the arbitrary constants and α≠1/2:(i)
ϕ=c1sinc0/sy+c2cosc0/syc3+c4sin−c0/qz+c5cos−c0/qz,(53)ux,y,z,t=a0−Γ1−αpΓ1−2αx+a1c1sinc0sy+c2cosc0syc3+c4sin−c0qz+c5cos−c0qzt−α.(ii)
ϕ=c1sinc0/sy+c2cosc0/syc3+c4ec0/qz+c5e−c0/qz,(54)ux,y,z,t=a0−Γ1−αpΓ1−2αx+a1c1sinc0sy+c2cosc0syc3+c4ec0/qz+c5e−c0/qzt−α.(iii)
ϕ=c1e−c0/sy+c2e−−c0/syc3+c4sin−c0/qz+c5cos−c0/qz,(55)ux,y,z,t=a0−Γ1−αpΓ1−2αx+a1c1e−c0/sy+c2e−−c0/syc3+c4sin−c0qz+c5cos−c0qzt−α.(iv)
ϕ=c1e−c0/qy+c2e−−c0/qyc3+c4ec0/qz+c5e−c0/qz,(56)ux,y,z,t=a0−Γ1−αpΓ1−2αx+a1c1e−c0/sy+c2e−−c0/syc3+c4ec0/qz+c5e−c0/qzt−α.Subcase IV. We can also assume the solution of (47) as a summary form of a separate variable ϕy,z=fy+gz, then(57)g′′′z=0.Thus,(58)gz=c22z22+c1z+c0,which leads to the exact solution by rewritingb2=a1c22,b1=a1c1anda=a0+a1c0:(59)ux,y,z,t=a−Γ1−αpΓ1−2αx+a1fy+b2z22+b1zt−α,α≠12,wherea,a1,b2,andb1 are the arbitrary constants and fy is an arbitrary function.Cases II. μ2−α=2μ2andμ1−α=μ1+μ2.By balancing equation (41), we have(60)μ1=0,μ2=−α.Comparing the coefficients oft−αandt−2α yields(61a)a0ω0+pa0a1vx+qa1vzzz+ra1vxxz+sa1vyyz=0,ω0=1Γ1−α,(61b)Ω0v+pa1vvx=0,Ω0=Γ1−αΓ1−2α.Solving (61b) leads to the following form:(62)v=−Ω0pa1x+ϕy,z.Substituting (62) into (61a) leads to(63)qa1ϕzzz+sa1ϕyyz=a0Ω0−ω0.(i)
Considering the complex transformationη=b1y+b2z+b0, we obtain(64)a1b2qb22+sb12ϕ′′′η=a0Ω0−ω0.Solving (81a) and (81b) leads to(65)ϕη=a0Ω0−ω06a1b2qb22+sb12η3+c22η2+c1η+c0.Then, we obtain the exact solution(66)ut,x,y,z=a0+−Γ1−αpΓ1−2αx+a0Γ1−α/Γ1−2α−1/Γ1−α6b2qb22+sb12b1y+b2z+b03+a1c22b1y+b2z+b02+a1c1b1y+b2z+b0+a1c0t−α,α≠12,wherea0,b1,b2,a1,c0,c1,andc2 are the arbitrary constants.(ii)
Note that (63) is a linear equation, and by the linear PDE method, we have(67)ϕy,z=f1y+12f2−sq/sy+zsqa1−24f1y/2−sz/2−sq−sqqa1+a0Ω0−ω0sz3+3qa0Ω0−ω0zy212sqa1,which leads to the exact solution(68)ut,x,y,z=a0+−Ω0px+a1f1y+12f2−sq/sy+zsqa1−24f1y/2−sz/2−sq−sqqa1+a0Ω0−ω0sz3+3qa0Ω0−ω0zy212sqt−α,α≠12,wherea0anda1 are constants and f1andf2 are the two arbitrary functions.(iii)
Using summary form of the separate variableϕ=fy+gz yields(69)qa1g′′′z=a0Ω0−ω0.Solving (69) directly, we arrive at(70)gz=a0Ω0−ω06qa1z3+b22z2+b1z+b0,which gives rise to the exact solution(71)ut,x,y,z=a0+−Γ1−αpΓ1−2αx+a1fy+a06qΓ1−αΓ1−2α−1Γ1−αz3+a1b22z2+a1b1z+a1b0t−α,α≠12,wherea0,a1,b2,b1,andb0 are the arbitrary constants and fy is an arbitrary function.(iv)
Whena0=0, using the product form of the separate variable ϕ=fygz yields(72)qfyg′′′z+sf′′′ygz=0,which give rises to(73)f′′y=c0fy,qg′′′z+sc0gz=0.From (73), we have(74)gz=c1e−sc0q23/qz+e−sc0q23/2qzc2cos3−sc0q232qz+c3sin3−sc0q232qz,and following three solutions (75)–(77) where a1,b1,b2,c1,c2,andc3 are constants and α≠1/2:(i)
Ifc0>0, then(75)ut,x,y,z=−Γ1−αpa1Γ1−2αx+b1e−c0y+b2ec0yc1e−sc0q23/qz+e−sc0q23/2qzc2cos3−sc0q232qz+c3sin3−sc0q232qzt−α.(ii)
Ifc0<0, then(76)ut,x,y,z=−Γ1−αpa1Γ1−2αx+b1cos−c0y+b2sin−c0yc1e−sc0q23/qz+e−sc0q23/2qzc2cos3−sc0q232qz+c3sin3−sc0q232qzt−α.(iii)
Ifc0=0, then(77)ut,x,y,z=−Γ1−αpa1Γ1−2αx+b1y+b2c1z2+c2z+c3t−α.Remark 5.
In case I, we do not concern the complex transformξ=ax+by+cz+d for (48) contains solutions of this type.Remark 6.
In case II, we do not consider the complex transformξ=ax+by+cz+d for omitting the trivial solution.Remark 7.
The singularity only appears whenα=1/2 in solutions (53)–(56), (66), (68), (71), and (76)–(78) in the Riemann–Liouville case.Remark 8.
If taking form (10) into the fractional ZK model (41), we see that the nonlinear term uux will lead to trivial solutions, thus the Caputo derivative is not considered.


#### 3.2.2. Dynamical Analysis of Exact Solution for Fractional ZK Model

(i)
As can be seen from Figure5, taking α=3/5,a=3,a0=c0=1,c1=1/2,b0=2,b1=1,b2=−2,c2=2,t∈0,10,y=1.5,z=2.5,x∈−10,10, x=−2.5,z=2.5,y∈−10,10, and x=−2.5,y=1.5,z∈−10,10, it is shown that solution (49) increases with increase in spatial variables and decreases as time increases. When x, y, and z are fixed, solution (49) tends to 0 at α decay rate as t⟶∞, the larger the α is, the faster the u decay is.
(ii)
From Figure6, by taking α=2/5,a0=1/2,a1=2,c1=1.5,c2=−2,c3=−3,c4=−2.5,c5=1,c0=2,p=3,s=q=1/2,t∈0,30,x=−2.5,z=2.5,y∈−30,30 and t∈0,10,andx=−2.5,y=1.5,z∈−10,10, we see solution (54) is periodic bounded and tends to 0 at α decay rate as y and t tend to ∞ if x and z are fixed. The larger the α is, the faster the u decay is. When x and y are fixed, solution (54) decreases to −∞ as z increases to +∞.
(iii)
As can be seen from Figure7, selecting α=2/3,p=3,a0=2,a1=1,b1=2,b2=−1,b0=3,c0=1,c1=1/2,c2=2,q=s=1/2,t∈0,10,y=1.5,z=2.5,x∈−10,10 , x=−2.5,z=2.5,y∈−10,10, and x=−2.5,y=1.5,z∈−10,10, while from Figure 8 by choosing α=2/3,p=3,a1=1,b1=1,b2=−1.5,c1=1.5,c2=2,c3=−2.5,c0=2,s=8,q=2,t∈0,10,y=1.5,z=2.5,x∈−10,10, x=−2.5,z=2.5,y∈−10,10, and x=−2.5,y=1.5,z∈−10,10, it is shown that solutions (66) and (75) increase to ∞ or decrease to −∞ as x, y, or z increases to ∞, and ones also tends to 0 at α decay rate as t⟶∞ if x, y, and z are fixed. The larger the α is, the faster the u decay is.
(iv)
From Figure9, as we choose α=2/3,p=1,a1=3,b1=1,b2=−1.5,c0=2,c1=1,c2=2,c3=−2.5,s=8,q=2,x=z=−2.5,y∈−30,30,t∈0,30, the solution (76) is periodic bounded and has α decay rate as t⟶∞ if x and z are fixed, and one has the same dynamical properties as (54) if else.


### 3.3. HBP for the (3 + 1)-Dimensional Fractional JM Equation

The (3 + 1)-dimensional fractional JM equation will be researched in the following contents:(78)uxxxy+puyuxx+quxuxy+r∂tαuy−suxz=0,0<α≤1,which is developed from the second members of integrable systems of the classical KP hierarchy.


#### 3.3.1. Exact Solutions under Caputo Derivatives

Assume(79)u=a0+a1vx,y,zEα,1λtα.Injecting (79) into (78), we arrive at(80a)vxxxy+rλvy−svxz=0,(80b)pvyvxx+qvxvxy=0.(i)
We assume thatv=vξ,ξ=ax+by+cz+d, then (80a) and (80b) becomes the following ODEs:(81a)a3bv4ξ+rλbv′ξ−sacv′′ξ=0,(81b)p+qv′ξv′′ξ=0.Whenq=−p, (81b) is a free equation, and then solving (81a), we have following solutions by noting Δ=327r2λ2b3−4s3c3/b−9rλb3b2,(82)vξ=c0bλr+c1e121/3Δ2+122/3scb/6abΔξ+e−121/3Δ2+122/3scb/12abΔξc2cos3121/3Δ2−122/3scb12abΔξ+c3sin3121/3Δ2−122/3scb12abΔξ,which leads to the exact solution:(83)ut,x,y,z=a0+a1c0bλr+c1e121/3Δ2+122/3scb/6abΔax+by+cz+d+e−121/3Δ2+122/3scb/12abΔax+by+cz+dc2cos3121/3Δ2−122/3scb12abΔax+by+cz+d+c3sin3121/3Δ2−122/3scb12abΔax+by+cz+dEα,1λtα,wherea0,a1,c0,a,b,c,andd are the arbitrary constants.Remark 9.
In case ofq≠−p, from v′′ξ=0, we obtain v′ξ=0 which only leads to a trivial solution, thus we do not consider this case.(ii)
Assumevx,y,z=vη,z,η=ax+by+c, and from (80a) and (80b), we have(84a)a3bvηηηη+rλbvη−savηz=0,(84b)p+qvηvηη=0.Then (84a) and (84b) yields the following two subcases(1)
Ifq=−p, we only consider (84a).(i)
Note that it is a linear equation, and we supposevη,z=fηgz, which arrives at(85)a3bf4ηgz+rλbf′ηgz−saf′ηg′z=0.Using the separate variable method, we have(86a)g′z=cgz,(86b)a3bf′′′η+rλb−sacfη=c0.Solving (86a) reads that(87)gz=b0ecz,and the solution of (86b) and the exact solution of (78) is described as the following two cases:Ifrλb≠sac, then(88)fη=c0rbλ−sac+c1esac−rbλb23/abη+e−sac−rbλb2/2ab3ηc2cos3sac−rbλb232abη+c3sin3sac−rbλb232abη,and the exact solution is(89)ut,x,y,z=a0+a1b0c0rbλ−sac+c1esac−rbλb23/abax+by+c+e−sac−rbλb23/2abax+by+cc2cos3sac−rbλb232abx+by+ca+c3sin3sac−rbλb232abax+by+ceczEα,1λtα,wherea0,a1,b0,c0,a,b,c,c1,c2,andc3 are the arbitrary constants.Ifc=rbλ/sa, then(90)fη=c06η3+c22η2+c1η+c3,and the exact solution is(91)ut,x,y,z=a0+ac06ax+by+rbλsa3+c22ax+by+rbλsa2+c1ax+by+rbλsa+c3erbλ/sazEα,1λtα,wherea0,c0,a,b,c1,c2,andc3 are the arbitrary constants.(ii)
Solving (84a) directly by the linear PDE method reads that(92)vη,z=c−1/3c1ec1/3η+ba2c/s+brλ/saz+4e−1/2c1/3η+ba2c/s+brλ/sazc2sin3c1/32η+c3cos3c1/32η+fz,ut,x,y,z=a0+a1c−1/3c1ec1/3ax+by+c+ba2c/s+brλ/saz+4e−1/2c1/3ax+by+c+ba2c/s+brλ/sazc2sin3c1/32ax+by+c+c3cos3c1/32ax+by+c+fzEα,1λtα,wherea0,a1,a,b,c,c1,c2,andc3 are the arbitrary constants and fz is an arbitrary function.(iii)
Assumevη,z=fη+gz, and from (84a), we have(93)a3f4η+rλf′η=0,and then(94)a3f′′′η+rλfη=c0.By solving (94), we have(95)fη=c0rλ+c1e−rλ3/aη+erλ3/2aηc2cos3rλ32aη+c3sin3rλ32aη,and the exact solution(96)ut,x,y,z=a0+a1c0rλ+c1e−rλ3/aax+by+c+erλ3/2aax+by+cc2cos3rλ32aax+by+c+c3sin3rλ32aax+by+c+gzEα,1λtα,wherea0,a1,a,b,c,c1,c2,andc3 are the arbitrary constants and gz is an arbitrary function.(2)
Ifq≠−p, then for the nontrivial case, we have vηη=0, and this leads to(97)v=f1zη+f2z.Then, substituting (97) into (84a), we have(98)rλbf1z−saf1′z=0,which yields(99)f1z=c0erλb/saz.This gives rise to the exact solution(100)ut,x,y,z=a0+a1c0erλb/sazax+by+c+f2zEα,1λtα,wherea0,a1,c0,a,b,andc are the arbitrary constants and f2z is an arbitrary function.(III)
Assumev=fx,z+gy,z.Then, from (84a) and (84b), we arrive at(101a)gyfxx=0,(101b)rλgy−sfxz=0.For the nontrivial case, we obtain(102)fxx=0⇒f=f1zx+f2z.Injecting (102) into (101b) yields(103)g=srλf1′zy+g0z,which leads to the exact solution(104)ut,x,y,z=a0+a1f1zx+srλf1′zy+ψzEα,1λtα,wherea0anda1 are the arbitrary constants and f1zandψz are the two arbitrary functions.


#### 3.3.2. Exact Solutions under Riemann–Liouville Derivatives

We assume(105)u=a0tμ1+a1vx,y,ztμ2,leads to μ2−α=2μ2 and μ2=−α.Substituting (105) into (78) by comparing the coefficients of t−αandt−2α, we have(106a)a1pvyvxx+qvxvxy+rΩ0vy=0,Ω0=Γ1−αΓ1−2α,(106b)vxxxy−svxz=0.(I)
Note that (106b) is a linear equation, and we assume the summary type of a separate variable v=fx,z+gy,z which yields(107a)pa1fxx+rΩ0=0,(107b)fxz=0.Then, from (107b), we have(108)fx,z=f1x+f2z.Plugging (108) into (107a), we arrive at(109)f1x=−rΩ02pa1x2+c1x+c0.This gives rise to the exact solution(110)ut,x,y,z=a0tμ1+−rΓ1−α2pΓ1−2αx2+b1x+a1f2z+a1gy,z+b0t−α,α≠12,μ1≠α−1,whereμ1,a0,a1,b1=a1c1andb0=a1c0 are the arbitrary constants and f2zandgy,z are the two arbitrary functions.(II)
Noticing equation (106a), we assume that v=vξ,z,ξ=ax+by+c, then (106a) and (106b) become(111a)a1a2p+qvξξ+rΩ0=0,Ω0=Γ1−αΓ1−2α,(111b)vz=−ψzsa.Solving (111a) leads to(112)v=−rΩ02a1a2p+qξ2+ϕ1zξ+ϕ0z.Substituting (112) into (111b) yields(113)ϕ1′zξ+ϕ0′z=−ψzsa.By comparing the coefficients, we have(114)ϕ1z=c1,ϕ0z=−1sa∫zψz¯dz¯=−1saϕz,which gives rise to the exact solution(115)ut,x,y,z=a0tμ1+−rΓ1−α2a2p+qΓ1−2αax+by+c2+a1c1ax+by+c−a1saϕzt−α,α≠12,μ1≠α−1,wherea0,a1,c1,a,b,andc are the arbitrary constants and ϕz is an arbitrary function.Remark 10.
We do not consider the transformationξ=ax+by+cz+d for omitting the trivial solution.Remark 11.
The singularity only appears whenα=1/2andμ1=α−1 in solutions (110) and (115) in the Riemann–Liouville case.


#### 3.3.3. Dynamical Analysis of Exact Solutions for Fractional JM Equation

Under the Caputo case:(i)
As can be seen from Figure10, taking α=1/2,μ1=α−1,a0=1,a1=1,c0=1,c1=2,c2=1,c3=−3,a=2,r=b=s=c=1,λ=−1,d=−3y=1.5,z=2,t∈0,10,x∈−10,10, it is shown that solution (83) increases with increase in the spatial variable and decreases as time increases on all (u, t, x), (u, t, y), and (u, t, z)-planes.(ii)
As Figures11 and 12 show, we see that by choosing a=b=1,s=7/2,c=2,r=2,α=1/2,λ=−1/2,a0=c0=1,a1=−2,b0=1,c1=2,c2=1,c3=−1.5,y=2,z=2,t∈0,10,x∈−10,10, and x=2.5,y=2,z∈−10,10, the solutions (89) and (91) increase as x and y increases and decrease as time increases on both (u, t, x) and (u, t, y)-planes; however, ones have exponential decay as z⟶∞ when x and y are fixed.Under the Riemann–Liouville case:As Figure13 shows, taking α=2/3,μ1=±3/2,a0=2,p=r=1,b1=−3,a1fz+a1gy,z+b0=10,x∈−100,100,t∈0,100, we see that(i)
Ifμ1>0, solutions (110) and (115) will blow up as time-space tends to infinity.(ii)
Ifμ1<=0, solutions (110) and (115) increase as x increases if y and z are fixed and decay at maxμ1,α rate as t⟶∞ if the spatial variables are fixed.


## 3.1. HBP for the (2 + 1)-Dimensional Time-Fractional KP Equation

We first begin with the (2 + 1)-dimensional fractional KP equation:(13)∂tαu−6uux+uxxxx+3uyy=0,0<α≤1,which developed from the classical KdV equation, describes long dispersive wave propagating in two dimension in shallow water.


### 3.1.1. Exact Solutions under Riemann–Liouville Derivatives

According to (9), we suppose the exact solutions as(14)u=a0tμ1+a1vx,ytμ2,where a0,a1,μ1,andμ2 are the undetermined constants.Balancing the coefficients oft-power yields the following two subcases:Subcase I.(15)μ2−α=μ1+μ2=2μ2.Thus,μ1=μ2=−α.Substituting (15) into (14) by comparing the coefficients of t−αandt−2α yields(16a)Ω0vx=6a0+a1vvxx+6a1vx2,Ω0=Γ1−αΓ1−2α,(16b)vxxxx+3vyy=0.Then, integrating (16a) corresponding to variable x once by selecting an integrate constant 0 leads to(17)dvdx=Ω0v6a0+a1v.This gives rise to the implicit relation:(18)6a0lnv+a1v=Ω0x+ϕy,whereϕy is an undetermined function.Calculating the derivatives of implicit function (18) reads that(19a)vxxxx=a0Ω046a12v2−8a0a1v+a02v1296a17v7+7a0a16v6+21a02a15v5+35a03a14v4+35a04a13v3+21a05a12v2+7a06a1v+a07,(19b)vyy=Ω0Ω0a0ϕ′y2+6a0+a1v2ϕ′′yv2v36a13v3+3a0a12v2+3a02a1v+a03.Plugging (19a) and (19b) into (15), we have(20)a0=0,ϕ′′y=0,which leads to the solution of (15)(21)vx,y=Ω06a1x+b1y+b2,with arbitrary constantsa1,b1,andb2  and the exact solution of (11)(22)ut,x,y=a0+Γ1−α6Γ1−2αx+b1y+b2t−α,α≠12,wherec1andc2 are the arbitrary constants.Subcase II. μ2−α=2μ2andμ1+μ2=μ2.Thus,μ1=0andμ2=−α.Similarly as above, substituting (13) into (11) by comparing the coefficients of t−αandt−2α yields(23a)Ω0vx=6a1vvxx+6a1vx2,Ω0=Γ1−αΓ1−2α,(23b)vxxxx+6a0vxx+3vyy=0.Then, integrating (23a) corresponding to variable x by choosing an integrate constant 0 once leads to(24)Ω0v=6a1vvx.Sincev≠0 (only nontrivial solutions are considered), we have(25)dvdx=Ω06a1,which give rises to(26)v=Ω06a1x+ϕy.Plugging (26) into (23b) yields(27)ϕy=c1y+c2.By rewritingb1=a1c1andb2=a2c2, we have the exact solution(28)ut,x,y=a0+Γ1−α6Γ1−2αx+b1y+b2t−α,α≠12,wherea0,b1,andb2 are the arbitrary constants.Remark 2.
We have the same type solutions under the complex transformationξ=c1x+c2y+c0, so this case is omitted here.


### 3.1.2. Exact Solutions under Caputo Derivatives

According to (10), we assume the exact solutions formed as(29)u=a0+a1vx,yEα,1λtα.Substituting (29) into (11) by comparing the coefficients of Eα,1λtα leads to(30a)λa1vx+a1vxxxx+a0vxx+3a1vyy=0,(30b)vvxx=vx2.We consider the following two cases:Case I. By using the total differential rule, from (30b), we have(31)vvxx−vx2v2=vxvx=0.Integrating (31) leads to(32)vx=Cyv.Thus,(33)v=ϕyeCyx,whereϕyandCy are the undetermined functions.Plugging (33) into (30a) yields(34)ϕ′′y+cλ+c4+a0c23ϕy=0.Then, we have following three type solutions:(i)
δ=cλ+c4+a0c2<0.We obtainϕy=c1e−−δ/3y+c2e−δ/3y and the exact solution(35)ut,x,y=a0+c1ecx−−δ/3y+c2ecx+−δ/3yEα,1λtα,wherea0,c,c1,andc2 are the arbitrary constants.(ii)
δ=cλ+c4+a0c2>0.We obtain the exact solution(36)ut,x,y=a0+c1cosδ3y+c2sinδ3yecxEα,1λtα,wherea0,c,c1,andc2 are the arbitrary constants.(iii)
δ=cλ+c4+a0c2=0.We obtainϕy=c1y+c2 and the exact solution(37)ut,x,y=a0+c1y+c2ecxEα,1λtα,wherea0,c,c1,andc2 are the arbitrary constants.Case II. Using the complex differential transformation ξ=c1x+c2y+c0 on (30a) and (30b), we have(38a)c14v4ξ+a0c12+3c22v′′ξ+λc1v′ξ=0,(38b)vξv′′ξ=v′2ξ.Solving (38b) similarly as (25) yields(39)vξ=recξ.Taking (39) into (38a), we obtain c14c3+a0c12+3c22c+λc1=c0 and the exact solution by rewriting a=a1recc0,b1=cc1,b2=cc2(40)ut,x,y=a0+aeb1x+b2yEα,1λtα,wherea0,a,b1,andb2 are the arbitrary constants.Remark 3.
(33) contains product forms of separate variable v=fxgy, and for omitting the trivial solution, we do not consider the summary form of a separate variable v=fxgy.Remark 4.
The singularity only appears whenα=1/2 in solutions (22) and (28) in the Riemann–Liouville case.


### 3.1.3. Dynamical Analysis of Exact Solution for Fractional KP Model

Under the Riemann–Liouville case:As can be seen from Figure1, taking parameters α=1/3 and α=2,a0=1,b1=2,b2=1,y=1, and interval x∈−10,10,t∈0,10, it is shown that solution (28) increases with increase in the spatial variable and decreases with increase in time. When x and y are fixed, the solution tends to 0 at an α decay rate as t⟶∞, the larger the α is, the faster the u decay is.Under the Caputo case:(i)
As can be seen from Figure2, choosing c1=1,c2=−2,c=1,λ=−2,α=1/2,δ=−1,y=2,x∈−10,10,t∈0,10, and x=2,y∈−10,10,t∈0,10, it is shown that solution (35) decreases with increase in time-space since slope ux<0 or uy<0; however, it is unbounded as x or y tends to ∞.(ii)
As we see from Figure3, selecting c1=1,c2=−2,c=2,λ=12,α=1/2,δ=−1,y=2,x∈−10,10,t∈0,10, and x=2,y∈−30,10,t∈0,10, solution (36) decreases as time increases, periodic bounded as y, t tends to ∞ if x is fixed while unbounded as x, and t tends to ∞ if y is fixed.(iii)
As we see from Figure4, choosing c1=1,c2=−2,c=1,λ=−1,α=1/2,δ=0,y=3,x∈−10,10,t∈0,10, and x=2,y∈−10,10,t∈0,10, solution (37) decreases as time increases and unbounded as time-space tend to ∞.


## 3.1.1. Exact Solutions under Riemann–Liouville Derivatives

According to (9), we suppose the exact solutions as(14)u=a0tμ1+a1vx,ytμ2,where a0,a1,μ1,andμ2 are the undetermined constants.Balancing the coefficients oft-power yields the following two subcases:Subcase I.(15)μ2−α=μ1+μ2=2μ2.Thus,μ1=μ2=−α.Substituting (15) into (14) by comparing the coefficients of t−αandt−2α yields(16a)Ω0vx=6a0+a1vvxx+6a1vx2,Ω0=Γ1−αΓ1−2α,(16b)vxxxx+3vyy=0.Then, integrating (16a) corresponding to variable x once by selecting an integrate constant 0 leads to(17)dvdx=Ω0v6a0+a1v.This gives rise to the implicit relation:(18)6a0lnv+a1v=Ω0x+ϕy,whereϕy is an undetermined function.Calculating the derivatives of implicit function (18) reads that(19a)vxxxx=a0Ω046a12v2−8a0a1v+a02v1296a17v7+7a0a16v6+21a02a15v5+35a03a14v4+35a04a13v3+21a05a12v2+7a06a1v+a07,(19b)vyy=Ω0Ω0a0ϕ′y2+6a0+a1v2ϕ′′yv2v36a13v3+3a0a12v2+3a02a1v+a03.Plugging (19a) and (19b) into (15), we have(20)a0=0,ϕ′′y=0,which leads to the solution of (15)(21)vx,y=Ω06a1x+b1y+b2,with arbitrary constantsa1,b1,andb2  and the exact solution of (11)(22)ut,x,y=a0+Γ1−α6Γ1−2αx+b1y+b2t−α,α≠12,wherec1andc2 are the arbitrary constants.Subcase II. μ2−α=2μ2andμ1+μ2=μ2.Thus,μ1=0andμ2=−α.Similarly as above, substituting (13) into (11) by comparing the coefficients of t−αandt−2α yields(23a)Ω0vx=6a1vvxx+6a1vx2,Ω0=Γ1−αΓ1−2α,(23b)vxxxx+6a0vxx+3vyy=0.Then, integrating (23a) corresponding to variable x by choosing an integrate constant 0 once leads to(24)Ω0v=6a1vvx.Sincev≠0 (only nontrivial solutions are considered), we have(25)dvdx=Ω06a1,which give rises to(26)v=Ω06a1x+ϕy.Plugging (26) into (23b) yields(27)ϕy=c1y+c2.By rewritingb1=a1c1andb2=a2c2, we have the exact solution(28)ut,x,y=a0+Γ1−α6Γ1−2αx+b1y+b2t−α,α≠12,wherea0,b1,andb2 are the arbitrary constants.Remark 2.
We have the same type solutions under the complex transformationξ=c1x+c2y+c0, so this case is omitted here.


## 3.1.2. Exact Solutions under Caputo Derivatives

According to (10), we assume the exact solutions formed as(29)u=a0+a1vx,yEα,1λtα.Substituting (29) into (11) by comparing the coefficients of Eα,1λtα leads to(30a)λa1vx+a1vxxxx+a0vxx+3a1vyy=0,(30b)vvxx=vx2.We consider the following two cases:Case I. By using the total differential rule, from (30b), we have(31)vvxx−vx2v2=vxvx=0.Integrating (31) leads to(32)vx=Cyv.Thus,(33)v=ϕyeCyx,whereϕyandCy are the undetermined functions.Plugging (33) into (30a) yields(34)ϕ′′y+cλ+c4+a0c23ϕy=0.Then, we have following three type solutions:(i)
δ=cλ+c4+a0c2<0.We obtainϕy=c1e−−δ/3y+c2e−δ/3y and the exact solution(35)ut,x,y=a0+c1ecx−−δ/3y+c2ecx+−δ/3yEα,1λtα,wherea0,c,c1,andc2 are the arbitrary constants.(ii)
δ=cλ+c4+a0c2>0.We obtain the exact solution(36)ut,x,y=a0+c1cosδ3y+c2sinδ3yecxEα,1λtα,wherea0,c,c1,andc2 are the arbitrary constants.(iii)
δ=cλ+c4+a0c2=0.We obtainϕy=c1y+c2 and the exact solution(37)ut,x,y=a0+c1y+c2ecxEα,1λtα,wherea0,c,c1,andc2 are the arbitrary constants.Case II. Using the complex differential transformation ξ=c1x+c2y+c0 on (30a) and (30b), we have(38a)c14v4ξ+a0c12+3c22v′′ξ+λc1v′ξ=0,(38b)vξv′′ξ=v′2ξ.Solving (38b) similarly as (25) yields(39)vξ=recξ.Taking (39) into (38a), we obtain c14c3+a0c12+3c22c+λc1=c0 and the exact solution by rewriting a=a1recc0,b1=cc1,b2=cc2(40)ut,x,y=a0+aeb1x+b2yEα,1λtα,wherea0,a,b1,andb2 are the arbitrary constants.Remark 3.
(33) contains product forms of separate variable v=fxgy, and for omitting the trivial solution, we do not consider the summary form of a separate variable v=fxgy.Remark 4.
The singularity only appears whenα=1/2 in solutions (22) and (28) in the Riemann–Liouville case.


## 3.1.3. Dynamical Analysis of Exact Solution for Fractional KP Model

Under the Riemann–Liouville case:As can be seen from Figure1, taking parameters α=1/3 and α=2,a0=1,b1=2,b2=1,y=1, and interval x∈−10,10,t∈0,10, it is shown that solution (28) increases with increase in the spatial variable and decreases with increase in time. When x and y are fixed, the solution tends to 0 at an α decay rate as t⟶∞, the larger the α is, the faster the u decay is.Under the Caputo case:(i)
As can be seen from Figure2, choosing c1=1,c2=−2,c=1,λ=−2,α=1/2,δ=−1,y=2,x∈−10,10,t∈0,10, and x=2,y∈−10,10,t∈0,10, it is shown that solution (35) decreases with increase in time-space since slope ux<0 or uy<0; however, it is unbounded as x or y tends to ∞.(ii)
As we see from Figure3, selecting c1=1,c2=−2,c=2,λ=12,α=1/2,δ=−1,y=2,x∈−10,10,t∈0,10, and x=2,y∈−30,10,t∈0,10, solution (36) decreases as time increases, periodic bounded as y, t tends to ∞ if x is fixed while unbounded as x, and t tends to ∞ if y is fixed.(iii)
As we see from Figure4, choosing c1=1,c2=−2,c=1,λ=−1,α=1/2,δ=0,y=3,x∈−10,10,t∈0,10, and x=2,y∈−10,10,t∈0,10, solution (37) decreases as time increases and unbounded as time-space tend to ∞.


## 3.2. HBP for the (3 + 1)-Dimensional Fractional ZK Equation

### 3.2.1. Exact Solutions under Riemann–Liouville Derivatives

The (3 + 1)-dimensional fractional ZK equation is investigated in this section.(41)∂tαu+puux+quzzz+ruxxz+suyyz=0,0<α≤1,which represented the acoustic dynamics in a magnetized plasma in three-dimensional space with a low pressure.The exact solution can be assumed as(42)u=a0tμ1+a1vx,y,ztμ2.Balancing the coefficients oft-power by plugging (42) into (41) yields the following two cases:Cases I. μ1−α=μ1+μ2=2μ2That is,(43)μ1=μ2=−α.By comparing the coefficients oft−αandt−2α, we have(44a)Ω0+pa1vx=0,Ω0=Γ1−αΓ1−2α,(44b)qvzzz+rvxxz+svyyz=0.By solving (44a) directly, we obtain(45)vx,y,z=−Ω0pa1x+ϕy,z.Substituting (45) into (44b) yields(46)qϕzzz+sϕyyz=0.Subcase I. Assume η=b1y+b2z+b0, then we have(47)qb22+sb12ϕ′′′η=0.Solving (47) leads to the following two solutions:(i)
Whenqb22+sb12=0,sq<0, then(48)ux,y,z,t=a0−Γ1−αpΓ1−2αx+ϕb1y+b2z+b0t−α,α≠12,wherea0andb0 are the arbitrary constants and ϕ is an arbitrary function.(ii)
Whenϕ=c2η2/2+c1η+c0, then(49)ux,y,z,t=a−Γ1−αpΓ1−2αx+c2b1y+b2z+b022+c1b1y+b2zt−α,α≠12,wherea=a0+c0+c1b0,c1,andc2 are the arbitrary constants.Subcase II. Note that (46) is a linear equation, and if sq<0, by the linear PDE method, we have(50)ϕy,z=f1y−2−sqsf1y2−sz2−sq+f2−sqsy+z,which gives rise to exact solution whensq<0:(51)ux,y,z,t=a0−Γ1−αpΓ1−2αx+a1f1y−2−sqsf1y2−sz2−sq+f2−sqsy+zt−α,α≠12,wherea0anda1 are the arbitrary constants and f1andf2 are the two arbitrary functions.Subcase III. According to the linear PDE method, we assume the solution of (46) as a product form of a separate variable ϕy,z=fygz, thus(52)qg′′′z=c0g′z,sf′′y+c0fy=0,wherec0 is an arbitrary constant.By solving them directly, we obtain the following four solutions (67)–(70) where a0,a1,c1,c2,c3,c4,andc5 are the arbitrary constants and α≠1/2:(i)
ϕ=c1sinc0/sy+c2cosc0/syc3+c4sin−c0/qz+c5cos−c0/qz,(53)ux,y,z,t=a0−Γ1−αpΓ1−2αx+a1c1sinc0sy+c2cosc0syc3+c4sin−c0qz+c5cos−c0qzt−α.(ii)
ϕ=c1sinc0/sy+c2cosc0/syc3+c4ec0/qz+c5e−c0/qz,(54)ux,y,z,t=a0−Γ1−αpΓ1−2αx+a1c1sinc0sy+c2cosc0syc3+c4ec0/qz+c5e−c0/qzt−α.(iii)
ϕ=c1e−c0/sy+c2e−−c0/syc3+c4sin−c0/qz+c5cos−c0/qz,(55)ux,y,z,t=a0−Γ1−αpΓ1−2αx+a1c1e−c0/sy+c2e−−c0/syc3+c4sin−c0qz+c5cos−c0qzt−α.(iv)
ϕ=c1e−c0/qy+c2e−−c0/qyc3+c4ec0/qz+c5e−c0/qz,(56)ux,y,z,t=a0−Γ1−αpΓ1−2αx+a1c1e−c0/sy+c2e−−c0/syc3+c4ec0/qz+c5e−c0/qzt−α.Subcase IV. We can also assume the solution of (47) as a summary form of a separate variable ϕy,z=fy+gz, then(57)g′′′z=0.Thus,(58)gz=c22z22+c1z+c0,which leads to the exact solution by rewritingb2=a1c22,b1=a1c1anda=a0+a1c0:(59)ux,y,z,t=a−Γ1−αpΓ1−2αx+a1fy+b2z22+b1zt−α,α≠12,wherea,a1,b2,andb1 are the arbitrary constants and fy is an arbitrary function.Cases II. μ2−α=2μ2andμ1−α=μ1+μ2.By balancing equation (41), we have(60)μ1=0,μ2=−α.Comparing the coefficients oft−αandt−2α yields(61a)a0ω0+pa0a1vx+qa1vzzz+ra1vxxz+sa1vyyz=0,ω0=1Γ1−α,(61b)Ω0v+pa1vvx=0,Ω0=Γ1−αΓ1−2α.Solving (61b) leads to the following form:(62)v=−Ω0pa1x+ϕy,z.Substituting (62) into (61a) leads to(63)qa1ϕzzz+sa1ϕyyz=a0Ω0−ω0.(i)
Considering the complex transformationη=b1y+b2z+b0, we obtain(64)a1b2qb22+sb12ϕ′′′η=a0Ω0−ω0.Solving (81a) and (81b) leads to(65)ϕη=a0Ω0−ω06a1b2qb22+sb12η3+c22η2+c1η+c0.Then, we obtain the exact solution(66)ut,x,y,z=a0+−Γ1−αpΓ1−2αx+a0Γ1−α/Γ1−2α−1/Γ1−α6b2qb22+sb12b1y+b2z+b03+a1c22b1y+b2z+b02+a1c1b1y+b2z+b0+a1c0t−α,α≠12,wherea0,b1,b2,a1,c0,c1,andc2 are the arbitrary constants.(ii)
Note that (63) is a linear equation, and by the linear PDE method, we have(67)ϕy,z=f1y+12f2−sq/sy+zsqa1−24f1y/2−sz/2−sq−sqqa1+a0Ω0−ω0sz3+3qa0Ω0−ω0zy212sqa1,which leads to the exact solution(68)ut,x,y,z=a0+−Ω0px+a1f1y+12f2−sq/sy+zsqa1−24f1y/2−sz/2−sq−sqqa1+a0Ω0−ω0sz3+3qa0Ω0−ω0zy212sqt−α,α≠12,wherea0anda1 are constants and f1andf2 are the two arbitrary functions.(iii)
Using summary form of the separate variableϕ=fy+gz yields(69)qa1g′′′z=a0Ω0−ω0.Solving (69) directly, we arrive at(70)gz=a0Ω0−ω06qa1z3+b22z2+b1z+b0,which gives rise to the exact solution(71)ut,x,y,z=a0+−Γ1−αpΓ1−2αx+a1fy+a06qΓ1−αΓ1−2α−1Γ1−αz3+a1b22z2+a1b1z+a1b0t−α,α≠12,wherea0,a1,b2,b1,andb0 are the arbitrary constants and fy is an arbitrary function.(iv)
Whena0=0, using the product form of the separate variable ϕ=fygz yields(72)qfyg′′′z+sf′′′ygz=0,which give rises to(73)f′′y=c0fy,qg′′′z+sc0gz=0.From (73), we have(74)gz=c1e−sc0q23/qz+e−sc0q23/2qzc2cos3−sc0q232qz+c3sin3−sc0q232qz,and following three solutions (75)–(77) where a1,b1,b2,c1,c2,andc3 are constants and α≠1/2:(i)
Ifc0>0, then(75)ut,x,y,z=−Γ1−αpa1Γ1−2αx+b1e−c0y+b2ec0yc1e−sc0q23/qz+e−sc0q23/2qzc2cos3−sc0q232qz+c3sin3−sc0q232qzt−α.(ii)
Ifc0<0, then(76)ut,x,y,z=−Γ1−αpa1Γ1−2αx+b1cos−c0y+b2sin−c0yc1e−sc0q23/qz+e−sc0q23/2qzc2cos3−sc0q232qz+c3sin3−sc0q232qzt−α.(iii)
Ifc0=0, then(77)ut,x,y,z=−Γ1−αpa1Γ1−2αx+b1y+b2c1z2+c2z+c3t−α.Remark 5.
In case I, we do not concern the complex transformξ=ax+by+cz+d for (48) contains solutions of this type.Remark 6.
In case II, we do not consider the complex transformξ=ax+by+cz+d for omitting the trivial solution.Remark 7.
The singularity only appears whenα=1/2 in solutions (53)–(56), (66), (68), (71), and (76)–(78) in the Riemann–Liouville case.Remark 8.
If taking form (10) into the fractional ZK model (41), we see that the nonlinear term uux will lead to trivial solutions, thus the Caputo derivative is not considered.


### 3.2.2. Dynamical Analysis of Exact Solution for Fractional ZK Model

(i)
As can be seen from Figure5, taking α=3/5,a=3,a0=c0=1,c1=1/2,b0=2,b1=1,b2=−2,c2=2,t∈0,10,y=1.5,z=2.5,x∈−10,10, x=−2.5,z=2.5,y∈−10,10, and x=−2.5,y=1.5,z∈−10,10, it is shown that solution (49) increases with increase in spatial variables and decreases as time increases. When x, y, and z are fixed, solution (49) tends to 0 at α decay rate as t⟶∞, the larger the α is, the faster the u decay is.
(ii)
From Figure6, by taking α=2/5,a0=1/2,a1=2,c1=1.5,c2=−2,c3=−3,c4=−2.5,c5=1,c0=2,p=3,s=q=1/2,t∈0,30,x=−2.5,z=2.5,y∈−30,30 and t∈0,10,andx=−2.5,y=1.5,z∈−10,10, we see solution (54) is periodic bounded and tends to 0 at α decay rate as y and t tend to ∞ if x and z are fixed. The larger the α is, the faster the u decay is. When x and y are fixed, solution (54) decreases to −∞ as z increases to +∞.
(iii)
As can be seen from Figure7, selecting α=2/3,p=3,a0=2,a1=1,b1=2,b2=−1,b0=3,c0=1,c1=1/2,c2=2,q=s=1/2,t∈0,10,y=1.5,z=2.5,x∈−10,10 , x=−2.5,z=2.5,y∈−10,10, and x=−2.5,y=1.5,z∈−10,10, while from Figure 8 by choosing α=2/3,p=3,a1=1,b1=1,b2=−1.5,c1=1.5,c2=2,c3=−2.5,c0=2,s=8,q=2,t∈0,10,y=1.5,z=2.5,x∈−10,10, x=−2.5,z=2.5,y∈−10,10, and x=−2.5,y=1.5,z∈−10,10, it is shown that solutions (66) and (75) increase to ∞ or decrease to −∞ as x, y, or z increases to ∞, and ones also tends to 0 at α decay rate as t⟶∞ if x, y, and z are fixed. The larger the α is, the faster the u decay is.
(iv)
From Figure9, as we choose α=2/3,p=1,a1=3,b1=1,b2=−1.5,c0=2,c1=1,c2=2,c3=−2.5,s=8,q=2,x=z=−2.5,y∈−30,30,t∈0,30, the solution (76) is periodic bounded and has α decay rate as t⟶∞ if x and z are fixed, and one has the same dynamical properties as (54) if else.


## 3.2.1. Exact Solutions under Riemann–Liouville Derivatives

The (3 + 1)-dimensional fractional ZK equation is investigated in this section.(41)∂tαu+puux+quzzz+ruxxz+suyyz=0,0<α≤1,which represented the acoustic dynamics in a magnetized plasma in three-dimensional space with a low pressure.The exact solution can be assumed as(42)u=a0tμ1+a1vx,y,ztμ2.Balancing the coefficients oft-power by plugging (42) into (41) yields the following two cases:Cases I. μ1−α=μ1+μ2=2μ2That is,(43)μ1=μ2=−α.By comparing the coefficients oft−αandt−2α, we have(44a)Ω0+pa1vx=0,Ω0=Γ1−αΓ1−2α,(44b)qvzzz+rvxxz+svyyz=0.By solving (44a) directly, we obtain(45)vx,y,z=−Ω0pa1x+ϕy,z.Substituting (45) into (44b) yields(46)qϕzzz+sϕyyz=0.Subcase I. Assume η=b1y+b2z+b0, then we have(47)qb22+sb12ϕ′′′η=0.Solving (47) leads to the following two solutions:(i)
Whenqb22+sb12=0,sq<0, then(48)ux,y,z,t=a0−Γ1−αpΓ1−2αx+ϕb1y+b2z+b0t−α,α≠12,wherea0andb0 are the arbitrary constants and ϕ is an arbitrary function.(ii)
Whenϕ=c2η2/2+c1η+c0, then(49)ux,y,z,t=a−Γ1−αpΓ1−2αx+c2b1y+b2z+b022+c1b1y+b2zt−α,α≠12,wherea=a0+c0+c1b0,c1,andc2 are the arbitrary constants.Subcase II. Note that (46) is a linear equation, and if sq<0, by the linear PDE method, we have(50)ϕy,z=f1y−2−sqsf1y2−sz2−sq+f2−sqsy+z,which gives rise to exact solution whensq<0:(51)ux,y,z,t=a0−Γ1−αpΓ1−2αx+a1f1y−2−sqsf1y2−sz2−sq+f2−sqsy+zt−α,α≠12,wherea0anda1 are the arbitrary constants and f1andf2 are the two arbitrary functions.Subcase III. According to the linear PDE method, we assume the solution of (46) as a product form of a separate variable ϕy,z=fygz, thus(52)qg′′′z=c0g′z,sf′′y+c0fy=0,wherec0 is an arbitrary constant.By solving them directly, we obtain the following four solutions (67)–(70) where a0,a1,c1,c2,c3,c4,andc5 are the arbitrary constants and α≠1/2:(i)
ϕ=c1sinc0/sy+c2cosc0/syc3+c4sin−c0/qz+c5cos−c0/qz,(53)ux,y,z,t=a0−Γ1−αpΓ1−2αx+a1c1sinc0sy+c2cosc0syc3+c4sin−c0qz+c5cos−c0qzt−α.(ii)
ϕ=c1sinc0/sy+c2cosc0/syc3+c4ec0/qz+c5e−c0/qz,(54)ux,y,z,t=a0−Γ1−αpΓ1−2αx+a1c1sinc0sy+c2cosc0syc3+c4ec0/qz+c5e−c0/qzt−α.(iii)
ϕ=c1e−c0/sy+c2e−−c0/syc3+c4sin−c0/qz+c5cos−c0/qz,(55)ux,y,z,t=a0−Γ1−αpΓ1−2αx+a1c1e−c0/sy+c2e−−c0/syc3+c4sin−c0qz+c5cos−c0qzt−α.(iv)
ϕ=c1e−c0/qy+c2e−−c0/qyc3+c4ec0/qz+c5e−c0/qz,(56)ux,y,z,t=a0−Γ1−αpΓ1−2αx+a1c1e−c0/sy+c2e−−c0/syc3+c4ec0/qz+c5e−c0/qzt−α.Subcase IV. We can also assume the solution of (47) as a summary form of a separate variable ϕy,z=fy+gz, then(57)g′′′z=0.Thus,(58)gz=c22z22+c1z+c0,which leads to the exact solution by rewritingb2=a1c22,b1=a1c1anda=a0+a1c0:(59)ux,y,z,t=a−Γ1−αpΓ1−2αx+a1fy+b2z22+b1zt−α,α≠12,wherea,a1,b2,andb1 are the arbitrary constants and fy is an arbitrary function.Cases II. μ2−α=2μ2andμ1−α=μ1+μ2.By balancing equation (41), we have(60)μ1=0,μ2=−α.Comparing the coefficients oft−αandt−2α yields(61a)a0ω0+pa0a1vx+qa1vzzz+ra1vxxz+sa1vyyz=0,ω0=1Γ1−α,(61b)Ω0v+pa1vvx=0,Ω0=Γ1−αΓ1−2α.Solving (61b) leads to the following form:(62)v=−Ω0pa1x+ϕy,z.Substituting (62) into (61a) leads to(63)qa1ϕzzz+sa1ϕyyz=a0Ω0−ω0.(i)
Considering the complex transformationη=b1y+b2z+b0, we obtain(64)a1b2qb22+sb12ϕ′′′η=a0Ω0−ω0.Solving (81a) and (81b) leads to(65)ϕη=a0Ω0−ω06a1b2qb22+sb12η3+c22η2+c1η+c0.Then, we obtain the exact solution(66)ut,x,y,z=a0+−Γ1−αpΓ1−2αx+a0Γ1−α/Γ1−2α−1/Γ1−α6b2qb22+sb12b1y+b2z+b03+a1c22b1y+b2z+b02+a1c1b1y+b2z+b0+a1c0t−α,α≠12,wherea0,b1,b2,a1,c0,c1,andc2 are the arbitrary constants.(ii)
Note that (63) is a linear equation, and by the linear PDE method, we have(67)ϕy,z=f1y+12f2−sq/sy+zsqa1−24f1y/2−sz/2−sq−sqqa1+a0Ω0−ω0sz3+3qa0Ω0−ω0zy212sqa1,which leads to the exact solution(68)ut,x,y,z=a0+−Ω0px+a1f1y+12f2−sq/sy+zsqa1−24f1y/2−sz/2−sq−sqqa1+a0Ω0−ω0sz3+3qa0Ω0−ω0zy212sqt−α,α≠12,wherea0anda1 are constants and f1andf2 are the two arbitrary functions.(iii)
Using summary form of the separate variableϕ=fy+gz yields(69)qa1g′′′z=a0Ω0−ω0.Solving (69) directly, we arrive at(70)gz=a0Ω0−ω06qa1z3+b22z2+b1z+b0,which gives rise to the exact solution(71)ut,x,y,z=a0+−Γ1−αpΓ1−2αx+a1fy+a06qΓ1−αΓ1−2α−1Γ1−αz3+a1b22z2+a1b1z+a1b0t−α,α≠12,wherea0,a1,b2,b1,andb0 are the arbitrary constants and fy is an arbitrary function.(iv)
Whena0=0, using the product form of the separate variable ϕ=fygz yields(72)qfyg′′′z+sf′′′ygz=0,which give rises to(73)f′′y=c0fy,qg′′′z+sc0gz=0.From (73), we have(74)gz=c1e−sc0q23/qz+e−sc0q23/2qzc2cos3−sc0q232qz+c3sin3−sc0q232qz,and following three solutions (75)–(77) where a1,b1,b2,c1,c2,andc3 are constants and α≠1/2:(i)
Ifc0>0, then(75)ut,x,y,z=−Γ1−αpa1Γ1−2αx+b1e−c0y+b2ec0yc1e−sc0q23/qz+e−sc0q23/2qzc2cos3−sc0q232qz+c3sin3−sc0q232qzt−α.(ii)
Ifc0<0, then(76)ut,x,y,z=−Γ1−αpa1Γ1−2αx+b1cos−c0y+b2sin−c0yc1e−sc0q23/qz+e−sc0q23/2qzc2cos3−sc0q232qz+c3sin3−sc0q232qzt−α.(iii)
Ifc0=0, then(77)ut,x,y,z=−Γ1−αpa1Γ1−2αx+b1y+b2c1z2+c2z+c3t−α.Remark 5.
In case I, we do not concern the complex transformξ=ax+by+cz+d for (48) contains solutions of this type.Remark 6.
In case II, we do not consider the complex transformξ=ax+by+cz+d for omitting the trivial solution.Remark 7.
The singularity only appears whenα=1/2 in solutions (53)–(56), (66), (68), (71), and (76)–(78) in the Riemann–Liouville case.Remark 8.
If taking form (10) into the fractional ZK model (41), we see that the nonlinear term uux will lead to trivial solutions, thus the Caputo derivative is not considered.


## 3.2.2. Dynamical Analysis of Exact Solution for Fractional ZK Model

(i)
As can be seen from Figure5, taking α=3/5,a=3,a0=c0=1,c1=1/2,b0=2,b1=1,b2=−2,c2=2,t∈0,10,y=1.5,z=2.5,x∈−10,10, x=−2.5,z=2.5,y∈−10,10, and x=−2.5,y=1.5,z∈−10,10, it is shown that solution (49) increases with increase in spatial variables and decreases as time increases. When x, y, and z are fixed, solution (49) tends to 0 at α decay rate as t⟶∞, the larger the α is, the faster the u decay is.
(ii)
From Figure6, by taking α=2/5,a0=1/2,a1=2,c1=1.5,c2=−2,c3=−3,c4=−2.5,c5=1,c0=2,p=3,s=q=1/2,t∈0,30,x=−2.5,z=2.5,y∈−30,30 and t∈0,10,andx=−2.5,y=1.5,z∈−10,10, we see solution (54) is periodic bounded and tends to 0 at α decay rate as y and t tend to ∞ if x and z are fixed. The larger the α is, the faster the u decay is. When x and y are fixed, solution (54) decreases to −∞ as z increases to +∞.
(iii)
As can be seen from Figure7, selecting α=2/3,p=3,a0=2,a1=1,b1=2,b2=−1,b0=3,c0=1,c1=1/2,c2=2,q=s=1/2,t∈0,10,y=1.5,z=2.5,x∈−10,10 , x=−2.5,z=2.5,y∈−10,10, and x=−2.5,y=1.5,z∈−10,10, while from Figure 8 by choosing α=2/3,p=3,a1=1,b1=1,b2=−1.5,c1=1.5,c2=2,c3=−2.5,c0=2,s=8,q=2,t∈0,10,y=1.5,z=2.5,x∈−10,10, x=−2.5,z=2.5,y∈−10,10, and x=−2.5,y=1.5,z∈−10,10, it is shown that solutions (66) and (75) increase to ∞ or decrease to −∞ as x, y, or z increases to ∞, and ones also tends to 0 at α decay rate as t⟶∞ if x, y, and z are fixed. The larger the α is, the faster the u decay is.
(iv)
From Figure9, as we choose α=2/3,p=1,a1=3,b1=1,b2=−1.5,c0=2,c1=1,c2=2,c3=−2.5,s=8,q=2,x=z=−2.5,y∈−30,30,t∈0,30, the solution (76) is periodic bounded and has α decay rate as t⟶∞ if x and z are fixed, and one has the same dynamical properties as (54) if else.


## 3.3. HBP for the (3 + 1)-Dimensional Fractional JM Equation

The (3 + 1)-dimensional fractional JM equation will be researched in the following contents:(78)uxxxy+puyuxx+quxuxy+r∂tαuy−suxz=0,0<α≤1,which is developed from the second members of integrable systems of the classical KP hierarchy.


### 3.3.1. Exact Solutions under Caputo Derivatives

Assume(79)u=a0+a1vx,y,zEα,1λtα.Injecting (79) into (78), we arrive at(80a)vxxxy+rλvy−svxz=0,(80b)pvyvxx+qvxvxy=0.(i)
We assume thatv=vξ,ξ=ax+by+cz+d, then (80a) and (80b) becomes the following ODEs:(81a)a3bv4ξ+rλbv′ξ−sacv′′ξ=0,(81b)p+qv′ξv′′ξ=0.Whenq=−p, (81b) is a free equation, and then solving (81a), we have following solutions by noting Δ=327r2λ2b3−4s3c3/b−9rλb3b2,(82)vξ=c0bλr+c1e121/3Δ2+122/3scb/6abΔξ+e−121/3Δ2+122/3scb/12abΔξc2cos3121/3Δ2−122/3scb12abΔξ+c3sin3121/3Δ2−122/3scb12abΔξ,which leads to the exact solution:(83)ut,x,y,z=a0+a1c0bλr+c1e121/3Δ2+122/3scb/6abΔax+by+cz+d+e−121/3Δ2+122/3scb/12abΔax+by+cz+dc2cos3121/3Δ2−122/3scb12abΔax+by+cz+d+c3sin3121/3Δ2−122/3scb12abΔax+by+cz+dEα,1λtα,wherea0,a1,c0,a,b,c,andd are the arbitrary constants.Remark 9.
In case ofq≠−p, from v′′ξ=0, we obtain v′ξ=0 which only leads to a trivial solution, thus we do not consider this case.(ii)
Assumevx,y,z=vη,z,η=ax+by+c, and from (80a) and (80b), we have(84a)a3bvηηηη+rλbvη−savηz=0,(84b)p+qvηvηη=0.Then (84a) and (84b) yields the following two subcases(1)
Ifq=−p, we only consider (84a).(i)
Note that it is a linear equation, and we supposevη,z=fηgz, which arrives at(85)a3bf4ηgz+rλbf′ηgz−saf′ηg′z=0.Using the separate variable method, we have(86a)g′z=cgz,(86b)a3bf′′′η+rλb−sacfη=c0.Solving (86a) reads that(87)gz=b0ecz,and the solution of (86b) and the exact solution of (78) is described as the following two cases:Ifrλb≠sac, then(88)fη=c0rbλ−sac+c1esac−rbλb23/abη+e−sac−rbλb2/2ab3ηc2cos3sac−rbλb232abη+c3sin3sac−rbλb232abη,and the exact solution is(89)ut,x,y,z=a0+a1b0c0rbλ−sac+c1esac−rbλb23/abax+by+c+e−sac−rbλb23/2abax+by+cc2cos3sac−rbλb232abx+by+ca+c3sin3sac−rbλb232abax+by+ceczEα,1λtα,wherea0,a1,b0,c0,a,b,c,c1,c2,andc3 are the arbitrary constants.Ifc=rbλ/sa, then(90)fη=c06η3+c22η2+c1η+c3,and the exact solution is(91)ut,x,y,z=a0+ac06ax+by+rbλsa3+c22ax+by+rbλsa2+c1ax+by+rbλsa+c3erbλ/sazEα,1λtα,wherea0,c0,a,b,c1,c2,andc3 are the arbitrary constants.(ii)
Solving (84a) directly by the linear PDE method reads that(92)vη,z=c−1/3c1ec1/3η+ba2c/s+brλ/saz+4e−1/2c1/3η+ba2c/s+brλ/sazc2sin3c1/32η+c3cos3c1/32η+fz,ut,x,y,z=a0+a1c−1/3c1ec1/3ax+by+c+ba2c/s+brλ/saz+4e−1/2c1/3ax+by+c+ba2c/s+brλ/sazc2sin3c1/32ax+by+c+c3cos3c1/32ax+by+c+fzEα,1λtα,wherea0,a1,a,b,c,c1,c2,andc3 are the arbitrary constants and fz is an arbitrary function.(iii)
Assumevη,z=fη+gz, and from (84a), we have(93)a3f4η+rλf′η=0,and then(94)a3f′′′η+rλfη=c0.By solving (94), we have(95)fη=c0rλ+c1e−rλ3/aη+erλ3/2aηc2cos3rλ32aη+c3sin3rλ32aη,and the exact solution(96)ut,x,y,z=a0+a1c0rλ+c1e−rλ3/aax+by+c+erλ3/2aax+by+cc2cos3rλ32aax+by+c+c3sin3rλ32aax+by+c+gzEα,1λtα,wherea0,a1,a,b,c,c1,c2,andc3 are the arbitrary constants and gz is an arbitrary function.(2)
Ifq≠−p, then for the nontrivial case, we have vηη=0, and this leads to(97)v=f1zη+f2z.Then, substituting (97) into (84a), we have(98)rλbf1z−saf1′z=0,which yields(99)f1z=c0erλb/saz.This gives rise to the exact solution(100)ut,x,y,z=a0+a1c0erλb/sazax+by+c+f2zEα,1λtα,wherea0,a1,c0,a,b,andc are the arbitrary constants and f2z is an arbitrary function.(III)
Assumev=fx,z+gy,z.Then, from (84a) and (84b), we arrive at(101a)gyfxx=0,(101b)rλgy−sfxz=0.For the nontrivial case, we obtain(102)fxx=0⇒f=f1zx+f2z.Injecting (102) into (101b) yields(103)g=srλf1′zy+g0z,which leads to the exact solution(104)ut,x,y,z=a0+a1f1zx+srλf1′zy+ψzEα,1λtα,wherea0anda1 are the arbitrary constants and f1zandψz are the two arbitrary functions.


### 3.3.2. Exact Solutions under Riemann–Liouville Derivatives

We assume(105)u=a0tμ1+a1vx,y,ztμ2,leads to μ2−α=2μ2 and μ2=−α.Substituting (105) into (78) by comparing the coefficients of t−αandt−2α, we have(106a)a1pvyvxx+qvxvxy+rΩ0vy=0,Ω0=Γ1−αΓ1−2α,(106b)vxxxy−svxz=0.(I)
Note that (106b) is a linear equation, and we assume the summary type of a separate variable v=fx,z+gy,z which yields(107a)pa1fxx+rΩ0=0,(107b)fxz=0.Then, from (107b), we have(108)fx,z=f1x+f2z.Plugging (108) into (107a), we arrive at(109)f1x=−rΩ02pa1x2+c1x+c0.This gives rise to the exact solution(110)ut,x,y,z=a0tμ1+−rΓ1−α2pΓ1−2αx2+b1x+a1f2z+a1gy,z+b0t−α,α≠12,μ1≠α−1,whereμ1,a0,a1,b1=a1c1andb0=a1c0 are the arbitrary constants and f2zandgy,z are the two arbitrary functions.(II)
Noticing equation (106a), we assume that v=vξ,z,ξ=ax+by+c, then (106a) and (106b) become(111a)a1a2p+qvξξ+rΩ0=0,Ω0=Γ1−αΓ1−2α,(111b)vz=−ψzsa.Solving (111a) leads to(112)v=−rΩ02a1a2p+qξ2+ϕ1zξ+ϕ0z.Substituting (112) into (111b) yields(113)ϕ1′zξ+ϕ0′z=−ψzsa.By comparing the coefficients, we have(114)ϕ1z=c1,ϕ0z=−1sa∫zψz¯dz¯=−1saϕz,which gives rise to the exact solution(115)ut,x,y,z=a0tμ1+−rΓ1−α2a2p+qΓ1−2αax+by+c2+a1c1ax+by+c−a1saϕzt−α,α≠12,μ1≠α−1,wherea0,a1,c1,a,b,andc are the arbitrary constants and ϕz is an arbitrary function.Remark 10.
We do not consider the transformationξ=ax+by+cz+d for omitting the trivial solution.Remark 11.
The singularity only appears whenα=1/2andμ1=α−1 in solutions (110) and (115) in the Riemann–Liouville case.


### 3.3.3. Dynamical Analysis of Exact Solutions for Fractional JM Equation

Under the Caputo case:(i)
As can be seen from Figure10, taking α=1/2,μ1=α−1,a0=1,a1=1,c0=1,c1=2,c2=1,c3=−3,a=2,r=b=s=c=1,λ=−1,d=−3y=1.5,z=2,t∈0,10,x∈−10,10, it is shown that solution (83) increases with increase in the spatial variable and decreases as time increases on all (u, t, x), (u, t, y), and (u, t, z)-planes.(ii)
As Figures11 and 12 show, we see that by choosing a=b=1,s=7/2,c=2,r=2,α=1/2,λ=−1/2,a0=c0=1,a1=−2,b0=1,c1=2,c2=1,c3=−1.5,y=2,z=2,t∈0,10,x∈−10,10, and x=2.5,y=2,z∈−10,10, the solutions (89) and (91) increase as x and y increases and decrease as time increases on both (u, t, x) and (u, t, y)-planes; however, ones have exponential decay as z⟶∞ when x and y are fixed.Under the Riemann–Liouville case:As Figure13 shows, taking α=2/3,μ1=±3/2,a0=2,p=r=1,b1=−3,a1fz+a1gy,z+b0=10,x∈−100,100,t∈0,100, we see that(i)
Ifμ1>0, solutions (110) and (115) will blow up as time-space tends to infinity.(ii)
Ifμ1<=0, solutions (110) and (115) increase as x increases if y and z are fixed and decay at maxμ1,α rate as t⟶∞ if the spatial variables are fixed.


## 3.3.1. Exact Solutions under Caputo Derivatives

Assume(79)u=a0+a1vx,y,zEα,1λtα.Injecting (79) into (78), we arrive at(80a)vxxxy+rλvy−svxz=0,(80b)pvyvxx+qvxvxy=0.(i)
We assume thatv=vξ,ξ=ax+by+cz+d, then (80a) and (80b) becomes the following ODEs:(81a)a3bv4ξ+rλbv′ξ−sacv′′ξ=0,(81b)p+qv′ξv′′ξ=0.Whenq=−p, (81b) is a free equation, and then solving (81a), we have following solutions by noting Δ=327r2λ2b3−4s3c3/b−9rλb3b2,(82)vξ=c0bλr+c1e121/3Δ2+122/3scb/6abΔξ+e−121/3Δ2+122/3scb/12abΔξc2cos3121/3Δ2−122/3scb12abΔξ+c3sin3121/3Δ2−122/3scb12abΔξ,which leads to the exact solution:(83)ut,x,y,z=a0+a1c0bλr+c1e121/3Δ2+122/3scb/6abΔax+by+cz+d+e−121/3Δ2+122/3scb/12abΔax+by+cz+dc2cos3121/3Δ2−122/3scb12abΔax+by+cz+d+c3sin3121/3Δ2−122/3scb12abΔax+by+cz+dEα,1λtα,wherea0,a1,c0,a,b,c,andd are the arbitrary constants.Remark 9.
In case ofq≠−p, from v′′ξ=0, we obtain v′ξ=0 which only leads to a trivial solution, thus we do not consider this case.(ii)
Assumevx,y,z=vη,z,η=ax+by+c, and from (80a) and (80b), we have(84a)a3bvηηηη+rλbvη−savηz=0,(84b)p+qvηvηη=0.Then (84a) and (84b) yields the following two subcases(1)
Ifq=−p, we only consider (84a).(i)
Note that it is a linear equation, and we supposevη,z=fηgz, which arrives at(85)a3bf4ηgz+rλbf′ηgz−saf′ηg′z=0.Using the separate variable method, we have(86a)g′z=cgz,(86b)a3bf′′′η+rλb−sacfη=c0.Solving (86a) reads that(87)gz=b0ecz,and the solution of (86b) and the exact solution of (78) is described as the following two cases:Ifrλb≠sac, then(88)fη=c0rbλ−sac+c1esac−rbλb23/abη+e−sac−rbλb2/2ab3ηc2cos3sac−rbλb232abη+c3sin3sac−rbλb232abη,and the exact solution is(89)ut,x,y,z=a0+a1b0c0rbλ−sac+c1esac−rbλb23/abax+by+c+e−sac−rbλb23/2abax+by+cc2cos3sac−rbλb232abx+by+ca+c3sin3sac−rbλb232abax+by+ceczEα,1λtα,wherea0,a1,b0,c0,a,b,c,c1,c2,andc3 are the arbitrary constants.Ifc=rbλ/sa, then(90)fη=c06η3+c22η2+c1η+c3,and the exact solution is(91)ut,x,y,z=a0+ac06ax+by+rbλsa3+c22ax+by+rbλsa2+c1ax+by+rbλsa+c3erbλ/sazEα,1λtα,wherea0,c0,a,b,c1,c2,andc3 are the arbitrary constants.(ii)
Solving (84a) directly by the linear PDE method reads that(92)vη,z=c−1/3c1ec1/3η+ba2c/s+brλ/saz+4e−1/2c1/3η+ba2c/s+brλ/sazc2sin3c1/32η+c3cos3c1/32η+fz,ut,x,y,z=a0+a1c−1/3c1ec1/3ax+by+c+ba2c/s+brλ/saz+4e−1/2c1/3ax+by+c+ba2c/s+brλ/sazc2sin3c1/32ax+by+c+c3cos3c1/32ax+by+c+fzEα,1λtα,wherea0,a1,a,b,c,c1,c2,andc3 are the arbitrary constants and fz is an arbitrary function.(iii)
Assumevη,z=fη+gz, and from (84a), we have(93)a3f4η+rλf′η=0,and then(94)a3f′′′η+rλfη=c0.By solving (94), we have(95)fη=c0rλ+c1e−rλ3/aη+erλ3/2aηc2cos3rλ32aη+c3sin3rλ32aη,and the exact solution(96)ut,x,y,z=a0+a1c0rλ+c1e−rλ3/aax+by+c+erλ3/2aax+by+cc2cos3rλ32aax+by+c+c3sin3rλ32aax+by+c+gzEα,1λtα,wherea0,a1,a,b,c,c1,c2,andc3 are the arbitrary constants and gz is an arbitrary function.(2)
Ifq≠−p, then for the nontrivial case, we have vηη=0, and this leads to(97)v=f1zη+f2z.Then, substituting (97) into (84a), we have(98)rλbf1z−saf1′z=0,which yields(99)f1z=c0erλb/saz.This gives rise to the exact solution(100)ut,x,y,z=a0+a1c0erλb/sazax+by+c+f2zEα,1λtα,wherea0,a1,c0,a,b,andc are the arbitrary constants and f2z is an arbitrary function.(III)
Assumev=fx,z+gy,z.Then, from (84a) and (84b), we arrive at(101a)gyfxx=0,(101b)rλgy−sfxz=0.For the nontrivial case, we obtain(102)fxx=0⇒f=f1zx+f2z.Injecting (102) into (101b) yields(103)g=srλf1′zy+g0z,which leads to the exact solution(104)ut,x,y,z=a0+a1f1zx+srλf1′zy+ψzEα,1λtα,wherea0anda1 are the arbitrary constants and f1zandψz are the two arbitrary functions.


## 3.3.2. Exact Solutions under Riemann–Liouville Derivatives

We assume(105)u=a0tμ1+a1vx,y,ztμ2,leads to μ2−α=2μ2 and μ2=−α.Substituting (105) into (78) by comparing the coefficients of t−αandt−2α, we have(106a)a1pvyvxx+qvxvxy+rΩ0vy=0,Ω0=Γ1−αΓ1−2α,(106b)vxxxy−svxz=0.(I)
Note that (106b) is a linear equation, and we assume the summary type of a separate variable v=fx,z+gy,z which yields(107a)pa1fxx+rΩ0=0,(107b)fxz=0.Then, from (107b), we have(108)fx,z=f1x+f2z.Plugging (108) into (107a), we arrive at(109)f1x=−rΩ02pa1x2+c1x+c0.This gives rise to the exact solution(110)ut,x,y,z=a0tμ1+−rΓ1−α2pΓ1−2αx2+b1x+a1f2z+a1gy,z+b0t−α,α≠12,μ1≠α−1,whereμ1,a0,a1,b1=a1c1andb0=a1c0 are the arbitrary constants and f2zandgy,z are the two arbitrary functions.(II)
Noticing equation (106a), we assume that v=vξ,z,ξ=ax+by+c, then (106a) and (106b) become(111a)a1a2p+qvξξ+rΩ0=0,Ω0=Γ1−αΓ1−2α,(111b)vz=−ψzsa.Solving (111a) leads to(112)v=−rΩ02a1a2p+qξ2+ϕ1zξ+ϕ0z.Substituting (112) into (111b) yields(113)ϕ1′zξ+ϕ0′z=−ψzsa.By comparing the coefficients, we have(114)ϕ1z=c1,ϕ0z=−1sa∫zψz¯dz¯=−1saϕz,which gives rise to the exact solution(115)ut,x,y,z=a0tμ1+−rΓ1−α2a2p+qΓ1−2αax+by+c2+a1c1ax+by+c−a1saϕzt−α,α≠12,μ1≠α−1,wherea0,a1,c1,a,b,andc are the arbitrary constants and ϕz is an arbitrary function.Remark 10.
We do not consider the transformationξ=ax+by+cz+d for omitting the trivial solution.Remark 11.
The singularity only appears whenα=1/2andμ1=α−1 in solutions (110) and (115) in the Riemann–Liouville case.


## 3.3.3. Dynamical Analysis of Exact Solutions for Fractional JM Equation

Under the Caputo case:(i)
As can be seen from Figure10, taking α=1/2,μ1=α−1,a0=1,a1=1,c0=1,c1=2,c2=1,c3=−3,a=2,r=b=s=c=1,λ=−1,d=−3y=1.5,z=2,t∈0,10,x∈−10,10, it is shown that solution (83) increases with increase in the spatial variable and decreases as time increases on all (u, t, x), (u, t, y), and (u, t, z)-planes.(ii)
As Figures11 and 12 show, we see that by choosing a=b=1,s=7/2,c=2,r=2,α=1/2,λ=−1/2,a0=c0=1,a1=−2,b0=1,c1=2,c2=1,c3=−1.5,y=2,z=2,t∈0,10,x∈−10,10, and x=2.5,y=2,z∈−10,10, the solutions (89) and (91) increase as x and y increases and decrease as time increases on both (u, t, x) and (u, t, y)-planes; however, ones have exponential decay as z⟶∞ when x and y are fixed.Under the Riemann–Liouville case:As Figure13 shows, taking α=2/3,μ1=±3/2,a0=2,p=r=1,b1=−3,a1fz+a1gy,z+b0=10,x∈−100,100,t∈0,100, we see that(i)
Ifμ1>0, solutions (110) and (115) will blow up as time-space tends to infinity.(ii)
Ifμ1<=0, solutions (110) and (115) increase as x increases if y and z are fixed and decay at maxμ1,α rate as t⟶∞ if the spatial variables are fixed.


## 4. Conclusions

In this paper, by taking advantage of the fractional calculus methods, we avoid the invalid chain rule and suggest the new types of HBPs to solve some multidimensional time-fractional PDEs whose exact solutions were hardly obtained before. We get some explicit solutions of the (2 + 1)-dimensional KP equation, (3 + 1)-dimensional ZK equation, and JM equation by solving the reduced PDE system in both Riemann–Liouville and Caputo cases. These solutions are all in the general separated variable forms of new type which even include arbitrary functions, and the singularity solutions only appear in the Riemann–Liouville case. Furthermore, the dynamical analysis and long-time behaviors of these solutions are also performed.However, to the best of our knowledge, that is far from enough, the HBP is only applicable to solve some special fractional nonlinear PDEs which satisfy the balance conditions (see Section3.1) and get a few special exact solutions but not more general solutions. For even more multidimensional fractional nonlinear partial differential systems, there is still no better way to acquire their exact solutions generally at present. Nevertheless, finding more effective methods for constructing more exact solutions of more general (N + 1)-dimensional fnPDEs will be quite a meaningful and challenging task in the future.



---

*Source: 9101982-2020-05-06.xml*