Electronic Journal of Differential Equations, Vol. 2021 (2021), No. 97, pp. 1–18. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu EXTENDING PUTZER’S REPRESENTATION TO ALL ANALYTIC MATRIX FUNCTIONS VIA OMEGA MATRIX CALCULUS ANTÔNIO FRANCISCO NETO Abstract. We show that Putzer’s method to calculate the matrix exponential in [28] can be generalized to compute an arbitrary matrix function defined by a convergent power series. The main technical tool for adapting Putzer’s formulation to the general setting is the omega matrix calculus; that is, an extension of MacMahon’s partition analysis to the realm of matrix calculus and the method in [8]. Several results in the literature are shown to be special cases of our general formalism, including the computation of the fractional matrix exponentials introduced by Rodrigo [30]. Our formulation is a much more general, direct, and conceptually simple method for computing analytic matrix functions. In our approach the recursive system of equations the base for Putzer’s method is explicitly solved, and all we need to determine is the analytic matrix functions. 1. Introduction Let φ(t) ∈ CN with N <∞, then the solution of the initial value problem φ′(t) = Aφ(t), φ(0) = φ0 (1.1) is φ(t) = exp(tA)φ0, where the matrix exponential exp(A) is defined by exp(A) = ∑ k≥0 Ak/k!, (1.2) which is a matrix valued convergent power series for any A ∈ CN×N . In general terms Putzer [28] constructed a representation of the matrix exponential in (1.2) avoiding the use of the Jordan canonical form and requiring [28, Theorem 2] or not [28, Theorem 1] the knowledge of the eigenvalues of A. Putzer’s method has the nice feature of being generic; that is, it holds for any square matrix even with repeated eigenvalues. Other papers searching for analogues of Putzer’s result also appeared in different contexts. We recall [1] where the role of the matrix exponential in (1.2) is replaced by the matrix logarithm and extensions to the discrete setting [9, 18] with the matrix power playing the role of the matrix exponential. 2010 Mathematics Subject Classification. 15A16, 26A33. Key words and phrases. Putzer’s method; omega matrix calculus; matrix valued convergent series; Mittag-Leffler function; fractional calculus. ©2021. This work is licensed under a CC BY 4.0 license. Submitted March 3, 2021. Published December 7, 2021. 1 2 A. F. NETO EJDE-2021/97 Recently, the fractional analogues of the IVP in (1.1) and their associated solu- tions were considered in an interesting article [30]. For a historical account on the origins of fractional calculus we refer the reader to the comprehensive work [24, 31] and for applications to [10, 30, 33, 35, 34] and references therein. We review some of the central results of [30] for clearness. In [30] two distinct and well-known frac- tional versions of the usual derivative were considered; that is, the Caputo fractional derivative C 0 D α t f(t) = 0D −(dαe−α) t Ddαef(t), t > 0 (1.3) and the Riemann-Liouville fractional derivative 0D α t f(t) = Ddαe0D −(dαe−α) t f(t), t > 0 (1.4) with α > 0 and dαe the least integer greater than or equal to α (0 ≤ dαe − α < 1). We remark that Ddαe is the ordinary differential operator of order dαe and the Riemann-Liouville fractional integral of order α is given by 0D −α t f(t) = 1 Γ(α) ∫ t 0 (t− s)α−1f(s)ds, t > 0, where Γ(α) is the Euler’s gamma function [26, Chapter 1]. Note that (1.3) and (1.4) agree with [26, (2.172) and (2.171)] upon setting p → α and n → dαe. We follow the notation in [30]: D−α, Dα, and Dα ∗ stand for 0D −α t , 0D α t , and C 0 D α t , respectively. In this way the solution of the IVP Dα ∗Φ(t) = AαΦ(t), DkΦ(0+) = Ak, k ∈ {0} ∪ [dαe − 1] (1.5) is given by the Caputo fractional exponential Exp∗(tA;α) = dαe−1∑ j=0 (tA)jEα,j+1((tA)α) (1.6) and the solution of the IVP DαΦ(t) = AαΦ(t), Dk−dαe+αΦ(0+) = Ak, k ∈ {0} ∪ [dαe − 1] (1.7) is given by the Riemann-Liouville fractional exponential Exp(tA;α) = tα−dαe dαe−1∑ j=0 (tA)jEα,α−dαe+j+1((tA)α) (1.8) with [n] = {1, . . . , n} (n a positive integer) and Eα,β(tA) = ∑ k≥0 (tA)k Γ(αk + β) (1.9) an entire function if α, β > 0 known as the matrix Mittag-Leffler function [26, Chapter 1]. We remark that with our choices Eα,j+1 in (1.6) and Eα,α−dαe+j+1 in (1.8) are entire (α − dαe + j + 1 > 0, ∀j ∈ {0} ∪ [dαe − 1]) and [30, Lemma 2.1] shows how to define Aα for any α > 0. More precisely, we have Cmk×mk 3 Aα k =  aαk ( α 1 ) aα−1 k · · · ( α mk−2 ) aα−mk+2 k ( α mk−1 ) aα−mk+1 k 0 aαk · · · ( α mk−3 ) aα−mk+3 k ( α mk−2 ) aα−mk+2 k ... ... . . . ... ... 0 0 · · · aαk ( α 1 ) aα−1 k 0 0 · · · 0 aαk  EJDE-2021/97 MATRIX FUNCTIONS VIA OMEGA MATRIX CALCULUS 3 with Aα = M(⊕rk=1A α k )M−1 and m1 + · · ·+mr = N . Here M is a nonsigular matrix such that M−1AM = ⊕rk=1Ak ≡ J (1.10) and ( α β ) = Γ(α+ 1) Γ(β + 1)Γ(α− β + 1) with α, β ∈ C. From now on, O and I stand for the null and the identity matrices, respectively. Note that Exp∗(tA; 1) = exp(tA) = Exp(tA; 1), Exp∗(O;α) = I = Exp(O;α), which follows by observing that E1,1(tA) = ∑ k≥0 (tA)k Γ(k + 1)︸ ︷︷ ︸ =k! (1.2) = exp(tA). Therefore, if we set α = 1 in (1.5) and (1.7) we recover the well-known IVP in (1.1) with (1.6) and (1.8) reducing to the corresponding solution given by (1.2). An adaptation of Putzer’s method to compute the fractional exponential functions in (1.6) and (1.8) as finite linear combinations of constant matrices with time-varying coefficients was obtained in [30, Theorems 5.1 and 5.2]. We also recall a comment taken from [30] highlighting the need for computational methods to determine the fractional matrix exponentials as close to the ordinary matrix exponential in (1.2) as possible and quoted verbatim here: “The numerical computation of these fractional matrix exponentials, akin to [18] for the usual matrix exponential, is of independent interest.” Note that [30, Ref. [18]] stands for [25] here. See also [14]. Even the extension of the basic properties of the ordinary exponential in (1.2) to the fractional setting is a subject of considerable interest [27, 32]. In this respect, [30] asks if the semigroup property of the matrix exponential in (1.2) is valid in the fractional setting with reference to (1.6) and (1.8). Summarizing, the extension of basic results from ODEs to the context of fractional calculus and the construction of computational procedures to determine the fractional matrix exponentials as close to the usual setting as possible are of general interest. The aim of this work is to introduce a general, direct, and conceptually simple method to compute analytic matrix functions, including the Mittag-Leffler function in (1.9) and the fractional matrix exponentials of [30] in (1.6) and (1.8). In our approach, there is no need to adapt Putzer’s formulation as in [1, Theorem 3] or [30, Theorems 5.1 and 5.2] dealing with the matrix logarithm and fractional matrix exponentials, respectively. More precisely, we obtain at once the solution of the recursive systems of equations in [28, Theorem 2], [9, Theorem 1], [1, Theorem 3], and [30, Theorems 5.1 and 5.2]. We also show that the determination of the analytic matrix functions depends on the same recursive system of equations in [28] which we explicitly solve. Furthermore, as our method relies on the usual matrix exponential it is more amenable to be treated by standard approaches available to compute (1.2). The main technical tool for our method is based on an extension of the usual Omega Calculus (i.e. MacMahon’s partition analysis [20]) to the context of Matrix Analysis introduced recently, the Omega Matrix Calculus (OMC for short) [11, 12], and an approach to compute the matrix exponential using the Jordan canonical 4 A. F. NETO EJDE-2021/97 form and properties of the minimal polynomial of a matrix [8]. We remark that OMC is a useful tool in representing a function defined by a convergent power series in terms of other functions under the action of the Omega operator. This feature comprises the starting point of [11] where the OMC was introduced and the inverse of a certain matrix function was used to obtain properties of the exponential in (1.2) (see [11, Lemma 2.3]). Therefore, motivated by the need to compute the fractional matrix exponentials in (1.6) and (1.8) akin to the exponential in (1.2) and with the aforementioned useful feature of OMC in mind, it is natural to explore OMC in the context of representing a matrix function defined by a convergent power series in terms of the exponential in (1.2) as we do here. This work is organized as follows. In Section 2 we state our main result Theorem 2.3. In Section 3 we give some auxiliary results to be used in Section 4 devoted to the proof of Theorem 2.3. In Section 5 we establish contact with previous results in the literature and we show the versatility of our main result. More precisely, Theorem 2.3 implies [28, Theorem 2], [9, Theorem 1], [1, Theorem 3], and [30, Theorems 5.1 and 5.2]. An example is also included in Section 5 in order to illustrate the simplicity of the proposed method. Finally, we summarize our findings in the conclusion. 2. Statement of main results First we introduce some notation and give a definition. We let fm = F (m)(0) and F (tA) = ∑ m≥0 fmt mAm/m! (2.1) be a convergent matrix valued function. We remark that there are other ways to define matrix functions and the connection between the several definitions is discussed in [29]. See also [16]. Of course, if we set F = exp in (2.1) we recover (1.2). Several other analytic matrix functions such as those considered in [1, 21, 22] and (1.9) are all special cases of (2.1) with appropriate domains of convergence. In this way, our access to OMC is based on the definition that follows. Throughout this article, 0n stands for the null vector in Cn. Definition 2.1. Let Xa ∈ CN×N for each a ∈ Zn and λa = λa11 · · ·λann . We define the linear operator acting on absolutely convergent matrix valued expansions in (2.1) by λ Ω = ∞∑ a1=−∞ · · · ∞∑ an=−∞ Xaλ a def = X0n in an open neighborhood of the complex circles |λi| = 1. Remark 2.2. Definition 2.1 is well-posed in the sense that we can ensure that all expressions considered here have no singularities in the λi variable in an open neighborhood of the circle |λi| = 1. As remarked in [2], this is an important ingre- dient leading to unique Laurent expansions (otherwise, ambiguous results appear as discussed in the introduction of [2]). In what follows, d stands for the degree of the minimal polynomial of A [4, 13, 17, 19]. We write the set of eigenvalues of A as S = {ai}ri=1 (2.2) EJDE-2021/97 MATRIX FUNCTIONS VIA OMEGA MATRIX CALCULUS 5 and the multiset S = {αi}di=1 (2.3) with α1, . . . , αn1 = a1, αn1+1, . . . , αn1+n2 = a2, and so on, until we obtain αn1+···+nr−1+1, . . . , αn1+···+nr−1+nr = ar. In other words, the elements of S are the distinct eigenvalues of A and the elements of S are the eigenvalues of A counted with multiplicity. From now on, F (tA;α) stands for (1.6), (1.8), and (2.1) with F (tA; 1) ≡ F (tA) for F in (2.1). We adopt the convention that ∑j k=i(· · · ) ≡ 0 and ∏j k=i(· · · ) ≡ 1 if i > j and write Cmn = m!/((m− n)!n!) throughout this article. Using Definition 2.1 we can now state our main result. Theorem 2.3. Let α > 0 and Pk(A) = { I if k = 0∏k j=1(A− αjI) if 1 ≤ k ≤ d− 1 . with αj ∈ S in (2.3) and j ∈ [d] = ∪rk=1{i + n1 + · · · + nk−1 + 1}nk−1 i=0 . Then we have F (tA;α) = dαe−1∑ j=0 d−1∑ k=0 yj+1,k+1(t)AjPk(Aα), (2.4) where yj+1,k+1(t) = lim l→∞ l∑ m=0 gjm(t;α) λ Ω = λmxk+1(tα/λ) (2.5) with gjm(t;α) given by gjm(t;α) =  m!tj Γ(αm+j+1) for F in (1.6), m!tj+α−dαe Γ(αm+α−dαe+j+1) for F in (1.8), g1m(t; 1) = fm for F in (2.1) and xj is determined recursively by xi+n1+···+nk−1+1(t) = ∑ i1+···+ik=i k−1∏ m=1 Cim+nm−1 nm−1 (−1)nm aim+nm m|k tik ik! exp(akt) − k−1∑ j=1 nj−1∑ l=0 alk|j ∑ ij+···+ik−1=i C ij+nj−l−1 nj−l−1 k−1∏ m=j+1 Cim+nm−1 nm−1 × k−1∏ m=j (−1)nm aim+nm m|k xl+n1+···+nj−1+1(t) (2.6) with ak ∈ S in (2.2) and ak|j = ak − aj. Before we prove Theorem 2.3 we introduce some auxiliary results. 6 A. F. NETO EJDE-2021/97 3. Auxiliary results We begin this section by recalling some basic results regarding OMC following [11]. See also [2, 20] for the introduction of the Omega package, a computer algebra package in MATHEMATICA that implements the Omega calculus, and for the original formulation in the scalar case (N = 1), respectively. A particular example of a matrix valued function defined by a convergent power series for any A ∈ CN×N and satisfying the requirements of Definition 2.1 is the matrix exponential in (1.2). It follows from [11, Lemma 2.3] that φ(t) = exp(tA)φ0 = λ Ω = exp(tλ)(I −A/λ)−1φ0. (3.1) Note that (3.1) is a generalization of the elimination rule proposed in the last equation of [3] dealing with the scalar case. We need a convergent Neumann series for (I −A/λ)−1 if one wants to use Definition 2.1 to prove (3.1). Following [11, Lemma 2.3], we introduce a rescaling λ→ λ/z with a small complex parameter z 6= 0 such that the Neumann series (I − zA/λ)−1 = ∑ n≥0 (zA/λ)n (3.2) converges if ‖A/λ‖ < 1/|z|. Since all matrix norms are equivalent if N < ∞, we write generically ‖ · ‖ without further specification of the norm used. We analyze now how the parameter z is used to prove (3.1). We have λ Ω = exp(tλ/z)(I − zA/λ)−1 (1.2),(3.2) = ∑ m,n≥0 tmzn m!zm An λ Ω = λm−n︸ ︷︷ ︸ =δm,n (1.2) = exp(tA) with δm,n the Kronecker delta. After the application of the Omega operator the parameter z cancels out in (3.1)! For this reason, we omit z from the notation in the right hand side of (3.1), but we assume this observation is used throughout this work whenever necessary to ensure convergence. We recall that [28, Theorems 1 and 2] gives a representation of the exponential function in (1.2) as a finite matrix sum. Other relevant papers in this direction employing distinct methods and including other matrix functions comprise [8] using the Jordan canonical form and properties of the minimal polynomial of a matrix, [23, 36] using the Horner polynomials, [5, 6, 7, 21, 22] concerning a combinatorial method based on generalized Fibonacci sequences, and [15] using path-sums. The representation of the analytic matrix function (2.1) as a finite sum is expected from the Cayley-Hamilton theorem which relates AN ∈ CN×N to lower powers of A. Note that the proof of [28, Theorems 1 and 2] uses the Cayley-Hamilton theorem, but the proof holds if the characteristic polynomial p(x) = det(xI −A) = r∏ i=1 (x− ai)mi (3.3) is replaced by any annihilating polynomial; that is, a polynomial f such that f(A) = O. In particular, it holds for the annihilating polynomial of the smallest possible EJDE-2021/97 MATRIX FUNCTIONS VIA OMEGA MATRIX CALCULUS 7 degree d called the minimal polynomial [13, 17, 19] q(x) = r∏ i=1 (x− ai)ni = xd + cd−1x d−1 + · · ·+ c1x+ c0 (3.4) with 1 ≤ ni ≤ mi. Therefore, it is clear that the expansion F (tA) = d−1∑ k=0 fk(t)Pk(A) (3.5) for F in (2.1) holds and the problem now becomes how to determine the coefficients fk(t) in (3.5). We write J = (⊕ri=1Ai)⊕B (3.6) for J in (1.10) with Cni×ni 3 Ai =  ai 1 · · · 0 0 0 ai · · · 0 0 ... ... . . . ... ... 0 0 · · · ai 1 0 0 · · · 0 ai  representing the block of J associated with the eigenvalue ai of A with i ∈ [r]. We also let B in (3.6) represent the remaining blocks, which might be of order zero if the minimal and the characteristic polynomials coincide as is the case if mi = 1 in (3.3) with i ∈ [r]. See [8, Theorem 2]. The Jordan canonical form in (3.6) is used in [8] as a key strategy to obtain the system of equations satisfied by the time-varying coefficients; that is, the coefficients fk(t) in (3.5) with F ≡ exp. As discussed in [8], this is not a restrictive procedure, because similar matrices have the same explicit expansion in (3.5). More precisely, we have the following lemma. Lemma 3.1. The following equivalence holds F (tA) = d−1∑ k=0 fk(t)Pk(A)⇔ F (tJ) = d−1∑ k=0 fk(t)Pk(J). The proof of the above lemma follows directly from (1.10). We use Lemma 3.1 in the proof of Theorem 2.3 in the next section. We also need two auxiliary lemmas. From now on, we write A = ((A)i,j); that is, the element in row i and column j of the matrix A is (A)i,j . Lemma 3.2. Let A =  a1 1 · · · 0 0 0 a2 · · · 0 0 ... ... . . . ... ... 0 0 · · · ad−1 1 0 0 · · · 0 ad  with ai 6= 0 for i ∈ [d]. Then (A−1)i,j = (−1)i−j∏j k=i ak with j ≥ i. 8 A. F. NETO EJDE-2021/97 Proof. The proof follows by induction on d. We need a formula for the inverse of a block triangular matrix A−1 = ( B1 B2 O B3 )−1 = ( B−1 1 −B−1 1 B2B −1 3 O B−1 3 ) (3.7) for non-singular blocks Bi with i = 1, 3. The base case is direct and gives B−1 3 = a−1 d . The induction step comprises assuming (B−1 1 )i,j = (−1)i−j∏j k=i ak . (3.8) Therefore, −(B−1 1 B2B −1 3 )i = − (B−1 1 )i,d−1 ad = (−1)i−d∏d k=i ak using (3.8) to obtain the second equality and the result follows from (3.7). � Lemma 3.3. The coefficients fk(t) in (3.5) are unique. Proof. The proof by contradiction is similar to the one given in [8, Proposition 2]. Recall that Pk(A) in Theorem 2.3 is a polynomial of degree at most d − 1. Next, observe that d−1∑ k=0 fk+1(t)Pk(A) = F (tA) = d−1∑ k=0 gk+1(t)Pk(A) ⇒ d−1∑ k=0 (fk+1(t)− gk+1(t))Pk(A) = O. Hence we have fk+1(t) ≡ gk+1(t) for all k ∈ [d − 1], for otherwise we would have an annihilator with smaller degree than q in (3.4), contradicting the fact that q is minimal. � In the next section we prove Theorem 2.3 using Lemmas 3.1 and 3.2. 4. Proof of Theorem 2.3 For clearness, since many blocks are involved in the decomposition (3.6), we use the notation Ik ≡ Ink ∈ Cnk×nk for the identity matrix associated with the block indexed by k in J . Similar considerations apply to the null matrix Ok. For simplicity of notation, we write Ai|j = Ai − ajIi. We are now ready to prove Theorem 2.3. Proof of Theorem 2.3. First, we show that exp(tA) = d−1∑ k=0 xk+1(t)Pk(A) (4.1) with xk+1(t) determined recursively by (2.6). The finite sum in (4.1) holds since the minimal polynomial in (3.4) is (by definition) an annihilator. Following Lemma 3.1, we restrict our attention to J in (3.6) to determine the coefficients xk+1(t) in (4.1). We have( ⊕rk=1 exp(tAk) ) ⊕ exp(tB) = exp(tJ) = d−1∑ l=0 xl+1(t) (( ⊕rk=1 Al k ) ⊕Bl ) . EJDE-2021/97 MATRIX FUNCTIONS VIA OMEGA MATRIX CALCULUS 9 Without loss of generality we consider the blocks Ak with k ∈ [r] except B, which results in redundant information as already noted in [8, Lemma 4]. Therefore, from now on, we focus on exp(tAk) = d−1∑ l=0 xl+1(t)Pl(Ak) with k ∈ [r]. Next, using (3.1) we obtain (exp(tAk))i,j = λ Ω = exp(tλ)((Ik −Ak/λ)−1)i,j . Using Lemma 3.2 we have( (Ik −Ak/λ)−1 ) i,j = −λ ( (Ak − λIk)−1 ) i,j = (−1)i−j−1λ (ak − λ)j−i+1 = 1 λj−i(1− ak/λ)j−i+1 . (4.2) Therefore, it follows that (exp(tAk))i,j = λ Ω = exp(λt) λj−i(1− ak/λ)j−i+1 = D (j−i) ak (j − i)! λ Ω = exp(λt) 1− ak/λ︸ ︷︷ ︸ =exp(akt) = tj−i exp(akt) (j − i)! with 1 ≤ i ≤ j ≤ nk. Next, we write [d] = ∪ri=1{l + n1 + · · · + ni−1 + 1}ni−1 l=0 to obtain ( exp(tAk) ) i,j = n1−1∑ l=0 xl+1(t) ( Al k|1︸︷︷︸ =Pl(Ak) ) i,j + n2−1∑ l=0 xl+n1+1(t) ( An1 k|1A l k|2︸ ︷︷ ︸ =Pl+n1 (Ak) ) i,j + · · ·+ nk−1∑ l=0 xl+n1+···+nk−1+1(t) ( k−1∏ m=1 Anm k|mAl k|k︸ ︷︷ ︸ =Pl+n1+···+nk−1 (Ak) ) i,j (4.3) using Ank k|k = Ok. By multiplying ∏k−1 m=1 A−nmk|m on both sides of (4.3) we find ( k−1∏ m=1 A−nmk|m exp(tAk) ) i,j = n1−1∑ l=0 xl+1(t) ( k−1∏ m=2 A−nmk|m A −(n1−l) k|1 ) i,j + n2−1∑ l=0 xl+n1+1(t)( k−1∏ m=3 A−nmk|m A −(n2−l) k|2 )i,j + · · ·+ nk−1∑ l=0 xl+n1+···+nk−1+1(t) ( Al k|k ) i,j . Equation (2.6) now follows by noting that( Al k|k ) i,j = λ Ω = λl((Ik −Ak/λ)−1)i,j ∣∣∣ ak=0 (4.2) = λ Ω = λl−j+i = δj,l+i 10 A. F. NETO EJDE-2021/97 and (A−mk|l )i,j = nk∑ i2,...,im=1 m∏ n=1 (A−1 k|l )in,in+1 = (−1)i−ja−mk|l ∑ j1+···+jm=j−i a−j1−···−jmk|l = Cj−i+m−1 m−1 (−1)m aj−i+ml|k , where i1 ≡ i and im+1 ≡ j. The change of variables jn = in+1− in ≥ 0 (recall that the matrix A−1 k|l is upper triangular) was used to obtain the second equality. The last equality follows from the well-known fact that the binomial coefficient Cj−i+m−1 m−1 counts the number of solutions of the linear diophantine equation j1+· · ·+jm = j−i in non-negative integers. Next, we show that F (tA) = lim l→∞ l∑ m=0 fm λ Ω = λm exp(tA/λ). (4.4) Indeed, starting with the right-hand side of (4.4) we obtain lim l→∞ l∑ m=0 fm λ Ω = λm exp(tA/λ) (1.2) = lim l→∞ l∑ m=0 λ Ω = (fmλ m ∑ k≥0 tkAk k!λk ) = lim l→∞ l∑ m=0 ∑ k≥0 fm tkAk k! λ Ω = λm−k︸ ︷︷ ︸ =δk,m (2.1) = F (tA). Equation (2.4) for F in (2.1) follows using (4.1) with the replacement t → t/λ in the right hand side of (4.4). Finally, the proof of (2.4) for F in (1.6) and (1.8) is similar, except that we use exp(tαAα/λ) instead of exp(tA/λ), to obtain tβ−1Eα,β((tA)α) = lim l→∞ l∑ m=0 gjm(t;α) λ Ω = λm exp(tαAα/λ), (4.5) where β is given by j + 1 and α−dαe+ j + 1 for F in (1.6) and (1.8), respectively. The required result now follows using (4.1) along with the replacement (t,A) → (tα/λ,Aα) in the right-hand side of (4.5) and multiplying by Aj before summing over j ∈ {0} ∪ [dαe − 1]. � From now on, we write y1,k+1(t) ≡ yk+1(t) in (2.4). Note that, as a consequence of Theorem 2.3, the computation of an analytic matrix function is reduced to the calculation of xk+1(t) using (2.6). Another nice feature of Theorem 2.3 is presented in the next section. More precisely, we show that Theorem 2.3 along with Lemma 3.3 imply several known results as special cases. EJDE-2021/97 MATRIX FUNCTIONS VIA OMEGA MATRIX CALCULUS 11 5. Connection with previous work In this section we show that Theorem 2.3 and Lemma 3.3 imply [28, Theorem 2], [9, Theorem 1], [1, Theorem 3], and [30, Theorems 5.1 and 5.2]. 5.1. Theorem 2.3, Lemma 3.3, and [28, Theorem 2]. Putzer’s original formu- lation is obtained from Theorem 2.3 by setting fm ≡ 1 in (2.4) and observing that yk+1(t) (2.5) = lim l→∞ l∑ m=0 λ Ω = λmxk+1(t/λ) = xk+1(t). We show that our approach implies [28, Theorem 2]. More precisely, we have x′1(t) = α1x1(t), x′k(t) = αkxk(t) + xk−1(t) if 2 ≤ k ≤ d (5.1) with x1(0) = 1 and xk(0) = 0 for 2 ≤ k ≤ d. (5.2) By taking the t-derivative on both sides of (4.1) we find that A exp(tA) = d−1∑ k=0 xk+1(t)APk(A) = d−1∑ k=0 x′k+1(t)Pk(A). (5.3) Next, observe that APk(A) = Pk+1(A) + αk+1Pk(A). (5.4) Therefore, d−1∑ k=0 xk+1(t)APk(A) = d−1∑ k=0 αk+1xk+1(t)Pk(A) + d−2∑ k=0 xk+1(t)Pk+1(A), (5.5) which follows from Pd(A) = q(A) = O using (3.4). It follows from (5.3) and (5.5) that d−1∑ k=0 x′k+1(t)Pk(A) = d−1∑ k=0 αk+1xk+1(t)Pk(A) + d−1∑ k=1 xk(t)Pk(A). (5.6) Thus, by comparing the corresponding coefficients of Pk(A) in (5.6) we obtain the desired result using Lemma 3.3. Finally, setting t = 0 in (4.1) we obtain exp(tA)|t=0 = I = d−1∑ k=0 xk+1(0)Pk(A) and using Lemma 3.3 once again we obtain the initial conditions stated in (5.2). 5.2. Theorem 2.3, Lemma 3.3, and [9, Theorem 1]. We show that the main result of [9, Theorem 1] follows from our approach using (2.1) with F (A) = An. In this case the analogue of the IVP in (1.1) is given by x(n+ 1) = Ax(n), x(0) = x0 ⇒ x(n) = Anx0 with the role of exp(A) in (1.2) played by An in the discrete setting. See [9, Theorems 1 and 2] and [18]. We have An = d−1∑ k=0 xk+1(n)Pk(A), (5.7) 12 A. F. NETO EJDE-2021/97 where x1(n+ 1) = α1x1(n), xk(n+ 1) = αkxk(n) + xk−1(n) if 2 ≤ k ≤ d with x1(0) = 1 and xk(0) = 0 for 2 ≤ k ≤ d. The proof of (5.7) follows from Theorem 2.3 (with t = 1) by observing that exp(A/λ) = d−1∑ k=0 xk+1(1/λ)Pk(A) to obtain An = n! λ Ω = λn exp(A/λ) = d−1∑ k=0 n! (λ Ω = λnxk+1(1/λ) ) Pk(A) and the result follows using (5.1). Indeed, we have n! λ Ω = λnx(1/λ) = n! λ Ω = λn exp(B/λ)x0 = Bnx0 = x(n), where x(1/λ) = (x1(1/λ), . . . , xd(1/λ))T , x(n) = (x1(n), . . . , xd(n))T , B =  α1 0 · · · 0 0 1 α2 · · · 0 0 ... ... . . . ... ... 0 0 · · · αd−1 0 0 0 · · · 1 αd  , and x0 = (1,0d−1)T . 5.3. Theorem 2.3, Lemma 3.3, and [1, Theorem 3]. In this case we have F (tA) ≡ ln(I + tA) with ‖tA‖ < 1 and we show that (1 + α1t)y ′ 1(t) = α1, (1 + α2t)y ′ 2(t) = −ty′1(t) + 1, (1 + αk+1t)y ′ k+1(t) = −ty′k(t) if 2 ≤ k ≤ d− 1 (5.8) with y1(0) = · · · = yd(0) = 0. (5.9) Our strategy here is to show that y′l(t) given in Theorem 2.3 with fm = m!(−1)m−1/m satisfying system (5.8). In this case we have ln(I + tA) = d−1∑ k=0 yk+1(t)Pk(A). (5.10) By taking the t-derivative on both sides of (5.10) we have A I + tA = d−1∑ k=0 y′k+1(t)Pk(A)⇒ A = (I + tA) d−1∑ k=0 y′k+1(t)Pk(A). EJDE-2021/97 MATRIX FUNCTIONS VIA OMEGA MATRIX CALCULUS 13 Thus, we find that α1P0(A) + P1(A) = A = d−1∑ k=0 (1 + αk+1t)y ′ k+1(t)Pk(A) + t d−2∑ k=0 y′k+1(t)Pk+1(A) using (5.4) and Pd(A) = O. By comparing the corresponding coefficients on both sides of the equation above we find the desired result using Lemma 3.3. Finally, setting t = 0 in (5.10) we obtain ln(I + tA)|t=0 = O = d−1∑ k=0 yk+1(0)Pk(A) and using Lemma 3.3 once again we obtain the initial conditions stated in (5.9). 5.4. Theorem 2.3, Lemma 3.3, and [30, Theorem 5.1]. We let j, l ∈ {0}∪[dαe− 1] and k ∈ [d− 1]. First, we recall some well-known facts from [24, 26]. We have Γ−1(−α) = 0 (5.11) if α ∈ N0. See, e.g., [26, (1.6)]. We also recall Dβtα = Γ(α+ 1) Γ(α− β + 1) tα−β (5.12) if α > −1. See [24, Appendix D, Ic] and [26, Chapter 2]. We show that the Caputo fractional matrix exponential admits the representa- tion Exp∗(tA;α) = dαe−1∑ j=0 d−1∑ k=0 yj+1,k+1(t)AjPk(Aα) with Dα ∗ yj+1,1 = α1yj+1,1, Dα ∗ yj+1,k+1 = αk+1yj+1,k+1 + yj+1,k if 1 ≤ k ≤ d− 1, (5.13) Dlyj+1,1(0+) = δj,l, Dlyj+1,k+1(0+) = 0. (5.14) If follows from (4.5) with β = j + 1 that tjEα,j+1(tαAα) = d−1∑ k=0 yj+1,k+1(t)Pk(Aα). (5.15) By taking the fractional derivative Dα ∗ in (1.3) on both sides we find that tjAαEα,j+1(tαAα) = d−1∑ k=0 Dα ∗ yj+1,k+1(t)Pk(Aα) (5.16) using Dα ∗ (tjEα,j+1(tαAα)) (1.3) = D−(dαe−α)Ddαe(tjEα,j+1(tαAα)) (1.9) = D−(dαe−α)( ∑ k≥0 tαk+j−dαe(Aα)k Γ(αk + j − dαe+ 1) ) (5.11) = D−(dαe−α)( ∑ k≥1 tαk+j−dαe(Aα)k Γ(αk + j − dαe+ 1) ) 14 A. F. NETO EJDE-2021/97 (5.12) = ∑ k≥1 tα(k−1)+j(Aα)k Γ(α(k − 1) + j) (1.9) = tjEα,j+1(tαAα). Note that to obtain the third equality we have used j ∈ {0} ∪ [dαe − 1]. We can now follow exactly what we have done in Subsection 5.1 by using (5.4) with the replacement A → Aα, (5.15), and (5.16) to obtain the system in (5.13). We now turn to the initial conditions in (5.14). We show that lim t↘0 Dl(tjEα,j+1(tαAα)) = δj,lI = d−1∑ k=0 Dlyj+1,k+1(0+)Pk(Aα), which implies (5.14). We have Dltαn+j = (αn+ j)(αn+ j − 1) · · · (αn+ j − l + 1)tαn+j−l with αn + j − l > 0 if n > 0. Here, the smallest possible exponent in tαn+j−l for fixed α and j occurs when l = dαe − 1 (recall that j, l ∈ {0} ∪ [dαe − 1]). In this case we have αn+ j − l = α(n− 1) + j︸ ︷︷ ︸ ≥0 +α− dαe+ 1︸ ︷︷ ︸ >0 > 0. Therefore, we obtain limt↘0 t αn+j−l = 0 if n > 0, which implies lim t↘0 Dl(tjEα,j+1(tαAα)) = ( lim t↘0 j(j − 1) · · · (j − l + 1)tj−l/Γ(j + 1) ) I. Next, note that j(j − 1) · · · (j − l + 1)/Γ(j + 1) = (j/Γ(j + 1))︸ ︷︷ ︸ =Γ−1(j) (j − 1) · · · (j − l + 1) = · · · = Γ−1(j − l + 1). Thus, we obtain lim t↘0 Dl(tjEα,j+1(tαAα)) = (lim t↘0 tj−l/Γ(j − l + 1))I. We consider three cases j− l < 0, j− l = 0, and j− l > 0. The first and third cases give zero using (5.11) and taking t↘ 0, respectively. We are left with j = l, which gives lim t↘0 Dl(tjEα,j+1(tαAα)) = I. 5.5. Theorem 2.3, Lemma 3.3, and [30, Theorem 5.2]. We show that the Riemann-Liouville fractional matrix exponential admits the representation Exp(tA;α) = dαe−1∑ j=0 d−1∑ k=0 yj+1,k+1(t)AjPk(Aα) with Dαyj+1,1 = α1yj+1,1, Dαyj+1,k+1 = αk+1yj+1,k+1 + yj+1,k if 1 ≤ k ≤ d− 1, (5.17) Dl−dαe+αyj+1,1(0+) = δl,j , D l−dαe+αyj+1,k+1(0+) = 0. (5.18) EJDE-2021/97 MATRIX FUNCTIONS VIA OMEGA MATRIX CALCULUS 15 If follows from (4.5) with β = α− dαe+ j + 1 that tα−dαe+jEα,α−dαe+j+1(tαAα) = d−1∑ k=0 yj+1,k+1(t)Pk(Aα). (5.19) By taking the fractional derivative Dα in (1.4) on both sides we have tα−dαe+jAαEα,α−dαe+j+1(tαAα) = d−1∑ k=0 Dαyj+1,k+1(t)Pk(Aα), (5.20) where we have used Dα(tα−dαe+jEα,α−dαe+j+1(tαAα)) (5.12) = t−dαe+jEα,−dαe+j+1(tαAα) (5.21) and Eα,−dαe+j+1(tαAα) (1.9) = ∑ k≥0 (tαAα)k Γ(αk − dαe+ j + 1) (5.11) = ∑ k≥1 (tαAα)k Γ(αk − dαe+ j + 1) (1.9) = tαAαEα,α−dαe+j+1(tαAα). Note that (5.21) is a special case of the scalar case in [26, (1.82)]. We can now follow exactly what we have done in Subsection 5.1 by using (5.4) with the replacement A → Aα, (5.19), and (5.20) to obtain the system in (5.17). We now turn to the initial conditions in (5.18). We show that lim t↘0 Dl−dαe+α(tα−dαe+jEα,α−dαe+j+1(tαAα)) = δj,lI = d−1∑ k=0 Dl−dαe+αyj+1,k+1(0+)Pk(Aα), which implies (5.18). Indeed, using Dl−dαe+α(tα−dαe+jEα,α−dαe+j+1(tαAα)) = tj−lEα,j−l+1(tαAα) and following the analogous calculations done in the previous subsection we find the desired result. The details are left to the reader. 5.6. Example. We revisit the example considered in [30, Section 6] to illustrate the simplicity of the approach proposed in Theorem 2.3. First, let us compute Exp∗(tA; 1/2) = E1/2,1(t1/2A1/2) using (1.6) with A = ( 5/2 3/2 3/2 5/2 ) ⇒ A1/2 = ( 3/2 1/2 1/2 3/2 ) . The eigenvalues of A1/2 can be easily seen to be a1 = 1 and a2 = 2, after correcting a minor typo presented in the expression for A1/2 in [30, Section 6]; that is, the off-diagonal elements must be positive. Using (2.6) we have x1(t) = et and x2(t) = e2t − et. (5.22) 16 A. F. NETO EJDE-2021/97 Thus, we obtain Exp∗(tA; 1/2) = y1,1(t)P0(A1/2) + y1,2(t)P1(A1/2) = ( y1,1(t) + y1,2(t)/2 y1,2(t)/2 y1,2(t)/2 y1,1(t) + y1,2(t)/2 ) , (5.23) where y1,1(t) (2.5) = lim l→∞ l∑ m=0 λ Ω = m!λm Γ(m/2 + 1) x1 ( t1/2/λ ) (5.22) = lim l→∞ l∑ m=0 λ Ω = m!λm Γ(m/2 + 1) exp ( t1/2/λ ) = lim l→∞ l∑ m=0 ∑ k≥0 m!tk/2 k!Γ(m/2 + 1) λ Ω = λm−k︸ ︷︷ ︸ =δm,k = E1/2,1 ( t1/2 ) and y1,2(t) (2.5) = lim l→∞ l∑ m=0 λ Ω = m!λm Γ(m/2 + 1) x2 ( t1/2/λ ) (5.22) = lim l→∞ l∑ m=0 λ Ω = m!λm Γ(m/2 + 1) ( exp ( 2t1/2/λ ) − exp ( t1/2/λ )) = E1/2,1 ( 2t1/2 ) − E1/2,1(t1/2) using Theorem 2.3. Note that (5.23) agrees with the example in [30, Section 6] apart from the minor typo previously mentioned and another one in [30, (6.1)]. The correct expression is given in [26, (1.107)] and reads as follows∫ t 0 τβ1−1Eα,β1 (z1τ α)(t− τ)β2−1Eα,β2 (z2(t− τ)α)dτ = z1Eα,β1+β2(z1t α)− z2Eα,β1+β2(z2t α) z1 − z2 tβ1+β2−1 (5.24) with β1, β2 > 0 and z1, z2 arbitrary complex numbers such that z1 6= z2. Note that the factor tβ1+β2−1 is missing in the right hand side of [30, (6.1)]. Using (5.24) and taking into account the corrected value for y1,2(t) we obtain y1,2(t) = 2t1/2E1/2,3/2 ( 2t1/2 ) − t1/2E1/2,3/2(t1/2) = E1/2,1 ( 2t1/2 ) − E1/2,1(t1/2) in agreement with [30, (6.2)]. Finally, we consider Exp(tA; 1/2) = t−1/2E1/2,1/2(t1/2A1/2) using (1.8). Following the procedure used to obtain (5.23) we find that Exp(tA; 1/2) = ( y1,1(t) + y1,2(t)/2 y1,2(t)/2 y1,2(t)/2 y1,1(t) + y1,2(t)/2 ) EJDE-2021/97 MATRIX FUNCTIONS VIA OMEGA MATRIX CALCULUS 17 with y1,1(t) = t−1/2E1/2,1/2(t1/2) and y1,2(t) = 2E1/2,1 ( 2t1/2 ) − E1/2,1(t1/2). Conclusion. We have shown how to generalize Putzer’s method in order to cover all analytic matrix functions using the OMC [11] and the approach described in [8] based on the Jordan canonical form and the minimal polynomial. Several results are shown to be special cases of the general approach introduced in this work. Indeed, we have shown that Theorem 2.3 along with Lemma 3.3 imply [28, Theorem 2], [9, Theorem 1], [1, Theorem 3], and [30, Theorems 5.1 and 5.2]. Furthermore, our ap- proach allows us to recursively compute xk(t) in Theorem 2.3 using (2.6) and, once xk(t) is determined, Putzer’s like representations of all analytic matrix functions follow at once using Theorem 2.3. An example related to the fractional matrix functions in (1.6) and (1.8) is included to illustrate the simplicity and versatility of the method described in Theorem 2.3. This work reinforces the message put forward hitherto in [11, 12] that OMC is a useful tool in generalizing and unifying previous work. Acknowledgements. The author gratefully acknowledges the helpful suggestions of the reviewers and Professor Michael O’Carroll leading to a much better readabil- ity of the article. References [1] C. D. Ahlbrandt, J. Ridenhour; Floquet theory for time scales and Putzer representations of matrix logarithms, J. Difference Equ. Appl., 9 No. 1 (2003), 77–92. [2] G. E. Andrews, P. Paule, A. Riese; MacMahon’s partition analysis: the Omega package, European J. Combin., 22 No. 7 (2001), 887–904. [3] G. E. Andrews, P. Paule, A. Riese; MacMahon’s partition analysis VI: a new reduction algorithm, emphAnn. Comb., 5 No. 3-4 (2001), 251–270. [4] S. Axler; Linear Algebra Done Right, Springer, New York, 2015. [5] R. Ben Taher, M. Mouline, M. Rachidi; Fibonacci-Horner decomposition of the matrix expo- nential and the fundamental system of solutions, Electron. J. Linear Algebra, 15 (2006). [6] R. Ben Taher, M. Rachidi; Linear recurrence relations in the algebra of matrices and appli- cations, Linear Algebra Appl., 330 No. 1-3 (2001), 15–24. [7] R. Ben Taher, M. Rachidi; Linear matrix differential equations of higher-order and applica- tions, Electron. J. Diff. Equ., 2008 No. 95 (2008), 1–12. [8] H.-W. Cheng, S. S.-T. Yau; More explicit formulas for the matrix exponential, Linear Algebra Appl., 262 (1997), 131–163. [9] S. N. Elaydi, W. A. Harris Jr; On the computation of An, SIAM Rev., 40 No. 4 (1998), 965–971. [10] G. Failla, M. Zingales; Advanced materials modelling via fractional calculus: challenges and perspectives, Philos. Trans. Roy. Soc. A, 378 (2020), 20200050. [11] A. Francisco Neto; Matrix analysis and Omega calculus, SIAM Rev., 62 No. 1 (2020), 264– 280. [12] A. Francisco Neto; An approach to isotropic tensor functions and their derivatives via Omega matrix calculus, J. Elasticity, 141 (2020), 165–180. [13] F. R. Gantmacher, The Theory of Matrices, AMS Chelsea Publishing, 1959. [14] R. Garrappa, M. Popolizio; Computing the matrix Mittag-Leffler function with applications to fractional calculus, J. Sci. Comput., 77 No. 1 (2018), 129–153. [15] P.-L. Giscard, S. J. Thwaite, D. Jaksch; Evaluating matrix functions by resummations on graphs: the method of path-sums, SIAM J. Matrix Anal. Appl., 34 No. 2 (2013), 445–469. [16] N. J. Higham; Functions of Matrices: Theory and Computation, SIAM, 2008. [17] R. A. Horn, C. R. Johnson; Matrix Analysis, Cambridge University Press, New York, 2013. [18] M. Kwapisz; Remarks on the calculation of the power of a matrix, J. Difference Equ. Appl., 10 No. 2 (2004), 139–149. 18 A. F. NETO EJDE-2021/97 [19] P. Lancaster, M. Tismenetsky; The Theory of Matrices: with Applications, Academic Press, San Diego, 1985. [20] P. A. MacMahon; Combinatory Analysis, Volumes I and II, 137, AMS Chelsea Publishing, Providence, 2001. [21] J. A. Marrero, R. Ben Taher, Y. El Khatabi, M. Rachidi; On explicit formulas of the principal matrix pth root by polynomial decompositions, Appl. Math. Comput., 242 (2014), 435–443. [22] J. A. Marrero, R. Ben Taher, M. Rachidi; On explicit formulas for the principal matrix logarithm, Appl. Math. Comput., 220 (2013), 142–148. [23] J. R. Mandujano, L. Verde-Star; Explicit expressions for the matrix exponential function obtained by means of an algebraic convolution formula, Electron. J. Diff. Equ., 2014 No. 79 (2014), 1–7. [24] K. S. Miller, B. Ross; An Introduction to the Fractional Calculus and Fractional Differential Equations, New York, Wiley, 1993. [25] C. Moler, C. Van Loan; Nineteen dubious ways to compute the exponential of a matrix, twenty-five years late, SIAM Rev., 45 No. 1 (2003), 3–49. [26] I. Podlubny; Fractional Differential Equations, 198, San Diego, Academic Press, 1999. [27] M. Popolizio; On the Matrix Mittag–Leffler function: theoretical properties and numerical computation, Mathematics, 27 No. 12 (2019), 1140. [28] E. J. Putzer; Avoiding the Jordan canonical form in the discussion of linear systems with constant coefficients, Amer. Math. Monthly, 73 No. 1 (1966), 2–7. [29] R. F. Rinehart; The equivalence of definitions of a matric function, Amer. Math. Monthly, 62 No. 6 (1955), 395–414. [30] M. R. Rodrigo; On fractional matrix exponentials and their explicit calculation, J. Differential Equations, 261 No. 7 (2016), 4223–4243. [31] B. Ross; Fractional calculus, Math. Mag., 50 No. 3 (1977), 115-122. [32] A. Sadeghi, J. R. Cardoso; Some notes on properties of the matrix Mittag-Leffler function, Appl. Math. Comput., 338 (2018), 733–738. [33] H. G. Sun, Y. Zhang, D. Baleanu, W. Chen, Y. Q. Chen; A new collection of real world ap- plications of fractional calculus in science and engineering, Commun. Nonlinear Sci. Numer. Simul., 64 (2018), 213–231. [34] V. E. Tarasov; Fractional dynamics: applications of fractional calculus to dynamics of par- ticles, fields and media, Springer Science & Business Media, 2011. [35] V. E. Tarasov; Mathematical Economics: Application of Fractional Calculus, Multidisci- plinary Digital Publishing Institute, Basel, 2020. [36] L. Verde-Star; Functions of matrices, Linear Algebra Appl., 406 (2005), 285–300. Antônio Francisco Neto DEPRO, Universidade Federal de Ouro Preto, CEP 35.400-000, Ouro Preto, MG, Brazil Email address: antonio.neto@ufop.edu.br 1. Introduction 2. Statement of main results 3. Auxiliary results 4. Proof of Theorem ?? 5. Connection with previous work 5.1. Theorem ??, Lemma ??, and [Theorem 2]putzer1966avoiding 5.2. Theorem ??, Lemma ??, and [Theorem 1]elaydi1998computation 5.3. Theorem ??, Lemma ??, and [Theorem 3]ahlbrandt2003floquet 5.4. Theorem ??, Lemma ??, and [Theorem 5.1]rodrigo2016fractional 5.5. Theorem ??, Lemma ??, and [Theorem 5.2]rodrigo2016fractional 5.6. Example Conclusion Acknowledgements References