А fractional analog of Crank-Nicholson method for the two sided space fractional partial equation with functional delay Текст научной статьи по специальности «Математика»

URAL MATHEMATICAL JOURNAL, Vol. 2, No. 1, 2016

A FRACTIONAL ANALOG OF CRANK-NICHOLSON METHOD FOR THE TWO SIDED SPACE FRACTIONAL PARTIAL EQUATION WITH FUNCTIONAL DELAY

Vladimir G. Pimenov

Department of Computational Mathematics, Ural Federal University;

Krasovskii Institute of Mathematics and Mechanics, Ural Branch of the Russian Academy of Sciences, Ekaterinburg, Russia,

v.g.pimenov@urfu.ru

Ahmed S. Hendy

Department of Computational Mathematics, Ural Federal University,

Ekaterinburg, Russia, ahmed.hendy@fsc.bu.edu.eg

Abstract: For two sided space fractional diffusion equation with time functional after-effect, an implicit numerical method is constructed and the order of its convergence is obtained. The method is a fractional analogue of the Crank—Nicholson method, and also uses interpolation and extrapolation of the prehistory of model with respect to time.

Key words: Fractional partial differential equation, Grunwald—Letnikov approximations, Grid schemes, Functional delay, Interpolation, Extrapolation, Convergence order.

Introduction

Diffusion equations with delay of general type, constant or variable, concentrated or distributed, are often applied in the simulation of dynamics processes. There are two effects often combined with these equations: distribution of parameters in space and heredity in time [1,2]. Numerical methods for solving such equations were considered in many papers, for example [3-7]. In paper [8], a technique of study on stability and convergence of numerical algorithms using the general theory of differential schemes was constructed for the heat conduction equation with delay [9] and the theory of numerical methods of the solution of the functional and differential equations were discussed in [10,11]. After that, this technique was applied to research of numerical methods of the solution of the equations of hyperbolic type with delay [12], various equations of parabolic type [13,14] and other types of the equations in partial derivatives with effect of heredity [15]. Here this technique is applied to the equations with partial derivatives of a fractional order effected with functional delay.

Fractional differential equations [16-18] gained a great interest in the past decades due to their accuracy in modeling a lot of problems in many fields of science. The equations in partial derivatives of a fractional order were subdivided into two big classes: with a fractional derivative on space and with a fractional derivative on time. Now there are many papers on numerical methods of the solution of such equations [19-21]. In this paper, we depend on results of work [22] in which explicit and purely implicit numerical methods for the solution of the equation with two sided space fractional derivative were investigated. In paper [23], an analog of Krank-Nikolson method was investigated for the equation which has a left-sided fractional derivative on space. In our paper, an analog of Krank-Nikolson method for the solution of the equation with two sided space fractional derivative and with functional time-delay is constructed and investigated. We need to extend this implicit scheme to work on the functional delayed form of these equations. The main purpose of

this approach is to test the possibility of embedding the constructed implicit difference scheme for the fractional diffusion equation with functional delay into a general difference scheme.

After that, the convergence order is deduced by using the technique of proving similar statements for functional differential equations and methods from the general theory of difference schemes. We note that in work [24] this technique was announced for the numerical method of the solution of the equation from a left-side fractional derivative on space and with delay on time.

We consider the following initial-boundary value problem with functional delay of fractional order diffusion equation 1 < a ^ 2 :

Qu d au dau

~dt = c+(x) d+Xa + c-(x) + f (x,t,u(x,t),u (x •))> (0A)

where x £ [0,X], t £ [t0,T], c+(x) > 0, c-(x) > 0, ut(x, •) = {u(x,t + s), —t ^ s < 0} is the prehistory function, t > 0 is the value of delay, with initial conditions:

u(x,t) = p(x,t), x £ [0,X], t £ [to — t, to],

and boundary conditions:

u(0, t) = 0, u(X, t) = 0, t £ [to, T]. The left-sides and the right-sides fractional derivatives are defined in Riemann-Liouville sense [17]

