Differential equations of fractional order have in recent years proven to be very useful tools for the mathematical modeling of many phenomena in science and engineering. Certain types of behavior that can be observed in reality can only be reproduced in a satisfactory way by such models, but not by their traditional counterparts. Unfortunately, the very property that is responsible for this high quality of fractional models, namely the nonlocality of the differential operators arising in these equations, also has the unpleasant consequence of substantially raising the computational complexity of associated numerical algorithms. This is particularly relevant in the case of a model that describes a process acting in a 2D or 3D environment. Therefore, it is not surprising to note that most papers devoted to such questions restrict their attention to models in one or, rarely, two space dimensions.

To overcome these limitations, Wang and Du have looked at a specific prototypical problem in three space dimensions, namely a space-fractional diffusion equation, and developed an efficient numerical scheme for its solution. This is a task that is of significant interest in its own right, and moreover it can be expected that a successful approach for handling this special problem will show a path to corresponding algorithms for related and more general classes of equations.

The fundamental concept for the fast solver is based on an innovative combination of well-known ideas. It starts with the shifted Grünwald-Letnikov method used for discretizing the fractional derivatives in each space dimension as suggested and proved to be stable and rapidly convergent by Meerschaert and Tadjeran [1]. These three 1D methods are then merged into a single 3D method via the alternating direction implicit (ADI) technique. In the final step, the authors then develop two different fast solvers for the matrix equation resulting from this procedure, an iterative method and one of multistep type. Given a discretization using N_x, N_y, and N_z points in the x, y, and z directions, respectively, and denoting N = N_x N_y N_z, the first approach gives rise to a method with a memory requirement of O(N) that uses only O(N log N) arithmetical operations in each time step; the second one has the same memory requirement and a slightly larger arithmetical cost of O(N log² N). In comparison, a traditional naive scheme would use O(N²) memory and an O(N²) operation count.