On the approximate solutions for system of fractional integro-differential equations using Chebyshev pseudo-spectral method

Accepted Manuscript On the approximate solutions for system of fractional integro-differential equa‐ tions using Chebyshev pseudo-spectral method M.M. Khader, N.H. Sweilam PII: DOI: Reference:

S0307-904X(13)00395-8 http://dx.doi.org/10.1016/j.apm.2013.06.010 APM 9564

To appear in:

Appl. Math. Modelling

Received Date: Revised Date: Accepted Date:

4 September 2011 4 November 2012 14 June 2013

Please cite this article as: M.M. Khader, N.H. Sweilam, On the approximate solutions for system of fractional integrodifferential equations using Chebyshev pseudo-spectral method, Appl. Math. Modelling (2013), doi: http:// dx.doi.org/10.1016/j.apm.2013.06.010

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

On the approximate solutions for system of fractional integro-differential equations using Chebyshev pseudo-spectral method M. M. Khader1,2 and N. H. Sweilam3 1

Department of Mathematics and Statistics, College of Science, Al-Imam Mohammed Ibn Saud Islamic University (IMSIU), Riyadh: 11566, Saudi Arabia

2

Permanent address: Department of Mathematics, Faculty of Science, Benha University, Benha, Egypt 3

Department of Mathematics, Faculty of Science, Cairo University, Giza, Egypt [email protected], and

n− [email protected]

Abstract In this paper, we implement Chebyshev pseudo-spectral method for solving numerically system of linear and non-linear fractional integro-differential equations of Volterra type. The proposed technique is based on the new derived formula of the Caputo fractional derivative. The suggested method reduces this type of systems to the solution of system of linear and non-linear algebraic equations. We give the convergence analysis and an upper bound of the error for the derived formula. To demonstrate the validity and applicability of the suggested method, some test examples are given. Also, we present a comparison with the previous work using the homotopy perturbation method.

Keywords: Chebyshev pseudo-spectral method; Systems of fractional integro-differential equations of Volterra type; Caputo fractional derivative; Convergence analysis.

1. Introduction The fractional integro-differential equations (FIDEs) play a pivotal role in the modeling of number of physical phenomenon from science and engineering. A particularly rich source is electric-circuit analysis and the activity of interacting inhibitory and excitatory neurons. Also, the control of systems governed by FIDEs with application to floating structures and viscoelastic material dynamics. In computational neuroscience, the Wilson-Cowan model describes the dynamics of interactions between populations of very simple excitatory and inhibitory model neurons. It was developed by H. R. Wilson and J. D. Cowan and extensions of the model have been widely used in modeling neuronal populations. The model is important historically because it uses phase plane methods and numerical solutions to describe the responses of neuronal populations to stimuli. Because the model neurons are simple, only elementary limit cycle behavior, i.e. neural oscillations, and stimulus-dependent evoked responses are predicted. The key findings include the existence of multiple stable states, and hysteresis, in the population response 1

