The Matrix Exponential via Eigenvalues and Eigenvectors 2.12 2002/06/21 2003/05/16 Steven Cox cox@rice.edu Liqun Wang liqun@rice.edu Steven Cox cox@rice.edu eigenvalue eigenvector exponential matrix This module introduces how to compute the matrix exponential using eigenvalues and eigenvectors. In this module we exploit the fact that the matrix exponential of a diagonal matrix is the diagonal matrix of element exponentials. In order to exploit it we need to recall that all matrices are almost diagonalizable. Let us begin with the clean case: if A is n-by-n and has n distinct eigenvalues, λj, and therefore n linear eigenvectors, sj, then we note that j j 1 … n A sj λj sj may be written A S Λ S where S s1 s2 … sn is the full matrix of eigenvectors and Λ diag λ1 λ2 … λn is the diagonal matrix of eigenvalues. One cool reason for writing A as in is that A 2 S Λ S S Λ S S Λ 2 S and, more generally A k S Λ k S If we now plug this into the definition in The Matrix Exponential as a Sum of Powers, we find A t S Λ t S where Λ t is simply diag λ1 t λ2 t … λn t Let us exercise this on our standard suite of examples. If A 1 0 0 2 then S I Λ A and so A t Λ t . This was too easy! As a second example let us suppose A 0 1 -1 0 and compute, in matlab, > [S, Lam] = eig(A) S = 0.7071 0.7071 0 + 0.7071i 0 - 0.7071i Lam = 0 + 1.0000i 0 0 0 - 1.0000i >> Si = inv(S) Si = 0.7071 0 - 0.7071i 0.7071 0 + 0.7071i >> simple(S*diag(exp(diag(Lam)*t))*Si) ans = [ cos(t), sin(t)] [-sin(t), cos(t)] ]]> If A 0 1 0 0 then matlab delivers > [S, Lam] = eig(A) S = 1.0000 -1.0000 0 0.0000 Lam = 0 0 0 0 ]]> So zero is a double eigenvalue with but one eigenvector. Hence S is not invertible and we can not invoke (). The generalization of () is often called the Jordan Canonical Form or the Spectral Representation. The latter reads A j 1 h λj Pj Dj where the λj are the distinct eigenvalues of A while, in terms of the resolvent R z z I A , Pj 1 2 z Cj R z is the associated eigen-projection and Dj 1 2 z Cj R z z λj is the associated eigen-nilpotent. In each case, Cj is a small circle enclosing only λj. Conversely we express the resolvent R z j 1 h 1 z λj Pj k 1 mj 1 1 z λj k 1 D j k where mj dim ℛ Pj with this preparation we recall Cauchy's integral formula for a smooth function f f a 1 2 z C a f z z a where C a is a curve enclosing the point a. The natural matrix analog is f A -1 2 z C r f z R z where C r encloses ALL of the eigenvalues of A. For f z z t we find A t j 1 h λj t Pj k 1 mj 1 t k k D j k with regard to our example we find, h 1 , λ 1 0 , P 1 I , m 1 2 , D 1 A so A t I t A Let us consider a slightly bigger example, if A 1 1 0 0 1 0 0 0 2 then > R = inv(s*eye(3)-A) R = [ 1/(s-1), 1/(s-1)^2, 0] [ 0, 1/(s-1), 0] [ 0, 0, 1/(s-2)] ]]> and so λ 1 1 and λ 2 2 while P1 1 0 0 0 1 0 0 0 0 and so m 1 2 D1 0 1 0 0 0 0 0 0 0 and P2 0 0 0 0 0 0 0 0 1 and m 2 1 and D 2 0 . Hence A t t P1 t D1 2 t P2 t t t 0 0 t 0 0 0 2 t