Matrices-Unit30.html

Matrices-Unit30.mws

MATRICES AND MATRIX OPERATIONS: Unit 30

Dr. Wlodzislaw Kostecki

The Papua New Guinea University of Technology (PNGUT)

Department of Electrical and Communication Engineering

Lae, Morobe Province

Papua New Guinea

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

(30) Application of the function exp([ A ] t ) in solving matrix equations

OBJECTIVES :

• To introduce the function exp([ A ] t ) for computing the exponential function of the product of a constant square matrix and a real variable t .

• To exemplify how this operation may be performed in alternative ways.

• To examine essential properties of matrix multiplication where matrices of the type exp([ A ] t ) are involved.

• To provide a brief introduction to, and necessary formula for, solving matrix equations.

• To state a condition necessary for the solution to a matrix equation to be unique.

• To provide a step-by-step example of solving a matrix equation and verify the solution.

• To provide a brief introduction to, and necessary formulae for, solving matrix differential equations.

• To provide a step-by-step example of solving a matrix differential equation and verify the solution.

• To provide a universal method for clearing fractions from matrix elements with Maple .

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

The function exp([ A ] t ) for a constant square matrix [ A ] and real variable t is of special importance.

The function was used in the form exp( [ A ( ) ] ) in Unit (24) where elements of a square matrix [ A ( ) ] contained only the proportional terms of linear function , for computing the exponential function of [ A ( ) ].

However, the major application of the function exp([ A ] t ) is in solving matrix equations and matrix differential equations.

The matrix function exp([ A ] t ) may be computed by setting [ B ] = [ A ] t and then calculating exp([ B ]), or exp([ A ] t ) may be computed directly, using the exponential function.

A few examples of computation of the matrix function exp([ A ] t ) for different constant matrices [ A ] follow.

Example 1

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

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

> B := evalm(A * t) : B = matrix(B) ;

The matrix function exp([ A ] t ) is the following matrix:

> `exp(B)` := exponential(B) : exp(A*t) = matrix(`exp(B)`) ;

The same result may be obtained in Maple directly using the exponential function and scalar multiplication [ A ] t , viz.

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

or using the exponential function and t as a scalar parameter, viz.

> `exp(At)` := exponential(A, t) : exp(A*t) = matrix(`exp(At)`) ;

* * *

N.B. The eigenvalues of [ B ] = [ A ] t are the eigenvalues of [ A ] multiplied by t .

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

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

* * *

Example 2

Compute exp([ U ] t ) for the ( × ) unit matrix [ U ]

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

> `exp(Ut)` := exponential(U, t) : exp(U*t) = matrix(`exp(Ut)`) ;

Example 3

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

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

> `exp(At)` := exponential(A, t) : exp(A*t) = matrix(`exp(At)`) ;

Example 4

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

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

> `exp(At)` := exponential(A, t) : exp(A*t) = matrix(`exp(At)`) ;

Example 5

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

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

> `exp(At)` := exponential(A, t) : exp(A*t) = matrix(`exp(At)`) ;

The above matrix may be displayed in a more compact form, viz.

> `exp(At)` := map(collect, `exp(At)`, exp(t)) : exp(A*t) = matrix(`exp(At)`) ;

* * *

N.B. Computation of the the matrix exponential of [ A ] t depends on finding eigenvalues for [ A ] t . A symbolic result can be achieved only if symbolic eigenvalues can be found.

* * *

N.B. If two matrices, [ A ] and [ B ], satisfy the commutative law of multiplication, then

exp([ A ] t ) exp([ B ] t ) = exp[( [ A ] + [ B ] ) t ]

Exemplarily, consider the two ( × ) commuting matrices [ A ] and [ B ] of Unit (4) given as

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

(a) The product exp([ A ] t ) exp([ B ] t ) is the following ( × ) matrix:

> `exp(At) + exp(Bt)` := multiply(exponential(A, t), exponential(B, t)) :

> `exp(At) + exp(Bt)` := map(simplify, `exp(At) + exp(Bt)`) :

