% Hill6.tex

\section{Further Directions}

\subsection{Function Spaces}
This section applies the ideas of Hill's Chapter 4 to the vector space 
$C[a,b]$ described in Examples~9 and 10 of Hill Section 3.3. In section 4.1,
the inner product
 $$ f\cdot g = \int_a^b f(x) g(x)\, dx $$
was introduced, and will be used throughout this section, together with its
Maple code:
\begin{verbatim}
> dot:= proc(fx,gx,lims)        int( expand(fx*gx),x=lims )        end:
\end{verbatim}

\subsubsection*{Fourier Series and Legendre Polynomials}
Hill discusses the Fourier series of functions in $C^1[0,2\pi]$.\quad
These are based on the following
\setcounter{theorem}{4}
\begin{theorem} The functions $\cos(nx)$ and $\sin(nx)$ are orthogonal in the
  sense that, for all integers $m,n\ge 0$,
  \begin{eqnarray*}
    & \cos (nx) \cdot \sin (mx) = 0,& \\
    \cos(nx)\cdot \cos(mx) &= 0,\qquad \sin(nx)\cdot \sin(mx) &= 0\qquad
                                                             {\rm for}\ n\ne m
  \end{eqnarray*}
\end{theorem}

\setcounter{theorem}{8}
\begin{theorem} Let $f\in C^1[0,2\pi]$. Then, for all $x\in [0,2\pi]$,
    $$ f(x) = a_0 +a_1\cos(x) +b_1\sin(x) +a_2\cos(2x) +b_2\sin(2x) +\cdots $$
  where
    $$ a_0 = {f\cdot 1\over 1\cdot 1}\qquad
       a_n = {f\cdot \cos(nx)\over \cos(nx)\cdot \cos(nx)}\qquad
       b_n = {f\cdot \sin(nx)\over \sin(nx)\cdot \sin(nx)}\quad {\rm for}~n>0.$$
\end{theorem}

\headq{EXERCISES} Use the Maple procedure\ \ \verb$dot$\ \ defined above
  (with\ \ \verb$lims:= 0..2*Pi$) to
  \begin{enumerate}
  \item prove the identities of theorem 6.5;
  \item find the coefficients $a_n, b_n$ in theorem 6.9 for the function
    $f(x)=x$ (see Section 3.3 for the Maple command for this function);
  \item (Harder) Write a Maple procedure to find the Fourier series for a
    general function.
  \end{enumerate}

Fourier series are very useful in many branches of Engineering. However, you
may be more accustomed to approximating functions by Power Series.
There are several ways of defining an inner product of polynomials
and so obtaining an orthonormal set of polynomials.  Hill gives the example
of Legendre polynomials which are orthogonal on the interval $[-1,1]$.
Maple includes these in the package {\tt orthopoly}\index{orthopoly}, and
the following Maple code gives the first three Legendre polynomials.
\begin{verbatim}
# Hill Section 6.1 Legendre Polynomials
> with(orthopoly):
# The first three Legendre polynomials (to within a constant factor) are
> f0:= P(0,x):            f1:= P(1,x):            f2:= P(2,x): 
> print( f0,f1,f2 );
                                             2
                            1,    x,    3/2 x  - 1/2

# Using the dot product given in Section 4.1, we can check their orthogonality:
> print( dot(f0,f1,-1..1), dot(f0,f2,-1..1), dot(f1,f2,-1..1) );
                                    0, 0, 0
\end{verbatim}

\subsubsection*{Least-Squares Approximation and Line of Best Fit}
In Section 4.4, Hill has shown that the line of best fit for the points
$(x_1,y_1),\ldots (x_n,y_n)$ can be found from the least-squares solution of
the equation\quad $A\b x=\b y$\quad where
  $$ A = \pmatrix{1&x_1\cr \vdots&\vdots\cr 1&x_n\cr}\quad
        \b x = \pmatrix{b\cr m}\quad
                  {\rm and}\quad \b y = \pmatrix{ y_1\cr \vdots\cr y_n\cr}. $$
