 # Highly Accurate Global Padé Approximations of Generalized Mittag-Leffler Function and its Inverse

The two-parametric Mittag-Leffler function (MLF), E_α,β, is fundamental to the study and simulation of fractional differential and integral equations. However, these functions are computationally expensive and their numerical implementations are challenging. In this paper, we present a unified framework for developing global rational approximants of E_α,β(-x), x>0, with { (α,β): 0 < α≤ 1, β≥α, (α,β) (1,1) }. This framework is based on the series definition and the asymptotic expansion at infinity. In particular, we develop three types of fourth-order global rational approximations and discuss how they could be used to approximate the inverse function. Unlike existing approximations which are either limited to MLF of one parameter or of low accuracy for the two-parametric MLF, our rational approximants are of fourth order accuracy and have low percentage error globally. For efficient utilization, we study the partial fraction decomposition and use them to approximate the two-parametric MLF with a matrix argument which arise in the solutions of fractional evolution differential and integral equations.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

In this paper, we consider the two-parametric MLF

 Eα,β(z)=∞∑k=0zkΓ(αk+β),Rα>0,β∈C,z∈C. (1)

This entire function generalizes the MLF of one-parameter, .

The function

plays a key role in the study and simulation of history-dependent evolution models that arise in many engineering and science areas such as flow in porous media, pattern recognition, rheology, anomalous diffusion, electric networks, etc. In particular, it is the cornerstone of the development of generalized exponential time differencing (GETD) schemes

Garrappa:2011 which extend the notion of exponential integrator Minchev:2005 to time-fractional problems.

Evaluation of with scalar arguments is very expensive and challenging. Although the series (1) converges analytically for all , it is not practical or may not be valid to use it computationally for . Consequently, different techniques for the evaluation of have been developed. Gorenflo, Loutchko and Luchko Gorenflo:2002 proposed an algorithm based on using appropriate techniques for different regions of . For small and large values, they used the series definition (1) and the asymptotic series at infinity, respectively. For the intermediate regions they used the integral representations. A similar approach has been followed by Hilfer and Seybold Hilfer:2006 . Garrappa Garrappa:2015 provided an approach based on the numerical inversion of the Laplace transform. For efficient implementation, he provided an algorithm for finding the optimal parabolic contour on the basis of the distance and strength of the singularities of the Laplace transform.

The evaluation of MLF with a matrix argument is still a tricky and tough task. Garrappa Garrappa:2018

developed an algorithm based on the similarity transform. This approach requires the evaluation of MLF and its derivatives for each eigenvalue, which is again obtained using the Laplace transform. Clearly, massive calculations will be required for large full matrices.

In summary, all existing algorithms for evaluating suffer from some drawbacks such as nontrivial software implementation, long CPU time especially when a fine error tolerance is imposed, overflow numbers, and catastrophic cancellations. Due to these computational complexities and the need for efficient matrix function evaluation, accurate and efficient approximations are imperative.

To the best of authors’ knowledge, there have been few studies about rational approximations of MLF. Freed et al.Freed:2002 developed a piecewise approximant for , , based on the truncated series representation for small values, the asymptotic expansion for large values, and a Padé type approximant for the intermediate values. For , , Starovoitov and Starovoitova Starovoitov:2007 analyzed Padé type approximants of the form , , and discussed their asymptotic rate of convergence on the compact unit disk as . Borhanifar and Valizadeh Borhanifar:2015 constructed a fourth order Padé approximant and used it to develop a numerical scheme for the time-space diffusion equation. Iyiola et al.Iyiola:2018real constructed a second order non-Padé type rational approximation for using real distinct poles (RDP). However, although the approximants in Borhanifar:2015 and Iyiola:2018real might be adequate for small values, they fail to account for the asymptotic power law behavior.

The global Padé approximation technique introduced by Winitzki Winitzki:2003 has been applied recently to construct rational approximations for the Mittag-Leffler function and its generalization. In this technique, rational approximations are constructed by matching them with selected combinations of the series definition and the asymptotic expansion. Atkinson and Osseiran Atkinson:2011 used this technique to construct a second-order rational approximation for . Later, Ingo et al.Ingo:2017 showed that the rational approximant in Atkinson:2011 is not satisfactory for close to one. Alternatively, they constructed a fourth-order global approximation for that behaves reasonably well for all values of . Zeng Zeng2015 extended this technique to construct a second-order global Padé approximant for . However, this approximation is not satisfactory for close to one, especially when and it is malfunctioning for , .

