Matrix exponential
|
|
In mathematics, particularly in the study of Lie groups and Lie algebras, the matrix exponential is the function on square matrices A defined by
- <math>e^A=\exp(A)=\sum_{n=0}^\infty{A^n \over n!}.<math>
| Contents |
Motivational ideas
This concept has applications to systems of linear differential equations. Recall that a differential equation of the form
- y′ = Cy
has solution eCx. If we consider the vector
- <math> \mathbf{y}(x) = \begin{pmatrix} y_1(x) \\ \vdots \\y_n(x) \end{pmatrix}<math>
we can express a system of coupled linear differential equations as
- <math> \mathbf{y}'(x) = A\mathbf{y}(x)+\mathbf{b}<math>
If we make an ansatz and use an integrating factor of e−Ax and multiply throughout, we obtain
- <math>e^{-Ax}\mathbf{y}'(x)-e^{-Ax}A\mathbf{y} = e^{-Ax}\mathbf{b}<math>
- <math>D (e^{-Ax}\mathbf{y}) = e^{-Ax}\mathbf{b}<math>
If we can calculate eAx, then we can obtain the solution to the system.
Consider then the power series expansion of the exponential function:
- <math>\sum_{n=0}^\infty {x^n \over n!}<math>
If we replace x with A, then we have a means of calculating the matrix exponential. The issue of convergence arises, which will be dealt with in functions of matrices.
Properties
Two properties are important in the calculation of the exponential of a matrix via the Jordan form. Suppose J is the Jordan form of A, with P the transition matrix. Then
- <math>e^{A}=Pe^{J}P^{-1}.\,<math>
Also, since
- <math>J=J_{a_1}(\lambda_1)\oplus J_{a_2}(\lambda_2)\oplus\cdots\oplus J_{a_n}(\lambda_n),<math>
- <math>e^{J} = e^{J_{a_1}(\lambda_1)\oplus J_{a_2}(\lambda_2)\oplus\cdots\oplus J_{a_n}(\lambda_n)} = e^{J_{a_1}(\lambda_1)}\oplus e^{J_{a_2}(\lambda_2)}\oplus\cdots\oplus e^{J_{a_k}(\lambda_k)}.<math>
In some cases we have
- <math>e^{A+B}\neq e^A e^B<math>
but when A and B commute with each other, the familiar identity holds.
If A is a skew-symmetric matrix then eA is an orthogonal matrix.
If a matrix is diagonal
- <math>A=\begin{bmatrix} a_1 & 0 & \ldots & 0 \\
0 & a_2 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & a_n \end{bmatrix}, <math>
then its exponential can be obtained by just exponentiating every entry:
- <math>e^A=\begin{bmatrix} e^{a_1} & 0 & \ldots & 0 \\
0 & e^{a_2} & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & e^{a_n} \end{bmatrix}. <math>
This also allows to exponentiate diagonalizable matrices. If A = U−1DU and D is diagonal, then eA = U−1eDU.
Calculations
Consider the matrix
- <math>B=\begin{bmatrix}
21 & 17 & 6 \\ -5 & -1 & -6 \\ 4 & 4 & 16 \end{bmatrix}<math>
which has Jordan form
- <math>J=\begin{bmatrix}
16 & 1 & 0 \\ 0 & 16 & 0 \\ 0 & 0 & 4 \end{bmatrix}<math>
and transition matrix
- <math>P=\begin{bmatrix}
-1 & 1 & {5 \over 8} \\ 1 & -1 & -{1\over 8} \\ 0 & 2 & 0 \end{bmatrix}<math>
Now,
- <math>J=J_2(16)\oplus J_1(4)<math>
and
- <math>e^B = P e^{J} P^{-1} = P (e^{J_2(16)} \oplus e^{J_1(4)} ) P^{-1}.<math>
Recall that each Jordan block can be written as λI+N where N is nilpotent (suppose nilpotent with index k). The fact that N is nilpotent will aid us greatly:
- <math>e^{\lambda I+N}=e^{\lambda I}e^{N} = e^\lambda \left(I+N+{N^2 \over 2!}+\cdots+{N^{k-1} \over (k-1)!}+0+\cdots\right)<math>
which is a calculation with only finitely many terms!
So,
- <math>e^{16I+\begin{bmatrix} 0 & 1 \\ 0 & 0 \end{bmatrix}} = e^{16}\left(\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} + \begin{bmatrix} 0 & 1 \\ 0 & 0 \end{bmatrix} + {1 \over 2!}\begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix}+\cdots\right)=\begin{bmatrix} e^{16} & e^{16} \\ 0 & e^{16} \end{bmatrix}<math>
The exponential calculation for a 1×1 matrix is clearly trivial, with eJ1(4)=e4 so,
- <math>e^B = P\left(\begin{bmatrix} e^{16} & e^{16} \\ 0 & e^{16} \end{bmatrix} \oplus [e^4]\right)P^{-1} = P\begin{bmatrix} e^{16} & e^{16} & 0 \\ 0 & e^{16} & 0 \\ 0 & 0 & e^4 \end{bmatrix}P^{-1}=\begin{bmatrix}
-1 & 1 & {5 \over 8} \\ 1 & -1 & -{1\over 8} \\ 0 & 2 & 0 \end{bmatrix}\begin{bmatrix} e^{16} & e^{16} & 0 \\ 0 & e^{16} & 0 \\ 0 & 0 & e^4 \end{bmatrix}\begin{bmatrix} -1 & 1 & {5 \over 8} \\ 1 & -1 & -{1\over 8} \\ 0 & 2 & 0 \end{bmatrix}^{-1}<math>
- <math> e^B = {1\over 4}\begin{bmatrix}
5e^4-e^{16} & 5e^4 - 5 e^{16} & -2e^{16} \\ -e^4 + e^{16} & -e^4 + 5e^{16} & 2e^{16} \\ 0 & 0 & 4e^{16} \end{bmatrix}<math>
Clearly, to calculate the Jordan form and to evaluate the exponential this way is very tedious. Often, it will often suffice to calculate the action of the exponential matrix upon some vector in applications, and there are other methods available to achieve this.
Column method
If we have a matrix A with one eigenvalue λ, we can examine the generalized eigenspace of λ, and we can determine a constant k such that
- <math>GE_\lambda = (A-\lambda I)^k.\,<math>
However, observe
- <math>e^A\mathbf{v} = \exp\left(\lambda I+(A-\lambda I)\right)\mathbf{v}=e^\lambda\left(I+(A-\lambda I)+{1\over 2!}(A-\lambda I)^2+\cdots + {1 \over (k-1)!}(A-\lambda I)^{k-1}\right)\mathbf{v}<math>
also
- <math>e^{tA}\mathbf{v} = \exp\left(\lambda tI+t(A-\lambda I)\right)\mathbf{v}=e^{\lambda t}\left(I+t(A-\lambda I)+{t^2\over 2!}(A-\lambda I)^2+\cdots + {t^{k-1} \over (k-1)!}(A-\lambda I)^{k-1}\right)\mathbf{v}<math>
which, again, is a calculation with only finitely many terms.
If we firstly determine a basis of generalized eigenvectors for A, {v1,v2,...,vn}, we can form a matrix B with these vectors as columns. Using the above we can calculate the action of eA upon each of the vi, say, eAvi = wi, and likewise,
- <math>C=e^{A}B=e^{A}(\mathbf{v}_1|...|\mathbf{v}_n)=(\mathbf{w}_1|...|\mathbf{w}_n), e^{A}=CB^{-1}.<math>
We have
- <math>C=e^{tA}B=e^{tA}(\mathbf{v}_1|...|\mathbf{v}_n)=(\mathbf{w}_1|...|\mathbf{w}_n), e^{tA}=CB^{-1}.<math>
also, which will be of use later in applications.
Example
Consider the matrix
- <math>M=\begin{bmatrix}
2 & -1 & 1 \\ 0 & 3 & -1 \\ 2 & 1 & 3 \end{bmatrix}<math>
We need the eigenvectors and generalized eigenvectors for this matrix. From the characteristic polynomial of M, we get (x-2)2(x-4), and so we have eigenvalues 2, 2, and 4.
One eigenvector is v1=(-1, 1, 1)T corresponding to one of the eigenvalues 2. Since
- <math>(A-2I)\mathbf{v}_1 = \mathbf{0},\quad e^{tM}\mathbf{v}_1=e^{2t}(I)\mathbf{v}_1 = e^{2t}\begin{bmatrix} -1 \\ 1 \\ 1 \end{bmatrix}<math>
Another is v2=(1, -1, 1)T corresponding to eigenvalue 4. Again,
- <math>(A-4I)\mathbf{v}_2 = \mathbf{0},\quad e^{tM}\mathbf{v}_2=e^{4t}(I)\mathbf{v}_2 = e^{4t}\begin{bmatrix} 1 \\ -1 \\ 1 \end{bmatrix}<math>
The generalized eigenspace corresponding to eigenvalue 2,
- <math>GE_{2} = \ker{(A-2I)^2} = \mathrm{span}\ \left\{ \begin{bmatrix} -1 \\ 0 \\ 1 \end{bmatrix},\begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix} \right\}<math>
Choose the second vector, ie., v3=(0,1,0)T. Then
- <math>(A-2I)^2\mathbf{v}_1 = \mathbf{0},\quad e^{tM}\mathbf{v}_3 = e^{2t}(I+t(A-2I))\mathbf{v}<math>
- <math> = e^{2t}\left(\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}+t\begin{bmatrix} 0 & -1 & 1 \\ 0 & 1 & -1 \\ 2 & 1 & 1 \end{bmatrix}\right)\mathbf{v}=e^{2t}\left(\begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix}+t\begin{bmatrix} 0 & -1 & 1 \\ 0 & 1 & -1 \\ 2 & 1 & 1 \end{bmatrix}\begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix}\right)=e^{2t}\left(\begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix}+t\begin{bmatrix} -1 \\ 1 \\ 1 \end{bmatrix}\right)<math>
- <math> =e^{2t}\begin{bmatrix} -t \\ t + 1 \\ t \end{bmatrix}<math>
Now,
- <math> C = e^{tM}\left(\mathbf{v}_1\left|\mathbf{v}_2\right|\mathbf{v}_3\right) = e^{tM}\begin{bmatrix}
-1 & 1 & 0\\
1 & -1 & 1\\
1 & 1 & 0\end{bmatrix}=\begin{bmatrix}
-e^t & e^t & -te^{2t} \\
e^t & -e^t & (t+1)e^{2t} \\
e^t & e^t & te^{2t} \end{bmatrix}<math>
and
- <math>\begin{bmatrix}
-1 & 1 & 0\\
1 & -1 & 1\\
1 & 1 & 0\end{bmatrix}^{-1}
={1\over 2}\begin{bmatrix} -1 & 0 & 1 \\ 1 & 0 & 1 \\ 2 & 2 & 0 \end{bmatrix}<math>
so
- <math>
e^{tM}\begin{bmatrix} -1 & 1 & 0\\
1 & -1 & 1\\
1 & 1 & 0\end{bmatrix}=\begin{bmatrix}
-e^t & e^t & -te^{2t} \\
e^t & -e^t & (t+1)e^{2t} \\
e^t & e^t & te^{2t} \end{bmatrix}<math>
- <math>e^{tM}={1\over 2}\begin{bmatrix}
-e^t & e^t & -te^{2t} \\
e^t & -e^t & (t+1)e^{2t} \\
e^t & e^t & te^{2t} \end{bmatrix}\begin{bmatrix}
-1 & 0 & 1 \\
1 & 0 & 1 \\
2 & 2 & 0 \end{bmatrix}<math>
- <math> = \begin{bmatrix} 2e^t - 2te^{2t} & -2te^{2t} & 0 \\
-2e^t + 2(t+1)e^{2t} & 2(t+1)e^{2t} & 0 \\
2te^{2t} & 2te^{2t} & 2e^t\end{bmatrix}<math>
Applications
One application of the matrix exponential is its use in solving systems of differential equations. Say we have a system of differential equations
- <math>\begin{matrix}
y_1' &=& a_{11} y_1 &+& a_{21} y_2 &+& \cdots &+& a_{n1} y_n + b_1 \\ y_2' &=& a_{12} y_1 &+& a_{22} y_2 &+& \cdots &+& a_{n2} y_n + b_2 \\
&\vdots& & & & & & & \\
y_n' &=& a_{1n} y_1 &+& a_{2n} y_2 &+& \cdots &+& a_{nn} y_n + b_n \end{matrix}<math>
Clearly, we can represent this system in matrix form:
- <math>\begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}' =
\begin{bmatrix} a_{11} & a_{21} & \cdots & a_{n1} \\ a_{12} & a_{22} & \cdots & a_{n2} \\ \vdots & \vdots & \ddots & \vdots \\ a_{1n} & a_{2n} & \cdots & a_{nn} \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} +\begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{bmatrix}<math>
- <math>\mathbf{y}' = A\mathbf{y}+\mathbf{b}<math>
defining the derivative of a vector to be the derivative of its components. If b=0, the system is termed homogeneous, otherwise, inhomogeneous.
The general solution of the homogeneous system is given by a linear combination of the columns of exp(tA). Say xi is one of these columns, then
- <math>D(\mathbf{x}_i)=D\left(\exp{(tA)}\mathbf{e}_i\right)=D\left(\exp{(tA)}\right)\mathbf{e}_i = A\exp{(tA)}\mathbf{e}_i = A\mathbf{x}_i<math>
Example (homogeneous)
Say we have the system
- <math>\begin{matrix}
x' &=& 2x&-y&+z \\ y' &=& &3y&-1z \\ z' &=& 2x&+y&+3z \end{matrix}<math>
We have the associated matrix
- <math>M=\begin{bmatrix}
2 & -1 & 1 \\ 0 & 3 & -1 \\ 2 & 1 & 3 \end{bmatrix}<math>
In the example above, we have calculated the matrix exponential
- <math>e^{tM}=\begin{bmatrix}
2e^t - 2te^{2t} & -2te^{2t} & 0 \\
-2e^t + 2(t+1)e^{2t} & 2(t+1)e^{2t} & 0 \\
2te^{2t} & 2te^{2t} & 2e^t\end{bmatrix}<math>
so the general solution of the system is
- <math>\begin{bmatrix}x \\y \\ z\end{bmatrix}=
C_1\begin{bmatrix}2e^t - 2te^{2t} \\-2e^t + 2(t+1)e^{2t}\\2te^{2t}\end{bmatrix} +C_2\begin{bmatrix}-2te^{2t}\\2(t+1)e^{2t}\\2te^{2t}\end{bmatrix} +C_3\begin{bmatrix}0\\0\\2e^t\end{bmatrix}<math>
that is,
- <math>\begin{matrix}
x &=& C_1(2e^t - 2te^{2t}) + C_2(-2te^{2t})\\ y &=& C_1(-2e^t + 2(t+1)e^{2t})+C_2(2(t+1)e^{2t})\\ z &=& (C_1+C_2)(2te^{2t})+2C_3e^t\end{matrix}<math>
Inhomogeneous case - variation of parameters
For the inhomogeneous case, we can use a method akin to variation of parameters. We seek a particular solution of the form yp(t)=exp(tA)z(t) :
- <math>\mathbf{y}_p' = (e^{tA})'\mathbf{z}(t)+e^{tA}\mathbf{z}'(t)<math>
- <math>= Ae^{tA}\mathbf{z}(t)+e^{tA}\mathbf{z}'(t)<math>
- <math>= A\mathbf{y}_p(t)+e^{tA}\mathbf{z}'(t)<math>
For yp to be a solution:
- <math>e^{tA}\mathbf{z}'(t) = \mathbf{b}(t)<math>
- <math>\mathbf{z}'(t) = (e^{tA})^{-1}\mathbf{b}(t)<math>
- <math>\mathbf{z}(t) = \int_0^t e^{-uA}\mathbf{b}(u)\,du+\mathbf{c}<math>
So,
- <math>\mathbf{y}_p = e^{tA}\int_0^t e^{-uA}\mathbf{b}(u)\,du+e^{tA}\mathbf{c}<math>
- <math>= \int_0^t e^{(t-u)A}\mathbf{b}(u)\,du+e^{tA}\mathbf{c}<math>
where c is determined by the initial conditions of the problem.
Example (inhomogeneous)
Say we have the system
- <math>\begin{matrix}
x' &=& 2x&-y&+z&+e^{2t} \\ y' &=& &3y&-1z& \\ z' &=& 2x&+y&+3z&+e^{2t} \end{matrix}<math>
So we then have
- <math>M=\begin{bmatrix}
2 & -1 & 1 \\ 0 & 3 & -1 \\ 2 & 1 & 3 \end{bmatrix}<math>
and
- <math>\mathbf{b}=e^{2t}\begin{bmatrix}1 \\0\\1\end{bmatrix}<math>
From before, we have the general solution to the homogeneous equation, Since the sum of the homogeneous and particular solutions give the general solution to the inhomogeneous problem, now we only need to find the particular solution (via variation of parameters).
We have, above:
- <math>\mathbf{y}_p = e^{t}\int_0^t e^{(-u)A}\begin{bmatrix}e^{2u} \\0\\e^{2u}\end{bmatrix}\,du+e^{tA}\mathbf{c}<math>
- <math>\mathbf{y}_p = e^{t}\int_0^t
\begin{bmatrix}
2e^u - 2ue^{2u} & -2ue^{2u} & 0 \\
-2e^u + 2(u+1)e^{2u} & 2(u+1)e^{2u} & 0 \\
2ue^{2u} & 2ue^{2u} & 2e^u\end{bmatrix}\begin{bmatrix}e^{2u} \\0\\e^{2u}\end{bmatrix}\,du+e^{tA}\mathbf{c}<math>
- <math>\mathbf{y}_p = e^{t}\int_0^t
\begin{bmatrix} e^{2u}( 2e^u - 2ue^{2u}) \\
e^{2u}(-2e^u + 2(1 + u)e^{2u}) \\
2e^{3u} + 2ue^{4u}\end{bmatrix}+e^{tA}\mathbf{c}<math>
- <math>\mathbf{y}_p = e^{t}\begin{bmatrix}
-{1 \over 24}e^{3t}(3e^t(4t-1)-16) \\ {1 \over 24}e^{3t}(3e^t(4t+4)-16) \\ {1 \over 24}e^{3t}(3e^t(4t-1)-16)\end{bmatrix}+ \begin{bmatrix}
2e^t - 2te^{2t} & -2te^{2t} & 0 \\
-2e^t + 2(t+1)e^{2t} & 2(t+1)e^{2t} & 0 \\
2te^{2t} & 2te^{2t} & 2e^t\end{bmatrix}\begin{bmatrix}c_1 \\c_2 \\c_3\end{bmatrix}<math>
which can be further simplified to get the requisite particular solution determined through variation of parameters.