3.8 Matrix exponentials

Definition

In this section we present a different way of finding the fundamental matrix solution of a system. Suppose that we have the constant coefficient equation

$${\overrightarrow{x}}^{\prime }=P\overrightarrow{x}$$

as usual. Now suppose that this was one equation ($P$ is a number or a $1\times 1$ matrix). Then the solution to this would be

$$\overrightarrow{x}=e^{Pt}.$$

That doesn’t make sense if $P$ is a larger matrix, but essentially the same computation that led to the above works for matrices when we define $e^{Pt}$ properly. First let us write down the Taylor series for $e^{at}$for some number $a$.

$$e^{at}=1+at+\frac{(at{)}^2}{2}+\frac{(at{)}^3}{6}+\frac{(at{)}^4}{24}+\cdots ={\sum }_{k=0}^{\infty }\frac{(at{)}^k}{k!}$$

Recall $k!=1·2·3\cdots k$ is the factorial, and $0!=1$. We differentiate this series term by term

$$\frac{d}{dt}\left(e^{at}\right)=a+a^{2}t+\frac{a^3{t}^2}{2}+\frac{a^4{t}^3}{6}+\cdots =a(1+at\frac{(at{)}^2}{2}+\frac{(at{)}^3}{6}+\cdots )=a{e}^{at}.$$

Maybe we can try the same trick with matrices. Suppose that for an $n\times n$ matrix $A$ we define the matrix exponential as

$$e^A\stackrel{\text{def}}{=}𝐼+A+\frac{1}{2}{A}^2+\frac{1}{6}{A}^3+\cdots +\frac{1}{k!}{A}^k+\cdots$$

Let us not worry about convergence. The series really does always converge. We usually write $Pt$as $tP$ by convention when $P$ is a matrix. With this small change and by the exact same calculation as above we have that

$$\frac{d}{dt}\left(e^{tP}\right)=P{e}^{tP}.$$

Now $P$ and hence $e^{tP}$ is an $n\times n$ matrix. What we are looking for is a vector. We note that in the $1\times 1$ case we would at this point multiply by an arbitrary constant to get the general solution. In the matrix case we multiply by a column vector $\overrightarrow{c}$.

Let us check.

$$\frac{d}{dt}\overrightarrow{x}=\frac{d}{dt}\left(e^{tP}\overrightarrow{c}\right)=P{e}^{tP}\overrightarrow{c}=P\overrightarrow{x}.$$

Hence $e^{tP}$is the fundamental matrix solution of the homogeneous system. If we find a way to compute the matrix exponential, we will have another method of solving constant coefficient homogeneous systems. It also makes it easy to solve for initial conditions. To solve ${\overrightarrow{x}}^{\prime }=A\overrightarrow{x}$, $\overrightarrow{x}(0)=\overrightarrow{b}$ we take the solution

$$\overrightarrow{x}=e^{tA}\overrightarrow{b}$$

This equation follows because $e^{0A}=𝐼$, so $\overrightarrow{x}(0)=e^{0A}\overrightarrow{b}=\overrightarrow{b}$.

We mention a drawback of matrix exponentials. In general $e^{A+B}\ne e^A{e}^B$. The trouble is that matrices do not commute, that is, in general $AB\ne BA$. If you try to prove $e^{A+B}\ne e^A{e}^B$ using the Taylor series, you will see why the lack of commutativity becomes a problem. However, it is still true that if $AB=BA$, that is, if $A$ and $B$ commute, then $e^{A+B}=e^A{e}^B$. We will find this fact useful. Let us restate this as a theorem to make a point.

Simple cases

In some instances it may work to just plug into the series definition. Suppose the matrix is diagonal. For example, $D=\left[\begin{array}{cc}a& 0\\ 0& b\end{array}\right]$. Then

$$D^k=\left[\begin{array}{cc}{a}^k& 0\\ 0& b^k\end{array}\right]$$

and