In Section 6.1, he discusses two approaches to finding the line\quad $y=mx+b$
\quad of best fit for for a function $f(x)$ defined on the interval $[0,1]$.
The first of these consists of adapting the above method by replacing
\quad $A$ by $(1,x)$\quad and\quad $\b y$ by $f(x)$.\quad The least-squares
solution becomes
\begin{eqnarray*}
    A^T A \b x &=& A^T \b y \\
  \pmatrix{1\cr x} (1,x) \pmatrix{b\cr m} &=& \pmatrix{1\cr x} f(x)
\end{eqnarray*}
where ``multiplication'' by the functions 1 and $x$ is replaced by the
inner product defined above.

The second approach is to use the first two terms of the power series for
$f(x)$:
 $$ f(x) = b + mx + \cdots . $$
This can be obtained by the same method as for a Fourier Series:
if $v_0,v_1,\ldots$ are orthogonal polynomials in $C[0,1]$ with $\deg(v_n)=n$,
then
 $$ f(x) = a_0 v_0 + a_1 v_1 + \cdots\quad {\rm where}\quad
     a_n = {f\cdot v_n\over v_n\cdot v_n}. $$

\noindent The following Maple code illustrates both these ideas.

\begin{verbatim}
# Hill Section 6.1 Example on p400-402.
# (See above for the definition of the procedure  dot )
> one_one:= dot(1,1,0..1):     one_x:= dot(1,x,0..1):     x_x:= dot(x,x,0..1):
> with(linalg):
> A:= matrix([ [one_one,one_x], [one_x,x_x] ]);
                                    [  1   1/2 ]
                               A := [          ]
                                    [ 1/2  1/3 ]

# The line of best fit for the function  sqrt(x)  over the range 0..1 is
> b_m:= linsolve( A, vector( [dot(1,sqrt(x),0..1),dot(x,sqrt(x),0..1)] ) ):
> y = b_m[1] + b_m[2] * x;
                                y = 4/15 + 4/5 x

# Alternatively, we find the orthogonal polynomials (using Gram Schmidt):
> v0:= 1:                                            v1:= x - one_x * v0;
                                 v1 := x - 1/2

# Since  v0, v1  are orthogonal, the line of best fit for any function f is
> best_fit:= proc(f)
>     dot(f,v0,0..1)/dot(v0,v0,0..1)*v0 + dot(f,v1,0..1)/dot(v1,v1,0..1)*v1
> end:
\end{verbatim}
\begin{verbatim}
# So the line of best fit for the function  sqrt(x)  over the range 0..1 is
> best_fit( sqrt(x) );
                                  4/15 + 4/5 x
\end{verbatim}

\subsection{Singular Value Decomposition}
\headq{Note} You should re-read Sections 4.3 and 5.4 before reading this
section.

\setcounter{theorem}{23}
\begin{theorem}[Singular Value Decomposition] Any $m\times n$ matrix $A$
  can be factored as
    $$ A = V\Sigma U^T $$
  where $V$ is an $m\times m$ orthogonal matrix, $U$ is an $n\times n$
  orthogonal matrix, and $\Sigma$ is an $m\times n$ matrix of the form
    $$ \Sigma = \pmatrix{ D&0\cr 0&0\cr }\qquad {\rm where}\qquad
                 D = \pmatrix{ \sigma_1&&\cr &\ddots&\cr &&\sigma_k\cr }. $$
  {\em The numbers $\rowi{\sigma}{k}$ are all positive real numbers called the
  {\bf singular values} of $A$.}
\end{theorem}
Maple has a function {\tt singularvals}\index{singularvals} in its {\tt linalg}
package to find singular values, and {\tt Svd}\index{Svd} not in its
{\tt linalg} package to find the matrices $V$, $\Sigma$ and $U$.
You should get help on these two functions.

\subsubsection*{The Pseudoinverse and Least-squares Problems}
The {\bf pseudoinverse} of the matrix $A$ is defined to be\qquad\qquad
$ A^+ = (A^T A)^{-1} A^T$.

\setcounter{theorem}{29}
\begin{theorem} If $A$ is an $m\times n$ matrix with singular value
  decomposition\quad $A=V\Sigma U^T$,\quad then
    $$ A^+ = U \Sigma^+ V^T\qquad {\rm where}\qquad
                                \Sigma^+ = \pmatrix{ D^{-1}&0\cr 0&0\cr }^T $$
