Acta Appl Math (2009) 105: 189–198 DOI 10.1007/s10440-008-9271-x

Series Solutions of Systems of Nonlinear Fractional Differential Equations A.S. Bataineh · A.K. Alomari · M.S.M. Noorani · I. Hashim · R. Nazar

Received: 22 April 2008 / Accepted: 27 June 2008 / Published online: 17 July 2008 © Springer Science+Business Media B.V. 2008

Abstract Differential equations of fractional order appear in many applications in physics, chemistry and engineering. An effective and easy-to-use method for solving such equations is needed. In this paper, series solutions of the FDEs are presented using the homotopy analysis method (HAM). The HAM provides a convenient way of controlling the convergence region and rate of the series solution. It is confirmed that the HAM series solutions contain the Adomian decomposition method (ADM) solution as special cases. Keywords Homotopy analysis method · System of nonlinear fractional differential equations · Caputo fractional derivative Mathematics Subject Classification (2000) 35Q53 · 41A58

1 Introduction The purpose of this paper is to find series solutions to systems of nonlinear fractional differential equations (FDEs) of the form: D αi yi = fi (t, y1 , . . . , yn ),

yik (0) = cki ,

0 ≤ k ≤ αi , i = 1, 2, . . . , n,

(1)

where D αi denotes Caputo derivative, fi are (linear or non-linear) real-valued functions and cki ∈ . In the last few decades, it has been shown that many phenomena cannot be described within the framework of the classical theory using integer order derivatives [1–3]. Fractional derivatives are very suitable, for example, for the description of mechanical and electrical properties of real materials. Other interesting applications involving fractional derivatives arise in dynamical systems theory [4, 5], modelling viscoelastic materials, modeling diffusion of atmospheric pollutants and heat and mass transfer problems. A.S. Bataineh · A.K. Alomari · M.S.M. Noorani · I. Hashim () · R. Nazar School of Mathematical Sciences, Universiti Kebangsaan Malaysia, 43600 UKM Bangi Selangor, Malaysia e-mail: [email protected]

190

A.S. Bataineh et al.

Finding accurate and efficient methods for solving FDEs has been an active research undertaking. Exact solutions of most of the FDEs cannot be found easily, thus analytical and numerical methods must be used. Some of the analytic methods for solving nonlinear problems include the Adomian decomposition method (ADM) [6–9], homotopy-perturbation method (HPM) [10–13] and variational iteration method (VIM) [14, 15]. Another analytical technique yielding series solutions is called the homotopy analysis method (HAM), initially proposed by Liao in his PhD thesis [16]. (For a systematic and clear exposition on the HAM, the reader is referred to [17, 18].) The HAM contains an auxiliary parameter  which provides us with a simple way to adjust and control the convergence region and rate of convergence of the series solution. We note that the HPM proposed in 1998 by He [19, 20] is exactly a copy of the early HAM (proposed in 1992) and is a special case of the recent HAM when  = −1. This fact has been confirmed by Abbasbandy [21], Sajid and Hayat [22] and Hayat and Sajid [23]. In recent years, the HAM has been successfully employed to solve many types of problems in science and engineering [21, 24–36]. Very recently, Song and Zhang [37] applied HAM to solve fractional KdVBurgers-Kuramoto equation. Hashim et al. [38] applied the HAM for fractional initial-value problems (IVPs). Cang et al. [39] presented the series solutions of non-linear Riccati differential equations with fractional order. Xu and Cang [40] used the HAM for time fractional wave-like equation. This paper investigates for the first time the applicability and effectiveness of the HAM on system of nonlinear fractional differential equations. Comparison with the results obtained by the ADM [6] shall be given. It is demonstrated that the solutions obtained by the ADM are special cases of the HAM solutions.

