MATRICES AND MATRIX OPERATIONS: Unit 23

Dr. Wlodzislaw Kostecki

The Papua New Guinea University of Technology (PNGUT)

Department of Electrical and Communication Engineering

Lae, Morobe Province

Papua New Guinea

Copyright © 2000 by Wlodzislaw Kostecki

All rights reserved

-------------------------------------------------------------------

(23) Functions of constant matrices

OBJECTIVES :

• To provide essential information concerning expansion of functions of one variable using Taylor’ s and Maclaurin’ s series.

• To introduce the concepts of an interval and circle of convergence of a series representing a function of a real and complex variable, respectively.

• To provide definitions of the matrix series for real-valued and complex-valued matrices, based on the Maclaurin series.

• To state the condition necessary for Maclaurin’ s series of a square matrix to be absolutely convergent.

• To introduce the concept of a well-defined function of a matrix.

• To provide a list of some common functions together with the interval of convergence of their corresponding Maclaurin’ s series for reference.

• To provide a reference example of computation of a trigonometric function of a real and complex variable using infinite Taylor’ s series.

• To introduce the only built-in function exponential for computation of the exponential function of matrices.

• To provide six step-by-step solved examples for computation of functions of matrices with real and complex elements using truncated Maclaurin’ s series.

• To show how to pre-process a matrix containing transcendental numbers to significantly reduce the computation time of its function.

• To examine six common functions of the zero matrix.

• To show how to obtain the exact result of the exponential function of the unit matrix and examine two other functions of this matrix.

• To examine seven cases of exponentiation of the imaginary unit matrix multiplied by some scalars involving the various multiplicities of the imaginary unit.

• To show how to make use of the Cayley-Hamilton theorem for computation of matrix functions.

• To provide three step-by-step solved examples for computation of functions of matrices with real elements using a method based on the Cayley-Hamilton theorem and verify the results using the truncated Maclaurin’ s series.

• To state the necessary conditions that a matrix must satisfy in order that its square root be unique .

• To provide two step-by-step solved examples for computation of the square root of real and complex matrices using the method based on the Cayley-Hamilton theorem and verify the results using the truncated Maclaurin’ s series.

• To determine the conditions that a matrix must satisfy in order that its unique square root could be computed in an alternative manner using the functions natural logarithm and exponential.

• To provide a number of variants of the various square matrices of second order that are each a non-unique square root of the unit matrix of the same order.

• To provide two step-by-step solved examples for computation of functions of diagonal matrices with real elements using similarity of matrices and verify the results using the truncated Maclaurin’ s series.

• To provide four step-by-step solved examples for computation of functions of non-diagonal real matrices with distinct eigenvalues using similarity of matrices and verify the results using the truncated Maclaurin’ s series.

> restart : with(linalg, augment, coldim, definite, diag, eigenvals, eigenvects, exponential, inverse, jordan, matadd, rowdim) :

Functions of matrices are defined only for square matrices.

This Unit deals with computation of functions for constant matrices whose elements are real or complex numbers.

Consider a function of a real variable, , and a function of a complex variable, , and assume that either function has a Taylor series expansion about the points and , respectively, given by

and

Suppose that there is an interval with the central point at and end points at and , within which the series representing is absolutely convergent for all values of , i.e. those which satisfy the inequality .

Suppose that there is a circle in the complex plane with the centre at , radius , and equation , within which the series representing is absolutely convergent for all values of , i.e. those which satisfy the inequality .

Setting, for convenience, and in the above Taylor’ s series yields the Maclaurin series, viz.

and

Similarly as for Taylor’ s series, it may be assumed that the Maclaurin series representing is absolutely convergent for and the series representing is absolutely convergent for .

The definitions of the matrix series for a real-valued matrix [ A ] and a complex-valued matrix [ Z ] are based on the Maclaurin series, viz.

and

Either series is absolutely convergent if the absolute value or modulus of every eigenvalue of [ A ] or the modulus of every eigenvalue of [ Z ] satisfies the respective condition

or

In such a case, the functions f( [ A ] ) and f( [ Z ] ) are said to be well defined for the matrix and may be computed using a suitable matrix series.

Maple has built-in procedures for computation of infinite Taylor’ s series of functions of one variable. For example, it uses the following series in computing :

