Solving Partial Differential Equation with Space- and Time-Fractional Derivatives via Homotopy Decomposition Method

Abstract

The analytical solution of the partial differential equation with time- and space-fractional derivatives was derived by means of the homotopy decomposition method (HDM). Some examples are given and comparisons are made. The evaluations show that the homotopy decomposition method is extremely successful and suitable. The achieved results make the steadfastness of the HDM and its wider applicability to fractional differential equation obvious. Additionally, the adding up implicated in HDM is exceptionally undemanding and uncomplicated. It is confirmed that HDM is an influential and professional apparatus for FPDEs. It was also established that HDM is supplementary well organized than the ADM, VIM, HAM, and HPM.

1. Introduction

Many observable fact in natural science, physics, chemistry, and other knowledge preserve to be illustrated incredibly fruitfully by representations using the supposition of derivatives and integrals with fractional order. Attention in the notion of differentiation and integration to noninteger order has existed since the progress of the conventional calculus [1–3]. By insinuation, mathematical modeling of many physical systems is governed by linear and nonlinear fractional differential equations in various applications in fluid mechanics, viscoelasticity, chemistry, physics, biology, and engineering. Since many fractional differential equations are nonlinear and do not have exact analytical solutions, various numerical and analytic methods have been used to solve these equations [4]. The Adomian decomposition method (ADM) [4], the homotopy perturbation method (HPM) [5], the variational iteration method (VIM) [6], homotopy analysis method [7, 8], and other methods have been used to provide analytical approximation to linear and nonlinear problems [9, 10]. However, the convergence region of the corresponding results is rather small [9, 10]. One of the popular FDE applications is anomalous diffusion; however, this problem can be effectively solved by radial basis function (RBF) collocation methods [11–13]. In this paper, we use homotopy decomposition method that was recently proposed in [14] to solve partial differential equation that arises in groundwater flow problem. The method was first used to solve time-fractional coupled-korteweg-de-vries equations [15]. The homotopy decomposition method is chosen to solve this nonlinear problem because of the following advantages that the method has over the existing methods. The method does not require the linearization or assumptions of weak nonlinearity [5]. The solutions are not generated in the form of general solution as in the Adomian decomposition method. With ADM, the recursive formula allows repetition of terms in the case of nonhomogeneous partial differential equation, this repetition leads to the noisy solution [16]. The solution obtained is noise-free compared to the variational iteration method [16]. No correctional function is required as in the case of the variational homotopy decomposition method [16]. No Lagrange multiplier is required in the case of the variational iteration method [6]. this method is more realistic compared to the method of simplifying the physical problems. If the exact solution of the partial differential equation exists, the approximated solution via the method converges to the exact solution [14]. A construction of a homotopy 𝑣(𝑟,𝑝)∶Ω×[0,1] is not needed as in the case of the homotopy perturbation method, because in this case one needs first to continuously deform a difficult problem into another one, which is easy to solve [6]. HDM provides us with a convenient way to control the convergence of approximation series without adapting ℎ, as in the case of [17], which is a fundamental qualitative difference in analysis between HDM and other methods.

In this paper, two cases of special interest such as the time-fractional foam drainage equation and the space-fractional foam drainage equation are discussed in detail. The paper has been arranged as follows. History of the fractional derivative order is presented in Section 2. In Section 3 the homotopy decomposition method for solving fractional derivative order is described. Complexity of the homotopy decomposition method is discussed in Section 4. In Section 5 the application of HDM to solve the space- and time-fractional partial differential equation is discussed. Conclusions are presented in Section 6.

2. History of the Fractional Derivative Order

There is in the literature numereous definitions about fractional derivatives [18–23]. The most popular ones are the Riemann-Liouville and the Caputo derivatives. For Caputo derivative we have ^𝐶₀𝐷^𝛼_𝑥𝑓(𝑥)=1Γ(𝑛−𝛼)^𝑥₀(𝑥−𝑡)^{𝑛−𝛼−1}𝑑^𝑛𝑓(𝑡)𝑑𝑡^𝑛𝑑𝑡. (1)