2 Preliminaries and Notations In this section, we give some definitions and properties of the fractional calculus [3]. Definition 1 A real function f (t), t > 0, is said to be in the space Cμ , μ ∈ , if there exists a real number p > μ, such that f (t) = t p f1 (t), where f1 (t) ∈ C(0, ∞), and it is said to be in the space Cμn if and only if h(n) ∈ Cμ , n ∈ N . Definition 2 The Riemann-Liouville fractional integral operator (J α ) of order α ≥ 0, of a function f ∈ Cμ , μ ≥ −1, is defined as 1 J f (t) = (α)

t (t − s)α−1 f (t)ds

α

J 0 f (t) = f (t),

(α > 0),

(2)

0

where (α) is well-known gamma function. Some of the properties of the operator J α , which we will need here, are as follows: For f ∈ Cμ , μ ≥ −1, α, β ≥ 0 and γ ≥ −1: 1. J α J β f (t) = J α+β f (t); 2. J α J β f (t) = J β J α f (t); (γ +1) α+γ 3. J α t γ = (α+γ t . +1)

Series Solutions of Systems of Nonlinear Fractional Differential

191

Definition 3 The fractional derivative (D α ) of f (t), in the Caputo sense is defined as 1 D f (t) = (n − α)

t (t − s)n−α−1 f (n) (t)ds,

α

(3)

0

n . for n − 1 < α < n, n ∈ N, t > 0, f ∈ C−1

The following are two basic properties of the Caputo fractional derivative [41]: n 1. Let f ∈ C−1 , n ∈ N . Than D α f, 0 ≤ α ≤ n is well defined and D α f ∈ C−1 . 2. Let n − 1 ≤ α ≤ n, n ∈ N and f ∈ Cμn , μ ≥ −1. Then

(J α D α )f (t) = f (t) −

n−1 

f (k) (0+ )

k=0

tk . k!

(4)

3 Solution Approach by HAM In HAM [17], system (1) is first written in the form, FiD [yi (t)] = 0,

i = 1, 2, . . . , n,

where FiD are a nonlinear fractional differential operators, t denotes the independent variable and yi (t) are an unknown function. By means of generalizing the traditional homotopy method, Liao [17] constructs the so-called zeroth-order deformation equation (1 − q)Li [φi (t; q) − yi,0 (t)] = qFiD [φi (t; q)],

(5)

where q ∈ [0, 1] is an embedding parameter,  is an auxiliary parameter, Li are an auxiliary linear operator, yi,0 (t) are initial guesses of yi (t) and φi (t; q) are unknown functions. It is important to note that, one has great freedom to choose auxiliary objects such as  and Li in HAM. We note that, in the frame of HAM, the solution yi (t) can be represented by many different base functions such as the polynomial functions, exponential functions, rational functions etc. Obviously, when q = 0 and q = 1, both φi (t; 0) = yi,0 (t)

and

φi (t; 1) = yi (t)

hold. Thus as q increases from 0 to 1, the solution φi (t; q) varies from the initial guess yi,0 (t) to the solution yi (t). Expanding φi (t; q) in Taylor series with respect to q, one has φi (t; q) = yi,0 (t) +

+∞ 

yi,m (t)q m ,

(6)

m=1

where  1 ∂ m φi (t; q)  . yi,m (t) = m! ∂q m q=0

(7)

192

A.S. Bataineh et al.

If the auxiliary linear operators, the initial guesses and the auxiliary parameters  are so properly chosen, then the series (6) converges at q = 1 and φi (t; 1) = yi,0 (t) +

+∞ 

yi,m (t),

(8)

m=1

which must be one of the solutions of the original nonlinear equations, as proved by Liao [17]. Note that as  = −1, (5) becomes (1 − q)Li [φi (t; q) − yi,0 (t)] + qFiD [φi (t; q)] = 0,

(9)

which is used mostly in the HPM [19, 20]. According to (7), the governing equation can be deduced from the zeroth-order deformation equation (5). Define the vector yi,n (t) = {yi,0 (t), yi,1 (t), . . . , yi,n (t)}. Differentiating (5) m times with respect to the embedding parameter q and then setting q = 0 and finally dividing them by m!, we have the so-called mth-order deformation equation Li [yi,m (t) − χm yi,m−1 (t)] = Ri,m (yi,m−1 ),

