Computing the matrix Mittag-Leffler function with applications to fractional calculus

1 Introduction

• Accurate scalar Mittag–Leffler evaluation is itself challenging and expensive, with efficient algorithms developed only relatively recently.The authors cite.Source: §1 · p. 2
• Matrix Mittag–Leffler functions have applications in fractional differential equations, multiterm equations, and control, motivating accurate evaluation.The authors aim for results close to machine precision.Source: §1 · p. 2
• Multiple or highly clustered eigenvalues make Schur–Parlett evaluation depend on high-order derivatives of the underlying scalar function.The required derivative order depends on the matrix spectrum.Source: §1 · p. 2
• Three derivative-evaluation methods work in limited complex-plane regions and parameter ranges, so the authors combine them to handle matrices with any eigenvalue location.The topic was largely unexplored beyond first derivatives.Source: §1 · p. 2
• The paper also analyzes matrix Mittag–Leffler conditioning to assess sensitivity to data perturbations, which are inevitable in numerical applications.Source: §1 · pp. 2–3

2 Applications of matrix Mittag-Leffler functions

• Matrix Mittag–Leffler evaluation is useful across applications, including fractional differential equations and control.The section offers a sampling of fields rather than a thorough review.Source: §2 · p. 3

2.1 Solution of systems of FDEs

• For linear systems of fractional differential equations, the exact solution can be evaluated directly at any time once matrix Mittag–Leffler computation is available.Step-by-step methods can incur considerable costs because fractional operators have persistent memory.Source: §2.1 · p. 3

2.2 Solution of time-fractional partial differential equations

• For time-fractional PDEs, spatial discretization yields a semi-linear fractional system; the matrix-function formulation avoids product-integration convergence limitations and addresses stiffness in the linear part.High-order methods are possible under the reasonable assumption of a sufficiently smooth source term.Source: §2.2 · p. 4

2.3 Solution of linear multiterm FDEs

• Linear multiterm fractional differential equations can be represented using matrix Mittag–Leffler functions, while an alternative uses scalar multivariate Mittag–Leffler functions.Comparing the approaches and their computational trade-offs is beyond this work's scope.Source: §2.3 · pp. 4–5

2.4 Controllability and observability of fractional linear systems

• Fractional-system controllability and observability can be characterized using positive definiteness of their respective Gramian matrices.Computing the Gramians can support control, state reconstruction, and transfer-function models derived by linearization.Source: §2.4 · p. 5

3.1 Theoretical background

• The matrix Mittag-Leffler series is useful for defining the function, but usually not for computation because it inherits the scalar series issues and adds matrix-computation cost.Source: §3.1 · p. 6
• The Jordan-form definition requires the values of f and its derivatives at each eigenvalue up to one less than that eigenvalue’s index.The index is the smallest k for which (A − λ_i I)^k = 0.Source: §3.1 · p. 6
• When multiple eigenvalues occur, defining a matrix function requires scalar-function derivatives, and no numerical method for arbitrary-order Mittag-Leffler derivatives had been considered beyond first order.Source: §3.1 · p. 6
• Jordan canonical form is rarely used in practice because its similarity transformation may be ill-conditioned and cannot be reliably computed in floating-point arithmetic.Source: §3.1 · p. 6

3.2 Schur-Parlett algorithm

• The Schur-Parlett algorithm uses a backward-stable Schur decomposition, reorders and blocks the triangular factor, and evaluates the matrix function through block computations.The diagonal blocks are evaluated by Taylor series about the mean eigenvalue, then the remaining blocks are computed with the block Parlett recurrence.Source: §3.2 · p. 7

4 Computation of derivatives of the ML function

• The Schur-Parlett algorithm requires scalar-function derivatives, so the section introduces three approaches and investigates their differing numerical behavior to combine them in an efficient algorithm.The methods are intended to compute derivatives for arguments throughout the complex plane and to orders determined by the spectrum of the matrix.Source: §4 · pp. 7–8

4.1 Series expansion

• Series truncation is constrained by Gamma overflow above 171.624 and can suffer severe cancellation when |z| > 1 and |arg(z)| > απ/2.Round-off estimates are only available a posteriori; the authors recommend sorting terms by modulus and using the mean of two error bounds to assess reliability.Source: §4.1 · pp. 8–9

4.2 Numerical inversion of the Laplace transform

• Laplace-transform inversion balances discretization and truncation errors, but contour placement must limit round-off and may require subtracting residues of poles lying to the contour’s right.For derivatives, the origin’s order-(k + 1) singularity must also enter the error balance; high-order derivatives can lose accuracy near the negative real axis or when poles constrain contour selection.Source: §4.2 · pp. 11, 13

4.3 Summation formulas

• The Prabhakar-type summation formula avoids the Djrbashian formula’s cancellation near z = 0 by avoiding division by z^k, making it more reliable for small |z|.The Djrbashian formula represents derivatives as linear combinations of Mittag-Leffler function values, but differences of nearly equal values become unreliable near the origin, especially at high orders.Source: §4.3 · pp. 14–15
• For small |z|, the Prabhakar summation formula is more accurate than the Djrbashian formula; at larger |z|, their performance is nearly similar.The comparisons use reference derivatives computed from the series expansion in 2000-digit arithmetic.Source: §4.3 · p. 15, Figure 1 · p. 16, Figure 2 · p. 16