For the case of Riemann-Liouville derivative we have the following definition: 𝐷^𝛼_𝑥𝑓(𝑥)=1Γ()𝑑𝑛−𝛼𝑑𝑥^𝑛𝑛𝑥₀(𝑥−𝑡)^{𝑛−𝛼−1}𝑓(𝑡)𝑑𝑡. (2)

Lately, Jumarie (see [21, 22]) proposed a simple alternative definition to the Riemann-Liouville derivative: 𝐷^𝛼_𝑥𝑓(𝑥)=1Γ()𝑑𝑛−𝛼𝑑𝑥^𝑛𝑛×^𝑥₀(𝑥−𝑡)^{𝑛−𝛼−1}𝑓(𝑡)−𝑓(0)𝑑𝑡. (3)

3. Basic Idea of the HDM [14, 15]

To demonstrate the fundamental design of this technique, we reflect on a universal nonlinear nonhomogeneous fractional partial differential equation with initial conditions of the following structure: 𝜕^𝛼𝐺(𝑥,𝑡)𝜕𝑡^𝛼=𝐿(𝐺(𝑥,𝑡))+𝑁(𝐺(𝑥,𝑡))+𝑙(𝑥,𝑡),𝛼>0, (4) subject to the initial condition 𝐷₀^𝛼−𝑘𝐺(𝑥,0)=𝑓_𝑘(𝑥),(),𝐷𝑘=0,…,𝑛−1_0𝛼−𝑛𝐺()[𝛼],𝐷𝑥,0=0,𝑛=_𝑘0𝐺(𝑥,0)=𝑔_𝑘(𝑥),(),𝐷𝑘=0,…,𝑛−1_𝑛0𝐺(𝑥,0)=0,𝑛=[𝛼], (5) where, 𝜕^𝛼/𝜕𝑡^𝛼 indicates the Caputo or Riemann-Liouville fractional derivative operator, 𝑙 is a known function, 𝑁 is the general nonlinear fractional differential operator, and 𝐿 represents a linear fractional differential operator. The process of the HDM primary pace here is to change the fractional partial differential equation to the fractional partial integral equation by applying the inverse operator 𝜕^𝛼/𝜕𝑡^𝛼 of on both sides of (4) to obtain [15]. In the case of Riemann-Liouville fractional derivative [15] 𝐺(𝑥,𝑡)=^𝑛−1_𝑗=1𝑓_𝑗(𝑥)Γ𝑡𝛼−𝑗+1^𝛼−𝑗+1Γ(𝛼)×^𝑡₀(𝑡−𝜏)^𝛼−1×[𝐿(𝐺(𝑥,𝜏))+𝑁(𝐺(𝑥,𝜏))+𝑙(𝑥,𝜏)]𝑑𝜏. (6)

In the case of Caputo fractional derivative 𝐺(𝑥,𝑡)=^𝑛−1_𝑗=1𝑔_𝑗(𝑥)Γ𝑡𝛼−𝑗+1^𝛼−𝑗+1Γ(𝛼)×^𝑡₀(𝑡−𝜏)^𝛼−1×[𝐿(𝐺(𝑥,𝜏))+𝑁(𝐺(𝑥,𝜏))+𝑙(𝑥,𝜏)]𝑑𝜏, (7) or in general by putting ^𝑛−1_𝑗=1𝑓_𝑗(𝑥)Γ𝑡𝛼−𝑗+1^𝛼−𝑗()𝑓()==𝑙𝑥,𝑡or𝑥,𝑡^𝑛−1_𝑗=1𝑔_𝑗(𝑥)Γ𝑡𝛼−𝑗+1^𝛼−𝑗, (8) we obtain 𝑈()()+1𝑥,𝑡=𝑇𝑥,𝑡Γ(𝛼)×^𝑡₀(𝑡−𝜏)^𝛼−1×𝐿(𝑈(𝑥,𝜏))+𝑁(𝑈(𝑥,𝜏))+𝑓(𝑥,𝜏)𝑑𝜏. (9)

