Separation method of semi-fixed variables together with dynamical system method for solving nonlinear time-fractional PDEs with higher-order terms

It is well known that methods for solving fractional-order PDEs are grossly inadequate compared with integer-order PDEs. In this paper, a new approach combined with the separation method of semi-fixed variables and dynamical system method is introduced. As an example, a time-fractional reaction-diffusion equation with higher-order terms is studied under two different kinds of fractional-order differential operators. In different parametric regions, phase portraits of systems derived from the reaction-diffusion equation are presented. Existence and dynamic properties of solutions of this nonlinear time-fractional model are investigated. In some special parametric conditions, some exact solutions of this time-fractional models are obtained. The dynamical properties of some exact solutions are discussed, and the graphs of them are illustrated.


Introduction
It is well known that the concept of fractional-order derivative was first appeared in a famous correspondence between L'Hôspital and Leibniz in 1695. Since then, an effective theory of fractional-order calculus was gradually established with the joint efforts of many mathematicians during the past three hundred years, so far, there are more than dozen definitions on fractional-order derivative. Although the concepts of both fractional-order calculus and integer-order calculus originated in the Leibniz's era, their developments were uneven. Obviously, the theory of fractionalorder calculus develops very slowly and difficultly compared with the theory of integer-order calculus. One reason which restricts development of the theory of fractional-order calculus is that it has been lacked application background for a long time in the past. Another reason is that the theory of fractional-order calculus has been lacked effective methods and analytical tools on investigating solutions of fractionalorder differential equations because many classical (effective) methods in the integer-order field cannot be directly applied to fractional-order field. Since 1960s, the two issues have been greatly improved and developed. On the one hand, more and more researchers have found that the fractional-order differential models can be more accurately described complex problems in various fields such as mathematical mechanics, control theory, signal processing, aerodynamics, chem-istry, biology and so forth. In particular, fractionalorder calculus has been recognized as one of the best tools to describe long-memory processes, viscoelastic phenomena and behavior of anomalous diffusion in the past decades. Fractional-order differential models are not only interesting for science and engineering but also for pure mathematics. In the study of pure mathematics, in order to make up for the deficiency of fractional-order differential models, researchers sometimes directly transform the classical integer-order differential models into the corresponding fractionalorder differential models to study them. On the other hand, there appeared many effective methods for solving nonlinear fractional-order differential equations in recent decades. The main methods include Adomian decomposition method [1,2], homotopy analysis method [3,4], invariant analysis method [5,6], fractional variational iteration method [7][8][9], invariant subspace method [10][11][12], method of fractional complex transformation [13][14][15] and the method of separating variables [16][17][18], etc.
Although some exact and approximate analytical solutions of fractional-order differential equations can be obtained by above methods, it is still far from enough to solve more complex nonlinear fractional-order partial differential equations (PDEs). What's more, we find that the method of fractional complex transformation appeared in Refs. [13][14][15] is based on a wrong fractional chain rule given by Jumarie in Refs. [19][20][21]. In fact, Jumarie's fractional chain rule has been verified that it is invalid in Refs. [22][23][24]. This means peoples need to redesign some new methods to specifically target more complex fractional-order PDEs such as nonlinear time-fractional PDEs. Recently, based on the separation method of variables and combined with other methods such as homogenous balanced principle, idea of invariant subspace, integral bifurcation method and dynamical system method, we introduced several new methods [24][25][26][27][28] for solving nonlinear time-fractional PDEs. Obviously, in these methods, unlike the traditional separation method of variables is that we set the part of time function as some specific special functions such as Mittag-Leffler function and power function. We call this improved separation method of variables a separation method of semi-fixed variables. In particular, in Ref. [28], used this kind of improved separation method of variables and combined with bifurcation theory of differential dynamical system (dynamical system method) [29][30][31][32], we investi-gated existence and dynamical property of solutions of a nonlinear time-fractional PDE. But can this method solve nonlinear time-fractional PDEs with high-order terms? In this paper, we will introduce applications on the separation method of semi-fixed variables combining with dynamical system method. As an example, by using this combination method, we will investigate exact solutions, existence and dynamical properties of solution of time-fractional reaction-diffusion model with higher-order terms formed as the parameters δ, σ are two nonzero constants and κ is an arbitrary constant, and the sign ∂ α ∂t α defines a fractional-order differential operator of Caputo type where F(u) and K (u) are two arbitrary and smooth functions of u, and the sign ∂ α ∂t α is fractional-order differential operator of Riemann-Liouville type or Caputo type. When α → 1, Eq. (1.2) becomes the following integer-order reaction-diffusion equation of Fisher-KPP type as model of biological population given by Malaguti and Marcelli [33]. In particular, when F(u) = u m and K (u) = λu n (1 − u), Eq. (1.3) becomes a reactiondiffusion model with higher-order terms as follows: studied by Pablo and Vazquez [34]. In the conditions m > 1, λ = 1, n = 1, Aronson investigated sharptype solutions of Eq. (1.4) in [35].
In this paper, we will investigate exact solutions, existence and dynamical property of solutions of Eq. (1.1) from mathematical point of view. Because there are distinct differences between the Riemann-Liouville operator and Caputo operator, the solutions of Eq. (1.1) are completely different, and the difficulty for searching exact solutions is also different. In fact, under the two definitions of differential operator (Riemann-Liouville type and Caputo type), the fractional-order derivatives of the constant, Mittag-Leffler function and power function are different, the details which can be seen below.
where γ, α are constants and 0 < α < 1. The organization of this paper is as follows: In Sect. 2, under the definition of fractional-order differential operator of Caputo type, we will investigate exact solutions, existence and dynamical properties of solutions of (1.1). In Sect. 3, under the definition of fractional-order differential operator of Riemann-Liouville type, we will investigate exact solutions, existence and dynamical properties of solutions of (1.1) in two different cases.