> `cos(x)` := Sum((-1)^n*x^(2*n)/(2*n)!, n=0..infinity) : cos(x) = `cos(x)` ; Cos(x) := `cos(x)` :

Maple evaluates this series to an internally adopted accuracy and returns the result with the number of digits determined by the environment variable Digits (the default is ). For example, set and and compute and to digits. Thus,

> Digits := 15 : x := 2.7 : `cos(x)` := evalf(`cos(x)`) : x := 'x' : cos(x) = `cos(x)` ; z := 2.7 + 1.2*I : `cos(z)` := evalf(subs(x=z, Cos(x))) : z := 'z' : cos(z) = `cos(z)` ; Digits := 10 :

No such built-in procedures exist in Maple for computation of matrix functions. The only exception is the exponential function, which is accessible through the linalg package.

Therefore, where computation of a matrix function is involved, the series must be truncated. The number of series terms adopted should result from a compromise between the computation time and required accuracy.

It is appropriate to list at this point some common functions together with the interval of convergence of their corresponding Maclaurin’ s series, which will be a good indication whether the sought-for function can be computed for a given matrix once its eigenvalues have been determined.

Although functions of a real variable are shown together with the convergence interval of their series, substitute with and the term " interval of convergence " with " circle of convergence " where complex-numbered matrices are involved. If a real-valued matrix has complex eigenvalues, the designation still holds but the circle of convergence comes into play.

(1) Functions whose series is convergent for every , i.e. having convergence interval . These functions are well defined for every square matrix:

• , , , , ,

(2) Functions with convergence interval

• , , ,

(3) Functions with convergence interval 0 <

• , , ,

(4) Functions with convergence interval

• , , , , ,

(5) Function with convergence interval

•

(6) Function with convergence interval

•

(7) Function with convergence interval <

• as given by the series

Note that there are also other series representations of with their respective, different convergence intervals .

Apart from computation of matrix functions based on the Maclaurin series representation of a given function, there are two other methods that may be used for this purpose. One of them is based on the Cayley-Hamilton theorem [refer to Unit (21)] and the other on similarity of matrices [refer to Unit (22)].

All the three methods together with illustrating examples are presented in Sections A , B , and C of this Unit.

A. Computation of matrix functions using the Maclaurin series

Below presented are examples for computation of several different matrix functions with the use of this method.

Example 1

Compute the natural logarithm of ( × ) matrices [ A ] and [ B ] given as

> A := matrix(2, 2, [2, 3, 0, 2]) : B := matrix(2, 2, [1, -2, 1, 2]) : A = matrix(A) ; B = matrix(B) ;

Step 1 . Find the eigenvalues of either matrix:

• for matrix [ A ]

> charroots(A) := eigenvals(A) : char_roots(A) = charroots(A) ;

Since neither of the characteristic roots is greater than , the function ln( [ A ] ) is well defined for this matrix and may be computed.

• for matrix [ B ]

> charroots(B) := eigenvals(B) : char_roots(B) = charroots(B) ;

> root1(B) := charroots(B)[1] : char_root[1](B) = root1(B) ;

> `abs(root1(B))` := abs(root1(B)) : Abs(char_root[1](B)) = `abs(root1(B))` ;

Since the modulus of neither characteristic root is greater than , the function ln( [ A ] ) is well defined for this matrix and may be computed.

Step 2 . Write the expression for a Maclaurin’ s series representing the function , substitute the matrix names for , and replace subtrahend with the name of the unit matrix [ U ]:

> `ln(A)` := subs((x-1)=(A-U), Sum((-1)^(n+1)*(x-1)^n/n, n=1..infinity)) : ln(A) = `ln(A)` ;

> `ln(B)` := subs((x-1)=(B-U), Sum((-1)^(n+1)*(x-1)^n/n, n=1..infinity)) : ln(B) = `ln(B)` ;

Step 3 . Bearing in mind that the unit matrix [ U ] appropriately sized for this case is

> U := diag(1, 1) : U = matrix(U) ;

substitute either matrix with its numerical values and the unit matrix:

> `ln(A)` := subs(A=matrix(A), U=matrix(U), `ln(A)`) : ln(A) = `ln(A)` ;

