\documentclass{article}
\input{../macros.tex}

\begin{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Lecture 23: Differential equations and $e^{\mb At}$}

The system of equations below describes how the values of
variables $u_1$ and $u_2$ affect each other over time:
\begin{eqnarray*}
\frac{du_1}{dt} = -u_1 + 2u_2\\
\frac{du_2}{dt} = u_1 - 2u_2.
\end{eqnarray*}
Just as we applied linear algebra to solve a difference equation, we can use it to solve this differential equation.  For example, the initial condition $u_1 = 1$, $u_2 = 0$ can be written $\vb u(0) =
\small \twobyone 10$.

\subsection*{Differential equations $\frac{d \vb u}{dt} = \mb A \vb u$}

By looking at the equations above, we might guess that over
time $u_1$ will decrease.
We can get the same sort of information more safely by looking at the eigenvalues
of the matrix $\mb A = \small \twobytwo {-1}21{-2}$ of our system $\displaystyle{\frac{d {\vb u}}{dt} = \mb A \vb u}$.  Because $\mb A$ is singular and
its trace is $-3$ we know that its eigenvalues are $\lambda_1 = 0$ and
$\lambda_2 = -3$.  The solution will turn out to include $e^{-3 t}$
and $e^{0t}$.  As $t$ increases, $e^{-3t}$ vanishes and $e^{0t}=1$
remains constant.  Eigenvalues equal to zero have eigenvectors that are {\em
  steady state} solutions.

$\vb x_1 = \twobyone 21$ is an eigenvector for which $\mb A \vb x_1 =
0 \vb x_1$.  To find an eigenvector corresponding to $\lambda_2 = -3$
we solve $(\mb A - \lambda_2 \mb I) \vb x_2 = \vb 0$:
$$\twobytwo 2211 \vb x_2 = 0 \quad \text{so} \quad \vb x_2 =
\twobyone 1{-1}$$ and we can check that $\mb A \vb x_2 = -3 \vb x_2$.
The general solution to this system of differential equations will be:
$$\vb u(t) = c_1 e^{\lambda_1 t}\vb x_1 + c_2 e^{\lambda_2 t}\vb x_2.$$ 

Is $e^{\lambda_1 t}\vb x_1$ really a solution to $\frac{d\vb u}{dt} =
\mb A \vb u$?  To find out, plug in $\vb u = e^{\lambda_1 t}\vb x_1$:
$$\frac{d \vb u}{dt} = \lambda_1 e^{\lambda_1 t} \vb x_1,$$
which agrees with:
$$\mb A \vb u = e^{\lambda_1t}\mb A \vb x_1 = \lambda_1
e^{\lambda_1 t}\vb x_1.$$ The two ``pure'' terms $e^{\lambda_1 t}\vb
x_1$ and $e^{\lambda_2 t}\vb x_2$ are analogous to the terms
$\lambda_i^k \vb x_i$ we saw in the solution $c_1 \lambda_1^k \vb x_1 +  c_2 \lambda_2^k \vb x_2 + \cdots +
 c_n \lambda_n^k\vb x_n$ to the difference equation
$\vb u_{k+1} = \mb A \vb u_k$.

Plugging in the values of the eigenvectors, we get:
$$\vb u(t) = c_1 e^{\lambda_1 t}\vb x_1 + c_2 e^{\lambda_2 t}\vb x_2 = c_1 \twobyone 21 + c_2 e^{-3t} \twobyone 1{-1}.$$
We know $\vb u(0) = \small\twobyone 10$, so at $t = 0$:
$$\twobyone 10 = c_1 \twobyone 21 + c_2 \twobyone 1{-1}.$$
$c_1 = c_2 = 1/3$ and 
$\vb u(t) = \frac{1}{3} \twobyone 21 + \frac{1}{3}e^{-3 t}\twobyone 1{-1}.$

This tells us that the system starts with $u_1 = 1$ and $u_2 = 0$ but
that as $t$ approaches infinity, $u_1$ decays to $2/3$ and $u_2$
increases to $1/3$.  This might describe stuff moving from $u_1$ to $u_2$.

The steady state of this system is $\vb u(\infty) = \small
\twobyone {2/3}{1/3}$.

\subsubsection*{Stability}

Not all systems have a steady state.  The eigenvalues of $\mb A$ will
tell us what sort of solutions to expect:
\begin{enumerate}
\item Stability:  $\vb u(t) \rightarrow 0$ when $\text{Re}
(\lambda) < 0$.
\item Steady state:  One eigenvalue is $0$ and all
other eigenvalues have negative real part.
\item Blow up:  if $\text{Re}(\lambda)>0$ for any eigenvalue $\lambda$.
\end{enumerate}

If a two by two matrix $\mb A = \small \twobytwo abcd$ has two
eigenvalues with negative real part, its trace $a+d$ is negative.  The
converse is not true: $\small \twobytwo {-2}001$ has negative trace
but one of its eigenvalues is $1$ and $e^{1t}$ blows up.  If $\mb A$
has a positive determinant and negative trace then the corresponding
solutions must be stable.

