\documentclass[leqno,11pt]{article}
 \setlength{\textheight}{10.4in}
 \setlength{\topmargin}{-1in}
 \setlength{\headheight}{0in}
 \setlength{\topsep}{0in}
 \setlength{\marginparwidth}{0in}
 \setlength{\marginparsep}{0in}
 \setlength{\oddsidemargin}{0in}
 \setlength{\evensidemargin}{0in}
 \setlength{\textwidth}{6.4in}
 \setlength{\parskip}{4pt}
 \setlength{\parindent}{0pt}
\usepackage{alltt}
%\usepackage{email}
\usepackage{hyperref}
\usepackage{graphicx}
\usepackage{quiz}
\newif\ifsolutions
 \solutionstrue
 \solutionsfalse
\begin{document}
\section*{Sample Quiz 13}

{\Large\bf Problem 1}.

\begin{tabular}{lr}
\includegraphics[scale=0.5]{figures/brine-tank-star-color-7tanks.jpg}
&
\begin{minipage}[b]{0.6\textwidth}
\begin{itemize}
\item[]
Flow through each pipe is $f$ gallons per unit time.
\item[]
Each tank has constant volume $V$.
\item[]
Symbols $x_1(t)$ to $x_7(t)$ are the salt amounts in tanks $T_1$ to $T_7$, respectively.
\end{itemize}

\vspace*{4pc}
\end{minipage}
\end{tabular}

The differential equations are obtained by the classical balance law, which says that the rate
of change in salt amount is the rate in minus the rate out.
Individual rates in/out are of the form (flow rate)(salt concentration), where
flow rate $f$ has units volume per unit time and $x_i(t)/V$ is
the concentration = amount/volume.

$$
\begin{array}{lcl}
x_1'(t) &=& \frac{f}{V}\left(x_2(t)+x_3(t)+x_4(t)+x_5(t)+x_6(t)+x_7(t)  -6x_1(t)\right) \\[1.5ex]
x_2'(t) &=& \frac{f}{V}\left(x_1(t) - x_2(t) \right),\\[1.5ex]
x_3'(t) &=&  \frac{f}{V}\left(x_1(t) - x_3(t) \right),\\[1.5ex]
x_4'(t) &=&  \frac{f}{V}\left(x_1(t) - x_4(t) \right),\\[1.5ex]
x_5'(t) &=&  \frac{f}{V}\left(x_1(t) - x_5(t) \right),\\[1.5ex]
x_6'(t) &=&  \frac{f}{V}\left(x_1(t) - x_6(t) \right),\\[1.5ex]
x_7'(t) &=&  \frac{f}{V}\left(x_1(t) - x_7(t) \right).
\end{array}
$$
\rule{\textwidth}{0.5ex}

{\large\bf Problem 1(a)}. Change variables $t=Vr/f$ to obtain the new system

$$
\begin{array}{lcl}
\dd\frac{dx_1}{dr} &=& x_2+x_3+x_4+x_5+x_6+x_7  -6x_1 \\[1.5ex]
\dd\frac{dx_2}{dr} &=& x_1 - x_2 ,\\[1.5ex]
\dd\frac{dx_3}{dr} &=&  x_1 - x_3 ,\\[1.5ex]
\dd\frac{dx_4}{dr} &=&  x_1 - x_4 ,\\[1.5ex]
\dd\frac{dx_5}{dr} &=&  x_1 - x_5 ,\\[1.5ex]
\dd\frac{dx_6}{dr} &=&  x_1 - x_6 ,\\[1.5ex]
\dd\frac{dx_7}{dr} &=&  x_1 - x_7 .
\end{array}
$$

\ifsolutions
{\bf Solution 1(a)}:

Because $\dd \frac{dx(t)}{dt}=\frac{dx}{dr}\frac{dr}{dt}=\frac{dx}{dr}\frac{f}{V}$, then
the fraction $f/V$ cancels in the equations, resulting in the new system.
\fi

\vfill\newpage
\rule{\textwidth}{0.5ex}

{\large\bf Problem 1(b)}. Formulate the equations in 1(a) in the system form
$\dd\frac{d}{dr}\vec{u} = A\vec{u}$.

{\bf Answer}:

$$ A= \left(
\begin{array}{rrrrrrrr}
-6 & 1 & 1 & 1 & 1 & 1 & 1 \\
1 & -1 & 0 & 0 & 0 & 0 & 0 \\
1 & 0 & -1 & 0 & 0 & 0 & 0 \\
1 & 0 & 0 & -1 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 & -1 & 0 & 0 \\
1 & 0 & 0 & 0 & 0 & -1 & 0 \\
1 & 0 & 0 & 0 & 0 & 0 & -1
\end{array}
\right),\quad \vec{u}=\colvector{x_1\\ x_2\\ x_3\\ x_4 \\ x_5 \\ x_6 \\ x_7}
$$

\rule{\textwidth}{0.5ex}

{\large\bf Problem 1(c)}. Find the eigenvalues of $A$.

{\bf Answer}: $\lambda=0,-1,-1,-1,-1,-1,-7$

\bigskip

\ifsolutions
{\bf Solution 1(c)}.