(10)

where Ri,m (yi,m−1 ) =

 ∂ m−1 FiD [φi (t; q)]  1  , (m − 1)! ∂q m−1 q=0

(11)

and  χm =

0, 1,

n ≤ 1, n > 1.

It should be emphasized that yi,m (t) (m ≥ 1) is governed by the linear equation (10) with the linear boundary conditions that come from the original problem, which can be easily solved by computer algebra systems such as Maple and Mathematica. The series solution (8) contains the auxiliary parameter  which was suggested to be called the convergence-control parameter [18]. The validity of the method is based on such an assumption that the series (6) converges at q = 1. It is the the convergence-control parameter  which ensures that this assumption can be satisfied. In general, by means of the so-called -curve, it is straightforward to choose a proper value of  which ensures that the solution series is convergent [18].

4 Illustrative Examples To demonstrate the effectiveness of the HAM algorithm discussed above, tow examples of nonlinear systems of FDEs will be studied.

Series Solutions of Systems of Nonlinear Fractional Differential

193

4.1 Example 1 We consider the following nonlinear system of FDEs: D α y1 = y12 + y2 ,

(12)

D y2 = y2 cos y1 ,

(13)

β

subject to the initial conditions y1 (0) = 0,

y2 (0) = 1,

(14)

where α, β ∈ (0, 1). We note that system (12)–(13) has unique solution as given in [42]. To solve (12)–(14) by means of the HAM, we choose the initial guess approximations y1,0 (t) = 0,

y2,0 (t) = 1,

(15)

L2 = D β ϕ2 (t; q),

(16)

and the linear operators L1 = D α ϕ1 (t; q), with the property Li [Ai ] = 0,

(17)

where Ai (i = 1, 2) are constants of integrations. Furthermore, (12)–(13) suggest that we define a system of nonlinear fractional differential operators as F1D (ϕi (t; q)) = D α ϕ1 (t; q) − ϕ12 (t; q) − ϕ2 (t; q),

(18)

F2D (ϕi (t; q)) = D β ϕ2 (t; q) − ϕ2 (t; q) cos[ϕ1 (t; q)].

(19)

According to (5), we construct the zeroth-order deformation equation (1 − q)Li [φi (t; q) − yi,0 (t)] = qFiD [φi (t; q)],

(20)

and the mth-order deformation equation Li [yi,m (t) − χm yi,m−1 (t)] = Ri,m (yi,m−1 ),

(21)

with the initial conditions yi,m (0) = 0, where R1,m (yi,m ) = D α y1,m−1 +

m−1 

y1,j y1,m−1−j + y2,m−1 ,

j =0

R2,1 (yi,0 ) = D β y2,0 + y1,0 cos[y2,0 ], R2,2 (yi,1 ) = D β y2,1 − 2y1,1 y2,0 sin[y1,0 ] + y2,1 cos[y1,0 ],

(22)

194

A.S. Bataineh et al.

1 2 R2,3 (yi,2 ) = D β y2,2 + y1,1 y2,0 cos[y1,0 ] − y1,2 y2,0 sin[y1,0 ] 2 − y1,1 y2,1 sin[y1,0 ] + y2,1 cos[y1,0 ], .. . Now, the solution of the mth-order deformation (10) for m ≥ 1 becomes y1,m (t) = χm y1,m−1 (t) + J α [R1,m (yi,m−1 )],

(23)

y2,m (t) = χm y2,m−1 (t) + J β [R2,m (yi,0 , yi,1 , yi,2 , . . .)].

(24)

We now successfully obtain by the HAM y1,1 = −

t α , (α + 1)

y1,2 = −

2 t α 2 t α+β t α − + , (α + 1) (α + 1) (α + β + 1)

y2,1 = −

t β , (β + 1)

y2,2 = −

2 t β 2 t 2β t β − + , (β + 1) (β + 1) (2β + 1)

etc. Then the series solutions expression by the HAM can be written in the form yi (t) = yi,0 (t) + yi,1 (t) + yi,2 (t) + · · · .