$$\begin{array}{c}{e}^D=𝐼+D+\frac{1}{2}{D}^2+\frac{1}{6}{D}^3+\cdots \\ =\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]+\left[\begin{array}{cc}a& 0\\ 0& b\end{array}\right]+\frac{1}{2}\left[\begin{array}{cc}{a}^2& 0\\ 0& b^2\end{array}\right]+\frac{1}{6}\left[\begin{array}{cc}{a}^3& 0\\ 0& b^3\end{array}\right]+\cdots =\left[\begin{array}{cc}{e}^a& 0\\ 0& e^b\end{array}\right]\end{array}$$

So by this rationale we have that

$$e^{𝐼}=\left[\begin{array}{cc}e& 0\\ 0& e\end{array}\right]\text{and}{e}^{a𝐼}=\left[\begin{array}{cc}{e}^a& 0\\ 0& e^a\end{array}\right]$$

This makes exponentials of certain other matrices easy to compute. Notice for example that the matrix $A=\left[\begin{array}{cc}5& 4\\ -1& 1\end{array}\right]$ can be written as $3\mathit{I}\mathit{+}\mathit{B}$ where $𝐵=\left[\begin{array}{cc}2& 4\\ -1& -2\end{array}\right]$. Notice that ${\mathit{B}}^{\mathit{2}}=\left[\begin{array}{cc}0& 0\\ 0& 0\end{array}\right]$. So ${\mathit{B}}^{\mathit{k}}=0$ for all $k\ge 2$. Therefore, $e^B=\mathit{I}\mathit{+}\mathit{B}$. Suppose we actually want to compute $e^{tA}$. The matrices $3\mathit{t}\mathit{I}$ and $tB$ commute (exercise: check this) and $e^{tB}=\mathit{I}\mathit{+}\mathit{t}\mathit{B}$, since ${\left(tB\right)}^2=t^2{B}^2=0$. We write

$$\begin{array}{c}{e}^{tA}=e^{\mathit{3}\mathit{t}\mathit{I}+t𝐵}=e^{3t𝐼}{e}^{t𝐵}=\left[\begin{array}{cc}{e}^{3t}& 0\\ 0& e^{3t}\end{array}\right](I+tB)\\ =\left[\begin{array}{cc}{e}^{3t}& 0\\ 0& e^{3t}\end{array}\right]\left[\begin{array}{cc}1+2t& 4t\\ -t& 1-2t\end{array}\right]=\left[\begin{array}{cc}(1+2t)e^{3t}& 4t{e}^{3t}\\ -t{e}^{3t}& (1-2t)e^{3t}\end{array}\right]\end{array}$$

So we have found the fundamental matrix solution for the system ${\overrightarrow{x}}^{\prime }=A\overrightarrow{x}$. Note that this matrix has a repeated eigenvalue with a defect; there is only one eigenvector for the eigenvalue 3. So we have found a perhaps easier way to handle this case. In fact, if a matrix $A$ is $2\times 2$ and has an eigenvalue $\lambda$ of multiplicity 2, then either $A$ is diagonal, or $A=\lambda 𝐼+B$ where $B^2=0$. This is a good exercise.

Matrices $B$ such that $B^k=0$ for some $k$ are called nilpotent. Computation of the matrix exponential for nilpotent matrices is easy by just writing down the first $k$ terms of the Taylor series.

General Matrices

In general, the exponential is not as easy to compute as above. We usually cannot write a matrix as a sum of commuting matrices where the exponential is simple for each one. But fear not, it is still not too difficult provided we can find enough eigenvectors. First we need the following interesting result about matrix exponentials. For two square matrices $A$ and $B$, with $B$ invertible, we have

$$e^{BA{B}^{-1}}=B{e}^A{B}^{-1}.$$

This can be seen by writing down the Taylor series. First note that

$${\left(BA{B}^{-1}\right)}^2=BA{B}^{-1}BA{B}^{-1}=BA𝐼A{B}^{-1}=B{A}^2{B}^{-1}$$

And hence by the same reasoning ${\left(BA{B}^{-1}\right)}^k=B{A}^k{B}^{-1}$. Now write down the Taylor series for $e^{BA{B}^{-1}}$.

