Approximate solution of linear and nonlinear fractional differential equations under m-point local and nonlocal boundary conditions

Abstract

This paper investigates a computational method to find an approximation to the solution of fractional differential equations subject to local and nonlocal m-point boundary conditions. The method that we will employ is a variant of the spectral method which is based on the normalized Bernstein polynomials and its operational matrices. Operational matrices that we will developed in this paper have the ability to convert fractional differential equations together with its nonlocal boundary conditions to a system of easily solvable algebraic equations. Some test problems are presented to illustrate the efficiency, accuracy, and applicability of the proposed method.

Khalil et al. Advances in Dierence Equations (2016) 2016:177 DOI 10.1186/s13662-016-0910-7

R E S E A R C H Open Access

http://crossmark.crossref.org/dialog/?doi=10.1186/s13662-016-0910-7&domain=pdf

Web End = Approximate solution of linear and nonlinear fractional differential equations under m-point local and nonlocal boundary conditions

Hammad Khalil1,2* , Rahmat Ali Khan3, Dumitru Baleanu4 and Samir H Saker5

*Correspondence: mailto:[email protected]

Web End [email protected]

1Department of Mathematics, University of Poonch Rawalakot, Rawalakot, 12350, Pakistan

2Department of Mathematics, University of Malakand, P.O. Box 18000, Chakdara, Dir Lower, Khybarpukhtunkhwa, Pakistan Full list of author information is available at the end of the article

Abstract

This paper investigates a computational method to nd an approximation to the solution of fractional dierential equations subject to local and nonlocal m-point boundary conditions. The method that we will employ is a variant of the spectral method which is based on the normalized Bernstein polynomials and its operational matrices. Operational matrices that we will developed in this paper have the ability to convert fractional dierential equations together with its nonlocal boundary conditions to a system of easily solvable algebraic equations. Some test problems are presented to illustrate the eciency, accuracy, and applicability of the proposed method.

MSC: 35C11; 65T99Keywords: Bernstein polynomials; operational matrices; m-point boundary conditions; fractional dierential equations

1 Introduction

Recently the studies of fractional dierential equations (FDEs) gained the attention of many scientists around the globe. This topic remains a central point in several special issues and books. Fractional-order operators are nonlocal in nature and due to this property, they are most nicely applicable to various systems of natural and physical phenomena. This property has motivated many scientists to develop fractional-order models by considering the ideas of fractional calculus. Examples of such systems can be found in many disciplines of science and engineering such as physics, biomathematics, chemistry, dynamics of earthquakes, dynamical processes in porous media, material viscoelastic theory, and control theory of dynamical systems. Furthermore, the outcome of certain observations indicates that fractional-order operators possess some properties related to systems having long memory. For details of applications and examples, we refer the reader to the work in [].

The qualitative study of FDEs which discusses analytical investigation of certain properties like existence and uniqueness of solutions has been considered by several authors. In Guptta [] studied the solvability of three point boundary value problems (BVPs).

2016 Khalil et al. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 2 of 28

Since then many researchers have been working in this area and provided many useful results which guarantee the solvability and existence of a unique solution of such problems. For the reader interested in the existence theory of such problem we refer to work presented by Ruyun Ma [] in which the author presents a detailed survey on the topic. In [] the author derived an analytic relation which guarantees the existence of positive solution of a general third order multi-point BVPs. Also some analytic properties of solutions of FDEs are discussed by El-Sayed in []. Often it is impossible to arrive at the exact solution when FDEs has to be solved under some constraints in the form of boundary conditions. Therefore the development of approximation techniques remains a central and active area of research.

The spectral methods, which belong to the approximation techniques, are often used to nd approximate solution of FDEs. The idea of the spectral method is to convert FDEs to a system of algebraic equations. However, dierent techniques are used for this conversion. Some of well-known techniques are the collocation method, the tau method, and the Galerkin method. The tau and Galerkin methods are analogous in the sense that the FDEs are enforced to satisfy some algebraic equations, then some supplementary set of equations are derived using the relations of boundary conditions (see, e.g., [] and the references therein). The collocation method [, ], which is an analog of the spectral method, consists of two steps. First, a discrete representation of the solution is chosen and then FDEs are discretized to obtain a system of algebraic equations.

These techniques are extensively used to solve many scientic problems. Doha et al. [] used collocation methods and employed Chebyshev polynomials to nd an approximate solution of initial value problems of FDEs. Similarly, Bhrawy et al. [] derived an explicit relation which relates the fractional-order derivatives of Legendre polynomials to its series representation, and they used it to solve some scientic problems. Saadatmandi and Dehghan [] and Doha et al. [] extended the operational matrices method and derived operational matrices of fractional derivatives for orthogonal polynomials and used it for solving dierent types of FDEs. Some recent and good results can be found in the articles like [, , ].

Some other methods have also been developed for the solution of FDEs. Among others, some of them are iterative techniques, reproducing kernel methods, nite dierence methods etc. Esmaeili and Shamsi [] developed a new procedure for obtaining an approximation to the solution of initial value problems of FDEs by employing a pseudo-spectral method, and Pedas and Tamme [] studied the application of spline functions for solving FDEs. In [], the author used a quadrature tau method for obtaining the numerical solution of multi-point boundary value problems. Also the authors in [] extended the spectral method to nd a smooth approximation to various classes of FDEs and FPDEs. Some recent results in which orthogonal polynomials are applied to solve various scientic problems can be found in [].

Multi-point nonlocal boundary value problems appears widely in many important scientic phenomena like in elastic stability and in wave propagation. For the solution of such a problem Rehman and Khan [, ] introduced an ecient numerical scheme, based on the Haar wavelet operational matrices of integration for solving linear multi-point boundary value problems for FDEs. FDEs subject to multi-point nonlocal boundary conditions are a little bit dicult. In this area of research a few articles are available. Some good results on solution of nonlocal boundary value problems can be found in [].

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 3 of 28

Bernstein polynomials are frequently used in many numerical methods. Bernstein polynomials enjoy many useful properties, but they lack the important property of orthogonality. As the orthogonality property is one of more important properties in approximation theory and numerical simulations, these polynomials cannot be directly implemented in the current technique of numerical approximations. To overcome this diculty, these non-orthogonal Bernstein polynomials are transformed into an orthogonal basis []. But as the degree of the polynomials increases the transformation matrix becomes ill conditioned [, ], which results some inaccuracies in numerical computations. Recently Bellucci introduced an explicit relation for normalized Bernstein polynomials []. One applied Gram-Schmidt orthonormalization process to some sets of Bernstein polynomials of dierent scale levels, identifying the pattern of polynomials, and generalizing the result. The main results presented in this article are based on these generalized bases.

In this article we present an approximation procedure to nd an approximate solution of the FDEs subject to local and nonlocal m-point BVPs. The method is designed to solve linear FDEs with constant coecients, linear FDEs with variable coecients, and nonlinear FDEs, under local and nonlocal m-point boundary conditions. In particular, we consider the following generalized class of FDEs:

D U(t) =

p

i=

p

i=

iDiU(t) = F U(t), U(t), . . . , U(t)p + f (t). ()

In equation () i R, t [, ], f (t) C([, ]). The orders of the derivatives are dened

