\documentclass[11pt]{article}
\textwidth 6.8in
\textheight 9.2in
\oddsidemargin -0.3in
\evensidemargin -0.3in
\usepackage{latexsym}
\usepackage{multicol}
%\usepackage{epsfig}
\usepackage{latexsym}
 \usepackage{epsfrag}
\usepackage{multicol}
\usepackage{fancybox}
\usepackage{pdfpages}
\usepackage{graphicx}


 \topmargin -.7in

\input{zmacros}
\def\be{\begin{enumerate}}
\def\ee{\end{enumerate}}
%\newcommand{\vthin}{\vskip 5pt}


\pagestyle{empty}
\begin{document}
\hfill\vbox{\hbox{MA/CSC 580, Homework \#4}
\hbox{\bf Due, October 31, 2018}}

\vskip .2in
\baselineskip 14pt

\be

\item Let $A$ be a symmetric positive definite matrix,  so it has the Cholesky decomposition $A= L L^T$. Show that
(a):  $0< l_{kk} \le \sqrt{a_{kk}}$, $k=1,2,\cdots, n$. (b): From (a) to derive $\dsp \max_{1\le j\le i \le n} |l_{ij}| \le \sqrt{ \max_{1\le i,j \le n} |a_{ij}| }$. That is, the Cholesky decomposition is a stable algorithm.
(c): Do the Jacobi, Gauss-Seidel, and SOR($\omega)$ iterative methods converge?

\vthin
\item Consider the Poisson equation
\eqmno
u_{xx} + u_{yy}  &=& xy, \quad (x,y) \in \Omega \\
\left. u(x,y) \right |_{\pa \Omega} &=& 0,
\enmno
where $\Omega$ is the unit square. Using the finite difference method, we
can get a linear system of equations
\eqm
\fr{U_{i-1,j}+U_{i+1,j}+U_{i,j-1}+U_{i,j+1}-4 U_{ij}}{h^2} = f(x_i,y_j),
\quad 1 \le i, j \le 3,
\enm
where $h=1/4$, $x_i=i h$, $y_j=jh$, $i,j=0,1,2,3,4$,
and $U_{ij}$ is  an approximation of $u(x_i,y_j)$.
Write down the coefficient matrix and the right hand side using the
{\em red-black} orderings given in the right diagram below.
What is the dimension of the coefficient matrix? How many nonzero
entries and how many zeros? Generalize your results to general case when
$0 \le i,j \le n$ and $h=1/n$. Write down the component form of the
$SOR(\omega)$ iterative method. Does the $SOR(\omega)$ iterative method depend
on the ordering? From your analysis, explain whether you
prefer to use Gaussian elimination method or an iterative method.

%\begin{figure}[hbtp]
\bc
%\epsfysize=2.0in
   %\centerline{\hbox{\protect\epsffile{mesh.eps}}}
   
 % \includepdf[width=10.1cm]{mesh.pdf}
  \includegraphics[width=3.6in]{mesh.pdf}
\ec
 %
%\caption{Digram for Exercise 1.} \label{fig1}

%\end{figure}

\vthin
\item Given the following linear system of equations:
\eqmno
3 x_1 -x_2 + x_3 &=& 3 \\
    2 x_2 + x_3 &=& 2 \\
-x_2 + 2 x_3 &=& 2
\enmno
\be \item  With $x^{(0)} = [1, -1,1]^T$, find the {\em first}
iteration of the Jacobi, Gauss-Seidel, and SOR ($\omega=1.5$) methods.
\item  Write down the Jacobi,  Gauss-Seidel, and  Seidel iteration matrices
 $R_J$ ,  $R_{GS}$, and SOR($\omega$).
\item  Do the Jacobi and Gauss-Seidel iterative methods converge? Why?
\ee


\vthin
\item Explain when we want to use {\em iterative} methods to solve linear
system of equations $Ax=b$ instead of {\em direct} \  methods.

 Also if $||R||=1/10,\,$ then the iterative method
$\; x^{(k+1)}=R\,\dsp x^{(k)} + c\,\;$ converges to the solution
$x^*$, \linebreak  $x^* = R x^* +c $.
How many iterations are required so that $||x^{(k)} - x^* || \leq
10^{-6}$? \ \ Suppose $\;||x^{(0)} - x^* || = O(1)$.



\vthin
\vthin
\item Judge whether the iterative method $\; x^{(k+1)}=R\,\dsp x^{(k)} +
c\,\;$ converges or not.
\eqmno
(a): \quad R =\left [\begin {array}{ccccc} e^{-1}&-e^1&-1&-1&-10
\\\noalign{\medskip}0&\sin \pi/4 &10^4&-
1&-1\\\noalign{\medskip}0&0&-0.1&-1&1\\\noalign{\medskip}0&0&0&1-e^{-2}&-1
\\\noalign{\medskip}0&0&0&0&1-\sin(\alf \pi)\end {array}\right ], \hspace{1cm} (b): \quad
\left [\begin {array}{ccc}
0.9&0&0\\\noalign{\medskip}0&0.3&-0.7\\\noalign{\medskip}0&0.69&0.2999
\end {array}\right ]
\enmno

\vthin
\vthin
\item Determine the convergence of the Jacobi and Gauss-Seidel method
applied to the system of equations $Ax=b$, where
\eqmno
(a):  \quad  A = \left [\begin {array}{ccc}
0.9&0&0\\\noalign{\medskip}0&1&2\\\noalign{\medskip}0&-2&1\end {array}\right
], \qquad (b): \quad A = \left [\begin {array}{ccccccc}
3 & -1 & 0   & 0 & \cdots &\cdots & 0 \\
2 & 3 & -1 & 0 & \cdots & \cdots& 0 \\
0 & 2 & 3 & -1 & \cdots &\cdots&0 \\
\cdots  &\cdots     & \cdots&\cdots&\cdots &\cdots & \cdots  \\
\cdots  &\cdots     & \cdots&\cdots&\cdots &\cdots & \cdots  \\
 0 &\cdots     & \cdots& 0 & 2 & 3 & -1 \\
0 &\cdots     & \cdots& \cdots& 0 &2 & 3  \\
\end {array}\right ]
\enmno

\vthin
\vthin
\item Modify the Matlab code {\em poisson\_drive.m} and {\em poisson\_sor.m}
to solve the following diffusion and convection equation:
\eqmno
u_{xx} + u_{yy} + a u_x - b u_y = f(x,y), \quad 0 \le x,y \le 1,
\enmno
Assume that solution at the boundary $x=0$, $x=1$, $y=0$, $y=1$ are given
(Dirichlet boundary conditions).  The central-upwinding
finite difference scheme is
\eqmno
\fr{U_{i-1,j} + U_{i+1,j} + U_{i,j-1} + U_{i,j+1}- 4 U_{ij} }{h^2}
+ a \fr{U_{i+1,j} - U_{ij} }{h} - b \fr{U_{i,j} - U_{i,j-1}}{h} = f_{ij}
\enmno
\be
\item Assume the exact solution is $u(x,y)=e^{2 y} \sin (\pi x)$,
find $f(x,y)$.
\vthin
\item Use the $u(x,y)$ above for the boundary condition and the $f(x,y)$
above for the partial differential equation.
Let $a=1$, $b=2$, and $a=100$, $b=2$, solve the problem
with $n=20$, $40$, $80$, and $n=160$.
Try $\omega=1$, the best $\omega$ for the Poisson equation discussed in the
class, the optimal $\omega$ by testing, for example $\omega=1.9, 1.8,
\cdots, 1$.
\item Tabulate the error, the number of iterations for $n=20$, $40$, $80$,
and $n=160$ with your tested optimal $\omega$, compare the number of
iterations with the Gauss-Seidel method.
\item Plot the solution and the error for $n=40$ with  your tested optimal
$\omega$. Label your plots as well.
\ee

\item {\bf Extra credit}.  Do some research to explain the  behavior using cond\_hw.m and possible ways of improving it. 

\ignore{
\vthin
\vthin
\item {\bf Extra Credit:} Consider the finite difference equations for 1D
Poisson equation
\eqmno
-u''(x) = f(x), \qquad u(0)=u(1)=0, \quad 0 < x < 1.
\enmno
\eqm
\fr{-u_{i-1} + 2 u_i - u_{i+1}}{h^2} = f(x_i), \quad
i=1,2,\cdots, n-1, n, \quad \mbox{with} \quad
u_0 = 0, \; u_{n+1} =0, \; h=\fr{1}{n+1}.
\enm
Write the system of equations as matrix-vector form $A {\bf u} = {\bf F}$,
where $A$ is an $n$ by $n$ symmetric positive definite (S.P.D) matrix.
\be
\item Show that $n$ eigenvalues of $A$ are given by
\eqm
\lambda_i = \fr{2}{h^2} \lp 1- \cos(i \pi h) \rp, \quad i=1,2,\cdots, n.
\enm
The eigenvector ${\bf u}^i$ corresponding to $\lambda_i$ is
\eqmno
&& \left [\begin {array}{ccccccc} \sin(i\pi h)& \sin(2 i\pi h)&
\cdots& \sin(ij\pi h) & \cdots& \sin( (n-1)i\pi h) &
\sin( n i\pi h) \end {array}\right ]^T, \\
&& \mbox{or in simpler notation} \quad {\bf u}^i = \sin( i \pi jh), \quad
j=1,2,\cdots n.
\enmno
\item Form (a), show that $cond_2(A)= \|A\|_2 \|A^{-1}\|_2 \sim
O(1/h^2)=O(n^2)$ as $n$ is large.
\item  Estimate the spectral radius of iterative matrices for Jacobi,
Gauss-Seidel, and SOR methods.
\item Find the optimal $\omega$ for the SOR iterative method.
\ee
{\bf Hint:} Refer to Sections 6.3.2, 6.3.3, 6.5.4, 6.5.5.

}
\ee

\end{document}

\newpage

{\bf Selected solutions}

\bi
 \item[1] The coefficient matrix of the red-black ordering is
 \eqmno
  A = \left(%
   \begin{array}{cc}
  D_1 & B \\
  B^T & D_2 \\
 \end{array}%
\right) \qquad D_i = - \frac{4}{h^2}\, I, \quad i=1, 2
 \enmno
$B$ is a sparse matrix in which each row has at most 4 non-zero
entries. The size of the matrix is $(n-1)^2$ by $(n-1)^2$ with
non-zero entries of $O(5 n^2)$. The $SOR(\omega)$ at $k$-th
iteration can be written as
\eqmno
   &&  u_{ij}^{k+1}=u_{ij}^{k}, \quad i,j=1,2,\cdots, n-1. \\ \eqsp
   && u_{ij}^{k+1}= (1-\omega)u_{ij}^{k} + \frac{\omega}{4}
   \left ( u_{i-1,j}^{k+1} + u_{i+1,j}^{k+1} + u_{i,j-1}^{k+1} +
   u_{i,j+1}^{k+1} - h^2  f_{ij} \right ), \quad i,j=1,2,\cdots, n-1.
\enmno
The iterative method does not depend on the ordering of the
equations and unknowns, but does depend on the index $i$ and $j$.

\item[2] (a) The results of one iteration are
 \eqmno
 \mathbf{x}_J\left(%
\begin{array}{c}
  {1}/{3} \\
  {1}/{2} \\
  {1}/{2} \\
\end{array}%
\right), \qquad \mathbf{x}_{GS} = \left(%
\begin{array}{c}
  {1}/{3} \\
  {1}/{2} \\
  {5}/{4} \\
\end{array}%
\right) \qquad \mathbf{x}_{SOR(1.5)} = \left(%
\begin{array}{c}
  0 \\
  {5}/{4} \\
  {31}/{16} \\
\end{array}%
\right)
 \enmno
(b) and (c), The iteration matrices are:
 \eqmno
R_J = \left [\begin {array}{ccc}
0&1/3&-1/3\\\noalign{\medskip}0&0&-1/2\\\noalign{\medskip}0&1/2&0\end
{array}\right ], \qquad R_{GS} = \left [\begin {array}{ccc}
0&1/3&-1/3\\\noalign{\medskip}0&0&-1/2\\\noalign{\medskip}0&0&-1/4\end
{array}\right ].
  \enmno
  Since $\|R_J\|_{\infty}=2/3 < 1$ and $\|R_{GS}\|_{\infty}=2/3 <
  1$, both iterative methods converge.


\item[3] For the first matrix, the eigenvalues of $R$ are the
diagonals. Notice that $|a_{ii}|<1$ for $i=1,2,3,4$. We just need
to check $|a_{55}|=1-\sin (\alpha \pi)$. Note that $0< \sin x < 1$
if $0 < x < \pi$ and $\sin x $ is a periodic function of $2 \pi$.
Thus, if $ 2k < \alpha < 2k +1 $, then the iterative method
converges, where $k$ is an integer.

For the second matrix we have $\|R\|_{1}= 0.9999<1$, the
iterative method converges.

\item[4] (a) The iteration matrices are:
 \eqmno
 R_J = \left [\begin {array}{ccc} 0&0&0\\\noalign{\medskip}0&0&-2\\
 \noalign{\medskip}0&2&0\end {array}\right ], \qquad R_{G-S} =
  \left [\begin {array}{ccc} 0&0&0\\\noalign{\medskip}0&0&-2\\
  \noalign{\medskip}0&0&-4\end {array}\right ].
 \enmno
Since $\rho(R_J) = 2 > 1$ and $\rho(R_{GS}) = 4  >
  1$, both iterative methods diverge.

  (b) The matrix is weakly diagonally dominant and irreducible.
  Both Jacobi and Gauss-Seidel iterative methods converge.

  \ei
\end{document}