Exact solutions and dynamical properties of (1.1) under the definition of Caputo operator
Under the definition of Caputo operator, Eq. (1.1) can be rewritten as where C 0 D α t defines a fractional-order differential operator of Caputo type, and parameters δ, σ, κ are nonzero constants. According to the classical separation method of variables, the form of solutions of Eq. (2.1) needs to be supposed as follows: where v(x) and T (t) are both undetermined functions.
Here, if we directly fix the time function T (t) to a Mittag-Leffler function, then (2.2) becomes where the function v = v(x) is undetermined function, the λ is undetermined coefficient, both of them can be determined later. Thus, this improved method will become simpler and the amount of computation is greatly reduced during the process for solving model (2.1). Obviously, this is appropriate since solutions of fractional differential equations often contain Mittag-Leffler function. We call this improved method the separation method of semi-fixed variables.
Equation (2.4) can be rewritten as Letting each coefficient of Mittag-Leffler functions in (2.5) to be zero, we obtain and Solving Eq. (2.6), it yields λ = κ δ . (2.8) Letting dv dx = y, the equation of (2.7) can be reduced to the following system dv dx = y, (2.9) Obviously, planar system (2.9) is a regular system when m = 1, but it is singular system when m ≥ 2. In other words, the derivative dy dx is not defined at v = 0 when m ≥ 2. However, v = 0 is a trivial solution of Eq. (2.7) when m ≥ 2. In other words, system (2.9) is not equivalent to Eq. (2.7) at v = 0 when m ≥ 2. In order to obtain a completely equivalent system of Eq. (2.7) no matter how the function v varies when m ≥ 2, we make a scalar transformation as follows: dx = mv m−1 dτ, (2.10) where τ is a parameter. Under transformation (2.10), singular system (2.9) can be reduced to a regular system as follows: (2.11) Obviously, both systems (2.9) and (2.11) have the same first integral as follows: where h is an integral constant. In order to facilitate discussion, Eq. (2.12) can be rewritten as From (2.11) and (2.13), we easily verify that ∂ H ∂ y = − dv dτ , ∂ H ∂v = dy dτ . Therefore, system (2.11) is not a Hamiltonian system. For convenience on discussion, we write Jacobian matrix and Jacobian determinant of (2.11) as follows: 14) (2.15) Obviously, system (2.11) has only one equilibrium points O(0, 0) at the v-axis. Substituting the equilibrium points O(0, 0) into (2.13) and (2.15), it yields According to theory of planar dynamical systems [29][30][31][32] and Lemma 1 in [28], we know that the equi- Obviously, the phase portraits of system (2.9) can be constructed by using the phase portraits of system (2.11) together the singular line v = 0. According to the above information, we draw graphs of phase portraits for system (2.9) in four cases, which can be shown in Fig. 1.
By means of phase portraits in Fig. 1, we obtain some information about distribution of orbits for system (2.9) as follows: When m = 1, σ < 0 and h > 0, system (2.11) has infinite numbers of closed orbits, which can be shown in Fig. 1a. When m = 1, σ > 0 and h = 0, system (2.11) has two orbits of line type marked by black color, which can be shown in Fig. 1b. When m = 1, σ > 0 and h < 0, system (2.11) has infinite numbers of hyperbolic orbits (two pairs of them marked by green color), which can be shown in Fig. 1b. When m = 1, σ > 0 and h > 0, system (2.11) has infinite numbers of hyperbolic orbits (two pairs of them marked by red color), which can be shown in Fig. 1b. In particular, when m ≥ 2, the equilibrium point O(0, 0) occurs degradation phenomenon due to the influence of the singular straight line, the orbits are divided into left and right parts by the singular line v = 0, and the original closed orbits are no longer closed. According to the above information about distribution of orbits, we get two theorems below.  Different kinds of solutions mentioned above theorems satisfy u → 0 in the condition of δκ < 0 when t → +∞. On the contrary, they satisfy u → ±∞ in the condition of δκ > 0 when t → +∞. (2.17) Obviously, each closed orbit in Fig. 1a (2.20) (ii) When m = 1 and σ > 0, h = 0, Eq. (2.12) can be reduced to Substituting (2.21) into the first equation dv dx = y of (2.9) to directly integrate it, we get two general solutions of exponential function type as follows: (2.26) Respectively, taking − h σ , 0 and − − h σ , 0 as initial points, then substituting (2.26) into the first equation dv dx = y of (2.9) and integrating them, we get two solutions of hyperbolic cosine function as follows: (2.28) When m = 1 and σ > 0, h > 0, Eq. (2.12) can be reduced to (2.29) Substituting (2.29) into the equation dv dx = y to directly integrate it, we get two general solutions as follows: where c is arbitrary constant. In particular, when c = √ σ , solutions (2.30) and (2.31) become the following two types of hyperbolic sine function.
(2.12) can be rewritten as Substituting (2.36) into the first equation dv dx = y of (2.9) to integrate it, we obtain two periodic solutions of Eq. (2.7) as follows: (ii) When m ≥ 2, σ > 0 and h = 0, Eq. (2.12) can be rewritten as Substituting (2.41) into the first equation dv dx = y of (2.9) to directly integrate it, we get two general solutions of exponential function type as follows: (2.45) (iii) When m ≥ 2, σ > 0 and h < 0, Eq. (2.12) can be rewritten as , 0 as initial points, then substituting (2.46) into the first equation dv dx = y of (2.9) and integrating them, we get two solutions of hyperbolic cosine function as follows: When m ≥ 2, σ > 0 and h > 0, Eq. (2.12) can be rewritten as  Solutions (2.20) and (2.39) and their profiles show that both of them have periodic property, but the former is smooth and the latter is nonsmooth. This is because system (2.9) increases a singular line v = 0 when m ≥ 2.