as

< < p < p p < p + .

In equation () i(t) C([, ]). In () F(U(t), U(t), . . . , U(t)p) is a nonlinear function

of U(t) and its fractional derivatives.

The main aim in this paper is to nd a smooth approximation to U(t), which satises a given set of m-point boundary conditions. We consider the following two types of boundary constraints.

Type : Multi-point local boundary conditions dened as

U() = u, U(i) = ui, i = , , p , U() = up. ()

Type : Non local m-point boundary conditions, dened as

Uj() = uj, j = , , . . . , p ,

m

i=

iU(i) = U(). ()

p

i=

iDiU(t) + f (t), ()

D U(t) =

i(t)DiU(t) + f (t), ()

D U(t) +

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 4 of 28

We use the normalized Bernstein polynomials for our investigation, which is based on the explicit relations presented in []. We use these polynomials to develop new operational matrices. We develop four operational matrices, two of them being operational matrices of integration and dierentiation, the formulation technique for these two operational matrices is the same as that used for traditional orthogonal polynomials. We introduce two more operational matrices to deal with the local and nonlocal boundary conditions. To the best of our knowledge such types of operational matrices are not designed for normalized Bernstein polynomials.

We organized the rest of article as follows. In Section , we recall some basic concepts denition from fractional calculus, approximation theory, and matrix theory. Also we present some properties of normalized Bernstein polynomials which are helpful in our further investigation. In Section , we present a detailed procedure for the construction of the required operational matrices. In Section , the developed matrices are employed to solve FDEs by introducing a new algorithm. In Section , the proposed algorithm is applied to some test problems to show the eciency of the proposed algorithm. The last section is devoted to a short conclusion.

2 Some basic denitions

In this section, we present some basic notation and denitions from fractional calculus and well-known results which are important for our further investigation. More details can be found in [, ].

Denition The Riemann-Liouville fractional-order integral of order R+ of a func

tion (t) (L[a, b], R) is dened by

aIt (t) =

()

t

(s) ds. ()

The integral on right hand side exists and is assumed to be convergent.

Denition For a given function (t) Cn[a, b], the fractional-order derivative in the

Caputo sense, of order is dened as

D (t) =

(n )

a (t s)

t

(n)(s)

(t s)+n ds, n

< n, n N.

The right side is assumed to be point wise dened on (a, ), where n = [] + in the case

that is not an integer.

This leads to I tk =

(+k)

(+k+) tk+ for > , k , D C = , for a constant C and

D tk =

a

( + k)

( + k )tk , for k [

]. ()

Throughout the paper, we will use the Brenstein polynomials of degree n. The analytic relation of the Brenstein polynomials of degree n dened on [, ] is given as

B(i,n)(t) =

ni

k=

()k n i

n i k

ti+k, i = , , . . . , n. ()

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 5 of 28

The set of polynomials dened by () have a lot of interesting properties. They are positive on [, ] and also approximate a smooth function on the domain [, ]. The polynomials dened in () are not orthogonal, after the application of the Gram-Schmidt process the explicit form of normalized Bernstein polynomials is obtained (as discussed in detail in []) by

j,n(t) = (n j) + ( t)nj

j

k=

(j,k)tjk, ()

where

(j,k) = ()k

(n k + ) (j )

(j k + ) (n j + ) (k + ) (j k + ).

The polynomials dened in () are not directly applicable in construction of operational matrices (it will be clear in the next section). Therefore, we further generalize the relation

j,n(t) = (n j) +

j

k=

nj

l=

(j,k,l) tl+jk, ()

where

(j,k,l) = (()(nj+k+l)

(n j + ) (n k + ) (j + ))

(l + ) (n j l + ) (j k + ) (n j + ) (k + ) (j k + ). ()

Note that in equation (), j,n(t) is of degree n for all choices of j. By analyzing we observe that the minimum power of t is at maximum value of k and minimum value of l, which is . Conversely the maximum power of t is n. These polynomials are orthogonal on the interval [, ]. To make them applicable on the interval [, ], we simply make substitution of t = t without loss of generality. So we can write the orthogonal Bernstein polynomials on the interval [, ] as follows:

j,n(t) = w(j,n)

j

k=

nj

l=

(j,k,l) tl+jk

l+jk , ()

where w(j,n) = (n j) + /. The orthogonality relation for these polynomials is dened as follows:

i,n(t)j,n(t) = (i,j). ()

As usual we can approximate any function f C[, ] in the normalized Bernstein poly

nomial as

f (t) =

n

j=

cjj,n(t), where cj =

j,n(t) dt, ()

which can always be written as

f (t) = HTN N(t), ()

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 6 of 28

where HTN and N(t) are N = n + terms column vector containing coecients and Bernstein polynomials, respectively, and one dened

HTN = [c, c, . . . , cn], () N(t) = ,n(t), ,n(t), . . . , j,n(t) T. ()

As N represents the size of the resulting algebraic equations it is considered as a scale level of the scheme.

The integral of the triple product of fractional-order Legendre polynomials over the domain of interest was recently used in []. There the author used this value to solve fractional dierential equations with variable coecients directly. We use a triple product for Bernstein polynomials to construct a new operational matrix which is of basic importance in solving FDEs with variable coecients. The following theorem is of basic importance.

Theorem The denite integral of the product of three Bernstein polynomials over the domain [, ] is constant and is dened as

a,n(t)b,n(t)c,n(t) dt = (a,b,c), ()

where

(a,b,c) = w((a,b,c),n)

a

k=

na

l=

b

r=

nb

s=

c

p=

nc

q=

(a,k,l)

(b,r,s)

(c,p,q)

(l + a + b + c + q p s r k + )

and w((a,b,c),n) = w(a,n)w(b,n)w(c,n).

Proof Consider the following expression:

a,n(t)b,n(t)c,n(t) dt = w((a,b,c),n)

a

k=

na

l=

b

r=

nb

s=

c

p=

nc

q=

(a,k,l)

(b,r,s)

(c,p,q) l+a+b+c+qpsrk

tl+a+b+c+qpsrk dt.

Evaluating the integral and using the notation (a,b,c) we can write

(a,k,l)

(b,r,s)

(c,p,q)

(l + a + b + c + q p + s r k + ).

Now, we present the inverse of the well-known Vandermonde matrix. The inverse of this matrix will be used when we use operational matrices to solve under local boundary conditions.

(a,b,c) = w((a,b,c),n)

a

k=

na

l=

b

r=

nb

s=

c

p=

nc

q=

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 7 of 28

Theorem [, ] Consider a matrix V dened as

V =

p

p

... ... ... ...p pp pp

, ()

where < < < p. The inverse of this matrix exists and in dened as

V = [b(j,i)], where i, j = , , . . . , p,

where the entries b(j,i) are dened by

b(j,i) =

j ()k

m < < mpk p

m,...,mnk =j

m mnk

m n

. ()

m =j (

m j)

3 Operational matrices of derivative and integral

Now we are in a position to construct new operational matrices. The operational matrices of derivatives and integrals are frequently used in the literature to solve fractional-order dierential equations. In this section, we present the proofs of constructions of four new operational matrices. These matrices act as building blocks in the proposed method.