(25)

The proper value of  that ensures convergence can be obtained from the -curves [18]. Figure 1 shows samples of the -curves obtained from the 5th-order HAM approximation solutions for (12)–(14) for the case α = 0.5 and β = 0.8. Based on this figure, the valid regions of  correspond to the line segments nearly parallel to the horizontal axis [18]. In the case presented in Fig. 1, we can take  in the range −1.4 <  < −0.4. If we take  = −1, then we obtain the ADM solution [6]. In Fig. 2 we plot the numerical results for the 5th-order HAM approximation solutions yi (t) (i = 1, 2) for various values of , α, β and γ . The results presented in Fig. 2 clearly show the good accuracy of HAM. Fig. 1 The -curves of yi (0.1) (i = 1, 2) given by (25) when α = 0.5 and β = 0.8: 5th-order approximation of yi (0.1)

Series Solutions of Systems of Nonlinear Fractional Differential

195

Fig. 2 The numerical results for yi (t) (i = 1, 2) by means of 5-term HAM solution (25) for different values of  when α = 0.5 and β = 0.3

4.2 Example 2 Let us consider the following nonlinear system of FDEs: D α y1 = 2y22 ,

(26)

D β y2 = ty1 ,

(27)

D y3 = y2 y3 ,

(28)

γ

subject to the initial conditions y1 (0) = 0,

y2 (0) = 1,

y3 (0) = 1,

(29)

where α, β and γ ∈ (0, 1]. In order to solve (26)–(29) by means of the HAM, we choose the initial guess approximations y1,0 (t) = 0,

y2,0 (t) = 1,

y3,0 (t) = 1,

(30)

and the linear operators L1 = D α ϕ1 (t; q),

L2 = D β ϕ2 (t; q), Li [Ai ] = 0,

L2 = D γ ϕ3 (t; q),

(31) (32)

with the property (32), where Ai (i = 1, 2, 3) are constants of integrations. Similar to the previous example and by using (30) and (31) with property (17) we successfully obtain by the HAM y1,1 = −

2t α ,  (α + 1)

y1,2 = −

2( + 1)t α , (α + 1)

y2,1 = 0, y2,2 =

22 t α+β+1 (α + 2) , (α + 1)(α + β + 2)

196

A.S. Bataineh et al.

y3,1 = −

t γ ,  (γ + 1)

y3,2 = −

2 t γ 2 t 2γ t γ − + ,  (γ + 1)  (γ + 1)  (2γ + 1)

etc. Then the series solutions expression by the HAM can be written in the form yi (t) = yi,0 (t) + yi,1 (t) + yi,2 (t) + · · · . Fig. 3 The -curve of yi (0.1) (i = 1, 2, 3) given by (33) when α = 0.5, β = 0.8 and γ = 0.7: 5th-order approximation of yi (0.1)

Fig. 4 The numerical results for yi (t) (i = 1, 2, 3) by means of the 5th-order HAM approximation solution (33) for different values of  when α = 0.5, β = 0.4 and γ = 0.3

Fig. 5 The numerical results for yi (t) (i = 1, 2, 3) by means of the 5th-order HAM approximation solution (33) for different values of  when α = β = γ = 1

(33)

Series Solutions of Systems of Nonlinear Fractional Differential

197

Figure 3 shows the -curves obtained from the 5th-order HAM approximation solutions for (26)–(29) for the case α = 0.5, β = 0.8 and γ = 0.7. From this figure, the valid values of  fall in the range −1.4 <  < −0.4. Again, if we take  = −1, then we obtain the ADM solution [6]. In Figs. 4 and 5 we plot the numerical results for the 5th-order HAM approximation solutions yi (t) (i = 1, 2, 3) for various values of , α, β and γ . Again, the results presented in Figs. 4 and 5 clearly show the good accuracy of HAM.

