# Series Solutions of Time-Fractional PDEs by Homotopy Analysis Method

Authors: O. Abdulaziz; I. Hashim; A. Saif
Journal: Differential Equations and Nonlinear Mechanics
 (2008)

---

## Abstract

The homotopy analysis method (HAM) is applied to solve linear and nonlinear fractional partial differential equations (fPDEs). The fractional derivatives are described by Caputo's sense. Series solutions of the fPDEs are obtained. A convergence theorem for the series solution is also given. The test examples, which include a variable coefficient, inhomogeneous and hyperbolic-type equations, demonstrate the capability of HAM for nonlinear fPDEs.


---

## Body

## 1. Introduction

Fractional
calculus has been given considerable popularity and importance during the past
three decades, due mainly to its   applications in numerous fields of science and
engineering. For example, phenomena in the areas of fluid flow, rheology,
electrical networks, probability and statistics, control theory of dynamical
systems, electrochemistry of corrosion, chemical physics, optics and signal
processing, and so on can be successfully modelled by linear or nonlinear
fractional differential equations (fDEs) [1–4].Finding accurate methods for solving nonlinear
differential equations has become important. Some of the analytical methods for
nonlinear differential equations are the Adomian decomposition method (ADM)
[5–14], the homotopy-perturbation method (HPM) [15–19], variational iteration method (VIM) [12, 20–24], and the EXP-function
method [25]. Another
analytical approach that can be applied to solve nonlinear differential
equations is to employ the homotopy analysis method (HAM) [26–29]. Some of the recent
applications of HAM can be found in [30–41]. An account of the recent
developments of HAM was given by Liao [42]. HAM has been successfully applied into engineering
fields. The method has been applied to give an explicit solution for the
Riemann problem of the nonlinear shallow-water equations [43]. The obtained Riemann
solver has been implemented into a numerical model to simulate long waves, such
as storm surge or tsunami, propagation and run-up.Very recently, Song and Zhang [44] applied HAM to solve
fractional KdV-Burgers-Kuramoto equation. Cang
et al.  [45] solved nonlinear Riccati differential equations of fractional order using HAM. Hashim et al.
[46] employed
HAM to solve fractional initial value problems (fIVPs) for ordinary
differential equations. In [47],
the applicability of the HAM was extended to construct numerical solution for
the fractional BBM-Burgers equation. The HAM solutions for systems of nonlinear fractional differential equations were presented by
Bataineh et al. [48].A specific linear, nonhomogeneous time fractional
partial differential equation (fPDE) with variable coefficients was first
transformed to two fractional ordinary differential equations which were then
solved by HAM in [49].
Recently, Xu  et al.  [50] applied the HAM to linear, homogeneous one- and
two-dimensional fractional heat-like PDEs subject to the Neumann boundary
conditions. Jafari and Seifi [51] applied HAM to linear and nonlinear homogeneous
fractional diffusion-wave equations. Very recently, the HAM was shown to be
capable of solving linear and nonlinear systems of fPDEs [52].In this paper, we shall consider linear and nonlinear
fPDEs of the form(1.1)Dtαu(x,t)=f(u,ux,uxx),n−1<α≤n,⁢t>0,subject to the initial
conditions(1.2)u(k)(x,0)=gk(x),k=0,1,2,…,n−1,where n is an integer, f is a linear/nonlinear function, and Dtα(⋅)=∂α(⋅)/∂tα is a fractional differential operator. We
shall demonstrate the applicability of HAM to fPDEs through several linear and
nonlinear test examples.


## 2. Preliminaries

