 # Numerical solution of time-dependent problems with fractional power elliptic operator

An unsteady problem is considered for a space-fractional equation in a bounded domain. A first-order evolutionary equation involves a fractional power of an elliptic operator of second order. Finite element approximation in space is employed. To construct approximation in time, standard two-level schemes are used. The approximate solution at a new time-level is obtained as a solution of a discrete problem with the fractional power of the elliptic operator. A Pade-type approximation is constructed on the basis of special quadrature formulas for an integral representation of the fractional power elliptic operator using explicit schemes. A similar approach is applied in the numerical implementation of implicit schemes. The results of numerical experiments are presented for a test two-dimensional problem.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Many applied mathematical models involve both sub-diffusion (fractional in time) and super-diffusion (fractional in space) operators (see, e.g., Podlubny (1998); Uchaikin (2013)). Super-diffusion 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:

 dudt+Aαu=f(t),0<α<1,0
 u(0)=u0,

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 self-adjoint 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 time-dependent problems of BVP for PDE, the general theory of stability (well-posedness) for operator-difference 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 three-level difference schemes considered in finite-dimensional Hilbert spaces. We emphasize a constructive nature of the general theory of stability for operator-difference schemes, where stability criteria are formulated in the form of operator inequalities, which are easy to verify. In particular, the backward Euler scheme and Crank-Nicolson scheme are unconditionally stable for a non-negative 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 function-vector 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 space-fractional reaction-diffusion 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 (

Bueno-Orovio 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 Dunford-Cauchy 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 pseudo-parabolic equation. For an auxiliary Cauchy problem, the standard two-level 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 one-dimensional problems for the space-fractional 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 Crank-Nicolson 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 one-dimensional transient problems of convection and space-fractional diffusion equation are considered.

In Vabishchevich (2016c), an unsteady problem is considered for a space-fractional diffusion equation in a bounded domain. To construct approximation in time, regularized two-level 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 pseudo-parabolic equation (Vabishchevich (2015)). Some more general unsteady problems are considered in Vabishchevich (2016a, b).

In the present work, standard two-level 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 time-level. When implementing the explicit scheme, the fractional power of the operator is approximated on the basis of Gauss-Jacobi quadrature formulas for the corresponding integral representation. In this case, we have (see Frommer et al. (2014)) a Pade-type 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

with coefficients , . The operator is defined on the set of functions that satisfy on the boundary the following conditions:

 k(x)∂u∂n+g(x)u=0,x∈∂Ω, (2)

where .

In the Hilbert space , we define the scalar product and norm in the standard way:

 (u,v)=∫Ωu(x)v(x)dx,∥u∥=(u,u)1/2.

For the spectral problem

 Aφk=λkφk,x∈Ω,
 k(x)∂φk∂n+g(x)φk=0,x∈∂Ω,

we have

 λ1≤λ2≤...,

and the eigenfunctions form a basis in . Therefore,

 u=∞∑k=1(u,φk)φk.

Let the operator be defined in the following domain:

 D(A)={u | u(x)∈L2(Ω), ∞∑k=0|(u,φk)|2λk<∞}.

The operator is self-adjoint and positive definite:

 A=A∗≥δI,δ>0, (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 :

 Aαu=∞∑k=0(u,φk)λαkφk,0<α<1.

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 first-order equation with the fractional power of the operator . The solution satisfies the equation

 dudt+Aαu=f(t),0

and the initial condition

 u(0)=u0. (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 right-hand 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

 a(u,u)≥δ∥u∥2.

Define the subspace of finite elements and the discrete elliptic operator as

 (Ay,v)=a(y,v),∀ y,v∈Vh.

The fractional power of the operator is defined similarly to . For the spectral problem

 A˜φk=˜λk

we have

 ˜λ1≤˜λ2≤...≤˜λMh,∥˜φk∥=1,k=1,2,...,Mh.

The domain of definition for the operator is

 D(A)={y | y∈Vh, Mh∑k=0|(y,˜φk)|2˜λk<∞}.

The operator acts on the finite dimensional space defined on the domain and, similarly to (3),

 A=A∗≥δ–hI,δ–h>0, (6)

where . For the fractional power of the operator , we suppose

 Aαy=Mh∑k=1(y,˜φk)˜λαk˜φk.

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 :

 dwdt+Aαw=ψ(t),0
 w(0)=w0, (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 right-hand side and initial conditions are sufficiently smooth.

Let us multiply equation (2) by and integrate it over the domain :

 (dwdt,w)+(Aαw,w)=(ψ,w).

In view of the self-adjointness and positive definiteness of the operator , we have

 (dwdt,w)≤(ψ,w).

The right-hand side can be evaluated by the inequality

 (ψ,w)≤∥ψ∥∥w∥.

By virtue of this, we have

 ddt∥w∥≤∥ψ∥.

The latter inequality leads us to the desired a priori estimate:

 ∥w(t)∥≤∥w0∥+∫t0∥ψ(θ)∥dθ. (9)

We will focus on the estimate (9) for the stability of the solution with respect to the initial data and the right-hand 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 two-level schemes. Let be a step of a uniform grid in time such that , . It seems reasonable to begin with the simplest explicit scheme

 wn+1−wnτ+Aαwn=ψn,n=0,1,...,N−1, (10)
 w0=w0. (11)

Advantages and disadvantages of explicit schemes for the standard parabolic problem () are well-known, 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 time-level is determined via (10) as

 wn+1=wn−τAαwn+τψn. (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:

 wn+1=−τArn+wn+τψn,rn=Aα−1wn.

We construct a numerical algorithm that employ the rational approximation of the operator

 A−β,β=1−α,0<β<1.

In this case, we solve standard problems that are related to the operator .

We use an approximation for based on integral representation of a self-adjoint and positive definite operator (see, e.g., Krasnoselskii et al. (1976); Carracedo et al. (2001)):

 A−β=sin(πβ)π∫∞0θ−β(A+θI)−1dθ,0<β<1. (13)

The approximation of is based on the use of one or another quadrature formulas for the right-hand side of (13). Various possibilities in this direction are discussed in Bonito and Pasciak (2015). One possibility is special Gauss-Jacobi 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 right-hand side (13), it is natural to focus on the use of Gauss quadrature formulas. Some possibilities of constructing quadratures directly for half-infinite 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

 θ=μ1−η1+η,μ>0.

From (13), we have

 A−β=2μ1−βsin(πβ)π∫1−1(1−η)−β(1+η)β−1(μ(1−η)I+(1+η)A)−1dη. (14)

To approximate the right-hand side of (14), we apply the Gauss-Jacobi quadrature formula with the weight (see Ralston and Rabinowitz (2001)):

 ∫1−1f(t)(1−η)~α(1+η)~βdη≈M∑m=1ωmf(ηm),α,β>−1.

Here are the roots of the Jacobi polynomial of degree . The weights are given by the formula

 ωm=−2M+~α+~β+2M+~α+~β+1Γ(M+~α+1)Γ(M+~β+1)Γ(M+~α+~β+1)(M+1)!2~α+~βJ′M(ηm;~α,~β)JM+1(ηm;~α,~β)>0,

where denotes the gamma function.

For the fractional power of the operator , from (14), we get

 A−β≈RM(A),RM(A)=M∑m=1dm(cmI+A)−1, (15)

where

 ~α=−β,~β=β−1,dm=2μ1−βsin(πβ)πωm1+ηm,cm=μ1−ηm1+ηm.

In view of (15), the approximate solution of the problem is associated with solving standard problems .

Instead of (12), we employ the scheme

 wn+1=−τARM(A)wn+wn+τψn,n=0,1,...,N−1. (16)

Let us consider the stability conditions for the scheme (11), (16). For a finite-dimensional self-adjoint operator , in addition to the lower bound (6), the following upper bound holds:

 A≤¯¯¯δhI, (17)

where . Thus

 δ–αhI≤Aα≤¯¯¯δαhI,0<α<1.

Similar estimates we have also for :

 γ––hI≤ARM(A)≤¯¯¯γhI,0<α<1, (18)

with some .

###### Theorem 1

If

 τ≤τ0=2¯¯¯γh, (19)

then the scheme (11), (16) is stable in and the solution satisfies the following estimate:

 ∥wn+1∥≤∥w0∥+τn∑j=0∥ψj∥,n=0,1,...,N−1. (20)

[Proof.] From (16), we directly obtain

 ∥wn+1∥≤∥I−τARM(A)∥∥wn∥+τ∥ψn∥,n=0,1,...,N−1.

By (6), (17), (18), we get

 ∥I−τARM(A)∥≤maxδ–h≤z≤¯δh|1−τzRM(z)|≤1

under the restrictions (19). In view of this, we have the level-wise estimate

 ∥wn+1∥≤∥wn∥+τ∥ψn∥,n=0,1,...,N−1,

that results in the estimate (20) for the stability of the solution on the right-hand 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 Pade-type 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 self-adjoint 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

 ¯¯¯γh

where

 ¯¯¯γ(M,α)=M∑m=1dm.

From theorem 1, it follows that for the explicit scheme (11), (16) the time-step 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 two-level schemes with weights (Samarskii (2001); Samarskii et al. (2002)). For a constant weight parameter , we define

 wn+σ=σwn+1+(1−σ)wn.

Instead of (10), let us consider the implicit scheme

 wn+1−wnτ+Aαwn+σ=ψn+σ,n=0,1,...,N−1. (21)

For , the difference scheme (21) is the symmetric scheme (the Crank-Nicolson 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

 wn+σ−wnστ+Aαwn+σ=ψn+σ,n=0,1,...,N−1.

In view of this, the transition to a new time-level involves the solution of the problem

 (νI+Aα)wn+σ=χn,ν=1στ.

For this problem, we construct the rational approximation of the operator

 (νI+Aα)−1,0<α<1.

The necessary approximation is based on the following integral representation:

 (νI+Aα)−1=sin(πα)π∫∞0θαθ2α+2θανcos(πα)+ν2(A+θI)−1dθ, (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

 (νI+Aα)−1=2μ1−αsin(πα)π∫1−1(1−η)−α(1+η)α−1g(η;ν,α)(μ(1−η)I+(1+η)A)−1dη, (23)

where

Further, the Gauss quadrature formula is used (see Gautschi (2004)) with the weight function

 (1−η)−α(1+η)α−1g(ηm;ν,α).

From (23), we get

 (νI+Aα)−1≈RM(A;ν),RM(A;ν)=M∑m=1dνm(cmI+A)−1. (24)

Thereby .

For , an approximate solution is obtained from

 R−1M(A;ν)wn+σ=νwn+ψn+σ,wn+1=1σ(wn+σ−(1−σ)wn),n=0,1,...,N−1. (25)
###### Theorem 2

The difference scheme (11), (25) for and

 R−1M(A;ν)≥νI (26)

is unconditionally stable in and its solution satisfies the a priori estimate

 ∥wn+1∥≤∥w0∥+τn∑j=0∥ψj+σ∥,n=0,1,...,N−1. (27)

[Proof.] The use of (25) means that instead of (21) we employ the scheme

 wn+1−wnτ+Dwn+σ=ψn+σ,n=0,1,...,N−1, (28)

where, in view of (26), we have

 D=D∗≥0,R−1M(A;ν)=νI+D.

Multiplying (29) scalarly in by , we get

 (wn+1−wn,wn+σ)≤τ(ψn+σ,wn+σ). (29)

Taking into account

 wn+σ=(2σ−1)wn+1+(1−σ)(wn+1+wn),

in view of , for the left-hand side of (29), we obtain

 (wn+1−wn,wn+σ)=(2σ−1)∥wn+1∥2−(2σ−1)(wn,wn+1)+(1−σ)(∥wn+1∥2−∥wn∥2)≥(∥wn+1∥−∥wn∥)((2σ−1)∥wn+1∥+(1−σ)(∥wn+1∥+∥wn∥))≥(∥wn+1∥−∥wn∥)∥wn+σ∥.

For the right-hand side of (29), we have

 (ψn+σ,wn+σ)≤∥ψn+σ∥∥wn+σ∥.

In view of this, from (29), we get the inequality

 ∥wn+1∥≤∥wn∥+τ∥ψn+σ∥,

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

 Au=−Δu,x∈Ω,

with the boundary conditions

 ∂u∂n=0,x∈Γ1,x∈Γ2,
 ∂u∂n+gu=0,x∈Γ3,g=const.

We study the case, where the solution depends only on , and . By virtue of this

 Au=−1rddr(rdudr),0
 dudr+gu=0,r=1.

The solution of the spectral problem

 −1rddr(rdφkdr)=λkφk,0
 dφkdr+gφk=0,r=1,

is well-known (see, e.g., Polyanin (2002); Carslaw and Jaeger (1986)). Eigenfunctions are represented as zero-order Bessel functions:

 φk(r)=J0(√λkr),

whereas eigenvalues , where are roots of the equation

 νJ′0(ν)+μJ0(ν)=0. (30)

The general solution of the homogeneous () Cauchy problem for equation (4) is

 u(r,t)=∞∑k=1akexp(−ν2αkt)J0(νkr).

To study the accuracy of the approximate solution of the time-dependent equation with the fractional power of an elliptic operator, we use the exact solution

 u(r,t)=exp(−ν2α1t)J0(ν1r)+1.5exp(−ν2α3t)J3(ν3r),r=(x21+x22)1/2. (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. Figure 2: The solution for different g (α=0.5): left: g=1; center: g=10; right: g=100. Figure 3: The solution for different α (g=10): left: α=0.25; center: α=0.5; right: α=0.75.

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 :

 ε2=∥wN(x)−u(x,T)∥,ε∞=maxx∈Ω|wN(x)−u(x,T)|. Figure 4: Grid: left: 1 — 123 vertices and 208 cells; center: 2 — 461 vertices and 848 cells; right: 3 — 1731 vertices and 3317 cells. Figure 5: Approximation error for β=0.5 (μ=δ, z−β0=0.458821546223).

The finite element approximation in space is based on the use of continuous

Lagrange element, namely, piecewise-linear 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 Krylov-Schur 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 right-hand 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.

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.

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 .

## Conclusion

1. 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 piecewise-linear approximations. An a priori estimate for stability with respect to the initial data and the right-hand side is provided.

2. The explicit scheme is implemented using a Pade-type 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.

3. 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 time-level.

4. Possibilities of the proposed algorithms were demonstrated through numerical solving a test two-dimensional 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 reaction-diffusion 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.
• 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.
• Bueno-Orovio et al. (2014) A. Bueno-Orovio, D. Kay, and K. Burrage. Fourier spectral methods for fractional-in-space reaction-diffusion 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 fractional-in-space reaction-diffusion 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 half-infinite 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 gauss-type 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. Data-sparse approximation to the operator-valued functions of elliptic operator. Mathematics of Computation, 73(247):1297–1324, 2004.
• Gavrilyuk et al. (2005) I. Gavrilyuk, W. Hackbusch, and B. Khoromskij. Data-sparse approximation to a class of operator-valued 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 fractional-in-space 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 space-fractional 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 978-3-642-23098-1.
• 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. Springer-Verlag, 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 second-order 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 space-fractional diffusion equation. International Journal of Numerical Analysis and Modeling, 13(2):296–309, 2016a.
• Vabishchevich (2014) P. N. Vabishchevich. Additive Operator-Difference 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 space-fractional 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 space-fractional 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.