INTRODUCTION
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 fluiddynamic 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 semiinfinite 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 RiemannLiouville and Caputo. Each definition uses RiemannLiouville fractional integration and derivatives of whole order. The difference between the two definitions is in the order of evaluation. RiemannLiouville fractional integration of order α is defined as:
The next two Equations define RiemannLiouville and Caputo fractional derivatives of order α, respectively,
where, m1<α<=m and m ∈0 N. For now, the Caputo fractional derivative
will be denoted by D_{*}^{α} to maintain a clear distinction
with the RiemannLiouville 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 RiemannLiouville 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
RiemannLiouville 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
roundoff 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.
FRACTIONAL DIFFERENCE METHOD
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ünwaldLetnikov sense as:
where, [t] means the integer part of t and h is the step size.
The definition of operator in the GrünwaldLetnikov sense (5)
is equivalent to the definition of operator in the RiemannLiouville sense (3)
(Podlubny, 1999). Nevertheless the GrünwaldLetnikov 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, t_{m} = mh, y_{m} = y (t_{m}) and c^{α}_{j} are GrünwaldLetnikov coefficients defined as:
Using the recurrence relationship (Podlubny, 1999).
We can compute the coefficients in a simple way. For j = 1 we have c_{1}^{α} = –α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).
DECOMPOSITION METHOD
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, A_{n} are the socalled 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 A_{n} (t) can be calculated for all forms of nonlinearity according to specific algorithms constructed by Adomian (1994). The general form of formula for A_{n} 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).
NUMERICAL EXPERIMENTS
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 y_{5},...,y_{18} and y_{19} 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 wellknown
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 y_{exact} 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 semiinfinite
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 t^{1/2} = x; then,
Calculating the [8/8] Padé approximants and recalling that x = t^{1/2}, we get
The obtained numerical results for both methods are shown in Table
2. From these results we conclude that the halforder 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 timefractional derivative.

Fig. 1: 
[8/8] Pade approximants of y (t): (—) α = 1, (....)
α = 0.75, (−−−α = 0.5) 
CONCLUDING REMARKS
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).