Spectral collocation method for linear fractional integro-differential equations

Accepted Manuscript Spectral collocation method for linear fractional integro-differential equations Xiaohua Ma, Chengming Huang PII: DOI: Reference:

S0307-904X(13)00520-9 http://dx.doi.org/10.1016/j.apm.2013.08.013 APM 9649

To appear in:

Appl. Math. Modelling

Received Date: Revised Date: Accepted Date:

14 June 2012 9 January 2013 8 August 2013

Please cite this article as: X. Ma, C. Huang, Spectral collocation method for linear fractional integro-differential equations, Appl. Math. Modelling (2013), doi: http://dx.doi.org/10.1016/j.apm.2013.08.013

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.

Spectral collocation method for linear fractional integro-differential equations? Xiaohua Maa∗ , Chengming Huang a∗ a School

of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan 430074, P. R. China

Abstract In this paper, we propose and analyze a spectral Jacobi-collocation method for the numerical solution of general linear fractional integro-differential equations. The fractional derivatives are described in the Caputo sense. First, we use some function and variable transformations to change the equation into a Volterra integral equation defined on the standard interval [−1, 1]. Then the Jacobi-Gauss points are used as collocation nodes and the Jacobi-Gauss quadrature formula is used to approximate the integral equation. Later, the convergence order of the proposed method is investigated in the infinity norm. Finally, some numerical results are given to demonstrate the effectiveness of the proposed method. Keywords: fractional differential equation, Caputo derivative, Volterra integral equation, Spectral method, Convergence analysis.

1. Introduction This paper considers linear fractional integro-differential equations (FIDEs) of the form α D∗ p y(t) +

p−1

∑

ai (t)D∗αi y(t) =

i=0

f (t) +

Z t 0

K(t, s)y(s)ds,

0 ≤ t ≤ T,

(1.1)

with the initial conditions y(i) (0) = γi ,

i = 0, · · · , n − 1.

(1.2)

Here we assume that the fractional derivative orders αi satisfy 0 = α0 < α1 < · · · < α p , n − 1 < α p ≤ n, n, p ∈ N, I This

work was supported by NSF of China(No.10971077 and 91130003). author . Email addresses: [email protected] (X. H. Ma),chengming [email protected] (C. M. Huang) ∗ Corresponding

Preprint submitted to Elsevier

(1.3)

that the functions ai (i = 0, · · · , p − 1) and f are continuous functions from [0, T ] to R, and that the kernel K is a continuous function from {(t, s) : 0 ≤ s ≤ t ≤ T } to R. In this paper, D∗α y(t) denotes the Caputo type fractional derivative of order α defined as Dα∗ y(t) =

1 Γ(r − α )

Z t 0

(t − s)r−α −1 y(r) (s)ds,

t > 0,