Theorem The fractional integration of order of the function vector N(t) (as dened in ()) is dened as

I N(t) = P(,)(NN)

N(t),

where P(,)(NN) is an operational matrix for fractional-order integration and is given as

P(,)(NN) =

(,) (,) (,n) (,) (,) (,n)

... ... ... ... (n,) (n,) (n,n)

, ()

where

(r,s) = w(r,n)w(s,n)

s

k =

ns

l =

r

k=

nr

l=

(s,k ,l ) (r,k,l,,)l+rk++

(l + l + s + r k k + + ),

where

(s,k ,l ) is as dened in (), and

(r,k,l,,) =

(r,k,l) (l + r k + ) (l + r k + + )l+rk .

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 8 of 28

Proof Consider the general element of () and apply fractional integral of order , consequently we will get

I r,n(t) = w(r,n)

r

k=

nr

l=

(r,k,l) I tl+rk

l+rk . ()

Using the denition of fractional-order integration we may write

I r,n(t) = w(r,n)

r

k=

nr

l=

(r,k,l,,)tl+rk+ , ()

where (r,k,l,,) =

(r,k,l) (l+rk+) (l+rk++)l+rk . We can approximate tl+rk+ with normalized Bernstein polynomials as follows:

tl+rk+ =

n

s=

c(r,s)s,n(t), where c(r,s) =

tl+rk+

s,n(t) dt. ()

Using equation (), we can write

c(r,s) = w(s,n)

s

k =

ns

l =

(s,k ,l )

l +sk

tl+l +s+rkk + dt. ()

On further simplications, we can get

(s,k ,l ) l+rk++

(l + l + s + r k k + + ). ()

Using () and () in () we get

I r,n(t) =

n

s=

c(r,s) = w(s,n)

s

k =

ns

l =

w(r,n)w(s,n)

s

k =

ns

l =

r

k=

nr

l=

(s,k ,l ) (r,k,l,,)l+rk++

(l + l + s + r k k + + )

s,n(t). ()

Using the notation

(s,k ,l ) (r,k,l,,)l+rk++

(l + l + s + r k k + + ),

and evaluating for r = , , . . . , n and s = , , . . . , n completes proof of the theorem.

Theorem The fractional derivative of order of the function vector N(t) (as dened in ()) is dened as

D N(t) = D(,)(NN)

N(t),

(r,s) = w(r,n)w(s,n)

s

k =

ns

l =

r

k=

nr

l=

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 9 of 28

where D(,)(NN) is dened as

D(,)(NN) =

(,) (,) (,n) (,) (,) (,n)

... ... ... ... (n,) (n,) (n,n)

, ()

where

(r,s) = w(r,n)w(s,n)

s

k =

ns

l =

r

k=

nr

l=

(s,k ,l ) (r,k,l,,)l+rk+

(l + l + s + r k k + )

and

(r,k,l,,) =

(r,k,l) (l+rk+) (l+rk+)l+rk if l + r k , if l + r k < .

()

Proof On application of the derivative of order to a general element of (), we may write

D r,n(t) = w(r,n)

r

k=

nr

l=

(r,k,l) D tl+rk

l+rk . ()

Using the denition of the fractional-order derivative we can easily write

D r,n(t) = w(r,n)

r

k=

nr

l=

(r,k,l,,)tl+rk , ()

where

(r,k,l,,) =

(r,k,l) (l+rk+) (l+rk+)l+rk if l + r k , if l + r k < .

()

We can approximate tl+rk with the Bernstein polynomials as follows:

tl+rk =

n

s=

c(r,s)s,n(t),

c(r,s) =

()

tl+rk

s,n(t) dt.

Using equation (), we can write

c(r,s) = w(s,n)

s

k =

ns

l =

(s,k ,l )

l +sk

tl+l +s+rkk dt. ()

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 10 of 28

On further simplication we can get

(s,k ,l ) l+rk+

(l + l + s + r k k + ). ()

Using () and () in () we get

D r,n(t) =

n

s=

c(r,s) = w(s,n)

s

k =

ns

l =

w(r,n)w(s,n)

s

k =

ns

l =

r

k=

nr

l=

(s,k ,l ) (r,k,l,,)l+rk+

(l + l + s + r k k + + )

s,n(t). ()

Using the notation

(s,k ,l ) (r,k,l,,)l+rk+

(l + l + s + r k k + ),

and evaluating for r = , , . . . , n, and s = , , . . . , n, we complete the proof of the theorem.

The operational matrices developed in the previous theorems can easily solve FDEs with initial conditions. Here we are interested in the approximate solution of FDEs under complicated types of boundary conditions. Therefore we need some more operational matrices such that we can easily handle the boundary conditions eectively.

The following matrix plays an important role in the numerical simulation of fractional dierential equations with variable coecients.

Theorem For a given function f C[, ], and u = HTN

(r,s) = w(r,n)w(s,n)

s

k =

ns

l =

r

k=

nr

l=

N(t), the product of f (t) and order fractional derivative of the function u(t) can be written in matrix form as

f (t)D u(t) = HTNQ(f,,)(NN)

N(t),

where Q(f,,)(NN) = D(,)(NN)R(f,)(NN), and D(,)(NN) is an operational matrix for fractional-order derivative and

R(f,)(NN) =

(,) (,) (,n) (,) (,) (,n)

... ... ... ... (n,) (n,) (n,n)

, ()

where

(r,s) =

n

q=

dq (q,r,s)r, s = , , . . . , n,

and the entries (q,r,s) are dened as in Theorem , and dq are the spectral coecients of the function f (t).

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 11 of 28

Proof Applying Theorem we can write

D u(t) = HTND(,)(NN)

N(t) ()

and

f (t)D u(t) = HTND(,)(NN)

N(t), ()

where

N(t) = f (t),n(t), f (t),n(t), . . . , f (t)n,n(t) T.

Consider the general element of

N(t), and approximate it with normalized Bernstein

polynomials as

f (t)r,n(t) =

n

s=

crss,n(t), ()

where cri can easily be calculated as

crs =

r,n(t)s,n(t) dt. ()

Now as f (t) C[, ], we can approximate it with normalized Bernstein polynomials as

f (t) =

n

q=

dqq,n(t),

dq =

()

q,n(t) dt.

Using equation () in () we get the following estimates:

crs =

n

q=

dq

q,n(t)r,n(t)s,n(t) dt. ()

In view of Theorem we obtain the following estimate:

crs =

n

q=

dq (q,r,s). ()

Using () in ()

f (t)r,n(t) =

n

s=

n

q=

dq (q,r,s)s,n(t). ()

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 12 of 28

Evaluating () for s = , , . . . , n and r = , , . . . , n we can write

f (t),n(t)

f (t),n(t) ...

f (t)n,n(t)

=

(,) (,) (,n) (,) (,) (,n)

... ... ... ... (n,) (n,) (n,n)

,n(t)

,n(t) ... n,n(t)

, ()

where (r,s) = n

q= dq (q,r,s). In simplied notation, we can write

N(t) = R(f,)(NN)

N(t). ()

Using () in () we get the desired result. The proof is complete.

Since one of our aims in this paper is to solve FDEs under dierent types of local and non-local boundary conditions, we have to face some complicated situations, so to handle these situations we will use the operational matrix developed in the next theorem.