> `ln(B)` := subs(B=matrix(B), U=matrix(U), `ln(B)`) : ln(B) = `ln(B)` ;

Step 4 . Evaluate the function of either matrix after truncating the series to the first terms:

> `ln(A)` := evalf(evalm(sum((-1)^(n+1)*(A - U)^n/n, n=1..50))) : ln(A) = matrix(`ln(A)`) ;

> `ln(B)` := evalf(evalm(sum((-1)^(n+1)*(B - U)^n/n, n=1..50))) : ln(B) = matrix(`ln(B)`) ;

Example 2

Compute the cosine of ( × ) matrices [ A ] and [ Z ] given as

> A := matrix(2, 2, [2, 4, 1, 2]) : Z := matrix(2, 2, [1+5*I, 2-4*I, 3+2*I, 1-2*I]) :

> A = matrix(A) ; Z = matrix(Z) ;

Step 1 . Find the eigenvalues of either matrix:

• for matrix [ A ]

> charroots(A) := eigenvals(A) : char_roots(A) = charroots(A) ;

• for matrix [ Z ]

> charroots(Z) := eigenvals(Z) : char_roots(Z) = charroots(Z) ;

> root1(Z) := charroots(Z)[1] : char_root[1](Z) = root1(Z) ;

> `abs(root1(Z))` := evalf(abs(root1(Z))) : Abs(char_root[1](Z)) = `abs(root1(Z))` ;

Note that this step is not necessary since the series of is convergent for all , which implies that cos( [ A ] ) is well defined for every square matrix. Similarly, the series of is convergent for all , which implies that cos( [ Z ] ) is well defined for every complex-valued square matrix.

Step 2 . Write the expression for a Maclaurin’ s series representing the function and and substitute the matrix names [ A ] and [ Z ] for and , respectively:

> `cos(A)` := subs(x=A, Sum((-1)^n*x^(2*n)/(2*n)!, n=0..infinity)) : cos(A) = `cos(A)` ;

> `cos(Z)` := subs(x=Z, Sum((-1)^n*x^(2*n)/(2*n)!, n=0..infinity)) : cos(Z) = `cos(Z)` ;

N.B. Maple returns the first term of either series as the scalar value , according to the result of the following computation performed exemplarily for the matrix [ A ]:

> `A^0` := evalm(A^0) : A ^ `0` = `A^0` ;

Therefore, still with [ A ] as an example, the corresponding series for to, e.g., would be of the form

> cos(A) = subs(x=A, sum((-1)^n*x^(2*n)/(2*n)!, n=0..5)) ;

in which the first term is the algebraic number .

However, adhering to the common literature convention that [ A ] ^ 0 = [ U ], the first term of the series should be the unit matrix [ U ]. To obtain this result, the original expression for the series is modified by

(a) replacing with in the sum function, and

(b) adding the unit matrix [ U ] as the first term of the series representation of the matrix function.

Consequently, the following expressions are adopted in this case:

> `cos(A)` := U + subs(x=A, Sum((-1)^n*x^(2*n)/(2*n)!, n=1..infinity)) : cos(A) = `cos(A)` ;

> `cos(Z)` := U + subs(x=Z, Sum((-1)^n*x^(2*n)/(2*n)!, n=1..infinity)) : cos(Z) = `cos(Z)` ;

Step 3 . With the unit matrix as in Example 1, substitute the matrix names with their corresponding numerical structures:

> `cos(A)` := subs(U=matrix(U), A=matrix(A), `cos(A)`) : cos(A) = `cos(A)` ;

> `cos(Z)` := subs(U=matrix(U), Z=matrix(Z), `cos(Z)`) : cos(Z) = `cos(Z)` ;

Step 4 . Evaluate the function of either matrix after truncating the series to the first terms:

> `cos(A)` := evalf(evalm(U + sum((-1)^n*A^(2*n)/(2*n)!, n=1..50))) : cos(A) = matrix(`cos(A)`) ;

> `cos(Z)` := evalf(evalm(U + sum((-1)^n*Z^(2*n)/(2*n)!, n=1..50))) : cos(Z) = matrix(`cos(Z)`) ;

Example 3

