Shifted Jacobi–Gauss-collocation with convergence analysis for fractional integro-differential equations

Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359

Contents lists available at ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

Research paper

Shifted Jacobi–Gauss-collocation with convergence analysis for fractional integro-differential equations E.H. Doha a, M.A. Abdelkawy b,c, A.Z.M. Amin c, António M. Lopes d,∗ a

Department of Mathematics, Faculty of Science, Cairo University, Giza, Egypt Department of Mathematics and Statistics, College of Science, Al-Imam Mohammad Ibn Saud Islamic University (IMSIU), Riyadh, Saudi Arabia c Department of Mathematics, Faculty of Science, Beni-Suef University, Beni-Suef, Egypt d UISPA–LAETA/INEGI, Faculty of Engineering, University of Porto, Porto, Portugal b

a r t i c l e

i n f o

Article history: Received 18 October 2018 Revised 18 December 2018 Accepted 11 January 2019 Available online 12 January 2019 Keywords: Fractional integro-differential equation Spectral collocation method Jacobi–Gauss quadrature Riemann–Liouville derivative

a b s t r a c t A new shifted Jacobi–Gauss-collocation (SJ-G-C) algorithm is presented for solving numerically several classes of fractional integro-differential equations (FI-DEs), namely Volterra, Fredholm and systems of Volterra FI-DEs, subject to initial and nonlocal boundary conditions. The new SJ-G-C method is also extended for calculating the solution of mixed Volterra–Fredholm FI-DEs. The shifted Jacobi–Gauss points are adopted for collocation nodes and the FI-DEs are reduced to systems of algebraic equations. Error analysis is performed and several numerical examples are given for illustrating the advantages of the new algorithm. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Fractional differential equations (FDEs) [1–5] are powerful tools for modeling phenomena in mathematical chemistry [6,7], biology [8], viscoelasticity [3], physics [4], and other areas [9–14]. The increasing applicability of FDEs has required efficient algorithms for calculating their solutions. However, since most FDEs cannot be solved analytically, numerical methods have been developed [10–12]. Despite the intense research that has been carried out in this topic, the problem is still challenging. Fractional integro-differential equations (FI-DEs) are widely used in science and engineering. For example, the activity of interacting inhibitory and excitatory neurons is well modeled by means of FI-DEs. Detailed numerical methods for solving one-dimensional FI-DEs were presented in [15–24], while several techniques for general FI-DEs were developed by many authors. For example, Nazari and Shahmorad [25] introduced the fractional differential transform method for FI-DEs with nonlocal boundary conditions. Jiang and Tian [26] improved the reproducing kernel scheme for nonlinear Volterra FI-DEs. Saeedi and Moghadam [27] used the CAS wavelets technique for solving nonlinear Volterra FI-DEs of arbitrary order. Susahab et al. [28] applied quadrature rules to a class of nonlinear FI-DEs of the Hammerstein type. Zhu and Fan [29] developed the second kind Chebyshev wavelet for solving nonlinear Fredholm FI-DEs. Other authors [30–35] developed and employed different numerical techniques for solving FI-DEs.

∗

Corresponding author. E-mail address: [email protected] (A.M. Lopes).

https://doi.org/10.1016/j.cnsns.2019.01.005 1007-5704/© 2019 Elsevier B.V. All rights reserved.

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359

343

Numerical methods are divided into local and global techniques. The finite-difference and finite-element methods are classified as local techniques, whilst the spectral method is global. In practice, finite-element methods are particularly well suited to problems in complex geometries, whereas spectral methods can provide superior accuracy, at the expense of domain flexibility. We emphasize that there are many numerical approaches, such as hp finite-elements and spectral-elements, which combine advantages of both the global and local methods. However in this paper, we shall restrict our attentions to the global spectral methods. These can be thought as a development of the so-called method of weighted residuals. Recently, spectral methods were recognized as efficient numerical schemes for solving FI-DEs [32,33]. Spectral methods are characterized by having faster convergence rates and better accuracy than the local methods. According to this method, the solution of FI-DEs is expressed in terms of a finite series of known functions, which are global in the sense that they are defined over the entire domain and are called trial/basis functions. After substituting this series in the FI-DEs, an inner product of the resulting equation with the so-called test functions is formed, which is used in order to guarantee that the equation is satisfied as closely as possible by the truncated series expansion. This is accomplished by minimizing the error in the differential equation produced by using the truncated series expansion instead of the exact solution, with respect to a suitable inner product. Regarding the methodology used, the spectral methods are divided into four categories, namely collocation [36–40], tau [41,42], Galerkin [43] and Petrov-Galerkin [44] methods. In this paper, we propose an accurate numerical algorithm for calculating the solutions of different classes of FI-DEs with initial and nonlocal boundary conditions. Using the shifted Jacobi–Gauss collocation (SJ-G-C) method with the Riemann– Liouville (R-L) fractional derivative of the shifted Jacobi polynomials, we reduce the FI-DEs to systems of algebraic equations. The solution of such equations is approximated by means of a finite expansion of shifted Jacobi polynomials for independent variables (for more details see Canuto et al. [45]). Then we evaluate the residuals of the mentioned problem at the shifted Jacobi–Gauss quadrature points. Substituting these approximations in the FI-DEs leads to a system of algebraic equations. This system may be solved numerically using the Newton’s iterative algorithm. This scheme is one of the most suitable methods for solving systems of algebraic equations. Indeed, with the freedom to select the shifted Jacobi indexes σ > −1, ρ > −1, the method can be calibrated for a wide variety of problems. Moreover, we develop and analyze spectral collocation methods based on Jacobi polynomials with general parameters σ and ρ . The main advantage of the proposed algorithm is that the Chebyshev, Legendre and ultraspherical collocation methods can be obtained as special cases from our method. Furthermore, an error analysis of the new method is developed and the results are discussed. The paper is organized as follow. Section 2 introduces the tools of fractional calculus and the shifted Jacobi polynomials. Section 3 applies the new SJ-G-C to one-dimensional linear Volterra FI-DEs subject to nonlocal conditions, and to nonlinear Volterra FI-DEs subject to initial conditions. Sections 4–6 extend the SJ-G-C method to solve one-dimensional nonlinear Fredholm FI-DEs, systems of Volterra FI-DEs, and mixed Volterra–Fredholm FI-DEs, respectively. Section 7 presents some useful lemmas and error analysis. Section 8 solves some numerical examples and, finally, Section 9 draws the main conclusions. 2. Mathematical preliminaries 2.1. Fractional calculus The fractional integral and derivative of order ν > 0 can be expressed by means different formulas. Often we use the R-L definitions. Definition 2.1. The R-L fractional integral of order ν > 0, Jν , is given by

1

J ν  (z ) =

(ν ) 0 J  ( z ) =  ( z ),

where

(ν ) =



∞



z 0

(z − ζ )ν −1 (ζ )dζ ,

ν > 0, z > 0, (2.1)

zν −1 e−z dz.

0

The operator Jν satisfies

J ν J μ  ( z ) = J ν +μ  ( z ),

J ν J μ  ( z ) = J μ J ν  ( z ), (ρ + 1 ) J ν zρ = z ρ +ν . (ρ + 1 + ν )

(2.2)

Definition 2.2. The R-L fractional derivative of order ν > 0, Dν , is given by

Dν  ( z ) =

1 dm (m − ν ) dzm



z 0

where m is the ceiling function of ν .

 (z − ζ )m−ν −1 (ζ )dζ ,

m − 1 < ν ≤ m, z > 0,

(2.3)

344

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359

2.2. Properties of the shifted Jacobi polynomials We recall that (σ ,ρ )

(σ ,ρ )

