SCW method for solving the fractional integro-differential equations with a weakly singular kernel

SCW method for solving the fractional integro-differential equations with a weakly singular kernel

Applied Mathematics and Computation 275 (2016) 72–80 Contents lists available at ScienceDirect Applied Mathematics and Computation journal homepage:...

410KB Sizes 0 Downloads 75 Views

Applied Mathematics and Computation 275 (2016) 72–80

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

SCW method for solving the fractional integro-differential equations with a weakly singular kernel Yanxin Wang a,∗, Li Zhu b a b

School of Science, Ningbo University of Technology, 315211 Ningbo, China School of Applied Mathematics, Xiamen University of Technology, 361024 Xiamen, China

a r t i c l e

i n f o

Keywords: Weakly singular integro-differential equations SCW Operational matrix Block pulse functions Fractional calculus

a b s t r a c t In this paper, based on the second Chebyshev wavelets (SCW) operational matrix of fractional order integration, a numerical method for solving a class of fractional integro-differential equations with a weakly singular kernel is proposed. By using the operational matrix, the fractional integro-differential equations with weakly singular kernel are transformed into a system of algebraic equations. The upper bound of the error of the second Chebyshev wavelets expansion is investigated. Finally, some numerical examples are shown to illustrate the efficiency and accuracy of the approach. © 2015 Elsevier Inc. All rights reserved.

1. Introduction In recent years, fractional calculus has attracted many researchers successfully in different disciplines of science and engineering. The advantage of this subject is that fractional derivative and integral are not a local property. Fractional differential equations and fractional integral equations are used to model a lot of practical problems, such as the field of electromagnetic waves, viscoelasticity, dielectric polarization and diffusion equations [1,2]. Fractional integro-differential equations with a weakly singular kernel have many applications in various areas. These equations appear in the radiative equilibrium [3], heat conduction problem [4], elasticity and fracture mechanics [5] and so on. There have been several numerical methods for the fractional integro-differential equations. For instance, Taylor expansion method [6], homotopy analysis method [7], hybrid collocation method [8], the variational iteration method [9], reproducing kernel method [10], pseudo spectral method [11] and wavelet methods [12–15]. However, only a few methods are proposed to solve the fractional integro-differential equations with a weakly singular kernel [16]. Wavelet analysis plays a prominent role in different areas of mathematics and engineering. Wavelets permit the accurate representation of a variety of functions and operators, and establish a connection with fast numerical algorithms [17]. The wavelets have been used in the solution of all kinds of equations. Using wavelet numerical method has several advantages: (a) The main advantage is that after discreting the coefficient matrix of algebraic equation is sparsity; (b) The wavelet method is computer oriented, thus solving higher order equation becomes a matter of dimension increasing; (c) The solution is a multi-resolution type; (d) The solution is convergence, even the size of increment may be large. Furthermore, SCW have been used to approximate the solution the integral equations of the second kind [18], fractional differential equations [19,20], fractional integro-differential equations [13,14], weakly singular Volterra integral equations [21] and the calculus of variations



Corresponding author. Tel.: +86 15258396709. E-mail addresses: [email protected] (Y. Wang), [email protected] (L. Zhu).

http://dx.doi.org/10.1016/j.amc.2015.11.057 0096-3003/© 2015 Elsevier Inc. All rights reserved.

Y. Wang, L. Zhu / Applied Mathematics and Computation 275 (2016) 72–80

73

question [22]. The main objective of the present paper is solving the fractional integro-differential equations with a weakly singular kernel base on SCW operational matrix. The structure of this paper is as follows: In Section 2, SCW are introduced, the convergence of SCW is given. The SCW operational matrix of fractional integration is introduced in Section 3. In Section 4, we summarize the process of solving the fractional integro-differential equations with a weakly singular kernel based on the SCW operational matrix method. The upper bound of the error of the SCW expansion is proposed in Section 5. Several examples are provided to clarify the approach in Section 6. Concluding remarks are given in the last section. 2. The SCW and their properties The SCW which defined on the interval [0, 1) have the following form [14,19]:

⎧ ⎨2 2k U˜ (2kt − 2n + 1 ), m ψnm (t ) = ⎩0,

n−1 n  t < k−1 , 2k−1 2 otherwise,

˜ m (t ) = where n = 1, . . . , 2k−1 and k is any positive integer, and U