Compute the exponential function of the ( × ) matrices [ A ] and [ Z ] given as before.

Step 1 . Find the eigenvalues of either matrix:

This step is not necessary since the series of is convergent for all , which implies that exp( [ A ] ) is well defined for every square matrix. Likewise, the series of is convergent for all , which implies that exp( [ Z ] ) is well defined for every square matrix with complex-numbered elements.

Step 2 . Write the expression for a Maclaurin’ s series representing the function and and substitute the matrix names [ A ] and [ Z ] for and , respectively:

> `exp(A)` := subs(x=A, Sum(x^n/n!, n=0..infinity)) : exp(A) = `exp(A)` ;

> `exp(Z)` := subs(x=Z, Sum(x^n/n!, n=0..infinity)) : exp(Z) = `exp(Z)` ;

N.B. Like in Example 2, Maple returns the first term of either series as the scalar value . With the matrix [ A ] as an example, the corresponding series for to, e.g., would be of the form

> exp(A) = subs(x=A, sum(x^n/n!, n=0..5)) ;

For the reasons set forth in Example 2, the following expressions are adopted in this case:

> `exp(A)` := U + subs(x=A, Sum(x^n/n!, n=1..infinity)) : exp(A) = `exp(A)` ;

> `exp(Z)` := U + subs(x=Z, Sum(x^n/n!, n=1..infinity)) : exp(Z) = `exp(Z)` ;

Step 3 . With the unit matrix as in Example 1, substitute the matrix names with their corresponding numerical structures:

> `exp(A)` := subs(U=matrix(U), A=matrix(A), `exp(A)`) : exp(A) = `exp(A)` ;

> `exp(Z)` := subs(U=matrix(U), Z=matrix(Z), `exp(Z)`) : exp(Z) = `exp(Z)` ;

Step 4 . Evaluate the function of either matrix after truncating the series to the first terms:

> `exp(A)` := evalf(evalm(U + sum(A^n/n!, n=1..50))) : exp(A) = matrix(`exp(A)`) ;

> `exp(Z)` := evalf(evalm(U + sum(Z^n/n!, n=1..50))) : exp(Z) = matrix(`exp(Z)`) ;

* * *

N.B. The exponential function of a square matrix may be computed in Maple directly , using the built-in exponential function. For the matrices of Example 3, it yields

> `exp(A)` := evalf(exponential(A)) : exp(A) = matrix(`exp(A)`) ;

> `exp(Z)` := evalf(exponential(Z)) : exp(Z) = matrix(`exp(Z)`) ;

Notice that computation using the exponential function is much faster.

[ Refer also to Units (24), (25), and (30) in which examples for the major applications of the exponential function are given. ]

* * *

Example 4

Compute the arc tangent of ( × ) matrices [ A ] and [ B ] given as

> A := matrix(2, 2, [1/3, -2/7, -1/6, 2/3]) : B := matrix(2, 2, [1/3, -2/7, 1/6, 2/3]) :

> A = matrix(A) ; B = matrix(B) ;

Step 1 . Find the eigenvalues of either matrix:

• for matrix [ A ]

> charroots(A) := eigenvals(A) : char_roots(A) = charroots(A) ;

> root1(A) := charroots(A)[1] : root2(A) := charroots(A)[2] :

> char_root[1](A) = root1(A) ; char_root[2](A) = root2(A) ;

> `abs(root1(A))` := evalf(abs(root1(A))) : `abs(root2(A))` := evalf(abs(root2(A))) :

> Abs(char_root[1](A)) = `abs(root1(A))` ; Abs(char_root[2](A)) = `abs(root2(A))` ;

Since either of the characteristic roots is smaller than , the function arctan( [ A ] ) is well defined for this matrix and may be computed.

• for matrix [ B ]

> charroots(B) := eigenvals(B) : char_roots(B) = charroots(B) ;

> root1(B) := charroots(B)[1] : root2(B) := charroots(B)[2] :

> char_root[1](B) = root1(B) ; char_root[2](B) = root2(B) ;

> `abs(root1(B))` := evalf(abs(root1(B))) : Abs(char_root[1](B)) = `abs(root1(B))` ;

Since the modulus of either characteristic root is less than , the function arctan( [ B ] ) is well defined for this matrix and may be computed.

