Subscribe Now Subscribe Today
Research Article

Numerical Approximations of a Dynamic System Containing Fractional Derivatives

Shaher Momani, Omar K. Jaradat and Rabha Ibrahim
Facebook Twitter Digg Reddit Linkedin StumbleUpon E-mail

This study presents numerical methods-fractional difference and Adomian decomposition-for solution of a dynamic system containing fractional derivative of order α, 0<α<=1. The fractional derivative is described in the Caputo sense. The Adomian decomposition method provides the solution in the form of a convergent power series with easily computable components. Then the diagonal Pade approximants are effectively used in the analysis to capture the essential behavior of the solution. The presented schemes are introduced in a general way so that they can be used in applied sciences.

Related Articles in ASCI
Similar Articles in this Journal
Search in Google Scholar
View Citation
Report Citation

  How to cite this article:

Shaher Momani, Omar K. Jaradat and Rabha Ibrahim , 2008. Numerical Approximations of a Dynamic System Containing Fractional Derivatives. Journal of Applied Sciences, 8: 1079-1084.

DOI: 10.3923/jas.2008.1079.1084



Although fractional derivatives have a long mathematical history, for many years they were not used in physics. One possible explanation of such unpopularity could be that there are multiple nonequivalent definitions of fractional derivatives (Podlubny, 1999). Another difficult is that fractional derivatives have no evident geometrical interpretation because of their nonlocal character (Podlubny, 2002). However, during the last ten years fractional calculus starts to attract much more attention of physicists and mathematicians. It was found that various, especially interdisciplinary applications can be elegantly modeled with the help of the fractional derivatives. For example, the nonlinear oscillation of earthquake can be modeled with fractional derivatives (He, 1998a) and the fluid-dynamic traffic model with fractional derivatives (He, 1998b) can eliminate the deficiency arising from the assumption of continuum traffic flow. Based on experimental data fractional partial differential equations for seepage flow in porous media are suggested in (He, 1998c) and differential equations with fractional order have recently proved to be valuable tools to the modeling of many physical phenomena (Grigorenko and Grigorenko, 2003; Podlubny, 1999). A review of some applications of fractional derivatives in continuum and statistical mechanics is given by (Mainardi, 1997).

In this study, we present numerical and analytical solutions for the nonlinear fractional ordinary differential Eq:


where, β is a real constant and α is a parameter describing the order of the fractional derivative. The general response expression contains a parameter describing the order of the fractional derivative that can be varied to obtain various responses. For the special case α = 0.5, (1) describes the cooling of a semi-infinite body by radiation and the initial value problem has been solved numerically using the fractional difference method (Podlubny, 1999).

There are several definitions of a fractional derivative of order α>0 (Oldham and Spanier, 1974; Podlubny, 1999) The two most commonly used definitions are the Riemann-Liouville and Caputo. Each definition uses Riemann-Liouville fractional integration and derivatives of whole order. The difference between the two definitions is in the order of evaluation. Riemann-Liouville fractional integration of order α is defined as:


The next two Equations define Riemann-Liouville and Caputo fractional derivatives of order α, respectively,



where, m-1<α<=m and m ∈0 N. For now, the Caputo fractional derivative will be denoted by D*α to maintain a clear distinction with the Riemann-Liouville fractional derivative. The Caputo fractional derivative first computes an ordinary derivative followed by a fractional integral to achieve the desired order of fractional derivative. The Riemann-Liouville fractional derivative is computed in the reverse order. We have chosen to use the Caputo fractional derivative because it allows traditional initial and boundary conditions to be included in the formulation of the problem, but for homogeneous initial condition assumption, these two operators coincide. For more details on the geometric and physical interpretation for fractional derivatives of both the Riemann-Liouville and Caputo types see (Podlubny, 2002).

