MATRICES AND MATRIX OPERATIONS: Unit 24

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

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

(24) Functions of matrices comprising linear functions

OBJECTIVES :

• To provide five step-by-step solved examples for computation of matrix functions with real and complex elements using in each case all the three methods of computation of matrix functions that are employed in Unit (23).

• To provide a comparison of the resultant matrix functions by computing the function of each matrix for a given numerical value of the variable.

• To include computation of the limit of the resultant matrix function having a discontinuity on the interval of convergence of the function.

• To point out a disadvantage of computing matrix functions using truncated Maclaurin’ s series if a display of the resultant matrix was essential.

• To reveal a difficulty when Maple does not return any symbolic solutions to a system of transcendental equations and show how to get around such a problem.

> restart : with(linalg, augment, coldim, diag, eigenvals, eigenvects, equal, exponential, rowdim) :

This Unit deals with computation of functions of square matrices whose elements are linear functions of type , where t is a real variable and and are constant coefficients.

The exponential function of matrices whose elements are linear functions of type , or exp([ A ] t ), is a special case. It is treated separately in more detail in Unit (30) together with specific applications of the function.

A function f of a matrix [ A ( ) ] whose elements are linear functions of a real variable t may be computed only if the function is well defined for the matrix. [ Refer to Unit (23) for the well-defined function. ]

A well-defined function of a matrix [ A ( ) ] may be computed using the Cayley-Hamilton theorem, similarity of matrices, or the Maclaurin series representation of the corresponding function .

However, it should be borne in mind that none of the methods is universal since each has limitations of different nature – refer to Units (22) and (23). Whether a given method is suitable for computation of a matrix function, may be assessed only when both the sought-for function and the matrix are known.

In the examples provided hereunder, the matrices have been carefully constructed such that all the three methods may be used for computation of functions f( [ A ( ) ] ) that are well defined for these matrices and can be determined without a great computational effort. This has been done for illustrating and comparative purposes.

Example 1

Compute cos( [ A ( ) ] ) for a ( × ) matrix [ A ( ) ] given as

> A(t) := matrix(2, 2, [t+2, 0, -t, 2*t-1]) : 'A(t)' = A(t) ;

The power series of is convergent for all [refer to Unit (23)], which implies that cos( [ A ] ) is well defined for every square matrix. This, in turn, means that the function cos( [ A ( ) ] ) is well defined for every value of t , so it may be computed for this matrix. See also the note at the end of Method 1.

For convenience in computations, let

> B(t) := A(t) : A(t) := 'A(t)' :

Method 1 . Using the Cayley-Hamilton theorem

Step 1 . Having in mind that for this matrix, write equation (1) of Unit (23) B for [ A ( ) ]:

> `cos(A(t))` := a[1] * A(t) + a[0]*U : cos(A(t)) = `cos(A(t))` ; Cos(A(t)) := `cos(A(t))` :

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

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

evaluate the matrix equation for cos( [ A ( ) ] ):

> `cos(A(t))` := evalm(subs(A(t)=B(t), `cos(A(t))`)) : cos(A(t)) = matrix(`cos(A(t))`) ;

Step 3 . Formulate equation (2) of Unit (23) B corresponding to this case:

> l := lambda : `r(l)` := subs(A(t)=l, U=1, Cos(A(t))) : r(l) = `r(l)` ;

Step 4 . Determine the eigenvalues of [ A ( ) ]:

> charroots(A(t)) := eigenvals(subs(A(t)=B(t), A(t))) : char_roots(A(t)) = charroots(A(t)) ;

Extracting the individual eigenvalues yields

> No_roots(A(t)) := nops([charroots(A(t))]) :

> for i to No_roots(A(t)) do ch_r[i](A(t)) := charroots(A(t))[i] : print(l[i] = ch_r[i](A(t))) : od :

Step 5 . Taking notice of the fact that the roots are distinct, formulate equation (3) of Unit (23) B in a general form:

Since the function corresponding to this case is

> `f(l)` := cos(l) : f(l) = `f(l)` ;

therefore, equation (3) assumes the form

> Eq3_f(l) := `f(l)` = `r(l)` : Eq3_f(l) ;