Pk+1 (z ) = (ak (σ ,ρ )

P0

(σ ,ρ )

Pk

( z ) = 1,

(σ ,ρ )

σ ,ρ ) )Pk(σ ,ρ ) (z ) − ck(σ ,ρ ) Pk(−1 ( z ), k ≥ 1, 1 1 ( z ) = 2 ( σ + ρ + 2 )z + 2 ( σ − ρ ),

z − bk

(σ ,ρ )

P1

(−z ) = (−1 )k Pk(σ ,ρ ) (z ),

(σ ,ρ )

Pk

(−1 ) =

(−1 )k (k+ρ +1 ) , k!(ρ +1 )

(2.4)

where σ , ρ > −1, z ∈ [−1, 1] and (σ ,ρ )

ρ +1 )(2k+σ +ρ +2 ) = (2k+2σ(k++1 , )(k+σ +ρ +1 )

(σ ,ρ )

(ρ −σ )(2k+σ +ρ +1 ) = 2(k+1 )(k+σ +ρ +1 )(2k+σ +ρ ) ,

(σ ,ρ )

k+σ )(k+ρ )(2k+σ +ρ +2 ) = (k(+1 )(k+σ +ρ +1 )(2k+σ +ρ ) .

ak

2

bk ck

2

( σ ,ρ )

Furthermore, the rth derivative of P j (σ ,ρ )

Dr P j

(z ), is given by

( j+σ +ρ +q+1 ) (σ +r,ρ +r ) P j−r ( z ), 2r ( j+σ +ρ +1 )

(z ) =

(2.5) ( σ ,ρ )

where r is an integer. For the shifted Jacobi polynomial PL,k (σ ,ρ )

PL,k

(z )

= =

k 

(z ) = Pk(σ ,ρ ) ( 2Lz − 1 ), L > 0 we have

k− j (k+ρ +1 )( j+k+σ +ρ +1 ) −1 zj ( j+ρ +1 )(k+σ +ρ +1 )(k− j )! j!L j j=0 k (k+σ +1 )(k+ j+σ +ρ +1 ) j j=0 j!(k− j )!( j+σ +1 )(k+σ +ρ +1 )L j z − L .

(

)

(

(2.6)

)

Thereby, we deduce that (σ ,ρ )

PL,k

k (k+ρ +1 ) = (−1 ) ( ρ +1 ) k! ,

(0 ) (L )

(σ ,ρ )

PL,k

(2.7)

(k+σ +1 ) = ( σ +1 ) k! ,

(σ ,ρ )

(0 ) =

(−1 )k−r (k + ρ + 1 )(k + σ + ρ + 1 )r , Lr (k − r + 1 )(r + ρ + 1 )

(2.8)

(σ ,ρ )

(L ) =

(k + σ + 1 )(k + σ + ρ + 1 )r , Lr (k − r + 1 )(r + σ + 1 )

(2.9)

(σ ,ρ )

(z ) =

(r + k + σ + ρ + 1 ) (σ +r,ρ +r ) P ( z ). Lr (k + σ + ρ + 1 ) L,k−r

Dr PL,k Dr PL,k Dr PL,k

( σ ,ρ )

Taking wL

(2.10)

(z ) = (L − z )σ zρ , we define the norm and inner product related to the weighted space L2 (σ ,ρ ) [0, L] as wL

(u, v )w(σ ,ρ ) = L



L 0

(σ ,ρ )

 ( z ) v ( z ) wL

(z ) dz,

1

vw(σ ,ρ ) = (v, v )w2 (σ ,ρ ) . L

(2.11)

L

A complete L2 (σ ,ρ ) [0, L]-orthogonal system is a set of shifted Jacobi polynomials, where wL

(σ ,ρ ) 2 PL,k w(σ ,ρ ) =

 L σ +ρ +1

L

2

( σ ,ρ )

( σ ,ρ )

(σ ,ρ )

L (σ ,ρ ) (z + 1 ), 2 N, j

(σ ,ρ )

hk

(σ ,ρ )

= hL,k

.

(2.12)

Using zN, j and N, j , 0  j  N, as the nodes and Christoffel numbers of the standard Jacobi–Gauss interpolation in the interval [−1, 1]. In the interval [0, L] the corresponding nodes and Christoffel numbers of the shifted Jacobi–Gauss interpolation are

zL,N, j =

L 2

(σ ,ρ ) (σ ,ρ ) L,N, = ( )σ +ρ +1 N, j , 0  j  N. j

For any positive integer N, φ ∈ S2N+1 [0, L] and, by means of Jacobi–Gauss quadrature property, we obtain

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359



L 0

 L σ +ρ +1 

(L − z )σ zρ φ (z )dz =

2

1

−1

 L σ +ρ +1  N

=

2

N,(σj,ρ ) φ

L

j=0

2

L 2

 (z + 1 ) dz

(σ ,ρ ) (zN, + 1) j



  (σ ,ρ ) (σ ,ρ ) L,N, φ z . j L,N, j

N 

=

( 1 − z )σ ( 1 + z )ρ φ

345

(2.13)

j=0

3. Volterra FI-DEs 3.1. Linear Volterra FI-DEs with nonlocal boundary condition The SJ-G-C method is applied to numerically solve the linear Volterra FI-DE with nonlocal conditions

Dν  ( z ) = f ( z ) +



z

0

 (0 ) + γ  (1 ) + λ

k(z, ζ )(ζ )dζ ,



b

0 < ν < 1,

φ ( ζ ) ( ζ )d ζ = d1 ,

a

(3.1)

(3.2)

where f(z), φ (ζ ) and k(z, ζ ) are given functions, γ and λ are constants, and (z) is an unknown function. The solution of Eq. (3.1) is approximated by

N ( z ) =

N 

(σ ,ρ )

a j PL, j

( z ),

(3.3)

j=0

and the fractional derivative of N (z) is estimated as

Dν N ( z ) =

N  j=0



Given the R-L derivative

1 ∂ (1 − ν )

Dν z k =

(σ ,ρ )

a j Dν PL, j



(z ) .

(3.4)

 χk dχ , ( z − χ )ν

z

0

zk−ν (1 + k ) = , (1 + k − ν )

0 < ν < 1,

(3.5)

then (σ ,ρ )

Dν PL, j

(z ) = L,(σj ,ρ ) (z ) =

j  (−1 ) j−k ( j + ρ + 1 ) ( j + k + σ + ρ + 1 ) ν k D z (k + ρ + 1 ) ( j + σ + ρ + 1 ) ( j − k )! k! Lk k=0

=

j 

(−1 )i−k ( j + ρ + 1 ) ( j + k + σ + ρ + 1 ) z k −ν . ( k + ρ + 1 )( j + σ + ρ + 1 )( j − k )!(k − ν + 1 ) Lk k= ν

(3.6)

Accordingly,

Dν N ( z ) =

N  j=0



(σ ,ρ )

a j Dν PL, j

N

 (z ) = a j L,(σj ,ρ ) (z ).

(3.7)

j=0

Using the above relations we can rewrite Eq. (3.1) as N  j=0

(σ ,ρ )

a j L, j

(z ) = f (z ) +

 0

z



k(z, ζ )

N 

(σ ,ρ )

a j PL, j

 (ζ ) d ζ .

(3.8)

j=0

In the new SJ-G-C method the residual of (3.8) is set to zero at N of the shifted Jacobi–Gauss points. Employing (3.3)–(3.8), then we write (3.1) in the form

346

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359 N 

(σ ,ρ )

a j L, j

(α ,ρ ) (σ ,ρ ) (zL,N,i ) = f (zL,N,i )+



(σ ,ρ )

zL,N,i



0

j=0

(σ ,ρ )

k(zL,N,i , ζ )

N 

(σ ,ρ )

a j PL, j

 (ζ ) dζ , i = 1, . . . , N.

(3.9)

j=0

We may rearrange the above equation, yielding N 

(σ ,ρ )

a j L, j

(σ ,ρ ) (zL,N,i )−



j=0

(σ ,ρ )

zL,N,i 0

(σ ,ρ )

(σ ,ρ )

k(zL,N,i , ζ )PL, j

(σ ,ρ ) (ζ )dζ ) = f (zL,N,i ), i = 1, . . . , N.