\end{theorem}
The following is a Maple procedure to find the pseudoinverse, and its
application to an example.
\begin{verbatim}
> with(linalg):
> pseudoinverse:= proc(A)
>     local U,V,Sigma_plus,singvals,i:
>     Sigma_plus:= array( sparse,1..coldim(A),1..rowdim(A) ):
>     singvals:= evalf(Svd(A,V,U));
>     for i from 1 to vectdim(singvals) do
>         if singvals[i]>0 then  Sigma_plus[i,i]:= evalf( 1/singvals[i] )  fi:
>     od:
>     map( evalf, evalm(U&*Sigma_plus&*transpose(V)) );
> end:
#
# Hill Section 4.3, Example 5, p264/5
> A:= matrix([ [1,3], [2,1], [0,1], [-2,1] ]):
# The matrix of the perpendicular projection function is P where
> P:= evalm( A&*pseudoinverse(A) ):
\end{verbatim}
\goodbreak
\begin{verbatim}
> map( evalf,evalm(33*P) );   # The previous answer had 33 in the denominators
            [ 25.00000000   9.999999992  8.000000003   5.999999996 ]
            [                                                      ]
            [ 9.999999996   15.00000001  .9999999999  -13.00000001 ]
            [                                                      ]
            [ 7.999999999   .9999999956  3.000000000   5.000000000 ]
            [                                                      ]
            [ 6.000000003  -13.00000002  5.000000000   23.00000000 ]
\end{verbatim}

\setcounter{theorem}{31}
\begin{theorem} The solution of the least-squares problem for the equation
  $A\b x=\b b$ is given by
    $$ \bar \b x = A^+\b b. $$
\end{theorem}
Thus the above procedure for pseudoinverses may be applied to Example 6 of Hill
Section 4.3, p265.
\begin{verbatim}
# Hill Section 4.3, Example 6, p265/6
> A:= matrix([ [40,1], [45,1], [50,1], [55,1], [60,1] ]):
> y:= vector( [482,467,452,433,421] ):
# The least-squares solution to the equation  A x = y  is given by
> map( evalf,evalm(pseudoinverse(A)&*y) );
                          [ -3.12000000, 607.0000012 ]
\end{verbatim}
Hill points out that this can be used for dealing with ill-conditioned systems
by setting the small singular values equal to zero and treating the system as a
least-squares problem.

\subsection{Iterative Methods}
The methods of finding eigenvalues and eigenvectors given in Hill Section 5.8
are all examples of {\em iterative methods}. If the number $x$ is the
solution of an equation which may be written in the form\quad $f(x)=x$,\quad
then the value of $x$ may be estimated iteratively as follows. If $x_n$ is an
estimate of the value of $x$, then a (hopefully better) estimate is
  $$ x_{n+1} = f(x_n). $$
The following piece of Maple code reproduces examples 1 to 4 of Hill
Section 6.3.  (Some blank lines have been omitted to save space.)

\begin{verbatim}
# Hill Section 6.3   Examples of Iterative Methods pp 414-418
# Example 1: Solve  cos(x) = x
> x:= 0.7:             # initial value
> for i from 1 to 4 do        x:= cos(x):        print(x);        od;
                                  .7648421873
                                  .7214916396
                                  .7508213288
                                  .7311287726
\end{verbatim}
\begin{verbatim}
# Example 2: Find the square root of 2.
> x:= 1:
> for i from 1 to 3 do     x:= evalf( (x+2/x)/2 ):      print(x);      od;
                                  1.500000000
                                  1.416666667
                                  1.414215686
                                  1.414213562
\end{verbatim}
\begin{verbatim}
> with(linalg):
# Example 3: Jacobi's method for solving the linear equations  A x = b  where
> A:= matrix([ [2,-1],
>              [-1,2] ]):
> b:= vector([3,6]):                                         m:= 6:
> x:= vector([0,0]):
> for i from 1 to m do
>     newx1:= ( b[1] - A[1,2]*x[2] )/A[1,1]:
>     newx2:= ( b[2] - A[2,1]*x[1] )/A[2,2]:
>     x[1]:= evalf(newx1):            x[2]:= evalf(newx2):            print(x);
> od:
                              [ 1.500000000, 3. ]
                          [ 3.000000000, 3.750000000 ]
                          [ 3.375000000, 4.500000000 ]
                          [ 3.750000000, 4.687500000 ]
                          [ 3.843750000, 4.875000000 ]
                          [ 3.937500000, 4.921875000 ]
