Approximate Analytical Solutions of Fractional Perturbed Diffusion Equation by Reduced Differential Transform Method and the Homotopy Perturbation Method

Abstract

The approximate analytical solutions of differential equations with fractional time derivative are obtained with the help of a general framework of the reduced differential transform method (RDTM) and the homotopy perturbation method (HPM). RDTM technique does not require any discretization, linearization, or small perturbations and therefore it reduces significantly the numerical computation. Comparing the methodology (RDTM) with some known technique (HPM) shows that the present approach is effective and powerful. The numerical calculations are carried out when the initial conditions in the form of periodic functions and the results are depicted through graphs. The two different cases have studied and proved that the method is extremely effective due to its simplistic approach and performance.

1. Introduction

Various fields of science and engineering deal with dynamical systems (see [1, 2]), some of which can be described by fractional-order equations (see [3–5] and the references therein). The last two decades have witnessed a great progress in fractional calculus and fractional-order dynamical systems. It has been found that fractional calculus is a mathematical tool that works adequately for anomalous social and physical systems with nonlocal, frequency- and history-dependent properties, and for intermediate states such as soft materials, which are neither idea solid nor idea fluid (see [3, 4, 6]). Differential equations with fractional-order derivatives/integrals are called fractional differential equations, and they have found many successful applications in viscoelasticity, heat conduction, electromagnetic wave, diffusion wave, control theory, and so on (see [7–9] and the references therein).

Bistable systems play an important role in the study of spatial patterns. A typical example, which appears in population dynamics, leads to 𝜕𝑢=𝜕𝜕𝑡²𝑢𝜕𝑥²+𝑢−𝑢³, (1) referred to as the classical Fisher-Kolmogorov (FK) equation. In investigating such spatial patterns, a key role is played by a model equation, which is simpler than the full equation describing the process. Recently, interest has turned to higher-order model equations involving bistable dynamics, such as the extended Fisher-Kolmogorov (EFK) equation 𝜕𝑢𝜕𝜕𝑡=−𝛾⁴𝑢𝜕𝑥⁴+𝜕²𝑢𝜕𝑥²+𝑢−𝑢³,𝛾>0, (2) proposed by Coullet et al. [10] in 1987 as well as by Dee and Van Saarloos [11] in 1988. The EFK equation has appeared in studies of phase transitions, for instance, near a Lifshitz point (cf. Zimmerman [12]). Another well-known equation of this type is the Swift-Hohenberg equation 𝜕𝑢𝜕𝜕𝑡=𝛼𝑢−1+²𝜕𝑥²²𝑢−𝑢³,𝛼>0, (3) proposed in 1977.

Whereas the EFK equation and the SH equation involve fourth-order spatial derivatives, certain phase-field models lead to even-higher-order spatial gradients. The following equation 𝜕𝑢=𝜕𝜕𝑡⁶𝑢𝜕𝑥⁶𝜕+𝐴⁴𝑢𝜕𝑥⁴𝜕+𝐵²𝑢𝜕𝑥²+𝑢−𝑢³, (4) studied by Gardner and Jones [13] and Caginalp and Fife [14], in which 𝐴 and 𝐵 are constants. In this paper, we consider the initial value problem (IVP) of perturbed diffusion equation with fractional time derivative 𝐷^𝛼_𝑡𝜕𝑢=⁶𝑢𝜕𝑥⁶𝜕+𝐴⁴𝑢𝜕𝑥⁴𝜕+𝐵²𝑢𝜕𝑥²+𝑢−𝑢²,𝑢(𝑥,0)=𝜑(𝑥), (5) where the fractional derivative in (5) is in the sense of Caputo and 0≤𝛼<1.

Historically, comparatively little was known about the extraordinary range of behavior presented by the solutions of nonlinear partial differential equations. Except for a few special cases, it is impossible to find a closed form solution for a fractional differential equation. Effective methods for solving such equations are needed. So approximate and numerical techniques must be used.

