\documentclass[twoside]{article}

\usepackage{palatino}
\usepackage{camnum}

\textheight=21.5cm
\textwidth=15cm
\topmargin=0cm
\oddsidemargin=0.5cm
\evensidemargin=0.5cm
\parindent=0pt
\parskip=6pt

\setcounter{section}{2}
\setcounter{subsection}{2}
\setcounter{equation}{0}

\begin{document}
\begin{center}
  {\large \textbf{Mathematical Tripos Part IB: Easter 2001}}

  {\Large \textbf{Numerical Analysis -- Lecture
      3\footnote{Corrections and suggestions to these notes should
        be emailed to \texttt{a.shadrin@damtp.cam.ac.uk}. All handouts
        are available on the WWW at the URL 
        \texttt{http://www.damtp.cam.ac.uk/user/na/PartIB/Handouts.html}.}}}
\end{center}

{\bf Band matrices.} The matrix $A$ is a {\em band matrix\/} if there
exists an integer $r<n$ such that $A_{i,j}=0$ for $|i-j|>r$,
$i,j=1,2,\ldots,n$. In other words, all the nonzero elements of $A$
reside in a band of width $2r+1$ along the main diagonal. In that case,
according to the statement from the end of the last lecture, $A=LU$
implies that $L_{i,j}=U_{i,j}=0$ $\forall$ $|i-j|>r$ and sparsity
structure is inherited by the factorization.\\
In general, the expense of calculating an LU factorization of an
$n\times n$ {\em dense\/} matrix $A$ is $\O{n^3}$ operations and the
expense of solving $A\MM{x}=\MM{b}$, provided that the
factorization is known, is $\O{n^2}$. However, in the case of a banded
$A$, we need just $\O{r^2n}$ operations to factorize and $\O{rn}$
operations to solve a linear system. If $r\ll n$ this represents a
very substantial saving!

\section{Iterative methods for linear systems}

\subsection{Basic iterative schemes}

Solution of $A\MM{x}=\MM{b}$ by factorization is frequently very
expensive for large $n$, even if we exploit sparsity. An alternative
is to use {\em iterative methods\/} which are very efficient.
An example of an iterative scheme
is to write $A=B-C$, where \textbf{\textit{(1)\/}} $B$ \& $C$ are
$n\times n$ matrices; \textbf{\textit{(2)\/}} $B$ is
nonsingular; \textbf{\textit{(3)\/}} the system $B\MM{x}=\MM{c}$ is
\emph{easy to solve\/} and \textbf{\textit{(4)\/}}
the matrix $C$ is somehow `small' in comparison with $B$.  We write
the original system in the form $B\MM{x}=C\MM{x}+\MM{b}$ and consider 
solving it by iteration. Choose an arbitrary $\MM{x}_0\in\BB{R}^n$ and
define $\MM{x}_{m+1}$, $m=0,1,\ldots$, by solving
\begin{equation}\label{3.1}
  B\MM{x}_{m+1}=C\MM{x}_m+\MM{b}.
\end{equation}
Provided that $B$ is, for example, banded, the solution of \R{3.1} is cheap
(and the LU factorization of $B$ can be re-used -- an example of why
the LU formalism is superior to Gaussian elimination). Often the
sequence $\{\MM{x}_m\}_{m=0}^\infty$ converges to the solution of
$A\MM{x}=\MM{b}$.

{\bf The Jacobi iteration.} We write
$A=A_{\CC{D}}+A_{\CC{L}}+A_{\CC{U}}$, where $A_{\CC{L}}$, 
$A_{\CC{D}}$ and $A_{\CC{U}}$ are portions of $A$: subdiagonal
(strictly lower triangular), diagonal and supdiagonal
(strictly upper triangular), respectively. 
Suppose that no diagonal element of $A$ is
zero. The {\em Jacobi iteration\/} is
\begin{equation}\label{3.2}
  A_{\CC{D}}\MM{x}_{m+1}= -(A_{\CC{L}}+A_{\CC{U}})\MM{x}_m+\MM{b},
  \qquad m=0,1,\ldots.
\end{equation}

{\bf The Gauss--Seidel iteration.} In the above notation, it takes the
form
\begin{equation}\label{3.3}
  (A_{\CC{D}}+A_{\CC{L}})\MM{x}_{m+1}= -A_{\CC{U}}\MM{x}_m+\MM{b},
  \qquad m=0,1,\ldots.
\end{equation}
Note that $A_{\CC{L}}+A_{\CC{D}}$ is lower triangular, hence the
solution of \R{3.3} is cheap.

\subsection{Necessary and sufficient conditions for convergence}

Suppose that $A$ is nonsingular and denote by $\MM{x}^*$ the solution
of $A\MM{x}=\MM{b}$. Having written $A=B-C$, we examine the
iterative scheme \R{3.1}. (Note that both \R{3.2} and \R{3.3} can be cast in
this form.) Our goal is to identify conditions so that
$\MM{x}_m\rightarrow\MM{x}^*$, regardless of the choice of
$\MM{x}_0\in\BB{R}^n$.

Subtract $B\MM{x}^*=C\MM{x}^*+\MM{b}$ from \R{3.1}. This gives
$B(\MM{x}_{m+1}-\MM{x}^*)=C(\MM{x}_m-\MM{x}^*)$, hence
$B\MM{\varepsilon}_{m+1}=C\MM{\varepsilon}_m$, where
$\MM{\varepsilon}_m:=\MM{x}_m-\MM{x}^*$ is 
the error in the $m$th iterate. Since $B$ is nonsingular (otherwise we
cannot execute \R{3.1} in the first place), it follows that
\begin{equation}\label{3.4}
  \MM{\varepsilon}_{m+1}=H\MM{\varepsilon}_m=\cdots=H^{m+1}
  \MM{\varepsilon}_0,\quad m=0,1,\ldots  
  \quad \mbox{where} \quad H=B^{-1}C\quad \mbox{is the
  \textit{iteration matrix.}}.
\end{equation}
This indicates that the errors tend to zero as $m\rightarrow\infty$
(regardless of the choice of $\MM{x}_0$) provided that
$\lim_{m\rightarrow\infty}H^m=O$. The next theorem establish a criterium 
for such a convergence.

We employ the notation $\rho(P)$ for the magnitude of the largest (in
absolute value) eigenvalue of the $n\times n$ matrix $P$. The quantity
$\rho(P)$ is called the \emph{spectral radius\/} of the matrix
$P$. ({\em Note:\/} Recall that, even if $P$ is real, its eigenvalues
might be complex.)

{\bf Theorem.} We have
$\lim\limits_{m\rightarrow\infty}H^m\MM{\varepsilon}_0=0$ 
for all
$\MM{\varepsilon}_0\in\BB{R}^n$ if and only if $\rho(H)<1$.\\
\hspace*{8pt}{\bf Proof.} We commence with the case $\rho(H)\geq1$ and
wish to demonstrate that, for some {\em real} $\MM{\varepsilon}_0$,  
the vector $\MM{\varepsilon}_m$ need not tend to $\MM{0}$. Let
$\lambda$ be an eigenvalue of (the real) $H$, real or complex,  
such that $|\lambda|=\rho(H)$, and
let $\MM{w}$ be a corresponding eigenvector, i.e., $H\MM{w}=\lambda\MM{w}$.
Denote the {\em Euclidean length\/} of
$\MM{p}\in\BB{C}^n$ by  
$\|\MM{p}\|=\left(\sum_{k=1}^n |p_k|^2\right)^{1/2}.$
Then the vector $H^m \MM{w} = \lambda^m\MM{w}$
is bounded away from zero, because 
%
\begin{equation}\label{3.5}
    \|H^m\MM{w}\| 
 =  |\lambda|^m\|\MM{w}\|
\ge \|\MM{w}\| \ge \max_k |\MM{w}_k| = \gamma > 0.
\end{equation}
%
If $\MM{w}$ is real, we choose $\MM{\varepsilon}_0=\MM{w}$, hence
$\MM{\varepsilon}_m=H^m\MM{\varepsilon}_0 \ge \gamma$, 
and this cannot tend to zero.


If $\MM{w}$ is complex, then $\MM{w} = \MM{u} + {\rm i}\MM{v}$
with  $\MM{u},\MM{v}$ some real vectors. But then at least one
of the sequences $(H^m\MM{u}),(H^m\MM{v})$ does not tend to zero.
For if both do, then also
$$
   H^m \MM{w} = H^m \MM{u} + {\rm i} H^m \MM{v} \to \MM{0},
$$
and this contradicts \R{3.5}.
This completes the proof of the `only if' part of the theorem.

Now we turn to the `if' case, i.e., $\rho(H)<1$, where we assume 
for simplicity that $H$ possesses $n$ linearly independent eigenvectors
$\MM{w}_1,\MM{w}_2,\ldots,\MM{w}_n$, say. Hence $H\MM{w}_j=\lambda_j
\MM{w}_j$, $|\lambda_j|<1$, $j=1,2,\ldots,n$. Linear 
independence means that every $\MM{\varepsilon}\in\BB{R}^n$ 
can be expressed as a linear combination of the eigenvectors. 
Therefore, given $\MM{\varepsilon}_0\in\BB{R}^n$, there exist
$\alpha_1,\alpha_2,\ldots,\alpha_n\in\BB{C}$ 
such that $\MM{\varepsilon}_0=\sum_{j=1}^n \alpha_j \MM{w}_j$. Thus, 
\begin{displaymath}
  \MM{\varepsilon}_1=H\MM{\varepsilon}_0=\sum_{j=1}^n
  \alpha_j\lambda_j\MM{w}_j \qquad\mbox{and, by induction,}\qquad
  \MM{\varepsilon}_m=\sum_{j=1}^n \alpha_j\lambda_j^m \MM{w}_j
\end{displaymath}
for all $m=0,1,\ldots$. Since $\rho(H)<1$, it follows that
$\lim_{m\rightarrow\infty} \MM{\varepsilon}_m=\MM{0}$, as required. \QED

{\bf Remark.} 
The complete proof of the `if' case exploits the so-called Jordan
normal form of the matrix $H$, 
$
  H = SJS^{-1},
$ 
where $J$ is a block diagonal matrix consisting of the Jordan blocks
$$
  J_i =
\left[\begin{array}{cccc}
           \lambda_i &      1    &           &   \\
                     &  \ddots   &   \ddots  &   \\
                     &           & \lambda_i & 1 \\
                     &           &           & \lambda_i
\end{array}\right], \qquad
  J_i \in \BB{R}^{n_i\times n_i}.
$$
To prove that $J_i^m \to 0$ if  $|\lambda_i| < 1$, 
one can split $J_i = \lambda_i I + P$, notice that $P^k = 0$
for $k \ge n_i$, and evaluate the terms of expansion $(\lambda_i I + P)^m$. 

{\bf Application 1.} If $A=A^{\rm T}$ is positive definite, then the 
Gauss--Seidel method converges [to a solution of $A\MM{x}=\MM{b}$
for any starting vector $\MM{x}_0$] (for then $\rho(H_{\rm GS})<1$). 

{\bf Application 2.} If both $A=A^{\rm T}$ and $2A_{\rm D} - A$ are positive
definite, then the Jacobi method converges (for then $\rho(H_{\rm J})<1$ 
as well). 





\end{document}