(1)

π Um (t ), here the coefficient 2



2/π is used for orthonormality;

Um (t) is the second Chebyshev polynomials of degree m which respect to the weight function ω (t ) = on the interval [−1, 1] by the recurrence:



1 − t 2 . They are defined

U0 (t ) = 1, U1 (t ) = 2t, Um+1 (t ) = 2tUm (t ) − Um−1 (t ), m = 1, 2, . . . The weight function ω ˜ (t ) = ω (2t − 1 ) has to be dilated and translated as ωn (t ) = ω (2k t − 2n + 1 ). A function f(t) defined over [0,1], may be expressed in terms of the SCW as

f (t ) =

∞  

cnm ψnm (t ),

n=0 m∈Z

where the coefficients cnm are given by

cnm = ( f (t ), ψnm (t ))ωn =



1

0

ωn (t )ψnm (t ) f (t ) dt.

We can approximate the function f(t) by the truncated series k−1

f (t ) 

M−1 2  

cnm ψnm (t ) = C T (t ),

(2)

n=1 m=0

where the coefficient vector C and SCW function vector  (t) are given by

C = [c10 , c11 , . . . , c1(M−1) , c20 , . . . , c2(M−1) , . . . , c2k−1 0 , . . . , c2k−1 (M−1) ]T ,

(3)

(t ) = [ψ10 , ψ11 , . . . , ψ1(M−1) , ψ20 , . . . , ψ2(M−1) , . . . , ψ2k−1 0 , . . . , ψ2k−1 (M−1) ]T .

(4)

Taking the collocation points as following

ti =

2i − 1 , i = 1, 2, . . . , 2k−1 M. 2k M

(5)

We define the SCW matrix m ×m as

m ×m

 

3

1 2m − 1 =  ,  ,..., ,    2m

2m

2m

where m = 2k−1 M. Also, a function k(t, s ) ∈ L2ωn ([0, 1] × [0, 1] ) can be approximated as

k(t, s ) = (t )T K (s ),

(6)

where K is matrix with ki j = (ψi (t ), (k(t, s ), ψ j (s ))). We use the wavelet collocation method to determine the coefficients kij . Taking collocation points as Eq. (5), we obtain the matrix form of Eq. (6) 2k−1 M

× 2k−1 M

U = Tm ×m K m ×m , where K = [ki j ]m ×m , U = [k(ti , s j )]m ×m . We investigate the convergence of the SCW expansion in the following lemma.

74

Y. Wang, L. Zhu / Applied Mathematics and Computation 275 (2016) 72–80 

Lemma 2.1 [20]. A square-integrable function with bounded second-order derivative, say f (t ) defined on [0, 1], is with bounded  ˜ can be expanded as an infinite sum of the second Chebyshev wavelets, and the series converges second derivative, say | f (t )| ≤ M, uniformly to the function f(t), that is

f (t ) =

∞  ∞ 

cnm ψnm (t ),

n=0 m=0

where the coefficients cnm are given by cnm = ( f (t ), ψnm (t ))ωn and (·, · )ωn is the inner product of f(t) and ψ nm (t). 3. Operational matrix of the fractional integration 3.1. Fractional calculus Before we introduce the SCW operational matrix of the fractional integration, we first review some basic definitions of fractional calculus, which have been given in [23,24]. Definition 3.1. The Riemann–Liouville fractional integral operator Jα of order α is given by

⎧ ⎨ 1  t (t − τ )α−1 f (τ )dτ , α J f (t ) = (α ) 0 ⎩ f (t ),

α > 0, t > 0,

(7)

α = 0.

The properties of the operator Jα are given as follows: (i) (ii) (iii)

Jα Jβ f (t ) = Jα +β f (t ); Jβ Jα f (t ) = Jα +β f (t ); Jα t γ = ((γ + 1 )/(γ + α + 1 ))t α +γ .

Definition 3.2. The Caputo definition of fractional differential operator is given by

⎧ r d f (t ) ⎪ ⎨ dt r ,

Dα∗ f (t ) =

⎪ ⎩

α = r ∈ N, 

1 (r − α )

t 0

f (r ) ( τ ) dτ , (t − τ )α−r+1

(8)

0 ≤ r − 1 < α < r.

