% FILE antonMO.tex
% AR Chapter 8-in-6th ed., 9-in-7th ed.
\section{Applications}

\subsection{Application to Differential Equations ([J] \S3.3)}

Constant coefficient homogenous systems of ordinary differential equations 
are treated in this section.
The general homogeneous system is
  $$ {d\b y\over dt} = A\b y. $$
If we know that $A$ is diagonalizable, then there is an invertible matrix $P$
such that $P^{-1} A P =\Lambda$ where $\Lambda$ is a diagonal matrix whose
diagonal entries are the eigenvalues of $A$. So if we substitute $\b y=P\b u$,
we have $\b u=P^{-1} \b y$ and so
  $$ {d\b u\over dt} = P^{-1}{d\b y\over dt} = P^{-1}A\b y
                     = P^{-1} AP\b u = \Lambda \b u. $$
It is easy to see that, if the eigenvalues of $A$ are $\rowi{\lambda}{n}$, then
the solution of this is
\begin{eqnarray*}
  \b u &=& \pmatrix{ u_1\cr \vdots \cr u_n } 
       = \pmatrix{ c_1 e^{\lambda_1 t}\cr \vdots \cr c_n e^{\lambda_n t} } \\
  {\rm and~so}\qquad\qquad \b y = P\b u &=& 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 $P$ (i.e. an eigenvector
associated with the eigenvalue $\lambda_i$).

\theoremNN{\it
Let $A$ be a diagonalizable matrix. Then the solutions of the
  differential equation\quad ${d\b y\over dt}=A\b y$\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 the eigenvector associated with $\lambda_i$.
}% \end{theorem}
\theoremNNEnd

For example, the following system (given in 
[AR6] Example 2 p370; [AR7] p445)
\begin{eqnarray*}
  {{d{y_1}}\over{dt}} &=& y_1 + y_2  \\
  {{d{y_2}}\over{dt}} &=&  4 y_1 - 2 y_2
\end{eqnarray*}
can be solved by re-writing it as\quad ${d\b y\over dt}=A\b y$\quad where
\quad $\b y=[y_1,y_2]^T$\quad and $A$ is the 2 by 2 matrix 
with coefficients read off the right-hand side.
The eigenvalues $\lambda_1=2$, $\lambda_2=-3$ are calculated
in the usual way. The columns of
  $$ P = \pmatrix{ 1&-{1\over4}\cr 1&1\cr } $$
are eigenvectors of $A$ associated with the eigenvalues 2, -3 respectively.
So the solution is
  $$\b y = c_1 e^{2t}\pmatrix{1\cr 1\cr } 
         + c_2 e^{-3t}\pmatrix{-{1\over4}\cr 1\cr} .
$$

\headq{Note.} By analogy with the equation\quad ${dy\over dt}=ay$,\quad we may
  define the exponential of a matrix by
  $$\Psi(t)
   = (e^{\lambda_1t}\b v_1, e^{\lambda_2t}\b v_2,\ldots e^{\lambda_nt}\b v_n)
  \qquad {\rm exponential}(At) = (\Psi(0))^{-1}\Psi(t) .$$
Maple has a function 
{\tt exponential}\index{exponential}
in its {\tt linalg} package.
Here is a code fragment.\quad
(Check that it agrees with the answer given above.)

\begin{maple}
# based on 
# Example 2 - AR6 Section 8.1 p370; AR7 Section 9.1 p445
> with(linalg,matrix,exponential,vector):
> A:= matrix([ [1, 1],
               [4,-2] ]):
> evalm( exponential(t*A) &* vector([c1-c2/4,c1+c2]) );
     [ - 1/4 exp(- 3 t) c2 + exp(2 t) c1, exp(- 3 t) c2 + exp(2 t) c1 ]
\end{maple}

For solving differential equations one also has an option of using {\tt dsolve} 
within 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{Geometry of Linear Transformations from
$\DoubleR^2$ to $\DoubleR^2$}

There is room for use of Maple {\tt plot} to illustrate the
items in this Section, and a Maple Worksheet seems the most
appropriate way of presenting the code and the graphics.

\subsection{Least Squares Fitting to Data}

