% Hill5.tex

\section{Eigenvectors and Eigenvalues}

\subsection{A Brief Introduction to Determinants}

\headq{Note} If you have already worked through Chapter 2, you
should skip this section and go straight on to Section 5.2.

\noindent If a matrix $A$ is known to Maple, the determinant of $A$ is obtained
by the function

\verb$ > det(A);$

\noindent Some examples are given below.

\begin{verbatim}
# Hill, Section 5.1, Determinants.
> with(linalg):
> A:= matrix([[a11,a12],
>             [a21,a22]]):
> det(A);                 # See Hill Section 5.1 Definition 5.3, p316
                               a11 a22 - a12 a21

# The formula for the determinant of a 3 by 3 matrix can also be displayed.
> A:= matrix( [[a11,a12,a13],
>              [a21,a22,a23],
>              [a31,a32,a33]]):
> det(A):  # This is left to the reader to run.
\end{verbatim}

The following is a definition by induction on the
dimension of the matrix. The determinant of $2\times 2$
(and of $3\times 3$) matrices has been given in the Maple code above.
\smallskip

\headq{Definitions 5.6.} Let $A_{ij}$ be the $(n-1)\times (n-1)$ submatrix
  obtained by deleting the $i$-th row and $j$-th column of the matrix $A$. Then
\begin{enumerate}\item[]
  \begin{enumerate}
  \setcounter{enumii}{1}
  \item the ($i,j$)'th {\bf cofactor} of $A$ is \qquad
    $C_{ij} = (-1)^{i+j} \det(A_{ij});$
  \item the {\bf determinant} of $A$ is defined by ``expansion by cofactors
    along the first row'':
      $$ \det(A) = a_{11}C_{11} + a_{12}C_{12} +\cdots + a_{1n}C_{1n}. $$
  \end{enumerate}
\end{enumerate}

\setcounter{theorem}{5}
\begin{theorem}[d] Let $A$ be an $n\times n$ matrix.
  Then $\det(A)$ can be computed using expansion by cofactors along {\em any}
  row:
    $$ \det(A) = a_{i1}C_{i1} + a_{i2}C_{i2} +\cdots + a_{in}C_{in} .$$
\end{theorem}
Example:
\begin{verbatim}
# Expansion of a 3 by 3 determinant along its third row
> with(linalg):
> A:= matrix(3,3):
> A31:= submatrix(A,[1,2],[2,3]):
> A32:= submatrix(A,[1,2],[1,3]):
> A33:= submatrix(A,[1,2],[1,2]):
> simplify( A[3,1]*det(A31)-A[3,2]*det(A32)+A[3,3]*det(A33) - det(A) );
                                       0
\end{verbatim}

\goodbreak
\noindent Other useful facts about determinants of matrices are:
\begin{verbatim}
> with(linalg):
# Hill Section 5.1 Theorem 5.7 p321
# If A is an n by n lower triangular matrix, then
# det(A) is the product of the diagonal entries of A.  For example
> L:= matrix( [[a11,  0,  0],
>              [a21,a22,  0],
>              [a31,a32,a33]]):
> det(L);
                                  a11 a22 a33

# Hill Section 2.3 Theorem 2.25 p106
# Let A be an n by n matrix.    Then det(transpose(A)) = det(A).  For example
> A:= matrix(3,3):
> simplify(det(transpose(A))-det(A));
                                       0

# COROLLARY Any true statement about det(A) is still true if the words `row'
# and `column' are interchanged wherever they occur.  Thus theorem 5.7 yields
> U:= matrix( [[a11,a12,a13],
>              [  0,a22,a23],
>              [  0,  0,a33]]):
> det(U);
                                  a11 a22 a33

# THEOREM. det(A)=0 if and only if the rows of A are linearly dependent
#    (This includes theorems 2.7, 2.12 of Hill 2.2).    For example
> A0:= matrix( [[a11,a12],
>               [a11,a12]]):
> det(A0);
                                       0

# Hill Section 5.1 Theorem 5.9 p323
# Let A and B be n by n matrices.    Then det( evalm(A&*B) ) = det(A)*det(B).
> B:= matrix(3,3):
> simplify( det(evalm(A&*B)) - det(A)*det(B) );
                                       0
\end{verbatim}
The last two theorems yield the following useful facts about matrices.

\setcounter{theorem}{0}
\begin{theorem} The $n\times n$ matrix $A$ is invertible
  if and only if $\det(A)\ne 0$.
\end{theorem}
\headq{COROLLARY 5.10} {\em If $A$ is invertible,
  then\qquad $\det(A^{-1})=\det(A)^{-1}$.}

\subsection{Eigenvalues and Eigenvectors}
\headq{Definition 5.11} Let $A$ be an $n\times n$ matrix,
$\b x\in \DoubleR^n$ with $\b x\ne \b 0$, and $\lambda$ a number such that
  $$ A \b x = \lambda \b x .$$
Then \b x is called an {\bf eigenvector} of $A$ and $\lambda$ is the
{\bf eigenvalue} corresponding to \b x (or $\lambda$ is called an
{\em eigenvalue} of $A$ and \b x is an {\em eigenvector} corresponding to
$\lambda$).

\setcounter{theorem}{11}
\begin{theorem} The number $\lambda$ is an eigenvalue
  of the square matrix $A$ if and only if
    $$\det(\lambda I-A) = 0 .$$
  {\em This is called the {\bf characteristic equation} of $A$, and the
 polynomial $\det(\lambda I-A)$ is called the {\bf characteristic polynomial} of $A$.}
\end{theorem}