Step 6 . Obtain a set of equations by substituting either eigenvalue into equation (3):

> for i to No_roots(A(t)) do Eq3_f[i](l) := subs(l=ch_r[i](A(t)), Eq3_f(l)) : print(Eq3_f[i](l)) : od :

Step 7 . Solve the simultaneous equations for the unknown coefficients:

> solution := solve({Eq3_f[1](l), Eq3_f[2](l)}, {a[0], a[1]}) : solution ;

Extracting either unknown from the solution set yields

> assign(solution) : for i to No_roots(A(t)) do a[i-1] := a[i-1] : print(evaln(a[i-1]) = a[i-1]) : od :

After some re-arrangements, the expressions for both coefficients assume more compact forms, viz.

> a[0] := ((2*t-1)*cos(t+2) - (t+2)*cos(2*t-1))/(t-3) : a[1] := numer(a[1])/denom(a[1]) :

> for i to No_roots(A(t)) do print(evaln(a[i-1]) = a[i-1]) : od :

Step 8 . Substitute the values of the coefficients into the matrix cos( [ A ( ) ] ), simplify the elements of the resultant matrix, and add a distinguishing subscript, say, C_H ( C ayley- H amilton) to its name for future comparative purposes:

> `cos(A(t))[C_H]` := map(simplify, map(x->eval(x), `cos(A(t))`)) :

> cos(A(t))[C_H] = matrix(`cos(A(t))[C_H]`) ;

> for i to No_roots(A(t)) do a[i-1] := evaln(a[i-1]) : od :

N.B. In general, the function cos( [ A ( ) ] ) is well defined for every value of t . This is true if the function elements of the matrix function are continuous functions for every t .

It can be easily noticed that the function element in the above matrix has a discontinuity at . Therefore, the function cos( [ A ( ) ] ) may be computed for only if a limiting value of the matrix element exists when t tends to 3 . Consequently, this must checked by computing the limit as follows:

> `limit(a[21](t))` := limit(`cos(A(t))[C_H]`[2,1], t=3) : Limit(a[21](t), t=3) = `limit(a[21](t))` ;

In computational practice, the above check may be omitted and an attempt made to compute the value of the function cos( [ A ( ) ] ) for directly as the limit of the matrix function when t tends to 3 . Thus,

> `lim cos(A(t))[C_H]` := map(limit, `cos(A(t))[C_H]`, t=3) :

> Limit(cos(A(t))[C_H], t=3) = matrix(`lim cos(A(t))[C_H]`) ;

Floating-point evaluation yields

> `lim cos(A(t))[C_H]` := evalf(matrix(`lim cos(A(t))[C_H]`)) :

> Limit(cos(A(t))[C_H], t=3) = matrix(`lim cos(A(t))[C_H]`) ;

[ For computation of limits of function matrices, refer to Unit (28). ]

Method 2 . Using similarity of matrices

Step 1 . Check if the eigenvalues of [ A ( ) ] are distinct:

See Step 4 above where the eigenvalues are found to be distinct.

Step 2 . Find the sequence of lists containing eigenvalues and eigenvectors of [ A ( ) ]:

> `roots&vectors(A(t))` := eigenvects(subs(A(t)=B(t), A(t))) :

> roots_and_vectors(A(t)) = `roots&vectors(A(t))` ;

Step 3 . Extract the eigenvectors corresponding to eigenvalues of [ A ( ) ]:

> for i to No_roots(A(t)) do e[i] := charroots(A(t))[i] : List[i] := `roots&vectors(A(t))`[i] : od :

> for j to No_roots(A(t)) do for i to No_roots(A(t)) do if List[i][1] = e[j] then Lst[j] := `roots&vectors(A(t))`[i] : fi : od : od :

> for i to No_roots(A(t)) do ch_v[i](A(t)) := op(Lst[i][3]) : print(v[lambda[i]] = eval(ch_v[i](A(t))), ` --> ` * lambda[i] = ch_r[i](A(t))) : od :

Conversion of the eigenvectors into column matrices yields

> for i to No_roots(A(t)) do X[l[i]] := convert(ch_v[i](A(t)), matrix) : print(X[l[i]] = matrix(X[l[i]])) : od :