(r−α ) Dr f (t ), where Dr is the usual integer differential operThe Caputo fractional derivatives of order is also defined as Dα ∗ f (t ) = J ator of order r. The relation between the Riemann–Liouville operator and Caputo operator is given by the following expressions:

Dα∗ Jα f (t ) = f (t ), Jα Dα∗ f (t ) = f (t ) −

(9) r−1 

f ( k ) ( 0+ )

k=0

tk , t > 0. k!

(10)

3.2. SCW operational matrix of the fractional integration In this part, we simply introduce the SCW operational matrix of the fractional integration, more detailed introduction can be found in [14,19]. Now we give the definition of block pulse functions (BPFs): an m−set of BPFs on [0, 1) is defined as



bi (t ) =

1,

i/m  t < (i + 1 )/m,

0,

otherwise,

(11)

where i = 0, 1, 2, . . . , m − 1. The BPFs have disjointness and orthogonality as following:



bi (t )b j (t ) =

0,

i = j,

bi (t ),

i = j,

and

 0

1

 bi ( τ )b j ( τ )d τ =

0,

i = j,

1/m,

i = j,

The block pulse operational matrix of the fractional integration Fα has been given in [25],

Jα (Bm (t )) ≈ F α Bm (t ),

(12)

Y. Wang, L. Zhu / Applied Mathematics and Computation 275 (2016) 72–80

where



ξ1

1

⎢0 ⎢ ⎢0 ⎢ 1 1 α ⎢. F = α m (α + 2 ) ⎢ ⎢ .. ⎢ ⎣0 0

1

ξ2 ξ1

0

1

ξ3 ξ2 ξ1

.. .

..

..

0

···

0

1

ξ1

0

0

···

0

1

.

··· ···

ξm −1 ξm −2 ξm −3

.

..

.. .

···

.

75

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

(13)

and

ξk = (k + 1 )α+1 − 2kα+1 + (k − 1 )α+1 . Note that for α = 1, Fα is BPF’s operational matrix of integration. There is a relation between BPFs and SCW, namely

m (t ) = m ×m Bm (t ).

(14)

If Jα is the fractional integration operator of SCW, we can get

Jα (m (t )) ≈ Pmα ×m m (t ),

(15)

α where matrix Pm  ×m is called SCW operational matrix of the fractional integration [13,14]. Using Eqs. (12) and (14), we have

Jα (m (t )) ≈ Jα (m ×m Bm (t )) = m ×m Jα (Bm (t )) ≈ m ×m F α Bm (t ).

(16)

From Eqs. (15) and (16) we can get

Pmα ×m m (t ) = Pmα ×m m ×m Bm (t ) = m ×m F α Bm (t ). Then, the matrix P α

m ×m



m ×m

(17)

is given by

= m ×m F α −1 . m ×m

(18)

4. Solving fractional integro-differential equations with weakly singular kernel Consider the fractional integro-differential equations with a weakly singular kernel:

Dα∗ y(t ) = λ1



t 0

y (s ) ds + λ2 ( x − t )β



1 0

k(t, s )y(s )ds + f (t ),

(19)

subject to initial condition

y(0 ) = 0,

(20)

k(t, s) and f(t) are known and λ1 , λ2 are real constants. Here, 0 < α , β < 1 where y(t) is an unknown function, and Dα ∗ denotes the Caputo fractional order derivative of order α . The different values of λ1 , λ2 , α , β determine the type of the equation. In the case when α = 1, β = 0, this equation is offered as the linear integro-differential equation. Particularly, if α = 0 and λ2 = 0, we have the Abel integral equation: L2 -functions

y(t ) = λ1



0

t

y (s ) ds + f (t ). ( x − t )β

The treatment of Eq. (19) is complicated at 0 < β < 1, because the solutions of weakly singular fractional integro-differential equation are well known to have a weak singularity at t = 0, even when the inhomogeneous term f(t) is smooth. Now we approximate Dα ∗ y (t ), f(t) and k(t, s) in terms of SCW as follows

Dα∗ y(t ) ≈ C T (t ),

(21)

f (t ) ≈ F T (t ),

(22)

k(t, s ) ≈ (t )T K (s ),

(23)

and where K is 2k−1 M × 2k−1 M matrix with ki j = (ψi (x ), (k(x, t ), ψ j (t ))). From Eqs. (10), (15) and (21), we obtain