In the homotopy perturbation method, the fundamental statement is that the solutions can be written as a power series in 𝑝: 𝐺𝑥,𝑡,𝑝=^∞_𝑛=0𝑝^𝑛𝐺_𝑛(),𝐺()𝑥,𝑡𝑥,𝑡=lim_𝑝→1𝐺𝑥,𝑡,𝑝. (10) The nonlinear term can be decomposed as 𝑁𝐺(𝑥,𝑡)=^∞_𝑛=0𝑝^𝑛ℋ_𝑛(𝐺), (11) where 𝑝∈(0,1] is an embedding parameter. ℋ_𝑛(𝐺) is the He’s polynomials [17] that can be generated by ℋ_𝑛𝐺₀,…,𝐺_𝑛=1𝜕𝑛!𝜕𝑝^𝑛𝑛𝑁^∞_𝑗=0𝑝^𝑗𝐺_𝑗(𝑥,𝑡),𝑛=0,1,2…. (12)

The homotopy decomposition method is obtained by the refined combination of homotopy technique with He’s polynomials [15] and is given by ^∞_𝑛=0𝑝^𝑛𝐺_𝑛()()=𝑝𝑥,𝑡−𝑇𝑥,𝑡Γ(𝛼)^𝑡₀(𝑡−𝜏)^𝛼−1𝑙(𝑥,𝜏)+𝐿^∞_𝑛=0𝑝^𝑛𝐺_𝑛(𝑥,𝜏)+𝑁^∞_𝑛=0𝑝^𝑛𝐺_𝑛(𝑥,𝜏)𝑑𝜏. (13)

Putting side by side expressions of identical powers of 𝑝 provides solutions of a variety of orders with the first term [15] 𝐺₀(𝑥,𝑡)=𝑇(𝑥,𝑡). (14)

4. Complexity of the Homotopy Decomposition Method

It is very important to test the computational complexity of a method or algorithm [15]. Complexity of an algorithm is the study of how long a program will take to run, depending on the size of its input and long of loops made inside the code [15]. We compute a numerical example which is solved by the homotopy decomposition method. The code has been presented with Mathematica 8 according to the following code [15].

Step 1. Set 𝑚←0.

Step 2. Calculate the recursive relation after the comparison of the terms of the same power is done.

Step 3. If ‖𝑈_𝑛+1(𝑥,𝑡)−𝑈_𝑛(𝑥,𝑡)‖<𝑟 with 𝑟 being the ratio of the neighbourhood of the exact solution [15], then go to Step 4, else 𝑚←𝑚+1 and go to Step 2.

Step 4. Print out 𝐺(𝑥,𝑡)=^∞_𝑛=0𝑈_𝑛(𝑥,𝑡) (15) as the approximate of the exact solution.

Lemma 1. If the exact solution of the fractional partial differential equation (4) exists, then 𝐺_𝑛+1(𝑥,𝑡)−𝐺_𝑛(𝑥,𝑡)<𝑟,∀(𝑥,𝑡)∈𝑋×𝑇. (16)

Proof. Let (𝑥,𝑡)∈𝑋×𝑇; since the exact solution exists, then we have the following: 𝐺_𝑛+1(𝑥,𝑡)−𝐺_𝑛(=𝐺𝑥,𝑡)_𝑛+1(𝑥,𝑡)−𝐺(𝑥,𝑡)+𝐺(𝑥,𝑡)−𝐺_𝑛(≤𝐺𝑥,𝑡)_𝑛+1(𝑥,𝑡)−𝐺(𝑥,𝑡)+−𝐺_𝑛()(≤𝑟𝑥,𝑡+𝐺𝑥,𝑡)2+𝑟2=𝑟. (17)

The last inequality follows from [15].

Lemma 2. The complexity of the homotopy decomposition method is of order 𝑂(𝑛).

