Progr. Fract. Differ. Appl. 2, No.

3, 227-232 (2016) 227

Progress in Fractional Differentiation and Applications

An International Journal

[Link]

A Highly Accurate Numerical Method for Solving

Time-Fractional Partial Differential Equation
Muhammad Khalid1,∗ , Fareeha Sami Khan1 , Husna Zehra2 and Muhammad Shoaib3
1 Department of Mathematical Sciences, Federal Urdu University of Arts, Sciences and Technology, Karachi-75300, Pakistan.
2 Department of Mathematics, NED University of Engineering and Technology, University Road, Karachi-75270, Pakistan.
3 Department of Physics, Federal Urdu University of Arts, Sciences and Technology, Karachi-75300, Pakistan.

Received: 2 Mar. 2016, Revised: 18 May 2016, Accepted: 20 May 2016

Published online: 1 Jul. 2016

Abstract: An efficient numerical method is employed to approximate the numerical solutions of some very functional, time-fractional
partial differential equations. Perturbation Iteration Algorithm applied on fractional PDEs can technically manipulate non-linear and
fractional terms pretty well. Similarly, the precision of its results are even better then that of different techniques. Explanatory figures
have been presented correlating the approximated and exact solutions and substantiating the precision of results.

Keywords: Perturbation iteration algorithm, advection diffusion equation, Fisher equation, hyperbolic partial differential equation.

1 Introduction

Fractional calculus is yet another subject of this era, attracting a large number of mathematical analysts who have been
dealing with this subject [1–3]. Different methods have been presented dealing with partial and ordinary differential
equations of fractional order. Some of them are proposed numerical solutions, e.g. Variational Iteration Method (VIM)
[4, 5], Homotopy Analysis method (HAM) [6, 7], and Adomian Decomposition Method (ADM) [8, 9], whereas, some
mathematical analysts even proposed analytical methods for time-fractional partial differential equation, such as Iteration
Method [10], Fourier transform method [11], Sumudu Transform Method [12], Greens Function Method [13] and Laplace
Transform Method [10–14].
Nonlinear fractional partial differential equations (FPDEs) are offshoots of established ordinary differential equations. If
first order derivative be replaced by a fractional order (single or multiple fractions) derivative, in basic PDE, a fractional
order PDE is acquired subsequently. Both, linear and non-linear FPDEs vitally contribute in the fields of social sciences,
engineering and many physical phenomena such as [15, 16].
In this paper, the time-fractional advection, hyperbolic and Fisher partial differential equations have been numerically
solved, that have been worked out by a good number of [Link] time-fractional diffusion equation as considered
by Wyss [17] through their solution in closed form in terms of H- function. In the work of Schneider and Wyss [18],
consideration of wave and fractional diffusion is found. Srivastava et al, numerically solved time fractional hyperbolic
telegraph equation by RDTM [19]. Fisher equation was initially designed by [20] for the breeding of a virile gene and
can be confronted by many chemical reactions such as the Brownian motion [21], chemical kinetics [22], auto catalytic
chemical reaction [23] etc. Neamaty [24] has analyzed time fractional PDEs by applying VIM, and also given comparison
of their results with different results of other numerical techniques. The analytical area has been segmented as:

[Link] definitions of fractional calculus.

[Link] Theory of Perturbation Iteration Algorithm on FPDEs.
[Link] Examples
[Link]
∗ Corresponding author e-mail: khalidsiddiqui@[Link]
c 2016 NSP
Natural Sciences Publishing Cor.
228 M. Khalid et al.: A Highly accurate numerical method...

2 Employed Definitions of Fractional Calculus

Most normally utilized definitions of fractional derivatives and integrals are Riemann-Liouville and Caputo sense,
therefore here is a brief introduction to these concepts.
Definition 2.1 The Riemann-Liouville fractional integral operator of order α > 0, of a function ∈ Cµ , µ > −1, is defined
as
Z t
1
J α y(t) = (t − s)α −1 y(s)ds; α > 0, (1)
Γ (α ) 0

J ◦ y(t) = y(t)
Some properties of the operator J α , used in this text are:

J α J β y(t) = J α +β y(t); α , β ≥ 0,