$$\begin{array}{c}{e}^{BA{B}^{-1}}=𝐼+BA{B}^{-1}+\frac{1}{2}{\left(BA{B}^{-1}\right)}^2+\frac{1}{6}{\left(BA{B}^{-1}\right)}^3+\cdots \\ =B{B}^{-1}+BA{B}^{-1}+\frac{1}{2}B{A}^2{B}^{-1}+\frac{1}{6}B{A}^3{B}^{-1}+\cdots \\ =B(𝐼+A+\frac{1}{2}{A}^2+\frac{1}{6}{A}^3+\cdots )B^{-1}\\ =B{e}^A{B}^{-1}.\end{array}$$

Given a square matrix $A$, we can sometimes write $A=ED{E}^{-1}$, where $D$ is diagonal and $E$ invertible. This procedure is called diagonalization. If we can do that, the computation of the exponential becomes easy. Adding $t$ into the mix we see that we can then easily compute the exponential

$$e^{tA}=E{e}^{tD}{E}^{-1}.$$

To diagonalize $A$ we will need $n$ linearly independent eigenvectors of $A$. Otherwise this method of computing the exponential does not work and we need to be trickier, but we will not get into such details. We let $E$ be the matrix with the eigenvectors as columns. Let ${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n$ be the eigenvalues and let ${\overrightarrow{v}}_1,\overrightarrow{v_2},\cdots ,{\overrightarrow{v}}_n$ be the eigenvectors, then $E=\left[\begin{array}{c}{\overrightarrow{v}}_1{\overrightarrow{v}}_2\cdots {\overrightarrow{v}}_n\end{array}\right]$. Let $D$ be the diagonal matrix with the eigenvalues on the main diagonal. That is

$$D=\left[\begin{array}{cccc}{\lambda }_1& 0& \cdots & 0\\ 0& {\lambda }_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\lambda }_n\end{array}\right]$$

We compute

$$\begin{array}{c}AE=A\left[\begin{array}{cccc}{\overrightarrow{v}}_1& {\overrightarrow{v}}_2& \cdots & {\overrightarrow{v}}_n\end{array}\right]\\ =\left[\begin{array}{cccc}A{\overrightarrow{v}}_1& A{\overrightarrow{v}}_2& \cdots & A{\overrightarrow{v}}_3\end{array}\right]\\ =\left[\begin{array}{cccc}{\lambda }_1{\overrightarrow{v}}_1& {\lambda }_2{\overrightarrow{v}}_2& \cdots & {\lambda }_n{\overrightarrow{v}}_n\end{array}\right]\\ =\left[\begin{array}{cccc}{\overrightarrow{v}}_1& {\overrightarrow{v}}_2& \cdots & {\overrightarrow{v}}_n\end{array}\right]D\\ =ED.\end{array}$$

The columns of $E$ are linearly independent as these are linearly independent eigenvectors of $A$. Hence $E$ is invertible. Since $AE=ED$, we right multiply by $E^{-1}$ and we get

$$A=ED{E}^{-1}.$$

This means that. $e^A=E{e}^D{E}^{-1}$ Multiplying the matrix by $t$ we obtain

The formula (3.8.1), therefore, gives the formula for computing the fundamental matrix solution $e^{tA}$ for the system ${\overrightarrow{x}}^{\prime }=A\overrightarrow{x}$, in the case where we have $n$ linearly independent eigenvectors.

Notice that this computation still works when the eigenvalues and eigenvectors are complex, though then you will have to compute with complex numbers. It is clear from the definition that if $A$ is real, then $e^{tA}$ is real. So you will only need complex numbers in the computation and you may need to apply Euler’s formula to simplify the result. If simplified properly the final matrix will not have any complex numbers in it.

Fundamental Matrix Solutions

We note that if you can compute the fundamental matrix solution in a different way, you can use this to find the matrix exponential $e^{tA}$. The fundamental matrix solution of a system of ODEs is not unique. The exponential is the fundamental matrix solution with the property that for $t=0$ we get the identity matrix. So we must find the right fundamental matrix solution. Let $X$ be any fundamental matrix solution to ${\overrightarrow{x}}^{\prime }=A\overrightarrow{x}$. Then we claim

