13.6. A Matrix Exponent and Systems of ODEs

We have already learned from Systems of Linear ODEs that ordinary differential equations (ODEs) of the form

\frac{dy(t)}{dt} = a\,y(t),

have the solution

y(t) = c\,e^{a\,t}.

The same principle applies to systems of ODEs, except that we use vectors and matrices to describe the equations.

\spalignsys{y_1' = a_{11}\,y_1 + a_{12}\,y_2 + \cdots + a_{1n}\,y_n; y_2' = a_{21}\,y_1 + a_{22}\,y_2 + \cdots + a_{2n}\,y_n; \. \vdots; y_n' = a_{n1}\,y_1 + a_{n2}\,y_2 + \cdots + a_{nn}\,y_n}

In matrix notation, this is \bm{y'} = \mathbf{A}\,\bm{y}, which has a solution of the form \bm{y}(t) = e^{\mathbf{A}\,t}\,\bm{k}.

We will show here that the matrix form of the solution leads to the same general solution given in Systems of Linear ODEs. First, let us put the matrix exponent e^{\bf{A}} into a more manageable form and then apply it to the ODE.

13.6.1. A Matrix in the Exponent

The equation e^{\bf{A}} seems unusual. What does it mean to have a matrix in the exponent of the natural number e?

We need to use a Maclaurin series expansion of e^{\bf{A}} to get the matrix out of the exponent and then we will use diagonalization from Diagonalization to simplify the powers of \bf{A} that the series expansion gives us.

From the Maclaurin series expansion we have,

e^{\mathbf{A}} = \sum_{n=0}^{\infty} \frac{\mathbf{A}^n}{n!} = \mathbf{I} + \mathbf{A} + \frac{\mathbf{A}^2}{2!} + \frac{\mathbf{A}^3}{3!} + \ldots.

From Diagonalization and Powers of A, we know that \mathbf{A}^k = \mathbf{X\,\Lambda}^k\,\mathbf{X}^{-1} Where \mathbf{X} and \mathbf{\Lambda} are the eigenvector and diagonal eigenvalue matrices of \bf{A}.

Thus,

\begin{array}{rl} e^{\mathbf{A}} &= \mathbf{I} + \mathbf{X\,\Lambda\,X}^{-1} + \frac{\mathbf{X\,\Lambda}^2\,\mathbf{X}^{-1}}{2!} + \frac{\mathbf{X\,\Lambda}^3\,\mathbf{X}^{-1}}{3!} + \ldots \\ &= \mathbf{X}\,\left( \mathbf{I} + \mathbf{\Lambda} + \frac{\mathbf{\Lambda}^2}{2!} + \frac{\mathbf{\Lambda}^3}{3!} + \ldots \right) \,\mathbf{X}^{-1} \\ &= \mathbf{X}\,e^{\mathbf{\Lambda}} \,\mathbf{X}^{-1} \end{array}

We will call the middle term \bf{\Gamma} (Gamma).

\mathbf{\Gamma} = e^{\mathbf{\Lambda}} = \left( \mathbf{I} + \mathbf{\Lambda} + \frac{\mathbf{\Lambda}^2}{2!} + \frac{\mathbf{\Lambda}^3}{3!} + \ldots \right)

If we add the elements in \bf{\Gamma}, we will find a Maclaurin series for e^{\lambda_i} in each diagonal element.

\begin{array}{rl} \bf{\Gamma} &= \spalignmat{ {\left( 1 + \lambda_1 + \frac{\lambda_1^2}{2!} + \frac{\lambda_1^3}{3!} + \ldots \right)} 0 \cdots{} 0; 0 {\left( 1 + \lambda_2 + \frac{\lambda_2^2}{2!} + \frac{\lambda_2^3}{3!} + \ldots \right)} \cdots{} 0; \vdots{} \vdots{} \ddots{} \vdots{}; 0 0 \cdots{} {\left( 1 + \lambda_n + \frac{\lambda_n^2}{2!} + \frac{\lambda_n^3}{3!} + \ldots \right)}} \\ \hfill \\ &= \spalignmat{e^{\lambda_1} 0 \cdots{} 0; 0 e^{\lambda_2} \cdots{} 0; \vdots{} \vdots{} \ddots{} \vdots{}; 0 0 \cdots{} e^{\lambda_n}} \end{array}

In MATLAB, if variable L (Lambda) is used for \Lambda, then Gamma = exp(L) .* eye(n), where n is the size of the matrix. Or if we have the variable e = exp(1), then Gamma = e.^L .* eye(n). The expression Gamma = expm(L) gives the same result using a built-in MATLAB function.

Finally,

\boxed{ e^{\mathbf{A}} = \mathbf{X\,\Gamma\,X}^{-1} }

13.6.2. Example Matrix Exponent

>> A = [10 0; -10 7]
A =
    10     0
   -10     7
>> [X, L] = eig(A)
X =
         0    0.2873
    1.0000   -0.9578
L =
    7     0
    0    10

>> e = exp(1)
e =
    2.7183
>> Gamma = e.^L .* eye(2)
Gamma =
   1.0e+04 *
    0.1097         0
         0    2.2026

>> % find e^A
>> X*Gamma*inv(X)
ans =
   1.0e+04 *
    2.2026         0
   -6.9766    0.1097

>> % verify with MATLAB's expm() function
>> expm(A)
ans =
   1.0e+04 *
    2.2026         0
   -6.9766    0.1097

13.6.3. Matrix Solution to a System of ODEs

Returning to the problem of systems of ODEs, we know that in matrix notation, the ODE \bm{y'} = \mathbf{A}\,\bm{y}, has a solution of the form

(13.5)\mathbf{y}(t) = e^{\mathbf{A}\,t}\bm{k} = \mathbf{X\,\Gamma^t\,X}^{-1}\bm{k},

where

\mathbf{\Gamma}^t = \spalignmat{e^{\lambda_1\,t} 0 \cdots{} 0; 0 e^{\lambda_2\,t} \cdots{} 0; \vdots{} \vdots{} \ddots{} \vdots{}; 0 0 \cdots{} e^{\lambda_n\,t}}.

Note

\begin{array}{ll} e^{\mathbf{A}} &= \mathbf{X\,\Gamma\,X}^{-1} \\ e^{\mathbf{A}\,2} &= (\mathbf{X\,\Gamma\,X}^{-1})\,(\mathbf{X\,\Gamma\,X}^{-1}) \\ &= \mathbf{X\,\Gamma^2\,X}^{-1} \\ e^{\mathbf{A}\,3} &= (\mathbf{X\,\Gamma^2\,X}^{-1})\,(\mathbf{X\,\Gamma\,X}^{-1}) \\ &= \mathbf{X\,\Gamma^3\,X}^{-1} \\ &\vdots \\ e^{\mathbf{A}\,t} &= \mathbf{X\,\Gamma^t\,X}^{-1} \end{array}

Using (13.5), we find \bm{k} from \bm{y}(t) when t = 0.

\bm{y}(0) = \mathbf{X\,\Gamma^0\,X}^{-1}\bm{k} = \mathbf{X\,I\,X}^{-1}\bm{k} = \bm{k}

Letting \bm{c} = \mathbf{X}^{-1}\,\bm{k} = \mathbf{X}^{-1}\,\bm{y}(0) gives the same solution as we used in Systems of Linear ODEs

\boxed{ \bm{y}(t) = \mathbf{X\,\Gamma}^t\,\mathbf{X}^{-1}\,\bm{y}(0) = \mathbf{X\, \Gamma}^t \, \bm{c} = c_1 e^{\lambda_1\,t}\bm{x}_1 + c_2 e^{\lambda_2\,t}\bm{x}_2 + \cdots + c_n e^{\lambda_n\,t}\bm{x}_n }

13.6.4. Example ODE Matrix Solution

In Systems of Linear ODEs, we solved the following ODE example with initial conditions,

\spalignsys{y_1(t)' = -2\,y_1(t) + y_2(t); y_2(t)' = y_1(t) - 2\,y_2(t)} \mbox{, \hspace{0.5in}} \spalignsys{y_1(0) = 6; y_2(0) = 2}.

Here we keep all of the variables as matrices and vectors in MATLAB and find the same solution.

>> A = [-2 1; 1 -2]
A =
    -2     1
     1    -2
>> [X, L] = eig(A)
X =
    0.7071    0.7071
   -0.7071    0.7071
L =
    -3     0
     0    -1
>> X = X*2/sqrt(2)    % scaling for appearance - not needed
X =
    1.0000    1.0000
   -1.0000    1.0000
>> y0 = [6; 2]
y0 =
    6
    2
>> c = X\y0
c =
    2.0000
    4.0000

\bm{y}(t) = \mathbf{X\, \Gamma}^t \, \bm{c} = \mat{1 1; -1 1}\,\mat{{e^{-3t}}, 0; 0, {e^{-t}}} \,\vector{2;4} = \mat{{e^{-3t}}, {e^{-t}}; {-e^{-3t}}, {e^{-t}}} \,\vector{2;4}

\spalignsys{y_1(t) = 2\,e^{-3t} + 4\,e^{-t}; y_2(t) = -2\,e^{-3t} + 4\,e^{-t}}