Proof. The number of computations including product, addition, subtraction, and division are the following.
In Step 2  𝐺₀: 0 because it is obtained directly from the initial guess [15],  𝐺₁: 3,  𝐺_𝑛: 3.
Now in Step 4 the total number of computations is equal to ∑^𝑛_𝑗=0𝐺_𝑗(𝑥,𝑡)=3𝑛=𝑂(𝑛).

5. Application

In this section we apply this method for solving partial differential equation with time- and space-fractional derivatives.

Example 3. Let us consider the fractional Riccati differential equation 𝐷^𝛼₀𝑦(𝑥)=−𝑦(𝑥)²+1,0<𝛼≤1,𝑥>0, (18) subject to the initial condition 𝑦(0)=0. (19)
Following the discussion presented in Section 3, we obtain the following: ^∞_𝑛=0𝑝^𝑛𝑦_𝑛(𝑥)(0)−𝑝=𝑦Γ(𝛼)^𝑥₀(𝑥−𝑡)^𝛼−1×^∞_𝑛=0𝑝^𝑛𝑦_𝑛(𝑡)²−1𝑑𝑡. (20)
Comparing the terms of the same power of 𝑝 and using (18), we obtain the following integral equations: 𝑝⁰:𝑦₀(𝑥)=𝑦(0)=0,𝑦₀(0)𝑝=0,₁:𝑦₁(𝑥)1=−Γ(𝛼)^𝑥₀(𝑥−𝑡)^𝛼−1𝑦²₀(𝑡)𝑦−1𝑑𝑡,₁(0)𝑝=0,_𝑛:𝑦_𝑛(𝑥)1=−Γ(𝛼)×^𝑥₀(𝑥−𝑡)^{𝛼−1𝑛−1}_𝑗=0𝑦_𝑗𝑦_{𝑛−𝑗−1}𝑦𝑑𝑡,_𝑛(0)=0,𝑛≥2. (21)
And the following solutions are obtained: 𝑦₀(𝑥)𝑦=0,₁(𝑥)=𝑥^𝛼Γ(),𝑦𝛼+1₂(𝑥)𝑦=0,₃(𝑥)𝑥=−^3𝛼Γ(1+2𝛼)Γ(1+3𝛼)Γ²(),𝑦1+𝛼₄(𝑥)𝑦=0,₅(𝑥)=2𝑥^5𝛼Γ(1+2𝛼)Γ(1+4𝛼)Γ³()Γ()Γ(),𝑦1+𝛼1+3𝛼1+5𝛼₆(𝑥)𝑦=0,₇(𝑥)=−𝑥^7𝛼Γ()×(()Γ()()Γ(Γ()×Γ1+2𝛼4Γ1+3𝛼1+4𝛼+Γ1+2𝛼1+5𝛼))1+6𝛼₃(1+𝛼)Γ²(1+3𝛼)Γ(1+5𝛼)Γ(1+7𝛼)⁻¹,𝑦₈(𝑥)𝑦=0,₉(𝑥)=𝑥^9𝛼Γ⁵(1+𝛼)Γ³()()Γ()Γ(1+3𝛼×Γ1+5𝛼1+7𝛼1+9𝛼)⁻¹×2𝑥^9𝛼Γ()×()Γ()Γ()(Γ()Γ()()Γ((⋮1+2𝛼4Γ1+3𝛼1+4𝛼1+6𝛼+Γ1+2𝛼)(1+4𝛼1+6𝛼+2Γ1+4𝛼1+7𝛼))×Γ1+8𝛼) (22)
Using the package Mathematica, in the same manner, one can obtain the rest of the components. But, here, 10 terms were computed and the asymptotic solution is given by 𝑦_𝑁=9(𝑥)=⁹_𝑛=0𝑦_𝑛(𝑥). (23)
Now to access the accuracy of HDM, we compare the approximated solution (23) when 𝛼=1 with the exact solution given as follows: 𝑦(𝑥)=tanh(𝑥). (24)
Notice that when 𝛼=1, 𝑦_𝑁=9(𝑥)=⁹_𝑛=0𝑦_𝑛(𝑥)𝑥=𝑥−³3+2𝑥⁵−1517𝑥⁷+31562𝑥⁹.2835 (25)
Therefore, for any 𝑁≥9, we have 𝑦_𝑁(𝑥)=^𝑁_𝑛=02^2𝑛2^2𝑛−1𝑥^2𝑛−1()!𝐵2𝑛^2𝑛,|𝑥|<𝜋2, (26) with 𝐵_2𝑛 being the Bernoulli number. Now also notice that the Taylor series of tanh(𝑥) of order 10 at the origin gives exactly 𝑦_𝑁=9(𝑥). Thus, it is worth concluding that when 𝛼=1 the approximated solution converges to the exact solution of Example 3.

