Int. J. Adv. Appl. Math. and Mech.

2(1) (2014) 83 - 91
ISSN: 2347-2529
Available online at www.ijaamm.com

International Journal of Advances in Applied Mathematics and Mechanics

Numerical solution of Bagley-Torvik equation using

Chebyshev wavelet operational matrix of fractional
derivative
Research Article

Fakhrodin Mohammadi ∗
Department of Mathematics, Hormozgan University, Bandarabbas, Iran

Received 04 July 2014; accepted (in revised version) 24 August 2014

Abstract: In this paper Chebyshev wavelet and their properties are employed for deriving Chebyshev wavelet operational
matrix of fractional derivatives and a general procedure for forming this matrix is introduced. Then Chebyshev
wavelet expansion along with this operational matrix are used for numerical solution of Bagley-Torvik boundary
value problems. The error analysis and convergence properties of the Chebyshev wavelet method are investi-
gated.

MSC: 34K37 • 42C40 • 65T60

Keywords: Chebyshev wavelet • Fractional derivatives • Operational matrix • Bagley-Torvik equation • Tau method

c 2014 IJAAMM all rights reserved.

1. Introduction
The idea of derivatives of noninteger order was initially appeared in a letter from Leibniz to L’Hospital in 1695.
For three centuries, studies of the theory of fractional order were mainly constraint to the field of pure theoreti-
cal mathematics, which were only useful for mathematicians. In the last several decades, many researchers found
that derivatives of non-integer order are very suitable for the description of various physical phenomena such as
damping laws, diffusion process, etc. These findings invoked the growing interest of studies of the fractional cal-
culus in various fields such as physics, chemistry and engineering. For these reasons we need reliable and efficient
techniques for the solution of fractional differential equations [1–5].
The Bagley-Torvik equation is a kind of fractional differential equation that apeare in the studies on behavior of real
material by use of fractional calculus [6, 7]. It has many applications in engineering and applied sciences fields,
for mor details see[3]. So, numerical solution of Bagley-Torvik fractional differential equation has attracted many
attention and has been studied by many authors. Several methods such as Adomian decomposition method [8, 9],
HeâĂŹs variational iteration method [10], Taylor collocation method [11] have been used to solve this fractional dif-
ferential equation. Diethelm [12] transformed this equation into first-order coupled fractional differential equation
and solved the problem with Adams predictor and corrector approach. Podlubny used successive approximation
method to solve the equation and recently applied the matrix approach to discretization of fractional derivatives for
the same problem [3, 13].
Wavelets theory is a relatively new and an emerging area in mathematical research. As a powerful tool, wavelets
have been extensively used in signal processing, numerical analysis, and many other areas. Wavelets permit the ac-
curate representation of a variety of functions and operators [14–16]. In this paper The Chebyshev wavelets are first
introduced, then by using shifted Chebyshev polynomial and their properties the operational matrix of derivative
and fractional derivative are derived. Then, applications of these operational matrices for solving fractional order
Bagley-Torvik boundary value problems are described. Illustrative examples are given to demonstrate the efficiency

∗ Corresponding author.
E-mail address: [email protected]

83
 Numerical solution of Bagley-Torvik equation using Chebyshev wavelet operational matrix of fractional derivative

and capability of the proposed method. The error analysis and convergence properties of the Chebyshev wavelet
method are investigated.
The article is organized as follows: We begin by introducing some necessary definitions and mathematical prelimi-
naries of the fractional calculus theory and Chebyshev polynomials. In Section 3 the Chebyshev wavelet and its oper-
ational matrix of derivatives and fractional derivative are derived and general procedure for forming these matrices
are introduced. Fractional order boundary value problems are introduced in Section 4 and then a method based
on Chebyshev wavelet and its operational matrices is establish for solving this fractional boundary value problems.
Numerical examples are included in section 5. Finally, a conclusion is given in Section 6.

2. Basic definitions
There are various definitions of fractional integration and differentiation, such as Grunwald-Letnikov’s definition
and Riemann-Liouville’s definition. The Riemann-Liouville derivative has certain disadvantages when trying to
model real world phenomena with fractional differential equations. Therefore, we shall introduce a modified frac-
tional differential operator D α proposed by Caputo in his work on the theory of viscoelasticity [3, 4].

Definition 2.1.
A real function f (t ), t > 0, is said to be in the space Cµ , µ ∈ R if there exists a real number p (> µ) and a function
f1 (t )C [0, ∞) such that f (t ) = t p f1 (t ), and it is said to be in the space Cµn , n ∈ N if f (n ) C .

