13.7. 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.

\[\begin{split}\left\{\begin{aligned} 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 \end{aligned}\right.\end{split}\]

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.7.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,

\[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!} + \cdots.\]

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{split}\begin{aligned} 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!} + \cdots \\ \hfill &= \mathbf{X}\,\left( \mathbf{I} + \mathbf{\Lambda} + \frac{\mathbf{\Lambda}^2}{2!} + \frac{\mathbf{\Lambda}^3}{3!} + \cdots \right) \,\mathbf{X}^{-1} \\ \hfill &= \mathbf{X}\,e^{\mathbf{\Lambda}} \,\mathbf{X}^{-1} \end{aligned}\end{split}\]

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!} + \cdots \right)\]

Adding the elements in the Maclaurin series of \(\bf{\Gamma}\) reveals a Maclaurin series for \(e^{\lambda_i}\) in each diagonal element.

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

In MATLAB, if variable L (for Lambda) is used for \(\bf{\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.

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

13.7.2. 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.8)\[\mathbf{y}(t) = e^{\mathbf{A}\,t}\bm{k} = \mathbf{X\,\Gamma^t\,X}^{-1}\bm{k},\]

where

\[\begin{split}\mathbf{\Gamma}^t = \begin{bmatrix} 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} \end{bmatrix}.\end{split}\]

Using equation (13.8), 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 solution used in Systems of Linear ODEs (equation (8.10)).

(13.9)\[\begin{split}\begin{aligned} \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 \end{aligned}\end{split}\]

13.7.3. Example ODE Matrix Solution

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

\[\begin{split}\left\{\begin{aligned} y_1(t)' &= -2\,y_1(t) + {y_2(t)}{,} \\ y_2(t)' &= y_1(t) - {2\,y_2(t)}{,} \end{aligned}\right. \qquad \begin{aligned} y_1(0) &= 6 \\ y_2(0) &= 2 \end{aligned}\end{split}\]

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
>> % scaling for appearance
>> X = X/X(1,1)
X =
    1.0000    1.0000
   -1.0000    1.0000
>> y0 = [6; 2]
y0 =
    6
    2
>> c = X\y0
c =
    2.0000
    4.0000

\[\begin{split}\bm{y}(t) = \mathbf{X\, \Gamma}^t \, \bm{c} = \begin{bmatrix} 1 & 1 \\ -1 & 1 \end{bmatrix}\,\begin{bmatrix} {e^{-3t}} & 0 \\ 0 & {e^{-t}} \end{bmatrix} \,\begin{bmatrix} 2 \\ 4 \end{bmatrix} = \begin{bmatrix} {e^{-3t}} & {e^{-t}} \\ {-e^{-3t}} & {e^{-t}} \end{bmatrix} \,\begin{bmatrix} 2 \\ 4 \end{bmatrix}\end{split}\]

\[\begin{split}\left\{\begin{aligned} y_1(t) &= 2\,e^{-3t} + 4\,e^{-t} \\ y_2(t) &= -2\,e^{-3t} + 4\,e^{-t} \end{aligned}\right.\end{split}\]