The homotopy perturbation method (HPM) is relatively new approach to provide an analytical approximation to linear or nonlinear problem. This method was first presented by He [15, 16] and applied to various nonlinear problems [17, 18]. The basic difference of this method from the other perturbation techniques is that it does not require small parameters in the equation, which overcomes the limitations of the traditional perturbation techniques. Recently, the application of this method is extended for fractional differential equations [5, 19–21]. Zhou [22] was the first one to use differential transform method (DTM) in engineering applications. He employed DTM in solution of initial boundary value problems in electric circuit analysis. This method constructs an analytical solution in the form of a polynomial. It is different from the traditional high order Taylor series method, which requires symbolic computation of necessary derivatives of the data functions. The differential transform is an iterative procedure for obtaining analytical differential equations. The concept of DTM has broadened to problems involving partial differential equations and systems of differential equations [23, 24]. Some researchers have lately applied DTM for analysis of uniform and non-uniform beams [25–27]. Recently, the application of reduced differential transform method is successfully extended to obtain analytical approximate solutions to linear and nonlinear ordinary differential equations of fractional order [28–30].

In this paper, we try to find an approximate analytical solution of (5) for fractional time derivatives with the help of powerful analytical method. We use the reduced differential transform method (RDTM) and homotopy perturbation method (HPM) to obtain the solutions and compare them with each other. We know that the HPM method is based on the use of homotopy parameter for classification of most favorable values of parameters in between [0,1]. While, RDTM technique does not require any parameter, discretization, linearization, or small perturbations and therefore it reduces significantly the numerical computation. For the standard cases, comparing the methodology with some known techniques shows that the present approach is effective and powerful.

2. Preliminaries

There are several approaches to define the fractional calculus, for example, Riemann-Liouville, Gruunwald-Letnikow, Caputo, and Generalized Functions approach. For the readers’ convenience, definitions of fractional integral/derivative and some preliminary results are given in this section.

Definition 1 (see [6]). The fractional integral of order 𝛼≥0 of a function 𝑢(𝑥,𝑡) is given by 𝐼^𝛼𝑢()=1𝑥,𝑡Γ(𝛼)^𝑡₀(𝑡−𝑠)^𝛼−1𝑢(𝑥,𝑠)𝑑𝑠, (6) provided that the right side is point-wise defined on (0,+∞), where Γ(⋅) is the well-known gamma function.

Definition 2 (see [6]). The Caputo derivative of order 𝑠>0 of a continuous function 𝑢(𝑥,𝑡) is defined to be ^𝐶𝐷^𝑠_𝑡𝑢()=1𝑥,𝑡Γ(𝑛−𝑠)^𝑡₀𝑢^(𝑛)(𝑥,𝜏)(𝑡−𝜏)^{𝑠−𝑛+1}𝑑𝜏, (7) where 𝑛=[𝑠]+1, provided that the right side is point-wise defined on (0,+∞).

Definition 3 (see [6]). The Riemann-Liouville derivative of order 𝑠>0 of a continuous function 𝑢(𝑥,𝑡) is defined to be ^𝑅𝐿𝐷^𝑠_𝑡𝑢()=1𝑥,𝑡Γ()𝑑𝑛−𝑠𝑑𝑡^𝑛𝑛𝑡₀𝑢(𝑥,𝜏)(𝑡−𝜏)^{𝑠−𝑛+1}𝑑𝜏, (8) where 𝑛=[𝑠]+1, provided that the right side is point-wise defined on (0,+∞).

Riemann-Liouville fractional derivative is mostly used by mathematicians but this approach is not suitable for real world physical problems since it requires the definition of fractional order initial conditions, which have no physically meaningful explanation yet. Caputo introduced an alternative definition, which has the advantage of defining integer order initial conditions for fractional order differential equations. We have chosen the Caputo fractional derivative because it allows traditional initial and boundary conditions to be included in the formulation of the problem. And some other properties of fractional derivative can be found in [4, 6].

The relation between the Riemann-Liouville operator and Caputo operator is given by ^𝐶𝐷^𝑠_𝑡𝑢(𝑥,𝑡)=^𝑅𝐿𝐷^𝑠_𝑡𝑢(𝑥,𝑡)−^𝑛−1_𝑖=1𝑢^(𝑖)()𝑡𝑥,0^𝑖.𝑖! (9) By using (9), from (8), we obtain fractional derivative in the Caputo sense as follows: ^𝐶𝐷^𝑠_𝑡𝑢()=1𝑥,𝑡Γ()𝑑𝑛−𝑠𝑑𝑡^𝑛𝑛×^𝑡₀𝜔(𝑥,𝜏)−∑^𝑛−1_𝑖=1𝑢^(𝑖)(𝑥,0)𝜏^𝑖/𝑖!(𝑡−𝜏)^{𝑠−𝑛+1}𝑑𝜏. (10)