In this study, we will use the Fractional Difference Method (FDM) (Podlubny, 1999) to find an approximate numerical solution of the nonlinear fractional ordinary differential Eq. 1. Also we will use the Adomian Decomposition Method (ADM) (Adomian, 1988, 1994) to provide analytic, verifiable, rapidly convergence approximation to Eq. 1. Unlike the fractional difference method, the decomposition method provides us with numerical solution without discretization of the given equation and therefore, it is not affected by computation round-off errors and one is not faced with necessity of large computer memory and time. The method introduces the solution in the form of a convergent power series with elegantly computable terms. Furthermore, the behavior of the model can be formally determined by using the Padé approximants of the series obtained. The combination of the series solution with the Padé approximants was successfully implemented in (Boyd, 1997; Momani, 2004; Wazwaz, 1999) and proved to be effective and promising. Moreover, it was found that the most accurate approximant is usually the diagonal approximant. The essential behavior of the solution y (t) of Eq. 1 will be addressed by using several diagonal Padé approximants of different degrees.


Here, we introduce another definition of the fractional derivative. With regard to (Oldham and Spanier, 1974; Podlubny, 1999) we define the fractional derivative in the Grünwald-Letnikov sense as:


where, [t] means the integer part of t and h is the step size.

The definition of operator in the Grünwald-Letnikov sense (5) is equivalent to the definition of operator in the Riemann-Liouville sense (3) (Podlubny, 1999). Nevertheless the Grünwald-Letnikov operator is more flexible and most straightforward in numerical calculations. Approximating the fractional derivative in 1 by 5, we obtain the following approximation for Eq. 1:


where, tm = mh, ym = y (tm) and cαj are Grünwald-Letnikov coefficients defined as:


Using the recurrence relationship (Podlubny, 1999).


We can compute the coefficients in a simple way. For j = 1 we have c1α = –αh.

Approximation (6) leads to the numerical solution algorithm described by:


This algorithm is simple for computational performance for all values of α, 0<α<=1. For details about the fractional difference method and its applications for solving fractional differential equations, we refer the reader to (Podlubny, 1999).


The decomposition method requires that the fractional differential Eq. 1 be expressed in terms of operator from as:


where, the fractional differential operator D*α is defined as in Eq. 4 denoted by:

Applying the operator Jα, the inverse of the operator D*α, to both sides of Eq. 10 and using the initial condition lead to:


The Adomian`s decomposition method (Adomian, 1988, 1994) suggests the solution y (t) be decomposed by the infinite series of components:


and the nonlinear function in Eq. 11 is decomposed as follows:


where, An are the so-called the Adomian polynomials.

Substitution the decomposition series 12 and 13 into both sides of 11 gives:


From Eq. 14, the iterates are determined by the following recursive way:


The Adomian polynomial An (t) can be calculated for all forms of nonlinearity according to specific algorithms constructed by Adomian (1994). The general form of formula for An Adomian polynomials is:


This formula is easy to compute by using Mathematica software or by setting a computer code to get as many polynomials as we need in the calculation of the numerical as well as explicit solutions.

Finally, we approximate the solution y (t) by the truncated series:


However, in many cases the exact solution in a closed form may be obtained. Moreover, the decomposition series solutions are generally converge very rapidly. The convergence of the decomposition series has investigated by several authors. The theoretical treatment of convergence of the decomposition method has been considered in the literature (Abbaoui and Cherruault, 1996a, b; Cherruault, 1989; Cherruault and Adomian, 1993; Répaci, 1990). They obtained some results about the speed of convergence of this method. In recent study of Abbaoui and Cherruault (1996a) have proposed a new approach of convergence of the decomposition series. The authors have given a new condition for obtaining convergence of the decomposition series to the classical presentation of the ADM in (Abbaoui and Cherruault, 1996b).


Here, we use the ADM to find an approximate solution of the fractional model Eq. 1 together with the initial condition y (0) = 0 and A = 1. Equation 1 is also solved numerically and the corresponding results are compared with the Adomian solution. The numerical method adopted in this study was the fractional difference method and the step size h is chosen to be 0.01.

To calculate the terms of the decomposition series 12 for y (t), we substitute the initial condition and the corresponding Adomian polynomials into 15 and using Mathematica consequently, we obtain


The components y5,...,y18 and y19 were also determined and will be used, but for brevity not listed. In this matter twenty components of the decomposition series (12) were obtained of which y (t) was evaluated to have the following expansion:


Our aim is to study the mathematical behaviour of the solution y (t) for different values of α. It was formally shown by Boyd (1997), Momani (2004) and Wazwaz (1999) that this goal can be achieved by forming Padé approximants (Baker, 1975) which have the advantage of manipulating the polynomial approximation into a rational function to gain more information about y (t). It is well-known that Padé approximants will converge on the entire real axis (Boyd, 1997) if y (t) is free of singularities on the real axis. It is of interest to note that Padé approximants give results with no greater error bounds (Burden and Faires, 1993) than approximation by polynomials. More importantly, the diagonal approximant is the most accurate approximant, therefore we will construct only the diagonal approximants in the following discussions. The computational work can be performed by using any manipulation language such as Maple and Mathematica.

To consider the behavior of solution of different values of α, we will take advantage of the explicit formula (19) available for 0<α<=1 and consider the following two special cases:

First order case: Setting α = 1 in Eq. 19, we obtain the approximate solution in a series form as:


The [4/4] Padé approximants gives:


Table 1 shows the values of y (t) by using the numerical solution obtained by the fractional difference method and by using the [4/4], [8/8] and [10/10] Padé approximants for specific value of α = 1. The exact values of yexact were evaluated using the exact solution of (1) given by:


and the values for both methods are in good agreement with the exact values.

An important conclusion can be drawn is that the error decreases dramatically with the increase of the degree of the Padé approximants.

Fractional order case: In this case we will examine the fractional nonlinear Eq. 1 when α = 1/2. It describes the cooling of a semi-infinite body by radiation (Podlubny, 1999). Setting α = 1/2 in Eq. 19 gives:

Table 1: Numerical values when α = 1 using the two methods

Table 2: Numerical values when α = 0.5 using the two methods


For simplicity, let t1/2 = x; then,


Calculating the [8/8] Padé approximants and recalling that x = t1/2, we get


The obtained numerical results for both methods are shown in Table 2. From these results we conclude that the half-order derivative system exhibits fast increase in the beginning and slow increase later. Similar behavior was observed in (Podlubny, 1999).

Figure 1 shows the [8/8] Padé approximants of y (t) for α = 1, α = 0.75 and α = 0.5. It can be seen from Fig. 1 that the solution continuously depends on the time-fractional derivative.

Fig. 1: [8/8] Pade approximants of y (t): (—) α = 1, (....) α = 0.75, (−−−α = 0.5)


The fundamental goal of this work has been to construct an approximation to the solution of a dynamic system containing fractional derivatives. The goal has been achieved by using the fractional difference method and the decomposition method. The methods were used in a direct way without using linearization, perturbation or restrictive assumptions.

There are three important points to make here. First, unlike the fractional difference techniques, where y (t) is defined at grid points only, the solution obtained by the decomposition method is given in a fractional series form. Second, the high agreement of the approximation of y (t) between the methods used at this study is clear and remarkable. Finally, the behavior of the model, in that it exhibits fast increase in the beginning and slow increase later, can be formally determined by using the Padé approximants. The combining of the series obtained with the Pade approximants provides a successful tool and promising scheme for identical applications.

In future researches we will focus on the implementation of a new analytical technique, He`s variational iteration method for solving the discussed problem. This method has many merits and has much advantages over the Adomian method. It can be introduced to overcome the difficulty arising in calculating Adomian polynomials. Applications of the variational iteration method can be found in Refs. (He, 1998a, b, 1999, 2000; He et al., 2004).

1:  Abbaoui, K. and Y. Cherruault, 1995. New ideas for proving convergence of decomposition methods. Comput. Math. Appl., 29: 103-108.
Direct Link  |  