As for the inverse of MLF, Hilfer and Seybold Hilfer:2006 introduced the inverse of as the solution of the equation

 Lα,β(Eα,β(z))=z,z∈C, (2)

where is evaluated by solving the functional equation (2) numerically. They discussed the principal branch of the function and showed that it reduces to the principal branch of the logarithm function as when . Hanneken and Achar Hanneken:2014 proposed a finite series representation for but only for some values of and . Lately, approximations of have been introduced to overcome the difficulty of solving the functional equation (2). Atkinson and Osseiran Atkinson:2011 ; and Ingo et al.Ingo:2017 discussed the approximation of based on the inversion of their global Padé approximants. Similarly, Zeng and Chen Zeng2015 ; and Iyiola et al.Iyiola:2018real inverted their second order approximants of to obtain an approximation of .

Consequently, based on the current state of the literature, more accurate rational approximations of and its inverse are needed. Such approximations are expected to ease computation cost and yield accurate values globally. In this paper, we introduce a framework that unifies the notion of global Padé approximation for the two-parametric MLF. Moreover, we develop different types of fourth-order global rational approximations for , , for . We also discuss analytically and numerically the approximation errors. Furthermore, we present the partial fraction decomposition of the rational approximants together with its advantage in efficient implementation for matrix arguments. An algorithm for computing based on the inversion of our fourth order approximants is presented. All along, we demonstrate through numerical experiments and comparisons, that the new developed fourth-order approximants provide superior global approximations for and its special case .

This paper is organized as follows. Section 2 contains the unified framework for the global Padé approximation and the error analysis. In section 3, we discuss the second order global Padé approximant constructed by Zeng and Chen in Zeng2015 and the need for more accurate approximants. Section 4 contains the construction of our fourth order approximants. The partial fraction decomposition and the algorithms to compute the poles and weights are discussed in section 5. In section 6 we discuss the inverse MLF and its approximation through the inversion of our rational approximants. Graphical and numerical demonstrations of the performance of our approximants are presented in section 7. In section 8, we apply our approximants to the solutions of fractional differential and integral equations and systems.

The computations in this paper are performed using Matlab software on a dell laptop with a core i5 processor.

## 2 Global rational approximation for Eα,β(−x), x>0

In this section, we introduce and outline the construction of the global Padé approximation for , , for the cases

 A={(α,β):0<α≤1,β≥α,(α,β)≠(1,1)}. (3)

Our approach is based on the technique proposed by Winitzki in Winitzki:2003 . This technique relies on the asymptotic expansion given by the following theorem (Podlubny:1999a , Theorem 1.4).

###### Theorem 2.1.

Let , and , . Then for ,

 Eα,β(z)=−n∑k=1(z)−kΓ(β−αk)+O(|z|−(n+1)),as |z|→∞,n≥1. (4)

In particular, when , the series in (4) takes the form

 Eα,α(z)=−n−1∑k=1z−(k+1)Γ(−αk)+O(|z|−(n+1)),as |z|→∞,n≥2. (5)

As an abbreviation for the rest of the paper, the cases and are to be understood as sub-cases of (3).

### 2.1 Definition

We proceed by considering the function

 Eα,β(x)=sα,β(x)Eα,β(−x), (6)

with chosen so that the first term in the asymptotic expansion of is 1. It follows from (4) and (5) that

 sα,β(x)=⎧⎨⎩Γ(β−α)x,β>α,−Γ(−α)x2,β=α. (7)