(3.10)

Combining Eqs. (3.2) and (3.3), we obtain N 

(σ ,ρ )

a j PL, j

N 

(0 ) + γ

j=0

(σ ,ρ )

a j PL, j

(1 ) + λ

j=0

N 

aj



b

a

j=0

φ (ζ )PL,(σj ,ρ ) (ζ ) = d1 .

(3.11)

Finally, a linear system of (N + 1 ) algebraic equations is generated and N (z) determined.

3.2. Non-linear Volterra FI-DEs with initial conditions We extend the SJ-G-C algorithm to the nonlinear Volterra FI-DE

Dν  ( z ) =



z

0

k(z, ζ )((ζ )) p dζ + f (z ),

1 < ν < 2,

(3.12)

with the initial conditions

u ( m ) ( 0 ) = dm

m = 0, 1.

(3.13)

The spectral solution of Eq. (3.12) is approximated by

N ( z ) =

N 

(σ ,ρ )

a j PL, j

( z ).

(3.14)

j=0

Similar steps to those followed in the previous subsection allow us to write the nonlinear Volterra FI-DE in the form N 

(σ ,ρ )

a j L, j

(z ) =

j=0

 0

z



k(z, ζ )

N 

(σ ,ρ )

a j PL, j

(ζ )

p

d ζ + f ( z ).

(3.15)

j=0

As a result of the above-mentioned relation, we get N − 1 algebraic nonlinear equations in the following form, N 

(σ ,ρ )

a j L, j



(σ ,ρ ) (σ ,ρ ) (zL,N,i ) = f (zL,N,i )+

j=0

(σ ,ρ )

zL,N,i 0



(σ ,ρ )

k(zL,N,i , ζ )

N 

(σ ,ρ )

a j PL, j

(ζ )

p

d ζ , i = 1, . . . , N − 1.

(3.16)

j=0

Combining Eqs. (3.13) and (3.14), we obtain N  j=0

(σ ,ρ )

a j Dm PL, j

( 0 ) = dm ,

m = 0, 1.

(3.17)

If we use (2.7) and (2.8), then Eq. (3.17) can be rewritten as N  j=0

( j+ρ +1 ) (−1 ) j ( ρ +1 ) j! a j = d0 ,

N  (−1 ) j−1 ( j+ρ +1 )( j+σ +ρ +1 ) a j = d1 . L( j−1 )!(ρ +2 )

(3.18)

(3.19)

j=0

The Eqs. (3.16), (3.18) and (3.19) are equivalent to a system of (N + 1 ) algebraic nonlinear equations in the unknowns a j , i = 0, . . . , N

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359 N 

(−1 ) j

j=0

( j + ρ + 1 ) a = d0 (ρ + 1 ) j! j

N  (−1 ) j−1 ( j + ρ + 1 )( j + σ + ρ + 1 ) a j = d1 L( j − 1 )!(ρ + 2 ) j=0

N 

347

(σ ,ρ )

a j L, j

(σ ,ρ )

(σ ,ρ )

(zL,N,i ) = f (zL,N,i ) +



j=0

(σ ,ρ )

zL,N,i



(σ ,ρ )

k(zL,N,i , ζ )

0

N 

(σ ,ρ )

a j PL, j

(ζ )

p

d ζ , i = 1, . . . , N − 1.

(3.20)

j=0

Using the Newton iterative scheme the system may be solved and the coefficients aj determined. Therefore, we compute the approximate solution N (z) at any value of z in the given domain. 4. Non-linear Fredholm FI-DEs with initial conditions The SJ-G-C method is applied to numerically solve the Fredholm FI-DE with initial conditions

Dν  ( z ) = f ( z ) +



L

0

k(z, ζ )((ζ )) p dζ ,

1 < ν < 2,

(4.1)

subject to

u ( m ) ( 0 ) = dm

m = 0, 1.

(4.2)

The solution of Eq. (4.1) is approximated by

N ( z ) =

N 

(σ ,ρ )

a j PL, j

( z ).

(4.3)

j=0

Based on the results in the last subsections, we obtain N 

(σ ,ρ )

a j L, j

(z ) = f (z ) +



L

0

j=0



k(z, ζ )

N 

(σ ,ρ )

a j PL, j

 (ζ ) d ζ ,

(4.4)

j=0

thus, we get N 



(σ ,ρ )

a j L, j

(σ ,ρ ) (zL,N,i )−



j=0

L 0

(σ ,ρ )

(σ ,ρ )

k(zL,N,i , t )PL, j

 (σ ,ρ ) (ζ )dζ = f (zL,N,i )

i = 1, . . . , N − 1.

(4.5)

Merging Eqs. (4.3) and (4.2), yields N 

(σ ,ρ )

a j Dm PL, j

( 0 ) = dm ,

m = 0, 1.

(4.6)

j=0

On the other hand, we may write N 

(−1 ) j

j=0

( j + ρ + 1 ) a = d0 , (ρ + 1 ) j! j

(4.7)

N  (−1 ) j−1 ( j + ρ + 1 )( j + σ + ρ + 1 ) a j = d1 . L( j − 1 )!(ρ + 2 )

(4.8)

j=0

The Eqs. (4.5), (4.7) and (4.8) are equivalent to a system of (N + 1 ) algebraic equations in the unknowns a j , i = 0, . . . , N N 

(−1 ) j

j=0

( j + ρ + 1 ) a = d0 , (ρ + 1 ) j! j

N  (−1 ) j−1 ( j + ρ + 1 )( j + σ + ρ + 1 ) a j = d1 , L( j − 1 )!(ρ + 2 ) j=0

N  j=0



(σ ,ρ )

a j L, j

(σ ,ρ ) (zL,N,i )−

 0

L

(σ ,ρ )

(σ ,ρ )

k(zL,N,i , ζ )PL, j

 (σ ,ρ ) (ζ )dζ = f (zL,N,i )

Finally, the system can be easily solved and N (z) computed.

i = 1, . . . , N − 1.

(4.9)

348

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359

5. System of Volterra FI-DEs with initial conditions In this Section, we extend the SJ-G-C technique to solve the system of Volterra FI-DEs

⎧  z ⎪ ⎨Dν (z ) = v(z ) + f (z ) + k1 (z, ζ )[(ζ ) + v(ζ )]dζ , ⎪ ⎩Dν v(z ) = (z ) + g(z ) +



0 z

0

0 < ν < 1, (5.1)

k2 (z, ζ )[(ζ ) + v(ζ )]dζ ,

subject to the conditions

 ( 0 ) = d1 ,

v ( 0 ) = d2 .

(5.2)

Here, we approximate (z) and v(z) using shifted Jacobi polynomials

N ( z ) =

N 

(σ ,ρ )

a j PL, j

( z ),

vN ( z ) =

N 

j=0

(σ ,ρ )

b j PL, j

( z ).

(5.3)

j=0

Using (5.3), we deduce that

⎧ N N ⎪  (σ ,ρ ) ⎪ ν (σ ,ρ ) ⎪ ⎪ ⎨ j=0 a j D PL, j (z ) = j=0 b j PL, j (z )+ ⎪ N N   ⎪ (σ ,ρ ) ⎪ ν (σ ,ρ ) ⎪ ⎩ b j D PL, j (z ) = a j PL, j (z )+ j=0

f (z ) +

 0