2:  Abbaoui, K. and Y. Cherruault, 1994. Convergence of Adomian's method applied to differential equations. Comput. Math. Applic., 28: 103-109.
CrossRef  |  Direct Link  |  

3:  Adomian, G., 1988. A review of decomposition method in applied mathematics. J. Math. Anal. Appl., 135: 501-544.
CrossRef  |  

4:  Adomian, G., 1994. Solving Frontier Problems of Physics: The Decomposition Method. 1st Edn., Kluwer Academic, Boston, Pages: 352.

5:  Baker, G.A., 1975. Essentials of Padé Approximants. 1st Edn., Academic Press, London.

6:  Boyd, J.P., 1997. Pade approximant algorithm for solving nonlinear ordinary differential equation boundary value problems on an unbounded domain. Comput. Phys., 11: 299-303.
CrossRef  |  Direct Link  |  

7:  Burden, R.L. and J.D. Faires, 1993. Numerical Analysis. 1st Edn., P.W.S., Boston.

8:  Cherruault, Y., 1989. Convergence of Adomian's method. Kybernetes, 8: 31-38.
CrossRef  |  

9:  Cherruault, Y. and G. Adomian, 1993. Decomposition methods: A new proof of convergence. Math. Comput. Modelling, 18: 103-106.
CrossRef  |  

10:  Grigorenko, I. and E. Grigorenko, 2003. Chaotic dynamics of the fractional Lorenz system. Phys. Rev. Lett., 91: 034101-034104.
Direct Link  |  

11:  He, J.H., 1998. Nonlinear oscillation with fractional derivative and its applications. Proceedings of the International Conference on Vibrating Engineering, (VE'98), Dalian, China, pp: 288-291.

12:  He, J.H., 1998. Approximate solution of nonlinear differential equations with convolution product nonlinearities. Comput. Methods Applied Mech. Eng., 167: 69-73.
CrossRef  |  Direct Link  |  

13:  He, J.H., 1998. Approximate analytical solution for seepage flow with fractional derivatives in porous media. J. Comput. Math. Applied Mech. Eng., 167: 57-68.
CrossRef  |  Direct Link  |  

14:  He, J.H., 1999. Variational iteration method-a kind of non-linear analytical technique: Some examples. Int. J. Nonlinear Mech., 34: 699-708.
CrossRef  |  Direct Link  |  

15:  He, J.H., 2000. Variational iteration method for autonomous ordinary differential systems. Applied Math. Comput., 114: 115-123.
CrossRef  |  Direct Link  |  

16:  He, J.H., Y.Q. Wan and Q. Guo, 2004. An iteration formulation for normalized diode characteristics. Int. J. Circuit Theor. Applied, 32: 629-632.
Direct Link  |  

17:  Mainardi, F., 1997. Fractional Calculus: Some Basic Problems in Continuum and Statistical Mechanics. In: Fractals and Fractional Calculus in Continuum Mechanics, Carpinteri, A. and F. Mainardi (Eds.). Springer-Verlag, New York, pp: 291-348.

18:  Momani, S., 2004. Analytical approximate solutions of nonlinear oscillators by the modified decomposition method. Int. J. Modern Phys. C., 15: 967-979.
Direct Link  |  

19:  Oldham, K.B. and J. Spanier, 1974. The Fractional Calculus. 1st Edn., Academic Press, New York.

20:  Podlubny, I., 1999. Fractional Differential Equations. 1st Edn., Academic Press, New York.

21:  Podlubny, I., 2002. Geometric and physical interpretation of fractional integration and fractional differentiation. Fract. Calculus Applied Anal., 5: 367-386.
Direct Link  |  

22:  Repaci, A., 1990. Nonlinear dynamical systems: On the accuracy of Adomian's decomposition method. Applied Math. Lett., 3: 35-39.
CrossRef  |  Direct Link  |  

23:  Wazwaz, A.M., 1999. Analytical approximations and Pade approximants for Volterra's population model. Applied Math. Comput., 100: 13-25.
CrossRef  |  

©  2021 Science Alert. All Rights Reserved