\documentclass[a4paper,11pt]{article}
\usepackage{a4wide}
\usepackage{subfigure}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{epsfig}
\usepackage{theorem}
\usepackage{verbatim}
\usepackage[english]{babel}

\def\bx#1{{\bf #1}}
\def\rref#1{(\ref{#1})}

\newtheorem{theorem}{Theorem}

{\theoremheaderfont{\bfseries}\theorembodyfont{\mdseries}
\newtheorem{exercise}{\underline{Exercise}}}
\parindent=0pt

\begin{document}

\begin{center}
\begin{tabular}{c}
Uppsala University\\
Graduate School in Mathematics and Computing\\
Institute for Information Technology\\
\rule{6in}{1pt}
\end{tabular}
\end{center}

\begin{center}
Numerical Linear Algebra FMB and MN1\\
Fall 2007 \\[1cm]
{\bf Mandatory Assignment 3a:}\\[5pt]
{\Large \bf \em Iterative solution methods and their rate of convergence}\\
\end{center}
 
 
\begin{exercise}[Implementation tasks]

Implement in {\tt Matlab} the following iterative schemes:

\begin{itemize}
\item[M1:] the (pointwise) Jacobi method;
\item[M2:] the second order Chebyshev iteration, described in Appendix A;
\item[M3:] the unpreconditioned conjugate gradient (CG) method (see the lecture notes);
\item[M4:] the Jacobi (diagonally) preconditioned CG method (see the lecture notes).
\end{itemize}


For methods M3 and M4, you have to \underline{write your own} {\tt Matlab} code.
For the numerical experiments you should use your implementation and 
not  the available {\tt Matlab} function {\tt pcg}. The latter can be used
to check the correctness of your implementation.
\end{exercise}


\begin{exercise}[Generation of test data]

Generate four test matrices which have different properties:

\begin{tabular}{cl}
A:&{\tt A=matgen\_disco(s,1)};\\
B:&{\tt B=matgen\_disco(s,0.001);}\\
C:&{\tt C=matgen\_anisot(s,1,1);}\\
D:&{\tt D=matgen\_anisot(s,1,0.001);}\\
\end{tabular}

Here $s$ determines the size of the problem, namely, the matrix size is $n=2^s$.

Check the sparsity of the above matrices ({\tt spy}). Compute the complete spectra of $A, B, C, D$ (say, for $s=3$) 
by using the {\tt Matlab} function {\tt eig}. Plot those and compare. 
Compute the condition numbers of these matrices. 
You are advised to repeat the experiment for a larger $s$ to see how the condition 
number grows with the matrix size $n$.
\end{exercise}


\begin{exercise}[Numerical experiments and method comparisons]

Test the four methods M1, M2, M3, M4 on the  matrices $A, B, C, D$.

Choose a right-side vector as {\tt b=rand(n,1);} and 
compute the corresponding 'exact' solution as 
{\tt sol\_A=A$\backslash$b}. Recall, that $n=2^s$ is the size of $A$.

For methods M1, M3 and M4 use as a stopping criterion 
$\|\bx{r}_k\|/\|\bx{r}_0\|/ < 10^{-6}$.
For methods M3 and M4  $\bx{r}_k$ should be the iteratively computed 
(not the true residual, computed as $\bx{r}_k = \bx{b}-A\bx{x}_k$) residual.
For method M2, first predict the number of iterations in advance using formula \rref{eps_k}, do that many iterations 
and compare the obtained error and residual reduction with the expected ones.


Plot the corresponding residual convergence curve for the four methods.
For the CG methods plot both the iteratively computed residual and the true residual obtained as {\tt r\_true=b-A*x\_k}. 

Perform a series of experiments for $s=3,4,5,6$ to see the the dependence of the methods behaviour when $n$ is increased.
\end{exercise}

\begin{exercise}[A theoretical task] 
Derive that the convergence of M4 is better than that of M3.

Here the intention is to try to show it yourself or to do a literature search to find a 
proper result to cite. A good source to check could be
Anne Greenbaum, \textit{Iterative methods for solving linear systems}, 
SIAM, 1997.
\end{exercise}