x

da u(x) 1 d2 f u(£)

d+xa r(2 — a) dx2 J (x — £)'

o

X

da u(x) =1 d2 f u(£) d-xa = r(2 — a) dx2 J (£ — x)a— 1 ^

x

We assume that the functions p(x,t), c+(x), c-(x) and the functional f are such that this problem has unique solution u(x,t).

We denote by Q = Q[—t, 0) the set of functions u(s) that are piecewise continuous on [—t, 0) with a finite number of points of discontinuity of the first kind and right continuous at the points of discontinuity. We define the norm of functions by the relation

II u(0 IIq= suP lu(s)|.

se[-r,o)

We additionally assume that the functional f (x,t,u,v()) is given on [0,X] x [t0,T] x R x Q and is Lipschitz in the last two arguments, i. e., there exists a constant Lf such that, for all x £ [0,X], t £ [t0,T], u1 £ R, u2 £ R, v*(•) £ Q and v2(^) £ Q the following inequality holds:

lf(x,t,u1,v1(•)) — f(x,t,u2,v2(»l < Lf {lu1 — u2|+ II v1« — v2(0 IIq ). (0.2)

We will use the approximations based on the shifted Griinwald definitions for the left-sides and right-sides fractional derivatives [22] respectively

dau(x) 1 , ,, . ,

^^ = (Ax+F L u(x—(k—1)Ax),

N-

da u (x) 1 _>

^^ = (Ax-a S 9a'k u(x+(k+1)Ax), V ' k=0

where N+, N- are positive integers, Ax+ = (x)/N+, ax- = (X — x)/N- and the normalized Grunwald weights are given by

. 1)k a(a — 1)...(a — k + 1) 9a, 0 = 1, 9a,k = ( —1) -k!-

for k = 1,2,3,....

1. Derivation of the fractional Crank—Nicholson method

Let h = X/N, introduce Xi = ih, i = 0,...,N, and A = (T - to)/M, tj = to + j a, j = 0,...,M, we assume that т/a = m is a positive integer. Denote by uj approximations of functions u(xi,tj) at the nodes.

Let us introduce a discrete prehistory for the time points tj, j = 0,... ,M : {uik}j = {uik ,j — m ^ k ^ j}. The mapping I : {uik}j — vi(t), t £ [tj — T,tj + A/2] will be called the interpolationextrapolation operator of discrete prehistory.

We will use piecewise-linear interpolation with extrapolation by continuation

( i(ui(t — t-i)+ ui-x(tl — t)), tl-1 < t < ti, 1 < l < j, vi(t) = I i(uj(t — tj-i) + uj-i(tj — t)), tj < t < tj + A/2,

[ p(Xi,t), to — T < t < to.

With the shifted Griinwald estimates for left and right fractional derivatives on space, applying usual approximation to a derivative on time and by the use of interpolation which is designed by extrapolation to the prehistory function, the discretized (0.1) at the nodes (xi,tj+1/2) will give us the analog of Crank-Nicholson method

ui -ui 1 / i+1 N-i+1 X

= gasj+1 + c- E МЙ-1) + j 2 , (1.1)

s=0 s=0

/;+! = f(xi,tj + | ,vi (tj + (•)), i = 1,..., N - 1, j = 0,...,M - 1,

C+ = C+(Xi), С_ = C-(Xi), uj+Tz/2 = 2(uj +

with the initial and boundary conditions

u0 = ^(xi,t0), i = 0,...,N, vj (t) = p(xi,t), t<to, i = 0,...,N,

u° = 0, uf = 0, j = 0,...,M.

Define

С i+! . ci N _i+1

$a,xuj : = h+a ga,su j + , $_,xu j : = J— 9a,su j+ ,

s=0 s=0

then (1.1) we can be written as follows

uj+i uj 1 / +i+i i i \ i

= 2 \^a,xuj + $a,xuj+1 + $_,xuj + $_,xuj+1J + fjj+2 . (1.2)

A

If a = 2, c(x) = c+(x)+c_(x), we receive the Crank-Nicholson scheme for the variable coefficient heat conductivity equation with functional delay, because

+ . . uj-1 - 2uj + uj+1

Ô2)Xu] = ô+uj + S- u = j j j

ruj = $a,xuj + $a,xuj J2

2. Stability and convergence of the method

Denote by ej = u(xi,tj) — uj the error of the method at the nodes. We will say that the method converges with order hp + aq if there exists a constant C independent of h and A such that lejl < C(hp + Aq) for all i = 0,1, • • • ,N and j = 0,1, • • • ,M.

The residual (without interpolation) of the method (1.2) is, by definition, the grid function

i = u(xi ,tj+1) — u(xi ,tj+1) Vj a

1

-Û,xU(Xi,tj) + Û,xU{Xi,tj+l) + Ôa,xU(Xi,tj) + Ôa,xU(Xi, j+l}) - f+ 1 , (2.1)

. i i+1 „ ci N-i+1 Ô+,xU(Xi ,tj ) = -J+^2 ga ,su(Xi-s+1,tj ), Ô-,xU(Xi,tj ) = h- ^ 9a Su(Xi+s-1 ,tj ),

s=0 s=0

fi+1 = f(xi,tj+1 ,u(xi,tj+1),ut.+ i (xi, •)).

Lemma 1. Suppose that the exact solution u(x,t) is twice continuously differentiable in t and

dau dau(x t)

four times continuously differentiable in x and also the fractional derivatives —- and —-—-—

d+xa d-xa

are twice continuously differentiable in x. Then, the residual (without interpolation) of the method has order a2 + h.

Proof. Based on Taylor series expansions, we have

u(xi,tj+1) — u(xi,tj+1) = du(xi,tj+ 1) A dt

Following [22], we find that 5+,xu(xi,tj) = ) + O2, S+xu(xt,tj+1) = dau(xj) + °3, lO2l< C2h, lO3l< Ch

+ O1, |01|< C1a2. (2.2)

and similarly

dau(Xi,tj ) - ( ) = dau(Xi,tj+1)

Sa ,xU(Xi,tj ) = -dX+ + Sa,xU(Xi,tj+1) = -d^-'--+ IO4I < C^k, |O5|< C5h,

then,

dau(Xi,tj+1 )

5+,xu(Xi,t3 ) + S+,xu(Xi,tj+1) = 2- aJ+2 + O2 + O3 + Ob, |O6| < Co A2, (2.3)

()X+

- - dau(Xi,tj+1 ) 5-,xu(Xi, tj) + 5-,xu(Xi, tj+1) = 2-J 2 + O4 + O5 + O7, |O7| < C7A2. (2.4)

Substituting (2.2), (2.3) and (2.4) in the definition of residual (without interpolation) (2.1) and consider that u(Xi,tj+1 ) is the exact solution of the equation (0.1), the statement of the lemma is received. □

For every fixed i = 0,..., N let's introduce the discrete prehistory of exact solution for the time points tj, j = 0,... ,M : {u(Xi,tk)}j = {u(Xi,tk),j — m ^ k ^ j}. We will use piecewise-linear interpolation with extrapolation by continuation of exact solution

- (u(Xi,ti)(t — ti-1)+ u(Xi,ti-1)(ti — t)), ti-1 < t < tl, 1 < l < j, W(t) = { -(u(Xi, tj)(t — tj-1) + u(Xi, tj-1)(tj — t)), tj < t < tj + a/2,

Xi,t), to — T < t < to.

The residual (with interpolation) of methods (1.2) can be defined using the following grid function

u(Xi,tj+l) — u(Xi,tj+l)

v. = .

jA

-S+,xu(xi,tj) + S+,xu(xi,tj+i) + S-,xu(xi,tj) + Sa,xu(xi,tj+i)^ - j2, (2.5)

f+ 1 = f (x.,tj+2 ,wi(tj + A),wj + f (.)).

Lemma 2. Assume that the conditions of Lemma 1 are satisfied, then the residual of the method with piecewise-linear interpolation and extrapolation by continuation has order A2 + h.

Proof. Since the piecewise linear interpolation operator with extrapolation by continuation

j + "2 > tj

has second order [8], therefore there exists a constant C, such that for all t & [tj + 4,tj — t]

\wi(t) — u(xi,t)\< CA2. (2.6)

The residual with interpolation (2.5) and the residual without interpolation (2.1) are connected with the ratio

vj = ^j + f (Xi,tj+i,u(Xi,-tj+2 W (Xi, •)) — f (Xi,tj+2,wi(tj + A),wl+a (0). (2.7) 2 2 2 2 j 2

From (0.2), (2.6), (2.7) and from Lemma 1 we get that the proof is completed. □

We can copy the method (1.2) as follows

(1 - T(S+>a: + S-J)uj+i = (1 + aô+x + ô~x))uj + Afi+1. (2.8)

In order to reduce it to the general scheme [8], introduce yj = (u], u2, • • • ,uN & Y, such that Y is a vector space of dimension N — 1 with norm

II yj ||y = max \uj\. 11 yj ly 1<i<N-11 j

Then (2.8) can be rewritten as

(E + A)yj+i = (E — A)yj + aFj+ 2, (2.9)

where elements of a matrix A with dimension N — 1 x N — 1 have the form

— (£i + Vi)9a,1, j = i,

— (tea,2 + Vi9a,0), j = i — 1,

Aij = < —((i9a,0 + Vi9a,2), j = i + 1,

^igai-j+U j<i — 1,

, ^igaj-i+U j>i + 1,

E is unit matrix, £i = "2+a", ni = 2ha , Fj+2 = (fj+1 , fj+1, • • • , fj+l ).

Lemma 3. [18,25] The coefficients ga,s satisfy the following properties

ga,o = 1, ga,i = -a, 1 > ga,2 > ga,3 > ... > 0,

œ m

^2,ga,s = 0, ^2,ga,s < 1, m > 1.

s=0 s=0

Lemma 4. The matrix E + A is reversible.

Proof. From the values of the elements of matrix A and from Lemma 3, we conclude that

N-1 i N-i

a* - E I aj 1= + ni)9a,1 - 6 £ ga,s - * E 9a,s >

j=1,i=j s=0,s=1 s=0,s=1

<x

>-(ti + Vi)ga,1 - (ti + Vi) E 9a,s = + ni)9a,1 + (ti + Vi)9a,1 = 0. (2.10)

s=0,s=1

This inequality means that the coefficient matrix E + A is a nonsingular (reversible), strictly diagonally dominant matrix.

We will note that if a is not integer number, then inequality (2.10) is strict and A is strictly diagonally dominant matrix. □

As the matrix of E + A is reversible, equation (2.8) can be copied as follows

Vj+1 = Syj + A$(tj ,I ({yk }j)), (2.11)

S = (E + A)-1(E - A), $(j,I({yk}j)) = (E + A)-1 Fj+2, I is the operator of piecewise-linear interpolation with extrapolation by continuation.

Lemma 5. If 1 < a < 2, then any eigenvalues X of matrix S = (E + A)-1(E - A) satisfies the condition | X | < 1.

P r o o f. Let X be an eigenvalue of matrix A. As a is not integer number, A is strictly diagonally dominant nonsingular matrix. Since ga,1 = -a and ti + n gives a positive value, then all diagonal elements of a matrix A are positive. According Gershgorin circle theorem [26, c. 135], we conclude that the real parts of its eigenvalues are all positive. This means that

I 1 + X |>| 1 - X |, | 1-X I< 1,

then IXI < 1. □

From this lemma, it follows there is such constant S, that

||Sra|| < S (2.12)

for any natural degree n. This means that the method (2.11) is (unconditionally) stable.

Now, the theorem of a convergence order [8] which is the main result of the general theory of numerical methods of the solution of evolutionary equations with heredity can be applied to the scheme (2.11). As the resulted scheme (2.11) has some different features in comparison with the scheme in [8], in particular, another determination of stability (2.12) is used in comparison with [8] then proof of the following theorem is needed.

Theorem 1. Assume that the condition in Lemma 3 is satisfied, the function $ satisfies Lip-schitz condition in the second argument, interpolation operator I satisfies Lipschitz condition, the residual with interpolation has order Aq+hp. Then, the method (2.11) converges with order Aq+hp.

Proof. Introduce the following vector of exact values

zn = (u(x1,tn) ,u(x2,tn), ■ ■ ■ , u(xN-1, tn)) and note that yn = zn - yn, n = -m,..., M. Then for n = 0,..., M - 1 we have

Yn+1 = SYn + A7n + Avn, (2.13)

where

Sn = $(tn, I({Zi}n)) - $(tn, I({yi}n)).

Here

Vn = (Zn+1 - Szn)/A - $(tn, I({Zi}n)), n = 0,..., M - 1

is the vector residual with interpolation vn = (vln, vn, ■ ■ ■ , v.N-n). The assumptions that mapping $ and I are Lipschitz implies that

||Sn|| < L ■ max {^H}, (2.14)

where L = L$LJ, L$ is Lipschitz constant of $ on second argument, Lj is Lipschitz constant of operator I.

We notice that in the case of the described method L$ = ||(E - A)-1HLf, The piecewise linear interpolation operator with extrapolation by continuation has Lipschitz constant Lj = 2. Also, we note that y0 = 0, as starting values of a method are the exact values, unlike the general scheme [8]. It follows from (2.13) that

nn

Yn+1 = Sn+17o + A E Sn-j Sj + A E Sn-j vj. (2.15)

j=0 j=0

From (2.15), (2.14) and from the definition of stability, we have

n

||Yn+1|| < SlA^ , max {||Yi||} + Stmax wiH}. (2.16)

j=0

We use the notation

R = max {HviH}, D = STR. (2.17)

0<i<N-1

Then we can write the estimate (2.16) in the form

n

pn+1|| < SLAV max {||Si||} + D. (2.18)

j-m^i^j

j=0

Depending on (2.18) and using the mathematical induction, let us prove the estimate

||Sn|| < D(1 + SLa)n, n = 1,..,M. (2.19)

Induction base. If we set n = 0, in (2.18), then

< slUsoU + D < (1 + sla)d.

Induction step. Let the estimate (2.19) be valid for all indices from 1 to n. We need to show that the estimate is also valid for n + 1. Fix j ^ n. Let i0 = i0(j) be an index for which {pj||}

is obtained. By induction assumption, we gain

max {||vi||} = IVi01| < D(1 + SLA)i0 < D(1 + sLv)j.

j

Thus, the following estimate is also valid

max {||vi||} < D(1 + SLv)j.

j

Using the previous inequality and (2.18), we have

|K+i|| < SLA E D(1 + SLa)j + D = D(1 + SLA)n+1. j=0

Thus, the estimate (2.19) is proved and this gives

IIvnII < D exp(SLT). (2.20)

Recall the notation of the value D from (2.17), then the inequality

D < C (Api + hp2),

holds and with the aid of (2.20) the proof is completed. □

Corollary 1. Under conditions of Lemmas 1 and 3, the method (1.2) with piecewise-linear interpolation and extrapolation by continuation converges with order A2 + h.

3. Numerical example

Let us consider the following two sided space fractional equation with varying delay in the time

^ = (X — 1/2)a + <3/2 — X)a ^ + (3.1)

f = ln(exp(—i/2)(x'- 1/2)3(3/2 — x)3) ((x 1/2)3(3/2 X'' +

exp(—t)

'(x — 1/2

+f(^(X —1/2)3(3/2 — X)3 — rgF——l)(X —1/2)4(3/2 — X)4+

+r3i"=6ik)(X — 1/2)5(3/2 — X)5 — rT—^1(X — 1/2)6(3/2 — x)6) ln(u(X, t — t/2)),

such that 1/2 < x < 3/2, 1 < t < 5, a is a constant, such that 1 < a < 2, with initial and boundary conditions on the form

u(x, r) = exp(—r)(x — 1/2)3(3/2 — x)3, 1/2 < r < 1, 1/2 < x < 3/2,

u(1/2, t) = 0, u(3/2, t) = 0, 1 < t < 5.

The exact solution of (3.1) is u(x,t) = exp(—t)(x — 1/2)3(3/2 — x)3. Define the maximum error by

E(A,h) = max \ u(xi,tj) — uj \,

0<i<N J

0<j<M

such that u(xi,tj) is the exact solution at the grid point (xi,tj). Let us test the spatial errors and their convergence orders by letting h vary from 1/20 to 1/320, fix the time step A = 4/1024. The convergence order with respect to the spatial step size is given by orders = log2 ^E^Ah)) which is clarified in table 1. For studying of temporal errors, let A varies from 1/16 to 1/256 and fix h = 1/4000. The convergence order with respect to the time step size is given by ordert = log2(Ejli?) and this is illustrated in table 2. From these tables, we can notice that the numerical results have a good agreement with the theory results.

h a = 1.1 a = 1.9

E(A, h) orders E(A, h) orders

1/20 0.00497 0.0027

1/40 0.00252 0.9752 0.0014 0.9865

1/80 0.00127 0.9887 0.0007 0.9984

1/160 0.00064 0.9996 0.0003 0.9997

1/320 0.00032 0.9999 0.0002 1.0002

Table 1. The spatial maximum norm errors and their convergence orders.

A a = 1.1 a = 1.9

E(A, h) ordert E(A, h) ordert

1/16 0.0054 0.0022

1/32 0.0014 1.979 0.0005 1.988

1/64 0.0003 1.985 0.0001 1.993

1/128 0.00008 1.996 0.00003 1.999

1/256 0.000002 1.999 8.7 x 10"6 2.001

Table 2. The temporal maximum norm errors and their convergence orders.

Conclusion

A fractional analog of Crank Nicholson method is constructed to deal with the one dimensional space two sided space fractional diffusion equation with functional delay. The method is based on the idea of separating the current state and the prehistory function. To deal with the prehistory function, we introduced a discrete prehistory for the time points and the piecewise-linear interpolation with extrapolation by continuation operator is used. Unconditional stability of the proposed scheme is achieved. A theorem is obtained on the order of convergence of the method, which used the technique of proving similar statements for functional differential equations and methods from the general theory of difference schemes. A numerical example with time varying delay is introduced to test the agreement between theoretical and numerical results.

Acknowledgements

This work was supported by Government of the Russian Federation program 02.A03.21.0006 on 27.08.2013 and by Russian Science Foundation 14-35-00005.

REFERENCES

1. Wu J. Theory and applications of partial functional differential equations. New York: Springer-Verlag, 1996. 438 p.

2. Zhang B., Zhou Y. Qualitative analysis of delay partial difference equations. New York: Hindawi Publishing Corporation, 2007. 375 p.

3. Tavernini L. Finite difference approximations for a class of semilinear Volterra evolution problems // SIAM J. Numer. Anal., 1977. Vol 14, no. 5. P. 931-949.

4. Van Der Houwen P.J., Sommeijer B.P., Baker C.T.H. On the stability of predictor-corrector methods for parabolic equations with delay // IMA J. Numer. Anal., 1986. Vol. 6. P. 1-23.

5. Zubik-Kowal B. The method of lines for parabolic differential-functional equations // IMA J. Numer. Anal., 1997. Vol 17. P. 103-123.

6. Kropielnicka K. Convergence of Implicit Difference Methods for Parabolic Functional Differential Equations // Int. Journal of Mat. Analysis, 2007. Vol. 1, no. 6. P. 257-277.

7. Garcia P., Castro M.A., Martin J.A., Sirvent A. Convergence of two implicit numerical schemes for diffusion mathematical models with delay // Mathematical and Computer Modelling. 2010. Vol. 52. P. 1279-1287.

8. Pimenov V.G., Lozhnikov A.B. Difference schemes for the numerical solution of the heat conduction equation with aftereffect // Proc. Steklov Inst. Math., 2011. Vol. 275. Suppl. 1. P. 137-148.

9. Samarskii A.A. Theory of difference schemes. Moscow: Nauka, 1989. 656 p. [In Russian]

10. Pimenov V.G. General linear methods for the numerical solution of functional-differential equations // Differential Equations, 2001, Vol. 37, no. 1. P. 116-127.

11. Kim A.V., Pimenov V.G. i-smooth calculus and numerical methods for functional differential equations. Moscow-Izhevsk: Regular and Chaotic Dynamics, 2004. 256 p. [In Russian]

12. Pimenov V.G., Tashirova E.E. Numerical methods for solving a hereditary equation of hyperbolic type // Proc. Steklov Inst. Math, 2013. Vol. 281. Suppl. 1. P. 126-136.

13. Lekomtsev A.V., Pimenov V.G. Convergence of the Alternating Direction Methods for the Numerical solution of a Heat Conduction Equation with Delay // Proc. Steklov Inst. Math., 2011. Vol. 272. Suppl. 1. P. 101-118.

14. Lekomtsev A, Pimenov V. Convergence of the scheme with weights for the numerical solution of a heat conduction equation with delay for the case of variable coefficient of heat conductivity // Appl.Math.Comput., 2015. Vol. 256. P. 83-93.

15. Pimenov V.G. Difference methods for the solution of the equations in partial derivatives with heredity. Ekaterinburg: Ural Federal University, 2014. 232 p. [In Russian]

16. Samko S.G., Kilbas A.A., Marichev O.I. Fractional Integrals and Derivatives: Theory and Applications. Boca Raton: CRC Press, 1993. 1006 p.

17. Miller K, Ross B. An Introduction to the Fractional Calculus and Fractional Differential Equations. New York: Wiley, 1993. 384 p.

18. Podlubny I. Fractional differential equations. San Diego: Acad. Press, 1999. 368 p.

19. Khader M.M., Danaf T.E., Hendy A.S. A computational matrix method for solving systems of high order fractional differential equations// Appl. Math. Modell., 2013. Vol. 37, no. 6. P. 4035-4050.

20. Pimenov V., Hendy V. Numerical studies for fractional functional differential equations with delay based on BDF-type shifted Chebyshev approximations // Abstract and Applied Analysis, 2015. Article ID 510875. P 1-12.

21. Alikhanov A.A. Numerical methods of solutions boundary value problems for multi-term veriable-distributed order diffusion equations // Appl.Math.Comput., 2015. Vol. 268. P. 12-22.

22. Meerschaert M.M., Tadjeran C. Finite difference approximations for two sided space fractional partial differential equations // Applied numerical mathematics, 2006. Vol. 65. P. 80-90.

23. Tadjeran C., Meerschaert M.M., Scheffler H.P. A second-order accurate numerical approximation for the fractional diffusion equation // Journal of Computational Physics, 2006. Vol. 213. P. 205-214.

24. Pimenov V.G., Hendy A.S. Numerical methods for the equation with fractional derivative on state and with functional delay on time // Bulletin of the Tambov university. Series: Natural and technical science. 2015. Vol. 20, no 5. P. 1358-1361.

25. Wang H., Wang K., Sircar T. A direct O(Nlog2N) finite difference method for fractional diffusion equations // Journal of Computational Physics, 2010. Vol. 229. P. 8095-8104.

26. Isaacson E., Keller H.B. Analysis of Numerical Methods. New York: Wiley, 1966. 541 p.