\documentclass[10pt]{article}

\usepackage{systeme}

\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{amssymb}
%\usepackage{palatino}
\usepackage{setdeck}

\usepackage{xcolor}
\usepackage{tcolorbox}

%\usepackage[textheight=8.5in, margin=1.6in]{geometry}
\usepackage[total={6in,8in}]{geometry}

\usepackage{tikz}
\usetikzlibrary{cd}
\tikzset{
      labr/.style={anchor=north, rotate=90, inner sep=.7mm}
    }
\tikzset{
      labl/.style={anchor=south, rotate=90, inner sep=.7mm}
    }

%\usepackage{pgfplots}
\tikzstyle{ball} = [circle,shading=ball, ball color=black,
    minimum size=1mm,inner sep=1.3pt]
%\usetikzlibrary{arrows.meta}
%\tikzset{>={Latex[length=2mm,width=1.5mm]}}

%\usepackage{graphicx}

\usepackage{needspace}  % <-- to prevent bad pagebreak over environment titles
% Use:  \Needspace{3\baselineskip} before opening an environment

\usepackage[shortlabels]{enumitem}
\setlist[enumerate]{leftmargin=*,label=(\alph*)}

%\usepackage[colorlinks=true,citecolor=black,linkcolor=black,urlcolor=blue]{hyperref}

