Alternative numerical computation of one-sided Levy and Mittag-Leffler distributions

We consider here the recently proposed closed form formula in terms of the Meijer G-functions for the probability density functions $g_\alpha(x)$ of one-sided L\'evy stable distributions with rational index $\alpha=l/k$, with $0<\alpha<1$. Since one-sided L\'evy and Mittag-Leffler distributions are known to be related, this formula could also be useful for calculating the probability density functions $\rho_\alpha(x)$ of the latter. We show, however, that the formula is computationally inviable for fractions with large denominators, being unpractical even for some modest values of $l$ and $k$. We present a fast and accurate numerical scheme, based on an early integral representation due to Mikusinski, for the evaluation of $g_\alpha(x)$ and $\rho_\alpha(x)$, their cumulative distribution function and their derivatives for any real index $\alpha\in (0,1)$. As an application, we explore some properties of these probability density functions. In particular, we determine the location and value of their maxima as functions of the index $\alpha$. We show that $\alpha \approx 0.567$ and $\alpha \approx 0.605$ correspond, respectively, to the one-sided L\'evy and Mittag-Leffler distributions with shortest maxima. We close by discussing how our results can elucidate some recently described dynamical behavior of intermittent systems.