In this section, we consider the situation where\quad $M\b x=\b y$\quad is
inconsistent (i.e. has no solution). 
(In \cite{AR} they use $\b v$ where we use $\b x$.)
In this case, it is natural to
look for an $\b x$ which makes $M\b x$ as close as possible to $\b y$.

\smallskip
\headq{Definition.} The {\em least squares problem} for the inconsistent
  system\quad $M\b x=\b y$\quad is to find a vector ${\b x}_*$
  (called the {\em least squares solution} to the system) such that
  $\Vert\b y - M{\b x}_*\Vert$ is a minimum.

\smallskip
The simplest, and possibly the most widely used application is the fitting of
a straight line to data points.
Maple has a function 
{\tt leastsqrs}\index{leastsqrs}
in the {\tt linalg} package (Get {\tt help}
on it). If $M$ is a matrix and \b y is a vector, then {\tt leastsqrs(M,y)}
finds the least squares solution of the system $M\b x=\b y$. For example

\begin{maple}
# AR6 S8.3; AR7 S9.3 Least squares######################################
# Maple has a function in the linalg package     leastsqrs
# (Get help on it.)
# leastsqrs(M,b)   M matrix   ,  b vector
# finds least squares solution of Mx = b

# Example 1 - AR6 Section 8.3 p391; AR7 Section 9.3 p463
# Find the equation  y = a + b x  for the line of best fit for the points
#   (0,1), (1,3), (2,4), (3,4)
with(linalg,leastsqrs):
M:=array([[1,0],[1,1],[1,2],[1,3]]):       y:=array([1,3,4,4]):
xStar:=leastsqrs(M,y);
                          xStar :=  [ 3/2, 1 ]
# So the line of best fit is  y = x + 3/2
\end{maple}
\smallskip
\begin{maple}
#
# The same problem can be done with functions from the stats package.
# CAUTION stats package revised in Maple Vr3. First is for Vr1 and Vr2:
with(stats):
dat := array([[y,x],[1,0],[3,1],[4,2],[4,3]]):
v:=regression(dat,y=a+b*x);
#                   v := {a = 1.500000000, b = 1.0000000000}
\end{maple}
% Maple Vr2 ONLY
% assign(v); # This assigns the above values to a and b   
% statplot(dat,y=a+b*x);
\goodbreak
\begin{maple}
# Other sorts of least-squares fits are possible with the
# regression function.
# (Apologies for the obsolete units used in the text.
# Think of it as a piece of history: "Galileo's data set".)
# Example 3 - AR6 S8.3; AR7 S9.3. With time scaled by factor of 10.
dat := array([[y,x],[-0.18,1],[0.31,2],[1.03,3],[2.48,4],[3.73,5]]):
v:=regression(dat,y=a0+a1*x+a2*x^2);
#  v := {a0 = -.3980000000, a1 = .03471428572, a2 = .1607142857}
\end{maple}
%#  Run the above for answers, below for a plot.
% assign(v); statplot(dat,y=a0+a1*x+a2*x^2);
\medskip
\begin{maple}
#
# CAUTION stats package revised in Maple Vr3. Here we use Vr3:
with(stats): with(fit):
dat := [[0,1,2,3],[1,3,4,4]]:
fit[leastsquare[[x,y],y=a+b*x, {a,b}]](dat);
# Giving       {a = 3/2, b = 1}
\end{maple}
\smallskip
\begin{maple}
# Again, other sorts of least-squares fits are possible 
# Example 3 - AR6 S8.3; AR7 S9.3. With time scaled by factor of 10.
dat := [[1,2,3,4,5],[-0.18,0.31,1.03,2.48,3.73]]:
fit[leastsquare[[x,y],y=a0+a1*x+a2*x^2 ]](dat);
# Giving as before {a0 = -.3980000000, a1 = .03471428572, a2 = .1607142857}
\end{maple}
\smallskip