5 Conclusions In this paper, the HAM was applied to systems of nonlinear fractional differential equations. The results of the test examples show that the ADM results are special cases of HAM. The HAM overcomes the difficulties arising in finding the Adomian polynomials and yields a family of solutions. The HAM is a promising tool for today’s numerous applications in physics, engineering and other fields. Acknowledgements The authors would like to acknowledge the financial support received from the MOSTI Sciencefund grants: 04-01-02-SF0177 and the SAGA Grant STGL-011-2006 (P24c). The referee’s comments which led to an improvement to the paper are gratefully acknowledged.

References 1. Oldham, K.B., Spanier, J.: The Fractional Calculus. Academic Press, New York (1974) 2. Miller, K.S., Ross, B.: An Introduction to the Fractional Calculus and Fractional Differential Equations. Wiley, New York (1993) 3. Podlubny, I.: Fractional Differential Equations. Academic Press, New York (1999) 4. Li, C.P., Peng, G.J.: Chaos in Chen’s system with fractional order. Chaos Solitons Fractals 22(2), 443– 450 (2004) 5. Wang, Y.H., Li, C.P.: Does the fractional Brusselator with efficient dimension less than 1 have a limit cycle? Phys. Lett. A 363(5–6), 414–419 (2007) 6. Jafari, H., Daftardar-Gejji, V.: Solving a system of nonlinear fractional differential equations using Adomian decomposition. J. Comput. Appl. Math. 196, 644–651 (2006) 7. Shawagfeh, N.T.: Analytical approximate solutions for nonlinear fractional differential equations. Appl. Math. Comput. 131, 517–29 (2002) 8. Abdulaziz, O., Hashim, I., Chowdhury, M.S.H., Zulkifle, A.K.: Assessment of decomposition method for linear and nonlinear fractional differential equations. Far East J. Appl. Math. 28, 95–112 (2007) 9. Abdulaziz, O., Hashim, I., Ismail, E.S.: Approximate analytical solutions to fractional modified KdV equations. Math. Comput. Model. (2008). doi:10.1016/j.mcm.2008.01.005 10. Abdulaziz, O., Hashim, I., Momani, S.: Solving systems of fractional differential equations by homotopy-perturbation method. Phys. Lett. A 372, 451–459 (2008) 11. Abdulaziz, O., Hashim, I., Momani, S.: Application of homotopy-perturbation method to fractional IVPs. J. Comput. Appl. Math. 216, 574–584 (2008) 12. Momani, S., Odibat, Z.: Homotopy perturbation method for nonlinear partial differential equations of fractional order. Phys. Lett. A 365, 345–350 (2007) 13. Odibat, Z., Momani, S.: Modified homotopy perturbation method: application to quadratic Riccati differential equation of fractional order. Chaos Solitons Fractals 36, 167–174 (2008) 14. Odibat, Z., Momani, S.: Application of variational iteration method to nonlinear differential equation of fractional order. Int. J. Nonlinear Sci. Numer. Simul. 1, 271–279 (2006) 15. Momani, S., Odibat, Z.: Numerical comparison of methods for solving linear differential equations of fractional order. Chaos Solitons Fractals 31, 1248–1255 (2007) 16. Liao, S.J.: The proposed homotopy analysis techniques for the solution of nonlinear problems. Ph.D. dissertation, Shanghai Jiao Tong University, Shanghai (1992) 17. Liao, S.J.: Beyond Perturbation: Introduction to the Homotopy Analysis Method. CRC/Chapman and Hall, Boca Raton (2003)

198

A.S. Bataineh et al.