4.4 Combining algorithms with derivatives balancing

• The Prabhakar formula gives accurate, satisfactory results in most tested cases, but higher derivative orders can increase errors through accumulated round-off.Numerical Laplace-transform inversion can also lose accuracy at high orders when the integration process encounters the strong singularity in Hk(s; z).Source: §4.4 · pp. 15–16
• For critical cases, the authors recommend mixed strategies, such as using the truncated series for small arguments and switching methods as |z| grows.Error estimates can support an automatic switching procedure.Source: §4.4 · p. 16
• Derivative balancing expresses kth-order ML derivatives through lower-order derivatives, helping avoid high-order Laplace-transform inversion and reduce terms in the Prabhakar formula.The combined approach uses the Prabhakar summation formula with derivatives evaluated by numerical Laplace-transform inversion.Source: §4.4 · pp. 16–17
• For a fifth derivative, balancing with p = 1 kept errors moderate across the interval and removed the inversion method’s peak on [0, 2].The summation formula lost accuracy at higher order, while numerical Laplace-transform inversion could produce a peak when its contour approached the singularity.Source: §4.4 · p. 17, Figure 3 · p. 17
• Across experiments on a wide set of arguments and parameters, derivative errors remained in the 10^-13–10^-15 range, which the authors considered satisfactory for applications.The plotted error was |E − Ẽ|/(1 + |E|), with E a reference value and Ẽ the evaluated derivative. Lower-order derivatives were sometimes more accurate, but no general trend could be inferred.Source: §4.4 · p. 17

5 Conditioning of the matrix ML function

• The Fréchet derivative provides a tool for sensitivity analysis, and the absolute and relative condition numbers are expressed through it.The section motivates conditioning by noting that rounding errors can affect computation even with exact data.Source: §5 · pp. 18–19
• For 25 Matlab-gallery matrices of dimension 10 × 10, conditioning is nearly independent of α except for chebspec; Pascal matrices are severely ill-conditioned.For the n = 10 Pascal matrix, κ(A) ≈ 8.13 × 10^9; chebspec conditioning is huge for small α.Source: §5 · p. 19
• For most test matrices, the matrix ML function is better conditioned than its argument; in the remaining cases, conditioning is not amplified much.The comparison uses κr(Eα,1, A)/κ(A).Source: §5 · p. 19
• The estimate κr(Eα,1, A)ϵ reasonably bounds computation error, suggesting the matrix ML evaluation is normwise forward stable.The section reports this behavior for the tested computations.Source: §5 · p. 19

6 Numerical experiments

• Experiments use Matlab's funm for Schur-Parlett and the combined derivative algorithm, with Maple 15 at 2000-digit precision for reference values.Errors are computed using the Frobenius norm as ∥E − Ẽ∥/(1 + ∥E∥).Source: §6 · pp. 19–20
• The conditioning plots cover 25 Matlab-gallery matrices of dimension 10 × 10, showing both κr(Eα,1, A) and its ratio to κ(A).Source: Figure 6 · p. 20, Figure 7 · p. 20

Example 1: the Redheffer matrix

• Errors for Redheffer matrices with n from 4 to 20 are close to machine precision, despite requiring derivatives up to order 24.The experiment uses β = 1; the matrices have highly clustered eigenvalues.Source: Example 1: the Redheffer matrix · pp. 20–21, Figure 8 · p. 21

Example 2: matrices with selected eigenvalues

• For four 40 × 40 matrices with clustered eigenvalues, evaluation is reasonably accurate, and κr(Eα,1, A)ϵ reasonably estimates the error.The maximum required derivative order is listed in Table 2; Figure 9 plots errors and conditioning bounds.Source: Example 2: matrices with selected… · p. 21, Table 2 · p. 21, Figure 9 · p. 22

Example 3: solution of a multiterm FDE

• For the multiterm FDE, matrix-ML evaluation gives excellent solution accuracy directly at any t, while trapezoidal PI needs very small steps and substantial computational cost.The equation is reformulated using a 16 × 16 companion matrix with simple and double eigenvalues.Source: Example 3: solution of a multiterm FDE · p. 22

• The paper discusses matrix Mittag–Leffler evaluation and illustrates applications of the function.Source: §7 · p. 22
• Figure 11 presents errors in the solution of multiterm FDE (34) for α = 0.8.Source: Figure 11 · p. 23
• Numerical experiments show matrix Mittag–Leffler functions can be evaluated with great accuracy even for not simple matrix arguments.Source: §7 · p. 23
• Different methods for evaluating Mittag–Leffler derivatives are combined to target high accuracy for matrices with any kind of spectrum.The paper notes that accurate evaluation of these derivatives had not previously been addressed beyond first derivatives.Source: §7 · p. 23
• The conditioning analysis indicates that errors can be kept under control when evaluating matrix Mittag–Leffler functions.Source: §7 · p. 23