Exact solutions and dynamical properties of (1.1) under the definition of Riemann-Liouville operator
In this section, under the definition of Riemann-Liouville operator, we will investigate exact solutions, existence and dynamical properties of solutions of Eq.
It is well known that the Riemann-Liouville operator has singularity compared with the Caputo operator, which can be seen by the derivatives of the constant and Mittag-Leffler function given in the previous section. In fact, we cannot obtain exact solutions of Eq.
We suppose that Eq. (3.2) has solutions formed as follows: where the function v = v(x) is undetermined function and γ is undetermined constant. Substituting (3.3) into (3.2), we get According to homogenous balanced principle, in Eq.
Letting dv dx = y, Eq. (3.7) can be reduced to the follow- Also, the derivative dy dx cannot be defined at v = 0. Similarly, we make a scalar transformation as follows: dx = mv m−1 dτ, (3.9) where τ is a parameter. Under transformation (3.9), system (3.8) can be reduced to dv dτ = mv m−1 y, (3.10) Obviously, systems (3.8) and (3.10) have same first integral as follows: where h is an integral constant. As in Sect. 2, we rewrite (3.11) as Similarly, the Jacobian matrix and Jacobian determinant of system (3.10) are written as follows: From (3.10) and (3.12), we easily verify − ∂ H ∂ y = dv dτ , ∂ H ∂v = dy dτ . So, system (3.10) is not a Hamiltonian system yet. In addition, we know that system (3.10) has two equilibrium points O(0, 0) and According to theory of planar dynamical systems [29][30][31][32] and Lemma 1 in [28], we easily obtain some information about equilibrium points of system (3.10) as follows: The original point O is a degenerative equilibrium point because J 0 = 0 and its Poincaré index does not equal zero; it is a degenerative saddle point if σ > 0 and it is a degenerative center point if σ < 0. When m is even number and σ < 0, the equilibrium point A is a center point; it is on the left side of the y-axle if δ 0 < 0, it is on the right side of the y-axle if δ 0 > 0. When m is even number and σ > 0, the equilibrium point A is a saddle point; it is on the left side of the y-axle if δ 0 > 0, it is on the right side of the y-axle if δ 0 < 0. When m is odd number and σ < 0, δ 0 > 0, the equilibrium points B 1,2 are two center points. When m is odd number and σ > 0, δ 0 < 0, the equilibrium points B 1,2 are two saddle points.
Obviously, the phase portraits of system (3.8) can be constructed by the phase portraits of system (3.10) together with the singular line v = 0. According to the above information, we draw graphs of the phase portraits for system (3.8), which can be shown in Figs. 4 and 5.
According to the above information about distribution of orbits, we get two theorems as follows:  We can only obtain exact expressions of these solutions mentioned in above theorems in some particular parameter conditions because the order number m of functions is too high. Therefore, we only prove some of the propositions in theorems 3.1 and 3.2.
Theoretically, these kinds of periodic solutions can be obtained by substituting (3.18) into (3.8), but the order number m is too high to obtain exact expressions of solutions by integrating, indeed. However, in the case of h = 0 or m = 2, we can easily obtain an exact expression for this kind of periodic solution, see the next discussion. If m is even number and h = 0, σ < 0, then Eq. (3.18) can be reduced to Taking (v 0 , 0) as initial point, substituting (3.19) into the first equation of (3.8) and then integrating it, we get Solving (3.20), we obtain a smooth periodic solution of (3.7) as follows: m . This implies that Eq. (3.2) has not the type of solution formed as (3.3) when α = m−1 m ; maybe it has other type of solution when α = m−1 m , but we cannot obtain its expression by current method at here.
(3.23) From Fig. 4a, b, we can see that all closed orbits have two intersection points with the v-axis. We suppose the two intersection points expressed by Q 1 (φ 1 , 0) and Q 2 (φ 2 , 0) and write φ 1 > φ 2 . In fact, v = φ 1 , φ 2 are two real roots of the following equation There is no doubt that Eq. (3.24) has other two conjugate complex roots, and we write as v = s,s. We always can obtain expressions of the four roots φ 1 , φ 2 , s,s by using Ferrari method or Descartes method, but their expressions are more complex (tediously long), so we omit them at here. Once the parameters h, δ, σ, 0 are evaluated, we can easily obtain expressions of these roots with the help of the computer software. For example, giving the values of parameters h = −1, α = 0.75, δ = 2, σ = −1, by using computer to solve Eq.
Taking (φ 2 , 0) as initial conditions, substituting (3.25) into the first equation of (3.10) and then integrating, it yields By using the formula of Jacobian elliptic integral in [36] to integrate (3.26), we get where cn −1 ( * , * ) is inverse function of the Jacobian elliptic function cn( * , * ), the parameters A = By the way, the Jacobian elliptic function cn( * , * ) is an even function, so Eq. (3.26) has only one solution final. After doing the inverse operation of the above equation, we obtain a general solution as follows: v = Aφ 2 + Bφ 1 A + B 1 + P 1 cn(ωτ,k) 1 + P 2 cn(ωτ,k) , (3.27) where τ is a parameter and ω = √ −σ AB,k = Substituting (3.27) into (3.9) and then integrating it by integral formula of Jacobian elliptic function in [36], we get where the am(ωτ,k), In order to show dynamical property of above solutions intuitively, as an example, we plot 3D-graphs of dynamical profiles of solution (3.22), which are shown in Fig. 6a, b.
Compared with solution (2.20), solution (3.22) always satisfies u → 0 when t → +∞, but the convergence rate of solution (3.22) is obviously not as fast as solution (2.20) during the process of t → +∞. This implies that the properties of the solutions of Eq. (1.1) are completely different under the definitions of two different fractional differential operators.
Here, we omit the proof of the proposition (ii) because the method and the result are very similar to the former.
(iii). If m is even number and h = 0, σ > 0, Eq. We rewrite (3.33) as follows: where a 1 = 0 δ, b 1 = σ, a = 2 0 δ 2 , b = −2 0 δσ, c = 3σ 2 . Substituting (3.34) into the first equation of (3.8) and then integrating it, we get vdv Completely integrating (3.35) and setting the integral constant as zero, we obtain two implicit solutions of (3.7) as follows: 3) and combining with (3.36), we get two exact solutions of (3.2) as follows: (3.37) If m is odd number and h = 0, σ > 0, δ 0 < 0, then Eq. (3.11) can be reduced to (3.43) Solving (3.43) and setting the integral constant is zero, we obtain two solutions of hyperbolic function type for (3.7) as follows: If m is odd number and h = 0, σ > 0, δ 0 > 0, then Eq. (3.11) can be reduced to Obviously, there is not any intersection point of the two open orbits and the v-axis. Directly substituting (3.46) into the first equation of (3.8), then integrating it, we get By using the improvising differentiation method to integrate (3.47) and setting the integral constant is zero, we obtain two solutions of hyperbolic function type for (3.7) as follows: According to traditional separation method of variables, we suppose that Eq. (3.50) has solutions formed as follows: where v(x) and T (t) are two undetermined functions. Substituting (3.51) into (3.50), it can be reduced to where On the other hand, it is easy to know that Eq. (3.55) has characteristic equation as follows: its characteristic root (eigenvalue) is a real double root Z = 0 if λ = κ − σ, its characteristic roots are two real roots Z = ± √ λ + σ − κ if λ > κ − σ, and its characteristic roots are two pure virtual roots (complex roots) Z = ±i √ κ − λ − σ if λ < κ − σ. According to the above information, we obtain three kinds of solutions of ODE (3.55) as follows: where C 1 , C 2 are arbitrary constants. Thus, Eq. (3.50) has three kinds of exact solutions as follows: where c 1 = C 0 C 1 , c 2 = C 0 C 2 are two arbitrary constants.
In order to show dynamical property of above solutions intuitively, as an example, we plot 3D-graphs of dynamical profiles of solution (3.63), which are shown in Fig. 8a, b.