> exp(A*t) * exp(B*t) = matrix(`exp(At) + exp(Bt)`) ;

(b) The exponential exp[( [ A ] + [ B ] ) t ] is the following ( × ) matrix:

> `exp(A + Bt)` := exponential((evalm(A + B)), t) : exp((A+B)*t) = matrix(`exp(A + Bt)`) ;

Both matrices of (a) and (b) are equal.

* * *

N.B. For any constant square matrix [ A ] and real variables t and s ,

exp([ A ] t ) exp(–[ A ] s ) = exp{[ A ] ( t – s )]}

Exemplarily, consider the above matrix [ A ].

(a) The product exp([ A ] t ) exp(–[ A ] s ) is the following matrix:

> `exp(At) exp(-As)` := multiply(exponential(A, t), exponential(-A, s)) :

> `exp(At) exp(-As)` := map(combine, map(expand, `exp(At) exp(-As)`)) :

> exp(A*t) * exp(-A*s) = matrix(`exp(At) exp(-As)`) ;

(b) The exponential exp{[ A ] ( t – s )]} is the following matrix:

> `exp(A(t-s))` := exponential(evalm(A*(t-s))) :

> `exp(A(t-s))` := map(combine, map(expand, `exp(A(t-s))`)) :

> exp(A*(t-s)) = matrix(`exp(A(t-s))`) ;

The equal function applied to the resultant matrices of (a) and (b), i.e.

> equal(`exp(At) exp(-As)`, `exp(A(t-s))`) ;

returns the Boolean value , which verifies that both matrices are equal.

* * *

The two major application domains of the function exp([ A ] t ) are illustrated hereunder in Sections A and B .

A. Solution of matrix equations

The equation

[ A ] [ X ] + [ X ] [ B ] = [ C ]

where [ A ], [ B ], and [ C ] denote constant square matrices of the same order, has a unique solution [ X ] if and only if [ A ] and [ B ] have no eigenvalues in common. This unique solution is given by

[ X ] = –

provided the improper integral exists.

Exemplarily, solve the matrix equation [ A ] [ X ] + [ X ] [ B ] = [ C ] for [ X ] if ( × ) matrices [ A ], [ B ], and [ C ] are given as

> A := matrix(2, 2, [-6, 0, 0, -8]) : B := matrix(2, 2, [3, 2, 4, 1]) : C := matrix(2, 2, [2, 0, -1, 0]) :

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

Step 1 . Find the eigenvalues of matrices [ A ] and [ B ]:

• for matrix [ A ]

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

• for matrix [ B ]

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

The two matrices have no eigenvalues in common. Thus, the equation has a unique solution [ X ].

Step 2 . Determine the exponential functions of [ A ] and [ B ] in variable t :

• for matrix [ A ]

> `exp(At)` := exponential(A, t) : exp(A*t) = matrix(`exp(At)`) ;

• for matrix [ B ]

> `exp(Bt)` := map(sort, exponential(B, t)) : exp(B*t) = matrix(`exp(Bt)`) ;

Step 3 . Determine the matrix integrand exp([ A ] t ) [ C ] exp([ B ] t ):

> M(t) := evalm(`exp(At)` &* C &* `exp(Bt)`) : exp(A*t) * C * exp(B*t) = M(t) ;

After some mathematical operations, the above matrix simplifies to the form

> M(t) := map(combine, map(expand, M(t), power), exp) : exp(A*t) * C * exp(B*t) = M(t) ;

Step 4 . Compute the integral :

> `int(M(t))` := evalm(map(-int, M(t), t=0..infinity)) : -Int(exp(A*t) * C * exp(B*t), t=0..infinity) = matrix(`int(M(t))`) ;

Step 5 . Obtain the solution matrix [ X ]:

> X := matrix(`int(M(t))`) : X = matrix(X) ;

