SQU Journal for Science, 2014, 19 (1), 74-86 © 2014 Sultan Qaboos University 74 Homotopy Perturbation Method For Frac- tional Shallow Water Equations Kamel Al-Khaled 1* and Mohamed K. Al-Safeen2 1 Department of Mathematics and Statistics, College of Science, Sultan Qaboos University, P.O. Box: 36, PC 123, Al-Khod, Muscat, Sultanate of Oman; 2 Department of Mathematics and Statistics, Jordan University of Science and Technology, Irbid 22110, Jordan. *Email: kamel@squ.edu.om. ABSTRACT: In this paper, the homotopy perturbation method is adopted to find explicit and numerical solutions for systems of non-linear fractional shallow water equations. The fractional derivatives are described in the Caputo sense. We apply both the homotopy perturbation method and the homotopy analysis method, to solve certain shallow water equations with time-fractional derivatives, and explicitly construct convergent power series solutions. The results obtained reveal that these methods are both very effective and simple for finding approximate solutions. Some numerical examples and plots are presented to illustrate the efficiency and reliability of these methods. Keywords: Homotopy perturbation; Fractional differential equation; Shallow water equations; Approximate solutions. لحل معادالت المياه الضحلة الكسرية طريقة هوموتوبي للتشتت محمد السعافينكامل الخالد و تمثلكسرٌة التً التفاضلٌة الجزئٌة ال المعادالتمن ألنظمةٌجاد الحلول العددٌة للتشتت إل هوموتوبًطرٌقة فً هذا البحث تم استخدامملخص : التحلٌلٌة وطرٌقة هوموتوبً للتشتت هوموتوبًقمنا بتطبٌق كال الطرٌقتٌن, طرٌقة حلة. المشتقات الكسرٌة معرفة بطرٌقة كابوتا. المٌاه الض النتائج بٌنت ، متقاربةالت هندسٌة الحلول على شكل متسلس بإٌجادقمنا ، حلة ذات المشتقات الكسرٌة الزمنٌةالض لمٌاهانظام الذي ٌمثل لحل ال حٌث تم توضٌح ذلك من خالل بعض االمثلة والرسومات. ،ان الطرٌقتٌن تتسمان بالسهولة والدقة اضطراب هوموتوبً، المعادالت التفاضلٌة الكسرٌة، معادالت المٌاه الضحلة و الحلول التقرٌبٌة.مفتاح الكلمات: 1. Introduction he flow of water in a wide channel, with rectangular cross-section, and a smoothly varying bottom surface, is governed by the one-dimensional shallow water equations (SWEs). The assumptions of hydrostatic pressure distribution and small bottom slope are used in deriving these equations [1]. In the one- dimensional case, the flow of a fluid in an infinitely wide channel with time-fractional derivative can be presented by using shallow water equations )(' 0)( xGvvuvD uvuD xxt xt     (1.1) T mailto:kamel@squ.edu.om HOMOTOPY PERTURBATION METHOD 75 subject to the conditions that the initial velocity and the initial height of the water be ),()0,( 1 xgxu  ),()0,( 2 xgxv  x .R The fractional differential operator  t D describes the fractional time-derivatives of order of the system (1.1), where 10  . The fractional derivatives are considered in the Caputa sense [2]. In the case 1 , the fractional differential equation reduces to the classical shallow water equations. In equation (1.1) the functions ),( txu and ),( txv represent the total height above the bottom and the flux velocity respectively, and 0)(' xG is the depth of a point from a fixed reference level of the water. Bermudez and Vazquez [3] failed to find an exact analytical solution for the shallow water equations when the bottom is not flat for 1 in a literature search. In last few years many papers have been devoted to the numerical solution of these equations using finite difference and finite element methods [3]. In [1], the authors consider the same problem with 1 , and they implement Adomian's decomposition method. In recent years, there has been a growing interest in the field of fractional calculus. Oldham and Spanier [4], Miller and Ross [5], Momani [6] and Podlubny [7,8] provide the history and a comprehensive treatment of this subject. Several fields of application of fractional differentiation and fractional integration are already well established, and some others are just starting. Many applications of fractional calculus can be found in the literature [9,10]. For this reason, it is important to solve time fractional partial differential equations. It was found that fractional time derivatives arise generally as infinitesimal generators of the time evolution when taking a long time scaling limit. Hence, the importance of investigating fractional equations arises from the necessity to sharpen the concepts of equilibrium, stability states, and time evolution in the long time limit. There has been some attempt to solve linear problems with multiple fractional derivatives. In this paper, we consider 10  , with 0)(' xG ,and with the same analysis as in [11,3], we implement the homotopy perturbation method (HPM), and homotopy analysis method (HAM) for obtaining explicit and numerical solutions of those fractional shallow water equations with source terms. We will also illustrate that how the homotopy methods help to achieve accurate solutions. 2. Basic principles of Fractional Calculus This section is devoted to a description of the operational properties for the purpose of acquainting the reader with sufficient fractional calculus theory. Many definitions and studies of fractional calculus have been proposed in the last two centuries. These definitions include those of, Riemman Liouville, Weyl, Reize, Campos, Caputa, and Nishimoto fractional operator. The Riemann Liouville definition of fractional derivative operator  a J is defined as follows: Definition 2.1 Let   R . The operator  a J defined on the usual Lebesque space ],[ 1 baL by         a a dttftxxfJ ,)()( )( 1 )( 1 and )()( 0 xfxfJ a  , for ,bxa  is called the Riemann- Liouville-fractional integral operator of order  . Properties of the operator  a J can be found in [4, 7], here we mention the following: For 0,],,[ 1  baLf and 1 1. )(xfJ a  exists for almost every ],[ bax  2. )()( xfJxfJJ aaa    3. )()( xfJJxfJJ aaaa   4.         )( )1( )1( axxJ a KAMEL AL-KHALED and MOHAMED K. AL-SAFEEN 76 As mentioned in [6,8], the Riemann-Liouville derivative has certain disadvantages when trying to model real- world phenomena with fractional differential equations. Therefore, we shall now introduce a modified fractional differentiation operator  D proposed by Caputo in his work on the theory of viscoelasticity [2]. Definition 2.2The fractional derivative of )(xf in the Caputo sense is defined as ,)()( )( 1 )()( )( 0 1 dttftx m xfDJxfD m x mm           where 0,,1  xNmmm  Lemma 2.1 If ,1 mm   and ],,[ 1 baLf  then )()( xfxfJD aa   , and       1 0 )( ! )( )0()()( m k k k aa k ax fxfxfDJ  The Caputo fractional derivative is considered in the Caputo sense. For more details on the geometric and physical interpretations of fractional derivatives of both Riemann-Liouville and Caputo types see [2, 6]. Definition 2.3 For m to be the smallest integer that exceeds  , the Caputo fractional derivatives of order 0 are defined as 1 0 ( , ) 1 ( , ) ( , ) ( ) , ( ) m t m m m u x t u x D u x t t d t m                   if mm  1 , and for Nm , we have ( , ) ( , ) ( , ) m m u x t u x t D u x t t t          For mathematical properties of fractional derivatives and integrals one can consult the above mentioned references. 3. Homotopy Methods Liao [12] and He [13] employed the basic ideas of homotopy in topology to propose general analysis methods for nonlinear problems, namely the homotopy analysis method (HAM), and homotopy perturbation method (HPM). These methods have been successfully applied to solve many types of nonlinear problems. To illustrate the basic ideas of these methods, we consider the general form of the partial differential equation 0)(][][  rfuAuD (3.1) with boundary conditions of the form    r n u uB ,0),( (3.2) where D is a nonlinear operator, A is a general differential operator, B is the boundary operator, )(rf a known analytical function, and  is the boundary of the domain  . Generally speaking, the operator A can be divided into two parts, namely L and N , where L is the linear operator, and N is the nonlinear operator. Therefore, equation (3.1) can be rewritten as 0)()()(  rfuNuL (3.3) By means of generalizing the traditional homotopy method, Liao [14] constructs a homotopy Rpr  ]1,0[:),( which satisfies the zero-order deformation equation 0 ( , ) (1 ) [ ( ) ( )] ( ) [ ] 0H p p L L u p H r D        (3.4) HOMOTOPY PERTURBATION METHOD 77 where ]1,0[p is the embedding parameter, 0 is a nonzero auxiliary parameter, 0)( rH is the an auxiliary function, 0 u is an initial guess of u , and  is an unknown function. It is important, that one has great freedom to choose auxiliary things in (HAM) (i.e., the initial guess, the auxiliary linear operator, the auxiliary parameter  , and the auxiliary function )(rH ). Obviously, from equation (3.4), we have 0 ( , 0) ( ) ( ) 0, ( ,1) [ ] [ ] ( ) 0H L L H D A f r            (3.5) The changing process of p from zero to unity is just that of ),( pr varies (or deforms) from the initial guess )( 0 ru to the solution u . In topology, this is called deformation, and ),()( 0 uLL  )()( rfA  are termed ‘homotopic’. Expanding );( pr in a Taylor series with respect to the embedding parameter p , we have m m m pruupr )(),( 1 0     (3.6) where 0 | ),( ! 1 )(     pm m m p pr m ru  . (3.7) As proved by Liao [14], if the auxiliary linear operator, the initial guess, the auxiliary parameter  , and the auxiliary function are properly chosen, the series (3.6) converges at 1p , and one has     1 0 )()( m m xuuru (3.8) which must be one of the solutions of the original nonlinear equations. For the sake of simplicity, the governing equation can be deduced from the zero-order deformation equation (3.4). Define the vector },...,{ 0 nm uuu   , and Differentiating equation (3.4) m times with respect to the embedding parameter p , and then setting 0p and finally dividing them by !m , we have the so-called mth -order deformation equation 1 1 [ ( ) ( )] ( ) ( ( )) m m m m m L u r u r H r R u r     (3.9) where 1 1 01 1 [ ( , )] ( ( )) | (1 ) ( ) ( 1)! m m m p mm D r p R u r f r m p            and       1,1 1,0 m m m  . (3.10) It should be emphasized that )(ru m for 1m is governed by the linear equation (3.9) with the linear boundary conditions that come from the original problem, which can be easily solved by symbolic computation software such as Maple or Mathematica. By directly substituting the series (3.6) into the zeroth- order deformation equation (3.4), and then equating the coefficients of the like powers of p , one can obtain exactly the same high-order deformation equation as (3.9), as proved by Sajid et al.[15]. As 1 and 1)( rH , equation (3.4) becomes 0 ( , ) (1 )[ ( ) ( )] [ ( )]H p p L L u p D      (3.11) KAMEL AL-KHALED and MOHAMED K. AL-SAFEEN 78 which is used mostly in the homotopy perturbation method [13], whereas the solution obtained directly, without using Taylor series [16]. For equation (3.11) the perturbation procedure is applied by assuming the p as a perturbation parameter, i.e. the solution can be written as a power series in p : ...)()()(),( 2 210  pruprurupr This method does not require an additional small parameter, this being its main advantage in comparison with the usual perturbation procedure. The approximate solution of equation (3.1), therefore, can be readily obtained: ...)()()()1;( 210  rururuH  (3.12) The convergence of the series (3.12) has been proved in [17]. In [18], the homotopy analysis method was used to derive the Adomian decomposition method, i.e., if the auxiliary parameter 1 , the auxiliary function 1)( rH , the auxiliary linear operator, and the initial guess are properly chosen, the Adomian decomposition is a special case of the homotopy analysis method, as Allan [18] shows. 4. Solutions of Fractional Shallow Water Equations We solve the fractional system (1.1), subject to the initial conditions 1 2 sec ( ) ( , 0) ( ) ( ); ( , 0) 0 ( ), 4 h x u x G x g x v x g x x R      (4.1) We take the reference point to be 1, and )exp(1 )exp( )( 2 2 x x xG    . According to [19], we choose the linear operators as: ),;,()];,([ ptxDptxL itii    with property 0][ cL i , for ,2,1i where c is a constant. Also, we choose the system of nonlinear operators as 1 1 2 1 2 [ ( , ; ), ( , ; ) [ ( , ; ) ( , ; )]N x t p x t p x t p x t p x        (4.2) And, 2 1 2 2 2 1 ( ( , ; ), ( , ; )] '( ) ( , , ) ( , , ) ( , , )N x t p x t p G x x t p x t p x t p x x              4.1 The HPM solution We implement the Homotopy Perturbation Method (HPM) to solve the system (1.1). The Homotopy equations are: For 2,1i , we have ,0 2 (1 ) [ ( , ; ) ( , )] [ [ ( , ; ) [ ( , , ), ( , , ) ( )] 0 i i i i i i i i p L x t p z x t p L x t p N x t p x t p f x         (4.3) where [0,1]p  is an embedding parameter, and ,0 ( , ), 1, 2 i z x t i  , are the initial approximations satisfying the given conditions, while );,( ptx i  , for 1, 2i  are approximate solutions. In view of the homotopy perturbation method, expanding the solutions of equations (4.3) as a power series in p as: For 1, 2i  , we have ...),(),();,( 1,0,  ptxztxzptx iii  (4.4) where the initial approximations are given by HOMOTOPY PERTURBATION METHOD 79 1 1,0 2 2,0 ( , ; ) ( , ) ( , 0), ( , ; ) ( , ) ( , 0)x t p z x t u x x t p z x t v x     Substituting the above, and arranging coefficients of the same power of p , we get Powers of p : 1,1 2,0 1,0 2,1 1,0 2,0 2,0 ( , ) [ ( , ) ( , )] 0 ( , ) [ ( , )] '( ) ( , ) [ ( , )] 0 t t D z x t z x t z x t x D z x t z x t G x z x t z x t x x               Powers of 2 p : 1,2 2,1 1,0 2,0 1,1 ( , ) [ ( , ) ( , ) ( , ) ( , )] 0 a D z x t z x t z x t z x t z x t x       And 0)],([),()],([ ),()],([),( 0,21,21,2 0,21,12,2           txz x txztxz x txztxz x txzD t  (4.5) and so on for the rest of the coefficients of the same powers of p . The method of finding the solution is based on applying the operator  J on both sides of equations (4.5) with initial approximations 0)0,( , xz mi , for 1, 2m  and 1, 2i  to obtain 2 1,0 1,12 0 1 3 2 1,2 1,2 1,2 1,2 1,2 sec ( ) exp( ) ( , ) , ( , ) 0 4 1 exp( ) 1 ( , ) [ ( , ) ( , )] [[ ( , ) ( , )] 4 (1 2 ) h x x z x t z x t x z x t S x t S x t S x t S x t              (4.6) And 2,0 2,2 2,1 sec ( ) tanh( ) ( , ) 0 , ( , ) 0 , ( , ) 4 (1 ) t h x x z x t z x t z x t        (4.7) Where 2 3 0 1,2 1,0 2 2 1 1,2 1,0 2 2 2 1,2 2 2 2 3 2 1,2 sec ( ) ( , ) [ ( , )] 4 (1 2 ) sec ( ) tanh ( ) ( , ) [ ( , )] 4 (1 2 ) 2 exp( 2 ) 2 exp( ) sec ( ) tanh( ) ( , ) (1 exp( )) 1 exp( ) 4 ( , ) sec ( ) tanh( ) t h x S x t z x t t h x x S x t z x t x x x x h x x S x t x x S x t t h x x                        and so on. In the same manner, the rest of the components of the homotopy perturbation solutions can be obtained using Maple Package. The third-order term approximate solutions for equation (1.1) are given by KAMEL AL-KHALED and MOHAMED K. AL-SAFEEN 80 1,0 1,1 1,2 2,0 2,1 2,2 ( , ) ( , ) ( , ) ( , ), ( , ) ( , ) ( , ) ( , )u x t z x t z x t z x t v x t z x t z x t z x t      (4.8) 4.2 The HAM solution We implement the homotopy analysis method to solve the system (1.1). Thus we construct the zeroth- order deformation equations ,0 1 2 (1 ) [ ( , ; ) ( , )] ( , )[ [ ( , ; )] [ ( , ; ), ( , ; )] ( )] i i i i i i i i i p L x t p z x t p H x t L x t p N x t p x t p f x          (4.9) For .2,1i When ,0p we get 1 ,0 2 ,0 ( , ; 0) ( , ) ( , 0), ( , ; 0) ( , ) ( , 0) i i x t z x t u x x t z x t v x     While for 1p , we have ).,()1;,(),,()1;,( 21 txvtxtxutx   Expanding );,( ptx i  in a Taylor series with respect to p , then for 1, 2i  one has m m miii ptxztxzptx ),(),();,( 1 ,0,     where 0, | );,( ! 1 ),(     pm i m mi x ptx m txz  . Assuming that the above series is convergent when 1p  , due to (4.9) for 1, 2, 3i  we have     1 ,0, ),(),()1;,( m miii txztxztx . For the sake of simplicity, we can deduce that governing equations (4.9) from the zeroth-order deformation equations to the mth -order deformation equations are as follows: For 1, 2i  , , 1 , 1, 1 2, 1 [ ( , ) ( , )] ( , ) [ , ] i i m m i m i i m m m L z x t z x t H x t R z z      (4.10) where         1 0 1,2,1 1,11,21,1,1 ),(),([ ),(],[ m n nmn mtmmm txztxz x txzDzzR  1 2, 1, 1 2, 1 2, 1 1, 1 2, 2, 1 0 [ , ] ( , ) ( , ) ( , ) ( , ) (1 ) '( ) m m m m t m m n m n m n R z z D z x t z x t z x t z x t G x x x                    Now, applying the operator  t J on both sides of , , 1 , 1, 1 2, 1 [ ( , ) ( , )] ( , ) [ , ] t i m m t i m i i i m m m D z x t D z x t H x t R z z         yields that for 1, 2i  , , , 1 , 1 , 1, 1 2, 1 ( , ) ( , 0) [ ( , ) ( , 0)] [ ( , ) [ , ]] i m i m m i m i m i t i i m m m z x t z x z x t z x J H x t R z z           (4.11) For simplicity, we assume , 21   and 1),(),( 21  txHtxH , and ,0)0,( , xz mi for 1, 2i  and 1, 2,...m  . According to equation (4.11), we now successively obtain HOMOTOPY PERTURBATION METHOD 81 2 1,0 1,12 0 1 3 2 1,2 1,2 1,2 1,2 1,2 exp( ) sec ( , ) , ( , ) 0, 1 exp( ) 4 ( , ) [ ( , ) ( , )] [ ( , ) ( , )] 4 (1 2 ) x hx z x t z x t x z x t S x t S x t S x t S x t             (4.12) and )1(4 tanhsec)1( ),( )1(4 tanhsec ),( 0),( 2,2 1,2 0,2            xhxth txz xhxt txz txz   (4.13) where 2 3 2 0 1,2 2 sec exp( ) sec ( , ) [ ] 4 (1 2 ) 1 exp( ) 4 t h x x hx S x t x           2 2 2 1 1,2 2 sec tanh exp( ) sec ( , ) [ ] 4 (1 2 ) 1 exp( ) 4 t hx x x hx S x t x          2 2 2 1,2 2 2 2 2 exp( 2 ) 2 exp( ) sec tanh ( , ) (1 exp( ) 1 exp( ) 4 x x x x hx x S x t x x          xhxttxS tanhsec),( 23 2,1  and so on. In the same manner, the rest of the components of the homotopy analysis solutions can be obtained using Maple Package. The third-order term approximate solutions for equation (1.1) are given by 1,0 1,1 1,2 2,0 2,1 2,2 ( , ) ( , ) ( , ) ( , ); ( , ) ( , ) ( , ) ( , )u x t z x t z x t z x t v x t z x t z x t z x t      5. Numerical Results In this section we shall illustrate the two techniques by different values of  for the fractional shallow water equations. All the results are calculated by using the symbolic calculus software Mathematica, using the (HPM). The result shown in Figure 1 represents the height of the water, when 1 , which agrees with the results in [1] in which the Adomian decomposition method was used. Figure 2 represents the velocity of the water when 1 , which begins with zero and, at any later instant of time, changes from positive to negative as x changes from positive to negative. Figure 3 shows the height of water when 25.0 , while Figure 4 shows the velocity of water for 25.0 . Figures 75  show two profiles of the height of the wave for different values of  when the time changes from 1 to 2 . Figures 108  show again profiles for the velocity. From the Figures 5-10, it is observed that the physical behavior of the approximate solutions at different values of  has the same behavior as that obtained when  =1. This shows the approximate solution is efficient. Comparison of Figures 5-10 shows that the solution continuously depends on the fractional derivatives. A clear conclusion can be drawn from Figures 5-10 that the solution tends to zero as |x| approaches infinity. It is also observed that as time changes from 1 to 2, the wave starts to bifurcate into two waves. We conclude that as the time-fractional derivative increases, the wave solution also changes continuously. Tables 1 and 2 show some numerical values for different values of  using (HPM) and (HAM), where only the third-order term of the approximate solutions were used. From the Tables, we can read that results presented by (HAM) when 1 can be considered as results of (HPM). KAMEL AL-KHALED and MOHAMED K. AL-SAFEEN 82 6. Conclusions The fundamental goal of this paper is to propose an efficient algorithm for the solution of fractional shallow water equations. This goal has been achieved by generalized mth -order deformation of HAM for approximate solutions of non-homogeneous systems. It should be emphasized that, within the frame of the homotopy perturbation method, we have great freedom to choose the initial guess, and the auxiliary linear operator L . The convergence of the series solution by HAM is dependent upon four factors [15,20], i.e., the initial guess, the auxiliary linear operator, the auxiliary function ),( txH , and the auxiliary parameter  , while the convergence of series solution that is obtained by HPM is only depend upon two factors: the auxiliary linear operator, and the initial guess. Using the HPM we overcome the difficulties arising in the calculation of HAM; moreover the calculations in the HPM are simple and straightforward. The reliability of the method and the reduction in the size of the computational domain give this method a wide applicability. Table 1. Numerical values of ),( txu using HPM. t x  =0.5  =0.75  =1 0.2 0.5 0.8 0.2 0.4 0.7 0.2 0.4 0.7 0.2 0.4 0.7 0.703295 0.673630 0.583692 0.655611 0.647070 0.590625 0.607928 0.620510 0.597558 0.724389 0.685380 0.580625 0.692810 0.667791 0.585216 0.649529 0.643682 0.591509 0.731904 0.689566 0.579532 0.715215 0.680270 0.581959 0.684221 0.663006 0.586465 Table 2. Numerical values of ),( txv using HPM. t x  =0.5  =0.75  =1 0.2 0.5 0.8 0.2 0.4 0.7 0.2 0.4 0.7 0.2 0.4 0.7 0.024410 0.044338 0.060744 0.038596 0.070105 0.096046 0.048820 0.088676 0.121490 0.015741 0.028591 0.039171 0.031295 0.056845 0.077879 0.044522 0.080860 0.110793 0.009674 0.017572 0.024075 0.024186 0.043932 0.060187 0.038698 0.070291 0.096300 HOMOTOPY PERTURBATION METHOD 83 Figure 1. Approximate numerical solution ),( txu with 1 Figure 2. Approximate numerical solution ),( txv with 1  Figure 3. Approximate numerical solution ),( txu with 0.25  KAMEL AL-KHALED and MOHAMED K. AL-SAFEEN 84 Figure 4. Approximate numerical solution ),( txv with 0.25  Figure 5. Approximate numerical solution ),( txu with 0.25  for 1t  to 2t  Figure 6. Approximate numerical solution ),( txu with 1  for 1t  to 2t  HOMOTOPY PERTURBATION METHOD 85 Figure 7. Approximate numerical solution ),( txu with 0.5  for 1t  to 2t  Figure 8. Approximate numerical solution ),( txv with 1  for 1t  to 2t  Figure 9. Approximate numerical solution ),( txv with 0.5  for 1t  to 2t  Figure 10. Approximate numerical solution ),( txv with 0.25  for 1t  to 2t  KAMEL AL-KHALED and MOHAMED K. AL-SAFEEN 86 References 1. Al-Khaled, K. and Allan, F. Construction of solutions for the shallow water equations by the decomposition method. Comput. Simul., 2004, 66, 479-486. 2. Caputo, M. Linear models of dissipation whose Q is almost frequency independent II, Geophysics. J. R. Astron. Soc., 1967, 13, 529-539. 3. Bermudez, A. and Vazquez, M. Upwind methods for hyperbolic conservation laws with source terms. Comput. Fluids, 1994, 23,1049-1071. 4. Oldham, K.B. and Spanier, J. The Fractional Calculus, Academic Press, New York, 1974. 5. Miller, K.S. and Ross, B. An Introduction to the Fractional Calculus and Fractional Differential Equations, John Wiley and Sons Inc. New York, 1993. 6. Momani, S. Non-perturbative analytical solutions of the space and time fractional Burgers' equations. Chaos Soliton Fractals., 2006, 28, 930-937. 7. Podlubny, I. Fractional Differential Equations, Academic Press, New York, 1999. 8. Podlubny, I. Gemmetric and physical interpretation of fractional integration and fractional differentiation. Frac. Calc. Appl. Anal. 2002, 5, 367-386. 9. Kisela, T. Power functions and essentials of fractional calculus on isolated time scales. Adv. Diff. Equ., 2013, 2013, 259. 10. Podlubny, I. What Euler could further write, or the unnoticed "big bang" of the fractional calculus. Frac. Calc. Appl. Anal. 2013, 16, 501-506. 11. Liao, S.J. Comparison between the homotopy analysis method and homotopy perturbation method. Appl. Math. Comput., 2005, 169, 1186-1194. 12. Liao, S.J. The proposed homotopy analysis technique for the solution of nonlinear problems, Ph.D. thesis, Shanghai Jiao Tong University, 1992. 13. He, J.H. Homotopy-Perturbation Method: a new nonlinear analytical technique. Appl. Math. Comput., 2003, 135, 73-79. 14. Liao, S.J. Beyond Perturbation: Introduction to the Homotopy Analysis Method, Chapman and Hall/CRC Press, Boca Raton, 2003. 15. Sahid, M., Hayat, T. and Asghar, S. Comparison between the HAM and HPM solutions of thin film flows of non-Newtonian fluids on a moving belt. Nonlin. Dynam., 2007, 50, 27-35. 16. Abbasbandy, S. The application of homotopyanalysis method to solve a genaralized Hiroto-Satsuma coupled KdV equations. Phys. Lett. A , 2007, 361, 478-483. 17. He, J.H. Homotopy-Perturbation Technique. Comput. Math. Appl. Mech. Eng., 1999, 178, 257-262. 18. Allan, F. Derivation of the Adomian decomposition method using the homotopy analysis method. Appl. Math. Comput., 2007, 190, 6-14. 19. Al-Saafeen, M.K. Homotopy Perturbation Method for Fractional Hirota-Satsuma Coupled KdV Equations, Master Thesis, Jordan University of Science and Technology, 2008. 20. Song, L.N. Application of homotopy analysis method to fractional KdV-Burgers-Kuramoto equation. Phys. Lett. A, 2007, 367(1-2), 88-94. Received 10 February 2014 Accepted 14 April 2014