Step 2 . Write the expression for a Maclaurin’ s series representing the function and substitute the matrix names [ A ] and [ B ] for :

> `arctan(A)` := subs(x=A, Sum((-1)^n*x^(2*n+1)/(2*n+1), n=0..infinity)) : arctan(A) = `arctan(A)` ;

> `arctan(B)` := subs(x=B, Sum((-1)^n*x^(2*n+1)/(2*n+1), n=0..infinity)) : arctan(B) = `arctan(B)` ;

Step 3 . Substitute the matrix names with their corresponding numerical structures:

> `arctan(A)` := subs(A=matrix(A), `arctan(A)`) : arctan(A) = `arctan(A)` ;

> `arctan(B)` := subs(B=matrix(B), `arctan(B)`) : arctan(B) = `arctan(B)` ;

Step 4 . Evaluate the function of either matrix after truncating the series to the first terms:

> `arctan(A)` := evalf(evalm(sum((-1)^n*A^(2*n+1)/(2*n+1), n=0..50))) :

> arctan(A) = matrix(`arctan(A)`) ;

> `arctan(B)` := evalf(evalm(sum((-1)^n*B^(2*n+1)/(2*n+1), n=0..50))) :

> arctan(B) = matrix(`arctan(B)`) ;

Example 5

Compute the arc tangent of ( × ) matrices [ A ] and [ Z ] given as

> A := matrix(2, 2, [0.2, -0.4, 0.1, 0.3]) : Z := matrix(2, 2, [0.1+0.5*I, 0.2-0.4*I, 0.3+0.2*I, 0.1-0.2*I]) :

> A = matrix(A) ; Z = matrix(Z) ;

Step 1 . Find the eigenvalues of either matrix:

• for matrix [ A ]

> charroots(A) := eigenvals(A) : char_roots(A) = charroots(A) ;

> root1(A) := charroots(A)[1] : root2(A) := charroots(A)[2] :

> char_root[1](A) = root1(A) ; char_root[2](A) = root2(A) ;

> `abs(root1(A))` := abs(root1(A)) : Abs(char_root[1](A)) = `abs(root1(A))` ;

Since the modulus of either characteristic root is less than , the function arctan( [ A ] ) is well defined for this matrix and may be computed.

• for matrix [ Z ]

> charroots(Z) := eigenvals(Z) : char_roots(Z) = charroots(Z) ;

> root1(Z) := charroots(Z)[1] : root2(Z) := charroots(Z)[2] :

> char_root[1](Z) = root1(Z) ; char_root[2](Z) = root2(Z) ;

> `abs(root1(Z))` := abs(root1(Z)) : `abs(root2(Z))` := abs(root2(Z)) :

> Abs(char_root[1](Z)) = `abs(root1(Z))` ; Abs(char_root[2](Z)) = `abs(root2(Z))` ;

Since the modulus of either characteristic root is less than , the function arctan( [ Z ] ) is well defined for this matrix and may be computed.

Step 2 . Write the expression for a Maclaurin’ s series representing the function and and substitute the matrix names [ A ] and [ Z ] for and , respectively:

> `arctan(A)` := subs(x=A, Sum((-1)^n*x^(2*n+1)/(2*n+1), n=0..infinity)) : arctan(A) = `arctan(A)` ;

> `arctan(Z)` := subs(x=Z, Sum((-1)^n*x^(2*n+1)/(2*n+1), n=0..infinity)) : arctan(Z) = `arctan(Z)` ;

Step 3 . Substitute the matrix names with their corresponding numerical structures:

> `arctan(A)` := subs(A=matrix(A), `arctan(A)`) : arctan(A) = `arctan(A)` ;

> `arctan(Z)` := subs(Z=matrix(Z), `arctan(Z)`) : arctan(Z) = `arctan(Z)` ;

Step 4 . Evaluate the function of either matrix after truncating the series to the first terms:

> `arctan(A)` := evalm(sum((-1)^n*A^(2*n+1)/(2*n+1), n=0..50)) : arctan(A) = matrix(`arctan(A)`) ;

> `arctan(Z)` := evalm(sum((-1)^n*Z^(2*n+1)/(2*n+1), n=0..50)) : arctan(Z) = matrix(`arctan(Z)`) ;

