Rational Approximations for the Oscillatory Two-Parameter Mittag–Leffler Function

: The two-parameter Mittag–Leffler function E α , β is of fundamental importance in fractional calculus, and it appears frequently in the solutions of fractional differential and integral equations. However, the expense of calculating this function often prompts efforts to devise accurate approximations that are more cost-effective. When α > 1, the monotonicity property is largely lost, resulting in the emergence of roots and oscillations. As a result, current rational approximants constructed mainly for α ∈ ( 0,1 ) often fail to capture this oscillatory behavior. In this paper, we develop computationally efficient rational approximants for E α , β ( − t ) , t ≥ 0, with α ∈ ( 1,2 ) . This process involves decomposing the Mittag–Leffler function with real roots into a weighted root-free Mittag–Leffler function and a polynomial. This provides approximants valid over extended intervals. These approximants are then extended to the matrix Mittag–Leffler function, and different implementation strategies are discussed, including using partial fraction decomposition. Numerical experiments are conducted to illustrate the performance of the proposed approximants.


Introduction
We are concerned with the approximation of the two-parameter Mittag-Leffler function (MLF) E α,β , defined by This entire function generalizes the MLF of one parameter, E α = E α,1 , and contains several well-known functions as special cases.In particular, E 1 is the exponential function, and E 2 (−z 2 ) and zE 2,2 (−z 2 ) are the cosine and sine functions, respectively, among others.For some surveys on the Mittag-Leffler functions, see, for example, [1,2].
The Mittag-Leffler function E α,β (z) arises frequently in the solutions of many physical problems described by differential and/or integral equations of fractional order.In the case α ∈ (1, 2), the MLF appears naturally in the solutions of fractional diffusion-wave equations, fractional differential equations for motion, and fractional plasma equations (see [3][4][5][6][7]).
Computing the MLF E α,β is usually challenging (devising and implementing suitable algorithms) and expensive (computation time).Although the series (1) converges for all z ∈ C, it is impractical or ineffective to use it computationally for |z| ≥ 1 because the series converges very slowly.As a result, various methods have been developed to evaluate the MLF.In [8], an algorithm based on the location of the argument z in the complex plane was developed, where for large |z| values, the asymptotic series as |z| → ∞ was used; for |z| < 1, the series definition (1) was used; and for values in intermediate regions, the integral representation was used.In [9], an approach based on numerical inversion of the Laplace transform was proposed.However, these approaches are computationally time-consuming (see [10,11] for some CPU time comparisons).
Various rational approximations have been sought to efficiently and accurately approximate E α,β .Some of them are based solely on the series definition and provide accurate approximations for small argument values (see, for example, [12][13][14]).Others are based on the global Padé approach [15], in which a hybrid of the local series definition and the asymptotic series representation is used.These global approximations lead to approximants that are accurate over a wide range of arguments (see [10,11,[16][17][18]).
The existing rational approximants are effective when the MLF is completely monotone, which is the case when α ∈ (0, 1) and β ≥ α.However, when α > 1, as discussed in [19][20][21], for some α-β combinations, E α,β is often oscillatory with multiple real zeros.In these situations, the rational approximants might fail to trace the oscillation profile and match the large zeros.
In this paper, we consider in the αβ-plane the strip The regions where E α,β is monotone and where it is oscillatory are classified by studying the zeros of its derivative.Additionally, we introduce new rational approximants that sufficiently capture more of the MLF roots and can effectively trace their oscillations.This is achieved by decomposing E α,β with real roots into another weighted Mittag-Leffler function devoid of roots, accompanied by a polynomial.Based on this decomposition, we construct new rational approximants as a sum of the global Padé approximant and a polynomial.Using this approach, we are able to sufficiently capture many roots of E α,β by choosing the appropriate degree of the polynomial.Further, we generalize these approximants to the matrix MLF and discuss different approaches in which the approximant can be implemented for matrix arguments.Apart from the computation cost, all approaches yield close approximations.Numerical experiments are presented to illustrate the performance of our approximants.Considering the MLF with scalar arguments, we observe the limitations of the existing global Padé approximants in approximating the oscillatory MLFs.It is then demonstrated that the new approximants introduced here are able to improve the situation, thereby providing sufficiently accurate approximations.For the matrix MLF, we consider four approaches for computing the rational approximants and compare them to the matrix MLF algorithm developed in [22].Using the values obtained from this algorithm as a reference, it is shown that our approximants provide close values with significantly less computation time.
The rest of this paper is organized as follows.A characterization of the oscillatory behavior of the MLF is given in Section 2. In Section 3, the derooting representation is discussed.Global Padé approximants of the MLF with scalar and matrix arguments are discussed in Sections 4 and 5, respectively.Some applications to the evaluation of solutions of fractional evolution equations are discussed in Section 6.

Monotonicity and Oscillatory Properties
In this section, we characterize the oscillatory behavior of For the discussion that follows, the following formula for the derivative [1] plays a fundamental role: The function E α,β (−t), t ≥ 0, is known to be completely monotone for 0 < α ≤ 1, β ≥ α.However, this monotonicity is not necessarily preserved when (α, β) ∈ Ω, since either the function or its derivative could have roots.On the other hand, it was established in [20] that E α,β has at most a finite number of roots when (α, β) ∈ Ω.Thus, E α,β (−t) is monotone for sufficiently large t.Furthermore, it was shown in [20] that there exists only one zero when α is sufficiently close to 1, while the number of zeros increases as α increases toward 2.
In Figure 1, the curve ϕ is the boundary given in Table 1 in [20] such that E α,β has a finite number of real roots when (α, β) is below it and none above it.However, the function in the region below ϕ(α) differs with respect to the number of roots.Typically, the number of roots increases as α approaches 2. Similarly, we constructed the boundary ψ for the derivative.The corresponding data are presented in Table 1.The monotonicity and oscillatory behavior of E α,β are characterized by the different regions marked in Figure 1 and described in Table 2.It follows from ( 4) that E α,β (−t) is monotonically decreasing for (α, β) in regions (D) and (F).In regions (B) and (C), although the function has no real roots, it could have a finite number of oscillations due to the roots of the derivative.Remark 1.Since E α,β (−t) → 0 as t → ∞, it is obvious that each root of it is followed by at least one root of its derivative.This is consistent with the inequality ϕ(α) < ψ(α).For the region below the line β = α + 1, which includes the regions (A), (B), and (D), the following decomposition is presented in [19], where the function f α,β (−t) is asymptotically approaching zero as t → ∞ and g α,β (−t) is an oscillatory function.
In region (B), although E α,β (−t) has no roots as it decays to zero, it undergoes oscillations due to the roots of its derivative, as demonstrated in Figure 2.This behavior is also consistent with the oscillatory nature of g α,β there.On the other hand, in region (D), in which E α,β is monotone, the behavior of the oscillatory function g α,β fades away and f α,β (−t) dominates.

Derooting Decomposition
A key step in accurately approximating E α,β (−t) when 1 < α < 2 is to capture its oscillatory behavior and roots.However, in region (A), in general, the number of roots could exceed those of the rational approximants.This makes these approximants valid only for small arguments.
A derooting approach can be used to extend the validity of the rational approximants to large intervals.For this purpose, we use the following recursive identity [24] where In this identity, E α,β is decomposed into another MLF with a shifted second parameter and a polynomial of degree r − 1.This allows us to replace the parameters (α, β) in region (A) with parameters (α, β + αr) above the boundaries ϕ(α) or above ψ(α).Consequently, we obtain a function E α,β+αr that has no roots or could even be monotone.Therefore, the approximation of E α,β (−t), 1 < α < 2, with a finite number of real roots can be replaced by approximating another one that has no roots or one that is monotone.

Rational Approximation
Based on the global rational approximation technique introduced for transcendental functions in [15], a variety of global Padé approximants for E α,β (−t), 0 < α < 1, have been developed and implemented [10,11,[16][17][18].These approximants are still valid for E α,β (−t), 1 < α < 2, over intervals in which the number of its roots does not exceed the number of approximant roots.This renders these approximants inadequate when approximating the MLF over an interval in which it has multiple roots, as is the case when (α, β) is in region (A) in Figure 1.
In general, the global Padé approximation R m,n α,β (t), where m and n are positive integers, for E α,β (−t) takes the form [10,11] The construction of the global Padé approximation depends on the series definition (1) of the MLF as well as its asymptotic expansion, given by the following theorem: It was shown in [10] that the approximant R m,n α,β has the asymptotic error In the rest of this section, we consider the global Padé approximants R 7,2 α,β and R 13,4 α,β , developed in [10,11], and examine their performance when 1 < α < 2.Then, we describe how the derooting decomposition leads to approximants with an increased number of roots, enabling them to better approximate the roots of the MLF.

Remark 2.
As mentioned in [11], it is important to highlight that systems (12) and (13), for determining the coefficients of the approximants R 7,2 α,β and R 13,4  α,β , exhibit relatively high condition numbers.Consequently, solving these systems necessitates careful attention to avoid accuracy loss due to possible ill-conditioning.In the implementations, we make use of the MATLAB function "linsolve" or the corresponding Python function "scipy.linalg.solve",which can handle this range of condition numbers.
Computationally, it is observed that R 13,4  α,β has at most three real roots, while R 7,2 α,β has exactly one real root.Consequently, in approximating some of the oscillatory MLFs, each approximant is only able to trace the function up to its largest root.To illustrate the accuracy and limitations of both approximants when 1 < α < 2, we compare them with the MATLAB routines ml [9] and ml_matrix [22] for the scalar and matrix arguments, respectively, as our reference values.Other implementations of the MLF, such as the MATLAB function mlf [26] and the Mathematica MittagLefflerE function, can also be used.However, these functions are not available for Matrix arguments, which are of primary interest in the current work.As such, we adopt those in the aforementioned references for our comparisons.
In Figure 3, two examples from region (A) are shown.In plot (a), the MLF has only one root, which is typically the case for 1 < α < 1.4.As observed, R 13,4  α,β provides a good approximation in this case.In plot (b), the MLF has many roots, which is typically the case as α approaches 2. As expected, each approximant fails beyond its largest root to capture the oscillations of the MLF.In Figure 4, typical cases from regions (B) and (C) are shown.The corresponding MLFs in both regions are root-free but are oscillatory.As shown, R 13,4  α,β is sufficiently accurate for an extended interval.However, eventually, it takes negative values, as shown in plot (b).
In Figure 5, typical cases from regions (D) and (F) are illustrated.In both regions, the MLF is monotone and globally positive.As can be clearly seen, R 13,4  α,β is a good approximation for an extended interval.As general guidelines, we have the following: • When (α, β) is in region (A), the existing global Padé approximants are not effective over extended intervals due to the existence of roots.• When (α, β) is in region (B) or (C), more accurate approximants should be used due to the oscillatory behavior of the MLF.• When (α, β) is in region (D) or (F), the approximant R 13,4 α,β is sufficiently accurate since the MLF is globally monotone.

Rational Approximation for Oscillatory MLFs
As observed in the above subsection, when (α, β) lies below the boundary ψ(α) in Figure 1, then more accurate and multi-root approximants are essential for extended intervals.Next, we propose a class of rational approximants that have the desirable performance.
As demonstrated in the previous subsection, this shifting of parameters would enhance the performance of R m,n α,β+αr since in this case, E α,β+αr (−t) is root-free and non-oscillatory.Furthermore, r could be chosen large enough so that the approximant R m,n,r α,β captures the desired number of roots and the oscillations of E α,β over an extended interval.

Asymptotic Behavior of the Approximation Error
We define the error e m,n,r α,β (t . By subtracting the decompositions ( 14) and ( 16), and then using the asymptotic error (9), we can determine the asymptotic behavior: Note that although the exponent (−n − 1 + r) in the error term (17) becomes positive as t → ∞ when r > n + 1, large values of r increase the number of roots in the modified approximant R m,n,r α,β , which improves its capability to capture the oscillations of the MLF over an extended interval.Consequently, the actual error starts to grow only after the final root of the approximant.
The performance of the generalized approximant R 13,4,r α,β when (α, β) is in regions (A) and (C) is demonstrated in Figures 6 and 7, respectively, where e m,n,r α,β (t As can be observed, the interval of the approximation can be extended by increasing r.
The generalized approximant R 7,2,r α,β has the same features, although it is not as accurate as R 14,3,r α,β , as shown in Figure 8.

Approximation of the Matrix Mittag-Leffler Function
In this section, we study the approximation of the matrix MLF, defined by [27] Using the definition of rational functions of matrices (see [28]), we define the global Padé approximation of E α,β (A) as for some appropriate matrix A, where p and q denote the numerator and denominator, respectively, of Γ(β − α)R m,n α,β .Furthermore, by extending the decomposition (6) to matrix arguments, we define the rational approximant where P r−1 α,β is the polynomial given by ( 7).Several approaches exist for extending a scalar rational function to matrix arguments.For a comprehensive overview, we recommend consulting [28].Next, we discuss some approaches for implementing (19) and compare their accuracies and computation times.

1.
Linear system approach A straightforward approach is to evaluate p(A) and q(A) using efficient methods for the evaluation of matrix polynomials, such as the Paterson-Stockmeyer method and Horner's method (nested multiplication) [28].Additionally, the powers A 2 , A 3 , . . ., can be precomputed and used in the computation of both p(A) and q(A) to minimize the overall cost.It is noteworthy that research in the field of matrix polynomial evaluation is increasingly active.In particular, a new family of methods for evaluating matrix polynomials, which are more efficient than the established Paterson-Stockmeyer method, was proposed in [29].This area could be a subject of future research for us, as this section concentrates on introducing general techniques for approximating the MLF matrix using rational approximation, aiming for a general comparison.Using this approach, the approximant R m,n α,β (A) is obtained by solving the matrix system q(A)R m,n α,β (A) = p(A).This requires solving N systems for an N × N matrix.

2.
Partial fraction approach Partial fraction decomposition is known to provide an efficient form for evaluating rational functions.For the global Padé approximants, it was discussed in [10,11] that these approximants have complex conjugate roots, which can contribute to efficient implementation.As an example, the approximant R 13,4 α,β (t), admits the partial fraction decomposition where {c 1 , c 2 , c 3 , c 4 } and {s 1 , s 2 , s 3 , s 4 } are the non-conjugate residues and poles, re- spectively.This expression can be simplified as So, for a matrix argument A, the approximant can be calculated as where I is the identity matrix.For example, in the implementation using the partial fraction approach, the partial fraction decompositions (22) are used to compute the matrix-vector products as outlined below.
For a given square matrix A and a vector v, the matrix-vector product E α,β (−A)v is computed using (20) as where the term R m,n α,β+αr (A)v, is computed using the partial fraction decomposition approach.From (23), we use the approximation

By solving the systems
we obtain where x j , are the solutions of the linear systems (24).As such, no matrix inversion is required.Moreover, using the partial fraction decomposition of the rational approximations is considered an efficient and cost-effective alternative.As mentioned in [11], the poles and residues in the partial fraction decomposition (22) rely on α and β.Therefore, their computation can be achieved independently of the argument A.
The Python or MATLAB "residue" function could be used for this purpose.

3.
Matrix diagonalization approach When the matrix argument A is diagonalizable, a scenario frequently encountered in matrices derived from the semi-discretization of partial differential equations, then the factorization A = ZDZ −1 could be considered, where D is the diagonal matrix containing the eigenvalues, and the columns of Z are the corresponding eigenvectors.In this case, the matrix MLF can be computed as [28,30] Accordingly, the approximant R m,n α,β can be computed as By employing this approach, we only need to calculate the rational approximation for scalar arguments and perform matrix multiplications, which considerably reduces the computational time.We note that using the diagonalization approach is advisable only under the condition that the matrix is guaranteed to be well conditioned [28].
Remark 3. The above approaches for computing the MLF matrix are discussed for completeness and general comparisons.However, they are mainly developed to compute the matrix-vector products involved in the implementation of numerical methods for solving fractional differential equations.To the best of our knowledge, computing matrix-vector products involving the MLF remains a fundamental challenge, as current algorithms for evaluating the MLF of scalar or matrix arguments [8,9,22] require considerable implementation cost.Therefore, developing such rational approximants is essential, as they have been proven to offer comparable accuracy and efficient approximations at a low cost.
For experimental purposes, we evaluate the performance of the approximant (16) when applied to the Redheffer matrix of size 100 × 100, comparing it with the reference values obtained by the ml_matrix function.In Figure 9, a color map shows the componentwise values and relative errors for R 13,4  α,β , with (α, β) = (1.9,1), calculated using partial fractions.Table 3 presents a comparison of the absolute errors, relative errors, and runtimes for the different techniques in computing R 13,4  α,β .The absolute and relative errors are given by , respectively.

Applications and Numerical Experiments
Here, we present some applications of the general Padé approximant (16) related to fractional oscillation equations, where MLFs arise naturally as solutions.Numerical experiments are included to highlight the efficiency and accuracy of these approximants.

Application: Fractional Plasma Oscillations
The fractional plasma oscillation model, described in [31], is given by where c D α denotes the Caputo fractional derivative of order α ∈ (1, 2), defined by The constant A is the fractional electron plasma frequency and f (t) is the electric field.We consider the model with a static electric field and with no electric field.

Fractional Plasma Oscillation Model with a Static Electric Field
Consider the special case of problem (25), The exact solution for this initial value problem is In Figure 10, a comparison of the approximations for (26) when α = 1.2 is provided.Considering the three terms in (26), the pair (α, β) with β = 1 falls within region (A) of the phase diagram in Figure 1, where the MLF is oscillatory.Meanwhile, the pairs (1.2, 2) and (1.2, 2.2) fall within regions (D) and (F), respectively, where the MLF has no real roots and no oscillations.A similar comparison for α = 1.9 is shown in Figure 11.In this case, the terms E 1.9,1 and E 1.9,2 correspond to (α, β) in parts of region (A) where the MLF is highly oscillatory and has many real roots.
To illustrate the computational efficiency, Table 4 shows that R 13,4 α,β and R 13,4,r α,β with r = 2, 5, 8 exhibit significantly lower computation costs (in terms of time) compared to the ml function.Moreover, while the modified approximant R 13,4,8   α,β yields more accurate approximations than R 13,4  α,β in comparison to the reference values, they have comparable computation times.
The exact solution of this problem is given by A comparison of the approximations of ( 27) when α = 1.9 is provided in Figure 12.Unlike the static electric field case, it can be observed that over an extended time interval, the solution of the fractional oscillation equation has more oscillations around zero.To sufficiently capture these oscillations on the interval [0, 20], the approximant R 13,4,15   α,β can be used.Clearly, as the values of r increase, the number of roots in the rational approximation also increases.As a result, the precision of the approximants improves, and their ability to capture a greater number of roots and oscillations across a relatively broader range is enhanced.Consequently, the accuracy of the approximants improves, enhancing their capability to capture more roots and oscillations over a relatively extended interval.27), with α = 1.9.

Concluding Remarks
In this paper, a characterization of the oscillatory and monotone behavior of the two-parameter Mittag-Leffler function is given.Furthermore, generalized rational approximants over extended intervals are developed.These approximants are based on decomposing a Mittag-Leffler function that possesses real roots into a combination of a weighted rootless Mittag-Leffler function and a polynomial.The approximants have good tracking capabilities of the roots and oscillations over extended intervals.To show the efficiency of these approximants, numerical experiments are presented alongside various applications.Furthermore, to demonstrate computational efficiency, comparisons with the ml and ml_matrix MATLAB functions are included.Such approximants play an effective role in the implementation of numerical methods for fractional oscillation equations.

Figure 2 .
Figure 2. Typical oscillations in region (B) about the roots (marked by the red lines) of the derivative.

Table 2 .
Number of roots in the regions in Figure1.