The fractional derivative is defined in the Caputo
sense as in [53],(2.1)Dαw(t)=Jn−αDnw(t).Here Dn is the usual integer differential operator of
order m and Jβ is the Riemann-Liouville fractional integral
operator of order β>0,
defined by(2.2)Jβw(t)=1Γ(β)∫0t(t−τ)β−1w(τ)dτ,(t>0),J0w(t)=w(t). Some of the
properties of the operator Jβ,
which we will need in our work, are as follows [2, 3]:(1)
JβJμw(t)=Jβ+μw(t),(2)
JβJμw(t)=JμJβw(t),(3)
Jβtγ=(Γ(γ+1)/Γ(β+γ+1))tβ+γ. Caputo’s
fractional derivative has a useful property [54](2.3)(JμDμ)w(t)=w(t)−∑k=0n−1w(k)(0+)tkk!,(n−1<μ≤n).The operator form of the
nonlinear fPDEs (1.1) can be written as follows:(2.4)Dtα(u(x,t))=A(u,ux,uxx)+B(u,ux,uxx)+C(x,t),n−1<α≤n,⁢t>0,subject to the initial
conditions(2.5)u(k)(x,0)=gk(x),k=0,1,2,…,n−1,where A is a linear operator which might include other
fractional derivatives of order less that α, B is a nonlinear operator which also might
include other fractional derivatives of order less that α and C is a known analytic function.Applying the operatorJα,
the inverse operator of Dα,
to both sides of (2.4) with considering the initial conditions (2.5) according to
(2.3), we obtain(2.6)u(x,t)=∑k=0n−1gk(x)tkk!+JαC(x,t)+JαA(u,ux,uxx)+JαB(u,ux,uxx),n−1<α≤n,⁢t>0.


## 3. Homotopy Analysis Method (HAM)

## 3.1. The Zeroth-Order Deformation Equation

Letℒ denote an auxiliary linear operator, u0(x,t) is an initial approximation of u(x,t) which satisfies the initial conditions (2.5).
Note that, in this paper, the auxiliary linear operator ℒ is not the same linear operator A of (2.4).Note that the original equation (1.1) contains the
linear operator Dtα.
So, it is straightforward for us to choose the auxiliary linear
operator(3.1)ℒ(ϕ)=Dtα(ϕ).According to (2.6), we can choose
the initial approximation to be(3.2)u0(x,t)=∑k=0m−1gk(x)tkk!+JαC(x,t).For simplicity, let us define, according to (2.4), the
nonlinear operator(3.3)N(ϕ)=Dtα(ϕ)−A(ϕ,ϕx,ϕxx)−B(ϕ,ϕx,ϕxx)−C(x,t).Hence, in the frame of HAM
[29], we can construct
the so-called zeroth-order deformation(3.4)(1−q)ℒ(U(x,t;q)−u0(x,t))=qℏN(U(x,t;q)),subject to the following initial
conditions:(3.5)U(k)(x,0;q)=gk(x),k=0,1,2,…,m−1,where q∈[0,1] is the embedding parameter, ℏ≠0 is an auxiliary parameter, and U(x,t;q) is an unknown function on the independent
variables x,t,
and q.Whenq=0,
since u0(x,t) satisfies all the initial conditions (2.5), and ϕ=0 is a solution of ℒϕ=0,
we have obviously(3.6)U(x,t;0)=u0(x,t),and when q=1,
the zeroth-order deformation equations (3.4) and (3.5) are equivalent to the
original equations (2.4) and (2.5), provided(3.7)U(x,t;1)=u(x,t).Using the parameter q,
we expand U(x,t;q) in Taylor series as follows:(3.8)U(x,t;q)=u0(x,t)+∑m=1∞um(x,t)qm,where(3.9)um(x,t)=1m!∂mU(x,t;q)∂mq|q=0.Assume that the auxiliary linear operatorℒ,
the initial guess u0 and the auxiliary parameter ℏ are properly chosen such that the series (3.8)
is convergent at q=1.
Thus, due to (3.7) we have(3.10)u(x,t)=u0(x,t)+∑m=1∞um(x,t).


## 3.2. Themth-Order Deformation Equation

