The \textit{LU decomposition} is not a direct method which solves a linear system, but it
allows to simplify the resolution by decomposing the $A$ matrix into a lower-triangular matrix $L$ and
an upper-triangular matrix $U$. The simplification is due to the major facility to invert the two matrices precedently
presented.
Once $A$ is decomposed, the process is straigh-forward:
\begin{align}
A \cdot \vec{x} = L \cdot U \cdot \vec{x} &= \vec{b} \nonumber \\
L \cdot \vec{y} &= \vec{b} \label{L_solve} \\
U \cdot \vec{x} &= \vec{y} \label{U_solve}
\end{align}
Both equations (\ref{L_solve}) and (\ref{U_solve}) can be solved sequentially using the
\textit{Gauss elimination} method.
\subsection{Diagonalization: introduction}
In case $A$ is a symmetric matrix, the spectral theorem \cite{spectral_thm} states that
such a matrix is equivalent (definition of equivalence here: \cite{matrix_equiv}) to a diagonal matrix $D$,
where the transition matrix $P$ is unitary ($P^{-1} = \bar{P}^T$), then:
\begin{equation} \label{diag_solve}
A = P \cdot D \cdot \bar{P}^T \implies \vec{x} = P \cdot D^{-1} \cdot \bar{P}^T \cdot \vec{b}
\end{equation}
Generally diagonalization is not used to solve general systems of linear equations,
but it's convenient when the problem is related to find the eigen-base related to the eigen-values.
%\newpage
\begin{homeworkProblem}
% When problems are long, it may be desirable to put a \newpage or a
% \clearpage before each homeworkProblem environment
\begin{homeworkSection}{(1) LU decomposition implementation}
This algorithm separes the input matrix $A$ into a lower triangular $L$ and an upper triangular $U$, garanteeing
that $A = L \cdot U$. Neverthless, not all the invertible square matrices are purely LU decomposable, then it may happen that the output can result ill formed.
The code (\ref{lu_decomposition}) shows at line 23 that
a division by the diagonal values is performed, causing eventually a singularity.
A possible work-around is to apply the partial pivoting technique in order to swap the problematic lines.
In listing (\ref{lu_decomposition}) is shown a full implementation with partial pivoting.
\matlabscript{lu_decomposition}{\textit{LU decomposition} implementation with partial pivoting}
\subsubsection{Partial pivoting}
The \textit{LU decomposition} algorithm (presented below in exercise 1.1) can easily run into
singularities, especially when $A$ presents zeros as diagonal terms.
In order to avoid divergent results, it would better select the rows of which element is not zero
in the requested columns and swap them with the current one.
More precisely, at the $k$-th step, select the $r$-th row such that $A_{rk} = \max\limits_{k \le i \le N} \abs*{A_{ik}}$,
then swap rows at the position $k$ and $r$.
If the pivoting is applied the resulting \textit{LU decomposition} won't be anymore like it was defined in the
previous section, but a correction to equation (\ref{L_solve}) must be applied:
\begin{equation} \label{pivot_correction}
P \cdot A = L \cdot U \implies L \cdot \vec{y} = P \cdot \vec{b}
\end{equation}
where $P$ is the orthogonal matrix that accumulated all row switching applications.
The rest of the solving method remains unchanged.
\subsubsection{(1.1) Solving a linear system}
A linear system can be solved applying the LU decomposition and then a gauss elimination process, as
shown in the equations (\ref{pivot_correction}) and (\ref{U_solve}).
For example, the system in equation (\ref{first_system}) is determined and can be solved using
the \texttt{solve.m} script. Addictionally the \texttt{test\_solve.m} script compares with the matlab
\texttt{x = A \symbol{92} b} verifying that the solution $\vec{x}$ is given correctly by the \texttt{solve.m} script.
\begin{equation} \label{first_system}
\begin{cases}
2 x_1 + x_2 - x_3 + 5 x_4 &= 13 \\
x_1 + 2 x_2 + 3 x_3 - x_4 &= 37 \\
x_1 + x_3 + 6 x_4 &= 30 \\
x_1 + 3 x_2 - x_3 + 5 x_4 &= 19 \\
\end{cases}
\quad \implies A =
\begin{pmatrix}
2 & 1 & -1 & 5 \\
1 & 2 & 3 & -1 \\
1 & 0 & 1 & 6 \\
1 & 3 & -1 & 5 \\
\end{pmatrix}, \quad \vec{b} =
\begin{pmatrix}
13 \\
37 \\
30 \\
19 \\
\end{pmatrix}
\quad \implies \vec{x} = A^{-1} \cdot \vec{b} =
\begin{pmatrix}
2 \\
4 \\
10 \\
3 \\
\end{pmatrix}
\end{equation}
\subsubsection{(1.2) Comparing with \texttt{matlab lu} built-in function}
\noindent
\begin{minipage}{\textwidth}
\begin{minipage}{0.45\textwidth}
The graph in figure (\ref{lu_compare}) underlines a comparison between
the previously implemented method and the built-in matlab \texttt{lu} function
within a range of sizes between $1$ and $100$.
The evident result is that they are the almost equal and
the matrices $U$ and $P$ are probably determined exactly in the same way
because of the completely null difference. However, there is a small difference, in the order
of $10^{-13}$ in the determination of $L$ and it seems to depend on the input matrix size and
the specific composition.
\end{minipage}
\hspace{0.8cm}
\begin{minipage}{0.57\textwidth}
\resizebox{\textwidth}{!}{
\input{graphs/lu_size.tex}
}
\captionof{figure}{Comparison of the output $L$, $U$ and $P$ between code (\ref{lu_decomposition}) and \texttt{matlab lu} function}
\label{lu_compare}
\end{minipage}
\end{minipage}
\subsubsection{(1.3) Decomposition of a matrix}
The example taken in equation (\ref{lu_example})
is a problematic case where a pure LU decomposition
doesn't exist. A necessary and sufficient condition to the existance of
a pure LU decomposition is that the matrix must be gauss reductible
without any row exchange (ref. \cite{exists_LU}), that's why if such a decomposition
exists, then pivoting matrix $P$ is the identity matrix.
So, the form $P \cdot A = L \cdot U$ is
obtainable using the pivoting described in the previous section.
\begin{equation} \label{lu_example}
A =
\begin{pmatrix}
1 & 2 & 3 \\
2 & 4 & 9 \\
4 & -3 & 1 \\
\end{pmatrix}
\implies
L =
\begin{pmatrix}
1 & 0 & 0 \\
0.5 & 1 & 0 \\
0.25 & 0.5 & 1 \\
\end{pmatrix}, \quad
U =
\begin{pmatrix}
4 & -3 & 1 \\
0 & 5.5 & 8.5 \\
0 & 0 & -1.5 \\
\end{pmatrix} \quad
P =
\begin{pmatrix}
0 & 0 & 1 \\
0 & 1 & 0 \\
1 & 0 & 0 \\
\end{pmatrix} \quad
\end{equation}
\end{homeworkSection}
\begin{homeworkSection}{(2) Puzzle board game}
The puzzle game problem is described by using the following formalization:
So, the output power of the second helicopter $P_2 = 3^{-\frac{7}{2}} \cdot P_1 \approx 0.021 P_1$.
\end{homeworkSection}
\end{homeworkProblem}
\section{The eigenvalue problem and diagonalization}
Let $\hat{A}$ be an operator defined over an hilbert space $\mathcal{H}$.
By solving an eigenvalue problem is meant to find all vectors (or functions) $x \in \mathcal{H}$
such that there exists a real (or complex) value $\lambda$ that satisfies the following condition:
\begin{equation} \label{eig_def}
\hat{A} \cdot x = \lambda \cdot x, \lambda \in \mathbb{K}
\end{equation}
In the case of this report, the interest is to computationally solve the eigenvalue problem for finite rank
operators, which can be expressed as square matrices.
So, let $N$ be rank of a square matrix $A$ and $\vec{v} \in \mathbb{K}^N$, then the equation (\ref{eig_def}) is equivalent to:
\begin{equation} \label{eig_vec}
A \cdot \vec{v} = \lambda \cdot \vec{v}, \lambda \in \mathbb{K}
\end{equation}
\subsection{Power method} \label{power_section}
The power method bases its functioning on the iterative application of a specific operation $T$.
The principle is that every iteration step tends to minimise of the distance between the old evaluated eigen value $\lambda_{k-1}$ and the current $\lambda_k$.
Given the unitary vector $\vec{v_k} \in \mathbb{K}^N | \|\vec{v_k}\| = 1$ at the iteration step $k$, the corresponding diagonal value, relative to a square matrix $A$, is given by the hermitian scalar product (see \cite{inner_product} for the notation):
The operation $T$ mentioned above variates depending on the specific method, which of there are three:
\begin{itemize}
\item \textit{Power} method: $T = A$.
\item \textit{Inverse power} method: $T = (A - 1 \cdot \tau)^{-1}$, $\tau \in \mathbb{K}$ is a fixed eigenvalue target.
\item \textit{Rayleigh quotient} method: $T = (A - 1 \cdot \lambda_{k-1})^{-1}$, $\lambda_{k-1} \in \mathbb{K}$ is the old evaluated eigenvalue, as defined in equation (\ref{herm_eig}).
\end{itemize}
Notice that for the \textit{Tayleigh quotient} method the application is adapting during each
iteration, garanteeing a faster convergence.
\subsection{Jacobi method}
Let $A$ be a real symmetric matrix, by the spectral theorem \cite{spectral_thm} such a matrix is diagonalizable in a form $A = PDP^T$,
where $P$ is an orthogonal matrix and $D$ a real diagonal matrix.
The idea at the base is that any orthogonal matrix $P$ can be decomposed into a series of \textit{axis aligned} rotation matrices $\{J^{(pq)}(\theta)\}, p,q \in \mathbb{N*}, p < q$:
an \textit{axis aligned} rotation matrix $J^{(pq)}$ is defined as:
$A$ in equation (\ref{secnd_deriv}) is the tridiagonal matrix containing $-1$ in the lower and the upper diagonal and $2$ on the diagonal,
all the terms divided by the squared step $\Delta x^2$. This matrix is symmetric, then by the spectral theorem \cite{spectral_thm}, it's eigenvalues
are real and for each pair of eigenvalues, the associated eigenspaces are orthogonal each other. Furthermore, the matrix is positive definite,
thus all the eigenvalues are strictly positive.
Notice that this setup automatically meets the boundary conditions, because the first and the last line of $A$ already discards the boundary terms (setting them implicitly to zero).
% TODO, units of measure!!!!
\subsubsection{(2.3) + (2.4), Find the first eigenvalues}
The first thing to notice is that the matrix $A$ is not degenerate (by its positivity), then
each eigenspace's dimension is $1$. This clearly means that the first four eigenvalues are different each other.
Theoretically, by the expression in equation (\ref{lambda_n}), the eigenvalues for $L = 1$ are