Definition 2.2.
The Riemann-Liouville fractional integration operator of order α ≥ 0 of a function f ∈ Cµ , ≥ −1, is defined as:
d n f (t )
α = n ∈ N,

α dtn ,
D f (t ) = 1 dn
Rt f (τ) (1)
Γ (n−α) d t n 0 (t −τ)α−n+1
d τ, 0 ≤ n − 1 < α < n.

Definition 2.3.
The fractional derivative of order α > 0 in the Caputo sense is defined as
¨ d n f (t )
α d t n ,R α = n ∈ N,
D∗ f (t ) = 1 t f (n) (τ) (2)
Γ (n −α) 0 (t −τ)α−n+1 d τ, t > 0, 0 ≤ n − 1 < α < n.

The useful relation between the Riemann-Liouvill operator and Caputo operator is given by the following expression

n −1
X tk
I α D∗α f (t ) = f (t ) − f (k ) (0+ ) , t > 0, (n − 1 < α ≤ n ). (3)
k =0
k!

Also for the Caputo’s derivative we have D∗α c = 0, in which c is a constant and

n < dαe, n ∈ Z+

α n 0,
D∗ t = Γ (n +1) , (4)
Γ (n+1−α) t
n −α
, n ≥ dαe, n ∈ Z+

where we use the ceiling function dαe to denote the smallest integer greater than or equal to α. For more details
about fractional calculus and its properties see [3].

2.1. Shifted Chebyshev polynomials

The sequence of Chebyshev polynomials {Tn (t )}∞ n =0 which are orthogonal with respect to the weight function w (t ) =
p 1 defined on the interval [−1, 1] and they obey the following recurrence relation [20]
1−t 2

T0 (t ) = 1, T1 (t ) = t ,

Tn (t ) = 2Tn−1 (t ) − Tn −2 (t ), n ≥ 2.

for useing Chebyshev polynomials on the interval [0, 1] we define the so-called shifted Chebyshev polynomials by
introducing the change of variable t = 2x − 1. So, the shifted Chebyshev polynomials Tn∗ (x ) on [0,1] can be obtained
as

Tn∗ (x ) = Tn (2x − 1).

84
 Fakhrodin Mohammadi / Int. J. Adv. Appl. Math. and Mech. 2(1) (2014) 83 - 91

The orthogonality condition for shifted polynomials is

Z1 § πγm
T ∗ (x )Tn∗ (x ) 4 , m = n,
pm dx = (5)
1 − (2x − 1)2 0, m 6= n ,
0

where
2, m = 0,
§
γm =
1, m ≥ 1.

Lemma 2.1.
The analytic form of the shifted Chebyshev polynomial Tm∗ (x ) of degree m is given by
m
X
Tm∗ (x ) = a m i x m−i , (6)
i =0

where
m (2m − i − 1)!
a m i = (−1)i 22m−2i . (7)
(i )!(2m − 2i )!

Proof. See [17].

3. Chebyshev wavelets
In this section, we introduce Chebyshev wavelet and its operational matrix. Chebyshev wavelets ψnm (x ) =
ψ(k , n , m, x ) are defined on the interval [0, 1) by
 k +1
2 2 T˜m (2k x − (2n + 1)), 2nk ≤ x ≤ n2+1
ψn m (x ) = k
, (8)
0, o t he r w i se

where
p1 m =0

,
T˜m (x ) = qπ ,
2
π Tm (x ), m >0

and n = 0, 1, . . . , 2k − 1, k ∈ N, m ∈ Z M = {0, 1, 2, ..., M − 1} for a fixed positive integer M . The Chebyshev wavelets
{ψn m (x )|n = 0, 1, . . . , 2k − 1, m ∈ Z M } forms an orthonormal basis for L w 2
n
[0, 1] with respect to the weight function
wn (x ) = w (2k +1 x − (2n + 1)), in which w (x ) = p1−x
1
2
.
Any square integralable function f (x ) defined over [0, 1] may be expanded by Chebyshev wavelets as
∞ X
X ∞
f (x ) = cn m ψn m (x ), (9)
n =1 m=0

If the infinite series in (9) is truncated, then it can be written as

k
2X −1 M
X −1
f (x ) ' cn m ψn m (x ) = C T Ψ(x ), (10)
n =0 m=0