Let us expand the analytical and continuous function 𝑢(𝑥,𝑡) in terms of a fractional power series as follows: 𝑢(𝑥,𝑡)=^∞_𝑘=0𝑈_𝑘(𝑥)𝑡^𝑘/𝛼, (11) where 𝛼 is the order of fraction and 𝑈_𝑘(𝑥) is the fractional differential transform of 𝑢(𝑥,𝑡).

Since the initial conditions are implemented to the integer order derivatives, the transformations of the initial conditions are defined as follows: 𝑈_𝑘(𝑥)=𝑘if𝛼∈𝑍⁺,1()!×𝑑𝑘/𝛼_𝑘/𝛼𝑑𝑡^𝑘/𝛼𝑢(𝑥,𝑡)_𝑡=0()𝑘for𝑘=0,1,2,…,𝛼𝑠−1if𝛼∈𝑅⧵𝑍⁺,0. (12)

The following theorems that can be deduced from (7) and (8) are given below.

Theorem 4. If 𝑤(𝑥,𝑡)=𝑥^𝑚𝑡^𝑛𝑢(𝑥,𝑡), then 𝑊_𝐾(𝑥)=𝑥^𝑚𝑈_𝑘−𝑛(𝑥).

Theorem 5. If 𝑤(𝑥,𝑡)=𝑢(𝑥,𝑡)𝑣(𝑥,𝑡), then 𝑊_𝐾(𝑥)=^𝑘_𝑟=0𝑈_𝑟(𝑥)𝑉_𝑘−𝑟(𝑥)=^𝑘_𝑟=0𝑈_𝑘−𝑟(𝑥)𝑉_𝑟(𝑥). (13)

Theorem 6. If 𝑤(𝑥,𝑡)=(𝜕^𝑚/𝜕𝑥^𝑚)𝑢(𝑥,𝑡), then 𝑊_𝐾(𝑥)=(𝜕^𝑚/𝜕𝑥^𝑚)𝑈_𝑘(𝑥).

Theorem 7. If 𝑤(𝑥,𝑡)=(𝜕^𝛽/𝜕𝑡^𝛽)𝑢(𝑥,𝑡), then 𝑊_𝐾(𝑥)=(Γ(𝛽+1+𝑘/𝛼)/Γ(1+𝑘/𝛼))𝑈_{𝑘+𝛼𝛽}(𝑥).

3. Solution of the Problem by the RDTM and HPM

To illustrate the basic ideas of HPM [15], we consider the following general nonlinear differential equation: 𝐴𝑦−𝑓(𝑟)=0,𝑟∈Ω, (14) with boundary conditions 𝐵𝑦,𝜕𝑦𝜕𝑛=0,𝑟∈Γ, (15) where 𝐴 is a general differential operator, 𝐵 is a boundary operator, 𝑓(𝑟) is a known analytic function, and Γ is the boundary of the domain Ω.

