1 Introduction
Many applied mathematical models involve both subdiffusion (fractional in time) and superdiffusion (fractional in space) operators (see, e.g., Podlubny (1998); Uchaikin (2013)). Superdiffusion problems are treated as evolutionary problems with a fractional power of an elliptic operator. For example, suppose that in a bounded domain on the set of functions , there is defined the operator : . We seek the solution of the Cauchy problem for the equation with the fractional power elliptic operator:
for a given , using the notation .
For approximation in space, we can apply finite volume or finite element methods oriented to using arbitrary domains and irregular computational grids (Knabner and Angermann (2003); Quarteroni and Valli (1994)). After this, we formulate the corresponding Cauchy problem with a fractional power of a selfadjoint positive definite discrete elliptic operator (see Bonito and Pasciak (2015); Szekeres and Izsák (2016)) — a fractional power of a symmetric positive definite matrix (Higham (2008)).
In the study of difference schemes for timedependent problems of BVP for PDE, the general theory of stability (wellposedness) for operatordifference schemes (Samarskii (2001); Samarskii et al. (2002)) is in common use. At the present time, the exact (matching necessary and sufficient) conditions for stability are obtained for a wide class of two and threelevel difference schemes considered in finitedimensional Hilbert spaces. We emphasize a constructive nature of the general theory of stability for operatordifference schemes, where stability criteria are formulated in the form of operator inequalities, which are easy to verify. In particular, the backward Euler scheme and CrankNicolson scheme are unconditionally stable for a nonnegative operator.
Problems in numerical solving unsteady problems with fractional powers of operators appear in using the simplest explicit approximations in time. A practical implementation of such approach requires the matrix functionvector multiplication. For such problems, different approaches (see
Higham (2008)) are available. Algorithms for solving systems of linear equations associated with fractional elliptic equations that are based on Krylov subspace methods with the Lanczos approximation are discussed, e.g., in Ilić et al. (2008). A comparative analysis of the contour integral method, the extended Krylov subspace method, and the preassigned poles and interpolation nodes method for solving spacefractional reactiondiffusion equations is presented in
Burrage et al. (2012). The simplest variant is associated with the explicit construction of the solution using the eigenvalues and eigenfunctions of the elliptic operator with diagonalization of the corresponding matrix (
BuenoOrovio et al. (2014); Ilic et al. (2006)). Unfortunately, all these approaches demonstrate very high computational complexity for multidimensional problems.There does exist a general approach to solve approximately equations involving a fractional power of operators based on an approximation of the original operator and then taking fractional power of its discrete variant. Using the DunfordCauchy formula the elliptic operator is represented as a contour integral in the complex plane. Further, applying appropriate quadratures with integration nodes in the complex plane, it is necessary to select a proper method that involves only inversion of the original operator. The approximate operator is treated as a sum of resolvents (Gavrilyuk et al. (2004, 2005)) ensuring the exponential convergence of quadrature approximations. In Bonito and Pasciak (2015), there was presented a more promising variant of using quadrature formulas with nodes on the real axis, which are constructed on the basis of the corresponding integral representation for the power operator (see Krasnoselskii et al. (1976); Carracedo et al. (2001)). In this case, the inverse operator of the problem has an additive representation, where each term is an inverse of the original elliptic operator. A similar rational approximation to the fractional Laplacian operator is studied in Aceto and Novati (2017).
We have proposed (Vabishchevich (2015)) a numerical algorithm to solve an equation for fractional power elliptic operators that is based on a transition to a pseudoparabolic equation. For an auxiliary Cauchy problem, the standard twolevel schemes are applied. The computational algorithm is simple for practical use, robust, and applicable to solving a wide class of problems. A small number of time steps is required to find a solution. This computational algorithm for solving equations with fractional power operators is promising for transient problems.
As for onedimensional problems for the spacefractional diffusion equation, an analysis of stability and convergence for this equation was conducted in Jin et al. (2014) using finite element approximation in space. A similar study for the CrankNicolson scheme was conducted earlier in Tadjeran et al. (2006) using finite difference approximations in space. We highlight separately the works Huang et al. (2008); Sousa (2012); Meerschaert and Tadjeran (2004), where numerical methods for solving onedimensional transient problems of convection and spacefractional diffusion equation are considered.
In Vabishchevich (2016c), an unsteady problem is considered for a spacefractional diffusion equation in a bounded domain. To construct approximation in time, regularized twolevel schemes are used (see Vabishchevich (2014)). The numerical implementation is based on solving the equation with the fractional power of the elliptic operator using an auxiliary Cauchy problem for a pseudoparabolic equation (Vabishchevich (2015)). Some more general unsteady problems are considered in Vabishchevich (2016a, b).
In the present work, standard twolevel schemes are applied to solve numerically a Cauchy problem for an evolutionary equation of first order with a fractional power of the operator. The numerical implementation is based on the rational approximation of the operator at a new timelevel. When implementing the explicit scheme, the fractional power of the operator is approximated on the basis of GaussJacobi quadrature formulas for the corresponding integral representation. In this case, we have (see Frommer et al. (2014)) a Padetype approximation of the power function with a fractional exponent. A similar approach is used when considering implicit schemes.
The paper is organized as follows. The formulation of an unsteady problem containing a fractional power of an elliptic operator is given in Section 2. Here finite element approximations in space are also discussed. In Section 3, we construct the explicit approximation in time and investigate its stability. The numerical implementation is based on the rational approximation of the fractional power operator. Implicit schemes are considered in Section 4. The results of numerical experiments are described in Section 5. At the end of the work the main results of our study are summarized.
2 Problem formulation
In a bounded polygonal domain , with the Lipschitz continuous boundary , we search the solution for the problem with a fractional power of an elliptic operator. Define the elliptic operator as
(1) 
with coefficients , . The operator is defined on the set of functions that satisfy on the boundary the following conditions:
(2) 
where .
In the Hilbert space , we define the scalar product and norm in the standard way:
For the spectral problem
we have
and the eigenfunctions form a basis in . Therefore,
Let the operator be defined in the following domain:
The operator is selfadjoint and positive definite:
(3) 
where is the identity operator in . For , we have . In applications, the value of is unknown (the spectral problem must be solved). Therefore, we assume that in (3). Let us assume for the fractional power of the operator :
More general and mathematically complete definition of fractional powers of elliptic operators is given in Carracedo et al. (2001); Yagi (2009).
We seek the solution of a Cauchy problem for the evolutionary firstorder equation with the fractional power of the operator . The solution satisfies the equation
(4) 
and the initial condition
(5) 
The key issue in the study of computational algorithms for solving the Cauchy problem (4), (5) is to establish the stability of the approximate solution with respect to small perturbations of the initial data and the righthand side in various norms.
To solve numerically the problem (4), (5), we employ finite element approximations in space (see, e.g., Brenner and Scott (2008); Thomée (2006)). For (1) and (2), we define the bilinear form
By (3), we have
Define the subspace of finite elements and the discrete elliptic operator as
The fractional power of the operator is defined similarly to . For the spectral problem
we have
The domain of definition for the operator is
The operator acts on the finite dimensional space defined on the domain and, similarly to (3),
(6) 
where . For the fractional power of the operator , we suppose
The use of finite element approximations for fractional power elliptic operators is discussed in detail, for instance, in the works Acosta and Borthagaray (2017); Szekeres and Izsák (2016).
For the problem (4), (5), we put into the correspondence the operator equation for :
(7) 
(8) 
where , with denoting projection onto .
Now we obtain an elementary a priori estimate for the solution of (
2), (3) assuming that the solution of the problem, coefficients of the elliptic operator, the righthand side and initial conditions are sufficiently smooth.Let us multiply equation (2) by and integrate it over the domain :
In view of the selfadjointness and positive definiteness of the operator , we have
The righthand side can be evaluated by the inequality
By virtue of this, we have
The latter inequality leads us to the desired a priori estimate:
(9) 
We will focus on the estimate (9) for the stability of the solution with respect to the initial data and the righthand side in constructing discrete analogs of the problem (7), (8).
3 Explicit scheme
To solve numerically the problem (7), (8), we use simplest explicit and implicit twolevel schemes. Let be a step of a uniform grid in time such that , . It seems reasonable to begin with the simplest explicit scheme
(10) 
(11) 
Advantages and disadvantages of explicit schemes for the standard parabolic problem () are wellknown, i.e., these are a simple computational implementation and a time step restriction (see, e.g., Samarskii (2001); Samarskii et al. (2002)). In our case (), the main drawback (conditional stability) remains, whereas the advantage in terms of implementation simplicity does not exist. The approximate solution at a new timelevel is determined via (10) as
(12) 
Thus, we must calculate . In view of these problems, considering the scheme (10), it is more correct to speak about the scheme with the explicit approximations in time in contrast to the standard fully explicit scheme.
The numerical implementation of (12) is based on the following representation:
We construct a numerical algorithm that employ the rational approximation of the operator
In this case, we solve standard problems that are related to the operator .
We use an approximation for based on integral representation of a selfadjoint and positive definite operator (see, e.g., Krasnoselskii et al. (1976); Carracedo et al. (2001)):
(13) 
The approximation of is based on the use of one or another quadrature formulas for the righthand side of (13). Various possibilities in this direction are discussed in Bonito and Pasciak (2015). One possibility is special GaussJacobi quadrature formulas considered in Frommer et al. (2014); Aceto and Novati (2017). Just this approximation of the operator is used in the present work.
To achieve higher accuracy in approximating the the righthand side (13), it is natural to focus on the use of Gauss quadrature formulas. Some possibilities of constructing quadratures directly for halfinfinite intervals are investigated, for example, in the work Gautschi (1991). The classical Gauss quadrature formulas can be used via introducing a new variable of integration.
Suppose (see Frommer et al. (2014)) that
From (13), we have
(14) 
To approximate the righthand side of (14), we apply the GaussJacobi quadrature formula with the weight (see Ralston and Rabinowitz (2001)):
Here are the roots of the Jacobi polynomial of degree . The weights are given by the formula
where denotes the gamma function.
For the fractional power of the operator , from (14), we get
(15) 
where
In view of (15), the approximate solution of the problem is associated with solving standard problems .
Instead of (12), we employ the scheme
(16) 
Let us consider the stability conditions for the scheme (11), (16). For a finitedimensional selfadjoint operator , in addition to the lower bound (6), the following upper bound holds:
(17) 
where . Thus
Similar estimates we have also for :
(18) 
with some .
Theorem 1
[Proof.] From (16), we directly obtain
under the restrictions (19). In view of this, we have the levelwise estimate
that results in the estimate (20) for the stability of the solution on the righthand side and the initial conditions.
It should be noted that the estimate (20) for the difference scheme (11), (16) is a discrete analog of the a priori estimate (9) for the problem (7), (8).
Special attention (see Frommer et al. (2014); Aceto and Novati (2017)) should be given to the problem of choosing the parameter in (14). Taking into account the definition of the operator , we are interested in the best approximation of for the smallest (principal) eigenvalue . In Frommer et al. (2014), there is established a remarkable fact that corresponds to a Padetype approximation for the function with expansion point . Thus, the optimal choice corresponds to the selection . In this case, in (18), we have . The computational complexity of finding (the principal eigenvalue of a discrete selfadjoint elliptic operator of second order) is not high. To this end, it is possible to use standard algorithms (see, e.g., Saad (2011)) and the corresponding software (see Hernandez et al. (2005)).
The function for is a positive and monotonically increasing function. In view of this, taking into account (15), we have
where
From theorem 1, it follows that for the explicit scheme (11), (16) the timestep restrictions do not depend on discretization in space, but depend on the power and the number of approximation nodes .
4 Implicit scheme
Unconditionally stable schemes are constructed on the basis of implicit approximations in time. Here we consider standard twolevel schemes with weights (Samarskii (2001); Samarskii et al. (2002)). For a constant weight parameter , we define
Instead of (10), let us consider the implicit scheme
(21) 
For , the difference scheme (21) is the symmetric scheme (the CrankNicolson scheme). It approximates the differential problem with the second order by , whereas for other values of , we have only the first order.
Rewrite the scheme (21) in the form
In view of this, the transition to a new timelevel involves the solution of the problem
For this problem, we construct the rational approximation of the operator
The necessary approximation is based on the following integral representation:
(22) 
taken from the work Carracedo et al. (2001) (see Proposition 5.3.2).
Using the new variable , from (22), we arrive at the representation
(23) 
where
Further, the Gauss quadrature formula is used (see Gautschi (2004)) with the weight function
From (23), we get
(24) 
Thereby .
For , an approximate solution is obtained from
(25) 
Theorem 2
[Proof.] The use of (25) means that instead of (21) we employ the scheme
(28) 
where, in view of (26), we have
Multiplying (29) scalarly in by , we get
(29) 
Taking into account
in view of , for the lefthand side of (29), we obtain
For the righthand side of (29), we have
In view of this, from (29), we get the inequality
which provides the estimate (27).
5 Numerical experiments
The test problem under the consideration is constructed using the exact solution of the problem in the unit circle (see Vabishchevich (2016b)). The computational domain is a quarter of the circle (see Fig. 1). Consider the equation
with the boundary conditions
We study the case, where the solution depends only on , and . By virtue of this
The solution of the spectral problem
is wellknown (see, e.g., Polyanin (2002); Carslaw and Jaeger (1986)). Eigenfunctions are represented as zeroorder Bessel functions:
whereas eigenvalues , where are roots of the equation
(30) 
The general solution of the homogeneous () Cauchy problem for equation (4) is
To study the accuracy of the approximate solution of the timedependent equation with the fractional power of an elliptic operator, we use the exact solution
(31) 
The values of the roots for different values of the boundary condition are given in Table 1. The exact solution for at different values of and is shown in Figs. 2 and 3.
1  1.25578371  2.17949660  2.38090166 