g( z ) +



z

z

k1 (z, ζ )





N 

(σ ,ρ )

a j PL, j

(ζ ) +

j=0

N 

k2 (z, ζ )

0

j=0



N 

 (σ ,ρ )

b j PL, j

(ζ )

(σ ,ρ )

a j PL, j

(ζ ) +

j=0

N 

dζ ,



j=0

(σ ,ρ )

b j PL, j

(ζ )

(5.4)

dζ ,

j=0

accordingly, we obtain

⎧ N N ⎪   ⎪ ⎪ a j L,(σj ,ρ ) (z ) = b j PL,(σj ,ρ ) (z )+ ⎪ ⎨ j=0 j=0 ⎪ N N   ⎪ (σ ,ρ ) (σ ,ρ ) ⎪ ⎪ ⎩ b j L, j (z ) = a j PL, j (z )+ j=0

f (z ) + g( z ) +

j=0





z 0



z

 k1 (z, ζ )



 k2 (z, ζ )

0

N 

(σ ,ρ )

a j PL, j

(ζ ) +

j=0

N 



N 

(σ ,ρ )

b j PL, j

(ζ )

(σ ,ρ )

a j PL, j

(ζ ) +

j=0

N 

dζ



j=0

(σ ,ρ )

b j PL, j

(ζ )

(5.5)

dζ .

j=0

In the proposed SJ-G-C method the residual of (5.5) is set to zero at N − 1 of the shifted Jacobi–Gauss points

⎧ N ⎪  ⎪ (σ ,ρ ) (σ ,ρ ) ⎪ a j L, j (zL,N,i ) = ⎪ ⎪ ⎪ j=0 ⎪ ⎪ ⎨ ⎪ N  ⎪ (σ ,ρ ) (σ ,ρ ) ⎪ b j L, j (zL,N,i ) = ⎪ ⎪ ⎪ j=0 ⎪ ⎪ ⎩

N  j=0

(σ ,ρ )

b j PL, j

(σ ,ρ )

(zL,N,i ) +



(σ ,ρ )

zL,N,i

N  j=0

(σ ,ρ )

a j PL, j

(σ ,ρ )

(σ ,ρ )

(zL,N,i ) +



 (σ ,ρ )

k1 (zL,N,i , ζ )

0

N 

(σ ,ρ )

a j PL, j

(ζ ) +

j=0

(σ ,ρ )

+ f (zL,N,i )



(σ ,ρ )

zL,N,i 0



 (σ ,ρ )

k2 (zL,N,i , ζ )

N  j=0

N 

(σ ,ρ )

b j PL, j

(ζ )

dζ

j=0

(σ ,ρ )

a j PL, j

(ζ ) +

N 

(σ ,ρ )

b j PL, j

(ζ )

(5.6)

dζ

j=0

+g(zL,N,i ), i = 1, . . . , N.

Using the conditions (5.2), we obtain N 

(−1 ) j

( j + ρ + 1 ) a = d1 , (ρ + 1 ) j! j

(5.7)

(−1 ) j

( j + ρ + 1 ) b = d2 . (ρ + 1 ) j! j

(5.8)

j=0 N  j=0

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359

349

Finally, from Eqs. (5.6)–(5.8), we get a system of algebraic equations which can be solved for the unknown coefficients N  j=0 N  j=0 N  j=0

( j+ρ +1 ) (−1 ) j ( ρ +1 ) j!

a j = d1 ,

( j+ρ +1 ) (−1 ) j ( ρ +1 ) j!

b j = d2 ,

(σ ,ρ )

a j L, j

N 

(σ ,ρ ) (zL,N,i )=

(σ ,ρ )

j=0

b j PL, j

(σ ,ρ ) (zL,N,i )+



(σ ,ρ )

zL,N,i



0

(σ ,ρ )

k1 (zL,N,i , t )

 N

+ f (zL,N,i ) j=0

(σ ,ρ )

b j L, j

N 

(σ ,ρ ) (zL,N,i )=

(σ ,ρ )

j=0

a j PL, j

(σ ,ρ ) (zL,N,i )+



(ζ ) +

N 

j=0

(σ ,ρ )

N 

(σ ,ρ )

a j PL, j

(σ ,ρ )

zL,N,i



0

(σ ,ρ )

k2 (zL,N,i , ζ )

 N

(ζ )



dζ

j=0

(σ ,ρ )

a j PL, j

(ζ ) +

j=0

(σ ,ρ )

(σ ,ρ )

b j PL, j

N 

(σ ,ρ )

b j PL, j

(ζ )



(5.9)

dζ

j=0

+g(zL,N,i ), i = 1, . . . , N. 6. Mixed Volterra–Fredholm FI-DEs with nonlocal boundary conditions We present the SJ-G-C method to numerically solve the linear fractional mixed Volterra–Fredholm FI-DE

Dν  ( z ) = f ( z ) +



z

0

k(z, ζ )(ζ )dζ +



L 0

k(z, ζ )(ζ )dζ ,

(6.1)

with the nonlocal boundary conditions

 (0 ) + γ  (1 ) + λ



b a

φ ( ζ ) ( ζ )d ζ = d1 .

(6.2)

Based on the results presented in the previous subsections, we obtain the following system of algebraic equations N 



(σ ,ρ )

a j PL, j

 b

 (0 ) + γ PL,(σj ,ρ ) (1 ) + λ φ (ζ )PL,(σj ,ρ ) (ζ ) = d1 , a

j=0 N 

  (σ ,ρ ) (σ ,ρ ) a j L, j (zL,N,i ) −

(σ ,ρ )

zL,N,i

0

j=0

(σ ,ρ )

(σ ,ρ )

k(zL,N,i , ζ )PL, j

(ζ )dζ −



L 0

(σ ,ρ )

(σ ,ρ )

