\documentclass[a4paper,12pt]{article}
\usepackage{graphicx} % Required for inserting images
\graphicspath{ {./figs/} }
\usepackage[margin=1.5cm]{caption}
\usepackage[T1]{fontenc}  % \ae\o\aa
\usepackage[utf8]{inputenc}
\usepackage{babel}          % ordelingsregler
\usepackage{amsmath, mathtools}
\usepackage{enumitem}
\setlist[itemize]{itemsep=0.3pt}
\usepackage{xcolor}
\usepackage{textgreek}
\usepackage{tabularx}
\usepackage{float}
\usepackage{hyperref}
\usepackage{subcaption} %needed for subfigures
\usepackage{placeins}
\usepackage{comment}
\usepackage{algorithm}
\usepackage{algpseudocode}

\usepackage{listings}
\usepackage{xcolor}

\definecolor{green}{RGB}{0,128,0}

\lstdefinestyle{mystyle}{
language=Python,
basicstyle=\small\ttfamily,
commentstyle=\color{green!40!black},
keywordstyle=\color{blue},
numberstyle=\tiny\color{gray},
numbers=left,
frame=single,
breaklines=true,
showstringspaces=false,
captionpos=b,
extendedchars=true,
inputencoding=utf8,
literate={"o}{{"o}}1 {"a}{{"a}}1 {"u}{{"u}}1
}

\lstset{style=mystyle}


\title{Project 2- Computational Physics}
\author{Vladislav Foss\\ Simen Husby Norrud\\ Sebastian Koranda}
\date{September 2025}
\begin{document}

\maketitle

\section*{Problem 1}
Our 2nd order differential function for this problem is given as:
\begin{align}
    \gamma \frac{d^2u(x)}{dx^2} = -Fu(x) \label{eq:org_func}
\end{align}
We will now be transforming this into a scaled equation with the following "limitation":
\begin{align}
    \hat{x} & \equiv x/L  \label{eq:hatxsub}\\
    \hat{x} & \in [0,1] \\
    \lambda & = F\frac{L^2}{\gamma}
\end{align}
We will firstly start off with a variable substitution in our equation with the relationship between $\hat{x}$ and $x$ given in Eq.(\ref{eq:hatxsub}).

From calculus we know that a variable substitution is given as $\frac{d}{dx} = \frac{d\hat{x}}{dx}\frac{d}{d\hat{x}}$. With this premise the variable change will give us:
\begin{align*}
    \frac{d}{dx} &= \frac{1}{L}\frac{d}{d\hat{x}} \\
    \text{meaning that:} \\
    \frac{d^2}{dx^2} &= \frac{d}{dx}(\frac{1}{L}\frac{d}{d\hat{x}^2}) = \frac{1}{L^2}\frac{d}{d\hat{x}}
\end{align*}
Now substituting into our Eq(\ref{eq:org_func}), we get:
\begin{align*}
    \gamma \frac{1}{L^2}\frac{d^2u(x)}{d\hat{x}^2} = -Fu(\hat{x}) \\
    \frac{d^2u(x)}{d\hat{x}^2} = -\frac{FL^2}{\gamma}u(\hat{x}) \\
    \frac{d^2u(x)}{d\hat{x}^2} = -\lambda u(\hat{x})
\end{align*}
Thereby we have shown how you can write the original function given in Eq(\ref{eq:org_func}) into $\frac{d^2u(x)}{d\hat{x}^2} = -\lambda u(\hat{x})$.


\section*{Problem 2}

We first computed the analytical eigenvalues and eigenvectors, and then verified our implementation by constructing the tridiagonal matrix and comparing its eigenvalues obtained with Armadillo’s \texttt{arma::eig\_sym} function to the analytical results. Table 1 compares the eigenvalues of the computed and analytical solutions. Table 2 presents the maximum differences between each of the Armadillo and analytic eigenvectors (allowing for sign flips), which are also on the order of $10^{-16}$, confirming that Armadillo's solution agrees with our analytical solution.

\begin{table}[H]
\centering
\begin{tabular}{c|c|c|c}
\hline
Mode $j$ & Armadillo $\lambda_j$ & Analytic $\lambda_j$ & $|\lambda_j^{\text{Arma}} - \lambda_j^{\text{Analytic}}|$ \\
\hline
1 & $9.7051 \times 10^{0}$  & $9.7051 \times 10^{0}$  & $3.5527 \times 10^{-15}$ \\
2 & $3.6898 \times 10^{1}$  & $3.6898 \times 10^{1}$  & $2.1316 \times 10^{-14}$ \\
3 & $7.6193 \times 10^{1}$  & $7.6193 \times 10^{1}$  & $1.4211 \times 10^{-14}$ \\
4 & $1.1981 \times 10^{2}$  & $1.1981 \times 10^{2}$  & $2.8422 \times 10^{-14}$ \\
5 & $1.5910 \times 10^{2}$  & $1.5910 \times 10^{2}$  & $0$ \\
6 & $1.8629 \times 10^{2}$  & $1.8629 \times 10^{2}$  & $5.6843 \times 10^{-14}$ \\
\hline
\end{tabular}
\caption{Difference of eigenvalues between Armadillo and analytic solution.}
\end{table}