\end{verbatim}
\begin{verbatim}
# Example 4: Gauss-Seidel method    (Note also how much neater the code is)
> x:= vector([0,0]):
> for i from 1 to m do
>     x[1]:= evalf( (b[1]-A[1,2]*x[2])/A[1,1] ):
>     x[2]:= evalf( (b[2]-A[2,1]*x[1])/A[2,2] ):         print(x);
> od:
                          [ 1.500000000, 3.750000000 ]
                          [ 3.375000000, 4.687500000 ]
                          [ 3.843750000, 4.921875000 ]
                          [ 3.960937500, 4.980468750 ]
                          [ 3.990234375, 4.995117188 ]
                          [ 3.997558594, 4.998779297 ]
\end{verbatim}
\begin{verbatim}
# The matrix form of the Gauss-Seidel method (written as a procedure) is
> GSsolve:= proc(A,b,m)     # m is the number of iterations
>     local n,DD,L,U,invLplusD,x,i,j:
>     n:= rowdim(A):
>     L:= array( sparse,1..n,1..n ):      U:= array( sparse,1..n,1..n ):
>     DD:= array( sparse,1..n,1..n ):
>     for i from 1 to n do
>         DD[i,i]:= A[i,i]:
>         for j from 1 to i-1 do    L[i,j]:= A[i,j]:    U[j,i]:= A[j,i]:    od:
>     od:
>     invLplusD:= inverse(L+DD):
>     x:= vector( [0$n] ):     #  0$n  denotes a sequence of n 0's
>     for i from 1 to m do
>         print(x);
>         x:= map( evalf,evalm( invLplusD &* (b-U&*x) ) ):
>     od;
> end:
#
\end{verbatim}
\goodbreak
\begin{verbatim}
# So Example 4 (above) becomes
> GSsolve(A,b,6);
                                    [ 0, 0 ]
                          [ 1.500000000, 3.750000000 ]
                          [ 3.375000000, 4.687500000 ]
                          [ 3.843750000, 4.921875000 ]
                          [ 3.960937500, 4.980468750 ]
                          [ 3.990234375, 4.995117188 ]
                          [ 3.997558594, 4.998779297 ]
\end{verbatim}

\subsection{Matrix Norms}\index{norm}
\headq{Definition 6.68} If $A$ is an $n\times n$ matrix, then the {\bf norm}
of $A$ is
  $$ \Vert A\Vert = \max \Vert A\b x\Vert / \Vert\b x\Vert $$
where $\b x\in \DoubleR^n$ and $\b x\ne \b 0$.

\smallskip
Maple has a function {\tt norm} which accepts as its first argument either a
vector or a square matrix. The second argument is a positive integer or
{\tt infinity} which determines which norm is to be used.

The norm used by Hill is the 2-norm (which will be denoted by $\Vert A\Vert_2$
here), although he also introduces the $\infty$-norm which is given by

\smallskip \headq{Definition 6.73}
$ \Vert \b x\Vert_\infty = max\{ \vert x_1\vert,\ldots \vert x_n\vert \}\qquad 
   \Vert A\Vert_\infty = max\{ \vert a_{i1}\vert + \cdots \vert a_{in}\vert\}$

\setcounter{theorem}{71}
\begin{theorem} If $A$ is an $n\times n$ matrix, then $\Vert A\Vert_2$ is the
  largest singular value of $A$.
\end{theorem}
\begin{verbatim}
# Hill Section 6.4 Example 5 p424
> with(linalg): 
> A:= matrix([ [1,1], [0,2] ]):
> singularvals(A);
                         1/2        1/2        1/2        1/2
                  [1/2 10    + 1/2 2   , 1/2 10    - 1/2 2   ]