\subsubsection*{Instructions for performing the numerical tests}
Download all the files from the course web-page.
The main files to call are the following:

\begin{itemize}
\item {\tt matgen\_disco.m}\\
The routine {\tt matgen\_disco.m} generates a finite element stiffness matrix for the Laplace equation 
\begin{equation}
 -\frac{\partial}{\partial x}
\left ( a\frac{\partial u}{\partial x} \right )
- \frac{\partial}{\partial y}
\left ( a\frac{\partial u}{\partial y} \right )
= f 
\label{problem2}
\end{equation}
in $\Omega \equiv (0,1)^{2}$
where $a$ is a constant and $a=\varepsilon\ll 1$ in a subset 
$\widetilde{\Omega}\subset \Omega$ and $a=1$ elsewhere with 
$\widetilde{\Omega}\equiv\{1/4\le x\le 3/4,\; 1/4\le y\le 3/4\}$.
The problem is then discretized using
regular isosceles triangles and piece-wise
linear basis functions. The mesh-size parameter $h$ is 
equal to $2^{-s}$ for some integer $s$ and is
chosen always such that 
there will be mesh-lines along the edges of $\widetilde{\Omega}$

The routine is called as
{\tt A=matgen\_disco(s,a);}. When {\tt a} is one, the generated matrix corresponds to  $-\Delta u = f$.

The above routine needs auxiliary files 
{\tt disco\_stiff.m}, {\tt disco\_rule.m}, {\tt Xdir.m} and {\tt Ydir.m}.

\item {\tt matgen\_anisot.m}\\
The routine {\tt matgen\_anisot.m} generates a finite difference (5-point)
discretization of the anisotropic Laplacian
$$- \varepsilon_x u_{xx} - \varepsilon_y u_{yy}.$$

The routine is called as
{\tt A=matgen\_anisot(s,epsx,epsy);}, where again mesh-size parameter $h$ is 
equal to $2^{-s}$.

If both {\tt epsx} and {\tt epsy} are equal to one, the generated matrix corresponds to  $-\Delta u = f$ discretized with central differences. 

\item {\tt Lanczos.m}\\
In order to use the Chebyshev iterative method one needs to estimate the 
extremal eigenvalues of $A$. The routine {\tt Lanczos} computes such
{\em approximations}. Note, however that there is no guarantee that we will obtain a lower bound for $\lambda_{min}$ and an upper bound for
$\lambda_{max}$!

The routine is called as
{\tt [lanmin,lanmax]=Lanczos(A);}
It performs internally a number of Lanczos steps until the following criteria are met:
$$|\lambda_1^{(k)}-\lambda_1^{(k-1)}|\le \epsilon \; \mbox{and} \;
|\lambda_n^{(k)}-\lambda_n^{(k-1)}|\le \epsilon $$
with a default value of $\epsilon=0.01$.


{\bf Remark:} You can estimate the extreme eigenvalues using some other technique, 
for example, using Gershgorin's theorem. 
In this case, you should include in the report a description how the bounds are found.

It is possible to use some other program (or program implementation of the Lanczos method) 
to compute approximations of the extreme eigenvalues of $A$.


\end{itemize}


\subsubsection*{Writing a report on the results}
The report has to have the following structure:

\begin{itemize}
\item[(i)] Brief problem description, namely, the solution of linear systems of equations 
with  symmetric positive definite matrices, using iterative solution methods.
\item[(ii)] Theory\\
Describe briefly the theory regarding the convergence of any of the methods M1-to M4,
including derivation of  Task 4.

\item[(iii)] Numerical experiments
\begin{itemize}
\item[(a)] Describe the experiments.
Include table(s) with iteration counts for various problem sizes.
Present some typical plots of the the relative residual 
($\|\bx{r}_k\|/\|\bx{r}_0\|$) convergence and the error convergence 
$\|\bx{x}^*-\bx{x}_k\|$. %You can use $\|\cdot \|$, for instance.