for r − 1 ≤ α < r, r ∈ N. Jtα denotes the Riemann-Liouville fractional integral operator of order α , defined by ( 1 Rt α −1 y(s)ds, α > 0, α Γ(α ) 0 (t − s) Jt y(t) = y(t), α = 0. For convenience in calculations, we set Q(t) =

n−1 γ

∑

j j

j!

j=0

t .

(1.4)

Then the function ye(t) = y(t) − Q(t) satisfies the equation α D∗ p ye(t) +

p−1

∑

ai (t)D∗αi ye(t) =

fe(t) +

Z t

i=0

0

K(t, s)e y(s)ds,

0 < t ≤ T,

(1.5)

with the homogeneous initial conditions ye(i) (0) = 0, where fe(t) = f (t) −

p−1

∑

i = 0, · · · , n − 1.

ai (t)D∗αi Q(t) +

i=0

Z t 0

K(t, s)Q(t)ds.

(1.6)

It is well-known that (see [1, 2]) ( γj j−αi , 0 ≤ α ≤ n − 1, ∑n−1 i αi Γ(1+ j−αi ) t j=d α e i D∗ Q(t) = 0, αi > n − 1, where the ceiling function dα e denotes the smallest integer greater than or equal to α . α α Set z = D∗ p ye, then ye = Jt p z. Problem (1.1)-(1.2) may be rewritten as p−1

z(t) = − ∑

α −α ai (t)Jt p i z(t) +

i=0

fe(t) +

Z t 0

K(t, s)e y(s)ds.

(1.7)

In order to characterize the behavior of high order derivatives of a solution of Eqs.(1.1)-(1.2), we introduce a weighted space of smooth functions Cr,µ (0, T ] (see [3, 4]). For given r ∈ N and 2

−∞ < µ < 1, Cr,µ (0, T ] is defined as the set of all r times continuously differentiable functions v : (0, T ] −→ R such that for all t ∈ (0, T ] and i = 1, · · · , r the following estimates hold:   i < 1 − µ, 1, (i) |v (t)| ≤ C 1 + | logt|, i = 1 − µ ,   1−µ −i t , i > 1 − µ. Remark 1.1. Note that Cr,µ (0, T ] is a Banach space and Cm [0, T ] ⊆ Cm,µ (0, T ] ⊆ Cr,ν (0, T ] for m ≥ r ≥ 1 and µ ≤ ν < 1. Conversely, if u ∈ Cm,v (0, T ] with m ≥ 2 and v < 1 − k, k ∈ {1, · · · , m − 1}, then u ∈ Ck [0, T ]. The existence and regularity of the solution of problem (1.1)-(1.2) is described in the following lemma which can be proved similarly to the proof of Theorem 2.1 in [3, 4]. Lemma 1.1. Let (1.3) be true and assume that ai ∈ Cr,µ (0, T ], i = 0, · · · , p − 1, f ∈ Cr,µ (0, T ] and K ∈ Cr ([0, T ] × [0, T ]), where r ∈ N and −∞ < µ < 1. Then problem (1.1)-(1.2) possesses a unique α solution y ∈ Cn−1 [0, T ] such that D∗ p y ∈ Cr,ν (0, T ], where ( i f Q , 0, max{µ , ν1 , ν2 } ν= (1.8) max{µ , ν1 } i f Q = 0, with

ν1 = max{1 − α p + αi : α p − αi < N, i = 0, · · · , p − 1}, ν2 = max{1 − dαi e + αi : αi < n − 1, αi < N0 , i = 0, · · · , p − 1},

(1.9)

Here N0 = N ∪ {0} and Q is defined by (1.4). If for all indices i = 0, · · · , p − 1 we have α p − αi ∈ N then ν = max{µ , ν2 } for Q , 0 and ν = µ for Q = 0. Analogously, if we have αi ∈ N0 for all indices i = 0, · · · , p − 1 such that αi < n − 1 then ν = max{µ , ν1 } or ν = µ if, in addition, α p − αi ∈ N, i = 0, · · · , p − 1. The problems of fractional differential equations (FDEs) have gained much interest in different research areas and engineering applications. In particular, fractional integro-differential equations have widely occurred in the mathematical modeling of various physical phenomena, such as heat conduction in materials with memory, combined conduction, convection and radiation problems [1, 2, 5]. The solution of FDEs is numerically challenging because of the singularity behavior at the origin. In recent years, there are several analytic methods which have been developed to solve linear and nonlinear problems [6, 7, 8, 9, 10]. In the last decade or so, extensive research has been carried out on the development of numerical methods which are numerically stable for both linear and nonlinear fractional differential equations. Diethelm et al. [11] presented a predictor-corrector method for numerical solution of FDEs. This method takes advantage of FDEs converting into Volterra integral equations. Fractional linear multi-step methods for solving nonlinear fractional differential equations has been studied in [12]. Furthermore, the authors in [13, 14] transformed the space fractional partial differential equation into a system of ordinary differential equations, 3

which is then solved by the method of lines. In [3, 4] Pedas and Tamme developed spline collocation methods for solving the linear multi-order FDEs. Spectral methods are a class of important method for the numerical solution of differential equations. Spectral methods have excellent error properties and they offer exponential rates of convergence for smooth problems. In [15], Chebyshev spectral methods were proposed to solve nonlinear Volterra-Hammerstein integral equations. In [16], a Legendre-collocation method was proposed for Volterra integral equations (VIEs) of the second kind and the spectral accuracy can be justified both theoretically and numerically. Also, the Legendre-collocation method has been successfully applied to pantograph-type delay differential equations [17, 18] and Volterra integrodifferential equations (VIDEs) [19]. In [20, 21, 22] the direct solution techniques for solving linear and nonlinear multi-order FDEs using the spectral tau and collocation methods are introduced. There has been considerable interest in solving integral and differential equations using JacobiGauss collocation spectral method. Bhrawy and Alofi [23] proposed a shifted Jacobi-Gauss collocation spectral method for solving the nonlinear Lane-Emden type equation. Doha et al.[24] used the shifted Jacobi-Gauss collocation method for solving nonlinear high order multi-point boundary value problems. Bhrawy and Alghamdi [25] constructed a Jacobi-Gauss-Lobatto collocation method for solving the nonlinear fractional Langevin equation with three-point boundary conditions. Recently, Chen and Tang [26, 27] developed a Jacobi-collocation spectral method to solve the weakly singular VIEs and then Wei and Chen [28] extended the method to the weakly singular VIDEs. The spectral accuracy of the approaches is verified both theoretically and numerically in [26, 27, 28]. So far, very few works have provided theoretical results to justify the high accuracy of the spectral approximations for FIDEs. In this work, we will consider the special case that the exact solutions of (1.1)-(1.2) are smooth (see also [20, 26, 28, 29]). The main purpose of this paper is to develop a Jacobi-collocation method for linear fractional integro-differential equations based on the work of [27, 28]. We will provide a rigorous error analysis for the proposed method. We mention that [28] has considered the Jacobi collocation method for solving weakly singular VIDEs. In that reference, the solution y and its first order derivative y0 are respectively approximated by two independent Lagrange interpolation polynomials uN and vN of degree N. This implies that vN differs from the exact first order derivative of uN . Consequently, the corresponding errors of y and y0 are in the same order. α However, in this paper we only approximate the highest order derivative D∗ p ye by Lagrange interpolation polynomial uN of degree N, then the approximation of the solution ye is obtained by vN = J α p uN . This means that vN is a polynomial of degree N + α p when α p ∈ N, and vN is not a polynomial when α p < N. As will be seen from numerical experiments in Section 5 this could lead to a better approximation to the solution y(t). This paper is organized as follows. In Section 2 we introduce the Jacobi-collocation spectral approaches for the linear problem (1.1)-(1.2). Some lemmas useful for the convergence analysis will be provided in Section 3. The convergence analysis in the infinity norm will be given in Section 4. The performance of the corresponding iterated collocation solution is also investigated and it is shown that the iterated solution has an improvement on the original collocation solution. Finally, in Section 5 the obtained theoretical results are verified by some numerical experiments.

4

2. Jacobi-collocation spectral method In this section we derive the Jacobi-collocation method and describe its implementation for the α ,β problem (1.1)-(1.2). For α , β > −1, as defined in [30, 31], the Jacobi polynomials Jn (x), n = 0, 1, . . . , are the eigenfunctions of the singular Sturm-Liouville problem ¶ µ d 2 α ,β d α ,β Jn (x) = λnα ,β ω α ,β Jnα ,β (x), (1 − x )ω − dx dx α ,β

where the weight function ω α ,β and eigenvalues λn

are given by

ω α ,β (x) = (1 − x)α (1 + x)β , λnα ,β = n(n + α + β + 1). We define the weighted space 2 Lω α ,β (−1, 1) = {u | u is measurable and kukω α ,β < ∞},

equipped with the following inner product and norm (u, v)ω α ,β =

Z 1 −1

1

u(x)v(x)ω α ,β (x)dx, kukω α ,β = (u, u)ω2 α ,β .

In order to apply the spectral method on the standard interval [−1, 1], take t=

2t T (1 + x), x = − 1, µi = α p − αi − dα p − αi e, i = 0, · · · , p, T 2

(2.1)

to rewrite (1.7) as Z

u(x) = g(x) +

0

T 2 (1+x)

p−1

T T T y(s)]ds, (2.2) [ ∑ ( (1 + x) − s)µi Ki ( (1 + x), s)z(s) + K( (1 + x), s)e 2 2 i=0 2

where ai (t)(t − s)dα p −αi e−1 T T , i = 0, · · · , p − 1. u(x) = z( (1 + x)), g(x) = fe( (1 + x)), Ki (t, s) = − Γ(α p − αi ) 2 2 Moreover, for converting the integral interval [0, T (1 + x)/2] to the interval [−1, x], we make a linear transformation T (2.3) s = (1 + τ ), τ ∈ [−1, x]. 2 Then (2.2) becomes ( Rx p−1 ei (x, τ )u(τ ) + K ep (x, τ )v(τ )]d τ , [∑i=0 (x − τ )µi K u(x) = g(x) + −1 R (2.4) x (x − τ )α p −1 u(τ )d τ , v(x) = ( T2 )α p Γ(α1 p ) −1 5

where ep (x, τ ) = ( T )K( T (1 + x), T (1 + τ )). ei (x, τ ) = ( T )1+µi Ki ( T (1 + x), T (1 + τ )), i = 0, · · · , p − 1, K K 2 2 2 2 2 2 (2.5) For a given positive integer N, we denote the collocation points by {x j }Nj=0 which is the set of (N + 1) Jacobi-Gauss points with respect to the Jacobi weight function ω α ,β . In this paper, let µ be an arbitrarily fixed number in (−1, 0]. We take the special collocation points {x j }Nj=0 which are corresponding to the case that α = β = µ . Assume that (2.4) holds at x j , j = 0, · · · , N, u(x j ) = g(x j ) +

Z x j p−1

ei (x j , τ )u(τ ) + K ep (x j , τ )v(τ )]d τ . [ ∑ (x j − τ )µi K

−1 i=0

(2.6)

For approximating the integral terms in (2.6) we will transfer the integration interval [−1, x j ] to the unit interval [−1, 1]. More precisely, we make a simple linear transformation

τ = τ (x, θ ) =

x−1 1+x , θ+ 2 2

θ ∈ [−1, 1].

Then Z x j p−1

ei (x j , τ )u(τ ) + K ep (x j , τ )v(τ )]d τ [ ∑ (x j − τ )µi K

−1 i=0 Z 1 p−1

=

(2.7)

[ ∑ (1 − θ )µi K i (x j , θ )u(τ (x j , θ )) + K p (x j , θ )v(τ (x j , θ ))]d θ ,

−1 i=0

where

1 + x 1+µi e ) Ki (x, τ (x, θ )), i = 0, 1, · · · , p. 2 By the (N + 1)-point Gauss quadrature formula relative to the Jacobi weight with α = µi and β = 0, we get K i (x, θ ) = (

Z 1 −1

N

µi

(1 − θ ) K i (x j , θ )u(τ (x j , θ ))d θ ≈

∑ K i(x j , θi,k )u(τ (x j , θi,k ))ωi,k , 0 ≤ i ≤ p,

k=0

N where {θi,k }N k=0 is the set of (N + 1) Jacobi Gauss points, and {ωi,k }k=0 denote the corresponding weights. From (2.4) we have

v(x) =

Z 1 −1

(1 − θ )µ0 K p+1 (x, θ )u(τ (x, θ ))d θ ≈

N

∑ K p+1(x, θ0,k )u(τ (x, θ0,k ))ω0,k ,

k=0 n−1

α p (1−θ ) where K p+1 (x, θ ) = ( T2 )α p ( 1+x 2 ) Γ(α p ) .

6

(2.8)

Let PN denote the space of all polynomials of degree not exceeding N. For any u(x) ∈ C[−1, 1], we define the Lagrange interpolation polynomial N

α ,β

IN u(x) =

∑ Fj (x)u(x j ),

j=0

N where {Fk (x)}N k=0 are the Lagrange basis functions corresponding to the non-uniform mesh {x j } j=0 . We use u j to indicate the approximate values for u(x j ), 0 ≤ j ≤ N. Then the Jacobi-collocation method is to seek an approximate solution uN (x) ∈ PN of the form

uN (x) =

N

∑ Fj (x)u j ,

(2.9)

j=0

where u j , j = 0, · · · , N are determined by the following discrete collocation conditions N

p−1

k=0

i=0

u j = g(x j ) + ∑ uk ( ∑ aijk + b jk ),

(2.10)

and j τi,k = τ (x j , θi,k ),

b jk =

aijk =

N

∑ K i(x j , θi,l ))Fk (τi,lj ))ωi,l ,

N

l=0 N

q=0

l=0

j j , θ0,l )Fk (τ (τ p,q , θ0,l ))ω0,l ). ∑ ω p,qK p(x j , θ p,q)( ∑ K p+1(τ p,q

Having determined the approximation uN (x) for problem (2.4), we can determine the approximations N 2t N z (t) = u ( − 1), v (x) = ∑ K p+1 (x, θ0,k )uN (τ (x, θ0,k ))ω0,k , T k=0 2t yeN (t) = vN ( − 1), yN (t) = Q(t) + yeN (t), T N

N

(2.11)

for the solutions of problems (1.7), (2.8), (1.5) and (1.1), respectively. It is well-known that the iterated solution may improve order of convergence for the original solution of the equations (see, e.g., [29, 32, 33]). For problem (1.7) a better approximation than zN (t) is an iterated approximation of the form zN it (t) = −

p−1

∑

α −α ai (t)Jt p i zN (t) +

i=0

7

fe(t) +

Z t 0

K(t, s)e yN (s)ds.

(2.12)

In the case that the coefficients ai , i = 1, · · · , p − 1 in (1.1) are constants, we can also construct α another discretization. In this case, if we denote fb(t) = Jt p fe(t), then (1.5) can be rewritten as ( α −α α p−1 ye(t) = − ∑i=0 ai Jt p i ye(t) + fb(t) + Jt p w(t), (2.13) Rt w(x) = 0 K(t, s)e y(s)ds. By using the linear transformation and tricks used before, (2.13) can be written as ( α −1 Rx p−1 ei (x, τ )u(τ ) + ( T )α p (x−τ ) p v(τ )]d τ , u(x) = g(x) + −1 [∑i=0 (x − τ )µi K 2 Γ(α p ) R ep (x, τ )u(τ )d τ , v(x) = x K

(2.14)

−1

where T T T u(x) = ye( (1 + x)), v(x) = w( (1 + x)), g(x) = fb( (1 + x)). 2 2 2 A similar treatment to the construction of the numerical scheme (2.10) yields N

p−1

k=0

i=0

u j = g(x j ) + ∑ uk ( ∑ aijk + c jk ),

0 ≤ j ≤ N,

(2.15)

where c jk =

N

N

q=0

l=0

j j , θ p,l )Fk (τ (τ0,q , θ p,l ))ω p,l ). ∑ ω0,qK p+1(x j , θ0,q)( ∑ K p(τ0,q

The approximated solution uN (x) of problem (2.14) can be obtained by scheme (2.15) together with the polynomial interpolation (2.9), and thus we define yN (t) = Q(t) + uN (

2t − 1) T

(2.16)

as the approximated solution of problem (1.1)-(1.2). 3. Preliminaries as

2 In this section, first we define the orthogonal projection operator πN,ω α ,β : Lω α ,β (−1, 1) −→ PN

(v − πN,ω α ,β v, φ )ω α ,β = 0,

2 ∀v ∈ Lω α ,β (−1, 1), φ ∈ PN .

For any nonnegative integer m, define 2 Hωmα ,β (−1, 1) = {v | ∂xk v ∈ Lω α ,β (−1, 1), 0 ≤ k ≤ m},

with the norm

m

|v|H m;N

ω α ,β

(−1,1) = (

∑

k=min(m,N+1)

8

1

k∂xk vk2ω α ,β ) 2 .

In order to study the properties of our method, we need some useful lemmas which play a significant role in the convergence analysis later. From Theorem 3.35, Theorem 3.41 of [31] and (5.5.28) in [30] we have the following Lemma 3.1 and Lemma 3.2. The results of Lemma 3.3 and Lemma 3.4 on Lagrange interpolation based on the zeros of the associated Jacobi polynomials and on suitable additional nodes are given in [34] and [35], respectively. Throughout this paper, C denotes a positive constant that is independent of N and may have different values in different occurrences. Lemma 3.1. Let α , β > −1. For any function v ∈ Hωmα ,β (−1, 1), we have k∂xk (v − πN,ω α ,β v)kω α +k,β +k ≤ CN k−m k∂xm vkω α +m,β +m , 0 ≤ k ≤ m ≤ N + 1.

(3.1)

Lemma 3.2. Let α , β > −1. For any function v ∈ Hωmα ,β (−1, 1), we have α ,β

kv − IN vkω α ,β ≤ CN −m k∂xm vkω α +m,β +m , and

− 1 ,− 12

kv − IN 2

(3.2)

1

vk∞ ≤ CN 2 −m |v|H m;N , c (−1,1)

(3.3)

ω

where ω c is the Chebyshev weight, i.e.,α = β = − 12 . Lemma 3.3. Assume that {Fj (x)}Nj=0 are the corresponding N-th Lagrange interpolation polynomials associated with the Gauss points of the Jacobi polynomials. Then ( N O(log N), − 1 ≤ α , β ≤ − 21 , α ,β kIN k∞ = max ∑ |Fj (x)| = 1 x∈[−1,1] j=0 O(N γ + 2 ), γ = max(α , β ), otherwise. Lemma 3.4. Assume that {Fj (x)}Nj=0 are the corresponding N-th Lagrange interpolation polynomials associated with the Jacobi collocation points {x j }Nj=0 . Then for every bounded function v(x), there exists a constant C independent of v such that N

sup k ∑ v(x j )Fj (x)kω α ,β ≤ Ckvk∞ . N

j=0

Lemma 3.5. If v, φ ∈ Hωmα ,β (−1, 1) for some m ≥ 1, then for the Jacobi-Gauss integration, we have α ,β

|(v, φ )ω α ,β − (v, φ )N,ω α ,β | ≤ CN −m k∂xm vkω α +m,β +m kφ kω α ,β +Ckvk∞ kφ − IN φ kω α ,β , where

(3.4)

N

(v, φ )N,ω α ,β =

∑ v(xk )φ (xk )ωk .

k=0

9

(3.5)

Proof. Note that the discrete inner product is based on the (N + 1)-degree Jacobi-Gauss points corresponding to the weight function ω α ,β . We have α ,β

α ,β

(v, φ )N,ω α ,β = (IN v, IN φ )ω α ,β . Consequently, we have α ,β

α ,β

α ,β

|(v, φ )ω α ,β − (v, φ )N,ω α ,β | = |(v − IN v, φ )ω α ,β + (IN v, φ − IN φ )ω α ,β | α ,β

α ,β

α ,β

≤ kv − IN vkω α ,β kφ kω α ,β + kIN vkω α ,β kφ − IN φ kω α ,β , which together with Lemma 3.2 and Lemma3.4 leads to the desired estimate (3.4). ¤ The following two lemmas of generalized Gronwall inequality and Hardy’s inequality for singular kernels can be found in [26, 27], which are essential for establishing our main results. Lemma 3.6. Suppose L > 0, 0 < µ < 1 and v(t) is a non-negative, locally integrable function defined on [0, T ] satisfying u(t) ≤ v(t) + L

Z t 0

(t − s)−µ u(s)ds.

Then there exists a constant C such that u(t) ≤ v(t) +CL

Z t 0

(t − s)−µ v(s)ds.

Lemma 3.7. For all measurable function f ≥ 0, the following generalized Hardy’s inequality Z b

(

a

q

1/q

|(K f )(x)| ω1 (x)dx)

Z b

≤ C(

a

| f (x)| p ω2 (x)dx)1/p

holds if and only if Z b

sup ( a
x

Z x

ω1 (t)dt)1/q (

a

0

0

ω21−p (t)dt)1/p < ∞, p0 =

p p−1

for the case 1 < p ≤ q < ∞. Here K is an operator of the form (K f )(x) =

Z x a

ρ (x,t) f (t)dt

with ρ (x,t) a given kernel, ω1 , ω2 weight functions, and −∞ ≤ a < b ≤ ∞. For given r ∈ N and 0 ≤ κ ≤ 1, we denote by Cr,κ [−1, 1] the set of continuous function v : [−1, 1] −→ R whose r-th derivatives are H older ¨ continuous with exponent κ , equipped with norm k.kr,κ (see, e.g.,[29, 31]). From[36] we know that there exists a constant Cr,κ ≥ 0 such that for any function v ∈ Cr,κ [−1, 1], there exists a polynomial function kN v ∈ PN such that kv − kN vk∞ ≤ Cr,κ N −r−κ kvkr,κ . 10

(3.6)

4. Convergence analysis In order to simplify the notations and without lose of generality, we consider the case p = 1 in the convergence analysis of the numerical scheme. We first give an error estimate for the numerical solution uN of the problem (2.4) in the sense of infinity norm. Theorem 4.1. Let u be the exact solution of the problem (2.4), which is assumed to be sufficiently smooth. Let the approximated solution uN be given by the spectral collocation scheme (2.10) together with the polynomial interpolation (2.9). If u ∈ Hωmµ ,µ (−1, 1), then for sufficiently large N the following error estimate holds  1 1 CN 2 −m+µ (K ∗ kuk∞ + k∂xm vkω m+µ1 ,m + N 21 |u| m;N , Hω c (−1,1) ), k1 ≤ − µ < 2 (4.1) ku − uN k∞ ≤ 1 1 CN −m log N(K ∗ kuk∞ + k∂xm vkω m+µ1 ,m + N 2 |u| m;N ≤ − µ < 1, ), 2 H c (−1,1) ω

where K∗ =

e j (xi , τ (xi , θ ))k m+µ j ,m , k1 = max{0, −µ0 − 1 }. k∂θm K ω 0≤i≤N, 0≤ j≤1 2 max

Proof. From (2.6)-(2.11) and the definition of the weighted inner product and discrete inner product (3.5) we get u(xi ) = g(xi ) + (K 0 (xi , θ ), u(τ (xi , θ )))ω µ0 ,0 + (K 1 (xi , θ ), v(τ (xi , θ )))ω µ1 ,0 , v(x) = (K 2 (x, θ ), u(τ (x, θ )))ω µ0 ,0 ,

(4.2)

and ui = g(xi ) + (K 0 (xi , θ ), uN (τ (xi , θ )))N,ω µ0 ,0 + (K 1 (xi , θ ), vN (τ (xi , θ )))N,ω µ1 ,0 , vN (x) = (K 2 (x, θ ), uN (τ (x, θ )))ω µ0 ,0 .

(4.3)

For convenience, we define integral operators Tei , 1 ≤ i ≤ 3, as (Tei u)(x) =

Z x −1

(x − s)

µi−1

ei−1 (x, s)u(s)ds, i = 1, 2, (Te3 u)(x) = K

Z x (x − τ )α1 −1 T u(s)ds. ( )α1 −1

2

Γ(α1 )

Subtracting (4.3) from (4.2) yields 1 2 u(xi ) − ui = Te1 (u − uN )(xi ) + Te2 (v − vN )(xi ) + Ii,2 + Ii,2 ,

(4.4)

where v − vN = (K 2 (x, θ ), u(τ (x, θ )) − uN (τ (x, θ )))ω µ0 ,0 , 1 + xi 1+µ0 e 1 e0 (xi , τ (xi , θ )), uN (τ (xi , θ ))) µ0 ,0 ), ) ((K0 (xi , τ (xi , θ )), uN (τ (xi , θ )))ω µ0 ,0 − (K Ii,2 =( N,ω 2 1 + xi 1+µ1 e 2 e1 (xi , τ (xi , θ )), vN (τ (xi , θ ))) µ1 ,0 ). ) ((K1 (xi , τ (xi , θ )), vN (τ (xi , θ )))ω µ1 ,0 − (K Ii,2 =( N,ω 2 11

Multiplying Fi (x) on both sides of (4.4) and summing up from i = 0 to i = N yields µ ,µ µ ,µ u − uN =u − IN u − (I − IN )(Te1 (u − uN ) + Te2 (v − vN )) N

1 2 + Te1 (u − uN ) + Te2 (v − vN ) + ∑ Fi (x)(Ii,2 + Ii,2 ). i=0

Due to (2.8) and using the Dirichlets formula we have |Te2 (v − vN )| ≤ C ≤C ≤C ≤C

Z x

|v − vN |ds,

Z−1 x Z s

−1 Z−1 x Z x

(ζ − s)α p −1 d ζ |u − uN |ds,

s Z−1 x −1

(s − ζ )α p −1 |u − uN |d ζ ds,

(x − s)α p |u − uN |ds.

It is clear that N

|u − u | ≤|I1 | + |I2 | + |I3 | + |I4 | +C

Z x −1

(x − s)µ0 |u − uN |ds,

(4.5)

where µ ,µ I1 = u − IN u, I2

N

=∑

1 Fi (x)Ii,2 , I3

N

2 = ∑ Fi (x)Ii,2 ,

i=0 i=0 µ ,µ e N e I4 = (I − IN )(T1 (u − u ) + T2 (v − vN )).

(4.6)

It follows from Lemma 3.6 that 4

N

ku − u k∞ ≤ C ∑ kIi k∞ .

(4.7)

i=1

By using Lemma 3.2 and Lemma 3.3 we obtain that µ ,µ

µ ,µ

− 1 ,− 12

kI1 k∞ = ku − IN uk∞ = k(I − IN )(u − IN 2 µ ,µ

u)k∞

− 12 ,− 12

≤ (1 + kIN k∞ )ku − IN uk∞  CN 1−m+µ |u| m;N Hω c (−1,1) , ≤ 1 CN 2 −m log N|u| m;N H c (−1,1) , ω

0 ≤ −µ < 12 , 1 2

(4.8)

≤ −µ < 1.

The implementation of Lemma 3.5 implies 1 e0 (xi , τ (xi , θ ))kω m+µ0 ,m kuN (τ (xi , θ ))k µ0 ,0 max |Ii,2 | ≤ CN −m max k∂θm K ω

0≤i≤N

≤ CN

−m

0≤i≤N ∗

K (kuk∞ + ku − uN k∞ ). 12

(4.9)

Next, it follows from Lemma 3.3 that N

µ ,µ

1 1 kI2 k∞ = k ∑ Fi (x)Ii,2 k∞ ≤ CkIN k∞ max |Ii,2 | 0≤i≤N

(i=0 1 CN 2 −m+µ K ∗ (kuk∞ + ku − uN k∞ ), ≤ CN −m log NK ∗ (kuk∞ + ku − uN k∞ ),

0 ≤ −µ < 12 , 1 2 ≤ − µ < 1.

(4.10)

Actually, by Lemma 3.5 2 e1 (xi , τ (xi , θ ))kω m+µ1 ,m kvN (τ (xi , θ ))k µ1 ,0 max |Ii,2 | ≤ CN −m max k∂θm K ω

0≤i≤N

0≤i≤N

1 + xi 1+µ1 e µ ,0 ) kK1 (xi , τ (xi , θ ))k∞ k(I − IN1 )vN (τ (xi , θ ))kω µ1 ,0 0≤i≤N 2 1 + xi 1+µ1 µ ,0 ≤ CN −m K ∗ kvN k∞ +C max ( ) k(I − IN1 )(v − v + vN )(τ (xi , θ ))kω µ1 ,0 0≤i≤N 2 1 + xi 1+µ1 m ≤ CN −m K ∗ kuN k∞ +CN −m max ( ) k∂θ v(τ (xi , θ ))kω m+µ1 ,m 0≤i≤N 2 +C max (

µ ,0

+C max k(I − IN1 )(v − vN )(τ (xi , θ ))kω µ1 ,0 . 0≤i≤N

(4.11) Direct calculations show that 1 + xi 1+µ1 m ) k∂θ v(τ (xi , θ )kω m+µ1 ,m 0≤i≤N 2 Z 1 1 + xi 1+µ1 xi 2 = max ( ) ( (xi − s)m+µ1 (1 + s)m (∂sm v)2 ds) 2 ≤ Ck∂xm vkω m+µ1 ,m . 0≤i≤N 2 −1 max (

(4.12)

From [37], we known that Tei , 1 ≤ i ≤ 3, are linear and compact operators from C[−1, 1]into C0,κ [−1, 1]. This implies that for any function u ∈ C[−1, 1], there exists a positive constant C such that kTei uk0,κ ≤ Ckuk∞ , 0 < κ < 1 + µ0 , 1 ≤ i ≤ 3. (4.13) Hence, by using (3.6) and Lemma 3.4 we have µ ,0

max k(I − IN1 )(v − vN )(τ (xi , θ ))kω µ1 ,0

0≤i≤N

µ ,0

= max k(I − IN1 )(I − kN )(v − vN )(τ (xi , θ ))kω µ1 ,0 0≤i≤N

N

−κ

(4.14)

N

≤ Ck(I − kN )(v − v )k∞ ≤ CN kv − v k0,κ ≤ CN −κ kTe3 (u − uN )k0,κ ≤ CN −κ ku − uN k∞ . A combination of (4.11), (4.12) and (4.14) yields 2 max |Ii,2 | ≤ CN −m (K ∗ kuN k∞ + k∂xm vkω m+µ1 ,m ) +CN −κ ku − uN k∞

0≤i≤N

≤ CN −m (K ∗ kuk∞ + k∂xm vkω m+µ1 ,m ) +CN −κ ku − uN k∞ . 13

(4.15)

Also, it follows from Lemma 3.3 that N

µ ,µ

2 2 kI3 k∞ = k ∑ Fi (x)Ii,2 k∞ ≤ CkIN k∞ max |Ii,2 | 0≤i≤N

(i=0 1 1 0 ≤ −µ < 21 , CN 2 −m+µ (K ∗ kuk∞ + k∂xm vkω m+µ1 ,m ) +CN 2 −κ +µ ku − uN k∞ , ≤ CN −m log N(K ∗ kuk∞ + k∂xm vkω m+µ1 ,m ) +CN −κ log Nku − uN k∞ , 21 ≤ −µ < 1. (4.16)

By using (3.6), (4.13) and Lemma 3.3 we have µ ,µ kI4 k∞ =k(I − I )(Te1 (u − uN ) + Te2 (v − vN ))k∞ N

µ ,µ ≤ k(I − IN )(Te1 (u − uN ) − kN Te1 (u − uN ) + Te2 (v − vN ) − kN Te2 (v − vN ))k∞ µ ,µ ≤ (1 + kI k∞ )(kTe1 (u − uN ) − kN Te1 (u − uN )k∞ + kTe2 (v − vN ) − kN Te2 (v − vN )k∞ ) N

µ ,µ

≤ C(1 + kIN k∞ )N −κ (kTe1 (u − uN )k0,κ + kTe2 (v − vN )k0,κ ) µ ,µ

≤ C(1 + kIN k∞ )N −κ (ku − uN k∞ + kv − vN k∞ ) ( 1 0 ≤ −µ < 21 , CN 2 −κ +µ ku − uN k∞ , ≤ 1 CN −κ log Nku − uN k∞ , 2 ≤ − µ < 1. (4.17) Now we can choose 12 + µ < κ < 1+ µ0 when k1 ≤ −µ < 21 , and 0 < κ < 1+ µ0 when 21 ≤ −µ < 1, then the estimate (4.1) follows from (4.7), (4.8), (4.10), (4.16) and (4.17). ¤ Remark 4.1. Note that the convergence analysis of the proposed method (2.15) can be established in a similar way as Theorem 4.1. Next, we analyze the errors of yN and zN it . The following result shows that they have higher N order accuracy than u . Theorem 4.2. Let u and uN be the same as those in Theorem 4.1. Let y be the exact solution of problem (1.1)-(1.2) and the approximated solution yN is given by the spectral collocation scheme (2.11). Let z be the exact solution of the problem (1.7) and zN it is obtained by the iterated spectral collocation scheme (2.12), then ( CN −m (N µ U1 + N 1−κ +µ U2 +U3 ), k1 ≤ −µ < 21 , (i) (i) (4.18) max ky − yN k∞ ≤ 1 1 0≤i≤q CN −m (N − 2 log NU1 + N 2 −κ log NU2 +U3 ), 21 ≤ −µ < 1, ( CN −m (N µ U1 + N 1−κ +µ U2 +U3 ), k2 ≤ −µ < 12 , (4.19) kz − zN k ≤ 1 1 it ∞ k3 ≤ −µ < 1, CN −m (N − 2 log NU1 + N 2 −κ log NU2 +U3 ), for sufficiently large N and for any κ ∈ (0, 1 + µ0 ), where 1 1+µ , k2 = max{0, −2µ0 − 1}, k3 = max{ , −2µ0 − 1}, q = α1 − 2 2 ∗ 0 c − 21 0 , U1 = K (ku kω + N kv kω 1,1 ), U2 = |u|H m;N c (−1,1) ω

U3 = k∂xm ukω m−µ ,m−µ

+ k∂xm vkω m,m + K ∗ kuk∞ . 14

Proof. For i = 0, 1, · · · , n − 1 we have (i)

ky(i) − yN k∞ = kJtα1 −i (z − zN )k∞ ≤ Ck

Z x −1

(x − s)α1 −i−1 (u − uN )k∞ . µ −1 2

Using the Cauchy-Schwartz inequality, we find that for α ≥ Z x

(

−1

(x − s)α v(s)ds)2 ≤ C

Z x −1

(1 + s)µ v2 (s)ds

≤ C(x + 1)2α −µ +1 ≤C

Z 1 −1

Z 1 −1

Z x −1

(4.20)

2 and v ∈ Lω µ ,µ (−1, 1)

(x − s)2α (1 + s)−µ ds

(1 + s)µ v2 (s)ds

(4.21)

2 (1 − s)µ (1 + s)µ v2 (s)ds ≤ Ckvkω µ ,µ .

If α1 − i ≥ 1+2 µ (which is equivalent to i ≤ α1 − 1+2 µ ), by using the Hardy inequality in Lemma 3.7 we have (i)

ky

(i) − yN k∞

4

N

≤ Cku − u kω µ ,µ ≤ C ∑ kIi kω µ ,µ .

(4.22)

i=1

Using (4.6) and Lemma 3.2 gives µ ,µ

kI1 kω µ ,µ = ku − IN ukω µ ,µ ≤ CN −m k∂xm ukω m+µ ,m+µ .

(4.23)

Next, by using (4.9), (4.15) and the convergence result in Theorem 4.1 we obtain N

N

1

1 1 2 2 kI2 kω µ ,µ = k ∑ Fi (x)Ii,2 kω µ ,µ = ( ∑ ωk (Ik,2 ) ) i=0

k=0 1 −m ∗ |Ik,2 | ≤ CN K (kuk∞ + ku − uN k∞ )

≤ C max 0≤k≤N ( 1 CN −m K ∗ (N µ ku0 kω c + N µ − 2 kv0 kω 1+µ1 ,1 + kuk∞ ), k1 ≤ −µ < 12 , ≤ 1 CN −m K ∗ (N − 2 log Nku0 kω c + N −1 log Nkv0 kω 1+µ1 ,1 + kuk∞ ), 12 ≤ −µ < 1, (4.24) and N

2 kI3 kω µ ,µ = k ∑ Fi (x)Ii,2 kω µ , µ i=0 −m

≤ CN (K ∗ kuk∞ + k∂xm vkω m+µ1 ,m ) +CN −κ ku − uN k∞  CN −m (K ∗ kuk∞ + k∂xm vkω m+µ1 ,m + N 1−κ +µ |u| m;N k1 ≤ −µ < 21 , Hω c (−1,1) ), ≤ 1 1 CN −m (K ∗ kuk∞ + k∂xm vkω m+µ1 ,m +CN 2 −κ log N|u| m;N H c (−1,1) ), 2 ≤ − µ < 1. ω

(4.25) 15

Finally, it follows from (3.6), (4.13), (4.1) and Lemma 3.4 that µ ,µ kI4 kω µ ,µ = k(I − IN )(Te1 (u − uN ) + Te2 (v − vN ))kω µ ,µ µ ,µ ≤ k(I − I )(Te1 (u − uN ) − kN Te1 (u − uN ) + Te2 (v − vN ) − kN Te2 (v − vN ))kω µ ,µ N

≤ C(kTe1 (u − uN ) − kN Te1 (u − uN )k∞ + kTe2 (v − vN ) − kN Te2 (v − vN )k∞ ) ≤ CN −κ (kTe1 (u − uN )k0,κ + kTe2 (v − vN )k0,κ ) ≤ CN −κ ku − uN k∞  1 CN 2 −m−κ +µ (K ∗ kuk∞ + k∂xm vkω m+µ1 ,m + N 12 |u| m;N Hω c (−1,1) ), ≤ 1 CN −m−κ log N(K ∗ kuk∞ + k∂xm vkω m+µ1 ,m + N 2 |u| m;N H c (−1,1) ), ω

k1 ≤ −µ < 12 , 1 2

≤ −µ < 1, (4.26)

which together with (4.22), (4.23), (4.24) and (4.25) yields the estimate (4.18). Applying the transformation (2.1) to (2.12) we have e N e N uN it = g + T1 u + T2 v , where

(4.27)

N T uN it (x) = yit ( (1 + x)). 2

Subtracting (4.27) from (2.4) yields N N e e u − uN it = T1 (u − u ) + T2 (v − v ) µ ,µ N e e e = Te1 u − Te1 I uN it + T1 I2 + T1 I3 + T2 (v − v )

N µ ,µ N e e e e e = T1 I1 + T1 IN (u − uN it ) + T1 I2 + T1 I3 + T2 (v − v ).

(4.28)

µ ,µ From (4.17) we have k(I − IN )Te1 k∞ → 0 as N → ∞. Since the operator I − Te1 is invertible, it follows from the Banach Lemma ([38]) that there exists an integer N0 such that for all N ≥ N0 the µ ,µ operator (I − IN Te1 ) is invertible and µ ,µ k(I − IN Te1 )−1 k∞ ≤ C.

The identity

µ ,µ µ ,µ µ ,µ (I − Te1 IN )−1 = I + Te1 (I − IN Te1 )−1 IN ,

µ ,µ then shows that I − Te1 IN is invertible for N sufficiently large, and moreover µ ,µ µ ,µ µ ,µ k(I − Te1 IN )−1 k∞ ≤ CkTe1 |IN |k∞ ≤ CkIN kω µ ,µ ≤ C

for µ ≤ 1 + 2µ0 , where we have used (4.21) and Lemma 3.4. Then N N kz − zN it k∞ = ku − uit k∞ ≤ C(kI1 kω µ ,µ + kI2 kω µ ,µ + kI3 kω µ ,µ + ky − y k∞ ),

µ ≤ 1 + 2 µ0 .

Since 1 + 2µ0 = 1 + 2(α1 − dα1 e) ≤ 2α1 − 1, we have q = α1 − 1+2 µ ≥ 0 when µ ≤ 1 + 2µ0 . Therefore, a combination of (4.18), (4.23), (4.24) and (4.25) yields the estimate (4.19). ¤ 16

5. Numerical experiments In this section some numerical examples are considered to demonstrate the efficiency and accuracy of the proposed method in the infinity norm. In all examples, we set µ = µ0 . Example 5.1. Consider the following multi-term fractional integro-differential equation:  Z t  y00 (t) + 3 √tD3/2 y(t) + (1 − t)y(t) = f (t) + s2 y(s)ds ∗ , 0 ≤ t ≤ 1, (5.1) 0  0 y (0) = y(0) = 1 with f (t) = 3E1,1/2 (t) −

3 + (t − t 2 )et + 2. Γ(1/2)

In this problem, the approximated solutions zN , zN it and yN are obtained by the numerical scheme (2.9)-(2.12). It can be verified that the exact solution of (5.1) is y(t) = et and thus the corresponding integral equation has a solution z = y00 = et . Table 1 shows the numerical errors of the problem (5.1) for 2 ≤ N ≤ 16 in infinity norm. We also plotted the numerical errors in Figure 1. As we can see from Table 1 and Figure 1 the errors decay rapidly, which is confirmed by spectral accuracy. This is in accordance with our theoretical results because the exact solution is sufficiently smooth. Table 1: The errors of Jacobi-collocation spectral method for example 5.1

N ky − yN k∞ kz − zN k∞ kz − zN it k∞ N ky − yN k∞ kz − zN k∞ kz − zN it k∞

2 4.4866e-04 9.0766e-03 1.0049e-02 10 5.5511e-16 2.2204e-14 1.3323e-14

4 3.8460e-07 3.3760e-05 2.3031e-05 12 2.7756e-16 3.1086e-15 3.5527e-15

6 3.3647e-10 5.3117e-08 2.9601e-08 14 2.7756e-16 3.1086e-15 2.4425e-15

8 1.7486e-13 3.6895e-11 2.2729e-11 16 3.3307e-16 2.6645e-15 2.6645e-15

Example 5.2. Consider the following linear fractional integro-differential equation:  Z t  D5/3 y(t) + y0 (t) + ty(t) = f (t) + 1 (t − s)2 y(s)ds ∗ 2 0 , 0 ≤ t ≤ 1.  0 y (0) = y(0) = 0

(5.2)

First we implement the numerical schemes (2.9)-(2.12) to solve this problem and get the apN proximated solutions zN , zN it and yN . Moreover, the approximated solutions y are obtained by the numerical scheme (2.15). 10

Case 1: If the function f (t) is selected so that the exact solution to (5.2) is y(t) = 11

2t 3 Γ(14/3)

2/3

4

t , then the corresponding integral equation has a solution z = D∗ y(t) = − 12

17

5

t3 Γ(13/3)

t3 Γ(8/3)

+

+ t2 −

0

10

||y−yN||∞ −2

||z−zN||∞

10

||z−zN || it ∞

−4

10

−6

10

−8

10

−10

10

−12

10

−14

10

−16

10

2

4

6

8

10 2≤ N≤16

12

14

16

Figure 1: The errors of example 5.1 versus the number of collocation points.

Table 2: The errors of Jacobi-collocation spectral method for example 5.2 (case 1)

N ky − yN k∞ ky − yN k∞ kz − zN k∞ kz − zN it k∞ N ky − yN k∞ ky − yN k∞ kz − zN k∞ kz − zN it k∞

2 6.7420e-03 7.6802e-04 1.7157e-02 2.9819e-03 10 1.3137e-07 1.0405e-07 9.4937e-05 1.9740e-06

4 6.7211e-05 1.3923e-05 1.6224e-03 1.0891e-04 12 4.0249e-08 4.0458e-08 5.3610e-05 8.0857e-07

18

6 3.7594e-06 1.6134e-06 4.6277e-04 1.7767e-05 14 1.4797e-08 1.8232e-08 3.2971e-05 2.1753e-07

8 5.6136e-07 3.3830e-07 1.9015e-04 4.8154e-06 16 6.2094e-09 9.0683e-09 2.1593e-05 1.1394e-07

−1

10

||y−yN||∞

−2

10

||y−yN||∞ ||z−zN||∞

−3

10

||z−zN || it ∞

−4

10

−5

10

−6

10

−7

10

−8

10

−9

10

2

4

6

8

10 2≤ N≤16

12

14

16

Figure 2: The errors of example 5.2 (case 1) versus the number of collocation points.

2t 7/3 Γ(10/3) .

In Table 2 the errors are given in the infinity norm for different values of N. Also the errors are depicted in Figure 2. It can be seen that for the numerical schemes (2.10) and (2.15) the errors decay algebraically as the exact solution for this example is not sufficiently smooth. But, the decay rate of the errors of the scheme (2.10) is slower than that of the scheme (2.15) because the original solution y(t) has better smoothness properties in comparison to its fractional derivative z(t). Case 2: If in this problem the function f (t) is selected so that the exact solution to (5.2) 5

is y(t) =

t3 Γ(8/3)

10

11

4

t 2t 3 t3 . As a result the corresponding integral equation has a − 12 + Γ(14/3) + Γ(13/3)

2/3

5

7/3

2t t3 . Similar to the previous examples, the errors are + t 2 − Γ(10/3) solution z = D∗ y(t) = 1 + Γ(8/3) demonstrated in Table 3 and Figure 3 for different values of N. Again we can see both schemes (2.10) and (2.15) give algebraical accuracy for not very large values of N. Although the original solution y(t) has the same smoothness properties as its fractional derivative z(t), however, the decay rate of the errors of the scheme (2.10) is faster than that of the scheme (2.15). From these tables and figures we can see that for the schemes (2.10) and (2.15) the error norms are quite small. Both methods give satisfactory numerical performance in terms of convergence. And each of them has its own advantages.

6. Concluding remarks In this work, a Jacobi collocation method for multi-order fractional integro-differential equations has been described. The spectral accuracy associated with infinity norm error estimates is demonstrated theoretically for the approximation methods. Several examples are given to support the theoretical results. The results show that the Jacobi collocation method is a competitive method for fractional integro-differential equations. Moreover, all the schemes presented also can be used for solving nonlinear fractional differential equations.

19

Table 3: The errors of Jacobi-collocation spectral method for example 5.2 (case 2)

N ky − yN k∞ ky − yN k∞ kz − zN k∞ kz − zN it k∞ N ky − yN k∞ ky − yN k∞ kz − zN k∞ kz − zN it k∞

2 4 6 8 3.7983e-03 1.2845e-03 3.9108e-04 1.6805e-04 7.6807e-04 1.3923e-05 1.6127e-06 3.3734e-07 1.7159e-02 1.6224e-03 4.6277e-04 1.9015e-04 2.9834e-03 1.0891e-04 1.7767e-05 4.8154e-06 10 12 14 16 8.6264e-05 4.96199e-05 3.0917e-05 2.0444e-05 1.0349e-07 4.0379e-08 1.8165e-08 9.0088e-09 9.4937e-05 5.3610e-05 3.2971e-05 2.1593e-05 1.9740e-06 8.0857e-07 2.1753e-07 1.1394e-07

−1

10

||y−yN||∞

−2

10

||y−yN||∞ ||z−zN||∞

−3

10

||z−zN || it ∞

−4

10

−5

10

−6

10

−7

10

−8

10

−9

10

2

4

6

8

10 2≤ N≤16

12

14

16

Figure 3: The errors of example 5.2 (case 2) versus the number of collocation points.

20

Acknowledgements The authors would like to thank the referees for their valuable comments and helpful suggestions which led to an improved version of the paper. References [1] I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. [2] K. Diethelm, The Analysis of Fractional Differential Equations, Springer-Verlag, Berlin, 2010. [3] A. Pedas, E. Tamme, On the convergence of spline collocation methods for solving fractional differential equations, J. Comput. Appl. Math. 235 (2011) 3502-3514. [4] A. Pedas, E. Tamme, Spline collocation methods for linear multi-term fractional differential equations, J. Comput. Appl. Math. 236(2011) 167-176. [5] K. Miller, B. Ross, An Introduction to the Fractional Calaulus and Fractional Differential Equations, John Wiley and Sons Inc. New York, 1993. [6] S. Momani, M. A. Noor, Numerical methods for fourth order fractional integro-differential equations. Appl.Math.Comput. 182(2006) 754-60. [7] A. Arikoglu, I. Ozkol, Solution of fractional integro-differential equations by using fractional differential transform method. Chaos, Solitons, & Fractals. 40(2009) 521-529. [8] D. Nazaria, S. Shahmorad, Application of the fractional differential transform method to fractional-order integro-differential equations with nonlocal boundary conditions. J.Comput.Appl.Math. 234(2010) 883-891. [9] E. A. Rawashdeh, Numerical solution of fractional integro-differential equations by collocation method. Appl.Math.Compt. 176(2006) 1-6. [10] L. Huang, X.F. Li, Y.L. Zhao and X.Y. Duan, Approximate solution of fractional integro-differential equations by Taylor expansion method. Comput. Math. Appl. 62(2011) 1127-1134. [11] K. Diethelm, Efficient solution of multi-term fractional differential equations using P(EC)mE methods, Computing 71 (2003) 305-319. [12] R. Lin, F. Liu, Fractional high order methods for the nonlinear fractional ordinary diffrential equation. Nonlinear Dynamics. 66(2007) 856-869. [13] F. Liu, V. Anh, I.Turner, Numerical Solution of the Space Fractional Fokker-Planck Equation. J. Comp. Appl. Math. 166(2004) 209-219. [14] Q. Yang, F. Liu, I. Turner, Numerical methods for fractional partial differential equations with riesz space fractional derivatives, Appl. Math. Model. 34(2010) 200-218. [15] G. Elnagar, M. Kazemi, Chebyshev spectral solution of nonlinear Volterra-Hammerstein integral equations, J. Comput. Appl. Math. 76 (1996) 147-158. [16] T. Tang, Xiang Xu, Jin Cheng, On spectral methods for Volterra type integral equations and the convergence analysis, J. Comput. Math. 26 (2008) 825-837. [17] I. Ali, H. Brunner, T. Tang, A spectral method for pantograph-type delay differential equations and its convergence analysis, J. Comput. Math. 27 (2009) 254-265. [18] I. Ali, H. Brunner, T. Tang, Spectral methods for pantograph-type differential and integral equations with multiple delays, Front. Math. China 4 (2009) 49-61. [19] Y. Jiang, On spectral methods for Volterra-type integro-differential equations, J. Comput. Appl. Math. 230 (2009) 333-340. [20] F. Ghoreishi, S. Yazdani, An extension of the spectral Tau method for numerical solution of multi-order fractional differential equations with convergence analysis, Comput. Math. Appl. 61 (2011) 30-43. [21] E.H. Doha, A.H. Bhrawy, S.S. Ezz-Eldien, Efficient Chebyshev spectral methods for solving multi-term fractional orders differential equations, Appl. Math. Model. 35 (2011) 5662-5672. [22] E.H. Doha, A.H. Bhrawy, S.S. Ezz-Eldien, A new Jacobi operational matrix: An application for solving fractional differential equations, Appl. Math. Modell. 36(2012) 4931-4943. [23] A.H. Bhrawy, A.S. Alofi, A Jacobi-Gauss collocation method for solving nonlinear Lane-Emden type equations, Commun Nonlinear Sci Numer Simulat 17 (2012) 62-70.

21

[24] E.H. Doha, A.H. Bhrawy, R.M. Hafez, On shifted Jacobi spectral method for high-order multi-point boundary value problems, Commun. Nonlinear Sci. Numer. Simul. 17 (2012) 3802-3810. [25] A.H. Bhrawy, M. A. Alghamdi, A shifted Jacobi-Gauss-Lobatto collocation method for solving nonlinear fractional Langevin equation involving two fractional orders in different intervals, Boundary Value Problems 2012 (2012) 62. [26] Y. Chen, T. Tang, Spectral methods for weakly singular Volterra integral equations with smooth solutions, J. Comput. Appl. Math. 233 (2009) 938-950. [27] Y. Chen, T. Tang, Convergence analysis for the Jacobi spectral-collocation methods for Volterra integral equations with a weakly singular kernel, Math. Comput. 79 (2010) 147-167. [28] Y.Wei,Y. Chen, Convergence Analysis of the Spectral Methods for Weakly Singular Volterra Integro-Differential Equations with Smooth Solutions. Adv. Appl. Math. Mech. 4 (2012) 1-20. [29] H. Brunner, Collocation Methods for Volterra Integral and Related Functional Differential Equations, Cambridge University Press, Cambridge MA 2004. [30] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods Fundamentals in Single Domains, Springer-Verlag, 2006. [31] J. Shen, T. Tang, L.L. Wang, Spectral Methods: Algorithms, Analyses and Applications, Springer, New York, 2011. [32] J. G. Blom, H. Brunner, The numerical solution of nonlinear Volterra integral equations of the second kind by collocation and iterated collocation methods, SIAM J. Sci. Statist. Comput. 8 (1987) 806-830. [33] H. Brunner, Iterated collocation methods and their discretizations for Volterra integral equations, SIAM J. Numer. Anal. 21 (1984) 1132-1145. [34] G. Mastroianni, D. Occorsio, Optimal systems of nodes for Lagrange interpolation on bounded intervals. A survey, J. Comput. Appl. Math. 134 (2001) 325-341. [35] P. Nevai, Mean convergence of Lagrange interpolation. III, Trans. Amer. Math. Soc. 282 (1984) 669-698. MR 85c:41009. [36] I.G. Graham, I.H. Sloan, Fully discrete spectral boundary integral methods for Helmholtz problems on smooth closed surfaces in R3, Numer. Math. 92 (2002) 289-323. [37] D. Colton, R. Kress, Inverse Coustic and Electromagnetic Scattering Theory, Applied Mathematical Sciences, Springer-Verlag, Heidelberg, 2ndEdition, 1998. [38] K.E. Atkinson, The Numerical Solution of Integral Equations of the Second Kind, Cambridge:Cambridge University Press, 1997.

22