k(zL,N,i , ζ )PL, j

 (σ ,ρ ) (ζ )dζ ) = f (zL,N,i ,

i = 1, . . . , N.

(6.3)

After the coefficients aj are determined, it is straightforward to compute the approximate solution N (z) at any value of z in the given domain. 7. Lemmas and error analysis In this Section useful lemmas and error analysis of the J-G-C algorithm presented in Section 3.1. 7.1. Lemmas Definition 7.1. Let PN : L2 (I) → ZN be the L2 orthogonal projection defined by (PN  − , v ) = 0, ∀v ∈ ZN , Definition 7.2. For a nonnegative integer m, define the weighted Hilbert space [45,47]

Hwm(σ ,ρ ) (−1, 1 ) = {υ : ∂zi υ ∈ L2w(σ ,ρ ) (−1, 1 ), 0 ≤ i ≤ m}. i where ∂zi υ (z ) = ∂ υ (iz ) , related to the semi-norm and norm

∂z

|υ|m,w(σ ,ρ ) = ∂zm υ w(σ ,ρ ) ,   υ m,w(σ ,ρ ) =

m 

12 ∂ υ i z

2 w(σ ,ρ )

.

i=0

Lemma 7.1. Assume that  ∈ H m (I ),

I ≡ (−1, 1 ). The interpolation of u



( σ ,ρ )

IN



 at any of the Jacobi–Gauss points (Gauss or

Gauss–Radau or Gauss–Lobatto points) satisfies [45]

  − IN(σ ,ρ )  L2 (σ ,ρ ) (I ) ≤ CN−m |  |Hm,N w

w (σ ,ρ )

(I ) ,

(7.1)

350

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359

  − IN(σ ,ρ )  L∞ (I ) ≤ CN 2 −m |  |Hm,N 1

.

(7.2)

( σ ,ρ )

 denote the interpolation of  based on (N + 1 ) degree Jacobi–Gauss points

w (σ ,ρ )

(I )

Lemma 7.2. Assume that u ∈ H m(σ ,ρ ) (I ) and IN w

corresponding to the weight function



(σ ,ρ ) σ Dσ  − I D  N

w(σ , ρ ) (z)

L2 (σ ,ρ ) (I ) ≤ CN

−m

with −1 < (σ , ρ ), then (see [45, 47])

| Dσ  |Hm,N

w (σ ,ρ )

w

(I )

.

(7.3)

7.2. Error analysis We carry out error analysis for the SJ-G-C method introduced in Section 3.1 and demonstrate its exponential rate of convergence, provided that the source and kernel functions are sufficiently smooth. In order to do that, some properties of Banach algebra and Sobolev inequality are considered. σ ,ρ

Theorem 7.3. Let (z) be the exact solution of the Volterra FI-DE (3.1) and assume that IN (z ) is the spectral approximation defined by Eq. (3.3), therefore we have

  e(z ) L2 (σ ,ρ ) (I ) ≤ N−m C1 | Dσ  |Hm,N

2 (I ) +C2 λN w (σ ,ρ )

w

1

|  |Hm,N

(I ) +C5 | k (z, ζ ) |H m,N (I )   L2w(σ ,ρ ) (I ) w (σ ,ρ ) w (σ ,ρ )



(7.4)

Proof. Using the exact solution (z), the Volterra FI-DE in (3.1) may be written as

Dσ  ( z ) − f ( z ) −



z

0

(k(z, ζ ))(ζ )dζ = 0,

(7.5)

while using the approximate solution, we have (σ ,ρ ) σ

IN

D  (z ) − f (z ) −



z

0

(σ ,ρ )

(σ ,ρ )

IN,N (k(z, ζ ))IN

 ( ζ )d ζ = 0,

(7.6)

subtracting (7.5) from (7.6), we get (σ ,ρ ) σ

e(z ) = (IN

D (z ) − Dσ (z )) +



z 0

(k(z, ζ ))(ζ )dζ −



z 0

(σ ,ρ )

(σ ,ρ )

IN,N (k(z, ζ ))IN

 ( ζ )d ζ ,

 z  z (σ ,ρ ) D (z ) − Dσ (z )) + eN,N (k(z, ζ ))IN  ( ζ )d ζ + (k(z, ζ ))(ζ )dζ 0 0  z − (k(z, ζ ))IN(σ ,ρ ) (ζ )dζ ,

(7.7)

(σ ,ρ ) σ

