Wrote a rough draft of problem 5

.gitignore

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

latex/assignment_1.tex

@@ -1,9 +1,7 @@
 \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
 
@@ -15,7 +13,7 @@
 %% I recommend downloading TeXMaker, because it includes a large library of the most common packages.
 
 \usepackage{physics,amssymb}  % mathematical symbols (physics imports amsmath)
-\usepackage{amsmath}
+\include{amsmath}
 \usepackage{graphicx}         % include graphics such as plots
 \usepackage{xcolor}           % set colors
 \usepackage{hyperref}         % automagic cross-referencing (this is GODLIKE)
@@ -74,9 +72,6 @@
 %%
 %% Don't ask me why, I don't know.
 
-% custom stuff
-\graphicspath{{./images/}}
-
 \begin{document}
 
 \title{Project 1}      % self-explanatory
@@ -89,22 +84,90 @@
     
 \textit{https://github.uio.no/FYS3150-G2-2023/Project-1}
     
-\input{problems/problem1}
+\section*{Problem 1}
 
-\input{problems/problem2}
+% 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/problem3}
+Using the boundary conditions, we can find $c_1$ and $c_2$ as shown below:
 
-\input{problems/problem4}
+\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/problem5}
+\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/problem6}
+Using the values that we found for $c_1$ and $c_2$, we get
 
-\input{problems/problem7}
+\begin{align*}
+    u(x) &= -e^{-10x} + (e^{-10} - 1) x + 1 \\
+         &= 1 - (1 - e^{-10}) - e^{-10x}
+\end{align*}
 
-\input{problems/problem8}
+\section*{Problem 2}
 
-\input{problems/problem9}
+% Write which .cpp/.hpp/.py (using a link?) files are relevant for this and show the plot generated.
+
+\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)}
+
+% Phrase it better
+$n = m - 2$, since $\textbf{A}$ is used to solve for all of the points in 
+between the end-points $0$ and $1$. For the complete solution, we need to add
+$u(0)$ and $u(1)$.
+
+\subsection*{b)}
+
+When solving for $\vec{v}$, we find the approximate solutions for $u(x)$
+that are in between the end-points, but not the end-points themselves.
+
+\section*{Problem 6}
+
+\subsection*{a)}
+
+% Use Gaussian elimination, and then use backwards substitution to solve the equation
+
+\subsection*{b)}
+
+% Figure it out
+
+% Linelevant files on gh and possibly add some comments
+
+\section*{Problem 8}
+
+%link to relvant files and show plots
+
+\section*{Problem 9}
+
+% Shon*{Proalgorithm, then calculate FLOPs, then link to relevant files
+
+\section*{Problem 10}
+
+% Time and show result, and link to relevant files
 
 \end{document}

latex/problems/problem1.tex

@@ -1,44 +0,0 @@
-\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}) - e^{-10x}
-\end{align*}

latex/problems/problem10.tex

@@ -1,3 +0,0 @@
-\section*{Problem 10}
-
-% Time and show result, and link to relevant files

latex/problems/problem2.tex

@@ -1,11 +0,0 @@
-\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
-
-Here is the plot of the analytical solution for $u(x)$.
-
-\includegraphics[scale=.5]{analytical_solution.pdf}

latex/problems/problem3.tex

@@ -1,37 +0,0 @@
-
-\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

@@ -1,44 +0,0 @@
-\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

@@ -1,6 +0,0 @@
-
-\section*{Problem 5}
-
-\subsection*{a)}
-
-\subsection*{b)}

latex/problems/problem6.tex

@@ -1,47 +0,0 @@
-\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

@@ -1,3 +0,0 @@
-\section*{Problem 7}
-
-% Link to relevant files on gh and possibly add some comments

latex/problems/problem8.tex

@@ -1,3 +0,0 @@
-\section*{Problem 8}
-
-%link to relvant files and show plots

latex/problems/problem9.tex