The operator 𝐴 can be generally divided into two parts 𝐿 and 𝑁, where 𝐿 is linear, while 𝑁 is nonlinear. Therefore (14) can be written as follows: 𝐿𝑦+𝑁𝑦−𝑓(𝑟)=0. (16) We construct a homotopy of (14) 𝑦(𝑟,𝑝)∶Ω×[0,1]→ℜ which satisfies 𝐻𝑦,𝑝=1−𝑝𝐿𝑦−𝐿𝑦₀𝐴𝑦(𝑟[]+𝑝−𝑓)=0,𝑝∈0,1,𝑟∈Ω, (17) which is equivalent to 𝐻𝑦,𝑝=𝐿𝑦−𝐿𝑦₀+𝑝𝐿𝑦₀+𝑝𝑁𝑦−𝑓(𝑟)=0, (18) where 𝑝∈[0,1] is an embedding parameter, and 𝑦₀ is an initial guess approximation of (14) which satisfies the boundary conditions. It follows from (17) and (18) that 𝐻𝑦,0=𝐿𝑦−𝐿𝑦₀𝐻𝑦(𝑟)=0,𝑦,1=𝐴−𝑓=0. (19) Thus, the changing process of 𝑝 from 0 to 1 is just that of 𝑦(𝑟,𝑝) from 𝑦₀(𝑟) to 𝑦(𝑟). In topology this is called deformation and 𝐿(𝑦)−𝐿(𝑦₀) and 𝐴(𝑦)−𝑓(𝑟) are called homotopic. Here the embedding parameter is introduced much more naturally, unaffected by artificial factors; further it can be considered as a small parameter for 0≤𝑝≤1. So it is very natural to assume that the solution of (18) and (19) can be expressed as 𝑦(𝑡)=𝑦₀(𝑡)+𝑝𝑦₁(𝑡)+𝑝²𝑦₂(𝑡)+⋯. (20) The approximate solutions of the original equations can be obtained by setting 𝑝=1; that is, 𝑦(𝑡)=lim^∞_{𝑝→1𝑛=0}𝑝^𝑛𝑦_𝑛(𝑡)=𝑦₀(𝑡)+𝑦₁(𝑡)+𝑦₂(𝑡)+⋯. (21) The convergence of series (21) has been proved by He in his paper [15].

Because of the knowledge of various perturbation methods that low-order approximate solution leads to high accuracy, there requires no infinite series. Then after a series of recurrent calculation by using Mathematica software, we will get approximate solutions of fractional model (5).

According to HPM, we construct the following homotopy of (5) 𝐷^𝛼_𝑡𝜕𝑢=𝑝⁶𝑢𝜕𝑥⁶𝜕+𝐴⁴𝑢𝜕𝑥⁴𝜕+𝐵²𝑢𝜕𝑥²+𝑢−𝑢², (22) where the homotopy parameter 𝑝 is considered to be small, 0≤𝑝≤1.

Assuming the solution of (22) has the form 𝑢(𝑥,𝑡)=𝑢₀+𝑝𝑢₁+𝑝²𝑢₂+𝑝³𝑢₃+⋯, (23) when 𝑝→1, (23) becomes the approximate solution of (5). Substituting (23) into (22) and equating the terms with identical powers of 𝑝, we obtain the following set of linear differential equations: 𝑝⁰∶𝐷^𝛼_𝑡𝑢₀𝑝=0,₁∶𝐷^𝛼_𝑡𝑢₁=𝐷_{𝑥𝑥𝑥𝑥𝑥𝑥}𝑢₀+𝐴𝐷_{𝑥𝑥𝑥𝑥}𝑢₀+𝐵𝐷_𝑥𝑥𝑢₀+𝑢₀−𝑢²₀,𝑝²∶𝐷^𝛼_𝑡𝑢₂=𝐷_{𝑥𝑥𝑥𝑥𝑥𝑥}𝑢₁+𝐴𝐷_{𝑥𝑥𝑥𝑥}𝑢₁+𝐵𝐷_𝑥𝑥𝑢₁+𝑢₁−2𝑢₀𝑢₁, (24) and so on. The method is based on applying the operators 𝐼^𝛼_𝑡 (the inverse operators of the Caputo derivative 𝐷^𝛼_𝑡) on both sides of the above linear differential equations. Using this selection we obtain the successive approximations and the solution may be obtained.

Finally, we approximate the analytical solutions of 𝑢(𝑥,𝑡) by the truncated series 𝑢(𝑥,𝑡)=lim_𝑁→∞^𝑁−1_𝑛=0𝑢_𝑛(𝑥,𝑡). (25)

For the illustration of the methodology of the reduced differential transform method, we write (5) in the standard operator form 𝐿(𝑢(𝑥,𝑡))−𝑅(𝑢(𝑥,𝑡))+𝑁(𝑢(𝑥,𝑡))=0, (26) with initial condition 𝑢(𝑥,0)=𝜑(𝑥), (27) where 𝐿(𝑢(𝑥,𝑡))=(𝜕^𝛼/𝜕𝑡^𝛼)𝑢(𝑥,𝑡) is the fractional time derivative operator and 𝑅(𝑢(𝑥,𝑡))=𝜕⁶𝑢/𝜕𝑥⁶+𝐴(𝜕⁴𝑢/𝜕𝑥⁴)+𝐵(𝜕²𝑢/𝜕𝑥²)+𝑢 is linear operator and 𝑁(𝑢(𝑥,𝑡))=𝑢² is nonlinear operator.