y(t ) = Jα Dα∗ y(t ) ≈ Jα (C T (t )) = C T Jα ((t )) = C T P α (t ). Then there is



t 0

y (s ) ds = C T P α (t − s )β



t 0

(s ) ds = (1 − β )C T P α P 1−β (t ). (t − s )β

(24)

(25)

76

Y. Wang, L. Zhu / Applied Mathematics and Computation 275 (2016) 72–80

1

From Eqs. (23), (24) and

 0

1

0

k(t, s )y(s )ds =

(s )(s )T ds = D (D is a m × m matrix), we have



1

0

(t )T K (s )(s )T Pα TCdt

= (t )T KDP α C T

= C T P α DT K T (t ).

(26)

By substituting Eqs. (21), (25) and (26) into Eq. (19), we obtain

C T (t ) = λ1 (1 − β )C T P α P 1−β (t ) + λ2C T P α DT K T (t ) + F T (t ).

(27)

Substituting Eq. (14) into Eq. (27), we have

C T m ×m Bm (t ) = λ1 (1 − β )C T P α P 1−β m ×m Bm (t ) + λ2C T P α DT K T m ×m Bm (t ) + F T m ×m Bm (t ).

(28)

Dispersing Eq. (28), we obtain

C T m ×m = λ1 (1 − β )C T P α P 1−β m ×m + λ2C T P α DT K T m ×m + F T m ×m .

(29)

By solving this system we can obtain the approximate solution of Eq. (19) according to Eq. (24). 5. Error analysis T In this section, in order to illustrate the effectiveness of Dα ∗ y (t ) ≈ C (t ) in Eq. (21), we have given the following theorem. α y(t ): Let Dα y ( t ) is the following approximation of D ∗ k,M ∗ k−1

Dα∗ yk,M (t ) =

M−1 2  

cnm ψnm (t ),

n=1 m=0

+∞

α then we have Dα ∗ y (t ) − D∗ yk,M (t ) =

n=2k−1 +1

+∞

m=M cnm ψnm (t ).

α α Theorem 5.1. Suppose that the function Dα ∗ yk,M (t ) obtained by using SCW are the approximation of D∗ y (t ), and D∗ y (t ) is with +2 y (t )| ≤ M ˜ ˜ for some constant M, then we have the following upper bound of error bounded second derivative, say |Dα ∗

Dα∗ y(t ) − Dα∗ yk,M (t ) E ≤ where y(t ) E = (

1 0



π M˜



23

+∞ 

n=2k−1 +1

+∞ 1  1 4 n5 m=M (m − 1 )

 12

,

1

y2 (t )ωn (t ) dt ) 2 .

Proof. The proof is similar to that of the theorem 1 in [20], so it is omitted here.  6. Numerical examples To show the efficiency and the accuracy of the proposed method, we consider the following four examples. Example 1. Consider the following fractional order integro-differential equation with weakly singular kernel [16]

D∗0.25 y(t ) =

1 2



t

0

1 y (s ) ds + 3 (t − s )1/2



1 0

(t − s )y(s )ds + f (t ),

(30)

with the condition y(0 ) = 0. The exact solution of this equation is y(t ) = t 2 + t 3 . In this problem

f (t ) =

(3 ) 1.75 (4 ) 2.75 t + t − (2.75 ) (3.75 )





3 π (3 )t 5/2 π (4 )t 7/2 7t − − + . 2(7/2 ) 2(9/2 ) 36 20

We apply SCW approach to solve Eq. (30) with various values of k. Table 1 shows the absolute errors obtained by SCW and CAS wavelets (m = 2k (2M + 1 )) respectively. The comparisons between the numerical solutions and the exact solution for M = 3 and some k are given in Fig. 1. Table 1 shows that the absolute errors become smaller and smaller with k increasing. From Table 1, we find that the SCW method can reach a higher degree of accuracy than the CAS wavelets method. From Fig. 1, we can see clearly that the numerical solution are more and more close to the exact solution when k increases. Example 2. Consider this equation [16]

D∗0.15 y(t ) =

1 4



t 0

1 y (s ) ds + 7 (t − s )1/2

 0

1

et+s y(s )ds + f (t ),

(31)

Y. Wang, L. Zhu / Applied Mathematics and Computation 275 (2016) 72–80

77

Table 1 The absolute errors of different k and M = 3 for Example 1. t

SCW k=2

CAS m = 6

