 Applied Mathematics, 2013, 4, 1260-1268 http://dx.doi.org/10.4236/am.2013.49170 Published Online September 2013 (http://www.scirp.org/journal/am) Matrix Functions of Exponential Order Mithat Idemen Engineering Faculty, Okan University, Istanbul, Turkey Email: midemen@gmail.com Received May 31, 2013; revised June 30, 2013; accepted July 7, 2013 Copyright © 2013 Mithat Idemen. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. ABSTRACT Both the theoretical and practical investigations of various dynamical systems need to extend the definitions of various functions defined on the real axis to the set of matrices. To this end one uses mainly three methods which are based on 1) the Jordan canonical forms, 2) the polynomial interpolation, and 3) the Cauchy integral formula. All these methods give the same result, say g(A), when they are applicable to given function g(t) and matrix A. But, unfortunately, each of them puts certain restrictions on g(t) and/or A, and needs tedious computations to find explicit exact expressions when the eigen-values of A are not simple. The aim of the present paper is to give an alternate method which is more logical, simple and applicable to all functions (continuous or discontinuous) of exponential order. It is based on the two-sided Laplace transform and analytical continuation concepts, and gives the result as a linear combination of certain n matri-ces determined only through A. Here n stands for the order of A. The coefficients taking place in the combination in question are given through the analytical continuation of g(t) (and its derivatives if A has multiple eigen-values) to the set of eigen-values of A (numerical computation of inverse transforms is not needed). Some illustrative examples show the effectiveness of the method. Keywords: Matrix; Matrix Functions; Analytical Continuation; Laplace Transform 1. Introduction For many theoretical and practical applications one needs to extend the definitions of functions defined on the real axis to the set of matrices. The history of the subject goes back to the second half of the nineteenth century when Cayley, who is the main instigator of the modern nota-tion and terminology, introduced the concept of the square-root of a matrix A . Since then a huge work has been devoted to the definitions, numerical computations and practical applications of the matrix functions. A rather detailed history (including a large reference list) and important results (especially those concerning the numerical computation techniques) are extensively dis-cussed in the book by Higham . So, we eschew here of making a review of the historical development and giv-ing a large reference list. To extend the definition of a scalar function g(t), de-fined for , to the set of matrices, one starts from an explicit expression of g(t), which can be continued analytically into the complex plane C, and replaces there t by A. If the result is meaningful as a matrix, then it is defined to be g(A). Before going into further detail, it is worthwhile to clarify the meaning of the word “de-fined” appearing in the expression of “defined for ttTh”. It is especially important when g(t) consists of a multi-valued inverse function. To this end consider, for example, the square-root function g(t) = t1/2. Its definition requires, first of all, a cut connecting the branch point t = 0 to the other branch point t =  in the complex plane C. en, by choosing one of its possible values at a given point, for example g(1), one defines it completely. The result consists of a well-defined branch of the square-root function. If one replaces t in this expression by A, then one gets a (unique) matrix to be denoted by A1/2. This matrix satisfies the equation X2 = A which may have many solutions denoted also by A1/2. The above-men- tioned function g(A), which consists merely of the exten-sion of the above-mentioned well-defined branch of the square-root function, can not permit us to find all these solutions. For example the equation X2 = I, where I de-notes the unit 2 × 2 matrix, has infinitely many solutions given by cos sinsin cosX, where  stands for any complex angle. All these matrices are defined to be the square-root I. But the above- mentioned matrix g(A) gives only one of them, namely I Copyright © 2013 SciRes. AM M. IDEMEN 1261or (−I). The known classical methods used in this context are grouped as follows (see , Sections 1, 2): 1) Methods based on the Jordan canonical formula; 2) Methods based on the Hermite interpolation for-mula; 3) Methods based on the Cauchy integral formula. All these methods are applicable when the function g(t), defined on the real axis, can be analytically contin-ued into a domain of the complex-plane, which involves the spectrum of the matrix A (see def. 1.2 and def. 1.4 in ). Consider, for example, the Heaviside unit step func-tion H(t) defined on the real axis by 1, 012, 00, 0.tHt tt (1a) It is obvious that the analytical continuation of H(t) into the complex z-plane, if it is exists, has the point z = 0 as a singular point. To reveal H(z), let us try to find its Taylor expansion about any point a > 0. This expansion is valid in the circle with center at the point z = a and radius equal to r = a. Since all the coefficients except the first one are equal to naught, one gets H(z)  1 at all points inside the circle in question. By letting a   one concludes that H(z) is regular in the right half-plane z > 0. If the above-mentioned Taylor expansion were made about a point a < 0, then one would get H(z)  0 for all z with z < 0. This shows that H(z) is a sectionally regular (holomorphic) function (see , Section 2.15). On the basis of the Plemelj-Sokhotiskii formulas (see , Sec-tion 2.17), for the points on the imaginary axis one writes H = 1/2, which yields 1, 012, 00, 0zHz zz. (1b) Notice that (1a) and (1b) can also be obtained by com-puting the improper integral 221d,2π2zHz z (1c) where the bar on the integral sign stands for the Cauchy principal value. From (1b) one concludes that the seemingly general and elegant method 3), which is based on the Cauchy integral  11d,2πCgiAIA (2) where C stands for a closed contour such that the domain bounded by C involves all the eigen-values of A and the function g(z) is regular there, can not be applicable to find H(A) (and other functions expressible through H(t)) when A has eigen-values having both positive and nega-tive real parts. As to the methods 1) and 2), they need, in general, some tedious and cumbersome computations if A has multiple eigen-values. In the present note we will consider the case when the function g(t), defined on the real axis, is of the exponen-tial order at both t = + and t = −, and give a new method which seems to be more logical and effective especially when the matrix A has multiple eigen-values. It gives the result as a linear combination of n matrices determined only by the matrix A. To this end we consider the Laplace transforms of g(t) on the right and left halves of the real axis, namely:   0ˆed,stgsggt  t (3a)  0ˆedstgsggt  t (3b) and write    11ˆˆ eded,2π2π,.ts tsLLgtgssgsiits (3c) If the orders of the function g(t) for t   and t  (−) are c+ and c−, respectively, then the function ˆgs is a regular function of s in the right-half plane sc and the integration path L+ appearing in (3c) consists of any vertical straight-line located in this half- plane (see Figure 1). Similarly, ˆgs is regular in the half-plane scL and the integration path L− is any vertical straight-line in this half-plane (if c+ < c−, then one can assume L). Furthermore, if g(t) as well as its derivatives up to the order (m-1) are all naught at t = 0, i.e. when  100 0mgg g0 , (4a) then one has c− O c+ ss L+ L− Figure 1. Regularity domains of ˆgs and the integration lines L when c− < c+. Copyright © 2013 SciRes. AM M. IDEMEN 1262 ˆmmgsg s (4b) and ˆ,mmgsg s (4c) which yield inversely  11ˆˆ ed ed,2π2π,.mmts mtsLLgss sgss sgtiit (4d) It is worthwhile to remark here that the formula (4d) permits us to compute mgt only at points on the real axis although g(t) and its derivatives are defined in (or can be continued analytically into) the complex t-plane. Therefore, when the point t is replaced by a complex C , in what follows we will replace the left hand side of (4d) by the analytical continuation of mgt to the point  and write  11ˆˆ ed ed,2π2π.mms msLLgsss gsssgiiC (4e) The formulas (3c), (4d) and (4e) will be the basis of our approach. Let A be a square matrix of dimensions n  n. We will define g(A) by replacing t in (3c) by A, namely:   11ˆˆ edsed.2π2πssLLggsgs siiAAA (5) Thus the computation of g(A) becomes reduced to the computation of exp{At}. As we will see later on, the lat-ter consists of a linear combination of certain constant matrices j ( = order of A). Hence g(A) will also be a linear combination of these j’s for every g(t). It is important to notice that to compute the coefficients in the combinations in question we will never need to compute the transform functions 1, ,jnˆgs as well as the integrals of the form (5) if the analytical continuation of mgt is known at the eigen-values of A (see the ex-amples to be given in Section 4). These points constitute the essential properties of the definition (5): 1) It unifies the definition of g(A) for all functions g(t) of exponential order; 2) It gives an expression of g(A) in terms of certain matrices which take place in the expression of exp(At) and are determined only by A; 3) It reduces the computation of g(A) to the computa-tion of exp(At) together with some scalar constants to be determined in terms of g(t) (and its derivatives when A has multiple eigen-values) at the eigen-values of A. The details are given in the theorems that follow. 2. Basic Results In what follows we will denote a square matrix A of en-tries ajk by A = [ajk]. Here the first and second indices show, respectively, the row and column where ajk is placed. The transpose of A will be indicated, as usual, by a super index T such as . The characteristic polynomial of A will be denoted by TTjkaAf, i.e. detfAI. Theorem-1. Let A = [ajk] be an n  n matrix with char-acteristic polynomial f. Then, 1) when all the zeros of f, say 1,,n, are dis-tinct, one has 01expeαnλtαtΓA (6a) with 0Γ given by  T01, jkfaf Γ (6b) 2) when f has p distinct zeros, say 1,,p, with multiplicities 1,,pmm, respectively, one has  1212 01expepmmmmttt t A  (6c) with the matrices 0, ,1kkmaΓ given by  T111! 1!d.dkmmkmkjkskm ksfsafss Γ (6d) Theorem-2. Let A = [ajk] be an nn matrix while g(z), defined in the complex z-plane, is regular at all the eigen- values of A and its restriction to the real axis is of expo-nential order at both t = + and t = −. If all the eigen- values of A, say 1,,n, are distinct, then one has  01.nggΓA (7) Here stands for the matrix taking place in the expression of exp{At}. 0 1,,aΓnTheorem-3. Let A = [ajk] be an n  n regular matrix which has p distinct eigen-values 1,,p with multi-plicities 1,,pmm1,,; 0ap, respectively, while g(z), defined in the complex z-plane, is regular at all the eigen-values of A and its restriction to the real axis is of exponential or-der at both t = + and t = −. Let the non-zero matrices taking place in the expression of exp{At} be . Then one has , ,kΓKk Copyright © 2013 SciRes. AM M. IDEMEN Copyright © 2013 SciRes. AM 1263G 1212 01,pmmNmmgGG ΓΓ ΓAA (8) where N stands for an integer such that NGtt gt and its derivatives up to the order (K-1) are all naught at t = 0. Let the characteristic polynomial of A be fs: 11 12121 22212............ ... ......nnnn nnas aaaasafsaaaBefore going into detail of proofs of the above-men- tioned theorems, it is worthwhile to draw the attention to the fact that theorems 2 and 3 give the matrix function g(A) as a combination of the n matrices kΓ which appear in the expression of exp(At). They are the same (invariant) for all g(t). s. (10) Then the entry of the inverse matrix 1sAI, which is placed at the k-th row and j-th column, can be computed through the polynomial fs as follows: Proof of theorem-1. Our basic matrix function exp{At} is defined, as usual, through the infinite series  1, ,1,2,,. jkfsjk nafs (11) 221exp1, 2!atatatt  Thus (9) yields by replacing there the scalar constant a by the square matrix A, namely:   T11exped.2πtsjkL fstsi afs A (12)  221exp, .2!tI tttt AAA X Here the integration line L is any vertical straight-line such that all the eigen-values of A are located in the left side of L. It is obvious that X(t) defined as above is the unique solution to the differential equation tXX'A t un- der the initial condition X(0) = I. Hence, by applying the Laplace transform to this equation one gets  1ˆ–ssXAI which permits us to write  1ts1exped.2πLtssi AAIIf the eigen-values are all simple, then the integral in (12) is computed by residues and gives (6a). When some of the eigen-values are multiple, as stated in theorem-1b, the residue method gives (9)   T1 11.1dexp e 1! dmmpstmjksstfamfsss A (13) Proof of theorem-2. When the eigen-values of A are all simple, in (5) one replaces exp(As) by its expression given in (6a) and obtains It is obvious that the derivatives in (13) yields a poly-nomial in t of degree . Hence the final expres-sion of exp(At) can be arranged as what is given in (6c). 1m   T1111ˆˆed ed2π2πnsjk LLsgfgss gsfa ii s A. (14) If all the eigen-values are real, then (3c) reduces (14) to (7). When some or all of the eigen-values are complex, we replace (3c) by (4e) with m = 0 and arrive again (7). Proof of theorem-3. Now consider the case when A has multiple non-zero eigen-values and define Gt Ntgt where the integer N will be determined appro-priately later on. If in (5) one replaces g(t) by G(t) and exp{As} by (6c), then we get 112 2001,pNmmm mgcc c   ΓΓ ΓAA (15a) where 11ˆˆdee2π2πsskkkLLcGsssGssii ds (15b) M. IDEMEN 1264 with 1,,pe coeffici and . Remark that some hents 0, ,1akmkof tΓ may be equalub-index k be Knsi to naught (see vativLet g(t) be the characteristic polynomial of A (i.e. g(t)  ) In order to show the application and effectiveness of the ider some simple One can easily check that A has a tple eigen-value  = 2. Therefore the theorems 1b and 3 are applicable di-re ties p and m mention= 3. On the other handex.-2)rgest s for which one has 0kΓ. Then, by codering the requirements in (4a), we will choose the integer N such that G(t) and its de-ries up to the order (K-1) are all naught for t = 0. In this case all terms existing in (15b) are computed through (4d) or (4e) and give (8). 3. A Corollary (Cayley-Hamilton Theorem) . Let the laf(t)). In this case all the terms taking place in (7) or (8are equal to zero. When the eigen-values are all simple, from (7) one gets directly f(A) = 0, which is valid for both regular and singular matrices. In the case of multi-ple eigen-values, (8) gives f(A) = 0 if A is not singular. We remark that the Cayley-Hamilton theorem is correct for all matrices. We will use this theorem to compute the factor A−N taking place in the formula (8) (see ex.-2). 4. Some Illustrative Examples method, in what follows we will consexamples. Ex.-1 As a first example consider the case where A is as follows: 812 2341122 A. rictly for all functions of exponential order. The quanti- ed in those theorems are: p = 1, m1 from  3812234f 1 2122   one computes Tf2222612204 25 31014224244jka    which gives (see (6c) and (6d)) 2210exp ettttΓΓΓA 2where 01612 2, 361120I and 224011202000Γ. Since 2  0, one has K = 2 which ows that the for-mula (8) is applicable with 0  N  2. For example, in or  shder to find the expressions of A and sin A (with  2), one can choose N = 0 while A, sin A, sin A, 3sin A, arcsinA etc. needs N = 1. To compute cosA, logA, cos A, signA and arcosA one has to choose N = 2. To check the formulas, we would like to compute first integer  2) through the formula (7) which gives An (n = 1112421122nnnnn nn 21 3 2210112 22211221 1242.22nnnnnnnnnnn nnn nnn nnnn  22nnn12n    ΓΓΓA Thus for n = 2 one gets 230 52813224484A. Notice that by a direct multiplication of A with itself one gets the sam result. Similarly, one gets also e 4343 4323334tt t ΓΓΓ ΓAA 1343 3210 2124, 293t ΓI gt tCopyright © 2013 SciRes. AM M. IDEMEN 1265  222221221log loglogloglog2log234log224log2,tgtxtxt xtt xtt xtxxx  ΓΓΓΓAAI 0Γ   22 22212221cos cos coscoscos72cos 2sin 24cos22sin 24cos2,22tgtxtxt xtt xtt xtxx0xxxxxx  ΓΓΓΓAAI Γ      212221sinsin sinsinsin3cos2sin 2cos2sin22sin 2,442 2tgtxtxt xttxtt xtxx xxx xx   ΓΓΓΓΓAAI 0x2 221030 528 signsign24413224sign.484gt t ΓΓΓAAA AI Remark that for different branches of t, 3t and lo one gets different expressions for 3Agxt , gloxA, cosxA and sinxA (sSection 5 ee and theoFinallylet us conse branchetrigonome arcsint and ut as shown in Figure 2 into the re taking place in the formula (8) can be computed rather easily by using the Cayley-Hamilton theorem as follows: rem-4). , ider ths of the inverse tric functions g(t) = h(t) = arccost which map the t-plane cgions in the g- and h-planes shown in Figures 3 and 4 (the so-called principal branches of these functions!). For the first function we have to choose N = 1 while the second one needs N = 2. The matrices A−1 and A−2 162041514 28244A, 256144 321368816 .64 1632 16A Thus, by starting fromd arcsinGt tt anHt 2arccostt one gets from (8) 816 02arcsin4 820 4010 202824 40004802 444086 2045 14 A  and 3264 064 1283256144 3283arccos16 32032 64163688166416 160001632 0163216       A  with arcsin2,arccos2, 1 3.  an interesting exercise to check that It isπarcsin arccos .2 0AAIII70483880 5A. In this case one has which shows again that the theorems 1b and 3 are appli- Ex.-2 Now consider the case where A is as follows:   13f , 2Copyright © 2013 SciRes. AM M. IDEMEN 1266 −1 1 B+ B− C2 A+ A− C1 t Figure 2. Complex plane-t cut along the lines (−) < t < −1 and 1 < t < . g B− C2 C1 A− −/2 /2 A+ B+  arcsin Figure 3. Mapping of the t-plane into the g-plane through the principal branch of the function g = arcsint. A− A+ C1 C2  O h B+ B− arccos  Figure 4. Mapping of the t-plane into the h-plane through the principal branch of the function h = arccost. cable with p = 2, 1 = −1, 2 = 3, m1 = 1 and m2 = 2. From the expression of f() one gets easily Therefore from (6c) and (8) we write .Since K = 0, (8) is applicable with N = 0 and ges, for example, 2Γand  12001012 0 1202,212 ΓΓ 20220 1and 21000000Γ0. 00012300ee ett tΓΓA iv 1210 21201200signsign 302414sign   ΓΓΓΓA sign40 3 12 221021200 HH HH  AΓΓ ΓΓ201212.20 1 By direct multiplication of these matrices by them-selves one gets are the basic properties for the original functions signt and H(t). Similarly, one gets also  22sign, HHAI AA, which 20arcsin2 22222 02  A and 20arccos222222 02  A with arcsin1 ,arcsin3, arcos1 ,arccos3. Here the functions and arcsin tgtarccos t ht consist of the functions consideredk that in the eove. One can easily checxample- 1 ab and πarcsin arccos 2AA, sin arcsinAAIh are the basic properties for the original functions t and arccost. Ex.-3 Finally, we want to give an exame with com-plex-valued eigen-values. To this end consider the case where whicarcsinplabbaA Copyright © 2013 SciRes. AM M. IDEMEN Copyright © 2013 SciRes. AM 1267 has with any real a and b. In this case one  120012001200sign signisigni if 0,if 0ab abaa ΓΓΓΓΓΓAI 12 iabi,ab  , Iand 101i12i1Γ, 201i1.2i1Γ ,if 0a0and Thus, for the following functions one writes     11arcsiniarcsinii arcsiniarcsini22arcsin 11iarcsin iarcsin iarcsinarcsin i22ab abab abab abab ab  A. Here arcsint stands for the function defined in ex.-1 above. It is interesting to check that A. Inverse Power of a Matrix m (= B) which cor-respond to a given matrix A through the relationm = A. The next lemma and theorem concern this case. nction 466gw, 477gw, 588gw. Theorem-4. Let A be an n  n regu matrix with p different eigen-values lar2sign if 0, sinarcsina AIA 1pn p differentwhile m is an integer. Then there exist at least m matrices B such that Bm = A. co5Proof. Let the eigen-values of A be 1, ,Consider a Riemann surface which the function jjp. rresponds toFrom the example 1 considered above one observes that there are many inverse powers A1/1mgtt and denote, one k-th sheet of this Riemann surface, the values of the function g(t) at th j by Bpoints1, 1,,kmkjjwk m1mgtt, Lemma. Let the values of the fu  . On the basis of the lemma, we can ange the cut line appropriately such that the values of tm at the points defined on a Riemann surface, at given p points 12,,,p be 0kjjgw, 1, ,km. Then there is a Riemann surface such that on one of its sheets g(t) takes previously chosen values at the points 1,,,ch1/,,,12p be-co me equal to previously chosen values. For example, 2p, namely: kjjgw with arbitrarily chosen jk2Proof. Let at the first r points (1  r < p) o1, 2,,, m1,,,j. pne has 1jjgw while at 1r one wants to have 11rrqgw, where ber. ne the ne which e branch pocirc2, ,qmint t = 0, en is any numw cut line as a spiral curveles the points Then, we defistarts from th,,,12 r (see betwthe Fi q-1) times passes gure 4) and (een the points r and 1r. Thus the analytic con-tinuation of 111gwto j 2, ,jr b ecomes 1jw while at 1r it is e1r. We continue this process to adjust also the values at the points 2,,rnqual to qw. Nowhen the super index kj + 1 in 11 tice that jkjw is smaller than kj, we can replace kj + 1 (kj + 1 + m) because one has 111jkjjww. At the end at a plane cut along an appropriate (spiral) curve joining the oint 0 to t =  such that at the given points ,, by mwe arrivept =  1jk12 ,p the function g(t) has desired values jkjw. xample on For ethe Riemann sheet shown in Figure 5 one has 111gw, 1122 55,, ,gwg w 1111mw, 1322mw, 1233mw, 14m 14w etc. Since these valun be arranged in mp dif-fees car (8)) can hav rent forms, the right-hand side of (7) (oe mp different values. This proves thorem-4. Ex. LetaC and bC any numbers which differ from zero. Then the matrix 2200abA 1 7 6 5 O t          8 Figure 5. A cut appropriate for the particular case when r = 5 and q = 4. M. IDEMEN 1268 has four square-roots given as follows: 0bREFERENCES  A. Cayley, “A Memoir on the Theory of Matrices,” Phi-losophical Transactions of the Royal Society of LondonVol. 148, 1858, pp. 17-37. doi:10.1098/rstl.1858.0002000, , , 00 0 0aaa abbb   . ,  N. J. Higham, “Functions of Matrices: Theory - putation,” Society for Industrial and Applied Mathe- matica (SIAM), Philadelphia, 2008.  N. I. Mushkhral Equations,” P. and Comelishvili, “Singular IntegNoordhoff Ltd., Holland, 1958.  E. C. Titchmarsh, “Introduction to the Theory of Fourier Integrals,” Chelsea Publishing Company, 1986, Chapter 1.3. Copyright © 2013 SciRes. AM