While we are considering regression, here is some Maple
deriving the widely-used formulae for least squares fitting of
a line $y=a + b x$ through data points $\{(x_i,y_i) \vert\ i=1..n\}$.
Define 
$$ {\bar x}={1\over{n}}\sum_{i=1}^n x_i, \qquad {\rm and}\qquad
{\bar y}={1\over{n}}\sum_{i=1}^n y_i .$$
The formulae we derive are:
$$a={{{\bar y} \sum_{i=1}^n x_i^2 -{\bar x}\sum_{i=1}^n x_i y_i}\over
{\sum_{i=1}^n x_i^2 - n{\bar x}^2}}, \qquad
b={{\sum_{i=1}^n x_i y_i - n {\bar x}{\bar y}}\over
{\sum_{i=1}^n x_i^2 - n{\bar x}^2}}.$$
(The familiar, partial, check is the identity: 
${\bar y} = a + b{\bar x}$.)
This linear regression problem is discussed in \cite{AR} % AR6 pp388-390,
and summarised in Theorem \ChMO.3.1.
We choose, for the sake of doing something a bit different
from the linear algebra in most of these notes, to use a
calculus approach.
\begin{maple}
# This approach is to minimize the sum of squares of
# the residuals, using calculus. 
# The quantities for which we solve are {a,b}. 
# Maple Vr1 and Vr2 solve commands do not like being asked
# solve({fa=0,fb=0},{a,b}) 
# but it is relatively straightforward to tell Maple that it is really 
# a simple set of linear equations to solve for a and b, 
# and to use linsolve.
> a:='a': b:='b': f:=sum((y[j]-a-b*x[j])^2,j=1..n):
  fa:=diff(f,a): fb:=diff(f,b):
  with(linalg,hessian,scalarmul,linsolve): 
  A:=hessian(f,[a,b]);
                        [                  n          ]
                        [                -----        ]
                        [                 \           ]
                        [      2 n         )   2 x[j] ]
                        [                 /           ]
                        [                -----        ]
                        [                j = 1        ]
                   A := [                             ]
                        [   n             n           ]
                        [ -----         -----         ]
                        [  \             \          2 ]
                        [   )   2 x[j]    )   2 x[j]  ]
                        [  /             /            ]
                        [ -----         -----         ]
                        [ j = 1         j = 1         ]



> fa0:=simplify(subs(a=0,b=0,fa));
                                   n
                                 -----
                                  \
                          fa0 :=   )   - 2 y[j]
                                  /
                                 -----
                                 j = 1

> fb0:=simplify(subs(a=0,b=0,fb));
                                 n
                               -----
                                \
                        fb0 :=   )   - 2 y[j] x[j]
                                /
                               -----
                               j = 1

> rhvec:=array([-fa0,-fb0]): linsolve(A,rhvec);
\end{maple}
Changing to Maple Vr2 output style to save paper, this gives:
% \b u &=& \pmatrix{ u_1\cr \vdots \cr u_n } 
%      = \pmatrix{ c_1 e^{\lambda_1 t}\cr \vdots \cr c_n e^{\lambda_n t} } \\
% {\rm and~so}\qquad\qquad \b y = P\b u &=& 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
\begin{eqnarray*}
[\qquad
&\ &{{
\bigl(\sum_{j=1}^n 2 x[j]\bigr)
\bigl(\sum_{j=1}^n -2 x[j] y[j]\bigr)
+ \bigl(\sum_{j=1}^n -2 y[j]\bigr)
\bigl(\sum_{j=1}^n 2 x[j]^2 y[j]\bigr)
}\over{
-\bigl(\sum_{j=1}^n 2 x[j]\bigr)^2
+ 2 \bigl(\sum_{j=1}^n 2 x[j]^2\bigr) n
}} ,\\
&\ &{{
\bigl(\sum_{j=1}^n 2 x[j]\bigr)
\bigl(\sum_{j=1}^n -2 y[j] y[j]\bigr)
-2 \bigl(\sum_{j=1}^n -2 y[j] x[j]\bigr) n
}\over{
-\bigl(\sum_{j=1}^n 2 x[j]\bigr)^2
+ 2 \bigl(\sum_{j=1}^n 2 x[j]^2\bigr) n
}}  \qquad ]
\end{eqnarray*}
Maple has been a bit obtuse about dividing out the 2s within 
the {\tt Sum} terms on both denom and numer.
\medskip

Regression, more general than the linear regression is sufficiently
important to be worth a few comments here.
The set-up is that we have data points, $\{(t_i,b_i)\}_{i=1}^m$
as before, and $n\leq m$ linearly independent functions $q_j$.
The hope is to approximate the data $b_i$ by
$\sum_{j=1}^n C_j q_j(t_i)$ and the least squares matrix is
$$M=\bigl( q_j(t_i) \bigr)_{m\times n} .$$