Let $D=|A-\lambda I|$. Replace $-1-\lambda$ by symbol $u$. Then
$$D=\left|
\begin{array}{crrrrrrr}
-5+u & 1 & 1 & 1 & 1 & 1 & 1 \\
1 & u & 0 & 0 & 0 & 0 & 0 \\
1 & 0 & u & 0 & 0 & 0 & 0 \\
1 & 0 & 0 & u & 0 & 0 & 0 \\
1 & 0 & 0 & 0 & u & 0 & 0 \\
1 & 0 & 0 & 0 & 0 & u & 0 \\
1 & 0 & 0 & 0 & 0 & 0 & u
\end{array}
\right|
$$
Add each of rows $2,3,4,5,6$ to row $1$. Then $1+u$ is a common factor of row $1$ and
the determinant multiply rule implies
$$D=(1+u)\left|
\begin{array}{rrrrrrrr}
1 & 1 & 1 & 1 & 1 & 1 & 1 \\
1 & u & 0 & 0 & 0 & 0 & 0 \\
1 & 0 & u & 0 & 0 & 0 & 0 \\
1 & 0 & 0 & u & 0 & 0 & 0 \\
1 & 0 & 0 & 0 & u & 0 & 0 \\
1 & 0 & 0 & 0 & 0 & u & 0 \\
1 & 0 & 0 & 0 & 0 & 0 & u
\end{array}
\right|
$$
Cofactor expansion along the last row, plus induction, gives the answer
$D=(u+1)(u-6)u^5=(-\lambda)(-\lambda-7)(-\lambda-1)^5$ with roots $\lambda=0,-7,-1,-1,-1,-1,-1$.
\fi

\rule{\textwidth}{0.5ex}

{\large\bf Problem 1(d)}. Find the eigenvectors of $A$.

\bigskip

\ifsolutions
{\bf Solution 1(d)}.

The root $\lambda=-1$ causes us to solve $(A+I)\vec{v}=\vec{0}$, which has coefficient matrix
$$ B= \left(
\begin{array}{rrrrrrrr}
-5 & 1 & 1 & 1 & 1 & 1 & 1 \\
1 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 & 0 & 0 & 0
\end{array}\right), \quad\mbox{with}\quad
\rref(B)= \left(
\begin{array}{rrrrrrrr}
1 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 1 & 1 & 1 & 1 & 1 & 1 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0
\end{array}\right).
$$
There are $2$ lead variables and $5$ free variables, hence $5$ basis vectors
$$
\colvector{0\\-1\\ 1 \\ 0\\ 0\\ 0\\ 0},
\colvector{0\\-1\\ 0 \\ 1\\ 0\\ 0\\ 0},
\colvector{0\\-1\\ 0 \\ 0\\ 1\\ 0\\ 0},
\colvector{0\\-1\\ 0 \\ 0\\ 0\\ 1\\ 0},
\colvector{0\\-1\\ 0 \\ 0\\ 0\\ 0\\ 1}.
$$

The eigenvector for $\lambda=0$ has all components equal to $1$. This fact is found from the
equation $(A-(0)I)\vec{v}=\vec{0}$, which has coefficient matrix $A$.

The eigenvector for $\lambda=-7$ has first component $-6$ and the remaining equal to $1$. The task begins with the equation $(A-(-7) I)\vec{v}=\vec{0}$, which has coefficient matrix
$$
\left(
\begin{array}{ccccccc}
7-6 & 1 & 1 & 1 & 1 & 1 & 1 \\
1 & 7-1 & 0 & 0 & 0 & 0 & 0 \\
1 & 0 & 7-1 & 0 & 0 & 0 & 0 \\
1 & 0 & 0 & 7-1 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 & 7-1 & 0 & 0 \\
1 & 0 & 0 & 0 & 0 & 7-1 & 0 \\
1 & 0 & 0 & 0 & 0 & 0 & 7-1
\end{array}
\right) =
\left(
\begin{array}{rrrrrrrr}
1 & 1 & 1 & 1 & 1 & 1 & 1 \\
1 & 6 & 0 & 0 & 0 & 0 & 0 \\
1 & 0 & 6 & 0 & 0 & 0 & 0 \\
1 & 0 & 0 & 6 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 & 6 & 0 & 0 \\
1 & 0 & 0 & 0 & 0 & 6 & 0 \\
1 & 0 & 0 & 0 & 0 & 0 & 6
\end{array}
\right)
$$
The eigenvectors for $\lambda=0$ and $\lambda=-7$ are respectively
$$
\colvector{1\\ 1\\ 1\\ 1\\ 1\\ 1\\ 1},\quad
\colvector{-6 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1}.
$$
\fi

\rule{\textwidth}{0.5ex}

{\large\bf Problem 1(e)}. Solve the differential equation $\dd
\frac{d\vec{u}}{dr} = A\vec{u}$ by the eigenanalysis method.

\bigskip

\rule{\textwidth}{0.5ex}


\begin{center}\Large\bf Three Methods for Solving $\frac{d}{dt}\vec{u}(t)=A\vec{u}(t)$\end{center}
\begin{itemize}
\item
{\bf Eigenanalysis Method}. The eigenpairs of matrix $A$ are
required. The matrix $A$ must be diagonalizable, meaning there are
$n$ eigenpairs
$(\lambda_1,\vec{v}_1),(\lambda_2,\vec{v}_2),\ldots,(\lambda_n,\vec{v}_n)$.
The main theorem says that the general solution of
$\vec{u}'=A\vec{u}$ is
$$\vec{u}(t) = c_1 e^{\lambda_1 t}\vec{v}_1 + c_2 e^{\lambda_2 t}\vec{v}_2
+\cdots +c_n e^{\lambda_n t}\vec{v}_n.$$