Theorem Let f be a function of the form f (t) = atn where a

R and n

N, then for

any function u(t) C[, ], u(t) = HTN

N(t), . Then we can generalize the product of

I u(t) dt and f (t) in matrix form as

f (t)I u(t) dt = HTNW(,,a,n

,) (NN)

N(t),

where

W(,,a,n

,) (NN) =

(,) (,) (,n) (,) (,) (,n)

... ... ... ... (n,) (n,) (n,n)

. ()

The entries of the matrix are dened by

(i,j) = a (i,,,)w(j,n)

j

p=

nj

q=

(j,p,q)

n +

(q + j p + n + ),

where

(l + i k + )l+ik+ (l + i k + + )l+ik .

Proof Consider the general term N(t), if we calculate the order denite integral from to we get

I u(t) dt = I HTN N(t) =

n

i=

(i,,,) = w(i,n)

i

k=

ni

l=

(i,k,l)

ciI i,n(x). ()

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 13 of 28

Using () we may write

I u(t) dt =

n

i=

ciw(i,n)

i

k=

ni

l=

(i,k,l) I xl+ik

l+ik

=

n

i=

ciw(i,n)

i

k=

ni

l=

(i,k,l)

(l + i k + )l+ik+ (l + i k + + )l+ik .

()

Using the notation

(i,,,) = w(i,n)

i

k=

ni

l=

(i,k,l)

(l + i k + )l+ik+ (l + i k + + )l+ik ,

we can write

atn I u(t) dt =

n

i=

cia (i,,,)tn . ()

Now a (i,,,)tn can be approximated with Bernstein polynomials as follows:

a (i,,,)tn =

n

j=

d(i,j)j,n(t), ()

where d(i,j) can be calculated as

d(i,j) = a (i,,,)w(j,n)

j

p=

nj

q=

(j,p,q)

xq+j+n p

q+jp

j

p=

nj

q=

(j,p,q)

n +

= a (i,,,)w(j,n)

(q + j p + n + ). ()

Using the notation (i,j) = d(i,j) and equation () in () we get the desired result. The proof is complete.

4 Application of operational matrices

The operational matrices derived in the previous section play a central role in the numerical simulation of FDEs under dierent types of boundary conditions. We start our discussion with the following class of FDEs:

D U(t) =

p

i=

iDiU(t) + f (t), ()

where i R, t [, ], f (t) C([, ]). The orders of the derivatives are dened as

< < p < p p < p + .

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 14 of 28

The solution to this equation can be assumed in terms of series representation of normalized Bernstein polynomials such that the following relation holds:

D U(t) = HTN N(t). ()

Now by the application of fractional integral of order , and use of Theorem , we can explicitly write U(t) as follows:

U(t)

p

i=

N(t). ()

Our main aim is to obtain HTN, which is an unknown vector. We will use this vector to get the solution to the problem.

Type : Suppose we have to solve () under the following condition:

U() = u, U(i) = ui, i = , , p , U() = up, ()

c can easily be obtained using the initial condition U() = . The solution in () must satisfy all the conditions, therefore using the conditions at the intermediate and boundary points we have

p

i=

citi = HTNP(,)(NN)

cii = (u u) HTNP(,)(NN)

N(),

p

i=

cii = (u u) HTNP(,)(NN)

N(),

... ...

p

i=

()

ciip = (up u) HTNP(,)(NN)

N(p),

p

i=

ciip = (up u) HTNP(,)(NN)

N().

In order to make the notation simple we use p = in last equation. Equation () can then also be written in matrix form as

p

p

... ... ... ...p pp pp

(u u) HTNP(,)(NN)

c c

... cp

N()

(u u) HTNP(,)(NN)

N() ...

(up u) HTNP(,)(NN)

N(p)

=

. ()

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 15 of 28

The matrix in the left is a Vandermonde matrix, in view of Theorem its inverse exists.

Using the inverse of V we get the values of ci:

ci =

p

j=

N(j)". ()

Using the values of ci in equation () we may write

U(t) u

p

i=

b(i, j)!(uj u) HTNP(,)(NN)

p

j=

b(i, j)!(uj u) HTNP(,)(NN)

N(j)"ti = HTNP(,)(NN)

N(t). ()

In simplied notation we can write

U(t) g(t) +

p

i=

p

j=

b(i, j)!tiHTNP(,)(NN)

N(j)" = HTNP(,)(NN)

N(t), ()

where g(t) = u + p

i=

pj= b(i,j)(uj u)ti. Approximating g(t) = GTN N(t) and making use

of Theorem , we can easily get

U(t) + HTN

p

i=

p

j=

!W(,j,b(i,j),i,)

(NN) " N(t) = HTNP(,)(NN)

N(t) + GTN N(t). ()

In simplied notation we can write

U(t) = HTNENN

N(t) + GTN N(t), ()

where

ENN = P(,)(NN)

!W(,j,b(i,j),i,)

(NN) "

is N N matrix related to type boundary conditions. For instance we stop the current

procedure here and discuss the procedure of obtaining analogous relation for U(t) under type boundary conditions.

Type : Suppose if we need to solve the problem under m-point non-local boundary conditions

Uj() = uj, j = , , . . . , p ,

m

i=

p

i=

p

j=

iU(i) = U(). ()

Then using the initial conditions we can nd the rst p constants of integrations in equation (),

U(t) = HTNP(,)(NN)

N(t) +

(p)

l=

ultl + cptp. ()

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 16 of 28

However, the last constant cp is unknown. In order to get a relation for cp we use nonlocal boundary conditions. Using the m-point boundary condition we get

U() = HTNP(,)(NN)

N() +

(p)

l=

ull + cnn, ()

m

i=

iU(i) = HTN

m

i=

iP(,)(NN)

N(i) +

m

i=

(p)

l=

iulli + cp

m

i=

ip. ()

From (), we see that left sides of () and () are equal, therefore for the sake of obtaining the value of cp, we can write

#

HTNP(,)(NN)

N()

m

i=

iP(,)(NN)

N(i)

$ + #(p)

l=

ull

m

i=

(p)

l=

iulli

$

$. ()

On further simplication we get

cp = #HTN

P(,)(NN)

N()

= cp#m

i=

ipi p

$ + , = , ()

where = { mi= ipi p} and = { (n)l= ull mi= (p)l= iulli}. Now, using ()

in () we get

U(t) = HTNP(,)(NN)

N(t) +

m

i=

i

P(,)(NN)

N(i)

(p)

l=

ultl

+ HTN tp

#

P(,)(NN)

N()

m

i=

P(,)(NN)

itp

N(i)

$ + tp. ()

Now, in view of Theorem , we can write () as

U(t) = HTNP(,)(NN)

(,,

,p,) (NN)

N(t) + HTNW

N(t)

m

i=

(,i,

HTNW

i

,p,) (NN)

N(t) + GTN N(t). ()

Here (p)

l= ultl + tp = GTN N(t). On further simplication we can write

U(t) = HTNENN

N(t) + GTN N(t), ()

where

(,,

,p,) (NN)

ENN = P(,)(NN) + W

m

i=

(,i,

W

i

,p,) (NN) .

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 17 of 28

