Fractional differential equations, which use derivatives of non-integer order, represent an appropriate mathematical model for describing memory and hereditary effects of various physical, biological or engineering complex systems. Due to the non-locality of the operators involved in them, analytical solutions can be found very rarely, therefore numerical methods are rapidly developed. In this article we collect classical discretisation methods and contemporary acceleration and optimisation strategies for providing an overview of the state-of-the-art numerical solvers for solving of the above mentioned type of equations. We describe time stepping schemes, which are based on finite differences and convolution-quadrature, finite elements and discontinuous Galerkin methods, spectral and operational matrix, meshless collocation, semi-analytical series and physics-informed neural networks and other learning-based surrogates. The main contribution of the presented work is the unified treatment of the deterministic and stochastic fractional cases, where the problem of the existing research gap in the scalability of the variable order solvers is addressed. A simple classification scheme is proposed in this work, which will allow researchers to classify existing solvers and methods by using their discretisation method and the way they manage the history of the equation. Therefore, researchers may avoid the known problems in implementing existing solvers and may have a better idea about the possible ways to develop new efficient and structure preserving fractional integrators.
Dozva et al. (Wed,) studied this question.