Γ (m + 1) α +m
Jα t m = t ; m ≥ −1.
Γ (m + α + 1)
Definition 2.2 A real valued function y(x), x > 0 is said to be in space Cµ , µ ∈ R, if there exists a real number p > µ , such
that y(t) = t p y1 (t), where y1 (t) ∈ C(0, inf), and it is said to be in the space Cµn if and only if yn ∈ Cµ , n ∈ N.
Definition 2.3 The fractional derivative of y(t) in the Caputo sense is defined as

Dα y(t) = J m−α Dm y(t); f or m − 1 < α ≤ m, m ∈ N, t > 0, and y ∈ C−1

m
. (2)

Firstly, Caputo fractional derivative evaluates only an ordinary derivative then through fractional integral obtains the
required fractional derivative. This approach Riemann- Liouville fractional integral operator resembles very much the
integer order integration so is a linear operation.

J α Σi=1 ci J α yi (t);
 
n
ci yi (t) = Σi=1
n
where {ci }ni=1 are constants. (3)

3 Perturbation Iteration Algorithm (PIA)

Step I
Consider an initial value problem such as

Dtα y + M(yxx , ytt , yxt , yx , yt , y) + H(yxx , ytt , yxt , yx , yt , y) = g(x,t); 0 < α ≤ 1, t > 0, x ∈ R (4)
∂k
with initial condition y(x, 0) = yk (x); k = 0, 1, 2, ..., m − 1. Where y = y(x,t), H is the linear operator, M is the
∂ tk
nonlinear operator and g(x,t) is the known analytic function.
Step II
Introducing ε with nonlinear term yield

Dα y + ε M + H − g(x,t) = 0 (5)
Here PIA(1, 1) will be considered, which means only one correction term in this expansion will be obtained by taking
n = 1, m = 1.
Step III
Consider the following fractional order differential equation.

F(Dtα y, yxx , ytt , yxt , yx , yt , y, ε ) = 0. (6)

By applying PIA on Eq.7 as [25, 26]. Only n correction terms in this perturbation expansion will be acknowledged

yn+1 = yn + ε (yn )c , (7)

where ε is the perturbation parameter. The developed Perturbation Iteration Algorithm is given here as PIA(n, m); here n
are the terms involved in the expansion , m is the mth order derivative in the Taylor’s Series expansion provided n ≤ m,

c 2016 NSP
Natural Sciences Publishing Cor.
Progr. Fract. Differ. Appl. 2, No. 3, 227-232 (2016) / [Link]/[Link] 229

w.r.t ε then comparing the coefficient of same power of ε renders the unknown correction terms. Back substitution of
these results in Eq. (7)thus yields an algorithm for the solution of Eq. (4). Substituting Eq. (8) in Eq. (7), expanding in a
Taylor’s Series with first derivative only yields

F(Dtα y, yxx , ytt , yxt , yx , yt , y, 0) + FDα y (Dtα y, yxx , ytt , yxt , yx , yt , y, 0)ε (Dα y)n +
 
c
α α
   
Fyxx (Dt y, yxx , ytt , yxt , yx , yt , y, 0)ε (yxx )n + Fytt (Dt y, yxx , ytt , yxt , yx , yt , y, 0)ε (ytt )n +
c  c
α α
 
Fyxt (Dt y, yxx , ytt , yxt , yx , yt , y, 0)ε (yxt )n + Fyt (Dt y, yxx , ytt , yxt , yx , yt , y, 0)ε (yt )n + (8)
 c  c
Fyx (Dtα y, yxx , ytt , yxt , yx , yt , y, 0)ε (yx )n + Fy (Dtα y, yxx , ytt , yxt , yx , yt , y, 0)ε (y)n +
 
c c
Fε (Dtα y, yxx , ytt , yxt , yx , yt , y, 0)ε = 0.

All above derivatives will be taken at ε = 0. First (y◦ )c has been calculated by using initial condition y◦ (x,t) and Eq.(8).
Then we substitute (y◦ )c into Eq. (7) to find y1 . Iteration process is repeated using Eq.(7) and Eq.(8) until we obtain a
satisfactory result.

4 Numerical Examples

4.1 Example:

Table 1: Result Comparison of PIA and other numerical methods provided in [24] for Eq. (9)
x t VIM ADM HPM VHPIM PIA Exact
0.25 0.0503090 0.0500000 0.0499876 0.0499876 0.0500001 0.0500000
0.50 0.1006190 0.1000000 0.0999780 0.0999746 0.1000002 0.1000000
0.2
0.75 0.1509280 0.1500010 0.1499680 0.1499620 0.1500004 0.1500000
1.00 0.2012370 0.2000010 0.1999570 0.1999510 0.2000005 0.2000000
0.25 0.1018940 0.1000230 0.0995290 0.0996450 0.1000158 0.1000000
0.50 0.2037870 0.2000460 0.1990590 0.1992900 0.2000316 0.2000000
0.4
0.75 0.3056810 0.3000690 0.2985880 0.2989350 0.3000475 0.3000000
1.00 0.4075750 0.4000920 0.3981180 0.3985800 0.4000633 0.4000000
0.25 0.1530940 0.1504110 0.1471580 0.1456900 0.1502739 0.1500000
0.50 0.3061880 0.3008230 0.2943170 0.2913800 0.3005478 0.3000000
0.6
0.75 0.4592820 0.4512340 0.4414750 0.4370700 0.4508218 0.4500000
1.00 0.6123760 0.6016460 0.5886340 0.5827590 0.6010957 0.6000000

Consider the time-fractional advection partial differential equation

Dtα y(x,t) + y(x,t)yx (x,t) = x(1 + t 2); t > 0, x ∈ R, 0 < α ≤ 1 (9)

with initial condition; y(x, 0) = 0. Applying perturbation parameter ε on non-linear and fractional terms, time-fractional
advection equation becomes
1
Dtα y(x,t) + ε y(x,t)yx (x,t) = ε x( + t 2 ).
ε
By applying PIA, the obtained corrected term is

yc (x,t) = J α − Dtα y(x,t) − y(x,t)yx (x,t) + x + xt 2

 

following iterations have been obtained

y◦ (x,t) = 0,

c 2016 NSP
Natural Sciences Publishing Cor.
230 M. Khalid et al.: A Highly accurate numerical method...

1 2t 2 
y1 (x,t) = xt α

+ ,
Γ (1 + α ) Γ (3 + α )
1 2t 2 t 2α Γ (1 + 2α ) 4t 2+2α Γ (3 + 2α )
y2 (x,t) = xt α
 
+ − − .
Γ (1 + α ) Γ (3 + α ) Γ (1 + α )2Γ (1 + 3α ) Γ (1 + α )Γ (3 + α )Γ (3 + 3α )

4.2 Example:

Table 2: Result Comparison of PIA and other numerical methods provided in [24] for Eq. (10)
t x VIM ADM HPM VHPIM PIA Exact
0.25 0.0434000 0.0433951 0.0434000 0.0432049 0.0434000 0.0434030
0.50 0.1736000 0.1735800 0.1736000 0.1728200 0.1735999 0.1736110
0.2
0.75 0.3906000 0.3905560 0.3906000 0.3888440 0.3905998 0.3906250
1.00 0.6944000 0.6943210 0.6944000 0.6912780 0.6943997 0.6944440
0.25 0.0317790 0.0315670 0.0317790 0.0299125 0.0317795 0.0318880
0.50 0.1271180 0.1262680 0.1271180 0.1196500 0.1271179 0.1275510
0.4
0.75 0.2860150 0.2841030 0.2860150 0.2692120 0.2860152 0.2869900
1.00 0.5084710 0.5050720 0.5084710 0.4786000 0.5084715 0.5084710
0.25 0.0236650 0.0220050 0.0236650 0.0188604 0.0236649 0.0244140
0.50 0.0946600 0.0880180 0.0946600 0.0754415 0.0946595 0.0976560
0.6
0.75 0.2129840 0.1980400 0.2129840 0.1697430 0.2129839 0.2197270
1.00 0.3786380 0.3520710 0.3786380 0.3017660 0.3786380 0.3906250

Consider the time fractional hyperbolic equations.