According to the RDTM and Theorems 6 and 7, we can construct the following iteration formula for (26): Γ𝛼+1+𝑘/𝑞Γ𝑈1+𝑘/𝑞^{𝑘+𝛼𝑞}(𝑥)=𝑅𝑈_𝑘(𝑥)−𝑁𝑈_𝑘(𝑥). (28) For the convenience of readers, we can give the first few nonlinear term are 𝑁𝑈₀(𝑥)=𝑈²₀(𝑥),𝑁𝑈₁(𝑥)=2𝑈₀(𝑥)𝑈₁(𝑥),𝑁𝑈₂(𝑥)=2𝑈₀(𝑥)𝑈₂(𝑥)+𝑈²₁(𝑥). (29) From initial condition, we write 𝑈₀(𝑥)=𝜑(𝑥). (30) Substituting (30) into (28) and by a straight forward iterative calculations, we get the following 𝑈_𝑘(𝑥) values: 𝑈₀()(𝑥),𝑈𝑥,𝑡=𝜑₁(𝑥,𝑡)=𝜑⁽⁶⁾(𝑥)+𝐴𝜑⁽⁴⁾(𝑥)+𝐵𝜑⁽²⁾(𝑥)+𝜑(𝑥)−𝜑(𝑥)²×𝑡^𝛼Γ(),𝛼+1 (31) and so on. Then the inverse transformation of the set of values {𝑈_𝑘(𝑥)}^𝑛_𝑘=0 gives approximation solution as 𝑢_𝑛(𝑥,𝑡)=^𝑛_𝑘=0𝑈_𝑘(𝑥)𝑡^𝑘, (32) where 𝑛 is order of approximation solution.

Therefore, the exact solution of problem is given by 𝑢(𝑥,𝑡)=lim_𝑛→∞𝑢_𝑛(𝑥,𝑡). (33)

4. Numerical Results and Discussion

Here, we take three different values of 𝜑(𝑥) and comparing the results of RDTM and HPM in the form of three- and two dimensional figures for each case, we would see that RDTM and HPM solutions are in excellent agreement. In this section we assume that 𝐴=𝐵=1.

Case Study 1. 𝜑(𝑥)=4sech(𝜇𝑥). See Figures 1 and 2.

(a)

(b)

Figure 1 

The comparison of the results of the RDTM and HPM for 𝛼=1/4 and 𝛼=3/4 at 𝑡=0.8 and 𝜇=1.

(a)

(b)

(c)

(d)

Figure 2 

Results of the HPM for 𝛼=1/4, 1/2, 3/4 and 1 and 𝜇=1.

Case Study 2. 𝜑(𝑥)=𝜋+𝜇cos(𝛽𝑥). See Figures 3 and 4.

(a)

(b)

Figure 3 

The comparison of the results of the RDTM and HPM for 𝛼=1/4 and 𝛼=3/4 at 𝑡=0.8 and 𝛽=𝜇=1.

(a)

(b)

(c)

(d)

Figure 4 

Results of the HPM for 𝛼=1/4, 1/2, 3/4 and 1 and 𝛽=𝜇=1.

5. Conclusions

In this paper, RDTM and HPM are successfully applied to find the solution of the fractional differential equation and standard motion with different initial conditions. Unlike the traditional methods, the solutions here are given in series 9 form. The approximate solution to the equation was computed without any need for special transformations, linearization, or discretization. In addition, we compare these two methods and show that the results of the RDTM method are in excellent agreement with results of the HPM method and the obtained numerical solutions are shown graphically. We use the mathematical simulation to calculate the functions obtained from the RDTM and HPM. It was shown that RDTM and HPM methods are powerful tools for solving analytically of nonlinear equations.

Acknowledgments

This project was supported by the National Natural Science Foundation of China (no. 41206163), Open Research Program of State Key Laboratory of Tropical Oceanography (South China Sea Institute of Oceanology Chinese Academy of Sciences) (no. LTO1205), and the Pre-research Foundation of PLA University of Science and Technology.