We discuss, study, and compare experimentally three methods for solving the system of algebraic equations $\mathbb{A}^\alpha \bf{u}=\bf{f}$, $0< \alpha <1$, where $\mathbb{A}$ is a symmetric and positive definite matrix obtained from finite difference or finite element approximations of second order elliptic problems in $\mathbb{R}^d$, $d=1,2,3$. The first method, introduced by Harizanov et.al, based on the best uniform rational approximation (BURA) $r_\alpha(t)$ of $t^{1-\alpha}$ for $0 \le t \le 1$, is used to get the rational approximation of $t^{-\alpha}$ in the form $t^{-1}r_\alpha(t)$. Here we develop another method, denoted by R-BURA, that is based on the best rational approximation $r_{1-\alpha}(t)$ of $t^\alpha$ on the interval $[0,1]$ and approximates $t^{-\alpha}$ via $r^{-1}_{1-\alpha}(t)$. The third method, introduced and studied by Bonito and Pasciak, is based on an exponentially convergent quadrature scheme for the Dundord-Taylor integral representation of the fractional powers of elliptic operators. All three methods reduce the solution of the system $\mathbb{A}^\alpha \bf{u}=\bf{f}$ to solving a number of equations of the type $(\mathbb{A} +c\mathbb{I})\bf{u}= \bf{f}$, $c \ge 0$. Comprehensive numerical experiments on model problems with $\mathbb A$ obtained by approximation of elliptic equations in one and two spatial dimensions are used to compare the efficiency of these three algorithms depending on the fractional power $\alpha$. The presented results prove the concept of the new R-BURA method, which performs well for $\alpha$ close to $1$ in contrast to BURA, which performs well for $\alpha $ close to $0$. As a result, we show theoretically and experimentally, that they have mutually complementary advantages.

Here we study theoretically and compare experimentally an efficient method for solving systems of algebraic equations, where the matrix comes from the discretization of a fractional diffusion operator. More specifically, we focus on matrices obtained from finite difference or finite element approximation of second order elliptic problems in $\mathbb R^d$, $d=1,2,3$. The proposed methods are based on the best uniform rational approximation (BURA) $r_{\alpha,k}(t)$ of $t^{\alpha}$ on $[0,1]$. Here $r_{\alpha,k}$ is a rational function of $t$ involving numerator and denominator polynomials of degree at most $k$. We show that the proposed method is exponentially convergent with respect to $k$ and has some attractive properties. First, it reduces the solution of the nonlocal system to solution of $k$ systems with matrix $(A +c_j I)$ and $c_j>0$, $j=1,2,\ldots,k$. Thus, good computational complexity can be achieved if fast solvers are available for such systems. Second, the original problem and its rational approximation in the finite difference case are positivity preserving. In the finite element case, positivity preserving results when mass lumping is employed under some mild conditions on the mesh. Further, we prove that in the mass lumping case, the scheme still leads to the expected rate of convergence, at times assuming additional regularity on the right hand side. Finally, we present comprehensive numerical experiments on a number of model problems for various $\alpha$ in one and two spatial dimensions. These illustrate the computational behavior of the proposed method and compare its accuracy and efficiency with that of other methods developed by Harizanov et. al. and Bonito and Pasciak.