Dtα y(x,t) = y(x,t)yx (x,t) ;

 
t > 0, x ∈ R, 1 < α ≤ 2 (10)
x

with initial condition; y(x, 0) = x2 , yt (x, 0) = −2x2 . Applying perturbation parameter ε on non-linear and fractional
terms, time-fractional hyperbolic equation becomes

Dtα y(x,t) − ε y(x,t)yx (x,t) = 0.

 
x
By applying PIA, the obtained corrected term is

yc (x,t) = J α − Dtα y(x,t) + y(x,t)yx (x,t)

  
− 2x2t + x2 .
x
Iterations obtained by adding initial condition in corrected term are as follows

y◦ (x,t) = x2 (1 − 2t),
 6t 24t 1+α 48t 2+α 
y1 (x,t) = x2 1 − 2t + − + ,
Γ (1 + α ) Γ (2 + α ) Γ (3 + α )
 6t 24t 1+α 48t 2+α 72t 2α 288t 1+2α 
y2 (x,t) = x2 1 − 2t + − + + − .
Γ (1 + α ) Γ (2 + α ) Γ (3 + α ) Γ (1 + 2α ) Γ (2 + 2α )

4.3 Example:
Consider the time-fractional Fisher’s equation

Dtα y(x,t) = yxx (x,t) + 6y(x,t) 1 − y(x,t) ;

 
t > 0, x ∈ R, 0 < α ≤ 1, (11)

c 2016 NSP
Natural Sciences Publishing Cor.
Progr. Fract. Differ. Appl. 2, No. 3, 227-232 (2016) / [Link]/[Link] 231

Table 3: Result Comparison of PIA and other numerical methods provided in [24] for Eq. (11)
x t VIM ADM HPM VHPIM PIA Exact
0.25 0.3159400 0.3179480 0.3159400 0.3280190 0.3159398 0.3160420
0.50 0.2499260 0.2505000 0.2499260 0.2565130 0.2499257 0.2500000
0.1
0.75 0.1916060 0.1909640 0.1916060 0.1943030 0.1916059 0.1916890
1.00 0.1424110 0.1409790 0.1424110 0.1427150 0.1424105 0.1425370
0.25 0.4593200 0.4811990 0.4593200 0.5121930 0.4593203 0.4612840
0.50 0.3864500 0.3969410 0.3864500 0.4146970 0.3864505 0.3874560
0.2
0.75 0.3154780 0.3152660 0.3154780 0.3247160 0.3154775 0.3160420
1.00 0.2490920 0.2411750 0.2490920 0.2458810 0.2490923 0.2500000
0.25 0.5911790 0.6814400 0.5911790 0.6302750 0.5911793 0.6041950
0.50 0.5276350 0.5818610 0.5276350 0.5076430 0.5276353 0.5344470
0.3
0.75 0.4597190 0.4758330 0.4597190 0.4882980 0.4597193 0.4612840
1.00 0.3870250 0.3729170 0.3870250 0.3784720 0.3870253 0.3874560

with initial condition; y(x, 0) = (1+e1 x )2 . Applying perturbation parameter ε on non-linear and fractional terms, then
time-fractional Fisher’s equation becomes

Dtα y(x,t) − ε yxx (x,t) − 6ε y(x,t) 1 − y(x,t) = 0.

 

By applying PIA, the obtained corrected term is

yc (x,t) = J α − Dtα y(x,t) + yxx (x,t) + 6y(x,t) − 6y(x,t)2 .

 

Iterations obtained by adding initial condition in corrected term are as follows

1
y◦ (x,t) = ,
(1 + ex)2

1 tα  6 6e2x 2ex 6 
y1 (x,t) = − + − + ,
(1 + ex )2 Γ (1 + α ) (1 + ex )4 (1 + ex )4 (1 + ex )3 (1 + ex )2

1 tα  6 6e2x 2ex 6 
y2 (x,t) = − + − + −
x
(1 + e )2 Γ (1 + α ) (1 + e )
x 4 x
(1 + e )4 x
(1 + e )3 x
(1 + e )2

t 2α   600t α Γ (1 + 2α )
50ex + 150e3x + 100e4x − .
(1 + e ) Γ (1 + 2α )
x 6 Γ (1 + α )2Γ (1 + 3α )