\medskip
\goodbreak
We now return to the general setting for linear least squares,
as discussed in Chapter \ChIPS,
but adapting to the notations in this section.
Approaches to the least-squares solution are given in the following
theorems.

\theoremNN{\it 
$\Vert\b y - M\b x\Vert$ is a minimum if and only
  if $\b y - M\b x\perp {\rm colSpace}(M)$.
} %\end{theorem}
\theoremNNEnd

%\begin{theorem} 
%\headq{THEOREM } 
\theoremNN{\it 
Let $M$ be an $m\times n$ matrix and $\b y\in \DoubleR^m$.
  A least-squares solution to the inconsistent system $M\b x=\b y$
  is a solution to the equation
    $$ M^T M {\bar\b x} = M^T \b y .$$
  If the columns of $M$ are linearly independent, then\quad
  ${\bar\b x} = (M^T M)^{-1} M^T \b y$  
} %\end{theorem}
\theoremNNEnd

$\b x$ satisfies the normal equations $M^T M\b x = M^T \b y$.
In fact, even Gaussian elimination presents problems here
as {\em the matrix $M^T M$ is almost always ill-conditioned}. 
See books 
on numerical linear algebra for details.
Some results related to the problem of calculating
the least squares solution of a system are repeated from
\ChIPS\ below.
The theory for least squares involves
a standard vector technique of projection.

\smallskip
\headq{Definition.} Suppose $W$ is a subspace of $\DoubleR^m$ and
  $\b y\in \DoubleR^m$.
  Then the {\em perpendicular projection} of $\b y$ into $W$ is the unique
  vector $\b p\in W$ such that $\b y - \b p\perp W$ and the function which
  maps each $\b y$ onto its perpendicular projection is called the
  {\em perpendicular projection function}.

\smallskip
%\headq{EXERCISE.} Show 
Earlier we saw that the perpendicular projection function is a linear
  transformation.

\theoremNN{\it 
  Assume $M$ has linearly independent columns.
  Then the matrix of the perpendicular projection function 
  into ${\rm colSpace}(M)$ is\quad
  $P=M (M^T M)^{-1} M^T$.
} %\end{theorem}
\theoremNNEnd

%\headq{EXERCISE} Show 
It can be shown that $P$ satisfies the following properties:

  \smallskip
  \qquad(a) $P^2 = P$.\qquad\qquad\qquad(b) $P^T = P$.

 \smallskip
 Conversely, any symmetric matrix with $P^2=P$ represents a projection function.

\begin{maple}
# Continuing the least squares example
> with(linalg,transpose,inverse,dotprod,col):
  Mtr:= transpose(M):           evalm( inverse(Mtr&*M) &* Mtr &* y - xStar );
                                    [ 0, 0 ]

> p:= evalm( M&*xStar):  # This is the projection of  y  into  colSpace(M).
> dotprod( y-p, col(M,1) );                 dotprod( y-p, col(M,2) );
                                       0

                                       0
\end{maple}

\noindent As discussed towards the end of \ChIPS, 
$QR$ decomposition is important in
in the practical solution of least squares.
%(and in an algorithm -- the $QR$-algorithm -- for numerical approximation of
%eigenvalues).

\theoremNN{\it 
  Suppose $M$ is an $m\times n$ matrix with linearly independent columns.  Then
  the $QR$-decomposition $M=QR$ replaces the equation $M^T M\b x=M^T \b y$ by
    $$ R \b x = Q^T \b y$$
} %\end{theorem}
\theoremNNEnd

\subsection{Approximation Problems; Fourier Series}

The use of Maple for this topic takes us further from the
{\tt linalg} package than seems appropriate for these notes.


\subsection{Quadratic Forms ([J] \S3.3)}

\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{maple}
# Here are some graphs of some conic sections
#
# The following is the circle    x^2 + y^2 = 1
> plot( [cos(t), sin(t), t=-Pi..Pi] );

# The following is the ellipse    9*x^2 + 4*y^2 = 36
# (which occurs in Example 4 - AR6 Section 8.6; AR7 Section 9.6)
> plot( [2*cos(t), 3*sin(t), t=-Pi..Pi] );

