  This study is motivated by the recent development in the fractional calculus
and its applications. During last few years, several different techniques are
proposed to localize the nonlocal fractional diffusion operator. They are based
on transformation of the original problem to a local elliptic or
pseudoparabolic problem, or to an integral representation of the solution, thus
increasing the dimension of the computational domain. More recently, an
alternative approach aimed at reducing the computational complexity was
developed. The linear algebraic system $\cal A^\alpha \bf u=\bf f$, $0< \alpha
<1$ is considered, where $\cal A$ is a properly normalized (scalded) symmetric
and positive definite matrix obtained from finite element or finite difference
approximation of second order elliptic problems in $\Omega\subset\mathbb{R}^d$,
$d=1,2,3$. The method is based on best uniform rational approximations (BURA)
of the function $t^{\beta-\alpha}$ for $0 < t \le 1$ and natural $\beta$.
  The maximum principles are among the major qualitative properties of linear
elliptic operators/PDEs. In many studies and applications, it is important that
such properties are preserved by the selected numerical solution method. In
this paper we present and analyze the properties of positive approximations of
$\cal A^{-\alpha}$ obtained by the BURA technique. Sufficient conditions for
positiveness are proven, complemented by sharp error estimates. The theoretical
results are supported by representative numerical tests.
