  We present and study a novel numerical algorithm to approximate the action of
$T^\beta:=L^{-\beta}$ where $L$ is a symmetric and positive definite unbounded
operator on a Hilbert space $H_0$. The numerical method is based on a
representation formula for $T^{-\beta}$ in terms of Bochner integrals involving
$(I+t^2L)^{-1}$ for $t\in(0,\infty)$.
  To develop an approximation to $T^\beta$, we introduce a finite element
approximation $L_h$ to $L$ and base our approximation to $T^\beta$ on
$T_h^\beta:= L_h^{-\beta}$. The direct evaluation of $T_h^{\beta}$ is extremely
expensive as it involves expansion in the basis of eigenfunctions for $L_h$.
The above mentioned representation formula holds for $T_h^{-\beta}$ and we
propose three quadrature approximations denoted generically by $Q_h^\beta$. The
two results of this paper bound the errors in the $H_0$ inner product of
$T^\beta-T_h^\beta\pi_h$ and $T_h^\beta-Q_h^\beta$ where $\pi_h$ is the $H_0$
orthogonal projection into the finite element space. We note that the
evaluation of $Q_h^\beta$ involves application of $(I+(t_i)^2L_h)^{-1}$ with
$t_i$ being either a quadrature point or its inverse. Efficient solution
algorithms for these problems are available and the problems at different
quadrature points can be straightforwardly solved in parallel. Numerical
experiments illustrating the theoretical estimates are provided for both the
quadrature error $T_h^\beta-Q_h^\beta$ and the finite element error
$T^\beta-T_h^\beta\pi_h$.