Example 6

Compute the third power of sine of a ( × ) matrix [ A ] given as

> A := matrix(2, 2, [-3, 5, 1, 9]) : A = matrix(A) ;

Step 1 . Find the eigenvalues of the matrix:

This step is not necessary since the series of is convergent for all , which implies that sin( [ A ] ) is well defined for every square matrix. So is {sin( [ A ] )}^3 and it may be computed.

Step 2 . Write the expression for a Maclaurin’ s series representing the function and substitute the matrix name [ A ] for :

> `sin(A)` := subs(x=A, Sum((-1)^n*x^(2*n+1)/(2*n+1)!, n=0..infinity)) : sin(A) = `sin(A)` ;

Step 3 . Substitute the matrix name with its corresponding numerical structure:

> `sin(A)` := subs(A=matrix(A), `sin(A)`) : sin(A) = `sin(A)` ;

Step 4 . Evaluate the matrix function after truncating the series to the first terms:

> `sin(A)` := evalf(evalm(sum((-1)^n*A^(2*n+1)/(2*n+1)!, n=0..50))) : sin(A) = matrix(`sin(A)`) ;

Step 5. Compute the third power of the matrix function:

> `(sin(A))^3` := evalm(`sin(A)` ^ 3) : (sin(A))^3 = matrix(`(sin(A))^3`) ;

N.B. Since Maclaurin’ s series representation of the function is known, viz.

> sin(x)^3 = (3/4)*Sum((-1)^n*(1-9^n)*x^(2*n+1)/(2*n+1)!, n=0..infinity) ;

it may be used for direct computation of the function {sin( [ A ] )}^3 by following similar steps as above.

Evaluation of the matrix function after truncating the series to the first terms yields the matrix

> `(sin(A))^3` := evalf(evalm((3/4)*sum((-1)^n*(1-9^n)*A^(2*n+1)/(2*n+1)!, n=0..50))) :

> (sin(A))^3 = matrix(`(sin(A))^3`) ;

which is practically equal to the matrix obtained by computing the third power of sin( [ A ] ) .

* * *

N.B. Various combinations of matrix functions may also be computed using the method based on Maclaurin’ s series representation of the functions concerned. Where each function in a combination has different convergence interval (or circle), the resultant function will be convergent within the narrowest interval.

Exemplarily, consider the product where the series for either function is convergent for every , i.e. has convergence interval . Both functions are well defined for every square matrix and so is their product. Using the same matrix [ A ] as before, the exponential function, and the expression for cos( [ A ] ) obtained in Example 2 yields

> `exp(A) cos(A)` := evalm(evalf(exponential(A)) &* evalf(evalm(U + sum((-1)^n*A^(2*n)/(2*n)!, n=1..50)))) : exp(A)*cos(A) = matrix(`exp(A) cos(A)`) ;

* * *

N.B. Transcendental numbers in a matrix are normally accepted and processed by the sum function, but the computation time is prohibitive . As an example, consider a ( × ) matrix [ A ] given as

> A := matrix(2, 2, [exp(1), Pi, sin(Pi/9), sqrt(3)]) : A = matrix(A) ;

and compute the cosine of the matrix using the relevant series expression obtained in Example 2.

Please wait. The following computation will take several tens of seconds.

> `cos(A)` := evalf(evalm(U + sum((-1)^n*A^(2*n)/(2*n)!, n=1..50))) : cos(A) = matrix(`cos(A)`) ;

Floating-point evaluation of the matrix elements before summation operations reduces the computation time significantly without any practical loss of accuracy, as demonstrated below.

> A := evalf(matrix(A)) : A = matrix(A) ;

> `cos(A)` := evalf(evalm(U + sum((-1)^n*A^(2*n)/(2*n)!, n=1..50))) : cos(A) = matrix(`cos(A)`) ;

Notice the insignificant differences in the values of the elements of both resultant matrices for . This is because of the floating-point evaluation of the elements of [ A ] prior to summation.

* * *

N.B. Notice the results of computation of some functions of the zero matrix [ 0 ]. Exemplarily, consider the ( × ) matrix:

> `0` := matrix(2, 2, [0, 0, 0, 0]) : 0 = matrix(`0`) ;