Conclusions
In this work, we improved the traditional separation method of variables, and we call this improved method the separation method of semi-fixed variables. Compared with the traditional separation method of variables, this improved method will become simpler and the amount of computation is greatly reduced during the process for solving a time-fractional PDE. By using this separation method of semi-fixed variables together with dynamical system method, we successfully studied a time-fractional reaction-diffusion equations with higher-order terms under two kinds of fractionalorder differential operators (Caputo differential operator and Riemann-Liouville differential operator). Different kinds of exact solutions of this time-fractional reaction-diffusion equation with higher-order terms are obtained. In addition, we find that the solutions and their dynamical properties of Eq. (1.1) are completely different under the definitions of two different fractional differential operators. All solutions of Eq. (1.1) can be divided into two major categories, one with power functions and another with Mittag-Leffler functions. Those solutions which contain power functions, all of them have attenuation properties and satisfy u → 0 as t → +∞ because their powers are negative. Those solutions contain Mittag-Leffler functions, and some of them have attenuation property and satisfy u → 0 in the condition of δκ < 0 as t → +∞ and others have incremental property and satisfy u → ±∞ in the condition of δκ > 0 as t → +∞. In order to show dynamical properties of these exact solutions intuitively, the 3D-graphs of dynamical profiles of some representative solutions are illustrated.