# The following is the right hand-side of the hyperbola 9*x^2 - 4*y^2 = 36
> plot( [2*sec(t), 3*tan(t), t=-Pi/4..Pi/4] );
# and an alternative parametrisation is as follows
> plot(  [2*cosh(t),3*sinh(t), t=-1 .. 1 ] );

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

\headq{Definition.} A {\em 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 {\em 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{maple}
# Moving the origin to the centre 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) ):
#
# 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{maple}

\theoremNN{\it 
  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}
\theoremNNEnd

\begin{maple}
# Example 5 - AR6 Section 8.5 p407; AR7 Section 9.5 p479
> with(linalg):
# a) The symmetric matrix for the quadratic form
#    q(x,y,z)  =  4 x^2 + 4 x y + 4 y^2 + 4 z^2 + 4 x z + 4 y z   is A where
> A:= matrix([ [ 4, 2, 2],
               [ 2, 4, 2],
               [ 2, 2, 4] ]):
# This can be checked as follows
> X:= vector( [x,y,z] ):                    evalm( transpose(X) &* A &* X );
                    2                      2              2
                 4 x  + 4 x y + 4 x z + 4 y  + 4 y z + 4 z

# b) Diagonalize A   
# Example 3 - AR6 Section 8.7 p426; AR7 Section 9.7 p498
> evals:= op(sort([eigenvals(A)]));
  evects:= eigenvects(A);

                              evals := 2, 2, 8

 evects := [2, 2, {[ -1, 1, 0 ], [ -1, 0, 1 ]}], [8, 1, {[ 1, 1, 1 ]}]

> if (nops(evects[1][3])=1) then j:=2 else j:=1 fi:
  v1:= evalm(evects[j][3][1]):
  u1:= evalm(v1/norm(v1,2)):
  v2:= evalm(evects[j][3][2]):
  u2:= GramSchmidt([v1,v2])[2]:
  u2:= evalm(u2/norm(u2,2)):
  if (j=1) then
    u3:= evalm(evects[2][3][1]/norm(evects[2][3][1],2)):
  else
    u3:= evalm(evects[1][3][1]/norm(evects[1][3][1],2)):
  fi:
  v2:= evalm(evects[1][3][2]):
  S:= augment( u1,u2,u3 );
                    [        1/2         1/2  1/2       1/2 ]
                    [ - 1/2 2     - 1/6 2    3     1/3 3    ]
                    [                                       ]
                    [       1/2          1/2  1/2       1/2 ]
               S := [  1/2 2      - 1/6 2    3     1/3 3    ]
                    [                                       ]
                    [                   1/2  1/2        1/2 ]
                    [      0       1/3 2    3      1/3 3    ]

> equal( evalm( transpose(S)&*A&*S ), diag(evals) );
                                      true

# Thus if the surface q(x,y,z) = C  is rotated so that the x'-, y'-, z'-axes
#  point along  u1, u2, u3  respectively, it becomes the ellipsoid
#
#                2 x'^2 + 2 y'^2 + 8 z'^2 = C .
\end{maple}

\subsection*{Positive definite matrices}

\subsection*{Minima, maxima and saddle points}

Given a function $f:\DoubleR^n \rightarrow\DoubleR$ its
first derivative is called the gradient
(Maple 
{\tt grad}\index{grad}) and its second derivative matrix
is called the hessian
(Maple {\tt hessian}\index{hessian}).
The {\it critical points} of $f$ are those points where the
gradient is zero.
See \cite{CH}.
The hessian determines a local quadratic approximation to
the function there.
Classifying whether a critical point is a
minimum, a maximum or a saddle point
depends on the whether the hessian there is
positive definite, negative definite or indefinite,
terms which we discuss soon.
Here is some Maple related to classifying
critical points, and we urge any reader trying it
to use Maple's {\tt plot3d} command to get a picture
of the surface
near each critical point.
\smallskip

\begin{maple}
with(linalg,grad,hessian,eigenvals):
x:='x': y:='y': f:=x^3 + x^4 - x^2*y^2 + (y^4)/2;
g:=grad(f,[x,y]);
cpts:=[solve({g[1]=0,g[2]=0},{x,y})];
h:=hessian(f,[x,y]);
for j from 1 to nops(cpts) do
 print(` ***************************************** `);
 h.j:=subs(cpts[j][1],cpts[j][2],op(h));
 eigenvals(h.j);
 f.j:=subs(cpts[j][1],cpts[j][2],f);
 od;
\end{maple}

Of course the second derivatives are not quite enough
for a complete analysis: zero eigenvalues of the hessian
require one to look at higher order terms.
The example above illustrates that.
\medskip

\par\noindent{\bf EXERCISE.}
Change the {\tt f} in the above to locate and classify the
critical points of $f$
\begin{maple}
f:=x^2*y^2  - 2*x -2*y;
\end{maple}

\medskip

\par\noindent{\bf EXERCISE.}
Consider more independent variables.
(You now don't have simple graphics like {\tt plot3d}
to help you.)
Adapt the code above and
change the {\tt f} in the above to
\begin{maple}
f:=2*x*y*z-4*x*z-2*y*z+x^2+y^2+z^2-2*x-4*y+4*z;
\end{maple}
Show that the critical points are
$(0,3,1)$, $(0,1,-1)$, $(1,2,0)$, $(2,1,1)$ and $(2,3,1)$
and classify them.

\medskip
\subsubsection*{
Tests for positive definiteness}

The following definition is in \cite{AR} \S\ChMO.5.

\headq{Definition.} A real symmetric matrix $A$ is said to be
{\it positive definite} if
${\b x}A{\b x} > 0$ for all ${\b x}\ne{\b 0}$.

% 6th ed. AR 8.5.2
\setcounter{theorem}{1}
\begin{theorem}
A real symmetric matrix $A$ is positive definite
if and only if
all the eigenvalues of $A$ are positive.
\end{theorem}
\smallskip

\headq{Definition.} Let $A$ be an $n\times n$ matrix.
The {\it principal submatrices} of $A$ are the submatrices
formed by the first $r$ rows and first $r$ columns of $A$
for $r=1,2,\ldots,n$.
\smallskip

% AR 8.5.3
\begin{theorem}
Each of the following tests is a
necessary and sufficient condition for the real symmetric matrix
$A$ to be positive definite: \newline
(a) $x^T A x>0$ for all nonzero vectors $x$.  \newline
(b) All the eigenvalues of $A$ satisfy $\lambda_i>0$.  \newline
(c) All the principal submatrices $A_k$ have positive
determinants.  \newline
(d) There is a matrix $R$ with independent columns such that
$A=R^T R$.
\end{theorem}
\smallskip

Condition (c) is called {\it Sylvester's condition}.
\smallskip

The Maple function
{\tt definite}\index{definite}
tests for positive definiteness of a matrix.
\smallskip


\medskip
\subsubsection*{
Semidefinite and indefinite matrices; $Ax=\lambda Mx$}

%\par\noindent{\bf
%THEOREM } 
\theoremNN{\it Each of the following tests is a
necessary and sufficient condition for the real symmetric matrix
$A$ to be positive semidefinite: \newline
(a) $x^T A x\geq 0$ for all vectors $x$.  \newline
(b) All the eigenvalues of $A$ satisfy $\lambda_i\geq 0$.  \newline
(c) All the principal submatrices $A_k$ have nonnegative
determinants.  \newline
(d) There is a matrix $R$, possibly with dependent columns, such that
$A=R^T R$.}
\theoremNNEnd

\subsubsection*{Congruence transformations and law of inertia}

\par\noindent{\bf Definition.} For symmetric matrix $A$ the
transformation
$$A\rightarrow C^T A C ,$$
where $C$ is nonsingular
is called a {\it congruence transformation}.
\medskip

%\par\noindent{\bf
%THEOREM }
\theoremNN{\bf (Sylvester's law of inertia).} {\it The matrix
$C^T AC$ has the same number of positive eigenvalues
as $A$, the same number of negative eigenvalues
and the same number of zero eigenvalues.}
\theoremNNEnd

Maple prefers to investigate polynomials and its
{\tt sturmseq} function is relevant here.
This, and {\tt charpoly} would facilitate illustrative
investigation in connection with this theorem.

\medskip

\subsubsection*{The generalized eigenvalue problem}

This is the problem
$$ A x =\lambda M x ,$$
where $M$ is positive definite.
The mass-spring example before,
with unequal masses, is an example.


%\vfill
%\eject  % Not for 7th ed. Numl methods in same Chapter
\endinput

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

\begin{maple}
#
# Show that matrix Q below is orthog
# and ask Maple to multiply transpose(Q) by Q.
with(linalg,orthog,transpose,multiply): s2:=1/sqrt(2): s3:=1/sqrt(3):
Q:=array([[s3,s2,s2*s3],[s3,0,-2*s2*s3],[s3,-s2,s2*s3]]):
orthog(Q);
#                                   true
multiply(transpose(Q),Q);
# gives the identity matrix verifying that transpose is inverse,
# as characterises real orthog matrices.
\end{maple}





\subsection{Difference equations}
Several applications of eigenvalues are given in 
Maple code dealing with them is given in the Worksheets.


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

\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$ satisfy $\vert\lambda_i\vert < 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}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Here is another exercise concerning least-squares fitting of
lines to data.
The source for the data, etc. is
L.C. Leinbach et al. (eds)
{\it Laboratory approach to teaching calculus} (1991).
(It can be argued that really nonlinear least squares,
which could be used to avoid the numerical differentiation in
this exercise would be a much better approach.
However, we begin with a simplistic approach 
involving linear least squares, and only discuss
nonlinear least squares at the end.)
\smallskip

\par\noindent{\bf EXERCISE.}
The first two statements of the Maple code below give the
{\tt year}, and the corresponding entry in the
{\tt popn} variable, an
estimate of the world's population, in millions, at that year.
(The estimates over the past millenium, i.e.
those given, are, I think, reasonably accurate,
say with an absolute error of $\pm 30$ millions prior to 1950.
The relative errors in recent figures will be less than
it is in the earlier ones.)
\smallskip

(a) Population increase in laboratory experiments
(e.g. with insects) in ordinary unrestricted
favourable circumstances for the species is often
governed by the differential equation,
$$ {{d P}\over{d t}} = k P^r ,\eqno{(DE)}$$
where $r=1$ and this has the solution,
as you learnt in first year,
$$P(t) = P(0) \exp(k t)\qquad {\rm or}\qquad
\log(P(t)) = \log(P(0)) + k t .$$
Run the Maple code fragment below to explore whether
the log of the human population data
is a reasonable fit to this latter
straight line formula.

\begin{verbatim}
with(linalg,augment,stack,transpose):
with(stats,regression,statplot):
years:=array([1000, 1650, 1750, 1800, 1850, 1900,
 1920, 1930, 1940, 1950, 1960, 1976, 1987, 1992]);
popn:=array([200, 545, 728, 906, 1171, 1608, 1834, 2070, 2295,
 2517, 3005, 4000, 5000, 5500]);
py:=map(evalf,stack(map(log,popn),years));
pt:=transpose(augment(array([[lp],[t]]),py));
lp0:='lp0': c:='c'; regression(pt, lp = lp0 + c*t);
assign(");
statplot(pt, lp = lp0 + c*t);
\end{verbatim}

What you should see from this {\tt statplot} is that
the straight-line fit is not particularly good, and that the
data seems to be such that a reasonably smooth curve would
fit close to it. The most important observation is that the
growth seems to be faster than exponential.
Attempts to fit things like
{\tt lp = lp0 + c*t +c2*t{\^{}}8} lead to
better fits, but lead to questions of why one should fit them,
and there is no reason why one should.
\smallskip

(b) Though it probably came from looking at the data
rather than any real `modelling',
population models have been suggested where,
at least for a certain range of time, the population
satisfies the (DE) with some $r>1$.
In first year you learnt that such DEs are `separable' and
can be solved exactly: see part (c).
The next thing is to check that the (DE) with an $r>1$
does indeed seem to fit the data reasonably,
 i.e. we will find $k$ and $r$.
We don't have the derivatives but we might estimate them, by
some method or other. A method we {\it don't} use is
$$ DP_m(t_{mid}) ={{P(t_{m+1}) - P(t_m)}\over
  {t_{m+1}-t_m}} \qquad {\rm at} \qquad
t_{mid}=(t_{m+1}+t_m)/2 .$$
It is a bit unfortunate that the data isn't in nicely spaced
intervals.
However, it isn't, so to estimate the populations at the same
midpoint of times as the above requires an interpolation.
Though more sophisticated interpolation and
numerical differentiation
is possible (see, for example, the treatment of cubic splines
%in Fisher, or 
in Anton and Rorres final
applications chapter \S17), here we just use quadratic
polynomial interpolation at the 3 nearest points.
Read the on-line help concerning {\tt interp}.
(We remark that numerical differentiation is always an
inaccurate process, and is best avoided, but let's ignore this.
More data would probably help, though not too much unless
we replace `interpolation' by `approximation'.)

\begin{verbatim}
m:=12: diffest:=array(1..m): pshift:=array(1..m):
for j from 1 to m do
   diffest[j]:=evalf(subs(tx=years[j+1], diff(
      interp([years[j],years[j+1],years[j+2]],
               [popn[j],popn[j+1],popn[j+2]],tx),tx)));
   pshift[j]:=popn[j+1];
   od;
pd:=map(evalf,stack(map(log,pshift),map(log,diffest)));
pl:=transpose(augment(array([[lps],[ld]]),pd);
cd:='cd': c0:='c0': regression(pl,ld=c0+cd*lps);
assign(");
statplot(pl,ld=c0 + cd*lps);
k:=evalf(exp(c0)); r:=cd - 1;
\end{verbatim}

We know no good reason why
the fit should be as good as it is,
and would like to have more information on the
quality of the data used, error estimates for the data
and so on.
If you find any good references for this, please tell him.
\smallskip

(c) The solution of (DE) is
$$P(t)^r = {1\over{r k (T - t)}} ,$$
where $T$ is a constant to be determined.
The value of $T$ has, flippantly, been called, `doomsday'
by some of the people who noticed that the data fitted this model.
Why?
Of course, the model just won't be accurate very close to $T$.
Use the solution to estimate $T$, for which
you might proceed as follows.

\begin{verbatim}
Tvec:=array(1..m+2):
for j from 1 to (m+2) do
   Tvec[j]:=evalf(years[j]+1/(r*k*popn[j]^r));
od;
with(stats,average,median): average(Tvec);
median(Tvec);
\end{verbatim}

The first entry, {\tt Tvec[1]} is clearly
worth omitting. Find the average of the remaining entries of
{\tt Tvec}.


(Observe that the idea of a long life in Old Testament
Biblical times was {\it three score years and ten} -
though the workers typically lived less than half as long as
this all the way through to about 1800AD.
For most of this class, the above `doomsday' estimate
is before you reach this age of 70.)

\smallskip
\par\noindent (Remark. One way to avoid the numerical
differentiation of part (b) would be to use {\it nonlinear} least
squares to simultaneously find $T$, $k$ and $r$ to fit the data.
This is harder numerically.
There are routines for it in the E04 chapter of the NAG library.)
\medskip


% Removed from end of derivation of linear least squares formula
\begin{maple}
[

- 2

  /  n         \ /  n                \   /  n           \ /  n
\
  |-----       | |-----              |   |-----         | |-----
|
  | \          | | \                 |   | \            | | \
2|
- |  )   2 x[j]| |  )   - 2 y[j] x[j]| + |  )   - 2 y[j]| |  )   2 x[j]
|
  | /          | | /                 |   | /            | | /
|
  |-----       | |-----              |   |-----         | |-----
|
  \j = 1       / \j = 1              /   \j = 1         / \j = 1
/
-------------------------------------------------------------------------
                   /  n         \2     /  n          \
                   |-----       |      |-----        |
                   | \          |      | \          2|
                 - |  )   2 x[j]|  + 2 |  )   2 x[j] | n
                   | /          |      | /           |
                   |-----       |      |-----        |
                   \j = 1       /      \j = 1        /

,

  /  n         \ /  n           \     /  n                \
  |-----       | |-----         |     |-----              |
  | \          | | \            |     | \                 |
  |  )   2 x[j]| |  )   - 2 y[j]| - 2 |  )   - 2 y[j] x[j]| n
  | /          | | /            |     | /                 |
  |-----       | |-----         |     |-----              |
  \j = 1       / \j = 1         /     \j = 1              /
2 ----------------------------------------------------------- ]
              /  n         \2     /  n          \
              |-----       |      |-----        |
              | \          |      | \          2|
            - |  )   2 x[j]|  + 2 |  )   2 x[j] | n
              | /          |      | /           |
              |-----       |      |-----        |
              \j = 1       /      \j = 1        /
\end{maple}