\begin{table}[H]
\centering
\begin{tabular}{c|c}
\hline
Mode $j$ & Max difference \\
\hline
1 & $2.22045 \times 10^{-16}$ \\
2 & $3.33067 \times 10^{-16}$ \\
3 & $5.82867 \times 10^{-16}$ \\
4 & $3.88578 \times 10^{-16}$ \\
5 & $7.77156 \times 10^{-16}$ \\
6 & $5.27356 \times 10^{-16}$ \\
\hline
\end{tabular}
\caption{Maximum difference between Armadillo and analytic eigenvectors.}
\end{table}

\section*{Problem 3}

We implemented the function \texttt{double max\_offdiag\_symmetric(const arma::mat\& A, int\& k, int\& l)} which, given a symmetric matrix \(A\), loops over the upper triangular part (since for symmetric matrices the lower part is identical) and finds the element with the largest absolute value. The function also updates the indices \(k\) and \(l\) to point to the location of this element and returns its value.

We test the function with the matrix:
\[
A =
\begin{bmatrix}
1.0 & 0.0 & 0.0 & 0.5 \\
0.0 & 1.0 & -0.7 & 0.0 \\
0.0 & -0.7 & 1.0 & 0.0 \\
0.5 & 0.0 & 0.0 & 1.0
\end{bmatrix}
\]
And correctly returns $-0.7$, at row 2, column 3.


\section*{Problem 4}
We want to solve the eigenvalue problem $ A\vec{v} = \lambda \vec{v}$, where A is going to be a tri-diagonal matrix.

If we reorder the equation we get:
\begin{align}
    \mathbf{A}\vec{x} = \lambda \vec{x} \\
    (\mathbf{A} - \lambda \mathbf{I})x = 0
\end{align}

