\magnification=\magstep1
\centerline{Mathematics 172 -- Numerical Analysis 2}
\centerline{Lab/Problem Set 5}
\centerline{April 8, 1999}
\bigskip
\noindent
{\bf Goals}
\bigskip
Let $A = -L + D - U$  be the decomposition of the matrix $A$ into
its diagonal, and negative upper- and lower-triangular parts.
Today's lab will illustrate the comparisons between the Jacobi:
$$x^{(n)} = D^{-1}(L+U)x^{(n-1)} + D^{-1}b$$
and Gauss-Seidel:
$$x^{(n)} = (D - L)^{-1}Ux^{(n-1)} + (D - L)^{-1}b$$
 methods we talked about in class on Tuesday.
We will also investigate the possibility of attaining even 
faster convergence for iterative methods for $Ax = b$, by 
considering the ``parametrized family'' of methods containing Gauss-Seidel:
$$x^{(n)} = (D - \omega L)^{-1}((1-\omega)D + \omega U)x^{(n-1)} + 
(D - \omega L)^{-1}\omega b\leqno(1)$$
(Gauss-Seidel is the case $\omega = 1$).  These methods go by the name
of SOR (``successive over-relaxation'') methods.  
\bigskip
\noindent
{\it Further Features of the {\tt linalg} package}
\bigskip
Maple's {\tt linalg} package contains a number of functions
for matrix arithmetic operations, including the following
operations that will be useful today:
\bigskip
\item{1)}  {\tt multiply} -- matrix multiplication
\item{2)}  {\tt scalarmul} -- multiply a matrix times a scalar
\item{3)}  {\tt inverse} -- compute the inverse of a square matrix
(This is done via the usual method: to find $A^{-1}$, we
reduce the $n \times 2n$ matrix $(A|I)$ to echelon form using 
Gauss-Jordan elimination.  In the resulting matrix $(I|B)$, the 
second block is $B = A^{-1}$.)
\item{4)}  {\tt norm} -- compute the matrix norm (Maple knows
the 1-, 2- , and $\infty$ norms, as well as 
the Frobenius norm as in the exercise from Problem Set 4. 
\item{5)} {\tt eigenvals} -- compute the eigenvalues of a matrix
(real and complex).  We will learn more about the numerical 
methods Maple uses for problems like this in the last unit of the course.
\bigskip
\noindent
Look at the on-line help page for the {\tt linalg} package
for information about the syntax of these commands.  
If you type in the entries of a matrix as real numbers in 
Maple, all arithmetic will be done in floating point.  {\it Be
sure you do that in all calculations today.}  
\bigskip
\noindent
{\it Important Note:}  Maple has special evaluation rules for
matrices.  This means that if you assign a array result to a variable {\tt A},
then tell Maple to print out the value of that variable, you will just
see the variable name as output, not the matrix you want.  To see the 
matrix, use {\tt eval(A);}

Here are sample commands for 
a Maple session illustrating their use.
Suppose we want to compare the rates of convergence for Jacobi
and Gauss-Seidel iteration for the system with $A$ given by
$$A = \pmatrix{6&-2&-1\cr
               -2&7&-2\cr
               -1&-2&8\cr}$$
We begin by defining the $L,D,U$ matrices.  First,
\medskip
\centerline{\tt DD:=array([[6.0,0,0],[0,7.0,0],[0,0,8.0]]);}
\medskip
\noindent
(don't use {\tt D} -- that's reserved for the differentiation operator
for functions in Maple!), and similarly for $L$ and $U$.  Then, for 
example, using the commands above, we can compute the Jacobi $G$-matrix
as follows:
\medskip
\centerline{\tt GJac := multiply(inverse(DD),L+U);}
\medskip
\noindent
The $1$-norm of this matrix can be computed as follows:
\medskip             
\centerline{\tt norm(GJac,1);}
\medskip
\noindent
You can also get the more refined information of the spectral radius
$\rho(G_{Jac})$ from the output of:
\medskip
\centerline{\tt eigenvals(GJac);}
\medskip
\noindent
(How?)  The Gauss-Seidel $G$-matrix is
\medskip
\centerline{\tt GGS := multiply(inverse(DD-L),U);}
\medskip
\noindent
You can compute its matrix norms and spectral radius for
$G_{GS}$ in the same way you did for $G_{Jac}$ above.
How do they compare to the norms and spectral radius for $G_{Jac}$ for 
this matrix?
\bigskip
\noindent
{\bf Lab Questions}
\bigskip
\item{A)}  Consider the matrix:
$$A = \pmatrix{10&-2&-1&0\cr
               -2&10&-5&0\cr
               -1 & -5 & 8 & -2\cr
               0 & 0 & -2 & 5\cr}$$
