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\mathrm{/}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.