\newcommand{\N}{\mathbb{N}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\sn}{\mathfrak{S}}
\newcommand{\PP}{\mathbb{P}}
\newcommand{\ve}{\varepsilon}
\newcommand{\ips}[2]{\langle{#1}, {#2}\rangle}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator{\Span}{\mathrm{Span}}
\DeclareMathOperator{\rk}{\mathrm{rank}}
\DeclareMathOperator{\rspace}{\mathrm{rowspace}}
\DeclareMathOperator{\cspace}{\mathrm{colspace}}
\DeclareMathOperator{\trace}{\mathrm{tr}}
\DeclareMathOperator{\im}{\mathrm{im}}
\DeclareMathOperator{\Sym}{\mathrm{Sym}}
\DeclareMathOperator{\Hom}{\mathrm{Hom}}
\DeclareMathOperator{\diag}{\mathrm{diag}}
\DeclareMathOperator{\glb}{\mathrm{glb}}
\DeclareMathOperator{\lub}{\mathrm{lub}}

\parskip = 0.05in
\parindent = 0.0in

%\pagestyle{empty}

\allowdisplaybreaks

\begin{document}
\centerline{Math 201 lecture for Monday, Week 13}
\bigskip

\centerline{\textbf{Systems of linear differential equations}}  
\bigskip
 
Suppose that
\[
  x(t) = \text{amount of yeast at time~$t$}
\]
and that rate of growth of yeast (at least in the time frame in which we are
interested) is proportional to the amount of yeast.  So there exists a
constant~$a$ such that
\[
  x'(t)=ax(t).
\]
Integrating, we get
\[
  \int\frac{x'(t)}{x(t)}\,dt=\int a\,dt\quad\Rightarrow\quad \ln(x(t))=at+b
\]
for some constant~$b$.  Exponentiating then yields
\[
  x(t)=e^{at}c
\]
where~$c=e^{b}$. Evaluating at~$t=0$ shows that~$c$ is the {\em initial
condition}: $x(0)=c$.

Now consider a two-dimensional system.  Let
\begin{align*}
  x_1(t)&=\ \text{population of frogs in a pond}\\
  x_2(t)&=\ \text{population of flies in a pond},
\end{align*}
and suppose the rate of change of these populations satisfies the following
system of differential equations:
\begin{align*}
  x_1'(t)&=ax_1(t)+bx_2(t)\\
  x_2'(t)&=cx_1(t)+dx_2(t).
\end{align*}
So we are assuming that the rate of growth of these populations depends linearly
on the sizes of the populations.  Letting
\[
  x(t):=
  \left(\begin{array}{c}
      x_1(t)\\x_2(t)
  \end{array} \right)\quad\text{and}\quad
  x'(t):=
  \left(\begin{array}{c}
      x_1'(t)\\x_2'(t)
  \end{array} \right),
\]
we can rewrite the system in matrix form:
\[
  x'(t)=  
  A x(t)
\]
where
\[
  A = 
  \left(\begin{array}{cc}
      a&b\\c&d
  \end{array} \right).
\]
Our problem is to find~$x(t)$.  The key to solving higher-dimensional systems
like this is the following:
\medskip

{\bf Theorem.} Let~$A$ be an~$n\times n$ matrix over the real or complex
numbers.  Then the solution to~$x'=Ax$ with initial condition~$x(0)=p$ is
\[
  x=e^{At}p.
\]
(Note that~$p$ is a column vector here.) 
\medskip

To make sense of this, we need to be able to exponentiate a matrix!  To do that,
recall that for a real or complex number~$a$, we have
\[
  e^a = \sum_{k=0}^{\infty}\frac{1}{k!}a^k,
\]
an infinite series that converges for all~$a$.  This formula generalizes: 
given any~$n\times n$ matrix~$A$ over the real or complex
numbers, define
\[
  e^{At}:=\sum_{k=0}^{\infty}\frac{1}{k!}(At)^k
  =\sum_{k=0}^{\infty}\frac{1}{k!}A^kt^k
  =I_n+At+\frac{1}{2}A^2t^2+\frac{1}{6}A^3t^2+\frac{1}{24}A^4t^4+\cdots.
\]
Each entry of~$e^{At}$ is a power series in~$t$, and that power series turns out to converge for
all~$t$.  To compute~$e^{At}$ though, we need to somehow compute all of the
powers of~$A$. As you might expect, diagonalization comes to the rescue.

{\bf Computing~$e^{At}$.} If~$A$ is diagonalizable, then we can write
\[
  P^{-1}AP=D=\diag(\lambda_1,\dots,\lambda_n)
\]
where the~$\lambda_i$ are the eigenvalues of~$A$. As we have seen earlier, it
follows that
\[
  A^k=(PDP^{-1})^k=PD^kP^{-1}=P\diag(\lambda_1^k,\dots,\lambda_n^k)P^{-1}.
\]
Therefore, modulo some technicalities involving convergence, we have
\[
  e^{At}
  =\sum_{k=0}^{\infty}\frac{1}{k!}A^kt^k
  =\sum_{k=0}^{\infty}\frac{1}{k!}(PD^kP^{-1})t^k
  =P\left(\sum_{k=0}^{\infty}\frac{1}{k!}D^kt^k\right)P^{-1}
  =Pe^{Dt}P^{-1}.
\]
Since~$D$ is diagonal, an easy calculation shows that
\[
  e^{Dt}=\diag(e^{\lambda_1t},\dots,e^{\lambda_nt}).
\]
So
\[
  e^{At}=P\diag(e^{\lambda_1t},\dots,e^{\lambda_nt})P^{-1}.
\]
\medskip

{\bf Example.} Consider the following two-dimensional system:
\begin{align*}
  x_1'&=x_2\\
  x_2'&=x_1.
\end{align*}
In matrix form,
\[
x'=Ax
\]
where
\[
  A=
\left(\begin{array}{cc}
0&1\\
1&0
\end{array} \right).
\]
Applying our algorithm to diagonalize~$A$, we find
\[
  P^{-1}AP=D=\diag(1,-1)
\]
where
\[
  P=
  \left(\begin{array}{rr}
      1&1\\
      1&-1
  \end{array} \right).
\]
Therefore,
\begin{align*}
  e^{At}
  &= Pe^{Dt}P^{-1}\\[10pt]
  &=
\left(\begin{array}{rr}
    1&1\\1&-1
\end{array} \right)
\left(\begin{array}{cc}
    e^t&0\\0&e^{-t}
\end{array} \right)
\left(\begin{array}{rr}
    \frac{1}{2}&\frac{1}{2}\\[5pt]
    \frac{1}{2}&-\frac{1}{2}
\end{array}\right) \\[10pt]
&=
\frac{1}{2}
\left(\begin{array}{cc}
    e^t+e^{-t}&e^t-e^{-t}\\
    e^t-e^{-t}&e^t+e^{-t}
\end{array} \right).
\end{align*}
So, for example, the solution with initial condition~$x(0)=(1,0)$ is
\[
  \left(\begin{array}{c}
      x_1(t)\\
      x_2(t)
  \end{array} \right)
  =
\frac{1}{2}
\left(\begin{array}{cc}
    e^t+e^{-t}&e^t-e^{-t}\\
    e^t-e^{-t}&e^t+e^{-t}
\end{array} \right)
\left(\begin{array}{c}
    1\\
    0
\end{array} \right)
=
\frac{1}{2}
\left(\begin{array}{c}
    e^t+e^{-t}\\
    e^t-e^{-t}
\end{array} \right),
\]
A plot of that solution~$(x_1(t),x_2(t))$ in the plane
appears in blue in the picture below.  The arrows indicate the following: at
each point~$(x_1,x_2)\in\R^2$, we attach the velocity vector~$(x_1',x_2')=(x_2,x_1)$.
\begin{center}
  \includegraphics[height=2in]{images/13mon1.png}
\end{center}
The solution in blue has velocity vector~$x'(0)=(x_2(0),x_1(0))=(0,1)$ at time~$t=0$.  To repeat:
geometrically, our solution is a parametrized curve in the
plane:
\begin{align*}
  x\colon\R&\to\R^2\\
  t&\mapsto x(t)=(x_1(t),x_2(t)).
\end{align*}
The differential equation specifies the tangent (velocity) vectors~$x'(t)$
at each time~$t$. It determines a ``flow'' as illustrated in the picture.
Specifying an initial condition is like dropping a speck into the flow.  We then
get a unique solution, which is the trajectory of that speck over time (shown in
blue, above).

{\bf Note:} the arrows determine new ``axes'' pointed in the directions of the
eigenvectors,~$(1,1)$ and~$(1,-1)$.  

{\bf Example.} Next consider the following two-dimensional system:
\begin{align*}
  x_1'&=x_2\\
  x_2'&=-x_1.
\end{align*}
It might approximate frog-fly populations since one would expect the frog
population~$x_1(t)$ to increase with the fly population and the fly population
to decrease with the frog population.
In matrix form,
\[
x'=Ax
\]
where
\[
  A=
\left(\begin{array}{rr}
0&1\\
-1&0
\end{array} \right).
\]
The characteristic polynomial is~$p_A(x)=x^2+1$, so~$A$ is not diagonalizable
over~$\R$.  However, it is diagonalizable over~$\C$.  So let's do that to see
where that goes.
Applying our algorithm to diagonalize~$A$, we find
\[
  P^{-1}AP=D=\diag(i,-i)
\]
where
\[
  P=
  \left(\begin{array}{rr}
      i&-i\\
      1&1
  \end{array} \right).
\]
\begin{align*}
  e^{At}
  &= Pe^{Dt}P^{-1}\\[10pt]
  &=
  \left(\begin{array}{rr}
      i&-i\\
      1&1
  \end{array} \right)
  \left(\begin{array}{rr}
      e^{i \, t} & 0 \\
      0 & e^{-i \, t}
  \end{array}\right)
  \left(\begin{array}{rr}
      -\frac{1}{2} i & \frac{1}{2} \\[5pt]
      \frac{1}{2} i & \frac{1}{2}
  \end{array}\right)\\[8pt]
&= 
\left(\begin{array}{rr}
    \frac{1}{2} \, e^{i \, t} + \frac{1}{2} \, e^{-i \,
    t} & -\frac{1}{2} i \, e^{i \, t} + \frac{1}{2} i \,
    e^{-i \, t} \\
    \frac{1}{2} i \, e^{i \, t} - \frac{1}{2} i \, e^{-i \,
    t} & \frac{1}{2} \, e^{i \, t} + \frac{1}{2} \,
    e^{-i \, t}
\end{array}\right)\\[8pt]
&= 
\left(\begin{array}{rr}
    \cos(t)&\sin(t)\\
    -\sin(t)&\cos(t)
\end{array} \right).
\end{align*}
So starting with equal populations of frogs and flies,~$x(0)=(1,1)$, we have
\[
  x(t)=
\left(\begin{array}{rr}
    \cos(t)&\sin(t)\\
    -\sin(t)&\cos(t)
\end{array} \right)
\left(\begin{array}{c}
    1\\1
\end{array} \right)
=(\cos(t)+\sin(t),-\sin(t)+\cos(t)).
\]
\begin{center}
  \includegraphics[height=2in]{images/13mon2.png}
\end{center}
Note how this system of equations is not a great model for frogs and flies:
starting at any initial population, the system evolves into one in which there
are negative amounts of frogs or flies.  One could hope that it applies locally,
say near times at which the populations for frogs and flies is nearly equal.  At
any rate, it raises the question as to whether any \emph{linear} system of
equations would make a good model.  Qualitatively, what are all of the
possibilities for a two-dimensional linear system?
\end{document}  
