Section 5.3

Two Applications

>

Warning, the protected names norm and trace have been redefined and unprotected

>

An Application: Exponential of a Matrix

Exponential Form of a Matrix

The differential equation dy(t)/dt = y'(t) = a*y(t) is a prototype equation that models

problems of growth and decay. For example, if y(t) is the count of bacteria at time t and

if the rate of change of bacteria is proportional to the count at time t, then the differential

equation models this phenomenon with a being the constant of proportionality.

The solution is y(t) = c*exp(a*t) and c is an arbitrary constant that can be determined

using the initial condition y(0) .

The exponential function is given by the expansion

> exp(a*t):=Sum((a*t)^i/i!,i=0..infinity);

exp(a*t) := Sum((a*t)^i/i!,i = 0 .. infinity)

>

Since this function is a solution of many important problems including the problem of

exponential growth or decay, a natural question to ask in this setting:

is there an equivalent expansion where the parameter a is replaced by a matrix A?

If so, how do we compute this expression?

This result is needed if we were to replace the scalar equation by the matrix equation

Y'(t) = A* Y(t),

where Y(t) = [ y[1](t), y[2](t) , ...., y[n](t) ] and solve the resulting system. The solution of this system is Y(t) = exp(A*t)*Y(0) where Y(0) is the initial condition.

To have an analogous exponential function for matrices

> exp(A*t):=Sum( (A*t)^i/i!,i=0..infinity);

exp(A*t) := Sum((A*t)^i/i!,i = 0 .. infinity)

we need to compute the powers of the matrix A for different type of matrices.

Case of a Diagonal Matrix

Let A be the 2x2 identity matrix

> A:=diag(1,1);

A := matrix([[1, 0], [0, 1]])

Since all powers of A are equal to A

> exp(A*t):=evalm(A)*Sum( (t)^i/i!, i =0..infinity);