We see that for both local and non-local boundary conditions we get U(t) in matrix form as

U(t) = HTNENN N(t) + GTN N(t). ()

Here HTN is an unknown vector to be determined. Note that we use E = E if type boundary conditions are given, otherwise use E = E. Now using () and using Theorem , we can write

p

i=

iDiU(t) = HTNENN

p

i=

iD(i,)(NN)

N(t) + GTN

p

i=

N(t). ()

Approximating f (t) = FTN N(t), and using (), () in () we may write

HTN N(t) = HTNENN

p

i=

iD(i,)(NN)

iD(i,)(NN)

N(t) + GTN

p

i=

iD(i,)(NN)

N(t) + FTN N(t). ()

After a long calculation and simplication we get

HTN HTNENN

p

i=

iD(i,)(NN) + GTN

p

i=

iD(i,)(NN) FTN = . ()

We see that () is an easily solvable matrix equation and it can easily be solved for HTN, by using HTN in () this will lead us to an approximate solution to the problem.

Next, we consider the linear FDEs with variable coecients,

D U(t) =

p

i=

i(t)DiU(t) + f (t), ()

where i(t) C([, ]), t [, ], f (t) C([, ]). We start our analysis with the initial

assumption

D U(t) = HTN N(t). ()

Using the integral of order we can write

U(t) = HTN N(t) +

p

i=

citi. ()

We can of course get U(t) as

U(t) = HTNENN N(t) + GTN N(t), ()

where ENN, GTN N(t) can be analogously derived as in the previous section depending on the type of boundary conditions to be used. Using () and Theorem , we can easily write

p

i=

i(t)DiU(t) = HTNENN

p

i=

Q(i(t),i,)(NN)

N(t) + GTN

p

i=

Q(i(t),i,)(NN)

N(t). ()

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 18 of 28

Approximating f (t) = FTN N(t), and using (), () in () we may write

HTN N(t) = HTNENN

p

i=

Q(i(t),i,)(NN)

N(t)

+ GTN

p

i=

Q(i(t),i,)(NN)

N(t) + FTN N(t). ()

On rearranging we get

HTN HTNENN

p

i=

Q(i(t),i,)(NN) + GTN

p

i=

Q(i(t),i,)(NN) FTN = . ()

Now, we consider the nonlinear FDEs. The nonlinear FDEs are often very dicult to solve with operational matrix techniques. One way is to use the collocation method, by doing so the FDEs result in nonlinear algebraic equations. These nonlinear algebraic equations are then solved iteratively using the Newton method or some other iterative method. The second approach is to linearize the nonlinear part of the fractional dierential equation using a quasilinearization technique. Doing so the nonlinear fractional dierential equations is converted to a recursively solvable linear dierential equations. This technique was introduced recently and has been used to nd an approximate solution of many scientic problem. The quasilinearization method was introduced by Bellman and Kalaba [] to solve nonlinear ordinary or partial dierential equations as a generalization of the Newton-Raphson method. The origin of this method lies in the theory of dynamic programming. In this method, the nonlinear equations are expressed as a sequence of linear equations and these equations are solved recursively. The main advantage of this method is that it converges monotonically and quadratically to the exact solution of the original equations []. Also some other interesting work in which quasilinearization method is applied to scientic problems is in []. Consider a nonlinear fractional-order dierential equation of the form

D U(t) +

p

i=

iDiU(t) = F U(t), U(t), . . . , U(t)p + f (t). ()

The general procedure of this method is as given now. First solve the linear part with a given set of local or nonlocal conditions

D U(t) +

p

j=

iDiU(t) = f (t). ()

This equation can easily be solved using the method developed in the previous discussion. Label the solution at this step as U(t).

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 19 of 28

The next step is to linearize the nonlinear part with a multivariate Taylor series expansion about U(t). So after linearizing and simplication we get

D U()(t) +

p

j=

jDjU()(t)

