There are several numerical methods to find approximate solutions and, as such, they provide something like an approximate analytical solution. My co-authors and I used Laguerre polynomials, if I remember correctly. Perhaps a more useful (for analysis) approach is to use Fourier transforms, i.e., Fourier integral operators.
It turns out that it matters a lot what type of fractional operator you want to approximate, whether it is in space or in time. Generally, the approximation schemes are much more limited in the sense that the regularity of these type of equations is also quite low (in most cases, only Holder regular solutions exist; and/or cont. Differentiability is not valid all the way to the origin, etc; the Hadamard notion of well-posedness is broken for equations with a fractional in time derivative, for instance). Most numerical schemes that I had seen in the engineering literature (and sometimes in the mathematical literature) are wrong, since they almost always assume infinite regularity of the solutions (and thus perform numerical truncations consistent with those assumptions). Thus, the question of optimal regularity is crucial for development of numerical schemes. See for instance:
Book Fractional-in-Time Semilinear Parabolic Equations and Applications