([15], [22]). Consequently, considerable attention has been given to the solutions of FIDEs of physical interest. Most FIDEs do not have exact solutions, so approximate and numerical techniques ([1], [2], [16], [24], [25], [27]), must be used. Recently, several numerical methods to solve the FIDEs have been given such as variational iteration method ([7], [21]), homotopy perturbation method ([3], [4], [8]), Adomian’s decomposition method [5], homotopy analysis method [6] and collocation method ([9]-[14], [18]). In the following we present some basic definitions for fractional derivatives which will be required in the present paper. Definition 1. The Caputo fractional derivative operator Dα of order α is defined in the following form Z x 1 f (m) (t) α D f (x) = dt, α > 0, x > 0, Γ(m − α) 0 (x − t)α−m+1 where m − 1 < α ≤ m, m ∈ N. Similar to integer-order differentiation, Caputo fractional derivative operator is a linear operation Dα (λ f (x) + µ g(x)) = λ Dα f (x) + µ Dα g(x), where λ and µ are constants. For the Caputo’s derivative we have [17] Dα C = 0, ( D α xn =

C is a constant,

(1)

0,

for n ∈ N0 and n < dαe;

Γ(n+1) xn−α , Γ(n+1−α)

for n ∈ N0 and n ≥ dαe.

(2)

We use the ceiling function dαe to denote the smallest integer greater than or equal to α and N0 = {0, 1, 2, ...}. Recall that for α ∈ N, the Caputo differential operator coincides with the usual differential operator of integer order. For more details on fractional derivatives definitions and its properties see ([17], [19]). The main aim of this paper is to introduce a numerical scheme for the solution of the system of fractional integro-differential equations of the general form ([15], [22])   Z x α D ui (x) = fi x; u1 , u2 , ..., un , ki (t, u1 , u2 , ..., un ) dt , i = 1, 2, ..., n,

(3)

0

with the following initial conditions (at 0 < α ≤ 1) ui (0) = i ,

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

(4)

where i , i = 1, 2, ..., n, are specified constants, α is a parameter describing the order of the space fractional derivative, and ui (x), the solution vector, is assumed to be a causal function of space, i.e., vanishing for x < 0. The fractional derivative is considered in the Caputo sense. The general response expression contains a parameter describing the order of the fractional derivative that can be varied to obtain various responses. When α = 1, the system of equations (3) reduces to the classical system of integro-differential equations of the form [26]   Z x dui = fi x; u1 , u2 , ..., un , ki (t, u1 , u2 , ..., un ) dt , i = 1, 2, ..., n. dx 0

(5)

Khader in [9] introduced a new approximate formula of the fractional derivative of the function u(x), where this function expanded by the Chebyshev polynomials. This formula is used to introduce an efficient numerical method for solving the fractional diffusion equation. The main aim of this work is to extend and apply the introduced formula in ([9], [23]) to solve numerically the system of FIDEs of Volterra type (3). Our technique uses the properties of Chebyshev polynomials to reduce the proposed models to linear/non-linear system of algebraic equations thus greatly simplifying the problem, and use a suitable method to solve the resulting system. Also, we present an estimation to the upper bound of the error of the introduced approximated formula. The organization of this paper is as follows: In the next section, the approximation of the fractional derivative Dα u(x) and the estimation of the upper bound of the error are obtained. In section 3, the application of Chebyshev collocation method to solve (3) is presented. As a result, a system of linear/non-linear of algebraic equations is formed and the solution of the considered problem is introduced. In section 4, some numerical examples are given to clarify the advantages of the method and a comparison with homotopy perturbation method is introduced. Also a conclusion is given in section 5. 2. Derivation an approximate formula for fractional derivatives using Chebyshev series expansion [9] Chebyshev polynomials are well known family of orthogonal polynomials on the interval [−1, 1] that have many applications [20]. In which have good widely used because of their good properties in the approximation of functions. However, with our best knowledge, very little work was done to adapt these polynomials to the solution of FIDEs. The well known Chebyshev polynomials [20] are defined on the interval [−1, 1] and can be determined with the aid of the following recurrence formula

3

Tn+1 (z) = 2 z Tn (z) − Tn−1 (z),

T0 (z) = 1, T1 (z) = z,

n = 1, 2, ... .

The analytic form of the Chebyshev polynomials Tk (z) of degree k is given by k

Tk (z) =

[2] X

(−1)i 2k−2 i−1

i=0

k (k − i − 1)! k−2 i z . (i)! (k − 2 i)!

Where [k/2] denotes the integer part of k/2. The orthogonality condition is   for i = j = 0;  Z 1  π, Ti (z) Tj (z) π √ dz = , for i = j = 6 0; 2  1 − z2 −1   0, for i 6= j.

(6)

(7)

In order to use these polynomials on the interval [0, 1] we define the so called shifted Chebyshev polynomials by introducing the change of variable z = 2x − 1.

√ The shifted Chebyshev polynomials are defined as Tk∗ (x) = Tk (2x − 1) = T2n ( x) and the analytic form is given by Tk∗ (x)

k X 22i (k + i − 1)! i =k (−1)k−i x, (2i)! (k − i)! i=0

k = 0, 1, ...

(8)

The function u(x), which belongs to the space of square integrable in [0, 1], may be expressed in terms of shifted Chebyshev polynomials as u(x) =

∞ X

ci Ti∗ (x),

(9)

i=0

where the coefficients ci , i = 0, 1, ... are given by Z 1 2 u(x) Ti∗ (x) √ ci = dx, h0 = 2, π hi 0 x − x2

hi = 1,

i = 1, 2, ....

(10)

In practice, only the first (m + 1)-terms of shifted Chebyshev polynomials are considered. Then we have um (x) =

m X

ci Ti∗ (x).

(11)

i=0

Theorem 1. (Chebyshev truncation theorem) The error in approximating u(x) by the sum of its first m terms is bounded by the sum of the absolute values of all the neglected coefficients. If um (x) =

m X k=0

4

ck Tk (x),

(12)

then ET (m) ≡ |u(x) − um (x)| ≤

∞ X

|ck |,

(13)

k=m+1

for all u(x), all m, and all x ∈ [−1, 1]. Proof. The Chebyshev polynomials are bounded by one, that is, |Tk (x)| ≤ 1 for all x ∈ [−1, 1] and for all k. This implies that the k-th term is bounded by |ck |. Subtracting the truncated series from the infinite series, bounding each term in the difference, and summing the bounds gives the theorem. The main approximate formula of the fractional derivative of u(x) is given in the following theorem. Theorem 2. Let u(x) be approximated by Chebyshev polynomials as (11) and suppose α > 0, then m i X X

α

D (um (x)) =

(α)

ci bi, k xk−α ,

(14)

i=dαe k=dαe (α)

where bi, k is given by (α)

bi, k = (−1)i−k

22k i (i + k − 1)!Γ(k + 1) . (i − k)! (2 k)! Γ(k + 1 − α)

(15)

Proof. Since the Caputo’s fractional differentiation is a linear operation we have Dα (um (x)) =

m X

ci Dα (Ti∗ (x)).

(16)

i=0

Employing Eqs.(1)-(2) in formula (8), we get Dα Ti∗ (x) = 0,

i = 0, 1, ..., dαe − 1,

α > 0.

(17)

Also, for i = dαe, dαe + 1, ..., m, and by using Eqs.(1)-(2) in formula (8), we get D

α

Ti∗ (x)

=i

i X

(−1)i−k

k=dαe

=i

i X k=dαe

22k (i + k − 1)! α k D x (i − k)! (2 k)! (18) 2k

(−1)i−k

2 (i + k − 1)!Γ(k + 1) k−α x . (i − k)! (2 k)! Γ(k + 1 − α)

A combination of Eqs.(15)-(18) leads to the desired result (14). 5

Theorem 3. The Caputo fractional derivative of order α for the shifted Chebyshev polynomials can be expressed in terms of the shifted Chebyshev polynomials themselves in the following form D

α

(Ti∗ (x))

k−dαe i X X

=

Θi,j,k Tj∗ (x),

(19)

k=dαe j=0

where Θi,j,k =

(−1)i−k 2i(i + k − 1)!Γ(k − α + 12 ) , hj Γ(k + 12 ) (i − k)! Γ(k − α − j + 1) Γ(k + j − α + 1)

j = 0, 1, ... .

Proof. Using the properties of the shifted Chebyshev polynomials [20], then xk−α in (18) can be expanded in the following form [9] k−dαe

x

k−α

=

X

ckj Tj∗ (x),

(20)

j=0

here ckj can be obtained using (10) where u(x) = xk−α , then we can claim the following Z 1 k−α ∗ x Tj (x) 2 √ ckj = dx, h0 = 2, hj = 1, j = 1, 2, ... . hj π 0 x − x2 Z 1 1 xk−α T0∗ (x) 1 Γ(k − α + 1/2) √ At j = 0 we have, ck0 = dx = √ , 2 π 0 π Γ(k − α + 1) x−x also, for any j, using the formula (10), we can claim j j X (j + r − 1)!22r+1 Γ(k + r − α + 1/2) , =√ (−1)j−r (j − r)!(2r)!Γ(k + r − α + 1) π r=0

ckj

j = 1, 2, ...,

employing Eqs.(18) and (20) gives D

α

(Ti∗ (x))

=

k−dαe i X X

Θi,j,k Tj∗ (x),

i = dαe, dαe + 1, ... ,

k=dαe j=0

where

Θi,j,k =

 (−1)i−k (i+k−1)! 22k k! Γ(k−α+ 12 )  √ i ,  (i−k)! (2k)! π (Γ(k+1−α))2    

(−1)i−k ij (i+k−1)! 22k+1 k! √ πΓ(k+1−α)(i−k)!(2k)!

×

Pj

r=0

j = 0; (−1)j−r (j+r−1)! 22r Γ(k+r−α+ 12 ) , (j−r)! (2r)! Γ(k+r−α+1)

j = 1, 2, ... .

After some lengthly manipulation Θi,j,k can put in the following form Θi,j,k

(−1)i−k 2i(i + k − 1)!Γ(k − α + 12 ) = , hj Γ(k + 12 ) (i − k)! Γ(k − α − j + 1) Γ(k + j − α + 1)

and this completes the proof of the theorem. 6

j = 0, 1, ... ,

(21)

Theorem 4. The error |ET (m)| = |Dα u(x) − Dα um (x)| in approximating Dα u(x) by Dα um (x) is bounded by k−dαe ∞ i X  X  X |ET (m)| ≤ ci Θi,j,k . i=m+1

(22)

k=dαe j=0

Proof. A combination of Eqs.(9), (11) and (19) leads to k−dαe ∞ i X  X  X α α ∗ |ET (m)| = D u(x) − D um (x) = ci Θi,j,k Tj (x) , i=m+1

k=dαe j=0

but |Tj∗ (x)| ≤ 1, so, we can obtain k−dαe ∞ i X  X X |ET (m)| ≤ ci Θi,j,k , i=m+1

k=dαe j=0

and subtracting the truncated series from the infinite series, bounding each term in the difference, and summing the bounds completes the proof of the theorem. 3. Procedure of solution for system of fractional integro-differential equations In this section, we introduce a discretization formula of (3) using the Chebyshev collocation method. To achieve this aim, we approximate ui (x) as uim (x) =

m X

cik Tk∗ (x),

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

(23)

k=0

From Eqs.(3), (23) and Theorem 2 we have m k X X

(α) ci k bk,r

xr−α = fi

x;

m X

ci k Tk∗ (x),

ki 0

k=0

k=dαe r=dαe

x

Z

t,

m X

!

!

ci k Tk∗ (t) dt .

(24)

k=0

We now collocate Eq.(24) at (m + 1 − dαe) points xp as m k X X

(α)

ci k bk,r xr−α = fi p

k=dαe r=dαe

xp ;

m X

ci k Tk∗ (xp ),

xp

Z

ki 0

k=0

t,

m X

!

!

ci k Tk∗ (t) dt .

(25)

k=0

∗ For suitable collocation points we use roots of shifted Chebyshev polynomial Tm+1−dαe (x). In

order to use the Gaussian integration formula for Eq.(25), we transfer the t-interval [0, xp ] into τ -interval [−1, 1] by means of the transformation τ = may then be restated as 7

xp t 2

− 1, Eq.(25), for p = 0, 1, ..., m − dαe,

m k X X

(α)

ci k bk,r xr−α = fi p

xp ;

m X

ci k Tk∗ (xp ),

k=0

k=dαe r=dαe

xp 2

!

m X

1

Z

xp xp (τ + 1), ci k Tk∗ ( (τ + 1)) dτ 2 2 k=0

ki −1

! ,

(26) by using the Gaussian integration formula, for p = 0, 1, ..., m − dαe, we get m k X X

 (α)

ci k bk,r xr−α ' fi xp ; p

m X k=0

k=dαe r=dαe

m xp X ∗ ci k Tk (xp ), ωq ki 2 q=0

xp (τq + 1), 2

m X k=0

!

ci k Tk∗ (

xp (τq + 1))  , 2 (27)

where τq are m + 1 zeros of Chebyshev polynomial Tm+1 (t) and wq are the corresponding weights given in [18]. The idea behind the above approximation is the exactness of the Gaussian integration formula for polynomials of degree not exceeding 2m + 1. Also, by substituting Eq.(11) in the given initial conditions (4) we can find ndαe equations as follows

m X (−1)k ci,k = i ,

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

(28)

k=0

Equations (27), together with ndαe equations of initial conditions (28), give n(m + 1) algebraic equations which can be solved for the unknowns ci k , k = 0, 1, ..., m, using a suitable method. Consequently ui (x), i = 1, 2, ..., n given in Eq.(23) can be calculated. In our computational we used the conjugate gradient method to solve the resulting linear system of algebraic equations, and the Newton’s iteration method to solve the resulting non-linear system of algebraic equations. 4. Numerical implementation with comparison and discussion In this section, in order to illustrate our theoretical results obtained in preceding sections and to compare these in the context of various applications of physics, we now present some numerical results. In the calculation process, we implement the proposed method to solve two examples of linear and non-linear systems of FIDEs. These examples can be used as a basis to illustrate the applicability of the proposed method. Example 1: Consider the following linear system of FIDEs (3) for i = 1, 2 and 0 < α ≤ 1 Z x α D u1 (x) = `(x) − u2 (x) − [u1 (t) + u2 (t)]dt,

(29)

0

Z

α

D u2 (x) = γ(x) + u1 (x) −

[u1 (t) − u2 (t)]dt, 0

8

x

(30)

where `(x) = 1 + x + x2 , γ(x) = −1 − x, with the initial conditions u2 (0) = −1.

u1 (0) = 1,

(31)

We implement the suggested method with m = 5, and approximate solution as u1 (x) '

5 X

c1 k Tk∗ (x),

u2 (x) '

k=0

5 X

c2 k Tk∗ (x).

(32)

k=0

Using Eq.(27) we have m k X X

(α) c1 k bk,r xr−α p

' `(xp ) −

5 X

k=dαe r=dαe

k=0

m k X X

5 X

(α)

c2 k bk,r xr−α ' γ(xp ) + p

k=0

k=dαe r=dαe

c2 k Tk∗ (xp )

" 5 # 5   X xp X ∗ xp − wq (c1 k + c2 k )Tk (τq + 1) , (33) 2 2 q=0

k=0

" 5 # 5 x  X X x p p c1 k Tk∗ (xp ) − wq (c1 k − c2 k )Tk∗ (τq + 1) , (34) 2 2 q=0

k=0

with p = 0, 1, 2, 3, 4 where xp are roots of the shifted Chebyshev polynomial T5∗ (x) and their values are x0 =

1 2

q q √ √ 1 1 − 2(5 − 5), x1 = 2 − 8 2(5 + 5), x2 = 12 , q q √ √ x3 = 12 + 18 2(5 + 5), x4 = 12 + 18 2(5 − 5). 1 8

Also τq are the roots of Chebyshev polynomials T6 (x) and their values are τ0 = 0.25882, τ1 = −0.25882, τ2 = 0.965926, τ3 = −0.96593, τ4 = 0.70711, τ5 = −0.70711. And wq are the corresponding weights and their values are w0 = 0.086716180, w3 = 0.454421768,

w1 = 0.287831394, w4 = 0.398241540,

w2 = 0.398241540, w5 = 0.287831394.

By using Eqs.(28) and (31) we get c10 − c11 + c12 − c13 + c14 − c15 = 1,

(35)

c20 − c21 + c22 − c23 + c24 − c25 = −1.

(36)

Now the solution of the linear system of equations (33)-(36) is given by c10 = 2.2534,

c11 = 1.3504,

c12 = 0.1052,

c13 = 0.0088,

c14 = 0.0005,

c15 = 0.0,

c20 = −1.2534, c21 = −0.3504, c22 = −0.1052, c23 = −0.0088, c24 = −0.0005, c25 = 0.0. 9

Thus using equation (32) we get u1 (x) ∼ = 1.00001 + 1.99951 x + 0.508008 x2 + 0.140811 x3 + 0.0699454 x4 , u2 (x) ∼ = −0.9999 − 0.0016 x − 0.4992 x2 − 0.1536 x3 − 0.064 x4 , which are the approximate solutions of the problem (29)-(30). The exact solution at α = 1 is given by u1 (x) = x + ex ,

u2 (x) = x − ex .

The obtained numerical results by means of the proposed method are shown in table 1 and figures 1 and 2. In table 1, the absolute error between the exact solution (u1ex , u2ex ) and the approximate solution (u1approx. , u2approx. ), at m = 3 (in columns 2, 3) and m = 5 (in columns 4, 5) respectively, are given. Figure l shows the evolution results for the system of fractional integro-differential equations (29)-(30) at m = 7 when α = 1. And figure 2 shows the behavior of the obtained approximate solution for the proposed system (29)-(30) at m = 7 with different values of α (α = 0.85 and 0.95). From these figures we can conclude that our approximate solution is in an excellent agreement with the exact values and with high accuracy if we compared with the approximate solution obtained in [22]. From Figures l and 2, it is easy to conclude that the solution continuously depends on the space-fractional derivative. This note confirms the physical meaning of the behavior of the solution for the proposed real problem. Table 1:

The absolute error between the exact solution and the approximate solution at m = 3 (in columns 2, 3) and m = 5 (in columns 4, 5).

x

|u1ex − u1approx. | |u2ex − u2approx. | |u1ex − u1approx. | |u2ex − u2approx. |

0.1

0.11e-02

0.60e-03

0.1059e-04

0.1584e-04

0.2

0.12e-02

0.50e-03

0.6446e-04

0.0632e-04

0.3

0.70e-03

0.20e-03

0.9289e-04

0.1294e-04

0.4

0.20e-03

0.10e-03

0.7228e-04

0.0496e-04

0.5

0.10e-03

0.20e-03

0.1503e-04

0.3809e-04

0.6

0.10e-03

0.20e-03

0.4634e-04

0.6712e-04

0.7

0.30e-03

0.70e-03

0.7659e-04

0.7251e-04

0.8

0.70e-03

0.13e-02

0.5678e-04

0.4584e-04

0.9

0.80e-03

0.13e-02

0.0514e-04

0.0312e-04

10

Figure 1. The behavior of the exact solution and the approximate solution at m = 7, α = 1.

Figure 2. The behavior of the approximate solution at m = 7, α = 0.85 and 0.95. 11

Example 2: Consider the following non-linear system of fractional differential equations for 0 < α ≤ 1 Dα u(x) = u(x)(a − b v(x)),

(37)

Dα v(x) = −v(x)(c − d u(x)),

(38)

with the initial conditions u(0) = u0 ,

v(0) = v 0 .

(39)

We implement the suggested method with m = 3, and approximate the solution as u(x) '

3 X

c1 k Tk∗ (x),

v(x) '

k=0

3 X

c2 k Tk∗ (x).

(40)

k=0

Using Eq.(27) we have 3 k X X

(α) c1 k bk,r

3 X

xr−α ' p

c1 k Tk∗ (xp )

a−b

k=0

k=dαe r=dαe 3 k X X

!

(α) c2 k bk,r

xr−α '− p

! c2 k Tk∗ (xp ) ,

(41)

k=0

3 X

! c2 k Tk∗ (xp )

c−d

k=0

k=dαe r=dαe

3 X

3 X

! c1 k Tk∗ (xp ) ,

(42)

k=0

with p = 0, 1, 2 where xp are roots of the shifted Chebyshev polynomial T3∗ (x) and their values are x0 =

3 , 50

x1 =

3 (2 100

+

√ 3),

x2 =

3 (2 100

−

√

3).

By using Eqs.(28) and (39) we get c10 − c11 + c12 − c13 = u0 ,

(43)

c20 − c21 + c22 − c23 = v 0 .

(44)

Now the solution of the nonlinear system of equations (41)-(44) at (u0 = 16, v 0 = 10, a = b = 1, c = 0.1, d = 1 and α = 1) is give by c10 = 5.639831,

c11 = −2.582608,

c20 = 19.775424,

c21 = 1.728723,

c12 = 0.126214, c22 = 0.212328,

c13 = −7.651282, c23 = 8.259036.

Thus using equation (40) we get u(x) ∼ = 15.9999 − 1199.15 x + 25574.4 x2 − 141690 x3 , v(x) ∼ = 9.99999 + 1253.51 x − 27412.2 x2 + 152945 x3 .

12

Figure 3. The behavior of the approximate solution at m = 3, α = 1.

Figure 4. The behavior of the approximate solution at m = 3, α = 0.9. Figures 3-6 show the evolution results for the system of differential equations (37)-(38) when m = 3 at α = 1, α = 0.9, α = 0.8 and α = 0.7, respectively. 13

Figure 5. The behavior of the approximate solution at m = 3, α = 0.8.

Figure 6. The behavior of the approximate solution at m = 3, α = 0.7. 14

From these figures we can conclude that our approximate solution is in an excellent agreement with the exact values and with high accuracy if we compared with the approximate solution obtained in [22]. From Figures 3-6, it is easy to conclude that the solution continuously depends on the space-fractional derivative. This note confirms the physical meaning of the behavior of the solution for the proposed real problem. 5. Conclusion In this article, the Chebyshev spectral method is implemented for solving the linear and nonlinear systems of fractional integro-differential equations. The fractional derivative is considered in the Caputo sense. The properties of the Chebyshev polynomials are used to reduce the proposed problem to the solution of a system of linear or non-linear algebraic equations which is solved by using conjugate gradient method or Newton iteration method, respectively. Special attention is given to study the convergence analysis and estimate an upper bound of the error of the derived formula. From the obtained numerical results we can see that the obtained solutions using the suggested method are in excellent agreement with the exact solution and with more accuracy than the solution obtained using homotopy perturbation method and show that this approach can be solved the problem effectively. It is evident that the overall errors can be made smaller by adding new terms in the series (32). Also, from the presented figures for the introduced models it is easy to conclude that the solution continuously depends on the space-fractional derivative and this note confirms the physical meaning of the behavior of the solution for the proposed real problem. All computations in this paper are done using Matlab 7.1.

References [1] S. Das, Functional Fractional Calculus for System Identification and Controls, Springer, New York, 2008. [2] Z. Z. Ganji, D. D. Ganji and M. Esmaeilpour, Study on nonlinear Jeffery-Hamel flow by He’s semi-analytical methods and comparison with numerical results, Computers and Mathematics with Applications, 58(11-12), (2107-2116), 2009. [3] D. D. Ganji and A. Rajabi, Assessment of homotopy-perturbation and perturbation methods in heat radiation equations, International Communications in Heat and Mass Transfer, 33(3), (391400), 2006. [4] D. D. Ganji and A. Sadighi. Application of He’s homotopy-perturbation method to nonlinear coupled systems of reaction-diffusion equations, International Journal of Nonlinear Sciences and Numerical Simulation, 7(4), (411-418), 2006.

15

[5] I. Hashim, Adomian decomposition method for solving BVPs for fourth-order integro-differential equations, Journal of Computational and Applied Mathematics, 193, (658-664), 2006. [6] I. Hashim, O. Abdulaziz and S. Momani, Homotopy analysis method for fractional IVPs, Communications in Nonlinear Science and Numerical Simulation, 14, (674-684), 2009. [7] J. H. He, Approximate analytical solution for seepage flow with fractional derivatives in porous media, Computer Methods in Applied Mechanics and Engineering, 167(1-2), (57-68), 1998. [8] J. H. He, Homotopy perturbation technique, Comput. Methods Appl. Mech. Engng., 178(3-4), (257-262), 1999. [9] M. M. Khader, On the numerical solutions for the fractional diffusion equation, Communications in Nonlinear Science and Numerical Simulation, 16, (2535-2542), 2011. [10] M. M. Khader, Numerical treatment for solving fractional Riccati differential equation, Journal of the Egyptian Mathematical Society, 21, (32-37), 2013. [11] M. M. Khader, Numerical treatment for solving the perturbed fractional PDEs using hybrid techniques, Journal of Computational Physics, 250, (565-573), (2013. [12] M. M. Khader and A. S. Hendy, A numerical technique for solving fractional variational problems, Mathematical Methods in Applied Sciences, 36(10), (1281-1289), 2013. [13] M. M. Khader, T. S. EL Danaf and A. S. Hendy, A computational matrix method for solving systems of high order fractional differential equations, Applied Mathematical Modelling, 37, (40354050), 2013. [14] M. M. Khader, N. H. Sweilam and A. M. S. Mahdy, Numerical study for the fractional differential equations generated by optimization problem using Chebyshev collocation method and FDM, Applied Mathematics and Information Science, 7(5), (2013-2020), 2013. [15] S. Momani and R. Qaralleh, An efficient method for solving systems of fractional integro-differential equations, Computers and Mathematics with Applications, 52, (459-470), 2006. [16] S. Momani and M. A. Noor, Numerical methods for fourth-order fractional integro-differential equations, Applied Mathematics Computation, 182, (754-760), 2006. [17] I. Podlubny, Fractional Differential Equations, Academic Press, New York, 1999. [18] E. A. Rawashdeh, Numerical solution of fractional integro-differential equations by collocation method, Applied Mathematics Computation, 176, (1-6), 2006.

16

[19] S. Samko, A. Kilbas and O. Marichev, Fractional Integrals and Derivatives: Theory and Applications, Gordon and Breach, London, 1993. [20] M. A. Snyder, Chebyshev Methods in Numerical Approximation, Prentice-Hall, Inc. Englewood Cliffs, N. J. 1966. [21] N. H. Sweilam, M. M. Khader and R. F. Al-Bar, Numerical studies for a multi-order fractional differential equation, Physics Letters A, 371, (26-33), 2007. [22] N. H. Sweilam, M. M. Khader and R. F. Al-Bar, Homotopy perturbation method for linear and nonlinear system of fractional integro-differential equations, International Journal of Computational Mathematics and Numerical Simulation, 1(1), (73-87), 2008. [23] N. H. Sweilam and M. M. Khader, A Chebyshev pseudo-spectral method for solving fractional order integro-differential equations, ANZIAM, 51 (464-475), 2010. [24] N. H. Sweilam, M. M. Khader and A. M. S. Mahdy, Numerical studies for fractional-order Logistic differential equation with two different delays, Journal of Applied Mathematics, 2012, 2012 Article ID 764894, 14 pages. [25] N. H. Sweilam, M. M. Khader and A. M. S. Mahdy, Numerical studies for solving fractional-order Logistic equation, Int. J. of Pure and Applied Mathematics, 78(8), (1199-1210), 2012. [26] N. H. Sweilam, M. M. Khader and W. Y. Kota, Numerical and analytical study for fourth-order integro-differential equations using a pseudo-spectral method, Mathematical Problems in Engineering, 2013, 2013 Article ID 434753, 7 pages. [27] N. H. Sweilam, M. M. Khader and A. M. Nagy, Numerical solution of two-sided space fractional wave equation using finite difference method, Computation and Applied Mathematics, 235, (28322841), 2011.

17