Neumann fractional diffusion problems: BURA solution methods and algorithms

Let us consider the non-local problem −Lαu=f, α∈(0,1), L is a second order self-adjoint elliptic operator in Ω⊂Rd with Neumann boundary conditions on ∂Ω. The problem is discretized by finite difference or finite element method, thus obtaining the linear system Aαu=f, A is sparse symmetric and positive semidefinite matrix. The proposed method is based on best uniform rational approximations (BURA) of degree k, rα,k, of the scalar function tα, t∈[0,1]. Then, the approximate solution ur of the fractional power linear system is defined as u≈ur=λ2−αrα,k(λ2A†)f, where λ2 is the first positive eigenvalue of A, and A† stands for the Moore–Penrose pseudo inverse of A. The BURA method reduces the non-local problem to solution of k linear systems with matrices A+diI, di>0, i=1,…,k. An exponential convergence rate with respect to k is proven. The error estimates are robust with respect to the condition number of A in the subspace orthogonal to the constant vectors. The algorithm has almost optimal computational complexity assuming that optimal iterative solvers are applied to the auxiliary sparse linear systems. The first group of numerical tests illustrates in detail the obtained theoretical results. Finite difference discretization of a model 2D problem is used for this purpose. The second part of numerical tests demonstrates the applicability of BURA methods for 3D problems in domains with general geometry. Linear finite elements on unstructured tetrahedral meshes with local mesh refinement are used in the presented large-scale experiments, confirming also the almost optimal computational complexity.