Remark 4. Example 3 has been solved using HAM [24], ADM [25], VIM [26], and HPM [27, 28], and recently it was solved in [29]. Comparison of the results of applying the HDM with others methods leads to significant consequences.
An Eton proverb says “an image is equivalent to one thousand words.” The following figure shows the comparison of the approximated solution and the exact solution.

Figure 1 shows that the fractional solution decrease in the fractional order 𝛼 causes an increase in the function. HDM solutions indicate that the present algorithm performs with extreme efficiency, simplicity, and reliability. The results obtained from HDM are fully compatible with the exact solution when 𝛼=1, and are more accurate than those obtained from existing methods.

Figure 1 

Comparison of approximated and exact solutions for 𝑁=15.

Example 5. Consider the following form of the time-fractional equation: 𝐷^𝛼_𝑡1𝑢=2𝑢𝑢_𝑥,𝑥−2𝑢²𝑢_𝑥+𝑢²_𝑥,0<𝛼≤1,𝑥>0, (27) with the initial condition 𝑢(𝑥,0)=−√𝑐tanh√𝑐𝑥. (28)
Following the steps of HDM, we obtain ^∞_𝑛=0𝑝^𝑛𝑢_𝑛()()+𝑥,𝑡=𝑢𝑥,0^𝑦₀12^∞_𝑛=0𝑝^𝑛𝑢_𝑛𝑥,𝑦^∞_𝑛=0𝑝^𝑛𝑢_𝑛𝑥,𝑦_𝑥,𝑥−2^∞_𝑛=0𝑝^𝑛𝑢_𝑛𝑥,𝑦^2∞_𝑛=0𝑝^𝑛𝑢_𝑛𝑥,𝑦_𝑥+^∞_𝑛=0𝑝^𝑛𝑢_𝑛𝑥,𝑦_𝑥𝑑𝑦². (29) Comparing the terms of the same power of 𝑝 and using (27), we obtain the following: 𝑝⁰:𝑢₀(𝑥,𝑡)=𝑢(𝑥,0)=−√𝑐tanh√,𝑝𝑐𝑥₁:𝑢₁(𝑥,𝑡)=^𝑡₀12𝑢₀𝑢_0𝑥,𝑥−2𝑢₀²𝑢_0𝑥+𝑢_0𝑥²𝑢𝑑𝜏,₁()𝑝𝑥,0=0,_𝑛:𝑢_𝑛(𝑥,𝑡)=^𝑡₀12^𝑛−1_𝑗=0𝑢_𝑗𝑢_{𝑛−1−𝑗𝑥,𝑥}−2^𝑛−1𝑗_𝑗=0𝑘=0𝑢_𝑗𝑢_𝑗−𝑘𝑢_{𝑛−1−𝑗𝑥}+^𝑛−1_𝑗=0𝑢_𝑗𝑥𝑢_{𝑛−𝑗−1𝑥}𝑢𝑑𝑦,_𝑛(𝑥,0)=0,𝑛≥2. (30) The following solutions are obtained: 𝑢₀(𝑥,𝑡)=𝑢(𝑥,0)=−√𝑐tanh√,𝑢𝑐𝑥₁(𝑥,𝑡)=√𝑐𝑐^3/2−𝑐^3/2tanh√𝑐𝑥²𝑡^𝛼Γ(),𝑢𝛼+1₂(𝑥,𝑡)=√𝑐12𝑐³tanh√𝑐𝑥−tanh√×−𝑐𝑐𝑥₃2+𝑐³tanh√𝑐𝑥²𝑡^2𝛼,𝑢₃(𝑥,𝑡)=√𝑐⁹√𝑐6−12⁹√𝑐tanh√𝑐𝑥²+√³𝑐−𝑐³2+𝑐³tanh√𝑐𝑥²−tanh√×−1𝑐𝑥3⁹√𝑐tanh√𝑐𝑥+√³𝑐tanh√×−𝑐𝑐𝑥₃2+𝑐³tanh√𝑐𝑥²𝑡^3𝛼,𝑢₄()=√𝑥,𝑡𝑐18𝑐⁶tanh√−1𝑐𝑥2𝑐³tanh√×−𝑐𝑐𝑥₃2+𝑐³tanh√𝑐𝑥²+√³𝑐−13⁹√𝑐tanh√𝑐𝑥+√³𝑐tanh√×−𝑐𝑐𝑥₃2+𝑐³tanh√𝑐𝑥²−tanh√×−𝑐𝑐𝑥₆+𝑐24⁶tanh√𝑐𝑥²6+12𝑐³−𝑐³2+𝑐³tanh√𝑐𝑥²−√³𝑐tanh√×𝑐𝑥₉√𝑐tanh√𝑐𝑥3³√𝑐tanh√−𝑐𝑐𝑥₃2+𝑐³tanh√𝑐𝑥²𝑡^4𝛼. (31)
Using the package Mathematica, in the same manner, one can obtain the rest of the components. But in this case, 5 terms were computed and the asymptotic solution is given by 𝑢_𝑁=5(𝑥,𝑡)=𝑢₀(𝑥,𝑡)+𝑢₁(𝑥,𝑡)+𝑢₂(𝑥,𝑡)+𝑢₃(𝑥,𝑡)+𝑢₄(𝑥,𝑡)+⋯. (32)
Now notice that if we set 𝛼=1, we recover the Taylor series of √𝑐tanh(√𝑐(𝑥−𝑐𝑡)) of order 4. Therefore, if 𝛼=1 using the package Mathematica, we recover ^∞_𝑛=0𝑢_𝑛(𝑥,𝑡)=−√𝑐tanh√𝑐(𝑥−𝑐𝑡). (33)
This is the exact solution of Example 5.

