\documentclass[11pt]{article}
\usepackage{epsfig}
\usepackage{graphicx}
\usepackage{color}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{xypic}
\usepackage{subcaption}
\usepackage{listings}
\lstset{language=Matlab}

\setlength{\textwidth}{6.5in}
\setlength{\topmargin}{0in}
\setlength{\headsep}{0in}
\setlength{\headheight}{0in}
\setlength{\oddsidemargin}{0in}
\setlength{\textheight}{9in}

\newcommand{\todo}[1]{\textcolor{red}{[#1]}}

\input xy
\xyoption{all}

\providecommand{\SolutionQuestionLULDMPartA}{}
\providecommand{\SolutionQuestionLULDMPartB}{}
\providecommand{\SolutionQuestionLULDMPartC}{}
\providecommand{\SolutionQuestionFlopCount}{}
\providecommand{\SolutionQuestionCholeskyPartA}{}
\providecommand{\SolutionQuestionCholeskyPartB}{}
\providecommand{\SolutionQuestionPositiveDefiniteNotSymmetricPartA}{}
\providecommand{\SolutionQuestionPositiveDefiniteNotSymmetricPartB}{}
\providecommand{\SolutionQuestionDifferenceEquationsPartA}{}
\providecommand{\SolutionQuestionDifferenceEquationsPartB}{}
\providecommand{\SolutionQuestionDifferenceEquationsPartC}{}

%Used to display the marks available in the left margin.
\newcommand{\rmp}{\reversemarginpar\marginpar}

\pagestyle{empty}

\begin{document}

\begin{center}
{\Large {\bf CS 475/675 Spring 2025: Crowdmark Assignment 2}}\\
\end{center}

\noindent
{\bf Due May 30, at 11:59 pm Eastern.} \\
Submit all components of your solutions (written/analytical work, code/scripts, figures,
plots, output, etc.) to CrowdMark in PDF form in the appropriate section for each question.\\

\noindent You must also separately submit a single zip file containing any and all code/scripts you
write to the Crowdmark Assignment 2 DropBox on LEARN, in runnable format (that is, .m).\\

\noindent
\textbf{For full marks, be sure to show all of your work!}

\begin{enumerate}

\item \label{QuestionLULDM}
Let $A = \begin{bmatrix} 2 & 4 & 10 \\ -2 & -5 & -8 \\ 4 & 6 & 27 \end{bmatrix}$.
\begin{enumerate}
\item \label{QuestionLULDMPartA}
Compute \rmp{[6]}the $LU$ factorization of $A$, in which $L$ is unit lower triangular.
Do this computation by hand, not in Matlab.
\SolutionQuestionLULDMPartA

\newpage

\item \label{QuestionLULDMPartB}
Compute \rmp{[3]}the $LDM^{T}$ factorization of $A$.
Do this computation by hand, not in Matlab.
You should use your answer from part \ref{QuestionLULDMPartA}.
\SolutionQuestionLULDMPartB

\newpage

\item \label{QuestionLULDMPartC}
Suppose \rmp{[6]}that an arbitrary square matrix, $B$, has two $LDM^{T}$ factorizations, i.e. suppose that
\begin{displaymath}
L_{1}D_{1}(M_{1})^{T} = B = L_{2}D_{2}(M_{2})^{T} ,
\end{displaymath}
where all of $L_{1}, L_{2}, M_{1}, M_{2}$ are unit lower triangular, and $D_{1}, D_{2}$ are diagonal, with all their diagonal entries non-zero.
Prove that $L_{1} = L_{2}, D_{1} = D_{2}$ and $M_{1} = M_{2}$.
(\textbf{Remark:} This proves that the $LDM^{T}$ factorization is unique.)

\noindent
\textbf{Hint:} Use the theorem that says that if $B$ has two $LU$ factorizations,
\begin{displaymath}
L_{1}U_{1} = B = L_{2}U_{2} ,
\end{displaymath}
in which $L_{1}, L_{2}$ are unit lower triangular and $U_{1}, U_{2}$ are upper triangular, with all their diagonal entries non-zero, then $L_{1} = L_{2}$ and $U_{1} = U_{2}$ (i.e. the $LU$ factorization, in which the $L$ is unit lower triangular and the $U$ has non-zero diagonal entries, is unique).
\SolutionQuestionLULDMPartC
\end{enumerate}

\newpage

\item
The \rmp{[8]}following algorithm, to carry out Gaussian Elimination on an $n \times n$ matrix $A$, is reproduced from the Lecture Notes and Slides.\\

\noindent
\textbf{GE Algorithm}
\begin{tabbing}
for $i$ = 1, 2, ..., n-1\\
\hspace{5mm}for $k$ = i+1, ..., n\\
\hspace{10mm}mult = $a_{ki}$ / $a_{ii}$ \hspace{20mm} \= \\
\hspace{10mm}$a_{ki} = 0$  \> not needed, but helpful for intuition\\
\hspace{10mm}for $j$ = i+1, ..., n\\
\hspace{15mm}$a_{kj}$ = $a_{kj}-$ mult $* a_{ij}$ \> update row k\\
\hspace{10mm}end\\
\hspace{10mm}$b_{k}$ = $b_{k}-$ mult $* b_{i}$ \> update RHS\\
\hspace{5mm}end\\
end
\end{tabbing}
Compute the count of the flops required to execute this algorithm, in terms of $n$.

\noindent
You may use, without proof, the sums:
\begin{eqnarray*}
\sum_{i=1}^{n} i & = & \frac{n(n+1)}{2} \textrm{, and}\\
\sum_{i=1}^{n} i^{2} & = & \frac{n(n+1)(2n+1)}{6} .
\end{eqnarray*}
\SolutionQuestionFlopCount

\newpage

\item 
Let $A = \begin{bmatrix} 9 & -12 & 3 \\ -12 & 17 & -4 \\ 3 & -4 & 5 \end{bmatrix}$.
It is clear that $A$ is symmetric.
\begin{enumerate}
\item 
Verify \rmp{[4]}that $A$ is positive definite.
You may use Matlab to do this.
Whichever way you do your verification, show all of your work.
\SolutionQuestionCholeskyPartA

\newpage

\item
Compute \rmp{[8]}the Cholesky factor, $G$, of $A$.
Do this computation by hand, not in Matlab.
\SolutionQuestionCholeskyPartB
\end{enumerate}

\newpage

\item
Let $A = \begin{bmatrix} 1 & -1 \\ 1 & 1 \end{bmatrix}$.
It is clear that $A$ is not symmetric.
\begin{enumerate}
\item 
Prove \rmp{[4]}that, for any $\vec{0} \neq x \in \mathbb{R}^{2}, \ x^{T}Ax > 0$.
\SolutionQuestionPositiveDefiniteNotSymmetricPartA

\newpage

\item
Prove \rmp{[4]}that there cannot exist any lower triangular $2 \times 2$ real matrix $G$, with strictly positive entries on its diagonal, such that $G G^{T} = A$.
\SolutionQuestionPositiveDefiniteNotSymmetricPartB
\end{enumerate}

\newpage

\item 
As in Lecture 04, in modelling (steady state) heat flow, the temperature $T = T(x,y)$ satisfies the Poisson equation
\begin{equation*}
-\frac{\partial^2T}{\partial x^2} - \frac{\partial^2T}{\partial y^2} = f(x,y)
\end{equation*}
where $f(x,y)$ represents the heat source function.

\noindent
We approximate $T(x,y)$ at discrete locations on a two dimensional grid with $m$ active grid points in the horizontal dimension and $2m+1$ active grid points in the vertical dimension (imagine a rectangular metal plate which is 1 unit wide and 2 units tall).
Let the $(i,j)$ grid point have location
$(x_i,y_j)$ where $x_i=ih$, $y_j=jh$, $h= \frac{1}{m+1}$ is the grid spacing, $0 \leq i \leq m + 1$ and $0 \leq j \leq 2m+2$.
If we let $T_{i,j} \approx T(x_i,y_j)$, then the finite difference approximation results in a
set of linear equations
\begin{equation}
\frac{1}{h^2} \left ( 4T_{i,j} - T_{i-1,j}- T_{i+1,j}- T_{i,j-1}- T_{i,j+1} \right ) = f_{i,j}.
\label{eq:fdstencil}
\end{equation}
We will assume that all the boundary temperatures along the sides of the grid are zero:
\begin{displaymath}
T_{0,j}=T_{i,0}=T_{m+1,j}=T_{i,2m+2}=0
\end{displaymath}
for all $i,j$ as above.
We want to analyze the heat flow with two heat sources as given by:
\begin{equation*}
f_{i,j} = \left\{
     \begin{array}{lr}
       1 \text{ if } \|(x_i,y_j) - (0.3,1.75)\| \leq 0.1\\
       2 \text{ if } \|(x_i,y_j) - (0.7,0.25)\| \leq 0.2\\
       0 \text{ otherwise.}
     \end{array}
   \right.
\end{equation*}

\begin{figure}[h]
    \centering
    \includegraphics[width=0.5\textwidth]{heat_results.pdf}
    \caption{Example plot for three heat sources}
    \label{fig:heatfig}
\end{figure}

An example of a final temperature distribution for a \textbf{different} set of three heat sources is shown in Figure \ref{fig:heatfig}.

Define a vector $x$ such that
\begin{displaymath}
x_k=T_{i,j},
\end{displaymath}
where $k=i+(j-1)m$.
Similarly, define the vector $b$ with $b_k = f_{i,j}$.
Then we can write equation (\ref{eq:fdstencil}) in matrix form as
\begin{equation}
Ax=b.
\label{eq:axb}
\end{equation}
The coefficient matrix $A$ is sparse with size $n \times n$, where $n = m (2m+1)$.
\begin{enumerate}

\item
Create \rmp{[10]}a Matlab function \texttt{[A,b] = Lap2D(m)}.
The input is a positive integer $m$, the number of active grid points in the horizontal dimension, as defined above.
The outputs are the matrix $A$ and the vector $b$ as in equation (\ref{eq:axb}).
\SolutionQuestionDifferenceEquationsPartA

\newpage

\item
Implement \rmp{[18]}the following numerical methods: Gaussian elimination, Cholesky factorization, and Banded Gaussian elimination.
Create the following Matlab functions:
 \begin{align*}
&\texttt{ x=GaussElim(A,b)}\\
&\texttt{ x=Cholesky(A,b)}\\
&\texttt{ x=BandGE(A,b,p,q)}
 \end{align*}
These Matlab functions take as inputs the matrix $A$ and right-hand side $b$ and compute the solution $x$ using the corresponding variant of Gaussian elimination.
The parameters \texttt{p} and \texttt{q} indicate the lower and upper matrix bandwidth, respectively.
(There is no need to implement pivoting or check for zero or small pivots.)
You can check your code by visualizing the solution $x$.
Use the following command:
\begin{equation*}
\texttt{>> mesh(reshape(x,m,2*m+1))}
\end{equation*}
This will convert the solution vector $x$ into a 2D array, and generate a 2D mesh plot of the result.
Submit a 2D mesh plot of the solution for the case when $m=20$.
\SolutionQuestionDifferenceEquationsPartB

\newpage

\item
Create \rmp{[12]}a Matlab script, \texttt{GETimes.m}, that solves (\ref{eq:axb}) using the three different variants of Gaussian elimination you wrote.
In particular, set up the matrix $A$ and right-hand side $b$ by calling the Matlab function \texttt{Lap2D}.
Then solve the equation by calling one of the Matlab functions in part (b).
Record the execution time using the Matlab commands \texttt{tic} and \texttt{toc}.
Construct a table of execution times for solving with each of the methods above.
Try the values $m = 8, 16, 24, 32$.
(Consider vectorizing the innermost loop of your implementations to keep the solve times more manageable.)
Compare and comment on the timing results for the three methods you implemented.
Include a mesh plot of $x$, a table of CPU times, and your comments on the timing results.
\SolutionQuestionDifferenceEquationsPartC
\end{enumerate}

\noindent
Submit to the LEARN Dropbox: \texttt{Lap2D.m}, \texttt{GaussElim.m}, \texttt{Cholesky.m}, \texttt{BandGE.m}, \texttt{GETimes.m}.
\end{enumerate}

\end{document}