$$e^{tA}=X(t)[X(0){]}^{-1}.$$

Clearly, if we plug $t=0$ into $X(t)[X(0){]}^{-1}$ we get the identity. We can multiply a fundamental matrix solution on the right by any constant invertible matrix and we still get a fundamental matrix solution. All we are doing is changing what the arbitrary constants are in the general solution $\overrightarrow{x}(t)=X(t)\overrightarrow{c}$.

Approximations

If you think about it, the computation of any fundamental matrix solution $X$ using the eigenvalue method is just as difficult as the computation of $e^{tA}$. So perhaps we did not gain much by this new tool. However, the Taylor series expansion actually gives us a very easy way to approximate solutions, which the eigenvalue method did not.

The simplest thing we can do is to just compute the series up to a certain number of terms. There are better ways to approximate the exponentialto the power 1. In many cases however, few terms of the Taylor series give a reasonable approximation for the exponential and may suffice for the application. For example, let us compute the first $4$ terms of the series for the matrix $$\begin{array}{c}12\\ 21\end{array}$$.

$$\begin{array}{c}{e}^{tA}\approx I+tA+\frac{t^2}{2}{A}^2+\frac{t^3}{6}{A}^3=I+t\left[\begin{array}{cc}1& 2\\ 2& 1\end{array}\right]+t^2\left[\begin{array}{cc}\frac{5}{2}& 2\\ 2& \frac{5}{2}\end{array}\right]+t^3\left[\begin{array}{cc}\frac{13}{6}& \frac{7}{3}\\ \frac{7}{3}& \frac{13}{6}\end{array}\right]=\\ =\left[\begin{array}{cc}1+t+\frac{5}{2}{t}^2+\frac{13}{6}{t}^3& 2t+2t^2+\frac{7}{3}{t}^3\\ 2t+2t^2+\frac{7}{3}{t}^3& 1+t+\frac{5}{2}{t}^2+\frac{13}{6}{t}^3\end{array}\right].\end{array}$$

Just like the scalar version of the Taylor series approximation, the approximation will be better for small $t$ and worse for larger $t$. For larger $t$, we will generally have to compute more terms. Let us see how we stack up against the real solution with $t=0.1$. The approximate solution is approximately (rounded to $8$ decimal places)

$$e^{0.1A}\approx 𝐼+0.1A+\frac{{0.1}^2}{2}+\frac{{0.1}^3}{6}{A}^3=\left[\begin{array}{cc}1.12716667& 0.22233333\\ 0.22233333& 1.12716667\end{array}\right]$$

And plugging $t=0.1$ into the real solution (rounded to $8$ decimal places) we get

$$e^{0.1A}=\left[\begin{array}{cc}1.12734811& 0.22251069\\ 0.22251069& 1.12734811\end{array}\right]$$

Not bad at all! Although if we take the same approximation for $t=1$ we get

$$𝐼+A+\frac{1}{2}{A}^2+\frac{1}{6}{A}^3=\left[\begin{array}{cc}6.66666667& 6.33333333\\ 6.33333333& 6.66666667\end{array}\right]$$

while the real value is (again rounded to $8$ decimal places)

$$e^A=\left[\begin{array}{cc}10.22670818& 9.85882874\\ 9.85882874& 10.22670818\end{array}\right]$$

So the approximation is not very good once we get up to $t=1$. To get a good approximation at $t=1$ (say up to $2$ decimal places) we would need to go up to the ${11}^{th}$ power (exercise).

Footnotes

[1] C. Moler and C.F. Van Loan, Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later, SIAM Review 45 (1), 2003, 3–49

Adapted from Differential Equations for Engineers by Jiří Lebl (Oklahoma State University), hosted on LibreTexts (math.libretexts.org) and licensed under CC BY-SA 4.0. Changes were made. License: CC-BY-SA-4.0.