In this section, to access the accuracy of the HDM, we compare the approximation (32) with the exact solution; this is depicted in Figures 2 and 3.

Figure 2 

Exact solution for 𝑐=1/2.

Figure 3 

Approximated solution for 𝑐=1/2 and 𝑁=6.

From the Figure 2, it is obvious that when 𝛼=1, the solution is nearly identical with the exact solution.

Example 6. Consider the fractional predator-prey equation 𝜕^𝛼𝑢𝜕𝑡^𝛼=𝜕²𝑢²𝜕𝑥²+𝜕²𝑢²𝜕𝑦²𝑢+𝑎𝑢−𝑏𝑢𝑣,𝑥,𝑦,0=𝑒^𝑥+𝑦,𝜕^𝛽𝑣𝜕𝑡^𝛽=𝜕²𝑣²𝜕𝑥²+𝜕²𝑣²𝜕𝑦²𝑣+𝑏𝑣𝑢−𝑐𝑣,𝑥,𝑦,0=𝑒^𝑥+𝑦. (34) Following the homotopy decomposition method, presented in Section 3, we obtain the following system: ^∞_𝑛=0𝑝^𝑛𝑢_𝑛+𝑝=𝑢𝑥,𝑦,0Γ(𝛼)^𝑡₀𝜕_𝑥𝑥^∞_𝑛=0𝑝^𝑛𝑢_𝑛²+𝜕_𝑦𝑦^∞_𝑛=0𝑝^𝑛𝑢_𝑛²+𝑎^∞_𝑛=0𝑝^𝑛𝑢_𝑛−𝑏^∞_𝑛=0𝑝^𝑛𝑢_𝑛^∞_𝑛=0𝑝^𝑛𝑣_𝑛×(𝑡−𝜏)𝑑𝜏,_∞𝑛=0𝑝^𝑛𝑣_𝑛+𝑝=𝑣𝑥,𝑦,0Γ𝛽×^𝑡₀𝜕_𝑥𝑥^∞_𝑛=0𝑝^𝑛𝑣_𝑛²+𝜕_𝑦𝑦^∞_𝑛=0𝑝^𝑛𝑣_𝑛²+𝑏^∞_𝑛=0𝑝^𝑛𝑢_𝑛^∞_𝑛=0𝑝^𝑛𝑣_𝑛−𝑐^∞_𝑛=0𝑝^𝑛𝑢_𝑛(𝑡−𝜏)^𝛽−1𝑑𝜏. (35)
Comparing the terms of the same power of 𝑝 yields the following two sets of linear equation: 𝑝⁰∶𝑢₀=𝑢𝑥,𝑦,0=𝑒^𝑥+𝑦,𝑝¹∶𝑢₁=1Γ(𝛼)^𝑡₀(𝑡−𝜏)^𝛼−1×𝜕_𝑥𝑥𝑢²₀+𝜕_𝑦,𝑦𝑢²₀+𝑎𝑢₀−𝑏𝑢₀𝑣₀𝑢𝑑𝜏,₁𝑝𝑥,𝑦,0=0,_𝑛:𝑢_𝑛=1𝑥,𝑦,𝑡Γ(𝛼)^𝑡₀(𝑡−𝜏)^𝛼−1×^𝑛−1_𝑗=0𝑢_𝑗𝑥𝑥𝑢_{𝑛−𝑗−1𝑥𝑥}+^𝑛−1_𝑗=0𝑢_{𝑗𝑦,𝑦}𝑢_{𝑛−𝑗−1𝑦𝑦}+𝑎𝑢_𝑛−1−𝑏^𝑛−1_𝑗=0𝑢_𝑗𝑣_{𝑛−𝑗−1}𝑢𝑑𝜏,_𝑛𝑝𝑥,𝑦,0=0,𝑛≥2,₀∶𝑣₀=𝑣𝑥,𝑦,0=𝑒^𝑥+𝑦,𝑝¹:𝑣₁=1Γ𝛽^𝑡₀(𝑡−𝜏)^𝛼−1×𝜕_𝑥,𝑥𝑣²₀+𝜕_𝑦,𝑦𝑣²₀+𝑏𝑢₀𝑣₀−𝑐𝑣₀𝑢𝑑𝜏,₁𝑝𝑥,𝑦,0=0,_𝑛:𝑣_𝑛=1𝑥,𝑦,𝑡Γ𝛽^𝑡₀^𝑛−1_𝑗=0𝑣_𝑗𝑥𝑥𝑣_{𝑛−𝑗−1𝑥𝑥}+^𝑛−1_𝑗=0𝑣_{𝑗𝑦,𝑦}𝑣_{𝑛−𝑗−1𝑦𝑦}+𝑏^𝑛−1_𝑗=0𝑢_𝑗𝑣_{𝑛−𝑗−1}−𝑐𝑣_𝑛−1,×(𝑡−𝜏)^𝛽−1𝑑𝜏,𝑣_𝑛𝑥,𝑦,0=0,𝑛≥2. (36)
The following solutions are obtained: 𝑢₀=𝑢𝑥,𝑦,0=𝑒^𝑥+𝑦,𝑣₀=𝑣𝑥,𝑦,0=𝑒^𝑥+𝑦,𝑢₁=𝑒^𝑥+𝑦𝑎−−8+𝑏𝑒^𝑥+𝑦𝑡^𝛼Γ(),𝑣1+𝛼₁=𝑒^𝑥+𝑦−𝑐+(8+𝑏)𝑒^𝑥+𝑦𝑡^𝛽Γ,𝑢1+𝛽₂=1Γ()ΓΓ()Γ×41+𝛼1+𝛽1/2+𝛽1+𝛼+𝛽_−𝛽𝑒^𝑥+𝑦𝑡^𝛽×𝑏𝑒^𝑥+𝑦4^𝛽𝑎+2^3+2𝛽−4^𝛽𝑏𝑒^𝑥+𝑦×𝑡^𝛼ΓΓ()Γ11+𝛽1+𝛼2+𝑐+𝛽₂−2(6+𝑏)𝑐𝑒^𝑥+𝑦+128+24𝑏+𝑏²√𝜋𝑡^𝛽Γ()11+𝛼×Γ2,𝑢+𝛽₂=1Γ()ΓΓ()Γ×41+𝛼1+𝛽1/2+𝛽1+𝛼+𝛽_−𝛼𝑒^𝑥+𝑦𝑡^𝛼×−𝑏𝑒^𝑥+𝑦−4^𝛼𝑐+2^3+2𝛽+4^𝛼𝑏𝑒^𝑥+𝑦×𝑡^𝛽ΓΓ()Γ11+𝛽1+𝛼2+𝑐+𝛽₂−2(−6+𝑏)𝑎𝑒^𝑥+𝑦+128−24𝑏+𝑏²√𝜋𝑡^𝛼Γ()11+𝛼×Γ2+𝛽. (37)
Using the package Mathematica, in the same manner, one can obtain the rest of the components. But in this case, 3 terms were computed and the asymptotic solution is given by 𝑢_𝑁=3(𝑥,𝑡)=𝑢₀(𝑥,𝑡)+𝑢₁(𝑥,𝑡)+𝑢₂()𝑣𝑥,𝑡+⋯,_𝑁=3(𝑥,𝑡)=𝑣₀(𝑥,𝑡)+𝑣₁(𝑥,𝑡)+𝑣₂(𝑥,𝑡)+⋯. (38)