\itemitem{1)}
Note that $A$ is a symmetric matrix, so by the spectral
theorem from linear algebra, all of its eigenvalues should be real.  Check
that in fact all of them are strictly positive.  We say a symmetric matrix 
$A$ is {\it positive-definite} if all of its eigenvalues are $> 0$.
\itemitem{2)}  Compute the Jacobi and Gauss-Seidel $G$-matrices for
this $A$ and compare their $1,2,\infty$-norms and spectral radii.  
Which iterative method will converge faster for a system $Ax = b$ for 
this coefficient matrix $A$?  (Note that the theorem we quoted in class
comparing the Jacobi and Gauss-Seidel methods {\it does\/} apply here -- 
this result is also Theorem 7.22 in the text.)
\itemitem{3)}  Find a $4\times 4$
matrix $A'$ with $a'_{ii} > 0$ all $i$, $a'_{ij} \le 0$
if $i \ne j$ for which $1 < \rho(G'_{Jac}) < \rho(G'_{GS})$.  (This is the 
second case in Theorem 7.22 from the text.)

\itemitem{4)}  Now consider the one-parameter family of matrices for 
SOR iteration:
$$G_\omega = (D - \omega L)^{-1}((1-\omega)D + \omega U)$$
Using Maple, determine the interval of $\omega$-values 
for which $\|G_\omega\|_\infty < 1$, 
and estimate the value of $\omega$ that {\it minimizes}
$\|G_\omega\|_\infty$.  (Suggestion:  You can {\it plot\/}
$\|G_\omega\|_\infty$ as a function of $\omega$ in Maple(!).)
\itemitem{5)}  Now, try to estimate the value of $\omega$ that
minimizes the {\it spectral radius\/} of $G_\omega$.  (It is harder, though
certainly not impossible, to generate a plot in this case; if you don't 
see how to do that, either ``trial and error'' or some kind of systematic
search will be necessary.)  Is this the same $\omega$ as in part 3?  Explain.
\itemitem{6)}  Perform 4 iterations of Gauss-Seidel and SOR (see equation (1))
above) using your optimal $\omega$ from part 4 to solve the system
$Ax = b$, where $b = (1,0,1,0)^t$.  Discuss your results.
\bigskip
\noindent
{\it Comment:}  As you might guess, the choice of the optimal $\omega$
for SOR iteration is something of a ``black art''.  There is one case,
though, in which the optimal $\omega$ can be determined exactly.  
If $A$ is {\it positive definite} and {\it tridiagonal}, then it can
be proved that the $\omega$ minimizing the spectral radius of 
$G_\omega$ is 
$$\omega_{opt} = {2 \over 1 + \sqrt{1 - (\rho(G_{Jac}))^2}}$$
where $\rho(G_{Jac})$ is the spectral radius of the corresponding 
Jacobi $G$-matrix.  This information may be useful for the next problem!
\bigskip
\item{B)}  Problem 12 from Section 7.3 in the text.  
Use Gauss-Seidel and/or SOR iteration to find approximate
solutions accurate to a tolerance of $10^{-5}$ in the $\infty$-norm.
Of course,  one of the main questions you will need to address is: When 
do you want to stop iterating?  As part of your solution explain the stopping
criteria you used.  Implementation note:  
You {\it DON'T} want to enter a $50 \times 50$ or $100 \times 100$ matrix entry 
by entry(!)  Since most of the entries are zero anyway, it will be 
simpler to set up a matrix of all zeroes, then fill in the non-zero
ones to get $A$ (or $L,D,U$).  The initialization:
\medskip
\centerline{\tt B := array(1..50,1..50,sparse);}
\medskip\item{}
sets up a $50\times 50$ matrix of all zeroes.  Then you can change 
an entry to something else when necessary.  {\tt B[i,j]} is the Maple notation
for the entry in row $i$ and column $j$ in $B$.
For instance, the command {\tt B[3,4] := .4;}
would set the entry in row 3 and column 4 in $B$ to $.4$.  Even now, 
{\it think} about how you can enter all the entries for $L$, $D$, $U$
before setting them manually one by one!!
\bigskip
\noindent
{\bf Assignment}
\bigskip
Individual write-ups due in class on April 15.  
\bye