= F( U()(t), U()(t), U()(t), . . . , Un()(t)

n

j=

()

The above equation is an FDE with variable coecients. It can easily be solved with the method developed in Section .. The solution at this stage will be labeled U(t) and is the solution of the problem at rst iteration. Again we have to linearize the problem about U(t) to obtain the solution at second iteration. The whole process can be seen as a recurrence relation like

D U(r+)(t) +

n

j=

+ Uj() Ui() f

Uj U

()(t), U()(t), U()(t), . . . , Un()(t) + f (t).

bj(t)DjU(r+)(t) =

F(t) .

And the boundary conditions become

Ujr+() = uj, Ur+() =

m

i=

iUr+(i),

where bj(t) = j

f U

j (U(r)(t), U(r)(t), U(r)(t), . . . , Un(r)(t)) and

F(t) = f U(r)(t), U(r)(t), U(r)(t), . . . , Un(r)(t)

n

j=

Uj(r)

f

Uj + g(t).

It can be easily noted that the above equation is fractional dierential equation with variable coecients. Using the initial solution X(t) we may start the iterations. The coecients bi(t) and the source term

F(t) can be updated at every iteration r to get the next solution at r + . At every step we may solve the problem at given nonlocal boundary conditions. We assume that the method converges to the exact solution of the problem if there is convergence at all.

5 Examples

To show the applicability and eciency of the proposed method, we solve some fractional dierential equations. The numerical simulation is carried out using MatLab. However, we believe that the algorithm can be simulated using any simulation tool kit.

Example As a rst problem we solve the following linear integer order boundary value problem

DU(t) =

DU(t) +

DU(t) +

U(t) + g(t), ()

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 20 of 28

where

& sin(t) cos(t) sin(t) + sin(t) .

We can easily see that the exact solution is U(t) = sin(t). We test our algorithm by solving this problem under the following sets of boundary conditions:

S = 'U() = , U%

& = ., U%

& = , U() = (,

S = !U() = , U(.) = , U(.) = , U() = ",

S = !U() = , U () = , U () = ,

U() + .U(.) + U(.) + .U(.) + .U(.) = U()".

We approximate the solution to this problem with dierent types of boundary conditions and observe that the approximate solution is very accurate. For illustration purposes we calculate the absolute dierence at a dierent scale levels. The results are displayed in Figure and Figure for boundary conditions S and S, respectively. We compare the exact solution and approximate solution under boundary conditions S at dierent scale levels, and the results are displayed in Table . It can easily be noted that with increase of scale level, that the approximate solution becomes more and more accurate, and at scale level N = , the approximate solution is accurate up to the seventh digit. The accuracy may be increased by using a high scale level. For instance we simulate the algorithm at high scale level and measure EN and EN at each scale level N. Table shows these

results at high scale levels under boundary conditions S and S. From this example, we conclude that the proposed method is convergent for integer order dierential equations (linear).

Figure 1 Absolute difference of exact and approximate solution for different values of N under boundary condition S1. Here we x = 3, and use notation EN = |Uexact UN|.

g(t) = %

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 21 of 28

Figure 2 Absolute difference of exact and approximate solution for different values of N under boundary condition S2. Here we x = 3, and use the notation EN = |Uexact UN|.

Table 1 Comparison of exact and approximate solution of Example 1 using the boundary condition as dened in S

t N

Exact U(t) N = 5 N = 7 N = 9 N = 11

t = 0.4 0.95105651629 0.9677391807 0.9499276887 0.9510565133 0.95105651554 t = 0.8 0.58778525229 0.5806739865 0.5877691649 0.5877785192 0.5877882463t = 1.2 0.58778525229 0.5879798137 0.5877896195 0.5877852801 0.58778528235 t = 1.6 0.95105651629 0.9318025727 0.9510687377 0.9510565637 0.9510565256t = 2.0 0.00000000000 0.0160927142 0.0000151536 0.0000001456 0.00000000975 t = 2.4 0.95105651629 0.9684572120 0.9510519403 0.9510565108 0.95105651183 t = 2.8 0.58778525229 0.5667990904 0.5879771600 0.5877852533 0.5877852523t = 3 0.00000000000 0.0634719308 0.0001774663 0.0000000045 0.00000000125

Table 2 Comparison of exact and approximate solution of Example 1 using the boundary condition as dened in S

S1 S2

EN 2 EN EN 2 EN

N = 5 7.23 102 6.53 101 4.72 101 3.68 101

N = 8 1.003 104 2.913 102 6.53 102 0.58 103

N = 10 4.82 106 3.21 105 0.587 103 0.778 104

N = 15 3.92 1011 7.84 109 9.51 105 9.510 108

N = 20 6.22 1016 5.61 1014 1.51 1010 1.456 1012

N = 30 9.36 1019 8.22 1018 9.510 1011 9.108 1013

N = 40 9.28 1021 0.781 1018 5.87 1016 5.33 1018

To show the eciency of proposed method in solving nonlocal m-point boundary problem, we solve Example under a -point nonlocal boundary condition as dened in S. We observe that the method works well, the absolute dierence is much less than , a very high accuracy for such complicated problems. We compare the approximate solution with the exact solution at dierent scale levels. We also calculate E at dierent scale levels.

The results are displayed in Figure .

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 22 of 28

Figure 3 Simulation and observation of Example 1. (a) Comparison of approximate solution at dierent scale levels with the exact solution. (b) Absolute dierence of exact and approximate solution.(c) Convergence of E 2 at dierent scale levels.

Example As a second example consider the following fractional-order dierential equation:

D.U(t) + D.U(t) + D.U(t) + U(t) = g(t). ()

Here we consider g(t) as

g(t) = ,,,,,t

(t ,t + ,t ,) ,,,,, t(t )

,,,,,t

(,,t ,,t

+ ,,t ,,t + ,,) /,,,,,,,

,,,,,t

(,t ,,t + ,,t

,,t + ,,) /,,,,,,,. ()

We consider the following two types of boundary conditions:

S = !U() = , U(.) = ., U(.) = ., U() = ",

S = !U() = , U () = , U () = ,

U() + .U(.) + U(.) + .U(.) .U(.) = U()".

It can easily be observed that the exact solution of the problem is U(t) = t(t ).

We solve this problem with the proposed method under boundary conditions S, and we observe that the approximate solution obtained via the new method is very accurate even for very small value of N. Figure (a) shows the comparison of the approximate solution with the exact solution. One can easily see that the approximate solution matches very well with the exact solution. In the same gure we also display the absolute error obtained with the new method. It is clear that the absolute error is less than at scale level

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 23 of 28

Figure 4 Example 2. (a) Comparison of approximate solution with the exact solution under boundary condition S. (b) Absolute error obtained with the new method at dierent scale levels. (c) E 2 obtained

with the new method at dierent values of N.

Figure 5 Example 2. (a) Comparison of approximate solution with the exact solution under boundary condition S. (b) Absolute error obtained with the new method at dierent scale levels. (c) E 2 obtained

with the new method at dierent values of N.

N = . We measure E at dierent values of N and observe that the method is highly

convergent. We also approximate solution of this problem under boundary condition S.

The results are displayed in Figure . We can easily observe that the approximate solution of the method converges to the exact solution as the value of N increases.

Example As a third example consider the following fractional dierential equation with variable coecients:

D.U(t) + t t D. + t + D.U(t) + tU(t) = g(t). ()

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 24 of 28

We consider the following types of boundary conditions:

S = !U() = , U () = , U(.) + .U(.) + U(.) .U(.) = U()",

S = !U() = , U(.) = ., U() = .".

We select a suitable g(t), such that the exact solution of the problem is et. We simulate the proposed algorithm to solve this problem under boundary conditions S and we observe that the method works well. The results are displayed in Figure . In this gure we display the comparison of the exact and the approximate solution, the absolute error of approximation and the square norm of the error. One can easily conclude that the method is highly ecient. We solve the problem under boundary condition S and the results are displayed in Table . One can easily see that the approximate solution is much more accurate.

Figure 6 Example 3. (a) Comparison of the approximate solution with the exact solution under boundary condition S. (b) Absolute error obtained with the new method at dierent scale levels. (c) E 2 obtained

with the new method at dierent values of N.

Table 3 Comparison of exact and approximate solution of Example 3 using the boundary condition as dened in S

t N

Exact U(t) N = 2 N = 3 N = 4 N = 5 N = 6

t = 0.0 1.000000000 1.0152591045 0.9989268650 1.0000584609 0.9999973751 1.000000102 t = 0.1 1.105170918 1.1048778054 1.1055690936 1.1051373901 1.1051724714 1.105170872 t = 0.2 1.221402758 1.2115236915 1.2222416889 1.2213694032 1.2214034018 1.221402755 t = 0.3 1.3498588075 1.3351967629 1.3506284171 1.3498396744 1.3498594914 1.349858762 t = 0.4 1.4918246976 1.4758970196 1.4924130441 1.4918009400 1.4918267043 1.491824610 t = 0.5 1.6487212707 1.6336244616 1.6492793361 1.6486734983 1.6487243484 1.648721203 t = 0.6 1.8221188003 1.8083790888 1.8229110590 1.8220452095 1.8221217781 1.822118764 t = 0.7 2.0137527074 2.0001609012 2.0149919788 2.0136714962 2.0137550986 2.013752651 t = 0.8 2.2255409284 2.2089698990 2.2272058616 2.2254753425 2.2255438696 2.225540838 t = 0.9 2.4596031111 2.4348060820 2.4612364735 2.4595472949 2.4596078085 2.459603059 t = 1.0 2.7182818284 2.6776694503 2.7187675805 2.7181454617 2.7182834949 2.718281628

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 25 of 28

Figure 7 Example 2, E 2 vs. iteration, for Example 4, obtained with the new method at different

values of N. Here we use boundary conditions S1.

Example As a last example consider the following nonlinear fractional dierential equation:

D.U(t) + D.U(t) + U(t) = D.UU + U + D.U + g(t), ()

with boundary conditions

S = !U() = , U(.) = ., U() = .",

S = 'U() = , U () = ,

.U%

& + .U%

& + .U%

& .U%

& = U()(.

Select a suitable g(t) such that the exact and unique solution of the above problem is U(t) = et/. We approximate the solution of this problem with the iterative method proposed in the paper under boundary condition S. We carry out the iteration at dierent scale levels N. We observe that the method converges more rapidly to the exact solution for high values of N. For instance, at some values of N we calculate E at each iteration.

The results are displayed in Figure , one can see that E falls below at the fth

iteration using scale level N = . On solving the problem under boundary conditions S we observe that the method provides very accurate estimate of solution. The results are displayed in Table .

6 Conclusion

From experimental results and analysis of the proposed method we conclude that the method is ecient in solving linear and nonlinear fractional-order dierential equations under dierent types of boundary conditions. The method has the ability to solve both

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 26 of 28

Table 4 Comparison of exact and approximate solution of Example 4 using the boundary condition as dened in S

t N

Exact U(t) First iteration Second iteration Third iteration Fourth iteration

t = 0.0 1.0000000000 0.9993612934 1.0001381646 1.0000260271 1.000000020 t = 0.1 1.0168063300 1.1131408265 1.0322239559 1.0171510593 1.016806465 t = 0.2 1.0338951135 1.2060058045 1.0595571295 1.0344373486 1.033895324 t = 0.3 1.0512710963 1.2803801650 1.0835040499 1.0519291558 1.051271352 t = 0.4 1.0689391057 1.3388134647 1.1052061945 1.0696628704 1.068939387 t = 0.5 1.0869040495 1.3839808784 1.1255801546 1.0876670104 1.086904346 t = 0.6 1.1051709180 1.4186831997 1.1453176345 1.1059622221 1.105171226 t = 0.7 1.1237447856 1.4458468410 1.1648854519 1.1245612803 1.123745105 t = 0.8 1.1426308117 1.4685238330 1.1845255377 1.1434690886 1.142631142 t = 0.9 1.1618342427 1.4898918253 1.2042549363 1.1626826787 1.161834578 t = 1.0 1.1813604128 1.5132540859 1.2238658054 1.1821912111 1.181360731

local and nonlocal boundary value problems. We use normalized Bernstein polynomials for our analysis. But the method can be used to generalize such types of operational matrices for almost all types of orthogonal polynomials. It is also possible to get a more approximate solution of such problems using other types of orthogonal polynomials like Legendre, Jacobi, Laguerre, Hermite etc. It is not clear to us which is the best set of orthogonal polynomials for this method. Further investigation is required to generalize the method to solve other types of scientic problems.

Competing interests

The authors declare that they have no competing interests.

Authors contributions

The authors have made equal contributions in preparing the proofs and numerical simulations. The authors have approved the nal manuscript.

Author details

1Department of Mathematics, University of Poonch Rawalakot, Rawalakot, 12350, Pakistan. 2Department of Mathematics, University of Malakand, P.O. Box 18000, Chakdara, Dir Lower, Khybarpukhtunkhwa, Pakistan. 3Dean Faculty of Science, University of Malakand, Chakdara, Dir Lower, Khybarpukhtunkhwa, Pakistan. 4Department of Mathematics and Computer Science, Cankaya University, Ankara, Turkey. 5Department of Mathematics, Mansoura University, Al Mansurah, Muhafazat ad Daqahliyah, Egypt.

Acknowledgements

The authors are thankful to the reviewers and the editor for carefully reading and useful suggestion which improved the quality of the article.

Received: 11 April 2016 Accepted: 24 June 2016

References

1. Baleanu, D, Diethelm, K, Scalas, E, Trujillo, JJ: Fractional Calculus Models and Numerical Methods, Series on Complexity, Nonlinearity and Chaos. World Scientic, Boston (2012)

2. Kilbas, AA, Srivastava, HM, Trujillo, JJ: Theory and Applications of Fractional Dierential Equations. North-Holland Mathematics Studies, vol. 204. Elsevier, Amsterdam (2006)

3. Lazarevic, MP, Spasic, AM: Finite-time stability analysis of fractional order time-delay systems: Gronwalls approach. Math. Comput. Model. 49(3), 475-481 (2009)

4. Magin, RL: Fractional Calculus in Bioengineering. Begell House, Redding (2006)5. Podlubny, I: Fractional Dierential Equations. Academic Press, San Diego (1999)6. Sabatier, J, Agrawal, OP, Machado, JAT (eds.): Advances in Fractional Calculus: Theoretical Developments and Applications in Physics and Engineering. Springer, Dordrecht (2007)

7. Zaslavsky, GM: Hamiltonian Chaos and Fractional Dynamics. Oxford University Press, Oxford (2008)8. Gupta, CP: Solvability of a three-point nonlinear boundary value problem for a second order ordinary dierential equation. J. Math. Anal. Appl. 168, 540-551 (1992)

9. Ma, R: A survey on nonlocal boundary value problems. Appl. Math. E-Notes 7, 257-279 (2007)10. Guezane-Lakoud, A, Zenkou, L: Existence of positive solutions for a third-order multi-point boundary value problem. Appl. Math. 3, 1008-1013 (2012)

11. El Sayed, AMA, Bin-Tahir, EO: Positive solutions for a nonlocal boundary-value problem of a class of arbitrary (fractional) orders dierential equations. Int. J. Nonlinear Sci. Numer. Simul. 14(4), 398-404 (2012)

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 27 of 28

12. Canuto, C, Hussaini, MY, Quarteroni, A, Zang, TA: Spectral Methods in Fluid Dynamics. Springer, New York (1988)13. Bhrawy, AH, Alo, AS: A Jacobi-Gauss collocation method for solving nonlinear Lane-Emden type equations. Commun. Nonlinear Sci. Numer. Simul. 17, 62-70 (2012)

14. Funaro, D: Polynomial Approximation of Dierential Equations. Lecturer Notes in Physics. Springer, Berlin (1992)15. Doha, EH, Bhrawy, AH, Ezz-Eldien, SS: Ecient Chebyshev spectral methods for solving multi-term fractional orders dierential equations. Appl. Math. Model. 35, 5662-5672 (2011)

16. Bhrawy, AH, Alo, AS, Ezz-Eldien, SS: A quadrature tau method for variable coecients fractional dierential equations. Appl. Math. Lett. 24, 2146-2152 (2011)

17. Saadatmandi, A, Dehghan, M: A new operational matrix for solving fractional-order dierential equations. Comput. Math. Appl. 59, 1326-1336 (2010)

18. Doha, EH, Bhrawy, AH, Ezz-Eldien, SS: A Chebyshev spectral method based on operational matrix for initial and boundary value problems of fractional order. Comput. Math. Appl. 62, 2364-2373 (2011)

19. Ghoreishi, F, Yazdani, S: An extension of the spectral Tau method for numerical solution of multi-order fractional dierential equations with convergence analysis. Comput. Math. Appl. 61, 30-43 (2011)

20. Vanani, SK, Aminataei, A: A Tau approximate solution of fractional partial dierential equations. Comput. Math. Appl. 62, 1075-1083 (2011)

21. Esmaeili, S, Shamsi, M: A pseudo-spectral scheme for the approximate solution of a family of fractional dierential equations. Commun. Nonlinear Sci. Numer. Simul. 16, 3646-3654 (2011)

22. Pedas, AA, Tamme, E: On the convergence of spline collocation methods for solving fractional dierential equations.J. Comput. Appl. Math. 235, 3502-3514 (2011)23. Bhrawy, AH, Al-Shomrani, MM: A shifted Legendre spectral method for fractional-order multi-point boundary value problems. Adv. Dier. Equ. 2012 8 (2012)

24. Khalil, H, Khan, RA: The use of Jacobi polynomials in the numerical solution of coupled system of fractional dierential equations. Int. J. Comput. Math. (2014). doi:10.1080/00207160.2014.945919

25. Shah, K, Ali, A, Khan, RA: Numerical solutions of fractional order system of Bagley-Torvik equation using operational matrices. Sindh Univ. Res. J. (Sci. Ser.) 47(4), 757-762 (2015)

26. Khalil, H, Khan, RA: A new method based on Legendre polynomials for solutions of the fractional two-dimensional heat conduction equation. Comput. Math. Appl. (2014). doi:http://dx.doi.org/10.1016/j.camwa.2014.03.008

Web End =10.1016/j.camwa.2014.03.008

27. Khalil, H, Khan, RA: New operational matrix of integration and coupled system of Fredholm integral equations. Chin. J. Math. 2014, Article ID 146013 (2014).

28. Khalil, H, Khan, RA: A new method based on Legendre polynomials for solution of system of fractional order partial dierential equation. Int. J. Comput. Math. 91, 2554-2567 (2014). doi:10.1080/00207160.2014.880781

29. Khalil, H, Khan, RA, Al Smadi, MH, Freihat, A: Approximation of solution of time fractional order three-dimensional heat conduction problems with Jacobi polynomials. J. Math. 47(1), 35-56 (2015)

30. Khalil, H, Rashidi, MM, Khan, RA: Application of fractional order Legendre polynomials: a new procedure for solution of linear and nonlinear fractional dierential equations under m-point nonlocal boundary conditions. Commun. Numer. Anal. 2016(2), 144-166 (2016)

31. Khalil, H, Khan, RA, Baleanu, D, Rashidi, MM: Some new operational matrices and its application to fractional order Poisson equations with integral type boundary constrains. Comput. Appl. Math. (2016). doi:http://dx.doi.org/10.1016/j.camwa.2016.04.014

Web End =10.1016/j.camwa.2016.04.014

32. Bhrawy, AH, Zaky, MA: A method based on the Jacobi tau approximation for solving multi-term time-space fractional partial dierential equations. J. Comp. Physiol. 281, 876-895 (2015)

33. Bhrawy, AH, Zaky, MA: Numerical simulation for two-dimensional variable-order fractional nonlinear cable equation. Nonlinear Dyn. 80(1-2), 101-116 (2015)

34. Bhrawy, AH, Zaky, MA: New numerical approximations for space-time fractional Burgers equations via a Legendre spectral-collocation method. Rom. Rep. Phys. 67(2), 340-349 (2015)

35. Zaky, MA, Bhrawy, AH, Van Gorder, RA: A space-time Legendre spectral tau method for the two-sided space-time Caputo fractional diusion-wave equation. Numer. Algorithms 71(1), 151-180 (2016)

36. Bhrawy, AH, Zaky, M: Shifted fractional-order Jacobi orthogonal functions: application to a system of fractional dierential equations. Appl. Math. Model. 40(2), 832-845 (2016)

37. Rehman, M, Khan, RA: A numerical method for solving boundary value problems for fractional dierential equations. Appl. Math. Model. 36, 894-907 (2012)

38. Rehman, M, Khan, RA: The Legendre wavelet method for solving fractional dierential equations. Commun. Nonlinear Sci. Numer. Simul. 16, 4163-4173 (2011)

39. Liu, Y: Numerical solution of the heat equation with nonlocal boundary condition. J. Comput. Appl. Math. 110, 115-127 (1999)

40. Ang, W: A method of solution for the one-dimensional heat equation subject to nonlocal condition. Southeast Asian Bull. Math. 26, 185-191 (2002)

41. Dehghan, M: The one-dimensional heat equation subject to a boundary integral specication. Chaos Solitons Fractals 32, 661-675 (2007)

42. Noye, KHBJ: Explicit two-level nite dierence methods for the two-dimensional diusion equation. Int. J. Comput. Math. 42, 223-236 (1992)

43. Avalishvili, G, Avalishvili, M, Gordeziani, D: On integral nonlocal boundary value problems for some partial dierential equations. Bull. Georgian. Natl. Acad. Sci. (N. S.) 5, 31-37 (2011)

44. Sajavicius, S: Optimization, conditioning and accuracy of radial basis function method for partial dierential equations with nonlocal boundary conditions, a case of two-dimensional Poisson equation. Eng. Anal. Bound. Elem. 37, 788-804 (2013)

45. Youse, SA, Behroozifar, M: Operational matrices of Bernstein polynomials and their applications. Int. J. Inf. Syst. Sci. 41(6), 709-716 (2010)

46. Doha, EH, Bhrawy, AH, Saker, MA: On the derivatives of Bernstein polynomials: an application for the solution of high even-order dierential equations. Bound. Value Probl. 2011, 829543 (2011)

47. Doha, EH, Bhrawy, AH, Saker, MA: Integrals of Bernstein polynomials: an application for the solution of high even-order dierential equations. Appl. Math. Lett. 24(1), 559-565 (2011)

Khalil et al. Advances in Dierence Equations (2016) 2016:177 Page 28 of 28

48. Juttler, B: The dual basis functions for the Bernstein polynomials. Adv. Comput. Math. 8(4), 345-352 (1998)49. Farouki, RT: Legendre Bernstein basis transformations. J. Comput. Math. 119(1), 145-160 (2000)50. Hermann, T: On the stability of polynomial transformations between Taylor, Bernstein, and Hermite forms. Numer. Algorithms 13(2), 307-320 (1996)

51. Bellucci, MA: On the explicit representation of orthonormal Bernstein polynomials. http://arxiv.org/abs/1404.2293v2

Web End =http://arxiv.org/abs/1404.2293v2 52. Chen, Y, Sun, Y, Liu, L: Numerical solution of fractional partial dierential equations with variable coecients using generalized fractional-order Legendre functions. Appl. Comput. Math. 244, 847-858 (2014)

53. Knuth, DE: The Art of Computer Programming. Fundamental Algorithms, vol. 1. Addison-Wesley, Reading (1968)54. https://proofwiki.org/wiki/Inverse_of_Vandermonde_Matrix

Web End =https://proofwiki.org/wiki/Inverse_of_Vandermonde_Matrix 55. Bellman, RE, Kalaba, RE: Quasilinearization and Non-linear Boundary Value Problems. Elsevier, New York (1965)56. Stanley, EL: Quasilinearization and Invariant Imbedding. Academic Press, New York (1968)57. Agarwal, RP, Chow, YM: Iterative methods for a fourth order boundary value problem. J. Comput. Appl. Math. 10(2), 203-217 (1984)

58. Akyuz Dascioglu, A, Isler, N: Bernstein collocation method for solving nonlinear dierential equations. Math. Comput. Appl. 18(3), 293-300 (2013)

59. Charles, A, Baird, J: Modied quasilinearization technique for the solution of boundary-value problems for ordinary dierential equations. J. Optim. Theory Appl. 3(4), 227-242 (1969)

60. Mandelzweig, VB, Tabakin, F: Quasilinearization approach to nonlinear problems in physics with application to nonlinear ODEs. Commun. Comput. Phys. 141(2), 268-281 (2001)