Numerical approximation of fractional powers of elliptic operators

Abstract

Abstract In this paper, we develop and study algorithms for approximately solving 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 method of Vabishchevich (2015, Numerically solving an equation for fractional powers of elliptic operators. J. Comput. Phys., 282, 289–302) 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é approximation to $(1+x)^{-\alpha }$. The first one uses geometrically graded meshes in order to compensate for the singular behaviour 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.