where C and Ψ(x ) are 2k M column vectors as

 T
C = c00 , . . . , c0(M −1) |c10 , . . . , c1(M −1) |, . . . , |c2k 0 , . . . , c2k (M −1) , (11)
 T
Ψ(x ) = ψ00 , . . . , ψ0(M −1) |ψ10 , . . . , ψ1(M −1) |, . . . , |ψ2k 0 , . . . , ψ2k (M −1) .

Lemma 3.1.
By using the shifted Chebyshev polynomials, any component Ψr (x ) of (11) can be written as:
v
k +1 t 2
Ψr (x ) = ψn m (x ) = 2 2 T˜m (2k x − n )χ[ nk , n+1 ],
πγm 2 2k

where r = n M + m + 1, m = 0, . . . , M − 1, n = 0, . . . , 2k − 1 and χ[ nk , n+1

k ]
is the characteristic function defined as:
2 2

1, x ∈ [ 2nk , n+1
2k ],
χ[ nk , n+1
k ]
(x ) =
2 2 0, t h e r w i s e .

85
 Numerical solution of Bagley-Torvik equation using Chebyshev wavelet operational matrix of fractional derivative

Now we present an useful theorem about operational matrix of derivative for Chebyshev wavelets.

Theorem 3.1.
Let Ψ(x ) be the Chebyshev wavelet vector defined in (11), and α > 0 (N − 1 < α ≤ N ) then we have
D∗α Ψ(x ) = D α Ψ(x ), (12)
α
where D is the (2 M ) × (2 M ) operational matrix of fractional derivative of order α, in the Caputo sense and its
k k

(p , q )-th component is
nM + 1 ≤ p ≤ nM + dαe,

α
0, Ç
[D ]p q = k +1 2 α (13)
2 2 πγm w mq , nM + dαe + 1 ≤ p ≤ (n + 1)M .
in which
m m−i  ‹
α
X X m −i
wm q = b j q am i (−1)m−i − j 2k j n m−i − j (14)
j
i =0 j =0

€ Š
b j q = D∗α (x j χ[ nk , n+1
k ]
), Ψq (x ) (15)
2 2 wn

and a m j is defined in (7).

Proof. For n M + 1 ≤ p ≤ nM + dαe, we have 0 ≤ m ≤ dαe − 1, and consequently

v
k +1 t 2
D∗α Ψp (x ) = Ψp (x ) = 2 2 D α T ∗ (2k x − n ) = 0, n = 0, . . . , 2k − 1, (16)
πγm ∗ m
for n M + dαe + 1 ≤ p ≤ (n + 1)M we have dαe ≤ m ≤ M − 1, so D∗α Ψp (x ) 6= 0. By using Lemma 3.1 we obtain
m
v
k +1 t 2 X
Ψp (x ) = 2 2 a m i (2k x − n )m −i χ[ nk , n+1 ]. (17)
πγm i =0 2 2k

this function is zero outside the interval χ[ 2nk , n2+1

k ]. So by substitute Newton expansion of (2 x − n )
k m−i
into (17) we
obtain
m m−i
v
2 XX
 ‹
k +1 t m −i
Ψp (x ) = 2 2 am i (−1)m−i − j 2k j n m−i − j x j χ[ nk , n+1 ], (18)
πγm j
i =0 j =0
2 2k

Now by apply D∗α

on both sides of (18) we have:
m m−i
v
2 XX
 ‹
α k +1 t m −i € Š
D∗ Ψp (x ) = 2 2 am i (−1)m−i − j 2k j n m−i − j D∗α x j χ[ nk , n+1 ] . (19)
πγm j
i =0 j =0
2 2 k

Eq. (19) can be rewritten as

m m−i
v
2 XX
 ‹
α k +1 t m −i
D∗ Ψp (x ) = 2 2 am i (−1)m−i − j 2k j n m−i − j f j (x ), (20)
πγm i =0 j =0 j

where
n+1
dn
tj
Z
2k
€ Š 1
f j (x ) = D∗α x j χ[ nk , n+1 ] = dtn
d t χ[ nk , n+1
k ]
k 2 2 Γ (n − α) n (x − t )α−n+1 2 2
2k

x dn
tj
Z
dtn
+ d t χ[ n+1
k ,1]
, j = 0, . . . , m − i . (21)
n+1 (x − t )α−n +1 2
2k

Here, by knowing that f j (x ) is zero outside the interval [ 2nk , n2+1

k ], it is consequence that Chebyshev wavelets expan-