SCW k=3

CAS m = 12

SCW k=4

CAS m = 24

0 1/6 2/6 3/6 4/6 5/6

8.3542e-03 1.2599e-03 9.3654e-03 2.2406e-02 1.9585e-02 3.2596e-02

3.0586e-02 4.4076e-02 3.8707e-02 1.5703e-02 2.8551e-02 9.8881e-02

1.0620e-03 1.1189e-03 1.8879e-03 4.7945e-03 5.9543e-03 8.0256e-03

1.4328e-02 2.2762e-02 1.9409e-02 6.4173e-03 1.8012e-02 5.6450e-02

1.4395e-04 2.2617e-04 5.9826e-04 1.1274e-03 1.4961e-03 2.2056e-03

6.3491e-03 1.1460e-02 9.6982e-03 2.9504e-03 9.6733e-03 2.9484e-02

k=2,M=3

k=3,M=3

1.8

k=4,M=3

2 Numerical solution Exact solution

1.6

2 Numerical solution Exact solution

1.8

1.4

Numerical solution Exact solution

1.8

1.6

1.6

1.4

1.4

1.2

1.2

1.2

y(t)

y(t)

y(t)

1 1

1

0.8 0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.6

0.4

0.2

0

0

0.5 t

1

0

0

0.5 t

0

1

0

0.5 t

1

Fig. 1. Comparisons of numerical solutions and exact solution of Example 1 for different k.

Table 2 The numerical solutions at different points for different k for Example 2. SCW k=2

CAS m = 6

SCW k=3

CAS m = 12

SCW k=4

CAS m = 24

Exact solution

t 0 1/6 2/6 3/6 4/6 5/6

−1.5509e-02 −1.4412e-01 −2.2290e-01 −2.4994e-01 −2.2128e-01 −1.3718e-01

6.1480e-03 −1.1115e-01 −1.8713e-01 −2.1160e-01 −1.8346e-01 −1.0222e-01

−7.2939e-03 −1.3954e-01 −2.2252e-01 −2.4997e-01 −2.2196e-01 −1.3843e-01

2.3540e-03 −1.2653e-01 −2.0680e-01 −2.3341e-01 −2.0590e-01 −1.2407e-01

−3.4910e-03 −1.3912e-01 −2.2230e-01 −2.4999e-01 −2.2216e-01 −1.3878e-01

9.0900e-04 −1.3345e-01 −2.1534e-01 −2.4266e-01 −2.1514e-01 −1.3269e-01

0 −1.3889e-01 −2.2222e-01 −2.5000e-01 −2.2222e-01 −1.3889e-01