It is important to describe what is plotted on the different coordinate axes.
Another suggestion for the residual plots is to use {\tt semilogy} in order to better see 
the convergence history.

OBS! The numerical tests  possible to performed are numerous.
Do not include everything (or rather too much) in the report. Choose a representative information to show the typical behaviour and comment on the rest of the tests run, if necessary. For instance, do not run inefficient (slowly converging) methods on very large problems. The behaviour can be seen on small-to-medium sized problems. Still, one needs, say, three consecutive 
problem sizes to see how does the iteration count grow with the problem size.


\item[(b)] Analyse the numerical results in comparison with the theoretical. 

Do you see iteraion counts for the CG as predicted from the theoretical estimates? 
If the condition
number of the matrices is proportional to $h^{-2}$, how does the iteration count change
when $h$ is decreased?

Does the convergence of the Chebyshev method improve significantly if you use the exact
eigenvalues instead of their approximations obtained by the Lanczos method? For which
problems?

\item[(c)] (Not compulsory) Upon your time, experience and interests, you could try another
preconditioner, for example, the incomplete Cholesky preconditioner, produced by 
{\tt Matlab}'s command
{\tt  U = cholinc(A,tol);},
where {\tt tol} can be chosen for instance as $0.1, 0.05, 0.01$.
How much does the preconditioner you have chosen improve the convergence of the CG method? 
How does the number of iterations depend on the size of the matrix, 
compared to the unpreconditioned CG? 
\end{itemize}
\item[(iv)] Give your conclusions. Which is you method of choice? Motivate your answer.

\end{itemize}

Printouts of the program codes should be attached to the report.

Working in pairs is recommended. 
However, all topics in the assignment should be covered by each participant.

For your convenience, the Latex source of this assignment is at your disposal and you can use it for the report if you wish.


\bigskip
{\bf Deadline:}
The solutions should be delivered to me no later than {\bf December 10, 2007}.\\
Success!\\

Maya ({\tt Maya.Neytcheva@it.uu.se}, $\;$ room 2307)

\vfill
\rule{6.1in}{1pt}

Any comments on the assignment will be highly appreciated  
and will be considered for further improvements. Thank you!

\newpage
\subsubsection*{Appendix:
The Second Order Chebyshev iterative solution method}

{\bf Description}\\
%
Let $A$ be a real symmetric positive definite (s.p.d.) matrix of order $n$.
Consider the solution of the linear system $A\bx{x}=\bx{b}$
using the following iterative scheme, known as 
the {\em Second Order Chebyshev iterative solution method}:

\begin{equation}
\begin{array}{rl}
\bx{x}_0 \;\mbox{given}, &\;\bx{x}_1=\bx{x}_0 + {1\over 2}\beta_0\bx{r}_0\\
\mbox{For} &k=0,1,\cdots \;\mbox{until convergence}\\
&\bx{x}_{k+1} = \alpha_k\bx{x}_{k}+(1-\alpha_k)\bx{x}_{k-1}+\beta_k\bx{r}_{k}.
\\
&\bx{r}_k=\bx{b}-A\bx{x}_k.
\end{array}
\label{chebrecurence}
\end{equation}

Let $\bx{x}^*$ be the exact solution of the above linear system.
Denote by $\bx{e}_k=\bx{x}^*-\bx{x}_k$ the iterative error at step $k$.
Clearly, $\bx{r}_k=\bx{b}-A\bx{x}_k = A(\bx{x}^*-\bx{x}_k) = A\bx{e}_k$.

It is seen from the recursive formula \rref{chebrecurence} that for each $k$,
the errors satisfy a relation of the form
$$
\bx{e}_{k+1} = \alpha_k\bx{e}_{k}+(1-\alpha_k)\bx{e}_{k-1}+\beta_kA\bx{e}_{k}
= Q_k(A)\bx{e}_0,
$$ 
where $Q_k(\cdot)$ is some polynomial of degree $k$. Furthermore, the polynomials
$Q_k(A)$ are related among themselves as follows:
\begin{equation}
Q_{k+1}(A) - \alpha_kQ_k(A)-\beta_kAQ_k(A)+(1-\alpha_k)Q_{k-1}(A) = 0,\;
k=1,2,\cdots
\label{polynrec}
\end{equation}
We compare the recurrence \rref{polynrec} with the recursive formula for the Chebyshev polynomials, namely,
\begin{equation}
T_0(z)=1,\quad T_1(z)=z, \quad T_{k+1}(z)-2T_{k}(z)+T_{k-1}(z) = 0.
\label{polyncheb}
\end{equation}

