A unified view of some numerical methods for fractional diffusion

In recent years, a number of numerical methods for the solution of fractional Laplace and, more generally, fractional diffusion problems have been proposed. The approaches are quite diverse and include, among others, the use of best uniform rational approximations, quadrature for Dunford-Taylor-like integrals, finite element approaches for a localized elliptic extension into a space of increased dimensions, and time stepping methods for a parabolic reformulation of the fractional differential equation. A systematic comparison, both theoretical and experimental, of these approaches has thus far been lacking. A main contribution of the present work is the observation that all approaches mentioned above can, in fact, be interpreted as realizing different rational approximations of a univariate function over the spectrum of the original (nonfractional) diffusion operator. While this is obvious for some of the methods, it is a new result in particular for extension-based and time stepping approaches. This observation allows us to cast all described methods into a unified theoretical and computational framework, which has a number of benefits. Theoretically, it enables us to develop new convergence proofs for several of the studied methods, clarifies similarities and differences between the approaches, suggests how to design new and improved methods, and allows a direct comparison of the relative performance of the various methods. Practically, it provides a single, simple to implement, efficient and fully parallel algorithm for the realization of all studied methods; for instance, this immediately yields a fast and memory-efficient way of realizing all tensor product extension methods and lets us parallelize the otherwise inherently sequential time stepping approach. Finally, we present a detailed numerical study comparing all investigated methods for various fractional exponents and draw conclusions from the results. The comparison is made fair by the central insight that the computational effort of all these methods depends only on a single parameter, the degree of the underlying rational approximation. As a point of comparison, we also test a simple rational approximation method based on a black-box direct rational approximation algorithm which performs very well in practice. (C) 2019 Elsevier Ltd. All rights reserved.