This function admits the following behavior:

 Eα,β(x)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩a(x)+O(xm), as % x→0,m≥{2,β>α,3,β=α,b(x−1)+O(x−n), as x→∞,n≥{1,β>α,2,β=α, (8)

where, from (1),

 a(x)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩Γ(β−α)xm−2∑k=0(−x)kΓ(β+αk),β>α,−Γ(−α)x2m−3∑k=0(−x)kΓ(αk+α),β=α, (9)

and from (4)–(5),

 b(x−1)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩−Γ(β−α)xn∑k=1(−x)−kΓ(β−αk),β>α,Γ(−α)x2n∑k=1(−x)−(k+1)Γ(−αk),β=α. (10)

Note that when , then . We will see later (equation (18)) that in this case the asymptotic expansion in (8) does not contribute to the rational approximation of . Therefore, for our purposes, we always take .

Next, we introduce the following definition.

###### Definition 2.2.

Consider with . Let and be positive integers such that

 n>1,m≥{2, if β>α,3, if β=α,,m+n is odd. (11)

Then, the global Padé approximation, , of type for is defined as

 Rm,nα,β(x)=1sα,β(x)p(x)q(x),0<α≤1,β≥α,(α,β)≠(1,1), (12)

where and are polynomials of degree ,

 ν:=m+n−12≥1, (13)

such that for and

 p(x)q(x)=⎧⎨⎩a(x)+O(xm−ν), as x→0,b(x−1)+O(x−n), as x→∞. (14)

Next, we present the procedure for constructing .

### 2.2 Construction of Rm,nα,β(x)

Let and be as in (11). We seek a rational approximation of the form

 Eα,β(x)≈p(x)q(x)=p0+p1x+⋯+pνxνq0+q1x+⋯+qνxν, (15)

where is as in (13). This means that coefficients are to be determined.

Since as and , we can set

 pν=qν=1.

To find the other unknowns , we solve the system of linear equations obtained by satisfying the requirement (14) which takes the form

 p(x)−q(x)a(x) =O(xm), as x→0, (16) x−νp(x)−x−νq(x)b(x−1) =O(x−n), as x→∞. (17)

By expanding the left-hand side of (16), it follows that the coefficients of , , must vanish. As such, we obtain linear equations. Similarly, by expanding the left hand side of (17), the coefficients of

 x−1,x−2,…,x−(n−1) (18)

must vanish and we obtain another system of linear equations. Collectively, (16) and (17) yield a linear system of (= ) equations which are then solved for the unknowns. By inspection, we have

 {p0=0,β>α,p0=p1=0,β=α. (19)

Hence, for we solve () equations with unknowns, while for we solve (= ) equations with unknowns.

We will see later (Remark 2.3) that for controlling the approximation error we must have . Table (1) provides the order of approximations for the types with and is odd.

### 2.3 Approximation error

Consider the pointwise error

 em,nα,β(x):=Eα,β(−x)−Rm,nα,β(x),x>0. (20)

Equations (6), (7), (8), (12), and (14) yield the following orders. As , we have

 em,nα,β(x)=Eα,β(−x)−Rm,nα,β(x) =1sα,β(x){Eα,β(x)−p(x)q(x)} =1sα,β(x){a(x)+O(xm)−a(x)−O(xm−ν)} =1sα,β(x){O(xm)+O(xm−ν)} ={O(xm−ν−1),β>α,O(xm−ν−2),β=α. (21)

As , we have

 em,nα,β(x)=Eα,β(−x)−Rm,nα,β(x) =1sα,β(x){Eα,β(x)−p(x)q(x)} =1sα,β(x){b(x−1)+O(x−n)−b(x−1)−O(x−n)} =1sα,β(x)O(x−n) ={O(x−n−1),n>1,β>α,O(x−n−2),n>1,β=α. (22)
###### Remark 2.3.

By inspection of (21) and (22), one can observe the following.

• For reliable approximations of at small values, one should consider when and when . This is why, for example, is not a good approximation when .

• For large values of , the approximation error can be made arbitrary small by taking sufficiently large. However, this may not be the case since the asymptotic series (4) of the Mittag-Leffler function is divergent.

###### Remark 2.4.

In all the numerical experiments and comparisons throughout this paper, the term ”exact” values of MLF refers to the values computed using the routines discussed in Garrappa:2015 and Garrappa:2018 .

## 3 Second-order global Padé approximant R3,2α,β

For completeness, we provide here an overview of the second order global Padé approximant constructed by Zeng and Chen in Zeng2015 . The approximant is given by

 R3,2α,β(x)=1Γ(β−α)p1+xq0+q1x+x2,β>α, (23)

with

 p1 =cα,β[Γ(β)Γ(β+α)−Γ(β+α)Γ2(β−α)Γ(β−2α)], (24) q0 =cα,β[ Γ2(β)Γ(β+α)Γ(β−α)−Γ(β)Γ(β+α)Γ(β−α)Γ(β−2α)], q1 =cα,β[Γ(β)Γ(β+α)−Γ2(β)Γ(β−α)Γ(β−2α)], cα,β =1Γ(β+α)Γ(β−α)−Γ2(β),

and

 R3,2α,α(x)=αΓ(1+α)+2Γ(1−α)2Γ(1−2α)x+Γ(1−α)x2,0<α<1. (25)

As shown in Figures 1 and 2, the approximation could be reasonable for small values of , however, it is not adequate otherwise. Figure 1: Plots of R3,2α,1 for different values of α Figure 2: Plots of R3,2α,β for β=α=0.25 (left) and β=α=0.5 (right)

## 4 Fourth-order global Padé approximants

Fourth order global Padé approximants () correspond to the types with . They include the types , , and . As discussed in subsection 2.2, the approximation for takes the form

 Rm,nα,β(x)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩1Γ(β−α)p1+p2x+p3x2+x3q0+q1x+q2x2+q3x3+x4,β>α,−1Γ(−α)^p2+^p3x+x2^q0+^q1x+^q2x2+^q3x3+x4,β=α. (26)

The unknown coefficients are obtained by applying (16) and (17). Below we present the systems for these coefficients for the different types.

### 4.1 Coefficients of R5,4α,β

For , the coefficients satisfy the system

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣100−Γ(β−α)Γ(β)000010Γ(β−α)Γ(β+α)−Γ(β−α)Γ(β)00001−Γ(β−α)Γ(β+2α)Γ(β−α)Γ(β+α)−Γ(β−α)Γ(β)0000Γ(β−α)Γ(β+3α)−Γ(β−α)Γ(β+2α)Γ(β−α)Γ(β+α)−Γ(β−α)Γ(β)1000−1Γ(β−α)Γ(β−2α)−Γ(β−α)Γ(β−3α)01000−1Γ(β−α)Γ(β−2α)001000−1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝p1p2p3q0q1q2q3⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝000−1−Γ(β−α)Γ(β−4α)Γ(β−α)Γ(β−3α)−Γ(β−α)Γ(β−2α)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (27)

For the coefficients satisfy the system

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣10Γ(−α)Γ(α)00001−Γ(−α)Γ(2α)Γ(−α)Γ(α)0000Γ(−α)Γ(3α)−Γ(−α)Γ(2α)Γ(−α)Γ(α)0000−1−Γ(−α)Γ(−2α)01000−1Γ(−α)Γ(−2α)01000−1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝^p2^p3^q0^q1^q2^q3⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝00−100−Γ(−α)Γ(−2α)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (28)

### 4.2 Coefficients of R6,3α,β

For , the coefficients satisfy the system

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣100−Γ(β−α)Γ(β)000010Γ(β−α)Γ(β+α)−Γ(β−α)Γ(β)00001−Γ(β−α)Γ(β+2α)Γ(β−α)Γ(β+α)−Γ(β−α)Γ(β)0000Γ(β−α)Γ(β+3α)−Γ(β−α)Γ(β+2α)Γ(β−α)Γ(β+α)−Γ(β−α)Γ(β)000−Γ(β−α)Γ(β+4α)Γ(β−α)Γ(β+3α)−Γ(β−α)Γ(β+2α)Γ(β−α)Γ(β+α)01000−1Γ(β−α)Γ(β−2α)001000−1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝p1p2p3q0q1q2q3⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝000−1Γ(β−α)Γ(β)Γ(β−α)Γ(β−3α)−Γ(β−α)Γ(β−2α)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (29)

For , the coefficients satisfy the system

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣10Γ(−α)Γ(α)00001−Γ(−α)Γ(2α)Γ(−α)Γ(α)0000Γ(−α)Γ(3α)−Γ(−α)Γ(2α)Γ(−α)Γ(α)000−Γ(−α)Γ(4α)Γ(−α)Γ(3α)−Γ(−α)Γ(2α)Γ(−α)Γ(α)1000−1Γ(−α)Γ(−2α)01000−1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝^p2^p3^q0^q1^q2^q3⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝00−10Γ(−α)Γ(−3α)−Γ(−α)Γ(−2α)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (30)

### 4.3 Coefficients of R7,2α,β

For , the coefficients satisfy the system

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣100−Γ(β−α)Γ(β)000010Γ(β−α)Γ(β+α)−Γ(β−α)Γ(β)00001−Γ(β−α)Γ(β+2α)Γ(β−α)Γ(β+α)−Γ(β−α)Γ(β)0000Γ(β−α)Γ(β+3α)−Γ(β−α)Γ(β+2α)Γ(β−α)Γ(β+α)−Γ(β−α)Γ(β)000−Γ(β−α)Γ(β+4α)Γ(β−α)Γ(β+3α)−Γ(β−α)Γ(β+2α)Γ(β−α)Γ(β+α)000Γ(β−α)Γ(β+5α)−Γ(β−α)Γ(β+4α)Γ(β−α)Γ(β+3α)−Γ(β−α)Γ(β+2α)001000−1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝p1p2p3q0q1q2q3⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝000−1Γ(β−α)Γ(β)−Γ(β−α)Γ(β+α)−Γ(β−α)Γ(β−2α)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (31)

For , the coefficients satisfy the system

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣10Γ(−α)Γ(α)00001−Γ(−α)Γ(2α)Γ(−α)Γ(α)0000Γ(−α)Γ(3α)−Γ(−α)Γ(2α)Γ(−α)Γ(α)000−Γ(−α)Γ(4α)Γ(−α)Γ(3α)−Γ(−α)Γ(2α)−Γ(−α)Γ(α)00Γ(−α)Γ(5α)−Γ(−α)Γ(4α)Γ(−α)Γ(3α)−Γ(−α)Γ(2α)01000−1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝^p2^p3^q0^q1^q2^q3⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝00−10Γ(−α)Γ(−3α)−Γ(−α)Γ(−2α)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (32)
###### Remark 4.1.

Although the type of the fourth order global Padé approximant of by Ingo et al.Ingo:2017 is not given, based on their construction, we expect the type to be a special case of one of the approximants above when .

## 5 Partial fraction decomposition

Partial fraction decomposition provides an efficient form for evaluating rational functions. In a recent work by Bertaccini:2019 , the efficiency of using partial fraction decomposition for computing functions of matrices is discussed. This efficiency is indisputable when the poles are complex conjugates and the argument is a matrix.

Unlike the the Padé approximations for the exponential function, the poles of are functions of and . Fortunately, through direct calculations, one can show that for most , the poles of are complex conjugates. Next, we explore the partial fraction decomposition of these approximations.

### 5.1 Decomposition of second-order global Padé approximant

 R3,2α,β(x)=c1x−r1+c2x−r2,

where

 r1=−q1+√q21−4q02, r2=−q1−√q21−4q02,

and

 c1=p1−r1r2−r1, c2=p1−r2r1−r2.

We can verify numerically that for we have which imply that and . As a result, we can write

 R3,2α,β(x)=2R[c1x−r1]. (33)

### 5.2 Decomposition of the fourth-order global Padé approximants

The partial fraction decomposition for , , takes the form

 Rm,nα,β(x)=c1x−r1+c2x−r2+c3x−r3+c4x−r4. (34)

Empirically, for , these poles are complex conjugates. If we let , , , and , then the partial fraction decomposition can be written as

 Rm,nα,β(x)=2R[c1x−r1]+2R[c2x−r2]. (35)

Computing of the poles and weights is outlined in the following algorithm.

## 6 Inverse Mittag-Leffler Function

The invertibility of , , follows from the complete monotonicity property of . As shown in Gorenflo:2014 , this function is completely monotone if and only if and . Since and , then for