exp(A*t) := matrix([[1, 0], [0, 1]])*Sum(t^i/i!,i =...

This is equivalent to

> exp(A*t):=evalm((A)*sum( t^i/i!, i =0..infinity));

exp(A*t) := matrix([[exp(t), 0], [0, exp(t)]])

>

Let A be the 2x2 diagonal matrix

> A:=diag(2,3);

A := matrix([[2, 0], [0, 3]])

Since A is a diagonal matrix, its powers are easy to compute. In particular, for the matrix A

> A^2=evalm(A^2);
A^3=evalm(A^3);
A^4=evalm(A^4);

A^2 = matrix([[4, 0], [0, 9]])

A^3 = matrix([[8, 0], [0, 27]])

A^4 = matrix([[16, 0], [0, 81]])

In general,

> A^i = matrix([[2^i,0],[0,3^i]]);

A^i = matrix([[2^i, 0], [0, 3^i]])

Therefore

> exp(A*t):=Sum(matrix([[2^i,0],[0,3^i]])*t^i/i!,
i=0..infinity);

exp(A*t) := Sum(matrix([[2^i, 0], [0, 3^i]])*t^i/i!...

That is

> exp(A*t):=Sum( matrix([[(2*t)^i/i!,0],[0,(3*t)^i/i!]]),
i =0..infinity);

exp(A*t) := Sum(matrix([[(2*t)^i/i!, 0], [0, (3*t)^...

or

> exp(A*t):= matrix([[Sum((2*t)^i/i!,
i=0..infinity),0],[0,Sum((3*t)^i/i!,
i=0..infinity)]]);

exp(A*t) := matrix([[Sum((2*t)^i/i!,i = 0 .. infini...

When the sums are evaluated , we get the matrix

> exp(A*t):= matrix([[sum((2*t)^i/i!,
i=0..infinity),0],
[0,sum((3*t)^i/i!,
i=0..infinity)]]);

exp(A*t) := matrix([[exp(2*t), 0], [0, exp(3*t)]])

>

>

>

Case of a Diagonalizable Matrix

Now consider the 2x2 matrix

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

A := matrix([[3, 4], [3, 2]])

>

The eigenvalues of the matrix A are: 6 and -1 and the corresponding eigenvectors are

> v1:=vector([4,3]); v2:=vector([1,-1]);

v1 := vector([4, 3])

v2 := vector([1, -1])

>

Thus the matrix A is similar to a diagonal matrix A = P*D[1]*P^`-1` with

> P:=matrix([[4,1],[3,-1]]):D1:=diag(6,-1): evalm(A)=evalm(P)*evalm(D1)*evalm(inverse(P));

matrix([[3, 4], [3, 2]]) = matrix([[4, 1], [3, -1]]...

>

Let us use the fact that A is diagonalizable to compute various powers of A:

> A^2=evalm(A^2);P*D1^2*P^(`-1`)=
multiply(P,multiply(D1^2,inverse(P)));

A^2 = matrix([[21, 20], [15, 16]])

P*D1^2*P^`-1` = matrix([[21, 20], [15, 16]])

> A^3=evalm(A^3);P*D1^3*P^(`-1`)=
multiply(P,multiply(D1^3,inverse(P)));

A^3 = matrix([[123, 124], [93, 92]])

P*D1^3*P^`-1` = matrix([[123, 124], [93, 92]])

> A^4=evalm(A^4);P*D1^4*P^(`-1`)=
multiply(P,multiply(D1^4,inverse(P)));

A^4 = matrix([[741, 740], [555, 556]])

P*D1^4*P^`-1` = matrix([[741, 740], [555, 556]])

>

Then

> A^n=P*D1^n*P^(`-1`);

A^n = P*D1^n*P^`-1`

>

Furthermore, the exponential of the matrix D[1] is:

> exp(D1*t):=diag(exp(6*t),exp(-t));

exp(D1*t) := matrix([[exp(6*t), 0], [0, exp(-t)]])

> exp(A*t):=evalm(P&*exp(D1*t)&*inverse(P));

exp(A*t) := matrix([[4/7*exp(6*t)+3/7*exp(-t), 4/7*...

>

******************************************************************

To get the exponential of a diagonalizable matrix A proceed as follows :

1. compute the eigenvalues and eigenvectors of the given matrix

2. construct the matrix P whose columns are the eigenvectors and form the product

A = P*D[1]*P^`-1` where D[1] is the diagonal matrix whose diagonal entries are the eigenvalues of A.

3. construct the exponential matrix exp(D[1]*t)

4. the exponential of the given matrix is P*exp(D1*t)*P^`-1` .

******************************************************************

Initial value problem

Solve the initial value problem

> A:=matrix([[3,2],[0,4]]): Y:=matrix([[y1],[y2]]):

> evalm((d/dt)*Y) = evalm(A)*evalm(Y);

matrix([[d/dt*y1], [d/dt*y2]]) = matrix([[3, 2], [0...

with initial data Y(0)

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

Y0 := matrix([[1], [2]])

1. Compute the eigenvalues of A (use the nostep mode)

> eigenvals(A);

3, 4

>

3

>

2. Compute the eigenvectors of matrix A (use the nostep mode)

> eigenvects(A);

[3, 1, {vector([1, 0])}], [4, 1, {vector([2, 1])}]

>

>

An eigenvector corresponding to the eigenvalue 3 is the vector [1,0] and an eigenvector corresponding to the eigenvalue 4 is the vector [2,1].

3.Construct the matrix P and the diagonal matrix D[1] .

> P:=matrix([[1,2],[0,1]]); D1:=diag(3,4);

P := matrix([[1, 2], [0, 1]])

D1 := matrix([[3, 0], [0, 4]])

>

4. construct the exponential of D[1] :

> exp(D1*t):=diag(exp(3*t),exp(4*t));

exp(D1*t) := matrix([[exp(3*t), 0], [0, exp(4*t)]])...

>

The exponential matrix of A is

> exp(A*t):=evalm(P&*exp(D1*t)&*inverse(P));

exp(A*t) := matrix([[exp(3*t), -2*exp(3*t)+2*exp(4*...

Therefore the solution to the initial value problem is

> Y:=evalm(exp(A*t)&*Y0);

Y := matrix([[-3*exp(3*t)+4*exp(4*t)], [2*exp(4*t)]...

>

Exercises

6, 7, 8.