\documentclass[12pt,a4paper]{article}

\usepackage{amstext}
\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{stmaryrd}
\usepackage{enumerate}
\usepackage{alltt}
\usepackage[latin1]{inputenc}
\usepackage{cancel}
\usepackage{fancyhdr}
\usepackage{listings}
\usepackage{pstricks}
\usepackage{pst-plot}
\usepackage{wasysym}
\usepackage{rotating}

\title{\textsc{Numerical Algorithms for Visual Computing II \\ \large Assignment 5}}
\author{\shortstack[c]{Jan Hendrik Dithmar\\ \scriptsize 2031259} \and \shortstack[c]{Andreas Steinel \\ \scriptsize
2029984}}
\date{}

\pagestyle{fancy}

\lhead{Assignment 5}
\chead{}
\rhead{Dithmar (2031259), Steinel (2029984)}
\lfoot{NAVC II}
\cfoot{Page \thepage}
\rfoot{WS 2008/2009}

\begin{document}

\renewcommand{\headrulewidth}{0.4pt}
\renewcommand{\footrulewidth}{0.4pt}

\newcommand{\exercise}[1]{\subsection{#1}}
\renewcommand{\theenumi}{\textbf{\alph{enumi}}}

\newcommand{\AND}{\ensuremath{\wedge}}
\newcommand{\OR}{\ensuremath{\vee}}
\newcommand{\SetR}{\ensuremath{\mathbb{R}}}
\newcommand{\SetN}{\ensuremath{\mathbb{N}}}
\newcommand{\SetC}{\ensuremath{\mathbb{C}}}
\newcommand{\SetZ}{\ensuremath{\mathbb{Z}}}
\newcommand{\SetQ}{\ensuremath{\mathbb{Q}}}
\newcommand{\SetB}{\ensuremath{\mathbb{B}}}
\newcommand{\QED}{\ensuremath{\hspace*{\fill}\square}}
\newcommand{\FLASH}{\ensuremath{\hspace*{\fill}\lightning}}
\newcommand{\OnotO}{\ensuremath{\mathcal{O}}}
\newcommand{\bild}[2]{
	\begin{figure}[htpb]
		\centering
			\includegraphics[#1]{#2}\\[0.2cm]
	\end{figure}
	}

\newcommand{\m}[1]{\begin{pmatrix}#1\end{pmatrix}}
\newcommand{\+}{\phantom{-}}
\newcommand{\f}{\dfrac}

\maketitle

\thispagestyle{empty}

\newpage
\setcounter{page}{1}
\setcounter{section}{5}

\exercise{How good is trivial?}
\begin{enumerate}
	\item \textit{linear:} We have $N = I$ and $I = diag(1,1)$.
		\begin{eqnarray*}
			\Phi(x,b) & = & Mx + Nb \\
			\Rightarrow \Phi(x,b) & = & M \left(\begin{array}{c}x_1 \\ x_2\end{array}\right) + \left(\begin{array}{cc}1 & 0 \\ 0 & 1\end{array}\right) \left(\begin{array}{c}0.3 \\ 0.3\end{array}\right) \\
			& = & M \left(\begin{array}{c}x_1 \\ x_2\end{array}\right) + \left(\begin{array}{c}0.3 \\ 0.3\end{array}\right)
		\end{eqnarray*}
		\textit{consistent:}
		\begin{eqnarray*}
			M & = & I - NA = I - IA = I - \left(\begin{array}{rr}0.7 & -0.4 \\ -0.2 & 0.5\end{array}\right) \\
			& = & \left(\begin{array}{cc}0.3 & 0.4 \\ 0.2 & 0.5\end{array}\right)
		\end{eqnarray*}
		$$\Rightarrow \Phi(x,b) = \left(\begin{array}{cc}0.3 & 0.4 \\ 0.2 & 0.5\end{array}\right) \left(\begin{array}{c}x_1 \\ x_2\end{array}\right) + \left(\begin{array}{c}0.3 \\ 0.3\end{array}\right)$$
	\item Eigenvalues of the arising iteration matrix
		$$M = \left(\begin{array}{cc}0.3 & 0.4 \\ 0.2 & 0.5\end{array}\right)$$
		have to be computed first: $Mx = \lambda x \Leftrightarrow (M - \lambda I)x \stackrel{!}{=} 0$
		\begin{eqnarray*}
			det(M - \lambda I) & = & 0 \\
			\left|\begin{array}{cc}0.3 - \lambda & 0.4 \\ 0.2 & 0.5 - \lambda\end{array}\right| & = & (0.3 - \lambda)(0.5 - \lambda) - 0.08 \stackrel{!}{=} 0 \\
			\Leftrightarrow \lambda^2 - 0.8 \lambda + 0.07 & = & 0 \\
			\lambda_{1,2} & = & \dfrac{0.8 \pm \sqrt{0.8^2 - 4 \cdot 0.07}}{2} \\
			\lambda_1 & = & 0.7 \\
			\lambda_2 & = & 0.1
		\end{eqnarray*}
		We have the set $\sigma(M) = \{0.1, 0.7\}$. This means that $\rho(M) = \max \{|\lambda| : \lambda \in \sigma(M)\} = 0.7$. This gives the spectral radius 0.7.
	\item \textit{Scilab.}
\end{enumerate}

\exercise{How good is Jacobi?}
\begin{enumerate}
	\item Let's define $D := diag(0.7, 0.5)$. It follows that $D^{-1} = diag\left(\dfrac{10}{7}, 2\right)$. The iteration matrix $M = D^{-1} (D - A)$ looks as follows:
		\begin{eqnarray*}
			M & = & \left(\begin{array}{cc}\dfrac{10}{7} & 0 \\[0.3cm] 0 & 2\end{array}\right) \left(\begin{array}{cc}0 & 0.4 \\ 0.2 & 0\end{array}\right) \\
			M & = & \left(\begin{array}{cc}0 & \dfrac{4}{7} \\[0.3cm] \dfrac{2}{5} & 0\end{array}\right)
		\end{eqnarray*}
		The eigenvalues of $M$ are:
		\begin{eqnarray*}
			det(M - \lambda I) & \stackrel{!}{=} & 0 \\
			\Leftrightarrow \left|\begin{array}{cc}0 - \lambda & \dfrac{4}{7} \\[0.3cm] \dfrac{2}{5} & 0 - \lambda\end{array}\right| & = & 0 \\
			\Leftrightarrow (0 - \lambda)^2 - \dfrac{8}{35} & = & 0 \\
			\lambda^2 - \dfrac{8}{35} & = & 0 \\
			\lambda_1 & = & \dfrac{2 \sqrt{70}}{35} \\
			\lambda_2 & = & -\dfrac{2 \sqrt{70}}{35}
		\end{eqnarray*}
		This means that $\rho(M) = \dfrac{2\sqrt{70}}{35} \approx 0.4781$.
	\item \textit{Scilab.}
	\item We see that the Jacobi method already converges after 30 instead of 50 iterations for the trivial case. So it needs a bit more than half of the iterations than the trivial method.
\end{enumerate}

\exercise{How good is Gau\ss-Seidel?}
\begin{enumerate}
	\item The Gau\ss-Seidel scheme decomposes $A$ into $A = L + D + U$ with
		$$L = \left(\begin{array}{rr}0 & 0 \\[0.3cm] -\dfrac{1}{5} & 0\end{array}\right) \quad D = \left(\begin{array}{cc}\dfrac{7}{10} & 0 \\[0.3cm] 0 & \dfrac{1}{2}\end{array}\right) \quad U = \left(\begin{array}{rr}0 & -\dfrac{2}{5} \\[0.3cm] 0 & 0\end{array}\right)$$
		We know that
		$$M = -(D + L)^{-1} U.$$
		Let's calculate the different components:
		\begin{eqnarray*}
			D + L & = & \left(\begin{array}{rr}\dfrac{7}{10} & 0 \\[0.3cm] -\dfrac{1}{5} & \dfrac{1}{2}\end{array}\right) \\
			-(D + L)^{-1} & = & \left(\begin{array}{cc}\dfrac{10}{7} & 0 \\[0.3cm] \dfrac{4}{7} & 2\end{array}\right) \\
			M & = & -(D + L)^{-1} U = -\left(\begin{array}{cc}\dfrac{10}{7} & 0 \\[0.3cm] \dfrac{4}{7} & 2\end{array}\right) \left(\begin{array}{rr}0 & -\dfrac{2}{5} \\[0.3cm] 0 & 0\end{array}\right) \\
%			& = & -\left(\begin{array}{rr}0 & -\dfrac{4}{7} \\[0.3cm] 0 & -\dfrac{8}{35}\end{array}\right) \\
			& = & \left(\begin{array}{cc}0 & \dfrac{4}{7} \\[0.3cm] 0 & \dfrac{8}{35}\end{array}\right)
		\end{eqnarray*}
		Eigenvalues of $M$:
		\begin{eqnarray*}
			det(M - \lambda I) & \stackrel{!}{=} 0 \\
			\Leftrightarrow \left|\begin{array}{cc}0 - \lambda & \dfrac{4}{7} \\[0.3cm] 0 & \dfrac{8}{35} - \lambda\end{array}\right| & = & 0 \\
			\Leftrightarrow (0 - \lambda) \left(\dfrac{8}{35} - \lambda\right) & = & 0 \\
			\lambda^2 - \dfrac{8}{35} \lambda & = & 0 \\
			\lambda_1 & = & 0 \\
			\lambda_2 & = & \dfrac{8}{35}
		\end{eqnarray*}
		This means that $\rho(M) = \dfrac{8}{35} = 0.2286$.
	\item \textit{Scilab.}
	\item The Gau\ss-Seidel method needs 20 iterations to converge which means that it needs 10 iterations less in comparison to the Jacobi method.
\end{enumerate}

\exercise{How good is SOR?}
\begin{enumerate}
	\item To compute the optimal relaxation parameter, let's consider
		$$M_J := D^{-1} (D-A) = \left(\begin{array}{cc}\dfrac{10}{7} & 0 \\[0.3cm] 0 & 2\end{array}\right) \left(\begin{array}{cc}0 & \dfrac{2}{5} \\[0.3cm] \dfrac{1}{5} & 0\end{array}\right) = \left(\begin{array}{cc}0 & \dfrac{4}{7} \\[0.3cm] \dfrac{2}{5} & 0\end{array}\right)$$
		The eigenvalues of $M_J$ are:
		\begin{eqnarray*}
			det(M_J - \lambda I) & \stackrel{!}{=} & 0 \\
			\left|\begin{array}{cc}0 - \lambda & \dfrac{4}{7} \\[0.3cm] \dfrac{2}{5} & 0 - \lambda\end{array}\right| & = & 0 \\
			\Leftrightarrow (0-\lambda)^2 - \dfrac{8}{35} & = & 0 \\
			\lambda^2 - \dfrac{8}{35} & = & 0 \\
			\lambda_1 & = & - \sqrt{\dfrac{8}{35}} = - \dfrac{2 \sqrt{70}}{35} \\
			\lambda_2 & = & \sqrt{\dfrac{8}{35}} = \dfrac{2 \sqrt{70}}{35}
		\end{eqnarray*}
		From this, it follows that $\rho := \rho(M_J) = \sqrt{\dfrac{8}{35}}$. According to (9.36), we get
		\begin{eqnarray*}
			\omega_{opt} & = & \dfrac{2}{1 + \sqrt{1 - \rho^2}} \\
			& = & \dfrac{2}{1 + \sqrt{1 - \dfrac{8}{35}}} \\
			& = & \dfrac{2}{1 + \dfrac{3 \sqrt{105}}{35}} \\
			& = & \dfrac{2}{\dfrac{35 + 3 \sqrt{105}}{35}} \\
			& = & \dfrac{70}{35 + 3 \sqrt{105}} \\
			& \approx & 1.065
		\end{eqnarray*}
	\item Now, let's compute the iteration matrix $M$ (using $\omega = \omega_{opt}$):
		$$M = (D + \omega L)^{-1} \left[(1 - \omega) D - \omega U\right]$$
		where $L$, $D$ and $U$ are known from Exercise 3.
		\begin{eqnarray*}
			D + \omega L & = & \left(\begin{array}{cc}\dfrac{7}{10} & 0 \\[0.3cm] 0 & \dfrac{1}{2}\end{array}\right) + \omega \left(\begin{array}{cc}0 & 0 \\[0.3cm] -\dfrac{1}{5} & 0\end{array}\right) \\
			& = & \left(\begin{array}{cc}\dfrac{7}{10} & 0 \\[0.3cm] - \dfrac{14}{35 + 3\sqrt{105}} & \dfrac{1}{2}\end{array}\right) \\
			(D + \omega L)^{-1} & = & \left(\begin{array}{cc}\dfrac{10}{7} & 0 \\[0.3cm] \dfrac{40}{35 + 3\sqrt{105}} & 2\end{array}\right)
		\end{eqnarray*}
		\begin{eqnarray*}
			(1 - \omega) D & = & \left(\begin{array}{cc}\dfrac{7}{10} - \dfrac{49}{35 + 3 \sqrt{105}} & 0 \\[0.3cm] 0 & \dfrac{1}{2} - \dfrac{35}{35 + 3 \sqrt{105}}\end{array}\right) \\
			\omega U & = & \left(\begin{array}{cc}0 & - \dfrac{28}{35 + 3\sqrt{105}} \\[0.3cm] 0 & 0\end{array}\right)
		\end{eqnarray*}
		This gives the final iteration matrix $M$:
		$$M = \left(\begin{array}{cc}\dfrac{-35 + 3\sqrt{105}}{35 + 3\sqrt{105}} & \dfrac{40}{35 + 3\sqrt{105}} \\[0.3cm] \dfrac{28 (-35 + 3\sqrt{105})}{(35 + 3\sqrt{105})^2} & \dfrac{840}{(35 + 3 \sqrt{105})^2}\end{array}\right) \approx \left(\begin{array}{cc}-0.0648 & 0.6084 \\ -0.0276 & 0.1944\end{array}\right)$$
		Let's compute the eigenvalues:
		\begin{eqnarray*}
			\left|\begin{array}{cc}-0.0648 - \lambda & 0.6084 \\ -0.0276 & 0.1944 - \lambda\end{array}\right| & \stackrel{!}{=} & 0 \\
			(-0.0648 - \lambda) (0.1944 - \lambda) - 0.6084 \cdot (-0.0276) & = & 0 \\
			-0.01259712 + 0.0648 \lambda - 0.1944 \lambda + \lambda^2 + 0.01679184 & = & 0 \\
			\lambda^2 - 0.1296 \lambda + 0.00419472 & = & 0
		\end{eqnarray*}
		\begin{eqnarray*}
			\lambda_{1,2} & = & \dfrac{0.1296 \pm \sqrt{0.1296^2 - 4 \cdot 0.00419472}}{2} \\
			\lambda_1 & = & 0.0669 \\
			\lambda_2 & = & 0.0627
		\end{eqnarray*}
		It follows that $\rho(M) = 0.0669$.
	\item \textit{Scilab.}
	\item \textit{Scilab.}
	\item SOR needs 10 iterations which means again that it needs 10 iterations less to converge in comparison to Gau\ss-Seidel.
\end{enumerate}

\end{document}