lars/FYS4150
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add project 1 with the corresponding source code and report.

Browse Source
This commit is contained in:
Lars Bogner
2025-08-31 14:03:13 +02:00
commit 49a519379c
18 changed files with 2421 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
projects/project1/include/citations.bib
+392
View File
@@ -0,0 +1,392 @@
@inreference{noauthor_tridiagonal_2025,
title = {Tridiagonal matrix algorithm},
rights = {Creative Commons Attribution-{ShareAlike} License},
url = {https://en.wikipedia.org/w/index.php?title=Tridiagonal_matrix_algorithm&oldid=1306920126},
abstract = {In numerical linear algebra, the tridiagonal matrix algorithm, also known as the Thomas algorithm (named after Llewellyn Thomas), is a simplified form of Gaussian elimination that can be used to solve tridiagonal systems of equations. A tridiagonal system for n unknowns may be written as
a
i
x
i
−
1
+
b
i
x
i
+
c
i
x
i
+
1
=
d
i
,
\{{\textbackslash}displaystyle a\_\{i\}x\_\{i-1\}+b\_\{i\}x\_\{i\}+c\_\{i\}x\_\{i+1\}=d\_\{i\},\}
where
a
1
=
0
\{{\textbackslash}displaystyle a\_\{1\}=0\}
and
c
n
=
0
\{{\textbackslash}displaystyle c\_\{n\}=0\}
.
[
b
1
c
1
0
a
2
b
2
c
2
a
3
b
3
⋱
⋱
⋱
c
n
−
1
0
a
n
b
n
]
[
x
1
x
2
x
3
⋮
x
n
]
=
[
d
1
d
2
d
3
⋮
d
n
]
.
\{{\textbackslash}displaystyle \{{\textbackslash}begin\{bmatrix\}b\_\{1\}\&c\_\{1\}\&\&\&0{\textbackslash}{\textbackslash}a\_\{2\}\&b\_\{2\}\&c\_\{2\}\&\&{\textbackslash}{\textbackslash}\&a\_\{3\}\&b\_\{3\}\&{\textbackslash}ddots \&{\textbackslash}{\textbackslash}\&\&{\textbackslash}ddots \&{\textbackslash}ddots \&c\_\{n-1\}{\textbackslash}{\textbackslash}0\&\&\&a\_\{n\}\&b\_\{n\}{\textbackslash}end\{bmatrix\}\}\{{\textbackslash}begin\{bmatrix\}x\_\{1\}{\textbackslash}{\textbackslash}x\_\{2\}{\textbackslash}{\textbackslash}x\_\{3\}{\textbackslash}{\textbackslash}{\textbackslash}vdots {\textbackslash}{\textbackslash}x\_\{n\}{\textbackslash}end\{bmatrix\}\}=\{{\textbackslash}begin\{bmatrix\}d\_\{1\}{\textbackslash}{\textbackslash}d\_\{2\}{\textbackslash}{\textbackslash}d\_\{3\}{\textbackslash}{\textbackslash}{\textbackslash}vdots {\textbackslash}{\textbackslash}d\_\{n\}{\textbackslash}end\{bmatrix\}\}.\}
For such systems, the solution can be obtained in
O
(
n
)
\{{\textbackslash}displaystyle O(n)\}
operations instead of
O
(
n
3
)
\{{\textbackslash}displaystyle O(n{\textasciicircum}\{3\})\}
required by Gaussian elimination. A first sweep eliminates the
a
i
\{{\textbackslash}displaystyle a\_\{i\}\}
's, and then an (abbreviated) backward substitution produces the solution. Examples of such matrices commonly arise from the discretization of 1D Poisson equation and natural cubic spline interpolation.
Thomas' algorithm is not stable in general, but is so in several special cases, such as when the matrix is diagonally dominant (either by rows or columns) or symmetric positive definite; for a more precise characterization of stability of Thomas' algorithm, see Higham Theorem 9.12. If stability is required in the general case, Gaussian elimination with partial pivoting ({GEPP}) is recommended instead.},
booktitle = {Wikipedia},
urldate = {2025-08-28},
date = {2025-08-20},
langid = {english},
note = {Page Version {ID}: 1306920126},
file = {Snapshot:/home/lars/Zotero/storage/LJCAC2BT/index.html:text/html},
}
projects/project1/main.pdf
BIN
View File
Binary file not shown.
projects/project1/main.tex
+334
View File
@@ -0,0 +1,334 @@
\documentclass[english,notitlepage]{revtex4-1} % defines the basic parameters of the document
%For preview: skriv i terminal: latexmk -pdf -pvc filnavn
\let\tablenum\relax % to avoid clash with siunitx
% if you want a single-column, remove reprint
% allows special characters (including æøå)
\usepackage[utf8]{inputenc}
%\usepackage[english]{babel}
%% note that you may need to download some of these packages manually, it depends on your setup.
%% I recommend downloading TeXMaker, because it includes a large library of the most common packages.
\usepackage{amssymb} % mathematical symbols (physics imports amsmath)
\usepackage{amsmath}
\usepackage{graphicx} % include graphics such as plots
\usepackage{xcolor} % set colors
\usepackage{hyperref} % automagic cross-referencing (this is GODLIKE)
\usepackage{listings} % display code
\usepackage{subfigure} % imports a lot of cool and useful figure commands
\usepackage{float}
%\usepackage[section]{placeins}
\usepackage{algorithm}
\usepackage[noend]{algpseudocode}
\usepackage{subfigure}
\usepackage{tikz}
\usepackage{siunitx}
\usetikzlibrary{quantikz}
% defines the color of hyperref objects
% Blending two colors: blue!80!black = 80% blue and 20% black
\hypersetup{ % this is just my personal choice, feel free to change things
colorlinks,
linkcolor={red!50!black},
citecolor={blue!50!black},
urlcolor={blue!80!black}}
%% Defines the style of the programming listing
%% This is actually my personal template, go ahead and change stuff if you want
%% USEFUL LINKS:
%%
%% UiO LaTeX guides: https://www.mn.uio.no/ifi/tjenester/it/hjelp/latex/
%% mathematics: https://en.wikibooks.org/wiki/LaTeX/Mathematics
%% PHYSICS ! https://mirror.hmc.edu/ctan/macros/latex/contrib/physics/physics.pdf
%% the basics of Tikz: https://en.wikibooks.org/wiki/LaTeX/PGF/Tikz
%% all the colors!: https://en.wikibooks.org/wiki/LaTeX/Colors
%% how to draw tables: https://en.wikibooks.org/wiki/LaTeX/Tables
%% code listing styles: https://en.wikibooks.org/wiki/LaTeX/Source_Code_Listings
%% \includegraphics https://en.wikibooks.org/wiki/LaTeX/Importing_Graphics
%% learn more about figures https://en.wikibooks.org/wiki/LaTeX/Floats,_Figures_and_Captions
%% automagic bibliography: https://en.wikibooks.org/wiki/LaTeX/Bibliography_Management (this one is kinda difficult the first time)
%% REVTeX Guide: http://www.physics.csbsju.edu/370/papers/Journal_Style_Manuals/auguide4-1.pdf
%%
%% (this document is of class "revtex4-1", the REVTeX Guide explains how the class works)
%% CREATING THE .pdf FILE USING LINUX IN THE TERMINAL
%%
%% [terminal]$ pdflatex template.tex
%%
%% Run the command twice, always.
%% If you want to use \footnote, you need to run these commands (IN THIS SPECIFIC ORDER)
%%
%% [terminal]$ pdflatex template.tex
%% [terminal]$ bibtex template
%% [terminal]$ pdflatex template.tex
%% [terminal]$ pdflatex template.tex
%%
%% Don't ask me why, I don't know.
\begin{document}
\title{Solutions to Project 1} % self-explanatory
\author{Lars Bogner} % self-explanatory
\date{\today} % self-explanatory
\noaffiliation % ignore this, but keep it.
\maketitle
\textit{List a link to your github repository here!}
\section*{Problem 1}
To prove that the given function is a solution to the one-dimensional Poisson equation, we are calculating the second derivative of $u(x)$ as
\begin{equation}
\frac{\partial u}{\partial x} = (1-e^{-10}) + 10 e^{-10x}
\end{equation}
and therefore
\begin{equation}
\frac{\partial^2 u}{\partial x^2} = -100 e^{-10x}
\end{equation}
Therefore the main equation $-\frac{\partial^2 u}{\partial x^2} = f(x)$ is satisfied with $f(x) = 100 e^{-10x}$. Now we have to check the remaining boundary conditions. There is no singularity within the $x$ range: $x \in [0, 1]$, therefore the Poisson equation is well-defined in this interval. Lastly we check the boundary conditions:
\begin{equation}
\begin{aligned}
u(0) &= 1 - 1 = 0\\
u(1) &= 1 - 1 + e^{-10} - e^{-10} = 0
\end{aligned}
\end{equation}
\section*{Problem 2}
The code for this problem can be found in the files \texttt{src/project1/poisson.cpp} and \texttt{src/\-project1/\-python/\-poisson\_plotter.py}. The resulting plot is displayed in \autoref{fig:poisson_plot}.
\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{include/poisson_solution.pdf}
\caption{Numerical values of the solution to the Poisson equation}
\label{fig:poisson_plot}
\end{figure}
\section*{Problem 3}
To discretize the Poisson equation, we start with discretizing derivatives. For this to work we will Taylor expand the function $u(x)$ around a point $x_i$:
\begin{equation} \label{eq:taylor}
u(x_j) = u(x_i) + (x_j - x_i) u'(x_i) + \frac{(x_j - x_i)^2}{2} u''(x_i) + \ldots
\end{equation}
If our $\{x_j\}$ are equally spaced, i.e. $x_{j+1} = x_j + \Delta, \quad \forall j \in \mathbb{Z}$, we can write the Taylor expansion in terms of the spacing $\Delta$:
\begin{equation} \label{eq:discrete_taylor}
u(x_j) = u(x_i) + (j - i) \Delta u'(x_i) + \frac{(j - i)^2}{2} \Delta^2 u''(x_i) + \ldots
\end{equation}
To find a discrete expression for the second derivative, we use the discrete version of derivatives:
\begin{equation}
\frac{\partial g(x)}{\partial x} = \lim_{\Delta \to 0} \frac{g(x+\Delta / 2) - g(x - \Delta / 2)}{\Delta}
\end{equation}
and thus follows
\begin{equation}
\frac{\partial}{\partial x} \frac{\partial g(x)}{\partial x} = \lim_{\Delta \to 0} \frac{1}{\Delta^2} \left( g(x + \Delta) - g(x) - g(x) + g(x - \Delta) \right).
\end{equation}
If we take this expression and apply it to $u(x_j)$, we get
\begin{equation}
\frac{\partial^2 u(x_j)}{\partial x_j^2} = \lim_{\Delta \to 0}\frac{u(x_{j+1}) - 2u(x_j) + u(x_{j-1})}{\Delta^2}
\end{equation}
The rest of the problem can be discretized by mapping the continous domain of $x \in [0, 1]$ to a discrete spectrum of points, i.e. $f(x) \to f_j = f(x_j)$ and $u(x) \to u_j = u(x_j)$. If we discard the limit of $\Delta \to 0$ we move to an approximation of $u_j$ which we call $v_j$. The problem can then be written as
\begin{equation}
- \frac{v_{j+1} - 2v_j + v_{j-1}}{\Delta^2} = f_j \quad \forall j \in \{1, \ldots, n-1\}
\end{equation}
\section*{Problem 4}
When slightly rearranging the discretized Poisson equation from Problem 3, we get
\begin{equation}
v_{j+1} - 2v_j + v_{j-1} = -\Delta^2 f_j \quad \forall j \in \{1, \ldots, n-1\}
\end{equation}
which can be directly interpreted as a matrix equation of the form $A \vec v = \vec g$. We need to pay special attention to the border regions of the matrix equation, as we can't compute $v_0$ and $v_n$. If we absorb those into $\vec g$, we get a correct matrix equation system. Therefore we define $g_j = -\Delta^2 f_j ,\,\forall j \in \{2, \ldots, n-2\}$ and $g_1 = -\Delta^2 f_1 - v_0$ and $g_{n-1} = -\Delta^2 f_{n-1} - v_n$. The matrix $A$ is a tridiagonal matrix with the following structure:
\begin{itemize}
\item The main diagonal consists of $-2$.
\item The first superdiagonal consists of $1$.
\item The first subdiagonal consists of $1$.
\end{itemize}
\section*{Problem 5}
\subsection*{Subproblem (a)}
As our discrete solution always depends on one value before and one value after the current value, our matrix has to be smaller in either direction by one row and column. If we put this into an equation we can say:
\begin{equation}
n = m - 2
\end{equation}
\subsection*{Subproblem (b)}
We can find the inner part of $\vec v^\ast$, i.e. $v_i = v_{i+1}^\ast \forall i \in \{1, \ldots, n\}$.
\section*{Problem 6}
\subsection*{Subproblem (a)}
To solve a general tridiagonal matrix equation we can use the algorithm described in \cite{noauthor_tridiagonal_2025}. The algorithm is structured as described in Algorithm \ref{algo:general}.
\begin{figure}
% NOTE: We only need \begin{figure} ... \end{figure} here because of a compatability issue between the 'revtex4-1' document class and the 'algorithm' environment.
\begin{algorithm}[H]
\caption{General Algorithm (Thomas Algorithm)\cite{noauthor_tridiagonal_2025} for solving $A \vec v = \vec g$ where $A$ is a tridiagonal matrix}
\label{algo:general}
\begin{algorithmic}
\Function{General Algorithm}{$\vec a, \vec b, \vec c, \vec g$}
\State $\vec {c'} \gets \vec 0$ \Comment{Initialize temporary vectors}
\State $\vec {g'} \gets \vec 0$
\State $c'_1 \gets c_1 / b_1$
\For{$i = 2, \ldots, n-1$} \Comment{Recursively compute $c'$}
\State $c'_i \gets c_i / (b_i - a_i c'_{i-1})$
\EndFor
\State $g'_1 \gets g_1 / b_1$
\For{$i = 2, \ldots, n$} \Comment{Recursively compute $g'$}
\State $g'_i \gets (g_i - a_i g'_{i-1}) / (b_i - a_i c'_{i-1})$
\EndFor
\State $\vec v \gets \vec 0$ \Comment{Initialize solution vector}
\State $v_n \gets g'_n$
\For{$i = n-1, n-2, \ldots, 1$} \Comment{Back substitution}
\State $v_i \gets g'_i - c'_i v_{i+1}$
\EndFor
\Return $\vec v$
\EndFunction
\end{algorithmic}
\end{algorithm}
\end{figure}
\subsection*{Subproblem (b)}
The floating point operations (FLOPs) can be counted as follows:
\begin{itemize}
\item Calculate first element of $c'$: 1 FLOP
\item Calculate rest of $c'$: 3 FLOPs per element ($3 (n-2)$ FLOPs)
\item Calculate first element of $g'$: 1 FLOP
\item Calculate rest of $g'$: 5 FLOPs per element ($5 (n-1)$ FLOPs)
\item Back substitution: 2 FLOPs per element ($2 (n-1)$ FLOPs)
\end{itemize}
In total this gives us
\begin{equation}
N = 1 + 3(n-2) + 1 + 5(n-1) + 2(n-1) = 10 \cdot n - 11 \sim \mathcal{O}(n)
\end{equation}
floating point operations (FLOPs).
\section*{Problem 7}
\subsection*{Subproblem (a)}
The implementation of the specialized algorithm can be found in the file \texttt{src/\-project1/\-src/\-solvers.cpp} with usage in \texttt{src/\-project1/\-poisson\-\_solver.cpp}.
\subsection*{Subproblem (b)}
Comparing the solutions for different values of $N$ we can clearly see the increase in accuracy as $N$ is increased. The numerical solutions converge towards the analytical solution. The results are summarized in \autoref{fig:big_ref_plot}. To assess the quality of the results, the solutions are plotted in comparison to the numerical solution in the upper plot, while the lower plot shows the relative error. The relative error is defined as
\begin{equation} \label{eq:relerror}
\frac{\delta v_i}{u_i} = \frac{v_i - u_i}{u_i}
\end{equation}
\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{include/big_ref_plot.pdf}
\caption{Numerical values of the solution to the Poisson equation for different values of $N$}
\label{fig:big_ref_plot}
\end{figure}
\section*{Problem 8}
\subsection*{Subproblem (a)}
Plotting the absolute error for different values of $N$ reveals the said convergence behavior. As $N$ increases, the absolute error decreases, indicating that the numerical solution is approaching the analytical solution. The results are summarized in \autoref{fig:abs_error_plot}.
\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{include/abs_error_plot.pdf}
\caption{Absolute error of the numerical solution for different values of $N$}
\label{fig:abs_error_plot}
\end{figure}
\subsection*{Subproblem (b)}
Plotting the relative error for different values of $N$ reveals the increase of relative error as $u$ approaches 0 in the vicinity of the boundaries. This behavior is prevalent in all solutions as the limit of small divisors is approached. This behavior is quite logical, but it also shows, that not all uncertainties can be correctly minimized with an increase in resolution. The results are shown in \autoref{fig:rel_error_plot}. To overcome the division by zero issue for the relative error, only the values of $u > \num{1e-10}$ are considered.
\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{include/rel_error_plot.pdf}
\caption{Relative error of the numerical solution for different values of $N$}
\label{fig:rel_error_plot}
\end{figure}
\subsection*{Subproblem (c)}
The maximum relative errors are calculated directly using C++ as the output size scales linearly with $N$ and gets overwhelmingly large for bigger $N$. The results are output using the \texttt{std::cout} command and retrieved using GNU/Linux command line utilities. The code for calculating the maximum relative error $\max_i \frac{\delta v_i}{u_i}$ (see \eqref{eq:relerror}) is included in the C++ implementation. The results are summarized in \autoref{tab:rel_errors}. We see a steady decrease of the relative error up to $N = \num{100000}$, after which the relative error increases again. This will most likely be due to the artifacts seen in \autoref{fig:rel_error_plot} starting at $N = \num{10000}$. The source of these artifacts is not entirely clear.
\begin{table}[H]
\centering
\begin{tabular}{
S[table-format=7]
S[scientific-notation=true,round-mode=places,round-precision=3]
}
\hline \\
{$N$} \relax & {Max Relative Errors} \\ \hline
10 \relax & 0.0793264 \\
100 \relax & 0.000832917 \\
1000 \relax & 8.33329e-06 \\
10000 \relax & 8.33301e-08 \\
100000 \relax & 1.4356e-09 \\
1000000 \relax & 8.40412e-07 \\
10000000 \relax & 2.9838e-06 \\
\hline
\end{tabular}
\caption{Summary of maximum relative errors for different values of $N$}
\label{tab:rel_errors}
\end{table}
\section*{Problem 9}
\subsection*{Subproblem (a)}
To specialize the general algorithm to our specific problem, we take a look at how our general algorithm modifies the input vectors during forward and backward substitution. As a first step, we can calculate the values $c_i'$ from Algorithm \ref{algo:general} explicitly as $c_i' = \frac{-i}{i+1}$. This explicit formulation also allows for a more straightforward definition of $g_i' \equiv \tilde{g}_i$. Combining these two results we can rewrite the algorithm as described in Algorithm \ref{algo:optimal}.
\begin{figure}
\begin{algorithm}[H]
\caption{Special Algorithm for computing $\vec v$ from $\vec g$ using constant tridiagonal elements}
\label{algo:optimal}
\begin{algorithmic}
\Function{Optimal Algorithm}{$\vec g$}
\State $\vec{v} \gets \vec{0}$ \Comment{Initialize solution vector}
\State $\vec{\tilde{g}} \gets \vec{0}$ \Comment{Initialize modified input vector}
\State $\tilde{g}_0 \gets g_0 / 2$
\For{$i = 1, \ldots, n-1$} \Comment{Forward substitution}
\State $\tilde{g}_i \gets (g_i + \tilde{g}_{i-1}) \cdot \frac{i+1}{i+2}$
\EndFor
\State $v_{n-1} \gets \tilde{g}_{n-1}$
\For{$i = n-2, n-3, \ldots, 0$} \Comment{Backward substitution}
\State $v_i \gets \tilde{g}_i + \frac{i+1}{i+2} \cdot v_{i+1}$
\EndFor
\Return $\vec{v}$
\EndFunction
\end{algorithmic}
\end{algorithm}
\end{figure}
\subsection*{Subproblem (b)}
The floating point operations (FLOPs) can be counted as follows:
\begin{itemize}
\item Calculate first element of $\tilde{g}$: 1 FLOP
\item Calculate remaining elements of $\tilde{g}$: $3(n-1)$ FLOPs
\item Calculate last element of $v$: 0 FLOPs
\item Calculate remaining elements of $v$: $3(n-1)$ FLOPs
\end{itemize}
If we dismiss integer additions and subtractions, the total number of FLOPs can be simplified to:
\begin{equation}
N =1 + 3(n-1) + 0 + 3(n-1) = 6\cdot n - 5 \sim \mathcal{O}(n)
\end{equation}
\subsection*{Subproblem (c)}
The implementation of the specialized algorithm can be found in the file \texttt{src/\-project1/\-src/\-solvers.cpp} with usage in \texttt{src/\-project1/\-poisson\_timing.cpp}. Secondly a variant with further optimizations on the memory access patterns is implemented in the same location. From now on we will reference this version as the \textit{Optimized Algorithm}.
\section*{Problem 10}
To evaluate the performance of the different algorithms, the algorithms were executed 100 times for varying problem sizes $N$. The execution times were recorded and analyzed to compare the efficiency of the algorithms. To compare the effect of different $N$ on the execution time, the execution times were normalized by dividing them by $N$. The results are summarized in \autoref{fig:timing_plot}. We can observe, that the specialized algorithm outperforms the general algorithm consistently across all problem sizes by more than a factor of 2. The optimized algorithm shows a further improvement over the specialized algorithm with very consistent execution times across all problem sizes.
\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{include/timing_results_100_iterations.pdf}
\caption{Normalized execution times for different algorithms and problem sizes}
\label{fig:timing_plot}
\end{figure}
\bibliography{include/citations}
\end{document}
projects/project1/mainNotes.bib
+2
View File
@@ -0,0 +1,2 @@
@CONTROL{REVTEX41Control}
@CONTROL{apsrev41Control,author="08",editor="1",pages="0",title="",year="1"}