sion of this function have only components of basis Chebyshev wavelets Ψ(x ) that are non-zero in this interval that
yield:
(n+1)M
X
f j (x ) = b j q Ψq (x ), j = 0, . . . , m − i , (22)
q =n M

where b j q = (f j (x ), Ψq (x )). Now by substitute (22) into (20), we obtain:

(nX
+1)M X
!
m m−i
v
2
 ‹
α k +1 t
X m − i m−i − j k j m−i − j
D∗ Ψp (x ) = 2 2 b j q a mi (−1) 2 n Ψq (x ), (23)
πγm q =n M
j
i =0 j =0

and this leads to desired results.

86
 Fakhrodin Mohammadi / Int. J. Adv. Appl. Math. and Mech. 2(1) (2014) 83 - 91

Lemma 3.2.
For the integer value α = 1 Chebyshev wavelet operational matrix of derivative can be expressed by

d Ψ(x )
= D Ψ(x ) (24)
dt

where D is the 2k M operational matrix of derivative defined as follow

F 0 ··· 0
 
 0 F ··· 0 
D =
 .. .. . . . .. 
 (25)
. . .
0 0 0 F

In which F is M × M matrix and its (r, s )-th element is defined as follow

 qγ
2k +2 m r −1
γs −1 r = 2, ..., M , s = 1, ..., r − 1 a nd (r + s ) o d d
Fr, s = (26)
0 o.w

Proof. It is an immediate consequence of Theorem 3.1.

Corollary 3.1.
By using Eq. (24) the operational matrix for n-th derivative can be derived as

d n Ψ(x )
= D n Ψ(x ), (27)
d xn
where D n is the n-th power of matrix D .

3.1. Error analysis

In this section we give the convergence properties and error bound for the Chebyshev wavelets expansion.

Theorem 3.2.
Any square integrable function f (x ) in [0, 1) with bounded second derivative | f 00 (x )| ≤ B can be expressed in terms of
Chebyshev wavelets, and the series converges uniformly to f (x ), that is
∞ X
X ∞
f (x ) = cn m ψn m (x ), (28)
n =1 m=0

where cn m = (f (x ), ψn m (x )) and (., .) denotes the inner product on L w

2
n
[0, 1].

Proof. Please see [18].

Theorem 3.3.
Suppose f (x ) be a continuous function defined on [0, 1), with bounded second derivative | f 00 (x )| ≤ B , and
P2k −1 PM −1
n =0 m=0 c n m ψn m (x ) be its approximated value using Chebyshev wavelets, then we have the following accuracy
estimation
p ‚ ∞ ∞
Π12
B π X 1 X 1
σM ,k = (29)
8 n =2k
n 5 m=M (m − 1)4

where
  2  12
1 k
Z 2X −1 M
X −1
σk,M =   f (x ) − cn m ψn m (x ) wn (x )d x 
0 n =0 m=0

Proof. Please see [18].

87
 Numerical solution of Bagley-Torvik equation using Chebyshev wavelet operational matrix of fractional derivative

4. Application and results

In this section, in order to show the high importance of operational matrix of derivative, we apply it to solve fractional
order Bagley-Torvikboundary value problems. These problems are considered because closed form solutions are
available for them, or they have also been solved using other numerical schemes. This allows one to compare the
results obtained using this scheme with the analytical solution or the solutions obtained using other schemes. The
general Bagley-Torvik boundary value problems of order 2 have the form [19]
3
€ Š
A 0 D 2 + A 1 D∗2 + A 2 y (t ) = f (t ) , t ∈ [0, T ], (30)

subject to initial-boundary conditions

0

γ0 y (0) + γ1 y (0) = α0 ,
0 (31)
γ3 y (T ) + γ4 y (T ) = α1 ,

where A 0 , A 1 , A 2 , α0 , γ0 , γ1 , γ3 , γ4 , α0 and α1 are constants with A 0 6= 0, and y ∈ L 1 [0, T ]. The existence and unique-
ness of the exact solution of the solution for these problems are discussed in [3]. Here we introduce a new method
based on Chebyshev wavelets expansion and their operational matrices of derivatives. To solving the Bagley-Torvik
boundary value problems of the form (30) subject to the conditions (31) we approximate the y (t ) and f (t ) by the
Chebyshev wavelets as

y (t ) ' C T Ψ(t ), f (t ) ' F T Ψ(t ) (32)