One can easily see that for the following  special choice of the method parameters 
$$\alpha_k={2c\,T_k(z)\over T_{k+1}(z)}=1+{T_{k-1}(z)\over T_{k+1}(z)}$$
and 
$$\beta_k={4\over b-a}T_k(z)/T_{k+1}(z), \quad 
\mbox{where}\; z={b+a\over b-a},\quad
0<a\le \lambda_{min}(A), \; \lambda_{max}(A)\le b,
$$
we get
$$
Q_k(A) = {T_k(z)\over T_k(z)}, \qquad Z={1\over b-a}\left[(b+a)I-2A\right].
$$
We now recall that the Chebyshev polynomials possess the following optimal
approximation  property
- among all normalized polynomials of degree $n$ defined in an interval $[a,b]$, the $n$th degree Chebyshev polynomial is the one which differs from zero least
(measured in local min and max in $[a,b]$).

Therefore, for this particular polynomial $Q_k(A)$ we have (due the 
above approximation properties of the Chebyshev polynomials)
$$
\max\limits_z|Q_k(A)z|\le \min\limits_{P\in\Pi_k}|P(A)z|,
$$
i.e., at each step we achieve the best possible error reduction.
(Here $\Pi_k$ is the set of all polynomials of degree $k$.)

In order to use the Chebyshev iteration method we need to estimate the extreme eigenvalues of $A$ and to determine an interval $[a,b]$ which contains the spectrum of $A$.

Having done that, 
one finds the following formulas to compute the method parameters recursively:

\begin{equation}
\alpha_k = {a+b\over 2}\beta_k,\quad
{1\over \beta_k} = {a+b\over 2}-\left({b-a\over 4}\right)^2\beta_{k-1},\quad
\beta_0={4\over a+b}.
\label{cheb_coef}
\end{equation}
Note that $\alpha_k>1, k\ge 1$. 

\begin{theorem}
The following results hold:
$$
\lim\limits_{k\to\infty}\beta_k = {4\over \left(\sqrt{a}+\sqrt{b}\right)^2},
$$
$$
{\|\bx{e}_k\|_A\over \|\bx{e}_0\|_A} \le 
{1\over T_k\left({b+a\over b-a}\right)}
\le 2{\sigma^k\over 1+\sigma^{2k}}, \qquad
\sigma={1-\sqrt{a/b}\over 1+\sqrt{a/b}}.
$$
It follows then that
${\|\bx{e}_k\|_A\over \|\bx{e}_0\|_A}\to 0$ monotonically.

(For the special choice $a=\lambda_{min}(A)$, $b=\lambda_{max}(A)$, we have
$\sigma={1-\sqrt{\varkappa(A)^{-1}}\over 1+\sqrt{\varkappa(A)^{-1}}}$.)


From the above convergence rate estimate one can determine 
\underline{a priory} the number of Chebyshev iterations needed to be 
performed in order to achieve an error reduction
\begin{equation}
{\|\bx{e}_k\|_A\over \|\bx{e}_0\|_A}\le\varepsilon.
\label{eps}
\end{equation}

Indeed, to insure \rref{eps}, it suffices to perform
\begin{equation}
k\ge \ln\left({1\over \varepsilon}+\sqrt{{1\over\varepsilon^2}-1}\right) / 
\ln(\sigma^{-1})\quad
\mbox{or} \quad
k^* = \left\lceil{1\over 2}\sqrt{{b\over a}}\ln{2\over\varepsilon}\right\rceil
\quad\mbox{iterations}.
\label{eps_k}
\end{equation}

\end{theorem}


\end{document}