\item
{\bf Laplace's Method}. Solve the scalar equations by the Laplace
transform method. The {\bf resolvent method} automates this process:
$\vec{u}(t)=\L^{-1}\left((sI-A)^{-1}\right)\vec{u}(0)$.

\item
{\bf Cayley-Hamilton-Ziebur Method}. The solution $\vec{u}(t)$ is a
vector linear combination of the Euler solution atoms
$f_1,\ldots,f_n$ found from the roots of the characteristic equation
$|A-\lambda I|=0$. The vectors $\vec{d}_1,\ldots,\vec{d}_n$ in the
linear combination
$$\vec{u}(t) = f_1(t)\vec{d}_1 + f_2(t)\vec{u}_2 + \cdots +
f_n(t)\vec{d}_n$$
 are determined by the explicit formula
$$<\vec{d}_1\,|\,\vec{d}_2\,|\cdots|\,\vec{d}_n> =
<\vec{u}_0\,|\,A\vec{u}_0\,|\cdots|\,A^{n-1}\vec{u}_0>\left(W(0)^T\right)^{-1},$$
where $W(t)$ is the Wronskian matrix of atoms $f_1,\ldots,f_n$ and
$\vec{u}_0$ is the initial data.

\end{itemize}

\ifsolutions
{\bf Solution 1(e)}.

The eigenvectors corresponding to $\lambda=0,-7,-1,-1,-1,-1,-1$ are respectively
$$
\colvector{1\\ 1\\ 1\\ 1\\ 1\\ 1\\ 1},
\colvector{-6 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1},
\colvector{0\\-1\\ 1 \\ 0\\ 0\\ 0\\ 0},
\colvector{0\\-1\\ 0 \\ 1\\ 0\\ 0\\ 0},
\colvector{0\\-1\\ 0 \\ 0\\ 1\\ 0\\ 0},
\colvector{0\\-1\\ 0 \\ 0\\ 0\\ 1\\ 0},
\colvector{0\\-1\\ 0 \\ 0\\ 0\\ 0\\ 1}.
$$
The Eigenanalysis method then implies the solution $\vec{x}(r)$ of
$\dd\frac{d\vec{x}}{dr} = A\vec{x}$ is given for arbitrary constants
$c_1,\ldots,c_7$ by the expression
$$
c_1e^{0r}\minicolvector{1\\ 1\\ 1\\ 1\\ 1\\ 1\\ 1}+
c_2e^{-7r}\minicolvector{-6 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1}+
c_3e^{-r}\minicolvector{0\\-1\\ 1 \\ 0\\ 0\\ 0\\ 0}+
c_4e^{-r}\minicolvector{0\\-1\\ 0 \\ 1\\ 0\\ 0\\ 0}+
c_5e^{-r}\minicolvector{0\\-1\\ 0 \\ 0\\ 1\\ 0\\ 0}+
c_6e^{-r}\minicolvector{0\\-1\\ 0 \\ 0\\ 0\\ 1\\ 0}+
c_7e^{-r}\minicolvector{0\\-1\\ 0 \\ 0\\ 0\\ 0\\ 1}.
$$
\rule{\textwidth}{0.5ex}
\fi

\vfill\newpage
\rule{\textwidth}{0.5ex}
{\Large\bf Problem 2}. {\bf Home Heating}

Consider a typical home with attic,
basement and insulated main floor.

\medskip

 \begin{center}
\includegraphics[scale=1.0]{figures/homeHeating.png}
\end{center}

\subsection*{Heating Assumptions and Variables}

\begin{itemize}
\item
It is usual to surround the main living area with insulation, but
the attic area has walls and ceiling without insulation.

\item
The walls and floor in the basement are insulated by earth.

\item
The basement ceiling is insulated by air space in the joists, a layer of
flooring on the main floor and a layer of drywall in the basement.
\end{itemize}

The changing temperatures in the three levels
is modeled by Newton's cooling law and the variables
\begin{center}\begin{tabular}{lcl}
$z(t)$ &=& Temperature in the attic,
\\
$y(t)$ &=& Temperature in the main living area,
\\
$x(t)$ &=& Temperature in the basement,
\\
$t$ &=& Time in hours.
\end{tabular}\end{center}

A typical mathematical model is the set of equations


$$
\begin{array}{lcl}
x' &=& \dd \frac{1}{2}(45-x) + \frac{1}{2}(y-x),\\[1.5ex]
y' &=& \dd
\frac{1}{2}(x-y)+\frac{1}{4}(35-y)+\frac{1}{4}(z-y)+20,\\[1.5ex]
z' &=& \dd \frac{1}{4}(y-z)+\frac{3}{4}(35-z).
\end{array}
$$
\rule{\textwidth}{0.5ex}
{\bf Problem 2(a)}. Formulate the system of differential equations as a matrix system
$\frac{d}{dt}\vec{u}(t) = A\vec{u}(t)+\vec{b}$. Show details.

{\bf Answer}.
$\dd\vec{u} = \colvector{x(t)\\ y(t)\\ z(t)},\quad \
\vec{b}= \colvectorC{\frac{45}{2} \\[1ex] 20 + \frac{35}{4} \\[1ex] \frac{105}{4} },
\quad
A=\left(\begin{array}{rrr}
-1 & \frac12 & 0\\[1.5ex]
\frac12 & -1 & \frac14\\[1.5ex]
0 & \frac14 & -1
\end{array}
\right)
$

\ifsolutions
\rule{\textwidth}{0.5ex}
{\bf Solution Details}.

Expand the right side of the system as follows.
$$
\begin{array}{lcl}
x' &=& \dd \frac{45}{2}-\frac{x}{2} + \frac{y}{2}-\frac{x}{2},\\[1.5ex]
y' &=& \dd
\frac{x}{2}-\frac{y}{2}+\frac{35}{4}-\frac{y}{4}+\frac{z}{4}-\frac{y}{4}+20\\[1.5ex]
z' &=& \dd \frac{y}{4}-\frac{z}{4}+\frac{105}{4}-\frac{3z}{4}.
\end{array}
$$
Then collect on the variables:
$$
\begin{array}{lcl}
x' &=& \dd -x + \frac{y}{2} + \frac{45}{2},\\[1.5ex]
y' &=& \dd \frac{x}{2}-y+\frac{z}{4}+20+\frac{35}{4}\\[1.5ex]
z' &=& \dd \frac{y}{4}-z+\frac{105}{4}.
\end{array}
$$
The right side of this system can be written as $A\vec{u}+\vec{b}$. Vector $\vec{b}$ is
obtained by formally setting $x=y=z=0$ on the right. This justifies the answer given.

The matrix $A$ has columns equal to the partial derivatives
$\partial_x,\partial_y,\partial_z$ of the right side of the scalar system.
This idea is important, because it allows the computation of matrix $A$
{\bf without any of the preceding details}.

\fi

\rule{\textwidth}{0.5ex}
{\bf Problem 2(b)}. The heating problem has
an {\bf equilibrium solution} $\vec{u}_p(t)$ which is a constant
vector of temperatures for the three floors. It is formally found by
setting $\frac{d}{dt}\vec{u}(t)=0$, and then $\vec{u}_p =
-A^{-1}\vec{b}$. Justify the algebra and explicitly find
$\vec{u}_p(t)$.

{\bf Answer 2(b)}. $\dd\vec{u}_p(t)=-A^{-1}\vec{b} = \colvectorC{\frac{620}{11} \\[1ex]
\frac{745}{11} \\[1ex] \frac{475}{11}}=\colvectorC{56.36 \\ 67.73 \\ 43.18}$.

% Vector(3, {(1) = 56.36363636, (2) = 67.72727273, (3) = 43.18181818})

\ifsolutions

{\bf Solution Details}.

The equation upon setting the derivative equal to zero becomes
$\vec{0}=A\vec{u} + \vec{b}$ which implies $A\vec{u}=-\vec{b}$ and finally
$\vec{u}= -A^{-1}\vec{b}$.

The calculation is done by technology. The maple code:
\begin{verbatim}
A:=<-1,1/2,0|1/2,-1,1/4|0,1/4,-1>^+; b:=<45/2,20+35/4,105/4>;
-A^(-1).b; evalf(%);
\end{verbatim}
The solution can also be obtained by hand from the augmented matrix of $A$ and $-\vec{b}$,
using the linear algebra toolkit of swap, combination and multiply.

\fi

\rule{\textwidth}{0.5ex}
{\bf Problem 2(c)}. The homogeneous problem is $\frac{d}{dt}\vec{u}(t) = A\vec{u}(t)$.
It can be solved by a variety of methods, three major methods enumerated below.
Choose a method and solve for $\vec{x}(t)$.

{\bf Answer 2(c)}: The homogenous scalar general solution is
$$
\begin{array}{lcl}
x_1(t) &=&\dd -\frac{1}{2} c_1 e^{-t} + 2c_2 e^{-at} + 2c_3 e^{-bt},\\[1.5ex]
x_2(t) &=&\dd -\sqrt{5}c_2e^{-at} + \sqrt{5}c_3 e^{-bt},\\[1.5ex]
x_3(t) &=&\dd c_1e^{-t} + c_2 e^{-at} + c_3e^{-bt}.
\end{array}
$$
%LinearAlgebra[CharacteristicPolynomial](A,r);solve(%,r);
%                      3      2   43     11
%                     r  + 3 r  + -- r + --
%                                 16     16
%                         1  (1/2)       1  (1/2)
%                -1, -1 - - 5     , -1 + - 5
%                         4              4
%a:=1+sqrt(5)/4;b:=1-sqrt(5)/4;c:=1;



\bigskip



\begin{center}\Large\bf Three Methods for Solving $\vec{u}'=A\vec{u}$\end{center}
\begin{itemize}
\item
{\bf Eigenanalysis Method}. Three eigenpairs of matrix $A$ are required. The
matrix $A$ must be diagonalizable, meaning there are 3 eigenpairs
$(\lambda_1,\vec{v}_1),(\lambda_2,\vec{v}_2),(\lambda_3,\vec{v}_3)$. The main
theorem says that the general solution of $\vec{u}'=A\vec{u}$ is
$$\vec{u}(t) = c_1 e^{\lambda_1 t}\vec{v}_1 + c_2 e^{\lambda_2 t}\vec{v}_2 +
c_3 e^{\lambda_3 t}\vec{v}_3.$$

\item
{\bf Laplace's Method}. Solve the scalar equations by the Laplace transform method.
The {\bf resolvent method} automates this process:
$\vec{u}(t)=\L^{-1}\left((sI-A)^{-1}\right)\vec{u}(0)$.

\item
{\bf Cayley-Hamilton-Ziebur Method}. The solution $\vec{u}(t)$ is a
vector linear combination of the Euler solution atoms found from the
roots of the characteristic equation $|A-\lambda I|=0$. The vectors
in the linear combination are determined by an explicit formula.
\end{itemize}


\ifsolutions
\subsubsection*{Solution Details for Problem 2(c)}.

The roots of the characteristic polynomial are used in all three methods. This is
the polynomial equation $|A-\lambda I|=0$, having $n$ roots real and complex, for
an $n\times n$ matrix $A$.

Subtract $\lambda$ from the diagonal of $A$ and form the determinant. Then
cofactor expansion on row 3 gives
$$
|A-\lambda I|=\left|\begin{array}{ccc}
-1-\lambda & \frac12 & 0\\[1.5ex]
\frac12 & -1 -\lambda& \frac14\\[1.5ex]
0 & \frac14 & -1-\lambda
\end{array}
\right| = (-1-\lambda)\left(-\frac{1}{16} + (-1-\lambda)^2-\frac14\right).
$$
The roots are $-1,-a,-b$ where $a=1+\sqrt5/4 = 1.56, b=1-\sqrt5/4=0.44$.
The three roots are distinct, real and negative.


\subsection*{Eigenanalysis Method}
 The eigenpairs must be found, in order to assemble the solution vector $\vec{u}(t)$.
 Technology can be used to find the answers, which are
% LinearAlgebra[Eigenvectors](A);
% Vector(3, {(1) = -1, (2) = -1-(1/4)*sqrt(5), (3) = -1+(1/4)*sqrt(5)})
% Matrix(3, 3, {(1, 1) = -1/2, (1, 2) = 2, (1, 3) = 2, (2, 1) = 0, (2, 2) = sqrt(5), (2, 3) = -sqrt(5), (3, 1) = 1, (3, 2) = 1, (3, 3) = 1})
$$
\left(-1,\minicolvector{-\frac12\\ 0 \\ 1}\right), \quad
\left(-a,\minicolvector{2\\ -\sqrt5 \\ 1}\right), \quad
\left(-b,\minicolvector{2\\ \sqrt5 \\ 1}\right).
$$
Without technology, there are three homogeneous problems to solve of the
form $B\vec{v}=\vec{0}$, for eigenvector $\vec{v}$. Enumerated, they are:

{\bf Case $\lambda=-1$}. Then $B=A+I=\left(\begin{array}{ccc}
0 & \frac12 & 0\\[1.5ex]
\frac12 & 0& \frac14\\[1.5ex]
0 & \frac14 & 0
\end{array}
\right)$


{\bf Case $\lambda=-a$}. Then $B=A+aI=\left(\begin{array}{ccc}
\frac{\sqrt5}{4} & \frac12 & 0\\[1.5ex]
\frac12 & \frac{\sqrt5}{4}& \frac14\\[1.5ex]
0 & \frac14 & \frac{\sqrt5}{4}
\end{array}
\right)$

{\bf Case $\lambda=-b$}. Then $B=A+bI=\left(\begin{array}{ccc}
-\frac{\sqrt5}{4} & \frac12 & 0\\[1.5ex]
\frac12 & -\frac{\sqrt5}{4}& \frac14\\[1.5ex]
0 & \frac14 & -\frac{\sqrt5}{4}
\end{array}
\right)$

In each case, the system $B\vec{v}=\vec{0}$ is solved using the last frame
algorithm (there is in each case one free variable). The eigenvector reported
is the partial derivative of the general solution on the invented symbol $t_1$,
which was assigned to the free variable.

\begin{center}\large\sf Application of the Theorem\end{center}
System $\vec{u}'=A\vec{u}$ is solved in the diagonalizable case
in terms of the eigenpairs of $A$, denoted as
$(\lambda_1,\vec{v}_1),(\lambda_2,\vec{v}_2),(\lambda_3,\vec{v}_3)$.
The solution of $\vec{u}'=A\vec{u}$  is given by the formula
$$\vec{u}(t) = c_1e^{\lambda_1 t}\vec{v}_1 + c_2e^{\lambda_2 t}\vec{v}_2 +
c_3e^{\lambda_3 t}\vec{v}_3.$$

In the present case, the solution is
$$\vec{u}(t) = c_1e^{-t}
\minicolvector{-\frac12\\ 0 \\ 1}+
c_2e^{-at}\minicolvector{2\\ -\sqrt5 \\ 1}+
c_3 e^{-bt}\minicolvector{2\\ \sqrt5 \\ 1}.
$$
Symbols $c_1,c_2,c_3$ in the solution are arbitrary constants, uniquely
determined by initial conditions.
In scalar form, the solution is
$$
\begin{array}{lcl}
x_1(t) &=&\dd -\frac{1}{2} c_1 e^{-t} + 2c_2 e^{-at} + 2c_3 e^{-bt},\\[1.5ex]
x_2(t) &=&\dd -\sqrt{5}c_2e^{-at} + \sqrt{5}c_3 e^{-bt},\\[1.5ex]
x_3(t) &=&\dd c_1e^{-t} + c_2 e^{-at} + c_3e^{-bt}.
\end{array}
$$

\subsection*{Laplace Transform Method}


The {\bf Laplace Method} for solving $\vec{u}'(t)=A\vec{u}(t)$ is based upon
transforming all differential equations into the frequency domain. Then time
variable $t$ no longer appears, being replaced by the frequency variable $s$.

The homogeneous system of differential equations is
$$
\begin{array}{lcl}
x' &=& \dd -x + \frac{y}{2},\\[1.5ex]
y' &=& \dd \frac{x}{2}-y+\frac{z}{4}+20\\[1.5ex]
z' &=& \dd \frac{y}{4}-z.
\end{array}
$$
Transforming to the $s$-domain uses the parts rule $\L(f'(t))=s\L(f(t)-f(0)$. Then
$$
\begin{array}{lcl}
s\L(x)-x(0) &=& \dd -\L(x) + \frac{1}{2}\L(y),\\[1.5ex]
s\L(y)-y(0) &=& \dd \frac{1}{2}\L(x)-\L(y)+\frac{1}{4}\L(z)+20\\[1.5ex]
s\L(z)-z(0) &=& \dd \frac{1}{4}\L(y)-\L(z).
\end{array}
$$
View these equations as linear algebraic equations for the symbols $\L(x),
\L(y),\L(z)$. Move terms left and right to re-write the scalar equations as a matrix system
$$
\left(\begin{array}{ccc}
s+1 & -\frac12 & 0 \\
\frac12 & s+1 & -\frac14 \\
0 & -\frac14 & s+1
\end{array}\right)
\colvector{\L(x) \\ \L(y) \\ \L(z)}
=
\colvector{x(0) \\ y(0) \\ z(0)}.
$$
The system is solved by inverting the coefficient matrix $C$ on the left, using the
adjugate formula $C^{-1}= \adj(C)/|C|$.
Write the answer  as
 $$C^{-1}(s)=\left(\begin{array}{ccc}
s+1 & -\frac12 & 0 \\
\frac12 & s+1 & -\frac14 \\
0 & -\frac14 & s+1
\end{array}\right)^{-1} = \frac{1}{\Delta(s)}
\left(\begin{array}{ccc}
s^2+2s+\frac{15}{16} & \frac{s+1}{2} & \frac18 \\[1.5ex]
\frac{s+1}{2} & (s+1)^2 & \frac{s+1}{4} \\[1.5ex]
\frac18 & \frac{s+1}{4} & s^2+2s+\frac34
\end{array}\right).
$$
%
Symbol $\Delta(s)=(s+1)(s+a)(s+b)$ is the determinant of $C(s)$.
%[[s^2+2*s+15/16, (1/2)*s+1/2, 1/8],
%[(1/2)*s+1/2, (s+1)^2, (1/4)*s+1/4],
%[1/8, (1/4)*s+1/4, s^2+2*s+3/4]]
Then
$$
\colvector{\L(x) \\ \L(y) \\ \L(z)}
=
C^{-1}(s)
\colvector{x(0) \\ y(0) \\ z(0)}.
$$
Backward table steps require solving nine equations of the form
$\L(f(t)) = \frac{p(s)}{\Delta(s)}$. A computer algebra system reduces the
effort, able to write
$C^{-1}(s)=\L(\Phi(t))$, using symbols  $f_1=e^{-t}, f_2=e^{-at}, f_3=e^{-bt}$, where
$$
\Phi(t) = \left(
\begin {array}{ccc}
  \frac15f_1+\frac25 f_2+ \frac25  f_3 &  \frac{1}{\sqrt5}(f_3- f_2 )
& \frac15 f_2+ \frac15 f_3- \frac25 f_1
\\
\noalign{\medskip}
 \frac{1}{\sqrt5} \left(  f_3- f_2 \right) & \frac12 f_2+ \frac12f_3 &
\frac{1}{2\sqrt5}
 \left(  f_3- f_2 \right)
\\
\noalign{\medskip}
\frac15f_2+ \frac15f_3-\frac25 f_1& \frac{1}{2\sqrt5} \left( f_3- f_2 \right)
& \frac45 f_1+ \frac{1}{10} f_2+ \frac{1}{10} f_3
\end {array}
 \right)
$$
Then
$\L(\vec{u}(t))=\L(\Phi(t)\vec{u}(0))$ implies by Lerch's cancelation law the formula
$$\colvector{x(t)\\[1.5ex]  y(t)\\[1.5ex] z(t)} = \left(
\begin {array}{ccc}
  \frac15f_1+\frac25 f_2+ \frac25  f_3 &  \frac{1}{\sqrt5}(f_3- f_2 )
& \frac15 f_2+ \frac15 f_3- \frac25 f_1
\\
\noalign{\medskip}
 \frac{1}{\sqrt5} \left(  f_3- f_2 \right) & \frac12 f_2+ \frac12f_3 &
\frac{1}{2\sqrt5}
 \left(  f_3- f_2 \right)
\\
\noalign{\medskip}
\frac15f_2+ \frac15f_3-\frac25 f_1& \frac{1}{2\sqrt5} \left( f_3- f_2 \right)
& \frac45 f_1+ \frac{1}{10} f_2+ \frac{1}{10} f_3
\end {array}
 \right)
 \colvector{x(0)\\[1.5ex]  y(0) \\[1.5ex]  z(0)}.
 $$


%[[(1/5)*f1+(2/5)*f3+(2/5)*f2, (1/5*(f3-f2))*sqrt(5), (1/5)*f3+(1/5)*f2-(2/5)*f1],
% [(1/5*(f3-f2))*sqrt(5), (1/2)*f3+(1/2)*f2, (1/10*(f3-f2))*sqrt(5)],
% [(1/5)*f3+(1/5)*f2-(2/5)*f1, (1/10*(f3-f2))*sqrt(5), (4/5)*f1+(1/10)*f3+(1/10)*f2]]

\bigskip

\begin{center} \large\sf The Resolvent Method \end{center}

The scalar system solved above is exactly
$$ (sI - A)\L(\vec{u}(t)) = \vec{u}(0),\quad I=\minithreebythree{1&0&0}{0&1&0}{0&0&1},\quad
A=\minithreebythree
{-1 & \frac12 & 0}
{\frac12 & -1 & \frac14}
{0 & \frac14 & -1}
,\quad
\vec{u}(t)=\minicolvector{x(t)\\ y(t)\\ z(t)}.$$

 The system $ (sI - A)\L(\vec{u}(t)) = \vec{u}(0)$ is called the {\bf resolvent equation}.
 The inverse of the coefficient matrix, $(sI-A)^{-1}$, is called the {\bf resolvent matrix},
 because $\L(\vec{u}(t))=(sI-A)^{-1}\vec{u}(0)$. If these statements make sense to you, then
 please use them to solve problems. Otherwise, ignore the information, and solve problems in
 the same manner as outlined earlier.

\begin{center} \large\sf Engineering and Laplace Transforms \end{center}

Both mechanical engineering and electrical engineering have rich support for Laplace theory.
Using Laplace theory has the advantage that many persons can help you through difficult times.
Independent persons prefer to choose the method from their own private toolbox.

\subsection*{Cayley-Hamilton-Ziebur Method}

The Ziebur Lemma implies that the solution of the system $\vec{u}'(t)=A\vec{u}(t)$ is given
by the formula
$$\vec{u}(t) = \vec{d}_1 e^{-t} + \vec{d}_2 e^{-at} + \vec{d}_3 e^{-bt}.$$


{\large\bf THEOREM}. Vectors $\vec{d}_1,\vec{d}_2,\vec{d}_3$ are uniquely determined by
initial condition $\vec{u}(0)$, which is a column vector of prescribed constants, by the
matrix equation.

$$
\langle\vec{d}_1|\vec{d}_2|\vec{d}_3\rangle
 = \langle\vec{u}(0)|A\vec{u}(0)|A^2\vec{u}(0)\rangle\, \left(W(0)^T\right)^{-1}
$$
Symbol $W(t)$ is the Wronskian matrix of the three Euler solution atoms. Notation
$\langle\vec{A}|\vec{B}|\vec{C}\rangle$ is the augmented matrix of the three columns
vectors $\vec{A},\vec{B},\vec{C}$.

\bigskip

{\Large\bf Illustration}. For the heating example, with
$a=1+\sqrt5/4 = 1.56, b=1-\sqrt5/4=0.44$, the Euler solution atoms are
$e^{-t},e^{-at},e^{-bt}$. The Wronskian matrix is
$$W(t)=\threebythree{e^{-t} & e^{-at} & e^{-bt}}{-e^{-t} & -ae^{-at} &
-be^{-bt}}{e^{-t} & a^2e^{-at} & b^2e^{-bt}}, \ \quad
W(0)=\threebythree{1&1&1}{-1&-a&-b}{1&a^2&b^2}.$$
$$
\mbox{Then~~} (W(0)^T)^{-1} =
\left( \begin {array}{ccc} -{\frac {11}{5}}&\frac85-\frac25\,\sqrt {5}&\frac85+
\frac25\,\sqrt {5}\\ \noalign{\medskip}-{\frac {32}{5}}&{\frac {16}{5}}-\frac25
\,\sqrt {5}&{\frac {16}{5}}+\frac25\,\sqrt {5}\\ \noalign{\medskip}-{
\frac {16}{5}}&\frac85&\frac85\end {array} \right).
$$
For initial state $\vec{u}(0)=\minicolvector{1\\ 0\\ 0}$,
$ \langle\vec{u}(0)|A\vec{u}(0)|A^2\vec{u}(0)\rangle\, =
\left(\begin {array}{ccc} 1&-1&5/4\\ \noalign{\medskip}0&1/2&-1
\\ \noalign{\medskip}0&0&1/8\end {array} \right)
$. Then
$$\vec{u}(t) = \langle\vec{d}_1|\vec{d}_2|\vec{d}_3\rangle
\colvector{e^{-t}\\ e^{-at} \\ e^{-bt}} =
\left(\begin {array}{ccc} 1&-1&5/4\\ \noalign{\medskip}0&1/2&-1
\\ \noalign{\medskip}0&0&1/8\end {array} \right)
\,
\left( \begin {array}{ccc} -{\frac {11}{5}}&\frac85-\frac25\,\sqrt {5}&\frac85+
\frac25\,\sqrt {5}\\ \noalign{\medskip}-{\frac {32}{5}}&{\frac {16}{5}}-\frac25
\,\sqrt {5}&{\frac {16}{5}}+\frac25\,\sqrt {5}\\ \noalign{\medskip}-{
\frac {16}{5}}&\frac85&\frac85\end {array} \right)
\colvector{e^{-t}\\ e^{-at} \\ e^{-bt}}$$

For general initial state $\vec{u}(0)=\minicolvector{c_1\\ c_2\\ c_3}$,
$$\dd
  \langle\vec{u}(0)|A\vec{u}(0)|A^2\vec{u}(0)\rangle\, =
  \left( \begin {array}{ccc} c_1&-c_1+\frac12c_2&\frac54\,c_1-c_2+\frac18c_3\\
  \noalign{\medskip}
c_2& \frac12c_1-c_2
+\frac14c_3&-c_1+{\frac {21}{16}c_2}-\frac12c_3
\\ \noalign{\medskip}
c_3&\frac14c_2-c_3&\frac18c_1-\frac12c_2+{\frac {17}{16}c_3}\end {array} \right)
.
$$
Then $\vec{u}(t)$ is this matrix times $\left( W(0)^T\right)^{-1}$
times the column vector of atoms $e^{-t},e^{-at},e^{-bt}$.

\subsubsection*{Details for the Theorem}

The idea for finding the three vectors is differentiation of Ziebur's equation,
two times, to get three equations

$$
\begin{array}{lcrrrrrr}
\vec{u}(t) &=& \vec{d}_1 e^{-t} &+& \vec{d}_2 e^{-at} &+& \vec{d}_3 e^{-bt}, \\
\vec{u}'(t) &=& -\vec{d}_1 e^{-t} &-& a \vec{d}_2 e^{-at} &-& b \vec{d}_3 e^{-bt}, \\
\vec{u}''(t) &=& \vec{d}_1 e^{-t} &+& a^2\vec{d}_2 e^{-at} &+& b^2 \vec{d}_3 e^{-bt}.
\end{array}
$$
Identities $\vec{u}'(t)=A\vec{u}(t)$ and $\vec{u}''(t)=A\vec{u}'(t)=AA\vec{u}(t)$ imply
that the left sides are simplified to
$$
\begin{array}{rcrrrrrr}
\vec{u}(t) &=& \vec{d}_1 e^{-t} &+& \vec{d}_2 e^{-at}& +& \vec{d}_3 e^{-bt}, \\
A\vec{u}(t) &=& -\vec{d}_1 e^{-t} &-& a \vec{d}_2 e^{-at} &-& b \vec{d}_3 e^{-bt}, \\
A^2\vec{u}(t) &=& \vec{d}_1 e^{-t} &+& a^2\vec{d}_2 e^{-at} &+& b^2 \vec{d}_3 e^{-bt}.
\end{array}
$$
The critical idea is to substitute $t=0$, which because of $e^0=1$ gives the following
three equations for unknowns $\vec{d}_1,\vec{d}_2,\vec{d}_3$:
$$
\begin{array}{rcrrrrrr}
\vec{u}(0) &=& \vec{d}_1 &+& \vec{d}_2 & +& \vec{d}_3, \\
A\vec{u}(0) &=& -\vec{d}_1 &-& a \vec{d}_2 &-& b \vec{d}_3, \\
A^2\vec{u}(0) &=& \vec{d}_1 &+& a^2\vec{d}_2 &+& b^2 \vec{d}_3.
\end{array}
$$
How to solve these vector equations for unknowns $\vec{d}_1,\vec{d}_2,\vec{d}_3$?
To begin,  solve the scalar system
$$
\threebythree{1&1&1}{-1&-a&-b}{1&a^2&b^2}\colvector{x\\ y\\ z} = \colvector{b_1\\ b_2\\ b_3}
$$
where variables $x,y,z$ are the first components of $\vec{d}_1,\vec{d}_2,\vec{d}_3$, and
similarly, $b_1,b_2,b_3$ are the first components of vector
$\vec{u}(0),A\vec{u}(0),A^2\vec{u}(0)$:
\begin{center}
\hspace*{4pc}
\begin{tabular}{l}
$x=\vec{d}_1\cdot \vec{v}$, \\
$y=\vec{d}_2\cdot \vec{v}$,\\
$z=\vec{d}_3\cdot \vec{v}$,
\end{tabular}
\hfill
\begin{tabular}{l}
$b_1=\vec{u}(0)\cdot \vec{v}$, \\
$b_2=A\vec{u}(0)\cdot \vec{v}$,\\
$b_3=A^2\vec{u}(0)\cdot \vec{v}$,
\end{tabular}
\hfill
$\vec{v}=\minicolvector{1\\ 0\\ 0}$.
\hspace*{4pc}
\end{center}
The equations also apply to find the second components, using
$\vec{v}=\minicolvector{0\\ 1\\ 0}$, then the third components
using $\vec{v}=\minicolvector{0\\ 0\\ 1}$. The three systems of
equations can be written as one huge matrix equation:
$$
\threebythree{1&1&1}{-1&-a&-b}{1&a^2&b^2}\,
\langle\vec{d}_1|\vec{d}_2|\vec{d}_3\rangle^T
 = \langle\vec{u}(0)|A\vec{u}(0)|A^2\vec{u}(0)\rangle^T
$$
Taking the transpose across the equation gives
$$
\langle\vec{d}_1|\vec{d}_2|\vec{d}_3\rangle\,\threebythree{1&1&1}{-1&-a&-b}{1&a^2&b^2}^T
 = \langle\vec{u}(0)|A\vec{u}(0)|A^2\vec{u}(0)\rangle
$$
Finally, invert the matrix $W(0)^T$ and multiply across the equation on the right to obtain
$$
\langle\vec{d}_1|\vec{d}_2|\vec{d}_3\rangle
 = \langle\vec{u}(0)|A\vec{u}(0)|A^2\vec{u}(0)\rangle\,
 \left(\threebythree{1&1&1}{-1&-a&-b}{1&a^2&b^2}^T\right)^{-1}
$$
This is exactly the equation reported in the theorem,
$$
\langle\vec{d}_1|\vec{d}_2|\vec{d}_3\rangle
 = \langle\vec{u}(0)|A\vec{u}(0)|A^2\vec{u}(0)\rangle\,
 \left(W(0)^T\right)^{-1}
$$
It has been observed that if $f_1=e^{-t}, f_2=e^{-at}, f_3=e^{-bt}$
are replaced by a new basis of solutions such that $W(0)=I$, then
$\vec{d}_1=\vec{u}(0),\vec{d}_2=A\vec{u}(0),\vec{d}_3=A^2\vec{u}(0)$.
The resulting solution in this case is
$$\vec{u}(t) = f_1(t)\vec{u}(0) + f_2(t)A\vec{u}(0) +
f_3(t)A^2\vec{u}(0).$$

\bigskip

\begin{center}\Large\bf Which Method is the Best?\end{center}

The eigenanalysis method seems to be the best, because it is a method for
simplifying coordinates, hence a shorter answer. Except in the case of
complex roots. Or in the case that the matrix $A$ fails to be diagonalizable.
In practice, the method used to solve the equation $\vec{u}'(t) = A\vec{u}(t)$
has to be tuned to the expected application. Dynamical systems are an important example.
For dynamical systems, the actual solution is less important than its formula, which
is a linear combination of Euler solution atoms, according to Cayley-Hamilton-Ziebur.

Laplace theory provides a simple formula for the solution of $\vec{u}'(t) = A\vec{u}(t)$.
It has the form
$$\vec{u}(t) = \Phi(t)\vec{u}(0).$$
The matrix $\Phi(t)$ in the Laplace formula is computed from
$\L(\Phi(t)) = (sI-A)^{-1}$. Although this computation is nontrivial
by hand, computer algebra system automation is possible.

Matrix $\Phi(t)$ is called the {\bf Exponential Matrix}, denoted by $e^{At}$. Computer
algebra system {\tt Maple} computes $\Phi(t)$ by the command
\fbox{\tt LinearAlgebra[MatrixExponential](A,t)}. Both {\tt Matlab} and {\tt Mathematica}
support symbolic computation of the matrix exponential $\Phi(t)$.
\fi


\end{document}
\endinput
