Compare commits

..

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

36 changed files with 233 additions and 643 deletions

3
.gitignore vendored
View File

@ -41,7 +41,4 @@
# C++ specifics
src/*
!src/Makefile
!src/*.cpp
!src/*.hpp
!src/*.py

View File

@ -3,36 +3,3 @@
## 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**.

Binary file not shown.

View File

@ -74,13 +74,10 @@
%%
%% Don't ask me why, I don't know.
% custom stuff
\graphicspath{{./images/}}
\begin{document}
\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
\noaffiliation % ignore this, but keep it.
@ -107,6 +104,4 @@
\input{problems/problem9}
\input{problems/problem10}
\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*}
u(x) &= -e^{-10x} + (e^{-10} - 1) x + 1 \\
&= 1 - (1 - e^{-10})x - e^{-10x}
&= 1 - (1 - e^{-10}) - e^{-10x}
\end{align*}

View File

@ -1,8 +1,3 @@
\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}

View File

@ -1,13 +1,3 @@
\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}
% Write which .cpp/.hpp/.py (using a link?) files are relevant for this and show the plot generated.

View File

@ -2,7 +2,7 @@
% 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*}
- v_{0} + 2 v_{1} - v_{2} &= h^{2} \cdot f_{1} \\
- 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} \\
\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*}
2 v_{1} - v_{2} &= h^{2} \cdot f_{1} + v_{0} \\
- 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]
\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}
\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}
\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$.
% Figure it out

View File

@ -1,14 +1,3 @@
\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}

View File

@ -1,27 +1,3 @@
\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}

View File

@ -1,55 +1,3 @@
\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$.
% Show the algorithm, then calculate FLOPs, then link to relevant files

View File

@ -1,30 +1,12 @@
CC=g++
CCFLAGS= -std=c++11
all: simpleFile
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 $@ $^
simpleFile: simpleFile.o
$(CC) -o $@ $^
%.o: %.cpp
$(CC) $(CCFLAGS) -c -o $@ $^
$(CC) -c $< -o $@
clean:
rm *.o
rm $(EXEC)
rm -r output
create_dirs:
mkdir -p output/general
mkdir -p output/special
mkdir -p output/error

View File

@ -6,36 +6,17 @@
#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);
int n = 1000;
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::string filename = "datapoints.txt";
std::ofstream outfile;
outfile.open(filename);
@ -51,5 +32,21 @@ void write_analytical_solution(std::string filename, int n) {
<< std::endl;
}
outfile.close();
}
return 0;
};
double u(double x) {
double result;
result = 1 - (1 - exp(-10))*x - exp(-10*x);
return result;
};
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;
}
}

17
src/analyticPlot.py Normal file
View File

@ -0,0 +1,17 @@
import numpy as np
import matplotlib.pyplot as plt
x = []
y = []
v = []
with open('testdata.txt') as f:
for line in f:
a, b, c = line.strip().split()
x.append(float(a))
# y.append(float(b))
v.append(float(c))
fig, ax = plt.subplots()
ax.plot(x, v)
plt.show()
# plt.savefig("main.png")

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 <cmath>
#include <ctime>
#include <fstream>
#include <iomanip>
#include <ios>
#include <string>
#include <iostream>
#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();
double f(double x) {
return 100. * std::exp(-10.*x);
}
int main()
{
timing();
general_algorithm_main();
general_algorithm_error();
special_algorithm_main();
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;
}

View File

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

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()