3  7.15579917  7.95688342  8.56783165 
Predictions were performed on a sequence of refining grids, which are shown in Fig. 4
. The numerical solution is compared with the exact one at the final time moment
. Error estimation is performed in and :grid  

1  1.57959231369  4225.51507674  
1  1.57699272630  2  1.57763558651  17104.1780271 
3  1.57715815735  74989.7519112  
1  4.76409956820  4252.23867499  
10  4.75020542941  2  4.75363764524  17143.1279728 
3  4.75108440807  74989.7519112  
1  5.68846224707  7310.80621520  
100  5.66869271459  2  5.67358161306  22017.7463507 
3  5.66994292109  74989.7519112 
The finite element approximation in space is based on the use of continuous
Lagrange element, namely, piecewiselinear elements. The calculations were performed using the computing platform FEniCS for solving partial differential equations (website
http://fenicsproject.org, Logg et al. (2012); Alnæs et al. (2015)). To solve spectral problems with symmetrical matrices, we use the SLEPc library (Scalable Library for Eigenvalue Problem Computations, http://slepc.upv.es, Hernandez et al. (2005)). We apply the KrylovSchur algorithm, a variation of the Arnoldi method, proposed by Stewart (2001).Table 2 presents the lower and upper bounds of the operator spectrum (see (6), (17)) on various grids for different values of the parameter in the boundary condition. A comparison with the exact minimum eigenvalue demonstrates good accuracy in evaluation of , which increases with the grid refinement. It is easy to see a significant dependence of the maximum eigenvalue on the grid size.
Peculiarities of the approximation (15) are illustrated by the accuracy of the approximation of the function for (see (3)). Figure 5 shows the absolute error arising from the approximation of by the function for and . In this case, (see Table 1) and . We see higher accuracy near the left boundary , whereas for large , the approximation accuracy decreases. The effect of increasing accuracy with increasing number of nodes of the quadrature formula is clearly observed.
Decreasing the approximation accuracy at , we can increase the accuracy for other values of . For example, Fig. 6 demonstrates the approximation accuracy for . In this case . We need to have good accuracy for small , and therefore, in calculations, we are guided by the choice of . The dependence of the approximation accuracy of the function on the value of is shown in Figs. 7 and 8. In these figures, we observe a significant drop in accuracy with decreasing .
The numerical implementation of the explicit scheme (11), (16) involves the approximation of the operator by the expression . Peculiarities of this approximation at , are shown in Fig. 9. It should be noted that the operator is bounded and the constant on the righthand side of (18) for the corresponding values of is presented via dotted lines.
The upper bounds of the operator are given in Table 3 for . Increasing with increasing the number of quadrature formula nodes results from increasing the accuracy of approximation of the unbounded operator . As decreases, the value of decreases drastically.
5  4.4602175  21.794966  142.00220 

10  6.3106349  43.589932  401.45610 
20  8.9256294  87.179864  1135.3565 
40  12.623116  174.35973  3211.1792 
Now we present numerical results obtained using the explicit scheme (11), (16). We confine ourselves to the case with the value of the boundary condition parameter . It is interesting to identify the dependence of accuracy on grids in space and time. In our case, we should also study the influence of the number of quadrature formula nodes . Table 4 demonstrates the numerical solution convergence for decreasing the time step and increasing the accuracy of approximation of the fractional power operator. Increasing the accuracy of an approximate solution due to spatial grid refinement is shown in Table 5.
25  50  100  200  

5  0.00436648  0.00207328  0.00094905  0.00041937  
0.01896717  0.00927482  0.00447550  0.00216564  
10  0.00507635  0.00277981  0.00164515  0.00108352  
0.02186657  0.01217982  0.00738292  0.00499609  
20  0.00507902  0.00278251  0.00164787  0.00108627  
0.02187724  0.01219044  0.00739354  0.00500671  
40  0.00507902  0.00278251  0.00164787  0.00108627  
0.02187723  0.01219043  0.00739352  0.00500669 
grid  25  50  100  200  

1  0.00641387  0.00419465  0.00310833  0.00257568  
0.02505950  0.01634587  0.01202969  0.00988175  
2  0.00507902  0.00278251  0.00164787  0.00108627  
0.02187724  0.01219044  0.00739354  0.00500671  
3  0.00472777  0.00241442  0.00126921  0.00069981  
0.02077071  0.01081355  0.00588308  0.00342987 
The numerical implementation of implicit schemes is associated (see (21)) with the function . It should be noted that and therefore, this parameter is large enough in numerical solving unsteady problems. The function , which corresponds to our test problem for , is shown in Fig. 10. As we noted earlier, if , then the optimal value is . This value is also used in our calculations for .
Figure 11 shows the approximating function for with . The approximation accuracy for various values of is presented in Figs. 12 and 13. Operator approximations were designed using the package ORTHPOL (see Gautschi (1994)) developed for constructing Gauss quadrature formulas with an arbitrary weight function.
The accuracy of the approximate solution of the test problem obtained using the implicit scheme (11), (25) was investigated for and . For the fully implicit scheme (), Table 6 demonstrates the dependence of the solution accuracy on the grid in time for various numbers of the quadrature formula nodes .
25  50  100  200  

5  0.00326281  0.00111138  0.00044581  0.00078676  
0.01672817  0.00699940  0.00208775  0.00313725  
10  0.00394436  0.00173870  0.00063823  0.00018398  
0.01975034  0.01002383  0.00511271  0.00264674  
10  0.00394696  0.00174126  0.00064056  0.00018400  
0.01976054  0.01003442  0.00512346  0.00265756  
10  0.00394696  0.00174126  0.00064056  0.00018399  
0.01976037  0.01003433  0.00512339  0.00265750 
Conclusion

There is considered a nonclassical problem with the initial data, which is described by an evolutionary equation of first order with a fractional power of an elliptic operator. The multidimensional problem is approximated in space using standard finite element piecewiselinear approximations. An a priori estimate for stability with respect to the initial data and the righthand side is provided.

The explicit scheme is implemented using a Padetype approximation for the fractional power elliptic operator. Sufficient conditions for the stability of the explicit scheme are formulated. They do not depend on spatial grid steps.

Rational approximation is employed to implement implicit schemes. It is based on a computational generation of Gauss quadrature formulas for an integral representation of the operator of transition to a new timelevel.

Possibilities of the proposed algorithms were demonstrated through numerical solving a test twodimensional problem.
Acknowledgements
The publication was financially supported by the Ministry of Education and Science of the Russian Federation (the Agreement # 02.a03.21.0008).
References
 Aceto and Novati (2017) L. Aceto and P. Novati. Rational approximation to the fractional Laplacian operator in reactiondiffusion problems. SIAM Journal on Scientific Computing, 39(1):A214–A228, 2017.
 Acosta and Borthagaray (2017) G. Acosta and J. P. Borthagaray. A fractional laplace equation: Regularity of solutions and finite element approximations. SIAM Journal on Numerical Analysis, 55(2):472–495, 2017.
 Alnæs et al. (2015) M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells. The fenics project version 1.5. Archive of Numerical Software, 3(100), 2015. doi: 10.11588/ans.2015.100.20553.
 Bonito and Pasciak (2015) A. Bonito and J. Pasciak. Numerical approximation of fractional powers of elliptic operators. Mathematics of Computation, 84(295):2083–2110, 2015.
 Brenner and Scott (2008) S. C. Brenner and L. R. Scott. The mathematical theory of finite element methods. Springer, New York, 2008.
 BuenoOrovio et al. (2014) A. BuenoOrovio, D. Kay, and K. Burrage. Fourier spectral methods for fractionalinspace reactiondiffusion equations. BIT Numerical Mathematics, 54(4):937–954, 2014.
 Burrage et al. (2012) K. Burrage, N. Hale, and D. Kay. An efficient implicit fem scheme for fractionalinspace reactiondiffusion equations. SIAM Journal on Scientific Computing, 34(4):A2145–A2172, 2012.
 Carracedo et al. (2001) C. M. Carracedo, M. S. Alix, and M. Sanz. The Theory of Fractional Powers of Operators. Elsevier, Amsterdam, 2001.
 Carslaw and Jaeger (1986) H. S. Carslaw and J. C. Jaeger. Conduction of Heat in Solids. Clarendon Press, 1986.
 Frommer et al. (2014) A. Frommer, S. Güttel, and M. Schweitzer. Efficient and stable arnoldi restarts for matrix functions based on quadrature. SIAM Journal on Matrix Analysis and Applications, 35(2):661–683, 2014.
 Gautschi (1991) W. Gautschi. Quadrature formulae on halfinfinite intervals. BIT Numerical Mathematics, 31(3):437–446, 1991.
 Gautschi (1994) W. Gautschi. Algorithm 726: ORTHPOL — a package of routines for generating orthogonal polynomials and gausstype quadrature rules. ACM Transactions on Mathematical Software, 20(1):21–62, 1994.
 Gautschi (2004) W. Gautschi. Orthogonal Polynomials: Computation and Approximation. Oxford University Press, 2004.
 Gavrilyuk et al. (2004) I. Gavrilyuk, W. Hackbusch, and B. Khoromskij. Datasparse approximation to the operatorvalued functions of elliptic operator. Mathematics of Computation, 73(247):1297–1324, 2004.
 Gavrilyuk et al. (2005) I. Gavrilyuk, W. Hackbusch, and B. Khoromskij. Datasparse approximation to a class of operatorvalued functions. Mathematics of Computation, 74(250):681–708, 2005.
 Hernandez et al. (2005) V. Hernandez, J. E. Roman, and V. Vidal. SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems. ACM Transactions on Mathematical Software (TOMS), 31(3):351–362, 2005.
 Higham (2008) N. J. Higham. Functions of Matrices: Theory and Computation. SIAM, Philadelphia, 2008.
 Huang et al. (2008) Q. Huang, G. Huang, and H. Zhan. A finite element solution for the fractional advection–dispersion equation. Advances in Water Resources, 31(12):1578–1589, 2008.
 Ilic et al. (2006) M. Ilic, F. Liu, I. Turner, and V. Anh. Numerical approximation of a fractionalinspace diffusion equation. II with nonhomogeneous boundary conditions. Fractional Calculus and Applied Analysis, 9(4):333–349, 2006.
 Ilić et al. (2008) M. Ilić, I. W. Turner, and V. Anh. A numerical solution using an adaptively preconditioned lanczos method for a class of linear systems related with the fractional poisson equation. International Journal of Stochastic Analysis, Article ID 104525:26 pages, 2008.
 Jin et al. (2014) B. Jin, R. Lazarov, J. Pasciak, and Z. Zhou. Error analysis of finite element methods for spacefractional parabolic equations. SIAM J. Numer. Anal., 52(5):2272–2294, 2014.
 Knabner and Angermann (2003) P. Knabner and L. Angermann. Numerical Methods for Elliptic and Parabolic Partial Differential Equations. Springer, New York, 2003.
 Krasnoselskii et al. (1976) M. A. Krasnoselskii, P. P. Zabreiko, E. I. Pustylnik, and S. P. E. Integral Operators in Spaces of Summable Functions. Noordhoff International Publishing, 1976.
 Logg et al. (2012) A. Logg, K.A. Mardal, G. N. Wells, et al. Automated Solution of Differential Equations by the Finite Element Method. Springer, 2012. ISBN 9783642230981. doi: 10.1007/9783642230998.
 Meerschaert and Tadjeran (2004) M. M. Meerschaert and C. Tadjeran. Finite difference approximations for fractional advection–dispersion flow equations. Journal of Computational and Applied Mathematics, 172(1):65–77, 2004.
 Podlubny (1998) I. Podlubny. Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications. Academic Press, 1998.
 Polyanin (2002) A. D. Polyanin. Handbook of Linear Partial Differential Equations for Engineers and Scientists. Chapman and Hall/CRC, 2002. ISBN 9781584882992,1584882999.
 Quarteroni and Valli (1994) A. Quarteroni and A. Valli. Numerical Approximation of Partial Differential Equations. SpringerVerlag, Berlin, 1994.
 Ralston and Rabinowitz (2001) A. Ralston and P. Rabinowitz. A First Course in Numerical Analysis. Dover Publications, 2001.
 Saad (2011) Y. Saad. Numerical Methods for Large Eigenvalue Problems. SIAM, 2011.
 Samarskii (2001) A. A. Samarskii. The Theory of Difference Schemes. Marcel Dekker, New York, 2001.
 Samarskii et al. (2002) A. A. Samarskii, P. P. Matus, and P. N. Vabishchevich. Difference Schemes with Operator Factors. Kluwer Academic, Dordrecht, 2002.
 Sousa (2012) E. Sousa. A second order explicit finite difference method for the fractional advection diffusion equation. Computers & Mathematics with Applications, 64(10):3141–3152, 2012.
 Stewart (2001) G. W. Stewart. A Krylov–Schur algorithm for large eigenproblems. SIAM Journal on Matrix Analysis and Applications, 23(3):601–614, 2001.
 Szekeres and Izsák (2016) B. J. Szekeres and F. Izsák. Finite element approximation of fractional order elliptic boundary value problems. Journal of Computational and Applied Mathematics, 292:553–561, 2016.
 Tadjeran et al. (2006) C. Tadjeran, M. M. Meerschaert, and H.P. Scheffler. A secondorder accurate numerical approximation for the fractional diffusion equation. Journal of Computational Physics, 213(1):205–213, 2006.
 Thomée (2006) V. Thomée. Galerkin finite element methods for parabolic problems. Springer Verlag, Berlin, 2006.
 Uchaikin (2013) V. V. Uchaikin. Fractional derivatives for physicists and engineers. Higher Education Press, 2013.
 Vabishchevich (2016a) P. Vabishchevich. Numerical solution of nonstationary problems for a convection and a spacefractional diffusion equation. International Journal of Numerical Analysis and Modeling, 13(2):296–309, 2016a.
 Vabishchevich (2014) P. N. Vabishchevich. Additive OperatorDifference Schemes: Splitting Schemes. de Gruyter, Berlin, 2014.
 Vabishchevich (2015) P. N. Vabishchevich. Numerically solving an equation for fractional powers of elliptic operators. Journal of Computational Physics, 282(1):289–302, 2015.
 Vabishchevich (2016b) P. N. Vabishchevich. Numerical solving unsteady spacefractional problems with the square root of an elliptic operator. Mathematical Modelling and Analysis, 21(2):220–238, 2016b.
 Vabishchevich (2016c) P. N. Vabishchevich. Numerical solution of nonstationary problems for a spacefractional diffusion equation. Fractional Calculus and Applied Analysis, 19(1):116–139, 2016c.
 Yagi (2009) A. Yagi. Abstract Parabolic Evolution Equations and Their Applications. Springer, Berlin, 2009.