3 ) 1.85 2 ) 0.85 (3 )t (2 )t such that y(0 ) = 0, f (t ) = (( t − (( t − π4( + π4( − e 7−3e . The exact solution is y(t ) = t 2 − t. 2.85 ) 1.85 ) 7/2 ) 5/2 ) Table 2 shows the numerical solutions by using the SCW and CAS wavelets, respectively. The comparisons of exact solution and numerical solutions for various k are shown in Fig. 2. From Table 2 and Fig. 2 we infer that the approximate solutions converge to the exact solution. 5/2

3/2

t+1

t

78

Y. Wang, L. Zhu / Applied Mathematics and Computation 275 (2016) 72–80

k=2,M=3

k=3,M=3

−0.06

k=4,M=3

0

0

Numerical solution Exact solution

−0.08

Numerical solution Exact solution

−0.1

Numerical solution Exact solution

−0.05

−0.05

−0.1

−0.1

−0.12

−0.18

y(t)

y(t)

−0.16

−0.15

−0.15

−0.2

−0.2

−0.2

−0.22

−0.24

−0.26

0

0.5 t

1

−0.25

0

0.5 t

1

−0.25

0

0.5 t

Fig. 2. Comparisons of numerical solutions and exact solution of Example 2 for different k.

k=5,M=2 3 Numerical solution of α=0.7 Numerical solution of α=0.8 Numerical solution of α=0.9 Numerical solution of α=1 Exact solution α=1

2.5

2

y(t)

y(t)

−0.14

1.5

1

0.5

0

0

0.2

0.4

0.6 t

Fig. 3. Numerical solutions and exact solution of α = 1.

0.8

1

1

Y. Wang, L. Zhu / Applied Mathematics and Computation 275 (2016) 72–80

79

Example 3. Consider the following equation

Dα∗ y(t ) =



t 0



y (s ) ds + (t − s )1/2



1 0

(t 2 + cos s )y(s )ds + f (t ),

(32)

(3 )t − t3 − 2 cos 1 + sin 1, with this supplementary condition y(0 ) = 0. The exact solution is y(t ) = t 2 where f (t ) = 2t − π2( 7/2 ) when α = 1. By setting k = 5 and M = 2, we obtain SCW solutions for various α . In Fig. 3, we show the SCW solutions and the exact solution for various values of α . It is evident from Fig. 3 that, as α is close to 1, the numerical solutions by SCW, converge to the exact solution, that is, the solution of the fractional integro-differential equation approaches the solution of the integer order integro-differential equation. 5/2

2

Let’s consider example with non-smooth and singular solutions.

k=4,M=2

k=5,M=2

6

k=6,M=2

8

14

Numerical solution Exact solution

Numerical solution Exact solution

Numerical solution Exact solution

5.5 7

12

6

10

5

8

5

4.5

3.5

u(t)

u(t)

y(t)

4

4

6

3

4

2

2

3

2.5

2

1.5

1

0

0.5 t

1

1

0

0.5 t

1

0

0

Fig. 4. Comparisons of numerical solutions and exact solution of Example 4 for different k.

Table 3 The absolute errors of different k and M = 2 for Example 4. t

k=4

k=5

k=6

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

6.7355e-02 9.4799e-02 1.1993e-01 1.1628e-01 1.0517e-01 9.9094e-02 8.9717e-02 7.9154e-02 6.8848e-02

9.7161e-02 5.6899e-02 9.2375e-02 9.7659e-02 9.0248e-02 8.3582e-02 7.4367e-02 6.5842e-02 5.8091e-02

9.8457e-02 4.1789e-02 6.9725e-02 7.3216e-02 6.9329e-02 6.3163e-02 5.6505e-02 5.0045e-02 4.4036e-02

0.5 t

1

80

Y. Wang, L. Zhu / Applied Mathematics and Computation 275 (2016) 72–80

Example 4. Consider the following equation

Dα∗ y(t ) = −



t 0

y (s ) ds − (t − s )1/2



1 0

tsy(s )ds + f (t ),

(33) √

−3/4

πt as the exact solution, with this supplementary condition y(0 ) = 0, where f (t ) = 23 t + ( + π . In this 1/4 ) case, there is a singularity at point t = 0. The solution around this point is not good (see Fig. 4 with M = 2 and some k). Due to this fact, the numerical results in Table 3 have been calculated over the interval [0.005: 1]. As it is observed, our method provides a reasonable estimate even in this case with singular solution.

which has y(t ) =

1 √ t

7. Conclusion In this paper, a numerical method is presented for numerical solutions of the fractional order integro-differential equation with a weakly singular kernel. Using the properties of SCW, we transform the initial problem into a linear algebraic system equations. By solving the linear system, numerical solutions are obtained. Moreover, The error analysis of SCW is proposed. Also, We have considered various types of solutions, with smooth, non-smooth and even singular behavior. In all the considered cases, the method provides reasonable estimates. The numerical results show that the approximations are in very good coincidence with the exact solution. The obtained results are compared with exact solution and with the solutions obtained by CAS wavelets. The present method is better than CAS wavelets method. The main character of our method is to combine the SCW method with the definition of Riemann–Liouville fractional integral operator dealing with the weakly singular integral, and the SCW basis functions are orthogonal and have small intervals of support, they can lead to a sparse matrix representation. Acknowledgments The authors are very thankful to the reviewers and the editor of this paper for their constructive comments and nice suggestions, which helped to improve the paper. This work was supported by the the Tian Yuan Special Funds of the National Natural Science Foundation of China (grant no. 11426192) and the Research Project of Education of Zhejiang Province (Y201431734). References [1] R.C. Koeller, Application of fractional calculus to the theory of viscoelasticity, J.Appl. Mech. 51 (1984) 299–307. [2] M. Ichise, Y. Nagayanagi, T. Kojima, An analog simulation of noninteger order transfer functions for analysis of electrode process, J. Electroanal. Chem. 33 (1971) 253–265. [3] P.K. Kythe, P. Puri, Computational Method for Linear Integral Equations, Birkhauser, Boston, 2002. [4] B.Q. Tang, X.F. Li, Solution of a class of Volterra integral equations with singular and weakly singular kernels, Appl. Math. Comput. 199 (2008) 406–413. [5] V.V. Zozulya, P.I. Gonzalez-Chi, Weakly singular, singular and hypersingular integrals in 3-D elasticity and fracture mechanics, J. Chin. Inst. Eng. 22 (1999) 763–775. [6] L. Huang, X.F. Li, Y.L. Zhao, X.Y. Duan, Approximate solution of fractional integro-differential equations by Taylor expansion method, Comput. Math. Appl. 62 (2011) 1127–1134. [7] X.D. Zhang, B. Tang, Y.N. He, Homotopy analysis method for higher-order fractional integro-differential equations, Comput. Math. Appl. 62 (2011) 3194–3203. [8] X.H. Ma, C.M. Huang, Numerical solution of fractional integro-differential equations by a hybrid collocation method, Appl. Math. Comput. 219 (2013) 6750– 6760. [9] K. Sayevand, Analytical treatment of Volterra integro-differential equations of fractional order, Appl. Math. Model. 39 (2015) 4330–4336. [10] W. Jiang, T. Tian, Numerical solution of nonlinear Volterra integro-differential equations of fractional order by the reproducing kernel method, Appl. Math. Model. 39 (2015) 4871–4876. [11] X.J. Tang, H.Y. Xu, Fractional pseudospectral integration matrices for solving fractional differential, integral, and integro-differential equations, Commun. Nonlinear Sci. Numer. Simulat. 30 (2016) 248–267. [12] H. Saeedi, M.M. Moghadam, N. Mollahasani, G.N. Chuev, A CAS wavelet method for solving nonlinear Fredholm integro-differential equations of fractional order, Commun. Nonlinear Sci. Numer. Simulat. 16 (2011) 1154–1163. [13] L. Zhu, Q.B. Fan, Solving fractional nonlinear Fredholm integro-differential equations by the second kind Chebyshev wavelet, Commun. Nonlinear Sci. Numer. Simul. 17 (6) (2012) 2333–2341. [14] L. Zhu, Q.B. Fan, Numerical solution of nonlinear fractional-order Volterra integro-differential equations by SCW, Commun. Nonlinear Sci. Numer. Simul. 18 (2013) 1203–1213. [15] Z. Meng, L. Wang, H. Li, W. Zhang, Legendre wavelets method for solving fractional integro-differential equations, Int. J. Comput. Math. 92 (6) (2015) 1275– 1291. [16] M.X. Yi, J. Huang, CAS wavelet method for solving the fractional integro-differential equation with a weakly singular kernel, Int. J. Comput. Math. 92 (8) (2015) 1715–1728. [17] G. Beylkin, R. Coifman, V. Rokhlin, Fast wavelet transform and numerical algorithms I, Commun. Pur. Appl. Math. 44 (1991) 141–183. [18] L. Zhu, Y.X. Wang, Q.B. Fan, Numerical computation method in solving integral equation by using the second Chebyshev wavelets, in: The 2011 International Conference on Scientific Computing, Las Vegas, USA, 2011, pp. 126–130. [19] Y.X. Wang, Q.B. Fan, The second kind Chebyshev wavelet method for solving fractional differential equation, Appl. Math. Comput. 218 (2012) 8592–8601. [20] F. Zhou, X. Xu, Numerical solution of the convection diffusion equations by the second kind Chebyshev wavelets, Appl. Math. Comput. 247 (2014) 353–367. [21] L. Zhu, Y.X. Wang, Numerical solutions of Volterra integral equation with weakly singular kernel using SCW method, Appl. Math. Comput. 260 (2015) 63–70. [22] L. Zhu, Y.X. Wang, SCW operational matrix of integration and its application in the calculus of variations, Int. J. Comput. Math. 90 (11) (2013) 2338–2352. [23] T. Kaczorek, Selected Problems of Fractional Systems Theory, Springer-Verlag, Berlin, 2011. [24] B.J. West, M. Bolognab, P. Grigolini, Physics of Fractal Operators, Springer, New York, 2003. [25] A. Kilicman, Z.A.A.A. Zhour, Kronecker operational matrices for fractional calculus and some applications, Appl. Math. Comput. 187 (2007) 250–265.