Let us define the vector(3.11)u→n={u0(x,t),u1(x,t),…,un(x,t)}.Following Liao [26–29], differentiating (3.4) m times with respect to the embedding parameter q,
then setting q=0,
and finally dividing them by m!,
we have the so-called mth-order deformation equation(3.12)ℒ[um(x,t)−χmum−1(x,t)]=ℏℛm(u→m−1),subject to the initial
conditions(3.13)um(k)(x,0)=0,k=0,1,2,…,m−1,where(3.14)ℛm(u→m−1)=1(m−1)!∂m−1N(U(x,t;q))∂m−1q|q=0,(3.15)χm={0,m≤1,1,m>1.Substituting (3.3) into (3.14), and
since A is a linear operator, ℛm(u→m−1) can be given by(3.16)ℛm(u→m−1)=Dtαum−1−A(u(m−1),u(m−1)x,u(m−1)xx)−1(m−1)!∂m−1B(U,Ux,Uxx)∂m−1q|q=0−(1−χm)C(x,t).According to (3.1), we can apply the operator Jα to both sides of (3.12) to
obtain(3.17)JαDα[um(x,t)−χmum−1(x,t)]=ℏJα[ℛm(u→m−1)].Using the property (2.3) and the
initial conditions (1.1), we have(3.18)um(x,t)=χmum−1(x,t)+ℏJα[ℛm(u→m−1)].Finally, for the purpose of computation, we will
approximate the HAM solution (3.10) by the following truncated
series:(3.19)ϕm(x,t)=∑k=0m−1uk(x,t).


## 3.3. Convergence Theorem

Theorem 3.1.
As long as the seriesu(x,t)=u0(x,t)+∑m=1+∞um(x,t) converges, where um(x,t) is governed by (3.12) under the definitions (3.14)
and (3.15), it must be a solution of (2.4).Proof.
If the series∑m=0+∞um(x,t) is convergent, we can write(3.20)S(x,t)=∑m=0+∞um(x,t),and it holds(3.21)limn→+∞un(x,t)=0.From (3.12) and by using (3.15),
it yields(3.22)ℏ∑m=1+∞ℛm(u→m−1)=∑m=1+∞ℒ[um(x,t)−χmum−1(x,t)]=limn→+∞∑m=1nℒ[um(x,t)−χmum−1(x,t)]=ℒ[limn→+∞∑m=1n(um(x,t)−χmum−1(x,t))]=ℒ[limn→+∞un(x,t)]=0.Since ℏ≠0,
then(3.23)∑m=1+∞ℛm(u→m−1)=0.Substituting (3.16) into the above
equation and simplifying it, due to the convergence of the series u(x,t)=u0(x,t)+∑m=1+∞um(x,t) and since A is a linear operator,
yield(3.24)∑m=1+∞ℛm(u→m−1)=∑m=1+∞[Dtαum−1−A(u(m−1),u(m−1)x,u(m−1)xx)]−∑m=1+∞[1(m−1)!∂m−1B(U,Ux,Uxx)∂m−1q|q=0−(1−χm)C(x,t)]=Dtα(∑m=0+∞um)−A(∑m=0+∞um,∑m=0+∞umx,∑m=0+∞umxx)−∑m=0+∞[1m!∂mB(U,Ux,Uxx)∂mq|q=0]−C(x,t).Now, expanding the nonlinear
term B(U(x,t;q),Ux(x,t;q),Uxx(x,t;q)) by using the general Taylor theorem at q=0 yields(3.25)B(U(x,t;q),Ux(x,t;q),Uxx(x,t;q))=∑m=0+∞[1m!∂mB(U,Ux,Uxx)∂mq|q=0qm].Setting q=1 in the above equation and using (3.8), we
obtain(3.26)B(∑m=0+∞um,∑m=0+∞umx,∑m=0+∞umxx)=∑m=0+∞[1m!∂mB(U,Ux,Uxx)∂mq|q=0].Then(3.27)∑m=1+∞ℛm(u→m−1)=Dtα(∑m=0+∞um)−A(∑m=0+∞um,∑m=0+∞umx,∑m=0+∞umxx)−B(∑m=0+∞um,∑m=0+∞umx,∑m=0+∞umxx)−C(x,t)=Dtα(S(x,t))−A(S(x,t),Sx(x,t),Sxx(x,t))−B(S(x,t),Sx(x,t),Sxx(x,t))−C(x,t).From the initial conditions (3.5)
and (3.13), it holds that(3.28)S(k)(x,0)=∑m=0+∞um(k)(x,0)=u0(k)(x,0)+∑m=1+∞um(k)(x,0)=gk(x).
Thus,S(x,t) is satisfied and
also must be the exact solution for (2.4).


## 4. Test Examples

In this section, we shall illustrate the applicability
of HAM to several linear and nonlinear fPDEs.


## 4.1. Problem 1

Let us consider the following linear time-fractional
wave-like equations:(4.1)Dtαu(x,t)=12x2uxx(x,t),t>0,x∈R,⁢1<α≤2,(4.2)u(x,0)=x,ut(x,0)=x2.We note that the heat-like
counterpart of (4.1) was solved by HAM in [50] without direct comparison with the result by the ADM.
According to (3.2), we can choose the initial guess to be(4.3)u0(x,t)=x+x2t.From (3.18), we
have(4.4)um=(χm+ℏ)um−1−ℏ2x2Jα[(um−1)xx].Consequently, the first few
terms of HAM series solutions are as follows:(4.5)u1(x,t)=−ℏtα+1Γ(α+2)x2,u2(x,t)=−ℏ(ℏ+1)tα+1Γ(α+2)x2+ℏ2t2α+1Γ(2α+2)x2,u3(x,t)=−ℏ(ℏ+1)2tα+1Γ(α+2)x2+2ℏ2(ℏ+1)t2α+1Γ(2α+2)x2−ℏ3t3α+1Γ(3α+2)x2,and so on. Hence, the HAM series
solution is(4.6)u(x,t)=u0+u1+u2+u3+⋯=x+x2t−ℏ[1+(ℏ+1)+(ℏ+1)2+⋯]tα+1Γ(α+2)x2+ℏ2[1+2(ℏ+1)+⋯]t2α+1Γ(2α+2)x2−ℏ3t3α+1Γ(3α+2)x2+⋯.Since we choose the initial
guess u0(x,t) to be the same initial guess used by ADM
[12], we can notice
that when ℏ=−1,
the above expression gives the same solution given by ADM. Table 1 shows the
HAM approximation solutions for (4.1)-(4.2) when α=1.5, 1.75,
and 2 with ℏ=−1 and −1.0453.
It is to be noted that the first four terms of the HAM series were used to
evaluate the approximate solutions in Table 1.Table 1
Approximate
solution of (4.1) for some values of α using the 4-term HAM approximation, ϕ4,
with ℏ=−1 and ℏ=−1.0453.
txα=1.5α=1.75α=2.0ℏ=−1ℏ=−1.045ℏ=−1ℏ=−1.045ℏ=−1ℏ=−1.045Exact0.20.250.262840620.262840620.262669890.262669910.262583500.262583500.262583500.500.551362460.551362500.550679590.550679640.550334000.550334020.550334000.750.865565540.865565620.864029090.864029190.863251500.863251560.863251501.01.205449851.205449991.202718391.202718561.201336001.201336101.201336000.40.250.276971140.276971090.276156680.276156690.275672020.275672050.275672020.500.607884550.607884370.604626750.604626790.602688080.602688200.602688080.750.992740250.992739830.985410190.985410270.981048180.981048460.981048181.01.431538211.431537481.418507021.418507151.410752321.410752821.410752320.60.250.293094810.293095010.291090090.291089970.289790840.289790840.289790840.500.672379240.672380060.664360390.664359900.659163380.659163390.659163390.751.137853281.137855131.119810881.119809781.108117621.108117641.108117641.01.689516941.689520241.657441561.657439611.636653551.636653581.63665358


## 4.2. Problem 2

In this example, we consider the following
one-dimensional linear inhomogeneous time-fractional equation:(4.7)Dαu(x,t)+xux(x,t)+uxx(x,t)=2t+2x2+2,t>0,x∈R,⁢0<α≤1,subject to the initial
condition(4.8)u(x,0)=x2.In Section 3, we chose the
initial guess to contain the initial conditions and the source term C(x,t).
In this example, due to the appearance of noise terms and also to get the exact
solution, we will modify the way we choose the initial guess. The initial guess
is set to contain only the initial condition (4.8),
and the source term, C(x,t)=2t+2x2+2,
will be added to u1(x,t).
The other terms are obtained the same as described in Section 3.Hence, the initial guess is given by(4.9)u0(x,t)=x2,and according to (3.18), we
have(4.10)u1=Jα[2t+2x2+2]+ℏJα[x(u0)x+(u0)xx],um=(χm+ℏ)um−1+ℏJα[x(um−1)xx+(um−1)xx],m≥2.The terms of the HAM solution
series can be given by(4.11)u1(x,t)=2tα+1Γ(α+2)+2(ℏ+1)(x2+1)tαΓ(α+1),u2(x,t)=2(ℏ+1)tα+1Γ(α+2)+2(ℏ+1)(x2+1)[(ℏ+1)tαΓ(α+1)+2ℏt2αΓ(2α+1)],and so on. Hence, the HAM series
solution is(4.12)u(x,t)=u0+u1+u2+⋯=x2+2[1+(ℏ+1)+⋯]tα+1Γ(α+2)+2(ℏ+1)(x2+1)[1+(ℏ+1)+⋯]tαΓ(α+1)+4ℏ(ℏ+1)(x2+1)t2α+1Γ(2α+1)x2+⋯.Taking ℏ=−1 in (4.12), we obtain the exact
solution,(4.13)u(x,t)=x2+2tαΓ(α+2).


## 4.3. Problem 3

Consider the following nonlinear time-fractional
hyperbolic equation:(4.14)Dtαu(x,t)=∂∂x(u(x,t)∂u(x,t)∂x),t>0,x∈R,⁢1<α≤2,subject to the initial
conditions(4.15)u(x,0)=x2,ut(x,0)=−2x2.Equation (4.14) can be rewritten
as follows:(4.16)Dtαu(x,t)=(∂u(x,t)∂x)2+u(x,t)∂2u(x,t)∂x2.From (3.4), construct the following zeroth-order
deformation:(4.17)(1−q)ℒ(U(x,t;q)−u0(x,t))=qℏN(U(x,t;q)),subject to the following initial
conditions:(4.18)U(x,0)=x2,Ut(x,0)=−2x2,where(4.19)N(θ)=Dtαθ−θx2−θθxx.The auxiliary linear operator
can be chosen as follows:(4.20)ℒ(θ)=Dtα(θ),with the
property(4.21)ℒ(θ)=0whenθ=0,while, the initial guess
is(4.22)u0(x,t)=x2−2x2t.Again from (3.12), the high-order
deformation equation can be given by(4.23)ℒ[um(x,t)−χmum−1(x,t)]=ℏℛm(u→m−1),subject to the initial
conditions(4.24)um(k)(x,0)=0,k=0,1,where(4.25)ℛm(u→m−1)=1(m−1)!∂m−1N(U(x,t;q))∂m−1q|q=0.Then, ℛm(u→m−1) can be given by(4.26)ℛm(u→m−1)=Dtαum−1−∑i=0m−1uixu(m−1−i)x−∑i=0m−1uiu(m−1−i)xx.Accordingly, the governing
equation is as follows:(4.27)um=χmum−1+ℏJα[Dtαum−1−∑i=0m−1uixu(m−1−i)x]−ℏJα[∑i=0m−1uiu(m−1−i)xx],m≥1.Consequently, the first few terms of HAM series solutions
are given by(4.28)u1(x,t)=−6ℏx2[tαΓ(α+1)−4tα+1Γ(α+2)+8tα+2Γ(α+3)],u2(x,t)=−6ℏx2[(ℏ+1)tαΓ(α+1)−4(ℏ+1)tα+1Γ(α+2)+8(ℏ+1)tα+2Γ(α+3)]+72ℏ2x2[t2αΓ(2α+1)−2(Γ(α+2)Γ(α+1)+2)t2α+1Γ(2α+2)]+576ℏ2x2[(Γ(α+3)Γ(α+2)+1)t2α+2Γ(2α+3)+2Γ(α+4)t2α+3Γ(2α+4)],and so on. Hence, the HAM series
solution is(4.29)u(x,t)=u0+u1+u2+u3+⋯=x2−2x2t−6ℏx2[1+(ℏ+1)+⋯]tαΓ(α+1)+24ℏx2[1+(ℏ+1)+⋯]tα+1Γ(α+2)−48ℏx2[1+(ℏ+1)+⋯]tα+2Γ(α+3)+72hℏ2x2t2αΓ(2α+1)+⋯.The four-term HAM approximate
solutions for (4.14)-(4.15), when α=1.5, 1.75,
and 2 with ℏ=−1 and −1.0966,
are shown in Table 2. Notice that the HAM approximate solution when α=2 with ℏ=−1.0966 is in good agreement with the exact solution, u(x,t)=(x/(t+1))2.Table 2
Approximate
solution of (4.14) for some values of α using the 4-term HAM approximation, ϕ4,
with ℏ=−1 and ℏ=−1.0966.
α=1.5α=1.75α=2.0txℏ=−1ℏ=−1.0966ℏ=−1ℏ=−1.0966ℏ=−1ℏ=−1.0966Exact0.20.250.0600490.0601990.0487800.0487880.0434030.0434040.0434030.500.2401950.2407950.1951190.1951520.1736100.1736170.1736110.750.5404390.5417890.4390180.4390920.3906230.3906380.3906251.00.9607810.9631800.7804770.7806070.6944410.6944680.6944440.40.250.0745360.0769760.0454340.0458210.0318530.0319040.0318880.500.2981420.3079030.1817360.1832850.1274120.1276160.1275510.750.6708210.6927830.4089070.4123920.2866770.2871370.2869901.01.1925701.2316140.7269460.7331410.5096490.5104650.5102040.60.250.0914690.0997170.0459070.0480480.0239200.0244140.0244140.500.3658780.3988670.1836260.1921910.0956800.0976560.0976560.750.8232250.8974500.4131590.4324300.2152800.2197270.2197271.01.4635111.5954670.7345040.7687650.3827200.3906260.390625


## 4.4. Problem 4

Consider the following nonlinear time-fractional
Fisher’s equation:(4.30)Dαu(x,t)=uxx(x,t)+6u(x,t)(1−u(x,t)),for t>0,x∈R,0<α≤1, subject to the initial
condition(4.31)u(x,0)=1(1+ex)2.According to (3.2), we can choose
the initial guess to be(4.32)u0(x,t)=1(1+ex)2,and according to (3.18), we
have(4.33)um=(χm+ℏ)um−1−ℏJα[(um−1)xx+6um−1−6∑i=0m−1uium−1−i].Consequently, the first few terms of HAM series
solutions are as follows:(4.34)u1(x,t)=−10ℏex(1+ex)3tαΓ(α+1),u2(x,t)=−10ℏex(1+ex)3[(ℏ+1)tαΓ(α+1)−5ℏ(2ex−1)(1+ex)t2αΓ(2α+1)],and so on. Hence, the HAM series
solution is(4.35)u(x,t)=u0+u1+u2+⋯=1(1+ex)2−10ℏex(1+ex)3(1+(ℏ+1)+⋯)tαΓ(α+1)+50ℏ2ex(2ex−1)(1+ex)4t2αΓ(2α+1)+⋯.Table 3 shows the 3-term HAM
approximate solutions for (4.30)-(4.31), ϕ3,
when α=0.5,  0.75, and 1 with ℏ=−1 and −1.05133.
We notice that the HAM approximate solution when α=2 with ℏ=−1.05133 is in good agreement with the exact solution, u(x,t)=1/(1+exp(x−5t))2.Table 3
Approximate
solution of (4.30) for some values of α using the 3-term HAM approximation, ϕ3,
with ℏ=−1 and −1.05133.
α=0.5α=0.75α=1.0txℏ=−1ℏ=−1.0513ℏ=−1ℏ=−1.0513ℏ=−1ℏ=−1.0513Exact0.10.250.9461290.9840660.4881950.4969110.3179480.3196120.3160420.500.8439080.8835900.4057400.4149260.2505000.2522920.2500000.750.7150130.7525400.3244530.3331800.1909640.1926890.1916891.00.5764660.6091850.2496830.2573150.1409790.1425010.1425370.20.251.4753181.5517850.7912500.8165330.4811990.4884210.4612840.501.3598271.4396790.6901460.7166460.3969410.4045760.3874560.751.1809811.2564210.5744050.5994980.3152660.3225350.3160421.00.9700781.0358070.4566470.4785430.2411750.2475400.2500000.30.251.9674482.0825431.1242301.1711230.6814400.6981160.6041950.501.8452341.9653361.0094821.0585310.5818600.5993910.5344470.751.6229081.7363240.8595090.9058960.4758330.4924640.4612841.01.3455051.4442920.6954790.7359220.3729170.3874460.387456


## 5. Conclusions

In this work, the homotopy analysis method (HAM) was
implemented to derive exact and approximate analytical solutions for both
linear and nonlinear partial differential equations of fractional order. The
convergence region of the series solution obtained by HAM can be controlled and
adjusted by the auxiliary parameterℏ.
We give some examples to show the efficiency and accuracy of the suggested
method. It was also demonstrated that the Adomian decomposition method (ADM) is
a special case of HAM for the first and second test examples.



---

*Source: 686512-2009-01-29.xml*