> charroots(`0`) := eigenvals(`0`) : char_roots(0) = charroots(`0`) ;

(1) The exponential function of [ 0 ] is the ( × ) unit matrix:

> `exp(0)` := eval(evalm(U + sum(`0`^n/n!, n=1..2))) : exp(`0`) = matrix(`exp(0)`) ;

or, using the exponential function,

> `exp(0)` := exponential(`0`) : exp(`0`) = matrix(`exp(0)`) ;

(2) The sine function of [ 0 ] is the same ( × ) zero matrix:

> `sin(0)` := eval(evalm(sum((-1)^n*`0`^(2*n+1)/(2*n+1)!, n=0..1))) : sin(`0`) = matrix(`sin(0)`) ;

(3) The cosine function of [ 0 ] is the ( × ) unit matrix:

> `cos(0)` := eval(evalm(U + sum((-1)^n*`0`^(2*n)/(2*n)!, n=1..2))) : cos(`0`) = matrix(`cos(0)`) ;

(4) The arc tangent function of [ 0 ] is the same ( × ) zero matrix:

> `arctan(0)`:=eval(evalm(sum((-1)^n*`0`^(2*n+1)/(2*n+1), n=0..1))) : arctan(`0`)=matrix(`arctan(0)`) ;

(5) The hyperbolic sine function of [ 0 ] is the same ( × ) zero matrix:

> `sinh(0)` := eval(evalm(sum(`0`^(2*n+1)/(2*n+1)!, n=0..1))) : sinh(`0`) = matrix(`sinh(0)`) ;

(6) The hyperbolic cosine function of [ 0 ] is the ( × ) unit matrix:

> `cosh(0)` := eval(evalm(U + sum(`0`^(2*n)/(2*n)!, n=1..2))) : cosh(`0`) = matrix(`cosh(0)`) ;

Notice the number of series terms adopted in computation of the above functions of matrix [ 0 ]: n = 0...1 or n = 1...2. These numbers are sufficient because every term for equals , as shown exemplarily below for the functions exp( [ 0 ] ) and sinh( [ 0 ] ):

> exp(`0`) = matrix(U) + sum(matrix(`0`)^n/n!, n=1..2) ;

> sinh(`0`) = sum(matrix(`0`)^(2*n+1)/(2*n+1)!, n=0..2) ;

* * *

N.B. Notice the result of computation of with the unit matrix given as

> U := diag(1, 1) : U = matrix(U) ;

> charroots(U) := eigenvals(U) : char_roots(U) = charroots(U) ;

> `exp(U)` := evalf(evalm(U + sum(U^n/n!, n=1..50))) : exp(U) = matrix(`exp(U)`) ;

The exact result of may be obtained only if the exponential function is used, i.e.

> `exp(U)` := exponential(U) : exp(U) = matrix(`exp(U)`) ;

* * *

N.B. Notice the result of computation of ln( [ U ] ) with the unit matrix given as before:

> `ln(U)` := evalm(sum((-1)^(n+1)*(matadd(U, -U))^n/n, n=1..2)) : ln(U) = matrix(`ln(U)`) ;

* * *

N.B. Notice the result of computation of [ U ] ^(1/2) for the same unit matrix:

> `sqrt(U)` := exponential(`ln(U)`, 1/2) : U^`1/2` = matrix(`sqrt(U)`) ;

[ For details on computation of the function square root of a matrix, refer to Section B , Example 4. ]

* * *

N.B. Notice the following results of exponentiation of the ( × ) imaginary unit matrix [ IU ] multiplied by some scalars.

> IU := diag(I, I) : IU = matrix(IU) ;

> charroots(IU) := eigenvals(IU) : char_roots(IU) = charroots(IU) ;

(1) The exponential function of [ IU ]:

> `exp(IU)` := exponential(IU) : exp(IU) = matrix(`exp(IU)`) ;

or, as a floating-point approximation,

> exp(IU) = evalf(matrix(`exp(IU)`)) ;

(2) The exponential function of the product of [ IU ] and the imaginary unit:

> `exp(IU I)` := exponential(IU*I) : exp(IU*I) = matrix(`exp(IU I)`) ;

or, as a floating-point approximation,

>