Raytcho Lazarov

The survey is devoted to numerical solution of the fractional equation $A^\alpha u=f$, $0 < \alpha <1$, where $A$ is a symmetric positive definite operator corresponding to a second order elliptic boundary value problem in a bounded domain $\Omega$ in $\mathbb R^d$. The operator fractional power is a non-local operator and is defined through the spectrum. Due to growing interest and demand in applications of sub-diffusion models to physics and engineering, in the last decade, several numerical approaches have been proposed, studied, and tested. We consider discretizations of the elliptic operator $A$ by using an $N$-dimensional finite element space $V_h$ or finite differences over a uniform mesh with $N$ grid points. The numerical solution of this equation is based on the following three equivalent representations of the solution: (1) Dunford-Taylor integral formula (or its equivalent Balakrishnan formula), (2) extension of the a second order elliptic problem in $\Omega \times (0,\infty)\subset \mathbb R^{d+1}$ (with a local operator) or as a pseudo-parabolic equation in the cylinder $(x,t) \in \Omega \times (0,1) $, (3) spectral representation and the best uniform rational approximation (BURA) of $z^\alpha$ on $[0,1]$. Though substantially different in origin and their analysis, these methods can be interpreted as some rational approximation of $A^{-\alpha}$. In this paper we present the main ideas of these methods and the corresponding algorithms, discuss their accuracy, computational complexity and compare their efficiency and robustness.

In this paper we consider one particular mathematical problem of this large area of fractional powers of self-adjoined elliptic operators, defined either by Dunford-Taylor-like integrals or by the representation through the spectrum of the elliptic operator. Due to the mathematical modeling of various non-local phenomena using such operators recently a number of numerical methods for solving equations involving operators of fractional order were introduced, studied, and tested. Here we consider the discrete counterpart of such problems obtained from finite difference or finite element approximations of the corresponding elliptic problems. In this report we provide all necessary information regarding the best uniform rational approximation (BURA) $r_{k,\alpha}(t) := P_k(t)/Q_k(t)$ of $t^{\alpha}$ on $[\delta, 1]$ for various $\alpha$, $\delta$, and $k$. The results are presented in 160 tables containing the coefficients of $P_k(t)$ and $Q_k(t)$, the zeros and the poles of $r_{k,\alpha}(t)$, the extremal point of the error $t^\alpha - r_{k,\alpha}(t)$, the representation of $r_{k,\alpha}(t)$ in terms of partial fractions, etc. Moreover, we provide links to the files with the data that characterize $r_{k,\alpha}(t)$ which are available with enough significant digits so one can use them in his/her own computations.

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.