# Note that, before Release 3, Maple did not simplify the answer to this
> norm(A,infinity);
                                       2

> norm(A,2);
                                    1/2        1/2
                              1/2 10    + 1/2 2

# Before Release 2, Maple could not decide that this was the larger one
> evalf(");
                                  2.288245611
\end{verbatim}

Hill points out that any sequence $\b x_k$ defined by
  $$ \b x_{k+1} = B \b x_k + \b c, $$
where $B$ is a matrix and $\b c$ is a vector, always converges when
$\Vert B\Vert<1$ (either norm).

\medskip
\headq{Definition 6.86} If $A$ is an invertible matrix, then the
{\bf condition number} of $A$ is
  $$ c(A) = \Vert A\Vert \times \Vert A^{-1}\Vert $$
where the definition depends on the choice of the norm.

\smallskip
Hill points out that the condition number of an invertible matrix is
approximately the same irrespective of the choice of the norm.  Also,
in the case of the 2-norm, $c(A)=|\lambda_1/\lambda_n|$ where $\lambda_1$ is
the largest and $\lambda_n$ the smallest eigenvalue of $A$ (in absolute terms).

\subsection{General Vector Spaces and Linear Transformations}
Hill defines a field, and gives the examples of $\DoubleQ$, $\DoubleR$,
$\DoubleC$ and $\DoubleZ_p$, where $p$ is a prime number. Maple's use of the
first three has already been considered. To calculate in $\DoubleZ_p$, or in
a vector space over $\DoubleZ_p$, you must use the {\tt modp}
(or {\tt mod}) function (not in {\tt linalg}). An example follows.

\begin{verbatim}
# Hill Section 6.5 Examples 7-10, pp 435-437
> with(linalg):
> alpha1:=vector([1,0,0]):  alpha2:=vector([1,1,0]):  alpha3:=vector([1,1,1]):
> ALPHA:= augment( alpha1,alpha2,alpha3 );
                                       [ 1  1  1 ]
                                       [         ]
                              ALPHA := [ 0  1  1 ]
                                       [         ]
                                       [ 0  0  1 ]

# Since the field is  Z2,  we must use the function  modp( - , 2)
# Example 7
# {alpha1,alpha2,alpha3}  is a basis if and only if  ALPHA  is invertible:
> modp( det(ALPHA),2 );
                                       1

# Example 8
> v:= map( modp, evalm(alpha1+alpha3), 2 );
                                v := [ 0, 1, 1 ]

> modp( dotprod(v,v),2 );
                                       0

# So  norm(v) = 0  although  v <> 0
#
# The coordinate vector of v relative to the basis ALPHA is  vAlpha  where
> vAlpha:= map( modp, linsolve(ALPHA,v), 2);
                             vAlpha := [ 1, 0, 1 ]

# Example 9
# Find the matrix of the linear transformation  T  where
> T:= proc(x)    vector( [x[1]+x[3], x[2]+x[3] ])        end:
#    relative to the bases  ALPHA and BETA  where
> beta1:= vector( [0,1] ):       beta2:= vector( [1,1] ):
> BETA:= augment( beta1,beta2 );
                                        [ 0  1 ]
                                BETA := [      ]
                                        [ 1  1 ]

# Use the coordinates of  T(alpha1), T(alpha2), T(alpha3)  as in Section 3.8:
> T_BetaAlpha:= augment( map( modp, linsolve(BETA,T(alpha1)), 2),
>                        map( modp, linsolve(BETA,T(alpha2)), 2),
>                        map( modp, linsolve(BETA,T(alpha3)), 2) );
\end{verbatim}
\goodbreak
\begin{verbatim}
                                          [ 1  0  0 ]
                           T_BetaAlpha := [         ]
                                          [ 1  1  0 ]

# Example 10
# The coordinates of T(v) relative to the basis  BETA  is  Tv_Beta  where
> Tv_Beta:= map( modp, linsolve(BETA,T(v)), 2);
                              Tv_Beta := [ 1, 1 ]

# By Theorem 6.97, this should be the same as the following
> map( modp, evalm( T_BetaAlpha &* vAlpha ), 2);
                                    [ 1, 1 ]
\end{verbatim}