Clearing of fractions may be applied to the elements of this matrix by multiplying them by the least common multiple of their denominators and dividing the new matrix by the same. If this to be done manually, it is intuitive and simple. If the task is to be performed with Maple , the procedure is quite elaborate and requires several operations as shown.

(a) Create a list of the elements of [ X ]:

> L := convert(convert(X, set), list) : elements(X) = L ;

(b) Extract denominators of the elements of the list:

> for i to nops(L) do den(el[i]('L')) := denom(op(i, L)) : print(evaln(denominator(element[i]('L'))) = den(el[i]('L'))) : od :

> `seq(den(el(L)))` := [seq(den(el[n]('L')), n=1..nops(L))] : den(el('L')) := op(`seq(den(el(L)))`) : denominators(elements(X)) = den(el('L')) ;

(d) Find the least common multiple of denominators of the elements of the list:

> LCM := lcm(den(el('L'))) : least_common_multiple(denominators(X)) = LCM ;

(e) Multiply the elements of [ X ] by the least common multiple of their denominators and divide the new matrix by the same:

> X := evalm(X*LCM)/LCM : 'X' = X ;

Step 6 . Verify the solution by computing the sum matrix [ A ] [ X ] + [ X ] [ B ]:

> X := evalm(X) :

> `AX + XB` := evalm(A &* X + X &* B) : A*X + X*B = matrix(`AX + XB`) ; X := 'X' :

The sum matrix [ A ] [ X ] + [ X ] [ B ] is equal to matrix [ C ], which verifies the solution.

* * *

B. Solution of matrix differential equations

The initial-value matrix differential inhomogeneous equation

= [ A ][ X ( t ) ] + [ F ( t ) ]

with initial condition

[ X ( ) ] = [ C ]

has the solution

[ X ( t ) ] = [ C ] + (1)

or, equivalently,

[ X ( t ) ] = [ C ] + (2)

NOTE 1: If the differential equation is homogeneous, i.e.,

[ F ( t ) ] = [ 0 ]

then formulae (1) and (2) reduce to

[ X ( t ) ] = [ C ]

NOTE 2: Obtaining [ X ( t ) ] from formula (2) requires one operation step less than using formula (1). However, the integrals arising in formula (2) are generally more difficult to evaluate than those in formula (1).

As an example, use formula (1) to find the solution [ X ( t ) ] to the differential equation = [ A ][ X ( t ) ] + [ F ( t ) ] with initial value [ X ( ) ] = [ C ] if , when matrices [ A ], [ F ( t ) ], and [ C ] are given as

> A := matrix(2, 2, [0, 1, 8, -2]) : F(t) := matrix(2, 1, [0, exp(t)]) : C := matrix(2, 1, [1, -4]) :

> A = matrix(A) ; 'F(t)' = F(t) ; C = matrix(C) ;

Step 1 . Determine the exponential function of [ A ] in variable :

> `exp(A(t-to))` := subs(t[0]=0, evalm(exponential(A, t) &* exponential(-A, t[0]))) :

> exp(A*(t-t[0])) = matrix(`exp(A(t-to))`) ;

After some mathematical operations, the above matrix simplifies to the form

> `exp(A(t-to))` := map(sort, map(combine, map(expand, `exp(A(t-to))`, power), exp)) :

> exp(A*(t-t[0])) = matrix(`exp(A(t-to))`) ;

Step 2 . Determine the product matrix [ C ]:

> `exp(A(t-to)) C` := multiply(`exp(A(t-to))`, C) : exp(A*(t-t[0])) * C = matrix(`exp(A(t-to)) C`) ;

Step 3 . Determine the exponential function of [ A ] in variable :

> `exp(At)` := map(sort, exponential(A, t)) : exp(A*t) = matrix(`exp(At)`) ;

Step 4 . Determine the exponential function of – [ A ] in variable :

> `exp(-As)` := map(sort, exponential(-A, s)) : exp(-A*s) = matrix(`exp(-As)`) ;

Step 5 . Determine the matrix integrand [ F ( ) ]:

> M(s) := evalm(`exp(-As)` &* subs(t=s, F(t))) : exp(-A*s) * F(s) = M(s) ;

After some mathematical operations, the above matrix simplifies to the form

> M(s) := map(combine, map(expand, M(s), power), exp) : exp(-A*s) * F(s) = M(s) ;

Step 6 . Compute the definite integral :

> `int(M(s))` := map(sort, evalm(subs(t[0]=0, map(int, M(s), s=t[0]..t)))) :

> Int(exp(-A*s)*F(s), s=t[0]..t) = matrix(`int(M(s))`) ;

After some mathematical operations, the above matrix simplifies to the form

> `int(M(s))` := map(combine, map(expand, `int(M(s))`)) :

> Int(exp(-A*s)*F(s), s=t[0]..t) = matrix(`int(M(s))`) ;

Step 7 . Determine the product matrix :

> `exp(At) int(M(s))` := evalm(`exp(At)` &* `int(M(s))`) :

> exp(A*t)*Int(exp(-A*s)*F(s), s=t[0]..t) = matrix(`exp(At) int(M(s))`) ;

After some mathematical operations, the above product matrix simplifies to the form

> `exp(At) int(M(s))` := map(sort, map(combine, map(expand, `exp(At) int(M(s))`))) :

> exp(A*t)*Int(exp(-A*s)*F(s), s=t[0]..t) = matrix(`exp(At) int(M(s))`) ;

Step 8 . Compute the sum matrix :

> `exp(A(t-to)) C + exp(At) int(M(s))` := evalm(`exp(A(t-to)) C` + `exp(At) int(M(s))`) :

> exp(A*(t-t[0])) * C + exp(A*t)*Int(exp(-A*s)*F(s), s=t[0]..t) = matrix(`exp(A(t-to)) C + exp(At) int(M(s))`) ;

Step 9 . Obtain the solution matrix [ X ( t ) ]:

> X(t) := evalm(`exp(A(t-to)) C + exp(At) int(M(s))`) : 'X(t)' = X(t) ;

Clearing of fractions may be applied to the elements of this matrix by multiplying them by the least common multiple of denominators of their numeric operands and dividing the new matrix by the same. If this to be done manually, it is intuitive and simple. If the task is to be performed with Maple , the procedure is quite elaborate and requires several operations. A suitable method proposed below has been developed with a view to providing a universal approach to tackling more complicated cases. The method also works when some of the numeric operands of matrix elements are natural numbers.

(a) Create a list of the elements of [ X ( t ) ]:

> L := convert(convert(X(t), set), list) : elements('X(t)') = L ;

(b) Extract denominators of the numeric operands of the elements of the list:

> for i to nops(L) do for n to nops(op(i, L)) do den(el[i,n]('L')) := denom(op(n, L[i])) : print(evaln(denominator(element[i,n]('L'))) = den(el[i,n]('L'))) : od : od :

> k := 0 : for i to nops(L) do for n to nops(op(i, L)) do k := k+1 : den(el[k]('L')) := den(el[i,n]('L')) : print(evaln(denominator(element[k]('L'))) = den(el[k]('L'))) : od : od :

(d) Count the total number of denominators of the numeric operands of the elements of the list:

> num_ops('L') := k : 'num_ops(L)' = num_ops('L') ;

(e) Create a sequence of denominators of the numeric operands of the elements of the list:

> `seq(den(el(L)))` := [seq(den(el[m]('L')), m=1..num_ops('L'))] : den(el('L')) := op(`seq(den(el(L)))`) : denominators(elements('X(t)')) = den(el('L')) ;

(f) Find the least common multiple of denominators of the numeric operands of the elements of [ X ( t ) ]:

> LCM := lcm(den(el('L'))) : least_common_multiple(denominators('X(t)')) = LCM ;

(g) Multiply the elements of [ X ( t ) ] by the least common multiple of denominators of their numeric operands and divide the new matrix by the same:

> X(t) := evalm(X(t) * LCM)/LCM : 'X(t)' = X(t) ;