@@ -1,55 +0,0 @@
-\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

@@ -1,17 +0,0 @@
-CC=g++
-
-.PHONY: clean
-
-all: simpleFile analyticPlot
-
-simpleFile: simpleFile.o
-	$(CC) -o $@ $^
-
-analyticPlot: analyticPlot.o
-	$(CC) -o $@ $^
-
-%.o: %.cpp
-	$(CC) -c $< -o $@
-
-clean:
-	rm *.o

src/analyticPlot.cpp

@@ -1,55 +0,0 @@
-#include <iostream>
-#include <cmath>
-#include <vector>
-#include <string>
-#include <numeric>
-#include <fstream>
-#include <iomanip>
-
-#define RANGE 1000
-#define FILENAME "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/analyticPlot.py

@@ -1,19 +0,0 @@
-import numpy as np 
-import matplotlib.pyplot as plt
-
-def main():
-    FILENAME = "analytical_solution.pdf"
-    x = []
-    v = []
-
-    with open('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/main.cpp

@@ -1,24 +0,0 @@
-#include "GeneralAlgorithm.hpp"
-#include <armadillo>
-#include <iostream>
-
-double f(double x) {
-    return 100. * std::exp(-10.*x);
-}
-
-double a_sol(double x) {
-    return 1. - (1. - std::exp(-10)) * x - std::exp(-10*x);
-}
-
-int main() {
-    arma::mat A = arma::eye(3,3);
-
-    GeneralAlgorithm ga(3, &A, f, a_sol, 0., 1.);
-
-    ga.solve();
-    std::cout << "Time: " << ga.time(5) << std::endl;
-    ga.error();
-
-    return 0;
-
-}

src/simpleFile.cpp

@@ -1,162 +0,0 @@
-#include <armadillo>
-#include <cmath>
-#include <ctime>
-#include <fstream>
-#include <iomanip>
-#include <ios>
-#include <string>
-
-#define TIMING_ITERATIONS 5
-
-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;
-}
-
-
-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 error(
-    std::string filename,
-    arma::vec* x_vec,
-    arma::vec* v_vec,
-    arma::vec* a_vec
-)
-{
-    std::ofstream ofile;
-    ofile.open(filename);
-
-    if (!ofile.is_open()) {
-        exit(1);
-    }
-
-    for (int i=0; i < a_vec->n_elem; i++) {
-        double sub = (*a_vec)(i) - (*v_vec)(i);
-        ofile << std::setprecision(8) << std::scientific << (*x_vec)(i)
-            << std::setprecision(8) << std::scientific << std::log10(std::abs(sub))
-            << std::setprecision(8) << std::scientific << std::log10(std::abs(sub/(*a_vec)(i)))
-            << std::endl;
-    }
-    
-    ofile.close();
-
-}
-
-double f(double x) {
-    return 100*std::exp(-10*x);
-}
-
-void build_array(
-    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);
-
-    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) = f((i+1)*step_size);
-    }
-
-}
-
-void timing() {
-    arma::vec sub_diag, main_diag, sup_diag, g_vec;
-    int n_steps;
-
-    std::ofstream ofile;
-    ofile.open("timing.txt");
-
-    // Timing
-    for (int i=1; i <= 8; i++) {
-        n_steps = std::pow(10, i);
-        clock_t g_1, g_2, s_1, s_2;
-        double g_res = 0, s_res = 0;
-      
-        for (int j=0; j < TIMING_ITERATIONS; j++) {
-            build_array(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;
-            build_array(n_steps, &sub_diag, &main_diag, &sup_diag, &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;
-
-        }
-        ofile 
-            << n_steps << ","
-            << g_res / (double) TIMING_ITERATIONS << ","
-            << s_res / (double) TIMING_ITERATIONS << std::endl;
-    }
-
-    ofile.close();
-}
-
-int main() 
-{
-    timing();
-}