18. Liao, S.J.: Notes on the homotopy analysis method: some definitions and theorems. Commun. Nonlinear Sci. Numer. Simul. (2008). doi:10.1016/j.cnsns.2008.04.013 19. He, J.H.: An approximate solution technique depending upon an artificial parameter. Commun. Nonlinear Sci. Numer. Simul. 3(2), 92–97 (1998) 20. He, J.H.: Newton-like iteration method for solving algebraic equations. Commun. Nonlinear Sci. Numer. Simul. 3, 106–109 (1998) 21. Abbasbandy, S.: The application of the homotopy analysis method to nonlinear equations arising in heat transfer. Phys. Lett. A 360, 109–113 (2006) 22. Sajid, M., Hayat, T.: Comparison of HAM and HPM methods for nonlinear heat conduction and convection equations. Nonlinear Anal. Real World Appl. (2007). doi:10.1016/j.nonrwa.2007.08.007 23. Hayat, T., Sajid, M.: On analytic solution for thin film flow of a fourth grade fluid down a vertical cylinder. Phys. Lett. A 361, 316–322 (2007) 24. Liao, S.J.: An approximate solution technique which does not depend upon small parameters (Part 2): an application in fluid mechanics. Int. J. Nonlinear Mech. 32, 815–822 (1997) 25. Liao, S.J.: An explicit totally analytic approximation of Blasius viscous flow problems. Int. J. Nonlinear Mech. 34, 759–778 (1999) 26. Liao, S.J.: On the homotopy analysis method for nonlinear problems. Appl. Math. Comput. 147, 499–513 (2004) 27. Liao, S.J., Pop, I.: Explicit analytic solution for similarity boundary layer equations. Int. J. Heat Mass Transfer 47, 75–78 (2004) 28. Liao, S.J.: Comparison between the homotopy analysis method and homotopy perturbation method. Appl. Math. Comput. 169, 1186–1194 (2005) 29. Liao, S.J.: A new branch of solutions of boundary-layer flows over an impermeable stretched plate. Int. J. Heat Mass Transfer 48, 2529–3259 (2005) 30. Ayub, M., Rasheed, A., Hayat, T.: Exact flow of a third grade fluid past a porous plate using homotopy analysis method. Int. J. Eng. Sci. 41, 2091–2103 (2003) 31. Hayat, T., Khan, M., Asghar, S.: Homotopy analysis of MHD flows of an Oldroyd 8-constant fluid. Acta Mech. 168, 213–232 (2004) 32. Hayat, T., Khan, M.: Homotopy solutions for a generalized second-grade fluid past a porous plate. Nonlinear Dyn. 42, 395–405 (2005) 33. Abbasbandy, S.: The application of homotopy analysis method to solve a generalized Hirota-Satsuma coupled KdV equation. Phys. Lett. A 15, 1–6 (2006) 34. Xu, H., Liao, S.J.: A series solution of the unsteady Von Karman swirling viscous flows. Acta Appl. Math. 94, 215–231 (2006) 35. Bataineh, A.S., Noorani, M.S.M., Hashim, I.: The homotopy analysis method for Cauchy reactiondiffusion problems. Phys. Lett. A 372, 613–618 (2008) 36. Bataineh, A.S., Noorani, M.S.M., Hashim, I.: Solving systems of ODEs by homotopy analysis method. Commun. Nonlinear Sci. Numer. Simul. 13, 2060–2070 (2008) 37. Song, L., Zhang, H.: Application of homotopy analysis method to fractional KdV-Burgers-Kuramoto equation. Phys. Lett. A 367, 88–94 (2007) 38. Hashim, I., Abdulaziz, O., Momani, S.: Homotopy analysis method for fractional IVPs. Commun. Nonlinear Sci. Numer. Simul. (2007). doi:10.1016/j.cnsns.2007.09.014 39. Cang, J., Tan, Y., Xu, H., Liao, S.J.: Series solutions of non-linear Riccati differential equations with fractional order. Chaos Solitons Fractals (2007). doi:10.1016/j.chaos.2007.04.018 40. Xu, H., Cang, J.: Analysis of a time fractional wave-like equation with the homotopy analysis method. Phys. Lett. A 372, 1250–1255 (2008) 41. Gorenflo, R., Mainardi, F.: Fractional calculus: Integral and differential equations of fractional order. In: Carpinteri, A., Mainardi, F. (eds.) Fractals and Fractional Calculus in Continuum Mechanics. Springer, New York (1997) 42. Daftardar-Gejji, V., Babakhani, A.: Analysis of a system of fractional differential equations. J. Math. Anal. Appl. 293, 511–522 (2004)