e(z ) = (IN

(7.8)

0

where (σ ,ρ )

eN,N (k(z, ζ )) = k(z, ζ ) − IN,N (k(z, ζ )), or

e(z ) = I1 + I2 + I3 ,

(7.9)

where (σ ,ρ ) σ

D  ( z ) − Dσ  ( z ) .

I1 = IN  I2 =

z 0

(σ ,ρ )

k(z, ζ )((ζ ) − IN



(ζ ))dζ ,

I3 =

z 0

(σ ,ρ )

eN,N k(z, ζ )IN

 ( ζ )d ζ .

From the Grönwall inequality we can write

 e(z ) L2 (σ ,ρ ) ≤ I1 L2 (σ ,ρ ) +  I2 L2 (σ ,ρ ) +  I3 L2 (σ ,ρ ) . w

w

w

(7.10)

w

Based on Lemma 7.2, we have

 I1 L2 (σ ,ρ ) (I ) = IN(σ ,ρ ) Dσ (z ) − Dσ (z )  w ≤ C1 N −m | Dσ  |H m,N (I ) .

(7.11)

w (σ ,ρ )

also

 I2 L2 (σ ,ρ ) (I ) = 



w

≤

z 0

(σ ,ρ )

k(z, ζ )((ζ ) − IN



z 0

(ζ ))dζ

(σ ,ρ )

k(z, ζ )((ζ ) − IN

L2 (σ ,ρ ) (I )

(ζ ))dζ

w

∞ ,

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359

351

 ((ζ ) − IN(σ ,ρ ) (ζ )) ∞ ,

≤λ

≤ λC2 N 2 −m |  |H m,N (I ) , w (σ ,ρ ) 1

(7.12)

where the value of λ is given as [48]

λ = max | k(z, ζ ) | . a≤z≤b

Finally, by using Cauchy–Schwartz inequality, we find

 I3 L2 (σ ,ρ ) (I ) ≤ eN,N k(z, ζ ) L2 ((σ ,ρ )) (I )  IN((σ ,ρ )) (ζ ) L2 (σ ,ρ ) (I ) . w

w

(7.13)

w

moreover, we have [19,45]

 eN,N k(z, ζ ) L2 (σ ,ρ ) (I ) ≤ k(z, ζ ) − IN(σ ,ρ ) k(z, ζ ) L2 (σ ,ρ ) (I ) ≤ C3 N−m | k(z, ζ ) |Hm,N w

w (σ ,ρ )

w

(I ) ,

(7.14)

thus, the desired upper bound for I4 is obtained as

 I3 L2 (σ ,ρ ) (I ) ≤ (C3 N−m | k(z, ζ ) |Hm,N

w (σ ,ρ )

w

(I ) )(C4

  L2 (σ ,ρ ) (I ) ),

(7.15)

w

or, in other words,

 I3 L2 (σ ,ρ ) (I ) ≤ C5 N−m | k(z, ζ ) |Hm,N

w (σ ,ρ )

w

(I ) 

 L2

w (σ ,ρ )

(I )

.

(7.16)

Finally, we have

C1 N −m | Dσ  |H m,N (I ) +λC2 N 2 −m |  |H m,N (I ) + w (σ ,ρ ) w (σ ,ρ ) C5 N −m | k(z, ζ ) |H m,N (I )   L2 (I ) ,

 e(z ) L2 (σ ,ρ ) (I ) ≤

1

w

or



N −m C1 | Dσ  |H m,N (I ) +λC2 N 2 |  |H m,N (I ) + w (σ ,ρ )

w(σ ,ρ ) C5 | k(z, ζ ) |H m,N (I )   L2 (I ) ,

 e(z ) L2 (σ ,ρ ) (I ) ≤ w

(7.17)

w (σ ,ρ )

w (σ ,ρ )

1

(7.18)

w (σ ,ρ )

w (σ ,ρ )

this completes the proof of the theorem. 8. Numerical results This Section presents several examples to illustrate the accuracy and effectiveness of the proposed SJ-G-C method. We define the absolute error, E(z), as

E ( z ) =|  ( z ) − N ( z ) |,

(8.1)

where (z) and N (z) are the exact and the approximate solutions at point z, respectively. The maximum absolute error, MAE, is given by

MAE = max{E (z )}.

(8.2)

Moreover, we define the root mean square error, RMSE, as

  em ( z ) 2 =

N  i=0

(σ ,ρ )

(σ ,ρ )

(zL,N,i ) − N (zL,N,i )

2

 12 .

N+1

(8.3)

Example 1. Firstly, we introduce the linear Volterra FI-DE [25] 2

1

D 3  (z ) =

3 z3 2 2 − 1 + ez − z 2 ez + 2 ( 23 )



z 0

z2 ezt (ζ )dζ ,

(8.4)

with the nonlocal condition

 ( 0 ) + 2 ( 1 ) + 3



1 0

ζ  ( ζ )d ζ = 3,

(8.5)

knowing that the exact solution is (z ) = z. Let

N ( z ) =

N  j=0

(σ ,ρ )

a j PL, j

( z ),

(8.6)

352

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359 Table 1 The MAE for Example 1. z

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Transform method

SJ-G-C method at N = 20

[25] at N = 20

σ = ρ = − 12

−5

1.425 × 10 1.125 × 10−5 8.282 × 10−6 5.369 × 10−6 2.527 × 10−6 2.490 × 10−7 3.0 0 0 × 10−6 5.813 × 10−6 8.854 × 10−6 1.239 × 10−5 1.681 × 10−5

σ =ρ=0

−16

σ = 0, ρ = −16

1.688 × 10 3.821 × 10−17 1.753 × 10−17 1.778 × 10−17 2.565 × 10−17 3.056 × 10−17 4.988 × 10−17 2.511 × 10−19 2.721 × 10−17 2.590 × 10−17 3.747 × 10−18

1.705 × 10 9.099 × 10−17 7.254 × 10−17 5.449 × 10−17 3.044 × 10−17 2.255 × 10−17 8.412 × 10−18 5.598 × 10−17 6.493 × 10−17 1.638 × 10−16 1.079 × 10−16

1 2

−16

1.662 × 10 7.089 × 10−17 1.002 × 10−16 7.244 × 10−17 2.645 × 10−17 1.026 × 10−18 3.949 × 10−17 4.729 × 10−17 4.729 × 10−17 1.085 × 10−16 1.006 × 10−16

Fig. 1. The absolute error, E, versus z of Example 1 for N = 20 and σ = ρ = − 12 .

Using Eqs. (3.7)–(3.11), we obtain



N   j j=0

k=0





a j χ ( j, k, σ , ρ ) (kk++12 ) zk− 3 − (−1 )k z1−k  k + 1, −z2 − (k + 1 ) ( 3) 1



=

2

3 z3 2 ( 2 ) 3

− 1 + ez

2

(8.7)

2

j N  

 a j χ ( j, k, σ , ρ )

j=0 k=0

(k + ρ + 1 ) 3 ( (−1 )k + 2 ) + (ρ + 1 ) k! k+2

−z2 ez ,

 = 3,

(8.8)

) j−k ( j+ρ +1 )( j+k+σ +ρ +1 ) where χ ( j, k, σ , ρ ) = (−1 . In the suggested technique, the residual of Eq. (8.7) can be zero at N of k!( j−k )!(k+ρ +1 )( j+σ +ρ +1 ) shifted Jacobi–Gauss points that together with the nonlocal condition (8.8) provide a linear system of algebraic equations, which may then be solved by any standard numerical technique. Table 1 compares the results obtained with the SJ-G-C and the transform methods [25], for σ = ρ = −1/2, σ = ρ = 0 and σ = 0, ρ = 1/2, showing that the SJ-G-C is more accurate. Fig. 1 depicts the absolute error, E, versus z, for N = 20 and σ = ρ = − 12 . Fig. 2 shows the approximate and exact solutions, N (z) and (z), versus z, illustrating the good match between them.

Example 2. In this example, we discuss the non-linear Volterra FI-DE with weakly singular kernel 1

D 2  (z ) = F (z ) +



z 0

 (ζ )



z−ζ

dζ ,

(8.9)

with the nonlocal condition

 (0 ) +  (1 ) − 3



1 0

ζ 2 (ζ )dζ = 0.5384615384615384,

the function F (z ) may be chosen such that the exact solution is a non-smooth (z ) = z3.5 .

(8.10)

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359

353

Fig. 2. The approximate, N (z), and the exact, (z), solutions of Example 1 for N = 20 and σ = ρ = − 12 . Table 2 The MAE of Example 2.

σ = ρ = − 12

σ =ρ=0

N

MAE

Time (s)

N

MAE

Time (s)

4 8 12 16 20

1.31812 × 10−3 6.17413 × 10−6 3.43769 × 10−7 5.02712 × 10−8 1.12635 × 10−8

4.877 9.89 21.25 43.889 100.108

4 8 12 16 20

1.64618 × 10−3 9.85145 × 10−6 5.33433 × 10−7 7.36237 × 10−8 1.70775 × 10−8

5.373 11.422 23.374 47.093 113.547

Table 3 The RMSE of Example 3.

σ, ρ

SJ-G-C method

σ = ρ = − 12 σ =ρ=0 σ = 0, ρ = 12

Method in [26]

N=2

N=4

N=8

N = 12

N = 24

N = 48

3.971 × 10−16 3.348 × 10−16 3.211 × 10−17

5.275 × 10−17 1.553 × 10−17 6.001 × 10−17

4.259 × 10−17 4.281 × 10−17 6.378 × 10−17

7.2839 × 10−5 – –

1.051 × 10−5 – –

1.513 × 10−6 – –

Table 2 summarizes the results obtained by the SJ-G-C method in terms of MAE, for different choices of N, σ and ρ . The computational time was also included in Table 2. Example 3. In this example, we discuss the non-linear Volterra FI-DE [26] 4

6

D 5  (z ) =

5z 5

2( 45 )

−

z9 + 252



z 0

( z − ζ ) 2 [  ( ζ )] 3 d ζ ,

(8.11)

with the initial condition 

 ( 0 ) =  ( 0 ) = 0,

(8.12)

knowing that the exact solution is (z ) =

z2 .

Table 3 summarizes the results obtained by the SJ-G-C method in terms of root mean square error, RMSE, for different choices of N, σ and ρ . The values are much lower than the best RMSE = 1.513 × 10−6 obtained in [26] for N = 48. Fig. 3 depicts E versus z, while Fig. 4 shows the exact and approximate solutions, when considering the values of parameters N = 8, σ = 0 and ρ = 12 . Example 4. We consider the following non-linear Fredholm FI-DE [29] 1

5

D 3  (z ) =

6z 3

( ) 1 3

−

z2 z 1 − − + 7 4 9



1 0

( z + ζ ) 2 [  ( ζ )] 3 d ζ ,

(8.13)

with the initial condition 

 ( 0 ) =  ( 0 ) = 0,

(8.14)

354

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359

Fig. 3. The E versus z in Example 2 for N = 8, σ = 0 and ρ =

1 . 2

Fig. 4. The approximate, N (z), and the exact, (z), solutions versus z of Example 2 for N = 8, σ = 0 and ρ =

1 . 2

Table 4 The RMSE for Example 4.

σ, ρ

SJ-G-C method

Method in [29]

N=2

σ =ρ= σ =ρ=0 σ = 0, ρ =

− 12 1 2

N=6 −16

N=8 −16

3.347 × 10 3.348 × 10−16 1.024 × 10−15

2.307 × 10 7.973 × 10−17 3.522 × 10−17

N=8 −16

1.462 × 10 6.030 × 10−17 6.244 × 10−17

3.186 × 10 – –

N = 16 −5

6.157 × 10 – –

N = 32 −6

2.490 × 10−7 – –

knowing that the exact solution is (z ) = z. Table 4 compares the RMSE for the new SJ-G-C and the second kind Chebyshev wavelet methods [29]. We verify that, for the later, the best RMSE = 2.490 × 10−7 , verified for N = 32, is much higher than those obtained for the SJ-G-C scheme. Fig. 5 depicts the E versus z for N = 8, σ = 0 and ρ = 12 , while Fig. 6 compares graphically the exact and the numerical solutions. Example 5. We consider the system of Volterra FI-DEs [30]

⎧  z ⎪ σ  (z ) = 1 + z + z2 −  (z ) + ⎪ D [  1 ( ζ ) +  2 ( z )] d ζ , 1 2 ⎨ 0

 z ⎪ ⎪ ⎩Dσ 2 (z ) = −1 − z − 1 (z ) + [1 (ζ ) − 2 (z )]dζ , 0

(8.15)

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359

Fig. 5. The E versus z in Example 3 for N = 8, σ = 0 and ρ =

355

1 . 2

Fig. 6. The approximate, N (z), and the exact, (z), solutions versus z of Example 3, for N = 8, σ = 0 and ρ =

1 . 2

Table 5 The absolute errors E1 and E2 for Example 5. z

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Method in [30]

SJ-G-C method

N=3 E1

N=3 E2

N=5 E1

N=5 E2

N=5 E1

N=5 E2

N = 13 E1

N = 13 E2

1.1 × 10−3 1.2 × 10−3 7.0 × 10−4 2.0 × 10−4 1.0 × 10−4 1.0 × 10−4 3.0 × 10−4 7.0 × 10−4 8.0 × 10−4

6.0 × 10−4 5.0 × 10−4 2.0 × 10−3 1.0 × 10−4 2.0 × 10−4 2.0 × 10−4 7.0 × 10−4 1.3 × 10−3 1.3 × 10−3

1.0 × 10−5 6.4 × 10−4 9.3 × 10−5 7.2 × 10−5 1.5 × 10−5 4.6 × 10−5 7.7 × 10−5 5.7 × 10−5 5.1 × 10−6

1.6 × 10−5 6.3 × 10−6 1.2 × 10−5 5.0 × 10−6 3.8 × 10−5 6.7 × 10−5 7.3 × 10−5 4.6 × 10−5 3.1 × 10−6

4.9 × 10−7 1.7 × 10−8 4.0 × 10−7 9.2 × 10−7 6.2 × 10−7 1.7 × 10−7 8.2 × 10−8 1.5 × 10−6 5.8 × 10−7

4.9 × 10−7 1.7 × 10−8 4.0 × 10−8 9.2 × 10−7 6.2 × 10−7 1.7 × 10−7 8.2 × 10−8 1.5 × 10−6 5.8 × 10−7

4.1 × 10−16 8.3 × 10−17 2.8 × 10−16 2.3 × 10−16 1.9 × 10−16 1.8 × 10−17 2.6 × 10−16 3.9 × 10−16 2.7 × 10−16

4.1 × 10−16 8.3 × 10−17 2.8 × 10−16 2.3 × 10−16 1.9 × 10−16 1.8 × 10−17 2.6 × 10−16 3.9 × 10−16 2.7 × 10−16

with the initial condition

1 (0 ) = 1, 2 (0 ) = −1, where, for σ = 1, we have 1 (z ) = z +

(8.16) ez ,

2 ( z ) = z

− ez .

Applying the technique developed in Section 5 with different values of N, we observe that the new SJ-G-C algorithm is more accurate than the Chebyshev pseudo-spectral method [30] (see Table 5, which summarizes the values of the absolute

356

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359 Table 6 The maximum absolute errors MAE1 and MAE2 of Example 5 for σ = 0 and ρ = N

3

MAE1 MAE2

5 −3

4.5 × 10 5.2 × 10−3

7 −5

9 −8

1.7 × 10 1.9 × 10−5

2.8 × 10 3.2 × 10−8

11 −11

2.8 × 10 3.0 × 10−11

1 . 2

13 −14

1.7 × 10 1.8 × 10−14

4.6 × 10−16 4.1 × 10−16

Fig. 7. The absolute error E1 versus z of Example 4, for N = 13, σ = 0 and ρ =

1 . 2

Fig. 8. The absolute error E2 versus z of Example 4, for N = 13, σ = 0 and ρ =

1 . 2

errors E1 and E2 ). The MAE values are listed in Table 6 for σ = 0 and ρ = 12 . Figs. 7 and 8 depict the absolute errors, E1 and E2 , versus z, respectively, for N = 13, σ = 0 and ρ = 12 . Fig. 9 illustrates the good matching between the exact and approximate solutions. In Fig. 10, we plot the log10 (MAE) at various values of N, demonstrating the exponential convergence of the method. Example 6. Finally, we discuss the linear mixed Volterra–Fredholm FI-DE [25] 1

D 2 (z ) = −z2

1 ez 1 1  (z ) − z2 + z 2 + ez 3 2 ( 32 )



z 0

t(ζ )dζ +



z 0

z2 (ζ )dζ ,

(8.17)

with the nonlocal condition

 (0 ) +  (1 ) − 3



1 0

t(ζ )dζ = 0,

(8.18)

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359

Fig. 9. The approximate, N (z), and exact, (z), solutions of Example 4 for N = 13, σ = 0 and ρ =

357

1 . 2

Fig. 10. The log10 (MAE) versus N for Example 4. Table 7 The absolute error E for Example 6.

z 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Method [25]

SJ-G-C method at N = 12

at N = 20

σ = ρ = − 12 −7

9.666 × 10 9.589 × 10−7 9.236 × 10−7 8.497 × 10−7 7.311 × 10−7 5.655 × 10−7 3.539 × 10−7 1.003 × 10−7 1.907 × 10−7 5.162 × 10−7 8.809 × 10−7

−27

6.070 × 10 2.610 × 10−28 7.373 × 10−28 7.953 × 10−28 5.548 × 10−28 1.326 × 10−27 6.024 × 10−28 7.538 × 10−28 7.020 × 10−28 1.425 × 10−28 2.032 × 10−27

σ =ρ=0

σ = 0, ρ = −27

6.069 × 10 2.102 × 10−28 7.763 × 10−28 7.661 × 10−28 5.868 × 10−28 1.281 × 10−27 6.257 × 10−28 7.711 × 10−28 7.061 × 10−28 1.441 × 10−28 2.035 × 10−27

1 2

−27

6.069 × 10 2.102 × 10−28 7.763 × 10−28 7.661 × 10−28 5.868 × 10−28 1.281 × 10−27 6.257 × 10−28 7.711 × 10−28 7.061 × 10−28 1.441 × 10−28 2.040 × 10−27

knowing that the exact solution is given by (z ) = z. The absolute errors are listed in Table 8.7, showing that the SJG-C method is more accurate than the transform method [25]. 9. Conclusion In this paper efficient numerical techniques based on the SJ-G-C method were developed for solving FI-DEs subject to initial and nonlocal conditions. The SG-G points were adopted for collocation nodes and the FI-DEs reduced to systems of algebraic equations. Spectral methods are promising candidates for solving many problems, since their global nature fits well with the nonlocal definition of fractional operators. Spectral methods can be used to solve linear and nonlinear differ-

358

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359

ential equations, variational and optimal control problems, and integral and integro-differential equations, both of integer and fixed or variable fractional order. It should be pointed out that the proposed collocation method can also accommodate other numerical methods. For instance, if the proposed collocation scheme suffered from the Runge phenomenon, we could prevent that by simply replacing the basis functions in space by low-order finite-element or radial basis functions, and deriving the corresponding scheme by following the procedure described. Error analysis was performed and numerical examples demonstrated that the method is simple and accurate, even when using a limited number of nodes. The codes used in this paper were developed using the MATHEMATICA program. References [1] Bhrawy AH, Alzaidy JF, Abdelkawy MA, Biswas A. Jacobi spectral collocation approximation for multi-dimensional time-fractional Schrödinger equations. Nonlinear Dyn 2016;84:1553–67. [2] Bhrawy AH, Abdelkawy MA, Ezz-Eldien SS. Efficient spectral collocation algorithm for a two-sided space fractional Boussinesq equation with non-local conditions. Mediterr J Math 2016;13:2483–506. [3] Podlubny I. Fractional differential equations. Mathematics in science and engineering. San Diego: Academic Press Inc.; 1999. [4] Hilfer R. Applications of fractional calculus in physics. Singapore: Word Scientific; 20 0 0. [5] Atangana A, Alabaraoye E. Solving a system of fractional partial differential equations arising in the model of HIV infection of CD4+ cells and attractor one-dimensional Keller-Segel equations. Adv Differ Equ 2013;2013:1–12. [6] Giona M, Roman HE. Fractional diffusion equation for transport phenomena in random media. Physica A 1992;185:87–97. [7] Kirchner JW, Feng X, Neal C. Fractal stream chemistry and its implications for contaminant transport in catchments. Nature 20 0 0;403:524–6. [8] Magin RL. Fractional calculus in bioengineering. Begell House Publishers; 2006. [9] Meerschaert MM, Tadjeran C. Finite difference approximations for two-sided space-fractional partial differential equations. Appl Numer Math 2006;56:80–90. [10] Moghaddam BP, Machado JAT. A stable three-level explicit spline finite difference scheme for a class of nonlinear time variable order fractional partial differential equations. Comput Math Appl 2017;73:1262–9. [11] Yang XJ, Machado JAT, Srivastava HM. A new numerical technique for solving the local fractional diffusion equation: two-dimensional extended differential transform approach. Appl Math Comput 2016;274:143–51. [12] Zhang Y, Yang XJ. An efficient analytical method for solving local fractional nonlinear PDEs arising in mathematical physics. Appl Math Model 2016;40:1793–9. [13] Doha EH, Bhrawy AH, Ezz-Eldien SS. Efficient Chebyshev spectral methods for solving multi-term fractional orders differential equations. Appl Math Model 2011;35:5662–72. [14] Zayernouri M, Ainsworth M, Karniadakis GE. A unified Petrov-Galerkin spectral method for fractional PDEs. Comput Methods Appl Mech Eng 2015;283:1545–69. [15] Maleknejad K, Mahmoudi Y. Taylor polynomial solution of high-order nonlinear Volterra–Fredholm integro-differential equations. Appl Math Comput 2003;145:641–53. [16] Maleknejad K, Asgari M. The construction of operational matrix of fractional integration using triangular functions. Appl Math Model 2015;39:1341–51. [17] Rawashdeh EA. Numerical solution of fractional integro-differential equations by collocation method. Appl Math Comput 2006;176:1–6. [18] Momani S, Noor MA. Numerical methods for fourth-order fractional integro-differential equations. Appl Math Comput 2006;182:754–60. [19] Mokhtary P, Ghoreishi F. The l2 -convergence of the Legendre spectral tau matrix formulation for nonlinear fractional integro-differential equations. Numer Algorithm 2011;58:475–96. [20] Zhu L, Fan Q. Numerical solution of nonlinear fractional-order Volterra integro-differential equations by SCW. Commun Nonlinear Sci Numer Simul 2013;18:1203–13. [21] Sayevand K. Analytical treatment of Volterra integro-differential equations of fractional order. Appl Math Model 2015;39:4330–6. [22] Eslahchi MR, Dehghanb M, Parvizi M. Application of the collocation method for solving nonlinear fractional integro-differential equations. J Comput Appl Math 2014;257:105–28. [23] Ma X, Huang C. Spectral collocation method for linear fractional integro-differential equations. Appl Math Model 2014;38:1434–48. [24] Heydari MH, Hooshmandasl MR, Mohammadi F, Cattani C. Wavelets method for solving systems of nonlinear singular fractional Volterra integro-differential equations. Commun Nonlinear Sci Numer Simul 2014;19:37–48. [25] Nazari D, Shahmorad S. Application of the fractional differential transform method to fractional-order integro-differential equations with nonlocal boundary conditions. J Comput Appl Math 2010;234:883–91. [26] Jiang W, Tian T. Numerical solution of nonlinear Volterra integro-differential equations of fractional order by the reproducing kernel method. Appl Math Model 2015;39:4871–6. [27] Saeedi H, Moghadam MM. Numerical solution of nonlinear Volterra integro-differential equations of arbitrary order by CAS wavelets. Commun Nonlinear Sci Numer Simul 2011;16:1216–26. [28] Nazari Susahab D, Shahmorad S, Jahanshahi M. Efficient quadrature rules for solving nonlinear fractional integro-differential equations of the Hammerstein type. Appl Math Model 2015;39:5452–8. [29] Zhu L, Fan Q. Solving fractional nonlinear Fredholm integro-differential equations by the second kind Chebyshev wavelet. Commun Nonlinear Sci Numer Simul 2012;17:2333–41. [30] Khader MM, Sweilam NH. On the approximate solutions for system of fractional integro-differential equations using Chebyshev pseudo-spectral method. Appl Math Model 2013;222:255–64. [31] Ma X, Huang C. Numerical solution of fractional integro-differential equations by a hybrid collocation method. Appl Math Comput 2010;217:480–7. [32] Doha EH, Abdelkawy MA, Amin AZM, Baleanu D. Spectral technique for solving variable-order fractional Volterra integro-differential equations. Numer Methods Partial Differ Equ 2018;34:1659–77. [33] Doha EH, Abdelkawy MA, Amin AZM, Lopes AM. On spectral methods for solving variable-order fractional integro-differential equations. Comp Appl Math 2017;37:3937–50. [34] Bhrawy AH, Abdelkawy MA, Baleanu D, Amin AZM. A spectral technique for solving two-dimensional fractional integral equations with weakly singular kernel. Hacet J Math Stat 2017;47:553–66. [35] Zhang X, Tang B, Hea Y. Homotopy analysis method for higher-order fractional integro-differential equations. Comput Math Appl 2011;62:3194–203. ˚ [36] Aahin N, YzbaA˚ S, Gülsu M. A collocation approach for solving systems of linear Volterra integral equations with variable coefficients. Comput Math Appl 2011;62:755–69. [37] Abdelkawy MA, Doha EH, Bhrawy AH, M Amin AZ. Efficient pseudospectral scheme for 3d integral equations. Prog Proc Roman Acad Ser A 2017;18:199–206. [38] Abdelkawy MA, Amin AZM, Bhrawy AH, Machado JAT, Lopes A. Jacobi collocation approximation for solving multi-dimensional Volterra integral equations. Int J Nonlinear Sci Num 2017;18:411–25. [39] Bhrawy AH, Abdelkawy MA, Machado JT, Amin AZM. Legendre-Gauss-Lobatto collocation method for solving multi-dimensional Fredholm integral equations. Comput Math Appl 2016. doi:10.1016/j.camwa.2016.04.011. [40] Ma X, Huang C. Spectral collocation method for linear fractional integro-differential equations. Appl Math Model 2014;38:1434–48.

E.H. Doha, M.A. Abdelkawy and A.Z.M. Amin et al. / Commun Nonlinear Sci Numer Simulat 72 (2019) 342–359

359

[41] Lau SR, Price RH. Sparse spectral-tau method for the three-dimensional helically reduced wave equation on two-center domains. J Comput Phys 2012;231:7695–714. [42] Ghoreishi F, Yazdani S. An extension of the spectral tau method for numerical solution of multi-order fractional differential equations with convergence analysis. Comput Math Appl 2011;61:30–43. [43] Doha EH, Bhrawy AH. An efficient direct solver for multidimensional elliptic robin boundary value problems using a Legendre spectral-Galerkin method. Comput Math Appl 2012;64:558–71. [44] Doha EH, Bhrawy AH, Hafez RM. A Jacobi-Jacobi dual-Petrov-Galerkin method for third- and fifth-order differential equations. Math Comput Model 2011;53:1820–32. [45] Canuto C, Hussaini MY, Quarteroni A, Zang TA. Spectral methods: fundamentals in single domains. Berlin, Heidelberg: Springer-Verlag; 2006. [46] Yang Y. Jacobi spectral Galerkin methods for fractional integro-differential equations. Calcolo 2015;52:519–42. [47] Bhrawy AH, Zaky MA, Gorder RAV. A space-time Legendre spectral tau method for the two-sided space-time Caputo fractional diffusion-wave equation. Numer Algorithm 2016;71:151–80. [48] Eslahchi MR, Dehghan M, Parvizi M. Application of the collocation method for solving nonlinear fractional integro-differential equations. J Comput Appl Math 2014;257:105–28.