Figures 4 and 5 show the numerical solutions for prey-predator population system with appropriate parameter. From the figures, we recognize that prey population density first increases with the spatial variables but will increase less than the predator population. However, the predator population density always increases with the spatial variables with the parameter we choose here. Analysis and results of prey-predator population system indicate that the fractional model matches the regular biological diffusion behavior observed in the field. From the figures, it is also clear to see the time evolution of prey-predator population density and we also know that the numerical solutions of fractional prey-predator population model are continuous with the parameters 𝛼 and 𝛽.

Figure 4 

Prey density for 𝑡=0.053,  𝛼=0.88,  𝛽=0.542,  𝑎=0.5,  and 𝑐=0.3.

Figure 5 

Predator density for 𝑡=0.053,  𝛼=0.88, 𝛽=0.542, 𝑎=0.5, and 𝑐=0.3.

To test the accuracy of the method used in this paper, we present the numerical result of the approximate solution of (27) and the numerical results of the exact solution in Table 1.

Table 1 

Comparison of the numerical results of approximate solutions obtained via HDM, VIM and HPM with exact solution of (27).

The numerical results show that the HPM and HDM are more accurate than the VIM in this case.

6. Conclusion

In this paper, we have productively developed HDM for solving partial differential equation of space- and time-fractional derivatives. the dependability of the HDM and its wider applicability to fractional differential equation. It is consequently, that the HDM makes available more practical series solutions that converge very speedily in real physical problems. Also, the amount of the computational effort has been abridged. In the bargain, the computations concerned in HDM are very straightforward. It is established that HDM is a prevailing and resourceful instrument for FPDEs.

Conflict of Interests

The authors declare that there is no conflict of interests for this paper.

Authors’ Contribution

The first draft of this paper was portrayed by Abdon Atangana; the revised version was carefully corrected by Samir Brahim Belhaouari, and both authors read and submitted the final version.

Acknowledgment

The authors will like to thank the reviewer for his valuable spare time to read this paper and for his valuable comments and suggestions toward the enhancement of this paper.