Mostly known as the man because of whom mathematicians are not entitled to the Nobel prize, the Swedish mathematician Gösta Mittag-Leffler, the man who – the legend has it – had an affair with Nobel’s wife, is little known for his namesake function (for the record, Nobel never had a wife). The Mittag-Leffler function arises in several applications in physics and mathematics, chiefly in the solution of fractional-order differential equations – a well-studied special type of integrodifferential equations [1]. The function is given in the form of an infinite convergent series. A straightforward way to compute then would be to sum up the terms of this series until, eventually, the sum does not change. This, however, would be neither fast not accurate. In this article we are discussing about numerical methods for the computation of the Mittag-Leffler function.

The (one-parameter) Mittag-Leffler function is defined as

Its two-parameter variant is given by

that is, , and the three-parameter ML function, also known as Prabhakar function, is given by

What is particularly problematic about this definition when it comes to computing the value of using the above formula is that often many terms of the sum would be required to achieve reasonable convergence.

At the same time, for large values of , the terms in the numerator of the -term in the definition becomes quite large as increases and even larger becomes the value of . For example, for and , the fiftieth term of the sum has numerator and denominator (computed using Octave). Despite the fact that we cannot be sure that such high numbers are accurate, it is clear that the direct summation approach is impractical. Additionally, it is not always clear when to terminate the summation so as to achieve some desired accuracy.

There are several alternative ways to compute the Mittag-Leffler function and, not surprisingly, commence from a fractional-order differential equation whose solution is given by a Mittag-Leffler function. To understand the analogy, this is like computing the exponential function by solving a linear system by a numerical method.

One of the first works to address these computational challenges was a 2005 article by Hilfer and Seybold [12].

The Gorenflo-Loutchko-Luchko method is the currently the most popular one for the computation of the two-parameter ML function [3]. According to this approach, for all with modulus smaller than 1 and for , we may compute the ML function as where can be easily computed explicitly so that the approximation error is, in modulus, smaller than a desired . In other cases, the ML function is written in as an indefinite integral which is approximated up to desirable accuracy with a definite one. Nevertheless, this approach is of questionable accuracy in certain regions of the complex plane [4].

A 2015 paper by R. Garrappa proposed an alternative method for faster and more accurate computations of the 3-parameter ML function based on numerical methods for the inverse Laplace transform (NILT) [4]. The paper is accompanied by a MATLAB implementation which seems to be very efficient. See also a 2013 paper on the same topic where the ML function is evaluated on the real-line only [11].

NILT methods have been quite popular for solving fractional-order equations; see for example [5-7] and for fractional PDEs in [8]. Naturally, any numerical method for solving fractional-order differential equations can be employed to compute the Mittag-Leffler function.

Another paper which attracted my attention was the work of Zheng and Chen who propose a global Padé approximation for the 2-parameter ML function [2]. The proposed approximation is a rational function, but is valid only for , however, this is not a severe limitation because of the recursive formula

where .

Another approximation has been proposed by Valério and Machado in [9] for computing the function which arises naturally as the solution of simple fractional-order differential equations. The authors use low and high time asymptotic expansions which then they combine in a common formula which exhibits approximation errors typically lower than 1%.

**List of related software**

- It seems that the most popular MATLAB implementation of the Mittag-Leffler function by I. Podlubny here. The implementation is based on [3].
- A MATLAB implementation of [4] is available on MATLAB’s file exchange. In contrast to #1 this allows for computations involving the one, two and three parameters ML functions and works with complex arguments as well.
- The one and two parameter functions are available in wolfram and Mathematica.
- Although I have not tried this one, there is an R-package called mittagleffler which you can download/clone from github. It is based on a F90 code by D. Verrotta which is, in turn, based on Podlubny’s MATLAB code.
- A well-documented and robust MATLAB implementation of a numerical method for the inverse Laplace method is available on MATLAB’s file exchange.

**References**

[1] H.J. Haubold, A.M. Mathai, R.K. Saxena, *Mittag-Leffler Functions and Their Applications*, 2009, arXiv: http://arxiv.org/abs/0909.0230

[2] C. Zheng and Y. Chen, *Global Padé approximations of the generalized Mittag-Leffler function and its inverse*, 2015, arXiv: http://arxiv.org/abs/1310.5592v3

[3] R. Gorenflo, J. Loutchko and Y. Luchko, *Computation of the Mittag-Leffler function and its derivative, *Fractional Calculus and Applied Analysis 5(4), 2002, available here.[4] R. Garrappa, Numerical evaluation of two and three parameter Mittag-Leffler functions, 2015, arXiv: http://arxiv.org/abs/1503.06569

[5] H. Sheng, Y. Li and Y. Cheng, *Application of numerical inverse Laplace transform algorithms in fractional calculus*, Journal of the Franklin Institute 348(2), 2001, pp. 315-330.

[6] S. Gupta, D. Kumar and J. Singh, *Numerical study for systems of fractional differential equations via Laplace transform*, Journal of the Egyptian Mathematical Society 23(2), 2015, pp.256–262

[7] D.W. Brzeziński and P. Ostalczyk, *Numerical calculations accuracy comparison of the Inverse Laplace Transform algorithms for solutions of fractional order differential equations*, Nonlinear dynamics 84(1), 2016, pp. 65-77.

[8] M. Javidi and B. Ahmad, *Numerical solution of fractional partial differential equations by numerical Laplace inversion technique*, Advances in difference equations 2013(375), 2013.

[9] D. Valério and J.T. Machado, *On the numerical computation of the Mittag-Leffler function*, Commun Nonlin Sci Numer Simulat 19, 2014, pp.3419-3424.

[10] A. Saa and R. Venegeroles, Alternative numerical computation of one-sided Lévy and Mittag-Leffler distributions, 2011, arXiv: http://arxiv.org/abs/1106.0531v2

[11] R. Garrappa, M. Popolizio, *Evaluation of generalized Mittag–Leffler functions on the real line*, Advances in Computational Mathematics 39(1), 2013, pp. 205-225.

[12] R. Hilfer, H.J. Seybold, Computation of the generalized Mittag-Leffler function and its inverse in the complex plane, Integral Transforms and Special Functions 17(9), 2006.