Eigenvalues and eigenvectors have many uses, some of which we will see
in later sections. We will begin with relevant Maple functions.
Maple has a set of linear algebra functions that allow you to determine, 
in those circumstances when it can be done in exact closed form solutions, 
the eigenvalues and eigenvectors of a square matrix. We will demonstrate this
using two matrices.

\label{eg4}
\begin{verbatim}
# Hill Section 5.2 Examples 4 and 5, pp. 328-330
> with(linalg):
> A:= matrix([ [5,-7,7],
>              [4,-3,4],
>              [4,-1,2] ]):
# The eigenvalues are the roots of the characteristic polynomial:
> charpoly(A,lambda);                         # This is  det( lambda * Id - A)
                                          2
                      (lambda - 5) (lambda  + lambda - 2)

> solve(",lambda);
                                    5, 1, -2

# Alternatively, use the function  eigenvals  in the  linalg package:
> eigenvals(A);
                                    5, 1, -2

> B:= matrix([ [0,-1],
>              [1, 0] ]):
> eigenvals(B);
                                     I, - I
\end{verbatim}
Notice that $A$ has three real eigenvalues and $B$ has two imaginary
eigenvalues (Maple uses {\tt I}\index{I=$\sqrt{-1}$} to denote the
imaginary number $\sqrt{-1}$).  These are the solutions to the
characteristic polynomial in each case.

\medskip\headq{EXERCISE} Find the eigenvalues of the following 
  three symmetric tridiagonal matrices:
\begin{verbatim}
  C1:= matrix([[1,1,0],[1,0,1],[0,1,1]]);
  C2:= matrix([[0,1,0],[1,0,1],[0,1,0]]);          
  C3:= matrix([[1,1,0],[1,0,1],[0,1,0]]);
\end{verbatim}

\medskip\headq{Notes}
\begin{enumerate}
\item Hill usually only discusses ``real vector spaces'' where the word
  ``number'' means ``real number''. However, when looking at
  eigenvalues and eigenvectors it is often important to look at complex
  vector spaces, where ``number'' now means ``complex number''. Maple's
  functions work the same when entries are complex as when entries are real.

\item The matrices {\tt C1}, {\tt C2}, {\tt C3} have cubic characteristic
  polynomials, and Cardan's formula for solving them can give lengthy answers.
  For example, all the solutions of {\tt C3} are real, but the
  number {\tt I} occurs in the answer. As stated earlier, exact solutions are
  sometimes both correct and confusing -- or horrible to look at -- whereas
  reasonably accurate floating-point numerical approximations can be easily
  comprehended. So try
  \begin{verbatim} 
  evalf(Eigenvals(C3));
  \end{verbatim} 
\end{enumerate}

\begin{theorem} Let $A$ be an $n\times n$ matrix. Then the following are
  equivalent:
  \begin{enumerate}\item[]
  \begin{enumerate}
  \item $\lambda$ is a (real) eigenvalue of $A$.
  \item There is a non-zero vector $\b x\in \DoubleR^n$ such that
    $A\b x=\lambda\b x$.
  \item The system of equations $(\lambda I-A)\b x = \b 0$ has a non-trivial
    solution.
  \item $\lambda$ is a (real) solution of the characteristic equation
    $\det(\lambda I-A)=0$.
  \end{enumerate}
  \end{enumerate}
\end{theorem}

\goodbreak
Maple also has an\quad {\tt eigenvects}\index{eigenvects}\quad function.
However, we will here find the eigenvectors using the functions you have learnt
in the earlier part of this document, and hence revising these.

\begin{theorem} The eigenvectors corresponding to a given
  eigenvalue $\lambda$ of a matrix $A$ are the non-zero vectors of the
  {\bf eigenspace} $W_\lambda = NS(A-\lambda I)$.
\end{theorem}

In the above example, you know that the eigenvalues of $A$ are 5, 1 and -2.
Thus, you are looking for the null space of the matrices $(A-5I)$, $(A-I)$ and
$(A+2I)$, where $I$ is the identity matrix.
\label{eg4.2}
\begin{verbatim}
# Continuing Hill Section 5.2 Example 4 (p328):
> Id:= diag(1,1,1):                        W5:= nullspace( evalm(A - 5*Id) );
                              W5 := {[ 1, 1, 1 ]}

> W1:= nullspace( evalm(A - Id) );
                              W1 := {[ 0, 1, 1 ]}

> W_2:= nullspace( evalm(A + 2*Id) );
                             W_2 := {[ -1, 0, 1 ]}

# For example,  W1  is a basis for  NS(A-Id),  and so is a set of eigenvectors.
# To CHECK that the vectors of this set satisfy the equation  A x = x,  find
> evalm(A &* W1[1]);  # The product of A with the first (and only) vector in W1
                                  [ 0, 1, 1 ]

# You should CHECK the other eigenvectors for yourself.
\end{verbatim}

\medskip\headq{EXERCISE} Find the eigenvectors of the matrices $C1$, $C2$, $C3$
  above. 

\medskip
\noindent [{\bf Note}\quad Prior to Maple V, the command {\tt nullspace(A)} or
  {\tt linsolve(A,b)} assumed that the entries of $A$ and $b$ were
  rational numbers, or similar. It did not want to deal with numbers like
  $\sqrt{2}$. Thus it was a harder problem to find the eigenvectors of the
  matrix $C2$ above.  If you have an early version of Maple, then one way is
  to find the rowspace of $C2-\lambda Id$ which is two dimensional.
  Then find a vector orthogonal to this with {\tt crossprod}. ]

\medskip
As well as the determinant, another useful function associated
with a matrix is the sum of its diagonal elements (called its {\bf trace}).
If $A=(a_{ij})$ is an $n\times n$ matrix then
  $${\rm trace}\,(A) = \sum_{i=1}^n a_{ii} = a_{11} + a_{22} +\cdots + a_{nn}.$$

\noindent Maple has a command {\tt trace}\index{trace} to find this.
\setcounter{theorem}{16}
\begin{theorem} Let $A$ be an $n\times n$ matrix with eigenvalues
  $\rowi{\lambda}{n}$. Then
  \begin{enumerate}\item[]
    \begin{enumerate}
    \item $\det(A)$ is the product
     $\lambda_1\times \lambda_2\times\cdots \lambda_n$ of the eigenvalues.
    \item ${\rm trace}\,(A)$ is the sum
     $\lambda_1 + \lambda_2 + \cdots \lambda_n$ of the eigenvalues.
    \end{enumerate}
  \end{enumerate}
\end{theorem}

\subsection{Diagonalization}
\setcounter{theorem}{19}
\begin{theorem} Suppose the $n\times n$ matrix $A$ has $n$ {\em linearly
  independent} eigenvectors $\rowi{\b v}{n}$ with associated eigenvalues
  $\rowi{\lambda}{n}$.  If $S$ is the matrix whose
  $i$'th column is $\b v_i$, then $S$ is invertible and
    $$ S^{-1} A S =\Lambda = \pmatrix{\lambda_1& & & \cr & \lambda_2 & & \cr
                                      & & \ddots & \cr & & & \lambda_n \cr }.$$
\end{theorem}

\headq{COROLLARY 5.22} If $A$, $\Lambda$ and $S$ are as above, then\qquad
  $S \Lambda S^{-1} = A$.
\medskip

\headq{Definition 5.23} An $n\times n$ matrix $A$ is {\bf diagonalizable} if
there is an invertible matrix $S$ such that $S^{-1} A S$ is a diagonal matrix.
\medskip

\headq{Note} At this point, it may be helpful to review Sections 4.1 and 4.7,
  noting that the matrix $S$ above is the transition matrix from the basis
  $\{\rowi{\b v}{n}\}$ to the standard basis, and $\Lambda$ is the matrix
  (relative to the basis $\{\rowi{\b v}{n}\}$\ ) of the induced function of $A$.

\begin{verbatim}
# Hill Section 5.3 Examples 1 and 2, p338.
> with(linalg):
> A:= matrix([ [2,3,0],
>              [4,3,0],
>              [0,0,6] ]):
> evects:= eigenvects(A);
    evects := [6, 2, {[ 1, 4/3, 0 ], [ 0, 0, 1 ]}], [-1, 1, {[ -1, 1, 0 ]}]

# The eigenvalues are  evects[1][1] and evects[2][1],  occurring
#   evects[1][2] and evects[2][2]  times respectively, and the corresponding
#   eigenvectors are the elements of the sets  evects[1][3] and evects[2][3].
#
# Since one of the eigenvalues is repeated, we force it to be the second one
> if evects[1][2] = 2 then    evects:= evects[2],evects[1]    fi;
    evects := [-1, 1, {[ -1, 1, 0 ]}], [6, 2, {[ 1, 4/3, 0 ], [ 0, 0, 1 ]}]
\end{verbatim}
\pagebreak
\begin{verbatim}
> S:= augment( evects[1][3][1],evects[2][3][1],evects[2][3][2] );
                                   [ -1   1   0 ]
                                   [            ]
                              S := [  1  4/3  0 ]
                                   [            ]
                                   [  0   0   1 ]

> evalsA:= evects[1][1], evects[2][1], evects[2][1]:
> equal( evalm( inverse(S)&*A&*S ), diag(evalsA) );
                                      true

# Hill Section 5.3 Example 3, p340.
> B:= matrix([ [a,1], [0,a] ]):
> evalsB:= eigenvals(B);
                                 evalsB := a, a

> nullspace( B-a*diag(1,1) );
                                   {[ 1, 0 ]}

# Since the dimension of this eigenspace is only 1,   B is not diagonalizable
\end{verbatim}

\goodbreak
\headq{THEOREMS 5.24 -- 5.26.} Let $A$ be an $n\times n$ matrix.
  Then the following are equivalent:
  \begin{enumerate}\item[]
    \begin{enumerate}
    \item $A$ is diagonalizable.
    \item $A$ has $n$ linearly independent eigenvectors.
    \item $\DoubleR^n$ has a basis consisting of eigenvectors of $A$.
    \item $A$ is similar to a diagonal matrix.
    \end{enumerate}
  \end{enumerate}

\goodbreak
\setcounter{theorem}{30}
\begin{theorem}[and 5.30] Let $A$ be an $n\times n$ matrix.\quad
  If the eigenvalues $\rowi{\lambda}{k}$ of $A$ are {\em distinct},
  then their associated eigenvectors are linearly independent.\quad
  In particular, if $A$ has $n$ real distinct eigenvalues,
  then $A$ is diagonalizable.
\end{theorem}
Thus the matrix $A$ in Example 4 of Hill Section 5.2 
(see page~\pageref{eg4} of these notes) is diagonalizable.

\smallskip
\headq{Note:} Example 1 of Hill Section 5.3 (above) shows that the converse
of theorem 5.30 is not true.
\medskip

\setcounter{theorem}{33}
\begin{theorem} Suppose $A$ is an $n\times n$ matrix, and
    $$ S^{-1} A S =\Lambda = \pmatrix{\lambda_1& & & \cr & \lambda_2 & & \cr
                                    & & \ddots & \cr & & & \lambda_n \cr }.
                   \qquad\qquad {\rm Then}\quad
      \Lambda^k = \pmatrix{\lambda_1^k& & & \cr & \lambda_2^k & & \cr
                                    & & \ddots & \cr & & & \lambda_n^k \cr }$$

  and\quad $A^k = S\Lambda^k S^{-1}$\quad for any positive integer $k$.\newline
  If all the eigenvalues $\rowi{\lambda}{n}$ are non-zero, this is also true
  for negative values of $k$.
\end{theorem}

\goodbreak
\begin{verbatim}
# An example of the power of a matrix (see Hill Section 5.3 Theorem 5.34)
#
# The k'th power of the diagonal form of the matrix A in example 1 (p338) is
> diag( evalsA[1]^k,evalsA[2]^k,evalsA[3]^k );
                               [     k         ]
                               [ (-1)    0   0 ]
                               [               ]
                               [         k     ]
                               [   0    6    0 ]
                               [               ]
                               [             k ]
                               [   0     0  6  ]

# and the  k'th  power of A is
> evalm( S &* " &* inverse(S) );
               [          k        k             k        k     ]
               [  4/7 (-1)  + 3/7 6    - 3/7 (-1)  + 3/7 6    0 ]
               [                                                ]
               [           k        k           k        k      ]
               [ - 4/7 (-1)  + 4/7 6    3/7 (-1)  + 4/7 6     0 ]
               [                                                ]
               [                                              k ]
               [           0                     0           6  ]
\end{verbatim}

\subsection{Symmetric Matrices}
\setcounter{theorem}{38}
\begin{theorem} Suppose $A$ is an $n\times n$ symmetric matrix. Then
\begin{enumerate}\item[]
  \begin{enumerate}
  \item all of the eigenvalues of $A$ are real;
  \item $A$ has $n$ linearly independent eigenvectors
    (and so $A$ is diagonalizable).
  \end{enumerate}
\end{enumerate}
\end{theorem}
\headq{Note:} The proof of this theorem uses not only the {\em idea} of complex
eigenvalues and eigenvectors (as mentioned in Section 5.2) but also a complex
{\bf inner product} $\b u.\b v = \b u^T\bar{\b v}$ for $\DoubleC^n$, where
$\bar{\b v}$ is the complex conjugate of \b v.  Maple's {\tt dotprod} function
uses this when its arguments are complex vectors, although it has a more general
function {\tt innerprod}\index{innerprod} used with real vectors and matrices.

\begin{verbatim}
> with(linalg):
> u:= vector([I,1+I]):                            v:= vector([-2*I,1+I]):
> simplify(dotprod(u,v)),                         innerprod(u,v);
                                   0, 2 + 2 I
\end{verbatim}

\setcounter{theorem}{42}
\begin{theorem}[and 5.45] Let $A$ be an $n\times n$ symmetric matrix. Then 
  \begin{enumerate} \begin{enumerate}
  \item eigenvectors of $A$ associated with different eigenvalues
    are orthogonal (and so linearly independent);
  \item $A$ has $n$ orthonormal eigenvectors, and so can be
    expressed as
      $$ A = Q \Lambda Q^T $$
    where $Q$ is orthogonal and $\Lambda$ is a diagonal matrix with the
    eigenvalues of $A$ on the diagonal.
  \end{enumerate} \end{enumerate}
\end{theorem}

\goodbreak
\noindent Finding the matrix $Q$ in theorem 5.45 depends, not only on these
theorems, but also on the Gram-Schmidt process (see Section 4.4).
This is illustrated in the following.
\begin{verbatim}
# Hill Section 5.4 Example 3, p352.
> with(linalg):
# Find an orthogonal matrix Q which diagonalizes the matrix A where
> A:= matrix([ [0, 2, 2],
>              [2, 0,-2],
>              [2,-2, 0] ]):
# Use the same ideas as in Example 1 of Hill Section 5.3:
> evects:= eigenvects(A);
     evects := [2, 2, {[ 1, 1, 0 ], [ 1, 0, 1 ]}], [-4, 1, {[ -1, 1, 1 ]}]

# One of the eigenvalues is repeated:  ensure that it is the second one
> if evects[1][2] = 2 then    evects:= evects[2],evects[1]    fi;
     evects := [-4, 1, {[ -1, 1, 1 ]}], [2, 2, {[ 1, 1, 0 ], [ 1, 0, 1 ]}]

#
> u1:= evects[1][3][1]:       u2:= evects[2][3][1]:      u3:= evects[2][3][2]:
> print( dotprod(u1,u2), dotprod(u1,u3), dotprod(u2,u3) );
                                    0, 0, 1

> v1:= evalm( u1/norm(u1,2) ):
# Since  u2, u3  are not orthogonal, use Gram-Schmidt:
> GS:= GramSchmidt( [u2,u3] ):
> v2:= evalm( GS[1]/norm(GS[1],2) ):    v3:= evalm( GS[2]/norm(GS[2],2) ):
> Q:= augment( v1,v2,v3 );
                        [        1/2       1/2        1/2  ]
                        [ - 1/3 3     1/2 2      1/6 6     ]
                        [                                  ]
                        [       1/2        1/2         1/2 ]
                   Q := [  1/3 3      1/2 2     - 1/6 6    ]
                        [                                  ]
                        [       1/2                   1/2  ]
                        [  1/3 3          0      1/3 6     ]

> orthog(Q);
                                      true

> evals:= evects[1][1], evects[2][1], evects[2][1]:
> equal( evalm( transpose(Q)&*A&*Q ), diag(evals) );
                                      true
\end{verbatim}

\goodbreak
\subsection*{Rotations and Reflections}
Although these are not dealt with extensively by Hill, they do appear
as examples of linear transformations used in Computer Aided Design.
In Example 5 of Section 4.1, Hill shows that
$$ Q(\theta) = \pmatrix{\cos (\theta)&-\sin (\theta)\cr
                        \sin (\theta)&\cos (\theta)\cr}$$
is the matrix of a rotation of $\DoubleR^2$ through an angle of $\theta$.
In figure 4.2 of Section 4.1, Hill also illustrates the fact that
$$
\pmatrix{\cos(\theta)&-\sin(\theta)&0\cr\sin(\theta)&\cos(\theta)&0\cr0&0&1\cr},
\qquad
\pmatrix{\cos(\theta)&0&-\sin(\theta)\cr0&1&0\cr\sin(\theta)&0&\cos(\theta)\cr},
\qquad
\pmatrix{1&0&0\cr0&\cos(\theta)&-\sin(\theta)\cr0&\sin(\theta)&\cos(\theta)\cr}
$$
are the matrices of rotations through an angle of $\theta$ in the
$xy$, $xz$ and $yz$ planes of $\DoubleR^3$ respectively (i.e. rotations about
axes orthogonal to those planes).
\medskip

\headq{THEOREM} {\em The matrix $Q$ induced by any rotation through an angle
  of $\theta$ in $\DoubleR^3$ is similar to each of the above three matrices.}
\medskip

\headq{COROLLARY} {\em Let $Q$ be the matrix in the above theorem. Then
  \begin{enumerate}\item[]
    \begin{enumerate}
    \item $Q$ is orthogonal and $\det(Q)=1$;
    \item $\lambda=1$ is an eigenvalue of $Q$,\quad and the associated
      eigenspace is the axis of $Q$\newline
      (i.e. is orthogonal to the plane of the rotation);
    \item ${\rm trace}\,(Q)=1+2cos(\theta)$,\quad and so the angle of the
      rotation is given by
       $$ cos(\theta) = ({\rm trace}\,(Q)-1)/2. $$
    \end{enumerate}
  \end{enumerate}}

\headq{THEOREM} {\em Let $Q$ be an orthogonal $3\times 3$ matrix.\quad Then
  $\det(Q)=\pm 1$.\newline
  If $\det(Q)=1$, then $Q$ is the matrix induced by a rotation,
  and if $\det(Q)= -1$, then
  \begin{enumerate}\item[]
    \begin{enumerate}
    \item $\lambda=-1$ is an eigenvalue of $Q$, and $Q$ reflects any associated
      eigenvector in the plane $P$ orthogonal to it;
    \item $Q$ is the matrix induced by a linear transformation which either
      fixes every point of $P$ or is a rotation in $P$ through an angle
      $\theta$ where
       $$ cos(\theta) = ({\rm trace}\,(Q)+1)/2. $$
    \end{enumerate}
  \end{enumerate}}

\begin{verbatim}
# Consider the matrix Q produced by the housetrans procedure of Hill 4.6
> with(linalg):
> Q:= matrix([ [ -2/3,  -1/3,   2/3 ],
>              [ -1/3, 14/15,  2/15 ],
>              [  2/3,  2/15, 11/15 ] ]):
> orthog(Q), det(Q);
                                    true, -1

# So  Q  reflects any vector  u  orthogonal to the plane P where
> u:= nullspace( Q+diag(1,1,1) )[1]:
> P:= dotprod( u,vector([x,y,z]) ) = 0;
                             P := 5 x + y - 2 z = 0

# It also rotates  P  through an angle  theta  where
> theta:= arccos( (trace(Q)+1)/2 );
                                   theta := 0

# So (as expected)  Q  is simply a reflection in the plane P.
\end{verbatim}

\goodbreak
\headq{EXERCISES}
\begin{enumerate}
  \item Show that the product and inverses of rotations are rotations.
  \item Show that the inverse of a reflection is a reflection.
  \item Show that the product of reflections through two intersecting planes is
    a rotation.
  \item Let $\b u\in\DoubleR^3$ with $\Vert \b u\Vert=1$.  Show that
  $I-2\b u\b u^T$ is a reflection in the plane orthogonal to \b u.
\end{enumerate}

\subsection{Difference Equations}
Several applications of eigenvalues are given in Hill Sections 5.5 to 5.7.\quad
Maple code dealing with these and other applications is given in the supplement.

Maple has a function {\tt fibonacci}\index{fibonacci} in the {\tt linalg}
package, where {\tt fibonacci(2)} yields the matrix $A$ of Hill's equation 5.49.
Maple also has a function {\tt fibonacci} in the {\tt combinat} package,
which gives the terms of the Fibonacci sequence,
and {\tt rsolve}\index{rsolve}, which solves difference equations.
You should get help on these if you are interested.
\medskip

\headq{Definition.} A {\bf Markov matrix} is a square matrix
with nonnegative entries, each of whose columns adds to 1.

\setcounter{theorem}{55}
\begin{theorem} If $A$ is a Markov matrix and any power of $A$ has all
   {\em positive} entries, then
  \begin{enumerate}\item[]
    \begin{enumerate}
    \item $\lambda=1$ is an eigenvalue of $A$ with multiplicity 1;
    \item the other eigenvalues of $A$ have absolute value less than 1, and
      the corresponding eigenvectors do not represent realistic distributions;
    \item for any initial population $\b u_0$, the populations
      $\b u_k=A^k \b u_0$ approach an eigenvector
      associated with the eigenvalue $\lambda=1$ as $k\to\infty$;
    \item the sum of the coordinates of this limit is equal to
      the sum of the coordinates of $\b u_0$.
    \end{enumerate}
  \end{enumerate}
\end{theorem}

\noindent As an example, we return to the first example in Chapter 1.
\begin{verbatim}
# Hill Section 1.1 Example 4, p4  (and pp 356-358)
> with(linalg):
> T:= matrix([ [8/10, 1/10],
>              [2/10, 9/10] ]):
> ES:= eigenvects(T);
                ES := [1, 1, {[ 1, 2 ]}], [7/10, 1, {[ -1, 1 ]}]

# The eigenvector associated with the eigenvalue  1  yields
> evalm( 90000000 * ES[1][3][1] );
                            [ 90000000, 180000000 ]

# which agrees with the long-term population given in Chapter 1.
\end{verbatim}

\goodbreak
\subsection{Differential Equations}
Constant coefficient homogenous systems of ordinary differential equations 
are treated in Hill \S5.6, who gives an example of a predator-prey problem.
The general homogeneous system is
  $$ {d\b x\over dt} = A\b x. $$
If we know that $A$ is diagonalizable, then there is an invertible matrix $S$
such that $S^{-1} A S =\Lambda$ where $\Lambda$ is a diagonal matrix whose
diagonal entries are the eigenvalues of $A$. So if we substitute $\b x=S\b y$,
we have $\b y=S^{-1} \b x$ and so
  $$ {d\b y\over dt} = S^{-1}{d\b x\over dt} = S^{-1}A\b x
                     = S^{-1} AS\b y = \Lambda \b y. $$
It is easy to see that, if the eigenvalues of $A$ are $\rowi{\lambda}{n}$, then
the solution of this is
\begin{eqnarray*}
  \b y &=& \pmatrix{ y_1\cr \vdots \cr y_n } 
       = \pmatrix{ c_1 e^{\lambda_1 t}\cr \vdots \cr c_n e^{\lambda_n t} } \\
  {\rm and~so}\qquad\qquad \b x = S\b y &=& c_1 e^{\lambda_1t}\b v_1
                 + c_2 e^{\lambda_2t}\b v_2 +\cdots + c_n e^{\lambda_nt}\b v_n
\end{eqnarray*}
where $\b v_i$ is the $i$'th column of $S$ (i.e. an eigenvector
associated with the eigenvalue $\lambda_i$).

\setcounter{theorem}{67}
\begin{theorem} Let $A$ be a diagonalizable matrix. Then the solutions of the
  differential equation\quad ${d\b x\over dt}=A\b x$\quad consist of the vector
  space spanned by the functions $e^{\lambda_i t}\b v_i$ where $\lambda_i$ is an
  eigenvalue of $A$ and $\b v_i$ is an eigenvector associated with $\lambda_i$.
\end{theorem}

For example, the following system (given in Hill Example 3 p366)
\begin{eqnarray*}
  {du\over dt} &=& 5u - 7v + 7w \\
  {dv\over dt} &=& 4u - 3v + 4w \\
  {dw\over dt} &=& 4u - v + 2w
\end{eqnarray*}
can be solved by re-writing it as\quad ${d\b x\over dt}=A\b x$\quad where
\quad $\b x=[u,v,w]^T$\quad and $A$ is the matrix in Example~4 of Hill
Section 5.2 (see page~\pageref{eg4} of these notes). From the continuation of
that example (see page~\pageref{eg4.2} of these notes), we know that the
columns of the matrix
  $$ S = \pmatrix{ 1&0&-1\cr 1&1&0\cr 1&1&1\cr } $$
are eigenvectors of $A$ associated with the eigenvalues 5, 1 -2 respectively.
So the solution is
  $$\b x = c_1 e^{5t}\pmatrix{1\cr 1\cr 1\cr} + c_2 e^t\pmatrix{0\cr 1\cr 1\cr}
                                 + c_3 e^{-2t}\pmatrix{-1\cr 0\cr 1\cr} $$

\headq{Note} By analogy with the equation\quad ${dx\over dt}=ax$,\quad we may
  define the {\bf exponential} of the matrix\quad $A = S\Lambda S^{-1}$\quad
  to be\quad $e^{At} = SE S^{-1}$\quad where $E$ is a diagonal matrix whose
  diagonal entries are\quad
  $e^{\lambda_1t}, e^{\lambda_2t},\ldots e^{\lambda_nt}.$\quad With this
  definition, the above solution may be written\quad $\b x = e^{At}\b c$.

Here is a code fragment using Maple's function
{\tt exponential}\index{exponential} in its
{\tt linalg} package (check that it agrees with the answer given above).

\begin{verbatim}
# Hill Section 5.6 Example 3, p366.
> with(linalg):
> A:= matrix([ [5,-7,7],
>              [4,-3,4],
>              [4,-1,2] ]):
> evalm( exponential(t*A) &* vector([c1-c3,c1+c2,c1+c2+c3]) );
            [ exp(5 t) c1 - exp(- 2 t) c3, exp(5 t) c1 + exp(t) c2,

                exp(5 t) c1 + exp(t) c2 + exp(- 2 t) c3 ]
\end{verbatim}

Differential equations can also be solved using {\tt dsolve}\index{dsolve}
in the main part of Maple.

\smallskip
\headq{EXERCISE.} (Optional, harder Maple.) Let $n=2$. Use Maple to set up
  $n\times n$ matrices $A$ and $B$ whose entries are general functions of $t$.
  Define matrices $dA/dt$, $dB/dt$ to be the derivative with respect to $t$ of
  matrices $A$ and $B$ respectively, and $d(AB)/dt$ to be the derivative 
  with respect to $t$ of the product $AB$. Calculate\quad
  $d(AB)/dt - A(dB/dt)-(dA/dt)B$\quad where the terms denote matrix products. 
  \quad
  (Hint: first consider $n=1$. You don't need Maple then. What is the answer?)

\subsection{Quadratic Forms}
\headq{Note} If you are not familiar with the graphs of conics and you have
graphing facilities on your terminal, then you should use Maple to draw the
following:
\begin{verbatim}
# Hill Section 5.7 Example 1, p369:    Graphs of some conic sections
#
# The following is the circle    4*x^2 + 4*y^2 = 1
> plot( [cos(t)/2, sin(t)/2, t=-Pi..Pi] );

# The following is the ellipse    9*x^2 + 4*y^2 = 1
> plot( [cos(t)/3, sin(t)/2, t=-Pi..Pi] );

# The following is the right hand-side of the hyperbola    9*x^2 - 4*y^2 = 1
> plot( [sec(t)/3, tan(t)/2, t=-Pi/4..Pi/4] );

# A more familiar example is the following (rotated) hyperbola
> plot( [x, 1/x, x=-infinity..infinity] );
\end{verbatim}

\headq{Definition} A {\bf quadratic form} in the variables $\rowi{x}{n}$ is
\begin{eqnarray*}
  q_n(\b x) &=& \sum_{i=1}^n\sum_{j=1}^n a_{ij} x^ix^j + \sum_{i=1}^n b_i x^i \\
            = a_{11} x_1^2 +\cdots + a_{nn} x_n^2
             &+& (a_{12}+a_{21})x_1x_2 +\cdots + (a_{n-1,n}+a_{n,n-1})x_{n-1}x_n
                                              + b_1x_1 +\cdots + b_nx_n.
\end{eqnarray*}
The first (double) sum is called the {\bf principal part} of the quadratic form.
\quad The following piece of Maple code illustrates how the second sum can be
removed by moving the origin to the centre of the corresponding locus\quad
$q_n(\b x)=c$.
\begin{verbatim}
# Moving the origin to the centre  (h,k)  of the quadratic locus  Q2  where
> Q2:= a11*X1^2 + a22*X2^2 + 2*a12*X1*X2 + 2*b1*X1 + 2*b2*X2 - c:
> q2:= expand( subs(X1=x1-h,X2=x2-k, Q2) ):
#
\end{verbatim}
\goodbreak
\begin{verbatim}
# We want the coefficients of x1,x2 to be 0: collect q2 as a polynomial in x1
> qx1:= collect( q2,x1 ):                        # Also as a polynomial in x2:
> qx2:= collect( q2,x2 ):
> coeff1:= coeff(qx1,x1);
                 coeff1 := 2 a12 x2 + 2 b1 - 2 a12 k - 2 a11 h

> coeff2:= coeff(qx2,x2);
                coeff2 := - 2 a22 k - 2 a12 h + 2 b2 + 2 a12 x1

# These include the unwanted term  2 a12 x1 x2:  extract the constant term
> lprint();   solve( {coeff(coeff1,x2,0)=0,coeff(coeff2,x1,0)=0}, {h,k} );

                        a11 b2 - a12 b1      b1 a22 - a12 b2
                   {k = ---------------, h = ---------------}
                                      2                    2
                         a11 a22 - a12        a11 a22 - a12

# CHECK that these values give a quadratic in  x1,x2  where the coefficients
#    of x1,x2 are 0.    (This is not as easy as it seems!)
\end{verbatim}

\setcounter{theorem}{78}
\begin{theorem} The principal part of any quadratic form can be rewritten as
  \quad $\b x^T A\b x$\quad where $A$ is a square symmetric matrix.
\end{theorem}

\begin{verbatim}
> with(linalg):
# Hill Section 5.7 Example 5, p374.
# a)
> A:= matrix([ [  3 ,-2/3,-2/3],
>              [-2/3, 8/3,  0 ],
>              [-2/3,  0 ,10/3] ]):
# is the symmetric matrix for the quadratic form
> lprint();    X:= vector( [x,y,z] ):   qxyz:= evalm( transpose(X) &* A &* X );

                         2                            2         2
              qxyz := 3 x  - 4/3 x y - 4/3 x z + 8/3 y  + 10/3 z

# b) Diagonalize A   (See Examples 1,3 of Hill Sections 5.3, 5.4 respectively)
> evects:= eigenvects(A);
evects := [2, 1, {[ 2, 2, 1 ]}], [4, 1, {[ -2, 1, 2 ]}], [3, 1, {[ 1, -2, 2 ]}]

> u1:= evalm(evects[1][3][1]/norm(evects[1][3][1],2)):
> u2:= evalm(evects[2][3][1]/norm(evects[2][3][1],2)):
> u3:= evalm(evects[3][3][1]/norm(evects[3][3][1],2)):
> S:= augment( u1,u2,u3 );
                                 [ 2/3  -2/3   1/3 ]
                                 [                 ]
                            S := [ 2/3   1/3  -2/3 ]
                                 [                 ]
                                 [ 1/3   2/3   2/3 ]

> evals:= evects[1][1], evects[2][1], evects[3][1]:
> B:= evalm( transpose(S)&*A&*S ):                      equal( B,diag(evals) );
                                      true

# Thus if the surface  qxyz = 12  is rotated so that the x1-, y1-, z1-axes
#  point along  u1, u2, u3  respectively, it becomes the ellipsoid
> X:= vector( [x1,y1,z1] ):               evalm( transpose(X) &* B &* X ) = 12;
                               2       2       2
                           2 x1  + 4 y1  + 3 z1  = 12
\end{verbatim}

\subsection{Solving the Eigenvalue Problem Numerically}
Hill gives an introduction to various numerical methods for finding
eigenvalues and eigenvectors which are illustrated by the following Maple code.

\headq{(Note:} Some of the answers are more accurate than those given in Hill
because Maple is using higher precision.)

\begin{verbatim}
# Examples from Hill Section 5.8
> with(linalg):
# The following matrix is used by Hill in all the examples of Section 5.8
> A:= matrix([ [4,-2],  [3,-1] ]);
                                      [ 4  -2 ]
                                 A := [       ]
                                      [ 3  -1 ]

> EV:= eigenvects(A);
                 EV := [2, 1, {[ 1, 1 ]}], [1, 1, {[ 1, 3/2 ]}]

# Make sure the eigenvectors are in the same order as Hill lists them:
> if EV[1][1] = 1 then      EV:= EV[2],EV[1]      fi;
> v1:= map( evalf, evalm(EV[1][3][1]/norm(EV[1][3][1],2)) );
                       v1 := [ .7071067810, .7071067810 ]

> v2:= map( evalf, evalm(EV[2][3][1]/norm(EV[2][3][1],2)) );
                       v2 := [ .5547001960, .8320502943 ]

# Example 1 (p378)  Direct iteration
> direct_it:= proc(A,w,m,lambda) #  w = initial value,  m = no. of iterations.
>     local j,wj,Aw:
>     wj:= w:
>     for j from 1 to m do
>         Aw:= map( evalf, evalm(A&*wj) ):
>         lambda:= evalf( norm(Aw,2) ):
>         wj:= map( evalf, evalm(Aw/lambda) ):
>         print( eval(lambda),wj );  # we need  eval  because lambda is a name
>     od:
>     wj:
> end:
> direct_it( A,vector([1,0]),5,'lambda' ):
\end{verbatim}
\goodbreak
\begin{verbatim}
                   5.000000000, [ .8000000000, .6000000000 ]

                   2.690724809, [ .7432941464, .6689647318 ]

                   2.260640891, [ .7233555443, .6904757467 ]

                   2.115834767, [ .7148340257, .6992941553 ]

                   2.054850372, [ .7108779363, .7033154052 ]

# this is converging (slowly) to the eigenvalue 2  and eigenvector v1.
#
# Example 3 (p381)  Inverse iteration
> inverse_it:= proc(A,w,m,lambda)  #  w,m  as for direct iteration.
>     local wj,y,j:         wj:= w:
>     for j from 1 to m do
>         y:= linsolve( A,wj ):
>         lambda:= evalf( 1/norm(y,2) ):
>         wj:= map( evalf, evalm(lambda*y) ):
>         print( 1/eval(lambda),wj );      # Hill lists the INVERSE of lambda
>     od:
>     wj:
> end:
> inverse_it( A,vector([1,0]),4,'lambda' ):
                  1.581138830, [ -.3162277660, -.9486832980 ]

                  1.627882058, [ -.4856429314, -.8741572763 ]

                  1.199449559, [ -.5263546146, -.8502651465 ]

                  1.083785109, [ -.5417013338, -.8405710350 ]

# lambda,y  are converging to the eigenvalue 1  and eigenvector v2.
#
# Example 4 (p382)  Shifted inverse iteration - this is faster than the above.
#   First find a reasonable approximation to use as the shift
> w:= direct_it( A,vector([1,0]),3,'a' ):
                   5.000000000, [ .8000000000, .6000000000 ]

                   2.690724809, [ .7432941464, .6689647318 ]

                   2.260640891, [ .7233555443, .6904757467 ]

> inverse_it( A-diag(a,a),w,3,'lambda' ):
                  4.189493804, [ -.7102126934, -.7039871657 ]

                   3.903648935, [ .7077390322, .7064739640 ]
\end{verbatim}
\nopagebreak
\begin{verbatim}
                  3.850307450, [ -.7072370847, -.7069764533 ]
\end{verbatim}
\goodbreak
\begin{verbatim}
> print(a-lambda);       # This is an approximation to the eigenvalue  2
                                  2.000921372

# Example 5 (p383/384)  The QR Algorithm
#  -  the diagonal entries converge to the eigenvalues
> Ak:= copy(A):
> for k from 1 to 4 do
>     Rk:= map( evalf, qrfac2(Ak,'Qk') ):         #See Section 4.6 for  qrfac2
>     Ak:= map( evalf, evalm(Rk&*Qk) ):                 print(Ak);
> od:
                         [ 2.680000000  -4.760000000 ]
                         [                           ]
                         [ .2400000000   .3200000000 ]

                         [  2.259668508  -4.933701658 ]
                         [                            ]
                         [ .06629834253   .7403314917 ]

                         [  2.115675674  -4.974054054 ]
                         [                            ]
                         [ .02594594595   .8843243244 ]

                         [  2.054817672  -4.988408594 ]
                         [                            ]
                         [ .01159140306   .9451823232 ]

# The Shifted QR Algorithm  (starting with the above results)
> n:= 2:                      a:= Ak[n,n]:
> for k from 1 to 3 do
>     Rk:= map( evalf, qrfac2(Ak-diag(a,a),'Qk') ):
>     Ak:= map( evalf, evalm(Rk&*Qk+diag(a,a)) ):                 print(Ak);
> od:
                        [   2.002713756   -4.999455714 ]
                        [                              ]
                        [ .0005442846406   .9972862380 ]

                       [   2.000140677    -4.999971865 ]
                       [                               ]
                       [ .00002814086360   .9998593169 ]

                       [    2.000007305    -4.999998546 ]
                       [                                ]
                       [               -5               ]
                       [ .1462058736*10     .9999926904 ]
\end{verbatim}

\headq{EXERCISE} Write a Maple procedure to ``pre-process'' a matrix for these
iterative processes by reducing it to upper-Hessenberg form
(this is not used in the above examples as the matrix given by Hill is
$2\times 2$ whereas this pre-processing is only useful for matrices of higher
dimension).

\vfill
\eject

