Compare commits

..

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

28 changed files with 193 additions and 488 deletions

1
.gitignore vendored
View File

@ -43,5 +43,4 @@
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.

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

@ -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

@ -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
\begin{figure}[H]
\centering
\includegraphics[width=0.8\linewidth]{images/analytical_solution.pdf}
\caption{Plot of the analytical solution $u(x)$.}
\end{figure}
Here is the plot of the analytical solution for $u(x)$.
\includegraphics[scale=.5]{analytical_solution.pdf}

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,30 +1,17 @@
CC=g++
CCFLAGS= -std=c++11
.PHONY: clean
OBJS=generalAlgorithm.o specialAlgorithm.o funcs.o
all: simpleFile analyticPlot
EXEC=main analyticPlot
.PHONY: clean create_dirs
all: create_dirs $(EXEC)
main: main.o $(OBJS)
$(CC) $(CCFLAGS) -o $@ $^
simpleFile: simpleFile.o
$(CC) -o $@ $^
analyticPlot: analyticPlot.o
$(CC) $(CCFLAGS) -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

@ -7,7 +7,7 @@
#include <iomanip>
#define RANGE 1000
#define FILENAME "output/analytical_solution.txt"
#define FILENAME "analytical_solution.txt"
double u(double x);
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
def main():
FILENAME = "../latex/images/analytical_solution.pdf"
FILENAME = "analytical_solution.pdf"
x = []
v = []
with open('output/analytical_solution.txt') as f:
with open('analytical_solution.txt') as f:
for line in f:
a, b = line.strip().split()
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 <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;
double f(double x) {
return 100. * std::exp(-10.*x);
}
ofile.close();
double a_sol(double x) {
return 1. - (1. - std::exp(-10)) * x - std::exp(-10*x);
}
int main()
{
timing();
general_algorithm_main();
general_algorithm_error();
special_algorithm_main();
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,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()

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