5 Conclusions

In this work a powerful and easily manageable numerical method PIA has been applied on three different time space
fractional partial differential equations. This method uses Riemann-Liouville and Caputo definitions for fractional
integration and differentiation. Results obtained by PIA in this work has been compared by the approximated results
given in [24] by different methods such as VIM, HPM , ADM and VHPIM. Also it can easily be observed that results of
this numerical method is more accurate and convergent than other numerical techniques especially comparison by VIM
shows the efficiency of PIA. It is recommended that this satisfactory method ought to be utilized vivaciously for other
complex dynamical systems.

Acknowledgement

The authors thank the reviewers for their thorough efforts in editing this publication and highly appreciate the comments
and constructive criticism that significantly contributed in improving its quality . The authors also appreciate Ms. Wishaal
Khalid for proofreading the research paper.

c 2016 NSP
Natural Sciences Publishing Cor.
232 M. Khalid et al.: A Highly accurate numerical method...

References
[1] K. B. Oldham and J. Spanier, The Fractional Calculus, Academic Press, New York and London, (1974).
[2] B. Ross, The Fractional Calculus and its Application, Springer-Verlag, Berlin, (1975).
[3] H. C. Torrey, Bloch equations with diffusion terms, Phys. Rev., 104, 563-565 (1956).
[4] K. Al-Khaled and S. Momani, An approximate solution for a fractional diffusion-wave equation using the decomposition method,
Appl. Math. Comput. 165, 473-483 (2005).
[5] N. H. Sweilam, M. M. Khader and R. F. Al-Bar, Numerical studies for a multi-order fractional differential equation, Phys. Lett. A
371, 26-33 (2007).
[6] L. Song and H. Zhang, Phys. Lett. A 367, 88-94 (2007).
[7] H. Jafari and V. Daftardar-Gejji, Appl. Math. Comput., 180, 488-497 (2006).
[8] C. Yang and J. Hou, J. Infor. Comput. Sci. 10, 213-222 (2013).
[9] M. Al-Refai and M. A. Hajji, Nonlinear Anal. A 74, 3531-3539 (2011).
[10] I. Podlubny, An ntroduction to Fractional Derivatives, Fractional Equation, to Methods of Their Solution and Some of Their
Application, Academic Press, New York, (1999).
[11] R. L. Magin and M. Ovadia, J. Vibr. Contr. 14, 1431-1442 (2008).
[12] V. G. Gupta and B. Sharma, Appl. Math. Sci. 4, 435-446 (2010).
[13] F. Mainardi, IUTAM symposium-nonlinear waves in solid Fairfield, (1995).
[14] Z. Odibat and S. Momani, Phys. Lett. A 365, 351-357 (2007).
[15] M. Javidi and B. Ahmad, [Link]. Equ. 375, 375-392 (2013).
[16] A. E. Mohamed, A. A. Ahmed, B. B. Dumitru, Rom. J. Phys. 55, 274-284 (2010).
[17] W. Wyss, J. Math. Phys. 27, 2782-2785 (1986).
[18] W. R. Schneider and W. Wyss, J. Math. Phys. 30, 134-144 (1989).
[19] V. K. Srivastava, M. K. Awasthi and M. Tamsir, AIP Adv. 3, Article ID 032142 (2013).
[20] R. A. Fisher, Ann. Eugene 7, 335-369 (1937).
[21] M. D. Bramson, Commun. Pure Appl. Math. 31, 531-581 (1978).
[22] W. Malflict, American J. Phys. 60, 650-654 (1992).
[23] D. J. Aronson and H. F. Weinberg, Nonlinear diffusion in population genetics combustion and never pulse propapgation, Springer,
New York, 1988.
[24] A. Neamaty, B. Agheli and R. Darzi, Progr. Fract Differ. Equ. 1, 47-55 (2015).
[25] Y. Aksoy and M. Pakdermirli, Comput. Math. Appl. 59, 2802-2808 (2010).
[26] M. Pakdermirli, Y. Aksoy and H. Boyaci, Math. Comput. Appl. 16, 890-899 (2011).

c 2016 NSP
Natural Sciences Publishing Cor.