where C and F are coefficient vector of Chebyshev wavelets expansion for functions y (t ) and f (t ). By using theo-
rems 3.1 we have

D y (x ) ' C T D 2 Ψ(x ),
 2
3 3 (33)
D 2 y (x ) ' C T D 2 Ψ(x ),

Substituting Eqs. (37)in (30) the residual R (x ) can be derived as

3
€ Š
R (x ) ' A 0 C T D 2 + A 1 C T D 2 + A 2 C T Ψ(x ). (34)

By using typical tau method [20] we generate 2k M − 2 linear equations by applying

R2k M (x ), Ψ j (x ) = 0, j = 0, 1, ..., 2k M − 2


(35)

Also, by substituting boundary conditions we get

γ0 C T Ψ(0) + γ1 C T D Ψ(0) = α0 ,
§
(36)
γ3 C T Ψ(T ) + γ4 C T D Ψ(T ) = α1 ,

Eqs. (35) and (36) together generate 2k M set of linear equations. These linear equations can be solved for unknown
coefficients of the vector C . By substituting the derived vector C in Eq. (37) approximation solution y (x ) can be
obtained.

5. Numerical experiments
In this section we will consider the three fractional order Bagley-Torvik boundary value problems. We used the
method described in the Section 4 for solving these problems. The algorithms are performed by Maple 12 with 30
digits precision.

Example 5.1.
Consider the following boundary value problem in the case of the inhomogeneous Bagley-Torvik equation [19]
 3
Æ
D 2 y (x ) + D 2 y (x ) + y (x ) = t 2 + 4 πt + 2
(37)
y (0) = 0, y (5) = 25.

Where the exact solution of this problem is y (x ) = x 2 . We solve this fractional boundary value problem by applying
the method described in Section 4 using Chebyshev wavelets expansion and its operational matrices of derivatives
with M=4, k=1. Using (35) we obtain four linear equations. And by applying boundary condition we have two linear
equations. By solving this linear system we get the unknown vector C . By substituting this vectorin Eq. (37) we have
the exact solution.

88
 Fakhrodin Mohammadi / Int. J. Adv. Appl. Math. and Mech. 2(1) (2014) 83 - 91

Table 1. The absolute errors for different values of M and k

x (M , k ) = (6, 1) (M , k ) = (6, 2) (M , k ) = (12, 1) (M , k ) = (12, 2)
4 1.1 × 10−3 2.4 × 10−4 3.9 × 10−6 6.1 × 10−7
8 4.4 × 10−4 3.5 × 10−4 7.6 × 10−7 8.8 × 10−9
12 2.6 × 10−4 5.4 × 10−5 8.7 × 10−8 5.4 × 10−9
16 2.0 × 10−4 6.5 × 10−6 l9.5 × 10−8 7.4 × 10−10

Example 5.2.
Consider the boundary value problem
 3
 D 2 y (x ) + y (x ) = x 5 − x 4 + 7128
p x 3.5 − 64
π
p x 2.5
5 π
(38)
y (0) = 0 , y (1) = 0.


where the exact solution of this problem is y (x ) = x 4 (1 − x ). We solve this fractional boundary value problem by
applying the method described in Section 4 by using Chebyshev wavelets expansion M = 6, k = 2. Similar to Example
1 by solving linear system derived for this problem we get the exact solution.

Example 5.3.
Consider the boundary value problem
 3
 D 2 y (x ) + 12 D 2 y (x ) + 21 y (x ) = 8ν(x ) − 8ν(x − 1)
(39)
 y (0) = 0 , y (20) = −1.48433,

In which ν(x ) is the Heaviside function. This example was solved theoretically by Podlubny [3] and the compact
form of the solution is given by
Z x
y (x ) = 8G (x − t )(ν(t ) − ν(t − 1))d t (40)
0

where G is the fractional Green function defined as

€ 1 Šj
X∞
(−1)k X
∞ ( j + k )! −t 2
G (t ) = (41)
2k k ! j =0 j !2 j Γ 12 j + 12 k + 32 k + 2


k =0

In Eq. (41) is the Gamma function. Here we solve this problem using Chebyshev wavelets method described in Sec-
tion 4 with (M , k ) = (6, 1), (M , k ) = (6, 2), (M , k ) = (12, 1) and (M , k ) = (12, 2). Fig. 1 shows the approximate solution and
theoretical solution derived by Podlubny[3] for (M , k ) = (12, 2) in the interval [0, 20]. The absolute errors for different
values in the interval are shown in the Table 1. From Table 1, we see that we can achieve a good approximation with
the exact solution by usinga few terms of Chebyshev wavelets.

