Numerical Approximation of Fractional Powers of Elliptic Operators

In this paper, we develop and study algorithms for approximately solving the linear algebraic systems: $\mathcal{A}_h^\alpha u_h = f_h$, $0< \alpha <1$, for $u_h, f_h \in V_h$ with $V_h$ a finite element approximation space. Such problems arise in finite element or finite difference approximations of the problem $\mathcal{A}^\alpha u=f$ with $\mathcal{A}$, for example, coming from a second order elliptic operator with homogeneous boundary conditions. The algorithms are motivated by the recent method of Vabishchevich, 2015, that relates the algebraic problem to a solution of a time-dependent initial value problem on the interval $[0,1]$. Here we develop and study two time stepping schemes based on diagonal Pad\'e approximation to $(1+x)^{-\alpha}$. The first one uses geometrically graded meshes in order to compensate for the singular behavior of the solution for $t$ close to $0$. The second algorithm uses uniform time stepping but requires smoothness of the data $f_h$ in discrete norms. For both methods, we estimate the error in terms of the number of time steps, with the regularity of $f_h$ playing a major role for the second method. Finally, we present numerical experiments for $\mathcal{A}_h$ coming from the finite element approximations of second order elliptic boundary value problems in one and two spatial dimensions.