\subsubsection*{Applying $\mb S$}

The final step of our solution to the system $\frac{d \vb u}{dt} = \mb A \vb u$ was to solve:
$$c_1 \twobyone 21 + c_2 \twobyone 1{-1} = \twobyone 10.$$
In matrix form:
$$\twobytwo 211{-1} \twobyone {c_1}{c_2} = \twobyone 10.$$ or $\mb S
\vb c = \vb u(0)$, where $\mb S$ is the eigenvector matrix.  The
components of $\vb c$ determine the contribution from each pure
exponential solution, based on the initial conditions of the system.

In the equation $\frac{d \vb u}{dt} = \mb A \vb u$, the matrix $\mb A$ couples the pure solutions.  We set $\vb u = \mb S \vb v$, where $\mb S$ is the matrix of eigenvectors of $\mb A$, to get:
$$\mb S \frac{d \vb v}{dt} = \mb A \mb S \vb v$$
or:
$$\frac{d \vb v}{dt} = \mb S^{-1}\mb A\mb S \vb v = \mb \Lambda \vb
v.$$ This diagonalizes the system: $\frac{d v_i}{dt} = \lambda_i v_i$.
The general solution is then:
\begin{eqnarray*}
\vb v(t) &=& e^{\mb \Lambda t}\vb v(0), \qquad \text{and}\\
\vb u(t) &=& \mb S e^{\mb \Lambda t} \mb S^{-1} \vb v(0) = e^{\mb A t} \vb u(0).
\end{eqnarray*}

\subsection*{Matrix exponential $e^{\mb A t}$}

What does $e^{\mb At}$ mean if $\mb A$ is a matrix?  We know that for
a real number $x$, $$e^x = \sum_{n=0}^\infty \frac{x^n}{n!} = 1 + x +
\frac{x^2}{2} + \frac{x^3}{6} + \cdots.$$ We can use the same formula
to define $e^{\mb At}$: $$e^{\mb At} = \mb I + \mb At + \frac{(\mb
  At)^2}{2} + \frac{(\mb At)^3}{6} + \cdots.$$

\smallskip

Similarly, if the eigenvalues of $\mb At$ are small, we can use the
geometric series $\displaystyle{\frac{1}{1-x} = \sum_{n=0}^\infty x^n}$ to estimate
$(\mb I - \mb A t)^{-1} = \mb I + \mb At + (\mb At)^2 + (\mb At)^3 +
\cdots$.

\smallskip

We've said that $e^{\mb A t} = \mb S e^{\mb \Lambda t}S^{-1}$.  If
$\mb A$ has $n$ independent eigenvectors we can prove this from the
definition of $e^{\mb At}$ by using the formula $\mb A = \mb S \mb \Lambda \mb
S^{-1}$:
\begin{eqnarray*}
e^{\mb At} &=& 
\mb I + \mb At + \frac{(\mb At)^2}{2} + \frac{(\mb At)^3}{6} + \cdots\\
 &=& 
\mb S \mb S^{-1} + \mb S \mb \Lambda \mb S^{-1}t + \frac{\mb S \mb \Lambda^{2} \mb S^{-1}}{2}t^2 + \frac{\mb S \mb \Lambda^3 \mb S^{-1}}{6}t^3 + \cdots\\
 &=& 
\mb S e^{\mb \Lambda t}\mb S^{-1}.
\end{eqnarray*}

It's impractical to add up infinitely many matrices.  Fortunately,
there is an easier way to compute $e^{\mb \Lambda t}$.  Remember that:
$$\mb \Lambda = \left[ \begin{array}{rrrr} 
\lambda_1 & 0 & \cdots & 0\\
0 & \lambda_2 &  & 0\\
\vdots &  & \ddots & \vdots\\
0 & \cdots & 0 & \lambda_n\end{array}\right].$$
When we plug this in to our formula for $e^{\mb At}$ we find that:
$$e^{\mb \Lambda t} = \left[ \begin{array}{rrrr} 
e^{\lambda_1 t}& 0 & \cdots & 0\\
0 & e^{\lambda_2 t} &  & 0\\
\vdots &  & \ddots & \vdots\\
0 & \cdots & 0 & e^{\lambda_nt}\end{array}\right].$$

This is another way to see the relationship between the stability of
$\vb u(t) = \mb S e^{\mb \Lambda t} \mb S^{-1} \vb v(0)$ and the
eigenvalues of $\mb A$.

\subsection*{Second order}

We can change the second order equation $y'' + by' + ky = 0$ into a
two by two first order system using a method similar to the one we
used to find a formula for the Fibonacci numbers.  If $\mb u = \small
\twobyone {y'}y$, then $$\mb u' = \small \twobyone {y''}{y'} =
\twobytwo {-b}{-k}10 \twobyone{y'}{y}.$$ We could use the methods we
just learned to solve this system, and that would give us a solution
to the second order scalar equation we started with.

If we start with a $k$th order equation
we get a $k$ by $k$ matrix with coefficients of the equation in the
first row and $1$'s on a diagonal below that; the rest of the entries
are $0$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document}