6. Conclusion
In this article a general formulation for deriving the Chebyshev wavelet operational matrix of derivatives has been
derived. Then a numerical method based on Chebyshev wavelets expansion, its operational matrix of fractional
order and tau method is introduced for approximate the solution of Bagley-Torvik fractional boundary value prob-
lems. Moreover, the convergence and error analysis for the proposed method is considered.

89
 Numerical solution of Bagley-Torvik equation using Chebyshev wavelet operational matrix of fractional derivative

Fig. 1. Approximate solution and exact solution for Example 3 in the interval [0, 20]

References

[1] K.S. Miller, B. Ross, An Introduction to The Fractional Calculus and Fractional Differential Equations, Wiley,
New York, 1993.
[2] K.B. Oldham, J. Spanier, The Fractional Calculus, Academic Press, New York, 1974.
[3] I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999.
[4] A. Saadatmandi, M. Dehghan, A new operational matrix for solving fractional-orderdifferential equa-
tions,Comput. Math. Appl, 59 (2010) 1326-1336.
[5] Y. Li, W. Zhao, Haar wavelet operational matrix of fractional order integration and itsapplications in solving the
fractional order differential equations,Appl Math Comput, 216 (2010) 2276-2285.
[6] P. J. Torvik and R. L. Bagley, On the appearance of the fractional derivative in the behavior of real materials,
Journal of Applied Mechanics 51.2 (1984): 294-298.
[7] R. L. Bagley and P. J. Torvik, Fractional calculus - A different approach to the analysis of viscoelastically damped
structures, AIAA Journal 21.5 (1983): 741-748.
[8] Y. Hu, Y. Luo, and Z. Lu, Analytical solution of the linear fractional differential equation by Adomian decompo-
sition method, Journal of Computational and Applied Mathematics 215.1 (2008): 220-29.
[9] A. M. A. El-Sayed, I. L. El-Kalla, and E. A. A. Ziada, Analytical and numerical solutions of multi-term nonlinear
fractional orders differential equations, Applied Numerical Mathematics 60.8 (2010): 788-97.
[10] A. Ghorbani and A. Alavi, Application of He’s variational iteration method to solve semidifferential equations of
nth order, Mathematical Problems in Engineering 2008 (2008): 1-10.
[11] Y. Cenesiz, Y. Keskin, and A. Kurnaz, The solution of the Bagley-Torvik equation with the generalized Taylor
collocation method, Journal of the Franklin Institute 347.2 (2010): 452-466.
[12] K. Diethelm, N.J. Ford, Numerical solution of the BagleyâĂŞTorvik equation, BIT 42 (2002) 490-507.
[13] Podlubny, Igor, Tomas Skovranek, and Vinagre Jara, Matrix approach to discretization of fractional derivatives
and to solution of fractional differential equations and their systems, Emerging Technologies & Factory Au-
tomation, 2009. ETFA 2009. IEEE Conference on. IEEE, 2009.
[14] A. Boggess, F. J. Narcowich, A first course in wavelets with Fourier analysis , John Wiley & Sons, 2001.
[15] Mallat S. A wavelet tour of signal processing. 2nd ed. Academic Press, 1999.
[16] M. Razzaghi, S. Yousefi, Legendre wavelets operational matrix of integration, Int. J. Syst. Sci. 32 (4) (2001) 495-
502.
[17] M. M. Khader, On the numerical solutions for the fractional diffusion equation, Commun Nonlinear Sci Numer
Simul., 16 (2011) 2535-2542.
[18] S. Sohrabi, Comparison Chebyshev wavelets method with BPFs method for solving AbelâĂŹs integral equation,
Ain Shams Engineering Journal (2011) 2, 249-254.
[19] Q. M. Al-Mdallal, M. I. Syam, M.N. Anwar, A collocation-shooting method for solving fractional boundaryvalue
problems, Commun Nonlinear SciNumerSimulat 15 (2010) 3814-3822.

90
 Fakhrodin Mohammadi / Int. J. Adv. Appl. Math. and Mech. 2(1) (2014) 83 - 91

[20] C. Canuto, M. Hussaini, A. Quarteroni, T. Zang, Spectral Methods in Fluid Dynamics, Springer, Berlin, 1988.

91