Step 4 . Construct the 'unique' modal matrix [ M ( ) ] associated with [ A ( ) ]:

> i := 'i' : `seq(X[l[i]])` := seq(X[l[i]], i=1..coldim(B(t))) : M(t) := augment(`seq(X[l[i]])`) : 'M(t)'=M(t) ;

Step 5 . Construct the unique Jordan form [ J ( ) ] of [ A ( ) ]:

> J_A := array(diagonal, 1..rowdim(B(t)), 1..coldim(B(t))) :

> for i to rowdim(B(t)) do for j to coldim(B(t)) do if j = i then J_A[i,j] := ch_r[i](A(t)) fi : od : od :

> J(t) := matrix(J_A) : 'J(t)' = J(t) ;

Step 6 . Compute the function cosine of matrix [ J ( ) ]:

> J_A := array(diagonal, 1..rowdim(J(t)), 1..coldim(J(t))) :

> for i to rowdim(J(t)) do for j to coldim(J(t)) do if j = i then J_A[i,j] := cos(J(t)[i,j]) fi : od : od :

> `cos(J(t))` := matrix(J_A) : cos('J(t)') = matrix(`cos(J(t))`) ;

Step 7 . Compute cos( [ A ( ) ] ) = [ M ( ) ] cos( [ J ( ) ] ) Inv [ M ( ) ] , simplify the elements of the resultant matrix, and add a distinguishing subscript, say, s_m ( s imilarity of m atrices) to its name for future comparative purposes:

> `cos(A(t))[s_m]` := map(simplify, evalm(M(t) &* `cos(J(t))` &* M(t)^(-1))) :

> cos(A(t))[s_m] = matrix(`cos(A(t))[s_m]`) ;

The symbolic forms of matrices and are equal.

Method 3 . Using truncated Maclaurin’ s series

Step 1 . Find the eigenvalues of [ A ( ) ]:

This step was already done earlier to check if the eigenvalues of [ A ( ) ] are distinct, which was essential for both preceding methods. It is not necessary for this method since the series of is convergent for all , which implies that cos( [ A ( ) ] ) is well defined for every square matrix.

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

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

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

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

Step 4 . Evaluate the matrix cos( [ A ( ) ] ) with its symbolic elements after truncating the series to its first terms.

At this point, a relative disadvantage of this method appears in comparison with the previous methods: displaying of the matrix for a large number of terms of the Maclaurin series becomes impractical.

Exemplarily, adopt at this point only , simplify the element polynomials, and sort the terms of the matrix elements. (Note that Maple sorts the terms of polynomials into descending order.) These operations yield the matrix

> `cos(A(t))` := map(sort, map(simplify, evalm(subs(A(t)=B(t), U + sum((-1)^n*A(t)^(2*n)/(2*n)!, n=1..2))))) : cos(A(t)) = matrix(`cos(A(t))`) ;

However, computational experiments show that the number of series terms adopted, n , should be about to ensure a satisfactory accuracy. On the other hand, it is rather unlikely that a display of such matrices with their symbolic elements would be of any use. Practically, the result of numerical computation of the matrix function for a given value of t is required only.

Consequently, is adopted and the resultant matrix is added a distinguishing subscript, say, M_s ( M aclaurin’ s s eries ) to its name for future comparative purposes. Moreover, the functions simplify and sort are dropped, which results in a significant saving of computation time.

Thus, the final expression for cos( [ A ( ) ] ) computed by this method is adopted as follows:

> `cos(A(t))[M_s]`:= evalm(subs(A(t)=B(t), U + sum((-1)^n*A(t)^(2*n)/(2*n)!, n=1..50))) :

> cos(A(t))[M_s] = subs(x=B(t), matrix(U) + Sum((-1)^n*x^(2*n)/(2*n)!, n=1..50)) ; B(t) := 'B(t)' :

Comparison of numerical results of computation of cos( [ A ( ) ] )

Compute cos( [ A ( ) ] ) if using all the three methods.

> t := 4 :

(a) For Method 1 using the Cayley-Hamilton theorem:

> `cos(A(t))[C_H]` := map(x->eval(x), `cos(A(t))[C_H]`) : cos(A('t'))[C_H] = matrix(`cos(A(t))[C_H]`) ;

Floating-point evaluation yields

> `cos(A(t))[C_H]` := evalf(matrix(`cos(A(t))[C_H]`)) : cos(A('t'))[C_H] = matrix(`cos(A(t))[C_H]`) ;

(b) For Method 2 using similarity of matrices:

> `cos(A(t))[s_m]` := map(x->eval(x), `cos(A(t))[s_m]`) : cos(A('t'))[s_m] = matrix(`cos(A(t))[s_m]`) ;

Floating-point evaluation yields

> `cos(A(t))[s_m]` := evalf(matrix(`cos(A(t))[s_m]`)) : cos(A('t'))[s_m] = matrix(`cos(A(t))[s_m]`) ;

(c) For Method 3 using truncated Maclaurin’ s series:

> `cos(A(t))[M_s]` := evalf(map(x->eval(x), `cos(A(t))[M_s]`)) :

> cos(A('t'))[M_s] = matrix(`cos(A(t))[M_s]`) ; t := 't' :

The matrices of (a) and (b) are precisely equal and the matrix of (c) may be considered practically equal to the former ones.

* * *

Example 2

Compute sinh( [ A ( ) ] ) for a ( × ) matrix [ A ( ) ] given as

> A(t) := matrix(2, 2, [6*t, -5*t, 2*t, 4*t]) : 'A(t)' = A(t) ;

The power series of is convergent for all [refer to Unit (23)], which implies that sinh( [ A ] ) is well defined for every square matrix. This, in turn, means that the function sinh( [ A ( ) ] ) is well defined for every value of t , so it may be computed for this matrix.

For convenience in computations, let

> B(t) := A(t) : A(t) := 'A(t)' :

Method 1 . Using the Cayley-Hamilton theorem

Step 1 . Having in mind that for this matrix, write equation (1) of Unit (23) B for [ A ( ) ]:

> `sinh(A(t))` := a[1] * A(t) + a[0]*U : sinh(A(t)) = `sinh(A(t))` ; Sinh(A(t)) := `sinh(A(t))` :

Step 2 . With the unit matrix [ U ] as in Example 1, evaluate the matrix equation for sinh( [ A ( ) ] ):

> `sinh(A(t))` := evalm(subs(A(t)=B(t), `sinh(A(t))`)) : sinh(A(t)) = matrix(`sinh(A(t))`) ;

Step 3 . Formulate equation (2) of Unit (23) B corresponding to this case:

> `r(l)` := subs(A(t)=l, U=1, Sinh(A(t))) : r(l) = `r(l)` ;

Step 4 . Determine the eigenvalues of [ A ( ) ]:

> charroots(A(t)) := eigenvals(subs(A(t)=B(t), A(t))) : char_roots(A(t)) = charroots(A(t)) ;

Extracting the individual eigenvalues yields

> No_roots(A(t)) := nops([charroots(A(t))]) :

> for i to No_roots(A(t)) do ch_r[i](A(t)) := charroots(A(t))[i] : print(l[i] = ch_r[i](A(t))) : od :

Note that the eigenvalues are complex numbers.

Step 5 . Taking notice of the fact that the roots are distinct, formulate equation (3) of Unit (23) B in a general form:

Since the function corresponding to this case is

> `f(l)` := sinh(l) : f(l) = `f(l)` ;

therefore, equation (3) assumes the form

> Eq3_f(l) := `f(l)` = `r(l)` : Eq3_f(l) ;

Step 6 . Obtain a set of equations by substituting either eigenvalue into equation (3):

> for i to No_roots(A(t)) do Eq3_f[i](l) := subs(l=ch_r[i](A(t)), Eq3_f(l)) : print(Eq3_f[i](l)) : od :

Step 7 . Solve the simultaneous equations for the unknown coefficients:

> solution := solve({Eq3_f[1](l), Eq3_f[2](l)}, {a[0], a[1]}) : solution ;

Extracting either unknown from the solution set yields

> assign(solution) : for i to No_roots(A(t)) do a[i-1] := a[i-1] : print(evaln(a[i-1]) = a[i-1]) : od :