The equation det($\mathbf{A}-\lambda \mathbf{I}$) will then be a polynomial of degree N. The problem arises when N gets very large and thereby becomes very inefficient to compute. Therefore we introduce \textbf{Jacobi's rotation Method}.

The motivation for Jacobi's rotation method is that if $\mathbf{A}$ can be transformed into a diagonal form, 
its diagonal entries are the eigenvalues, and the corresponding eigenvectors are given by the columns of the transformation matrix. 
To achieve this we apply an orthogonal rotation matrix $\mathbf{R}$ such that
$$
    \mathbf{R}^T \mathbf{R} = \mathbf{I}, \qquad \mathbf{A}' = \mathbf{R}^T \mathbf{A} \mathbf{R}.
$$
Orthogonal transformations preserve eigenvalues and keep the matrix symmetric. 
At each step, $\mathbf{R}$ is chosen to eliminate the largest off-diagonal element in $\mathbf{A}$. 
By repeating this process, the matrix becomes increasingly diagonal, 
until all off-diagonal terms are below a given tolerance. 
The diagonal elements then approximate the eigenvalues, 
and the accumulated product of rotations $\mathbf{R}$ contains the eigenvectors.


In our code we have a tolerance of $tol = 1e^{-15}$. The function \texttt{max\_offdiag\_symmetric} scans through the upper triangular part of the symmetric matrix $\mathbf{A}$. 

It returns the indices $(k,l)$ of the largest off-diagonal element in absolute value, as well as its magnitude. 
This step identifies the element that contributes most to the "non-diagonality" of the matrix and is therefore the best candidate to eliminate in the next rotation.

The function \texttt{jacobi\_rotate} eliminates the chosen off-diagonal element $A_{kl}$. 
Given the indices $(k,l)$ from the previous step, it computes the rotation parameters
$$
    \tau = \frac{A_{ll} - A_{kk}}{2A_{kl}}, \qquad
    t = -\tau \pm \sqrt{1 + \tau ^2} , \qquad
    c_{\theta} = \frac{1}{\sqrt{1+t^2}}, \qquad
    s = c_{\theta} t_{\theta},
$$
and applies the orthogonal similarity transformation 
$$
    \mathbf{A}' = \mathbf{R}_{kl}^T \mathbf{A} \mathbf{R}_{kl}.
$$
This sets $A'_{kl}=A'_{lk}=0$, updates the diagonal entries $A_{kk}$ and $A_{ll}$, and mixes the remaining elements in rows/columns $k$ and $l$. 
At the same time, the eigenvector matrix $\mathbf{R}$ is updated with the same rotation so that its columns converge to the eigenvectors.

\section*{Problem 5}




\subsection*{Scaling behaviour}

For the tridiagonal matrix of  A, most elements will be zero and Jacobi's method only needs to eliminate the non-zero elements off the diagonal.
In this case, the number of required similarity transformations scales somewhat like $\mathcal{O}(N^2)$. 


If $\mathbf{A}$ is instead a dense symmetric matrix, then there are $\tfrac{N(N-1)}{2}$ off-diagonal elements that must be reduced. 
Each Jacobi rotation eliminates only one element at a time, while simultaneously altering an entire row and column. 
As a result, the number of required transformations scales more poorly, on the order of $\mathcal{O}(N^3)$, and since each rotation requires $\mathcal{O}(N)$ work, 
the total computational complexity of Jacobi's method for a dense matrix is $\mathcal{O}(N^4)$. 


For a dense symmetric matrix $\mathbf{A}$, the number of required Jacobi rotations still scales approximately like 
$$
\mathcal{O}(N^2).
$$
This is because one sweep involves applying a rotation for each of the 
$$
\frac{N(N-1)}{2} = \mathcal{O}(N^2)
$$
off-diagonal elements. In practice, convergence is usually reached in only a few sweeps, so the total number of similarity transformations grows quadratically with $N$, similar to the tridiagonal case.  

However, there is a crucial difference in terms of computational cost.  
Each Jacobi rotation in a dense matrix requires updating an entire row and column of $\mathbf{A}$, which is an 
$$
\mathcal{O}(N)
$$ 
operation. Thus, while the number of transformations is $\mathcal{O}(N^2)$, the \emph{total computational cost} is  
$$
\mathcal{O}(N^2) \times \mathcal{O}(N) = \mathcal{O}(N^3).
$$
This scaling matches that of more advanced algorithms (such as the QR method), but with a much larger constant factor, which makes Jacobi inefficient for large dense matrices.  

Table~\ref{tab:dense_sparse} and Figure~\ref{fig:njacobi} compares the number of Jacobi rotations required for convergence between the sparse tridiagonal case and a dense random symmetric matrix of the same size. The results illustrate that the number of rotations is very similar in both cases, since Jacobi is not sparsity-preserving: each rotation introduces additional nonzero elements, and the matrix quickly behaves like a dense one.  

\begin{table}[h!]
\centering
\begin{tabular}{c c c}
\hline
$N$ & Jacobi rotations (Sparse) & Jacobi rotations (Dense) \\
\hline
6   & 42    & 58    \\
10  & 168   & 179   \\
20  & 758   & 795   \\
40  & 3288  & 3265  \\
100 & 20778 & 20620 \\
200 & 83733 & 83025 \\
\hline
\end{tabular}
\caption{Number of Jacobi rotations required for convergence for sparse (tridiagonal) vs dense symmetric matrices.}
\label{tab:dense_sparse}
\end{table}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.8\linewidth]{iterations.pdf}
    \caption{Number of Jacobi rotations for sparse and dense matrix}
    \label{fig:njacobi}
\end{figure}

\section*{Problem 6}
Figure~\ref{fig:eigen10} shows the three lowest eigenvectors obtained using our Jacobi algorithm for a discretization with n $=10$ steps, together with the corresponding analytical eigenvectors (including the boundary points). We see that the numerical and analytical eigenvectors overlap almost perfectly, confirming the accuracy of our implementation for small n.

Figure~\ref{fig:eigen100} presents the same comparison for n$=100$ steps. Again, the numerical eigenvectors follow the analytical ones closely, and the higher resolution gives smoother curves. This demonstrates that our method converges to the analytical solutions as n increases.

\begin{figure}[H]
    \centering
    \includegraphics[width=0.8\linewidth]{eigvecs_10n.pdf}
    \caption{Eigenvectors corresponding to the three lowest eigenvalues for N = 10. Analytic and computed values.}
    \label{fig:eigen10}
\end{figure}
\begin{figure}[H]
    \centering
    \includegraphics[width=0.8\linewidth]{eigvecs_100n.pdf}
    \caption{Eigenvectors corresponding to the three lowest eigenvalues for N = 100. Analytic and computed values.}
    \label{fig:eigen100}
\end{figure}

\appendix
\section{Source code}
GitHub repository with code: \url{https://github.uio.no/sebasmk/fys4150.git}


\end{document}
