Compare commits

..

No commits in common. "main" and "coryab/implement-problem-2" have entirely different histories.

35 changed files with 204 additions and 591 deletions

1
.gitignore vendored
View File

@ -43,5 +43,4 @@
src/* src/*
!src/Makefile !src/Makefile
!src/*.cpp !src/*.cpp
!src/*.hpp
!src/*.py !src/*.py

View File

@ -3,36 +3,3 @@
## Practical information ## Practical information
- [Project](https://anderkve.github.io/FYS3150/book/projects/project1.html) - [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**.

Binary file not shown.

View File

@ -80,7 +80,7 @@
\begin{document} \begin{document}
\title{Project 1} % self-explanatory \title{Project 1} % self-explanatory
\author{Cory Alexander Balaton \& Janita Ovidie Sandtrøen Willumsen} % self-explanatory \author{Cory Balaton \& Janita Willumsen} % self-explanatory
\date{\today} % self-explanatory \date{\today} % self-explanatory
\noaffiliation % ignore this, but keep it. \noaffiliation % ignore this, but keep it.
@ -107,6 +107,4 @@
\input{problems/problem9} \input{problems/problem9}
\input{problems/problem10}
\end{document} \end{document}

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

View File

@ -40,5 +40,5 @@ Using the values that we found for $c_1$ and $c_2$, we get
\begin{align*} \begin{align*}
u(x) &= -e^{-10x} + (e^{-10} - 1) x + 1 \\ u(x) &= -e^{-10x} + (e^{-10} - 1) x + 1 \\
&= 1 - (1 - e^{-10})x - e^{-10x} &= 1 - (1 - e^{-10}) - e^{-10x}
\end{align*} \end{align*}

View File

@ -1,8 +1,3 @@
\section*{Problem 10} \section*{Problem 10}
% Time and show result, and link to relevant files % 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}

View File

@ -6,8 +6,6 @@ Point generator code: https://github.uio.no/FYS3150-G2-203/Project-1/blob/main/s
Plotting code: https://github.uio.no/FYS3150-G2-2023/Project-1/blob/main/src/analyticPlot.py Plotting code: https://github.uio.no/FYS3150-G2-2023/Project-1/blob/main/src/analyticPlot.py
\begin{figure}[H] Here is the plot of the analytical solution for $u(x)$.
\centering
\includegraphics[width=0.8\linewidth]{images/analytical_solution.pdf} \includegraphics[scale=.5]{analytical_solution.pdf}
\caption{Plot of the analytical solution $u(x)$.}
\end{figure}

View File

@ -2,7 +2,7 @@
% Show that each iteration of the discretized version naturally creates a matrix equation. % 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 The value of $u(x_{0})$ and $u(x_{1})$ is known, using the discretized equation we can approximate the value of $f(x_{i}) = f_{i}$. This will result in a set of equations
\begin{align*} \begin{align*}
- v_{0} + 2 v_{1} - v_{2} &= h^{2} \cdot f_{1} \\ - v_{0} + 2 v_{1} - v_{2} &= h^{2} \cdot f_{1} \\
- v_{1} + 2 v_{2} - v_{3} &= h^{2} \cdot f_{2} \\ - v_{1} + 2 v_{2} - v_{3} &= h^{2} \cdot f_{2} \\
@ -10,7 +10,7 @@ The value of $u(x_{0})$ and $u(x_{1})$ is known, using the discretized equation
- v_{m-2} + 2 v_{m-1} - v_{m} &= h^{2} \cdot f_{m-1} \\ - v_{m-2} + 2 v_{m-1} - v_{m} &= h^{2} \cdot f_{m-1} \\
\end{align*} \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 Rearranging the first and last equation, moving terms of known boundary values to the RHS
\begin{align*} \begin{align*}
2 v_{1} - v_{2} &= h^{2} \cdot f_{1} + v_{0} \\ 2 v_{1} - v_{2} &= h^{2} \cdot f_{1} + v_{0} \\
- v_{1} + 2 v_{2} - v_{3} &= h^{2} \cdot f_{2} \\ - v_{1} + 2 v_{2} - v_{3} &= h^{2} \cdot f_{2} \\
@ -40,5 +40,5 @@ We now have a number of linear eqations, corresponding to the number of unknown
\right. \right.
\right] \right]
\end{align*} \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$. where $g_{i} = h^{2} f_{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$.

View File

@ -1,6 +1,6 @@
\section*{Problem 5} \section*{Problem 5}
\subsection*{a \& b)} \subsection*{a)}
$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. \subsection*{b)}

View File

@ -1,47 +1,9 @@
\section*{Problem 6} \section*{Problem 6}
\subsection*{a)} \subsection*{a)}
% Use Gaussian elimination, and then use backwards substitution to solve the equation % 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)} \subsection*{b)}
% Figure out FLOPs
Counting the number of FLOPs for the general algorithm by looking at one procedure at a time. % Figure it out
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$.

View File

@ -1,14 +1,3 @@
\section*{Problem 7} \section*{Problem 7}
\subsection*{a)}
% Link to relevant files on gh and possibly add some comments % 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}

View File

@ -1,27 +1,3 @@
\section*{Problem 8} \section*{Problem 8}
%link to relvant files and show plots %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}

View File

@ -1,55 +1,3 @@
\section*{Problem 9} \section*{Problem 9}
\subsection*{a)} % Show the algorithm, then calculate FLOPs, then link to relevant files
% 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$.

View File

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

View File

@ -7,7 +7,7 @@
#include <iomanip> #include <iomanip>
#define RANGE 1000 #define RANGE 1000
#define FILENAME "output/analytical_solution.txt" #define FILENAME "analytical_solution.txt"
double u(double x); double u(double x);
void generate_range(std::vector<double> &vec, double start, double stop, int n); void generate_range(std::vector<double> &vec, double start, double stop, int n);

View File

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

View File

@ -1,38 +0,0 @@
#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);
}

View File

@ -1,24 +0,0 @@
#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

View File

@ -1,84 +0,0 @@
#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();
}
}

View File

@ -1,18 +0,0 @@
#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

View File

@ -1,68 +1,24 @@
#include "GeneralAlgorithm.hpp"
#include <armadillo> #include <armadillo>
#include <cmath> #include <iostream>
#include <ctime>
#include <fstream>
#include <iomanip>
#include <ios>
#include <string>
#include "funcs.hpp" double f(double x) {
#include "generalAlgorithm.hpp" return 100. * std::exp(-10.*x);
#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() double a_sol(double x) {
{ return 1. - (1. - std::exp(-10)) * x - std::exp(-10*x);
timing(); }
general_algorithm_main();
general_algorithm_error(); int main() {
special_algorithm_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;
} }

View File

@ -1,30 +0,0 @@
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()

View File

@ -1,40 +0,0 @@
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()

0
src/problem2.cpp Normal file
View File

162
src/simpleFile.cpp Normal file
View File

@ -0,0 +1,162 @@
#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();
}

View File

@ -1,52 +0,0 @@
#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();
}
}

View File

@ -1,16 +0,0 @@
#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

View File

@ -1,22 +0,0 @@
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()