After some re-arrangements, the expressions for both coefficients are written in more compact forms, viz.

> a[0] := (1/2+5/6*I)*sinh((5+3*I)*t) + (1/2-5/6*I)*sinh((5-3*I)*t) : a[1] := numer(a[1])/denom(a[1]) :

> for i to No_roots(A(t)) do print(evaln(a[i-1]) = a[i-1]) : od :

Step 8 . Substitute the values of the coefficients into the matrix sinh( [ A ( ) ]:

> `sinh(A(t))` := map(x->eval(x), `sinh(A(t))`) : sinh(A(t)) = matrix(`sinh(A(t))`) ;

Each element of this matrix may be written in a more compact form after some extra manual re-arrangements, viz.

> `sinh(A(t))`[1,1] := (1/2-I/6)*sinh((5+3*I)*t) + (1/2+I/6)*sinh((5-3*I)*t) :

> `sinh(A(t))`[1,2] := 5*I/6*(sinh((5+3*I)*t) - sinh((5-3*I)*t)) :

> `sinh(A(t))`[2,1] := -I/3*(sinh((5+3*I)*t) - sinh((5-3*I)*t)) :

> `sinh(A(t))`[2,2] := (1/2+I/6)*sinh((5+3*I)*t) + (1/2-I/6)*sinh((5-3*I)*t) :

> for i to rowdim(`sinh(A(t))`) do for j to coldim(`sinh(A(t))`) do print(sinh(A(t))[i,j] = `sinh(A(t))`[i,j]) : od : od :

Substitute the above elements into the matrix sinh( [ A ( ) ] ) and add a distinguishing subscript, say, C_H ( C ayley- H amilton) to the name of the resultant matrix for future comparative purposes:

> `sinh(A(t))[C_H]` := matrix(`sinh(A(t))`) : sinh(A(t))[C_H] = matrix(`sinh(A(t))[C_H]`) ;

> for i to No_roots(A(t)) do a[i-1] := evaln(a[i-1]) : od :

Method 2 . Using similarity of matrices

Step 1 . Check if the eigenvalues of [ A ( ) ] are distinct:

See Step 4 above where the eigenvalues are found to be distinct.

Step 2 . Find the sequence of lists containing eigenvalues and eigenvectors of [ A ( ) ]:

> `roots&vectors(A(t))` := eigenvects(subs(A(t)=B(t), A(t))) :

> roots_and_vectors(A(t)) = `roots&vectors(A(t))` ;

Step 3 . Extract the eigenvectors corresponding to eigenvalues of [ A ( ) ]:

> for i to No_roots(A(t)) do e[i] := charroots(A(t))[i] : List[i] := `roots&vectors(A(t))`[i] : od :

> for j to No_roots(A(t)) do for i to No_roots(A(t)) do if List[i][1] = e[j] then Lst[j] := `roots&vectors(A(t))`[i] : fi : od : od :

> for i to No_roots(A(t)) do ch_v[i](A(t)) := op(Lst[i][3]) : print(v[lambda[i]] = eval(ch_v[i](A(t))), ` --> ` * lambda[i] = ch_r[i](A(t))) : od :

Conversion of the eigenvectors into column matrices yields

> for i to No_roots(A(t)) do X[l[i]] := convert(ch_v[i](A(t)), matrix) : print(X[l[i]] = matrix(X[l[i]])) : od :

Step 4 . Construct the 'unique' modal matrix [ M ( ) ] associated with [ A ( ) ]:

> i := 'i' : `seq(X[l[i]])` := seq(X[l[i]], i=1..coldim(B(t))) : M(t) := augment(`seq(X[l[i]])`) : 'M(t)'=M(t) ;

Step 5 . Construct the unique Jordan form [ J ( ) ] of [ A ( ) ]:

> J_A := array(diagonal, 1..rowdim(B(t)), 1..coldim(B(t))) :

> for i to rowdim(B(t)) do for j to coldim(B(t)) do if j = i then J_A[i,j] := ch_r[i](A(t)) fi : od : od :

> J(t) := matrix(J_A) : 'J(t)' = J(t) ;

Step 6 . Compute the function hyperbolic sine of matrix [ J ( ) ]:

> J_A := array(diagonal, 1..rowdim(J(t)), 1..coldim(J(t))) :

> for i to rowdim(J(t)) do for j to coldim(J(t)) do if j = i then J_A[i,j] := sinh(J(t)[i,j]) fi : od : od :

> `sinh(J(t))` := matrix(J_A) : sinh('J(t)') = matrix(`sinh(J(t))`) ;

Step 7 . Compute sinh( [ A ( ) ] ) = [ M ( ) ] sinh( [ J ( ) ] ) Inv [ M ( ) ] and simplify the elements of the resultant matrix:

> `sinh(A(t))`:=map(simplify, evalm(M(t) &* `sinh(J(t))` &* M(t)^(-1))) : sinh(A(t))=matrix(`sinh(A(t))`) ;

Each element of this matrix may be written in a more compact form, which appears to be identical to the corresponding element of the resultant matrix of Method 1. Thus, the above matrix may be rebuilt by replacing its original elements with the re-arranged elements of matrix . Consequently, the same matrix elements are used as in Method 1, which are assigned values in the following manner:

> for i to rowdim(`sinh(A(t))`) do for j to coldim(`sinh(A(t))`) do `sinh(A(t))`[i,j] := `sinh(A(t))[C_H]`[i,j] : od : od :

Substitute the above elements into the matrix sinh( [ A ( ) ] ) and add a distinguishing subscript, say, s_m ( s imilarity of m atrices) to the name of the resultant matrix for future comparative purposes:

> `sinh(A(t))[s_m]` := matrix(`sinh(A(t))`) : sinh(A(t))[s_m] = matrix(`sinh(A(t))[s_m]`) ;

Therefore, the symbolic forms of matrices and are equal.

Method 3 . Using truncated Maclaurin’ s series

Step 1 . Find the eigenvalues of [ A ( ) ]:

This step was already done earlier to check if the eigenvalues of [ A ( ) ] are distinct, which was essential for both preceding methods. It is not necessary for this method since the series of is convergent for all , which implies that sinh( [ A ( ) ] ) is well defined for every square matrix.

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

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

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

> `sinh(A(t))` := subs(A(t)=B(t), `sinh(A(t))`) : sinh(A(t)) = `sinh(A(t))` ;

Step 4 . Evaluate the matrix sinh( [ A ( ) ] ) with its symbolic elements after truncating the series to its first terms.

Since displaying of the matrix for a large number of terms of the Maclaurin series is impractical, only is adopted at this point and the matrix thus obtained is displayed for instructive purposes upon simplification and sorting of its elements:

> `sinh(A(t))`:=map(sort, map(simplify, evalm(subs(A(t)=B(t), sum(A(t)^(2*n+1)/(2*n+1)!, n=0..3))))) :

> sinh(A(t)) = matrix(`sinh(A(t))`) ;

In final computation, is adopted and the resultant matrix is added a distinguishing subscript, say, M_s ( M aclaurin’ s s eries ) to its name for future comparative purposes. Moreover, the functions simplify and sort are dropped, which results in a significant saving of computation time.

Thus, the final expression for sinh( [ A ( ) ] ) computed by this method is adopted as follows:

> `sinh(A(t))[M_s]`:= evalm(subs(A(t)=B(t), sum(A(t)^(2*n+1)/(2*n+1)!, n=0..50))) :

> sinh(A(t))[M_s] = subs(x=B(t), Sum(x^(2*n+1)/(2*n+1)!, n=0..50)) ; B(t) := 'B(t)' :

Comparison of numerical results of computation of sinh( [ A ( ) ] )

Compute sinh( [ A ( ) ] ) if using all the three methods.

> t := 0.5 :

(a) For Method 1 using the Cayley-Hamilton theorem:

> `sinh(A(t))[C_H]` := map(x->eval(x), `sinh(A(t))[C_H]`) : sinh(A('t'))[C_H] = matrix(`sinh(A(t))[C_H]`) ;

Disregarding the meaningless imaginary parts simplifies the above to

> `sinh(A(t))[C_H]` := evalm(Re(`sinh(A(t))[C_H]`)) : sinh(A('t'))[C_H] = matrix(`sinh(A(t))[C_H]`) ;

(b) For Method 2 using similarity of matrices:

> `sinh(A(t))[s_m]`:=map(x->eval(x), `sinh(A(t))[s_m]`) : sinh(A('t'))[s_m]=matrix(`sinh(A(t))[s_m]`) ;

Disregarding the meaningless imaginary parts simplifies the above to

> `sinh(A(t))[s_m]` := evalm(Re(`sinh(A(t))[s_m]`)) : sinh(A('t'))[s_m] = matrix(`sinh(A(t))[s_m]`) ;

(c) For Method 3 using truncated Maclaurin’ s series:

> `sinh(A(t))[M_s]` := evalf(map(x->eval(x), `sinh(A(t))[M_s]`)) :

> sinh(A('t'))[M_s] = matrix(`sinh(A(t))[M_s]`) ; t := 't' :

The matrices of (a) and (b) are precisely equal and the matrix of (c) may be considered practically equal to the former ones.

* * *

Example 3

Compute exp( [ A ( ) ] ) for a ( × ) matrix [ A ( ) ] given as

> A(t) := matrix(2, 2, [0, t, -t, 0]) : 'A(t)' = A(t) ;

The power series of is convergent for all [refer to Unit (23)], which implies that exp( [ A ] ) is well defined for every square matrix. This, in turn, means that the function exp( [ A ( ) ] ) is well defined for every value of t , so it may be computed for this matrix.

* * *

N.B. The above matrix [ A ( ) ] may be written equivalently as the product of a constant matrix [ A ] given as

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

and the scalar t , viz.

> A*t = matrix(A) * t ;

since the product [ A ] t evaluates to the matrix

> `At` := evalm(A * t) : A*t = matrix(`At`) ;

This easy observation has a deeper meaning when special applications of the function exp([ A ] t ) are examined in Unit (30).

> A := 'A' : A(t) := matrix(`At`) :

* * *

The exponential function of a square matrix may be computed in Maple directly , using the exponential function – refer also to Units (23) and (24). Therefore, computation of exp( [ A ( ) ] ) using this function is also included in this example as Method 4.

For convenience in computations, let

> B(t) := A(t) : A(t) := 'A(t)' :

Method 1 . Using the Cayley-Hamilton theorem

Step 1 . Having in mind that for this matrix, write equation (1) of Unit (23) B for [ A ( ) ]:

> `exp(A(t))` := a[1] * A(t) + a[0]*U : exp(A(t)) = `exp(A(t))` ; Exp(A(t)) := `exp(A(t))` :

Step 2 . With the unit matrix [ U ] as in Example 1, evaluate the matrix equation for exp( [ A ( ) ] ):

> `exp(A(t))` := evalm(subs(A(t)=B(t), `exp(A(t))`)) : exp(A(t)) = matrix(`exp(A(t))`) ;

Step 3 . Formulate equation (2) of Unit (23) B corresponding to this case:

> `r(l)` := subs(A(t)=l, U=1, Exp(A(t))) : r(l) = `r(l)` ;

Step 4 . Determine the eigenvalues of [ A ( ) ]:

> charroots(A(t)) := eigenvals(subs(A(t)=B(t), A(t))) : char_roots(A(t)) = charroots(A(t)) ;

Extracting the individual eigenvalues yields

> No_roots(A(t)) := nops([charroots(A(t))]) :

> for i to No_roots(A(t)) do ch_r[i](A(t)) := charroots(A(t))[i] : print(l[i] = ch_r[i](A(t))) : od :

Step 5 . Taking notice of the fact that the roots are distinct, formulate equation (3) of Unit (23) B in a general form:

Since the function corresponding to this case is

> `f(l)` := exp(l) : f(l) = `f(l)` ;

therefore, equation (3) assumes the form

> Eq3_f(l) := `f(l)` = `r(l)` : Eq3_f(l) ;

Step 6 . Obtain a set of equations by substituting either eigenvalue into equation (3):

> for i to No_roots(A(t)) do Eq3_f[i](l) := subs(l=ch_r[i](A(t)), Eq3_f(l)) : print(Eq3_f[i](l)) : od :

Step 7 . Solve the simultaneous equations for the unknown coefficients:

> solution := solve({Eq3_f[1](l), Eq3_f[2](l)}, {a[0], a[1]}) : solution ;

Extracting either unknown from the solution set yields

> assign(solution) : for i to No_roots(A(t)) do a[i-1] := a[i-1] : print(evaln(a[i-1]) = a[i-1]) : od :

Simplification of either coefficient gives

> for i to No_roots(A(t)) do a[i-1] := simplify(a[i-1]) : print(evaln(a[i-1]) = a[i-1]) : od :

Step 8 . Substitute the values of the coefficients into the matrix exp( [ A ( ) ] and add a distinguishing subscript, say, C_H ( C ayley- H amilton) to the name of the resultant matrix for future comparative purposes:

> `exp(A(t))[C_H]` := map(x->eval(x), `exp(A(t))`) : exp(A(t))[C_H] = matrix(`exp(A(t))[C_H]`) ;

> for i to No_roots(A(t)) do a[i-1] := evaln(a[i-1]) : od :

Method 2 . Using similarity of matrices

Step 1 . Check if the eigenvalues of [ A ( ) ] are distinct:

See Step 4 above where the eigenvalues are found to be distinct.

Step 2 . Find the sequence of lists containing eigenvalues and eigenvectors of [ A ( ) ]:

> `roots&vectors(A(t))` := eigenvects(subs(A(t)=B(t), A(t))) :

> roots_and_vectors(A(t)) = `roots&vectors(A(t))` ;

Step 3 . Extract the eigenvectors corresponding to eigenvalues of [ A ( ) ]:

> for i to No_roots(A(t)) do e[i] := charroots(A(t))[i] : List[i] := `roots&vectors(A(t))`[i] : od :

> for j to No_roots(A(t)) do for i to No_roots(A(t)) do if List[i][1] = e[j] then Lst[j] := `roots&vectors(A(t))`[i] : fi : od : od :

> for i to No_roots(A(t)) do ch_v[i](A(t)) := op(Lst[i][3]) : print(v[lambda[i]] = eval(ch_v[i](A(t))), ` --> ` * lambda[i] = ch_r[i](A(t))) : od :

Conversion of the eigenvectors into column matrices yields

> for i to No_roots(A(t)) do X[l[i]] := convert(ch_v[i](A(t)), matrix) : print(X[l[i]] = matrix(X[l[i]])) : od :

Step 4 . Construct the 'unique' modal matrix [ M ( ) ] associated with [ A ( ) ]:

> i := 'i' : `seq(X[l[i]])` := seq(X[l[i]], i=1..coldim(B(t))) : M(t) := augment(`seq(X[l[i]])`) : 'M(t)'=M(t) ;

Step 5 . Construct the unique Jordan form [ J ( ) ] of [ A ( ) ]:

> J_A := array(diagonal, 1..rowdim(B(t)), 1..coldim(B(t))) :

> for i to rowdim(B(t)) do for j to coldim(B(t)) do if j = i then J_A[i,j] := ch_r[i](A(t)) fi : od : od :

> J(t) := matrix(J_A) : 'J(t)' = J(t) ;

Step 6 . Compute the function exponential of matrix [ J ( ) ]:

> J_A := array(diagonal, 1..rowdim(J(t)), 1..coldim(J(t))) :

> for i to rowdim(J(t)) do for j to coldim(J(t)) do if j = i then J_A[i,j] := exp(J(t)[i,j]) fi : od : od :

> `exp(J(t))` := matrix(J_A) : exp('J(t)') = matrix(`exp(J(t))`) ;

Step 7 . Compute exp( [ A ( ) ] ) = [ M ( ) ] exp( [ J ( ) ] ) Inv [ M ( ) ] :

> `exp(A(t))`:= evalm(M(t) &* `exp(J(t))` &* M(t)^(-1)) : exp(A(t))=matrix(`exp(A(t))`) ;

Simplify the above matrix and add a distinguishing subscript, say, s_m ( s imilarity of m atrices) to the name of the resultant matrix for f