Change front ticks to back ticks.

Janitaws/plot and latex

.gitignore

@@ -36,3 +36,12 @@
 *.log
 *.out
 *.bib
+*.synctex.gz
+*.bbl 
+
+# C++ specifics
+src/*
+!src/Makefile
+!src/*.cpp
+!src/*.hpp
+!src/*.py

README.md

@@ -3,3 +3,36 @@
 ## Practical information
 
 - [Project](https://anderkve.github.io/FYS3150/book/projects/project1.html)
+
+## How to compile C++ code
+
+Make sure that you are inside the **src** directory before compiling the code.
+
+Now you can execute the command shown under to compile:
+
+```
+    make
+```
+
+This will create object files and link them together into 2 executable files.
+These files are called **main** and **analyticPlot**. 
+To run them, you can simply use the commands below:
+
+```
+    ./main
+```
+
+```
+    ./analyticPlot
+```
+
+## How to generate plots
+
+For generating the plots, there are 4 Python scripts. 
+You can run each one of them by using this command:
+
+```
+    python <PythonFile>
+```
+
+The plots will be saved inside **latex/images**.

latex/assignment_1.tex

@@ -1,7 +1,9 @@
 \documentclass[english,notitlepage]{revtex4-1}  % defines the basic parameters of the document
 %For preview: skriv i terminal: latexmk -pdf -pvc filnavn
 
-
+% Silence warning of revtex4-1
+\usepackage{silence}
+\WarningFilter{revtex4-1}{Repair the float}
 
 % if you want a single-column, remove reprint
 
@@ -13,7 +15,7 @@
 %% I recommend downloading TeXMaker, because it includes a large library of the most common packages.
 
 \usepackage{physics,amssymb}  % mathematical symbols (physics imports amsmath)
-\include{amsmath}
+\usepackage{amsmath}
 \usepackage{graphicx}         % include graphics such as plots
 \usepackage{xcolor}           % set colors
 \usepackage{hyperref}         % automagic cross-referencing (this is GODLIKE)
@@ -72,10 +74,13 @@
 %%
 %% Don't ask me why, I don't know.
 
+% custom stuff
+\graphicspath{{./images/}}
+
 \begin{document}
 
 \title{Project 1}      % self-explanatory
-\author{Cory Balaton \& Janita Willumsen}          % self-explanatory
+\author{Cory Alexander Balaton \& Janita Ovidie Sandtrøen Willumsen}          % self-explanatory
 \date{\today}                             % self-explanatory
 \noaffiliation                            % ignore this, but keep it.
 
@@ -84,84 +89,24 @@
     
 \textit{https://github.uio.no/FYS3150-G2-2023/Project-1}
     
-\section*{Problem 1}
+\input{problems/problem1}
 
-% Do the double integral
-\begin{align*}
-    u(x) &= \int \int \frac{d^2 u}{dx^2} dx^2\\
-         &= \int \int -100 e^{-10x} dx^2 \\
-         &= \int \frac{-100 e^{-10x}}{-10} + c_1 dx \\
-         &= \int 10 e^{-10x} + c_1 dx \\
-         &= \frac{10 e^{-10x}}{-10} + c_1 x + c_2 \\
-         &= -e^{-10x} + c_1 x + c_2
-\end{align*}
+\input{problems/problem2}
 
-Using the boundary conditions, we can find $c_1$ and $c_2$ as shown below:
+\input{problems/problem3}
 
-\begin{align*}
-    u(0) &= 0 \\
-    -e^{-10 \cdot 0} + c_1 \cdot 0 + c_2 &= 0 \\
-    -1 + c_2 &= 0 \\
-    c_2 &= 1
-\end{align*}
+\input{problems/problem4}
 
-\begin{align*}
-    u(1) &= 0 \\
-    -e^{-10 \cdot 1} + c_1 \cdot 1 + c_2 &= 0 \\
-    -e^{-10} + c_1 + c_2 &= 0 \\
-    c_1 &= e^{-10} - c_2\\
-    c_1 &= e^{-10} - 1\\
-\end{align*}
+\input{problems/problem5}
 
-Using the values that we found for $c_1$ and $c_2$, we get
+\input{problems/problem6}
 
-\begin{align*}
-    u(x) &= -e^{-10x} + (e^{-10} - 1) x + 1 \\
-         &= 1 - (1 - e^{-10}) - e^{-10x}
-\end{align*}
+\input{problems/problem7}
 
-\section*{Problem 2}
+\input{problems/problem8}
 
-% Write which .cpp/.hpp/.py (using a link?) files are relevant for this and show the plot generated.
+\input{problems/problem9}
 
-\section*{Problem 3}
-
-% Show how it's derived and where we found the derivation.
-
-\section*{Problem 4}
-
-% Show that each iteration of the discretized version naturally creates a matrix equation.
-
-\section*{Problem 5}
-
-\subsection*{a)}
-
-\subsection*{b)}
-
-\section*{Problem 6}
-
-\subsection*{a)}
-
-% Use Gaussian elimination, and then use backwards substitution to solve the equation
-
-\subsection*{b)}
-
-% Figure it out
-
-\section*{Problem 7}
-
-% Link to relevant files on gh and possibly add some comments
-
-\section*{Problem 8}
-
-%link to relvant files and show plots
-
-\section*{Problem 9}
-
-% Show the algorithm, then calculate FLOPs, then link to relevant files
-
-\section*{Problem 10}
-
-% Time and show result, and link to relevant files
+\input{problems/problem10}
 
 \end{document}

latex/problems/problem1.tex

@@ -0,0 +1,44 @@
+\section*{Problem 1}
+
+First, we rearrange the equation.
+
+\begin{align*}
+    - \frac{d^2u}{dx^2} &= 100 e^{-10x} \\
+    \frac{d^2u}{dx^2} &= -100 e^{-10x} \\
+\end{align*}
+
+Now we find $u(x)$.
+
+% Do the double integral
+\begin{align*}
+    u(x) &= \int \int \frac{d^2 u}{dx^2} dx^2 \\
+         &= \int \int -100 e^{-10x} dx^2 \\
+         &= \int \frac{-100 e^{-10x}}{-10} + c_1 dx \\
+         &= \int 10 e^{-10x} + c_1 dx \\
+         &= \frac{10 e^{-10x}}{-10} + c_1 x + c_2 \\
+         &= -e^{-10x} + c_1 x + c_2
+\end{align*}
+
+Using the boundary conditions, we can find $c_1$ and $c_2$
+
+\begin{align*}
+    u(0) &= 0 \\
+    -e^{-10 \cdot 0} + c_1 \cdot 0 + c_2 &= 0 \\
+    -1 + c_2 &= 0 \\
+    c_2 &= 1
+\end{align*}
+
+\begin{align*}
+    u(1) &= 0 \\
+    -e^{-10 \cdot 1} + c_1 \cdot 1 + c_2 &= 0 \\
+    -e^{-10} + c_1 + c_2 &= 0 \\
+    c_1 &= e^{-10} - c_2\\
+    c_1 &= e^{-10} - 1\\
+\end{align*}
+
+Using the values that we found for $c_1$ and $c_2$, we get
+
+\begin{align*}
+    u(x) &= -e^{-10x} + (e^{-10} - 1) x + 1 \\
+         &= 1 - (1 - e^{-10})x - e^{-10x}
+\end{align*}

latex/problems/problem10.tex

@@ -0,0 +1,8 @@
+\section*{Problem 10}
+
+% Time and show result, and link to relevant files
+\begin{figure}[H]
+    \centering 
+    \includegraphics[width=0.7\linewidth]{images/problem10.pdf}
+    \caption{Timing of general algorithm vs. special for step sizes $n_{steps}$}
+\end{figure}

latex/problems/problem2.tex

@@ -0,0 +1,13 @@
+\section*{Problem 2}
+
+The code for generating the points and plotting them can be found under.
+
+Point generator code: https://github.uio.no/FYS3150-G2-203/Project-1/blob/main/src/analyticPlot.cpp
+
+Plotting code: https://github.uio.no/FYS3150-G2-2023/Project-1/blob/main/src/analyticPlot.py
+
+\begin{figure}[H]
+    \centering
+    \includegraphics[width=0.8\linewidth]{images/analytical_solution.pdf}
+    \caption{Plot of the analytical solution $u(x)$.}
+\end{figure}

latex/problems/problem3.tex

@@ -0,0 +1,37 @@
+
+\section*{Problem 3}
+
+To derive the discretized version of the Poisson equation, we first need
+the Taylor expansion for $u(x)$ around $x$ for $x + h$ and $x - h$.
+
+\begin{align*}
+    u(x+h) &= u(x) + u'(x) h + \frac{1}{2} u''(x) h^2 + \frac{1}{6} u'''(x) h^3 + \mathcal{O}(h^4)
+\end{align*}
+
+\begin{align*}
+    u(x-h) &= u(x) - u'(x) h + \frac{1}{2} u''(x) h^2 - \frac{1}{6} u'''(x) h^3 + \mathcal{O}(h^4)
+\end{align*}
+
+If we add the equations above, we get this new equation:
+
+\begin{align*}
+    u(x+h) + u(x-h) &= 2 u(x) + u''(x) h^2 + \mathcal{O}(h^4) \\
+    u(x+h) - 2 u(x) + u(x-h) + \mathcal{O}(h^4) &= u''(x) h^2 \\
+    u''(x) &= \frac{u(x+h) - 2 u(x) + u(x-h)}{h^2} + \mathcal{O}(h^2) \\
+    u_i''(x) &= \frac{u_{i+1} - 2 u_i + u_{i-1}}{h^2} + \mathcal{O}(h^2) \\
+\end{align*}
+
+We can then replace $\frac{d^2u}{dx^2}$ with the RHS (right-hand side) of the equation:
+
+\begin{align*}
+    - \frac{d^2u}{dx^2} &= f(x) \\
+    \frac{ - u_{i+1} + 2 u_i - u_{i-1}}{h^2} + \mathcal{O}(h^2) &= f_i \\
+\end{align*}
+
+And lastly, we leave out $\mathcal{O}(h^2)$ and change $u_i$ to $v_i$ to 
+differentiate between the exact solution and the approximate solution, 
+and get the discretized version of the equation:
+
+\begin{align*}
+    \frac{ - v_{i+1} + 2 v_i - v_{i-1}}{h^2} &= 100 e^{-10x_i} \\
+\end{align*}

latex/problems/problem4.tex

@@ -0,0 +1,44 @@
+\section*{Problem 4}
+
+% Show that each iteration of the discretized version naturally creates a matrix equation.
+
+The value of $u(x_{0})$ and $u(x_{1})$ is known, using the discretized equation we solve for $\vec{v}$. This will result in a set of equations 
+\begin{align*}
+    - v_{0} + 2 v_{1} - v_{2} &= h^{2} \cdot f_{1} \\
+    - v_{1} + 2 v_{2} - v_{3} &= h^{2} \cdot f_{2} \\
+    \vdots & \\
+    - v_{m-2} + 2 v_{m-1} - v_{m} &= h^{2} \cdot f_{m-1} \\
+\end{align*}
+
+where $v_{i} = v(x_{i})$ and $f_{i} = f(x_{i})$. Rearranging the first and last equation, moving terms of known boundary values to the RHS
+\begin{align*}
+    2 v_{1} - v_{2} &= h^{2} \cdot f_{1} + v_{0} \\
+    - v_{1} + 2 v_{2} - v_{3} &= h^{2} \cdot f_{2} \\
+    \vdots & \\
+    - v_{m-2} + 2 v_{m-1} &= h^{2} \cdot f_{m-1} + v_{m} \\
+\end{align*}
+
+We now have a number of linear eqations, corresponding to the number of unknown values, which can be represented as an augmented matrix 
+\begin{align*}
+    \left[
+      \begin{matrix}
+        2v_{1} & -v_{2} & 0 & \dots & 0 \\
+        -v_{1} & 2v_{2} & -v_{3} & 0 & \\
+        0 & -v_{2} & 2v_{3} & -v_{4} & \\
+        \vdots  & & & \ddots & \vdots \\
+        0 & & & -v_{m-2} & 2v_{m-1} \\
+      \end{matrix}
+      \left|
+        \,
+        \begin{matrix}
+          g_{1}  \\
+          g_{2}  \\
+          g_{2}  \\
+          \vdots  \\
+          g_{m-1}  \\
+        \end{matrix}
+      \right.
+    \right]
+    \end{align*}
+    Since the boundary values are equal to $0$ the RHS can be renamed $g_{i} = h^{2} f_{i}$ for all $i$. An augmented matrix can be represented as $\boldsymbol{A} \vec{x} = \vec{b}$. In this case $\boldsymbol{A}$ is the coefficient matrix with a tridiagonal signature $(-1, 2, -1)$ and dimension $n \cross n$, where $n=m-2$.
+

latex/problems/problem5.tex

@@ -0,0 +1,6 @@
+
+\section*{Problem 5}
+
+\subsection*{a \& b)}
+
+$n = m - 2$ since when solving for $\vec{v}$, we are finding the solutions for all the points that are in between the boundaries and not the boundaries themselves. $\vec{v}^*$ on the other hand includes the boundary points.

latex/problems/problem6.tex

@@ -0,0 +1,47 @@
+\section*{Problem 6}
+
+\subsection*{a)}
+% Use Gaussian elimination, and then use backwards substitution to solve the equation
+Renaming the sub-, main-, and supdiagonal of matrix $\boldsymbol{A}$ 
+\begin{align*}
+    \vec{a} &= [a_{2}, a_{3}, ..., a_{n-1}, a_{n}] \\
+    \vec{b} &= [b_{1}, b_{2}, b_{3}, ..., b_{n-1}, b_{n}] \\
+    \vec{c} &= [c_{1}, c_{2}, c_{3}, ..., c_{n-1}] \\
+\end{align*}
+
+Following Thomas algorithm for gaussian elimination, we first perform a forward sweep followed by a backward sweep to obtain $\vec{v}$
+\begin{algorithm}[H]
+    \caption{General algorithm}\label{algo:general}
+    \begin{algorithmic}
+        \Procedure{Forward sweep}{$\vec{a}$, $\vec{b}$, $\vec{c}$}
+        \State $n \leftarrow$ length of $\vec{b}$
+        \State $\vec{\hat{b}}$, $\vec{\hat{g}} \leftarrow$ vectors of length $n$.
+        \State $\hat{b}_{1} \leftarrow b_{1}$ \Comment{Handle first element in main diagonal outside loop}
+        \State $\hat{g}_{1} \leftarrow g_{1}$
+        \For{$i = 2, 3, ..., n$}
+        \State $d \leftarrow \frac{a_{i}}{\hat{b}_{i-1}}$ \Comment{Calculating common expression}
+        \State $\hat{b}_{i} \leftarrow b_{i} - d \cdot c_{i-1}$
+        \State $\hat{g}_{i} \leftarrow g_{i} - d \cdot \hat{g}_{i-1}$
+        \EndFor
+        \Return $\vec{\hat{b}}$, $\vec{\hat{g}}$
+        \EndProcedure
+
+        \Procedure{Backward sweep}{$\vec{\hat{b}}$, $\vec{\hat{g}}$}
+        \State $n \leftarrow$ length of $\vec{\hat{b}}$
+        \State $\vec{v} \leftarrow$ vector of length $n$.
+        \State $v_{n} \leftarrow \frac{\hat{g}_{n}}{\hat{b}_{n}}$
+        \For{$i = n-1, n-2, ..., 1$}
+        \State $v_{i} \leftarrow \frac{\hat{g}_{i} - c_{i} \cdot v_{i+1}}{\hat{b}_{i}}$ 
+        \EndFor
+        \Return $\vec{v}$
+        \EndProcedure
+    \end{algorithmic}
+\end{algorithm}
+
+
+\subsection*{b)}
+% Figure out FLOPs
+Counting the number of FLOPs for the general algorithm by looking at one procedure at a time. 
+For every iteration of i in forward sweep we have 1 division, 2 multiplications, and 2 subtractions, resulting in $5(n-1)$ FLOPs. 
+For backward sweep we have 1 division, and for every iteration of i we have 1 subtraction, 1 multiplication, and 1 division, resulting in $3(n-1)+1$ FLOPs. 
+Total FLOPs for the general algorithm is $8(n-1)+1$.

latex/problems/problem7.tex

@@ -0,0 +1,14 @@
+\section*{Problem 7}
+
+\subsection*{a)}
+% Link to relevant files on gh and possibly add some comments
+The code can be found at https://github.uio.no/FYS3150-G2-2023/Project-1/blob/main/src/generalAlgorithm.cpp
+
+\subsection*{b)}
+Increasing the number of steps results in an approximation close to the exact solution.
+\begin{figure}[H]
+    \centering 
+    \includegraphics[width=0.8\linewidth]{images/problem7.pdf}
+    \caption{Plot showing the numeric solution of $u_{approx}$ for $n_{steps}$ and the exact solution $u_{exact}$.}
+\end{figure}
+

latex/problems/problem8.tex

@@ -0,0 +1,27 @@
+\section*{Problem 8}
+
+%link to relvant files and show plots
+\subsection*{a)}
+Increasing number of steps result in a decrease of absolute error.
+\begin{figure}[H]
+    \centering 
+    \includegraphics[width=0.8\linewidth]{images/problem8_a.pdf}
+    \caption{Absolute error for different step sizes $n_{steps}$.}
+\end{figure}
+
+\subsection*{b)}
+Increasing number of steps also result in a decrease of absolute error.
+\begin{figure}[H]
+    \centering 
+    \includegraphics[width=0.8\linewidth]{images/problem8_b.pdf}
+    \caption{Relative error for different step sizes $n_{steps}$.}
+\end{figure}
+
+\subsection*{c)}
+Increasing number of steps result in a decrease of maximum relative error, up to a certain number of steps. At $n_{steps} \approx 10^{5}$ the maximumrelative error increase. 
+This can be related to loss of numerical precicion when step size is small.
+\begin{figure}[H]
+    \centering 
+    \includegraphics[width=0.7\linewidth]{images/problem8_c.pdf}
+    \caption{Maximum relative error for each step sizes $n_{steps}$.}
+\end{figure}

latex/problems/problem9.tex

@@ -0,0 +1,55 @@
+\section*{Problem 9}
+
+\subsection*{a)}
+% Specialize algorithm
+The special algorithm does not require the values of all $a_{i}$, $b_{i}$, $c_{i}$. 
+We find the values of $\hat{b}_{i}$ from simplifying the general case
+\begin{align*}
+    \hat{b}_{i} &= b_{i} - \frac{a_{i} \cdot c_{i-1}}{\hat{b}_{i-1}} \\
+    \hat{b}_{i} &= 2 - \frac{1}{\hat{b}_{i-1}}
+\end{align*} 
+Calculating the first values to see a pattern 
+\begin{align*}
+    \hat{b}_{1} &= 2 \\
+    \hat{b}_{2} &= 2 - \frac{1}{2} = \frac{3}{2} \\
+    \hat{b}_{3} &= 2 - \frac{1}{\frac{3}{2}} = \frac{4}{3} \\
+    \hat{b}_{4} &= 2 - \frac{1}{\frac{4}{3}} = \frac{5}{4} \\
+    \vdots & \\
+    \hat{b}_{i} &= \frac{i+1}{i} && \text{for $i = 1, 2, ..., n$}
+\end{align*}
+
+
+\begin{algorithm}[H]
+    \caption{Special algorithm}\label{algo:special}
+    \begin{algorithmic}
+        \Procedure{Forward sweep}{$\vec{b}$}
+        \State $n \leftarrow$ length of $\vec{b}$
+        \State $\vec{\hat{b}}$, $\vec{\hat{g}} \leftarrow$ vectors of length $n$.
+        \State $\hat{b}_{1} \leftarrow 2$ \Comment{Handle first element in main diagonal outside loop}
+        \State $\hat{g}_{1} \leftarrow g_{1}$
+        \For{$i = 2, 3, ..., n$}
+        \State $\hat{b}_{i} \leftarrow \frac{i+1}{i}$
+        \State $\hat{g}_{i} \leftarrow g_{i} + \frac{\hat{g}_{i-1}}{\hat{b}_{i-1}}$
+        \EndFor
+        \Return $\vec{\hat{b}}$, $\vec{\hat{g}}$
+        \EndProcedure
+
+        \Procedure{Backward sweep}{$\vec{\hat{b}}$, $\vec{\hat{g}}$}
+        \State $n \leftarrow$ length of $\vec{\hat{b}}$
+        \State $\vec{v} \leftarrow$ vector of length $n$.
+        \State $v_{n} \leftarrow \frac{\hat{g}_{n}}{\hat{b}_{n}}$
+        \For{$i = n-1, n-2, ..., 1$}
+        \State $v_{i} \leftarrow \frac{\hat{g}_{i} + v_{i+1}}{\hat{b}_{i}}$ 
+        \EndFor
+        \Return $\vec{v}$
+        \EndProcedure
+    \end{algorithmic}
+\end{algorithm}
+
+
+\subsection*{b)}
+% Find FLOPs
+For every iteration of i in forward sweep we have 2 divisions, and 2 additions, resulting in $4(n-1)$ FLOPs. 
+For backward sweep we have 1 division, and for every iteration of i we have 1 addition, and 1 division, resulting in $2(n-1)+1$ FLOPs. 
+Total FLOPs for the special algorithm is $6(n-1)+1$.
+

src/Makefile

@@ -0,0 +1,30 @@
+CC=g++
+
+CCFLAGS= -std=c++11
+
+OBJS=generalAlgorithm.o specialAlgorithm.o funcs.o
+
+EXEC=main analyticPlot
+
+.PHONY: clean create_dirs
+
+all: create_dirs $(EXEC)
+
+main: main.o $(OBJS)
+	$(CC) $(CCFLAGS) -o $@ $^
+
+analyticPlot: analyticPlot.o
+	$(CC) $(CCFLAGS) -o $@ $^
+
+%.o: %.cpp
+	$(CC) $(CCFLAGS) -c -o $@ $^
+
+clean:
+	rm *.o
+	rm $(EXEC)
+	rm -r output
+
+create_dirs:
+	mkdir -p output/general
+	mkdir -p output/special
+	mkdir -p output/error

src/analyticPlot.cpp

@@ -0,0 +1,55 @@
+#include <iostream>
+#include <cmath>
+#include <vector>
+#include <string>
+#include <numeric>
+#include <fstream>
+#include <iomanip>
+
+#define RANGE 1000
+#define FILENAME "output/analytical_solution.txt"
+
+double u(double x);
+void generate_range(std::vector<double> &vec, double start, double stop, int n);
+void write_analytical_solution(std::string filename, int n);
+
+int main() {
+    write_analytical_solution(FILENAME, RANGE);
+
+    return 0;
+};
+
+double u(double x) {
+    return 1 - (1 - exp(-10))*x - exp(-10*x);
+};
+
+void generate_range(std::vector<double> &vec, double start, double stop, int n) {
+    double step = (stop - start) / n;
+
+    for (int i = 0; i <= vec.size(); i++) {
+        vec[i] = i * step;
+    }
+}
+
+void write_analytical_solution(std::string filename, int n) {
+    std::vector<double> x(n), y(n);
+    generate_range(x, 0.0, 1.0, n);
+
+    // Set up output file and strem
+    std::ofstream outfile;
+    outfile.open(filename);
+
+    // Parameters for formatting 
+    int width = 12;
+    int prec = 4;
+
+    // Calculate u(x) and write to file
+    for (int i = 0; i <= x.size(); i++) {
+        y[i] = u(x[i]);
+        outfile << std::setw(width) << std::setprecision(prec) << std::scientific << x[i] 
+            << std::setw(width) << std::setprecision(prec) << std::scientific << y[i] 
+            << std::endl;
+    }
+    outfile.close();
+}
+

src/funcs.cpp

@@ -0,0 +1,38 @@
+#include "funcs.hpp"
+
+double f(double x) {
+    return 100*std::exp(-10*x);
+}
+
+double u(double x) {
+    return 1. - (1. - std::exp(-10.))*x - std::exp(-10.*x);
+}
+
+void build_g_vec(int n_steps, arma::vec& g_vec) {
+    g_vec.resize(n_steps-1);
+
+    double step_size = 1./ (double) n_steps;
+    for (int i=0; i < n_steps-1; i++) {
+        g_vec(i) = step_size*step_size*f((i+1)*step_size);
+    }
+}
+
+void build_arrays(
+    int n_steps, 
+    arma::vec& sub_diag, 
+    arma::vec& main_diag, 
+    arma::vec& sup_diag,
+    arma::vec& g_vec
+) 
+{
+    sub_diag.resize(n_steps-2);
+    main_diag.resize(n_steps-1);
+    sup_diag.resize(n_steps-2);
+
+    sub_diag.fill(-1);
+    main_diag.fill(2);
+    sup_diag.fill(-1);
+
+    build_g_vec(n_steps, g_vec);
+}
+

src/funcs.hpp

@@ -0,0 +1,24 @@
+#ifndef __FUNCS__
+#define __FUNCS__
+
+#include <armadillo>
+#include <cmath>
+
+#define PRECISION 8
+
+#define N_STEPS_EXP 7
+
+double f(double x);
+
+double u(double x);
+
+void build_g_vec(int n_steps, arma::vec& g_vec);
+
+void build_arrays(
+    int n_steps, 
+    arma::vec& sub_diag, 
+    arma::vec& main_diag, 
+    arma::vec& sup_diag,
+    arma::vec& g_vec
+);
+#endif

src/generalAlgorithm.cpp

@@ -0,0 +1,84 @@
+#include "funcs.hpp"
+#include "generalAlgorithm.hpp"
+#include <cmath>
+
+arma::vec& general_algorithm(
+    arma::vec& sub_diag,
+    arma::vec& main_diag,
+    arma::vec& sup_diag,
+    arma::vec& g_vec
+)
+{
+    int n = g_vec.n_elem;
+    double d;
+
+    for (int i = 1; i < n; i++) {
+        d = sub_diag(i-1) / main_diag(i-1);
+        main_diag(i) -= d*sup_diag(i-1);
+        g_vec(i) -= d*g_vec(i-1);
+    }
+
+    g_vec(n-1) /= main_diag(n-1);
+
+    for (int i = n-2; i >= 0; i--) {
+        g_vec(i) = (g_vec(i) - sup_diag(i) * g_vec(i+1)) / main_diag(i);
+    }
+    return g_vec;
+}
+
+void general_algorithm_main() 
+{
+    arma::vec sub_diag, main_diag, sup_diag, g_vec, v_vec;
+    std::ofstream ofile;
+    int steps;
+    double step_size;
+
+    for (int i = 0; i < N_STEPS_EXP; i++) {
+        steps = std::pow(10, i+1);
+        step_size = 1./(double) steps;
+
+        build_arrays(steps, sub_diag, main_diag, sup_diag, g_vec);
+
+        v_vec = general_algorithm(sub_diag, main_diag, sup_diag, g_vec);
+
+        ofile.open("output/general/out_" + std::to_string(steps) + ".txt");
+
+        for (int j=0; j < v_vec.n_elem; j++) {
+            ofile << std::setprecision(PRECISION) << std::scientific << step_size*(j+1) << ","
+                << std::setprecision(PRECISION) << std::scientific << v_vec(j) << std::endl;
+        }
+        ofile.close();
+    }
+}
+
+void general_algorithm_error() 
+{
+    arma::vec sub_diag, main_diag, sup_diag, g_vec, v_vec;
+    std::ofstream ofile;
+    int steps;
+    double step_size, abs_err, rel_err, u_i, v_i;
+
+    for (int i=0; i < N_STEPS_EXP; i++) {
+        steps = std::pow(10, i+1);
+        step_size = 1./(double) steps;
+
+        build_arrays(steps, sub_diag, main_diag, sup_diag, g_vec);
+
+        v_vec = general_algorithm(sub_diag, main_diag, sup_diag, g_vec);
+
+        ofile.open("output/error/out_" + std::to_string(steps) + ".txt");
+
+        for (int j=0; j < v_vec.n_elem; j++) {
+            u_i = u(step_size*(j+1));
+            v_i = v_vec(j);
+            abs_err = u_i - v_i;
+            ofile << std::setprecision(PRECISION) << std::scientific 
+                << step_size*(j+1) << ","
+                << std::setprecision(PRECISION) << std::scientific 
+                << std::log10(std::abs(abs_err)) << ","
+                << std::setprecision(PRECISION) << std::scientific 
+                << std::log10(std::abs(abs_err/u_i)) << std::endl;
+        }
+        ofile.close();
+    }
+}

src/generalAlgorithm.hpp

@@ -0,0 +1,18 @@
+#ifndef __GENERAL_ALG__
+#define __GENERAL_ALG__
+
+#include <armadillo>
+#include <iomanip>
+
+arma::vec& general_algorithm(
+    arma::vec& sub_diag,
+    arma::vec& main_diag,
+    arma::vec& sup_diag,
+    arma::vec& g_vec
+);
+
+void general_algorithm_main();
+
+void general_algorithm_error();
+
+#endif

src/main.cpp

@@ -0,0 +1,68 @@
+#include <armadillo>
+#include <cmath>
+#include <ctime>
+#include <fstream>
+#include <iomanip>
+#include <ios>
+#include <string>
+
+#include "funcs.hpp"
+#include "generalAlgorithm.hpp"
+#include "specialAlgorithm.hpp"
+
+#define TIMING_ITERATIONS 5
+
+void timing() {
+    arma::vec sub_diag, main_diag, sup_diag, g_vec;
+    int n_steps;
+
+    std::ofstream ofile;
+    ofile.open("output/timing.txt");
+
+    // Timing
+    for (int i=1; i < N_STEPS_EXP; i++) {
+        n_steps = std::pow(10, i);
+        clock_t g_1, g_2, s_1, s_2;
+        double g_res = 0, s_res = 0;
+      
+        // Repeat a number of times to take an average
+        for (int j=0; j < TIMING_ITERATIONS; j++) {
+            
+            build_arrays(n_steps, sub_diag, main_diag, sup_diag, g_vec);
+
+            g_1 = clock();
+
+            general_algorithm(sub_diag, main_diag, sup_diag, g_vec);
+
+            g_2 = clock();
+
+            g_res += (double) (g_2 - g_1) / CLOCKS_PER_SEC;
+            // Rebuild g_vec for the special alg
+            build_g_vec(n_steps, g_vec);
+
+            s_1 = clock();
+
+            special_algorithm(-1., 2., -1., g_vec);
+
+            s_2 = clock();
+
+            s_res += (double) (s_2 - s_1) / CLOCKS_PER_SEC;
+
+        }
+        // Write the average time to file
+        ofile 
+            << n_steps << ","
+            << g_res / (double) TIMING_ITERATIONS << ","
+            << s_res / (double) TIMING_ITERATIONS << std::endl;
+    }
+
+    ofile.close();
+}
+
+int main() 
+{
+    timing();
+    general_algorithm_main();
+    general_algorithm_error();
+    special_algorithm_main();
+}

src/plot_analytic_solution.py

@@ -0,0 +1,19 @@
+import numpy as np 
+import matplotlib.pyplot as plt
+
+def main():
+    FILENAME = "../latex/images/analytical_solution.pdf"
+    x = []
+    v = []
+
+    with open('output/analytical_solution.txt') as f:
+        for line in f:
+            a, b = line.strip().split()
+            x.append(float(a))
+            v.append(float(b))
+
+    plt.plot(x, v)
+    plt.savefig(FILENAME)
+
+if __name__ == "__main__":
+    main()

src/plot_general_alg.py

@@ -0,0 +1,30 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+analytical_func = lambda x: 1 - (1 - np.exp(-10))*x - np.exp(-10*x)
+
+def main():
+    for i in range(7):
+        x = []
+        y = []
+        x.append(0.)
+        y.append(0.)
+        with open(f"output/general/out_{10**(i+1)}.txt", "r") as f:
+            lines = f.readlines()
+            for line in lines:
+                x_i, y_i = line.strip().split(",")
+                x.append(float(x_i))
+                y.append(float(y_i))
+
+        x.append(1.)
+        y.append(0.)
+
+        plt.plot(x, y, label=f"n$_{{steps}} = 10^{i+1}$")
+
+    x = np.linspace(0, 1, 1001)
+    plt.plot(x, analytical_func(x), label=f"u$_{{exact}}$")
+    plt.legend()
+    plt.savefig("../latex/images/problem7.pdf")
+
+if __name__ == "__main__":
+    main()

src/plot_general_alg_error.py

@@ -0,0 +1,40 @@
+import matplotlib.pyplot as plt
+
+# plt.rc('text', usetex=True)
+# plt.rc('font', family='serif')
+
+def main():
+    for i in range(7):
+        x = []
+        abs_err = []
+        rel_err = []
+        with open(f"output/error/out_{10**(i+1)}.txt", "r") as f:
+            lines = f.readlines()
+            for line in lines:
+                x_i, abs_err_i, rel_err_i = line.strip().split(",")
+                x.append(float(x_i))
+                abs_err.append(float(abs_err_i))
+                rel_err.append(float(rel_err_i))
+
+        plt.figure(1)
+        plt.plot(x, abs_err, label=f"n$_{{steps}} = 10^{i+1}$")
+        plt.figure(2)
+        plt.plot(x, rel_err, label=f"n$_{{steps}} = 10^{i+1}$")
+
+        plt.figure(3)
+        plt.plot(i+1, max(rel_err), marker="o", markersize=10)
+    
+    plt.figure(1)
+    plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
+    plt.figure(2)
+    plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
+    
+    plt.figure(1)
+    plt.savefig("../latex/images/problem8_a.pdf", bbox_inches="tight")
+    plt.figure(2)
+    plt.savefig("../latex/images/problem8_b.pdf", bbox_inches="tight")
+    plt.figure(3)
+    plt.savefig("../latex/images/problem8_c.pdf", bbox_inches="tight")
+
+if __name__ == "__main__":
+    main()

src/specialAlgorithm.cpp

@@ -0,0 +1,52 @@
+#include "funcs.hpp"
+#include "specialAlgorithm.hpp"
+
+arma::vec& special_algorithm(
+    double sub_sig,
+    double main_sig,
+    double sup_sig,
+    arma::vec& g_vec
+)
+{
+    int n = g_vec.n_elem;
+    arma::vec diag = arma::vec(n);
+
+    for (int i = 1; i < n; i++) {
+        // Calculate values for main diagonal based on indices
+        diag(i-1) = (double)(i+1) / i;
+        g_vec(i) += g_vec(i-1) / diag(i-1);
+    }
+    // The last element in main diagonal has value (i+1)/i = (n+1)/n
+    g_vec(n-1) /= (double)(n+1) / (n);
+
+    for (int i = n-2; i >= 0; i--) {
+        g_vec(i) = (g_vec(i) + g_vec(i+1))/ diag(i);
+    }
+
+    return g_vec;
+}
+
+void special_algorithm_main() 
+{
+    arma::vec g_vec, v_vec;
+    std::ofstream ofile;
+    int steps;
+    double sub_sig, main_sig, sup_sig, step_size;
+
+    for (int i = 0; i < N_STEPS_EXP; i++) {
+        steps = std::pow(10, i+1);
+        step_size = 1./(double) steps;
+        build_g_vec(steps, g_vec);
+
+        v_vec = special_algorithm(sub_sig, main_sig, sup_sig, g_vec);
+
+        ofile.open("output/special/out_" + std::to_string(steps) + ".txt");
+
+        for (int j=0; j < v_vec.n_elem; j++) {
+            ofile << std::setprecision(PRECISION) << std::scientific << step_size*(j+1) << ","
+                << std::setprecision(PRECISION) << std::scientific << v_vec(j) << std::endl;
+        }
+        ofile.close();
+    }
+}
+

src/specialAlgorithm.hpp

@@ -0,0 +1,16 @@
+#ifndef __SPECIAL_ALG__
+#define __SPECIAL_ALG__
+
+#include <armadillo>
+#include <iomanip>
+
+arma::vec& special_algorithm(
+    double sub_sig,
+    double main_sig,
+    double sup_sig,
+    arma::vec& g_vec
+);
+
+void special_algorithm_main();
+
+#endif

src/timing.py

@@ -0,0 +1,22 @@
+import matplotlib.pyplot as plt
+
+def main():
+    x = []
+    gen_alg = []
+    spec_alg = []
+    with open(f"output/timing.txt", "r") as f:
+        lines = f.readlines()
+        for line in lines:
+            x_i, gen_i, spec_i = line.strip().split(",")
+            x.append(float(x_i))
+            gen_alg.append(float(gen_i))
+            spec_alg.append(float(spec_i))
+
+        plt.plot(x, gen_alg, label=f"General algorithm")
+        plt.plot(x, spec_alg, label=f"Special algorithm")
+
+    plt.legend()
+    plt.savefig("../latex/images/problem10.pdf")
+
+if __name__ == "__main__":
+    main()