From a15f8a1da38a1a245cdd1463f99c0f05fc86877f Mon Sep 17 00:00:00 2001 From: Cory Date: Mon, 11 Sep 2023 18:27:03 +0200 Subject: [PATCH] Make some changes --- src/funcs.cpp | 34 +++++++++++++++++ src/funcs.hpp | 20 ++++++++++ src/generalAlgorithm.cpp | 79 ++++++++++++++++++++++++++-------------- src/generalAlgorithm.hpp | 13 ++++--- src/specialAlgorithm.cpp | 26 ++++++------- src/specialAlgorithm.hpp | 5 ++- 6 files changed, 129 insertions(+), 48 deletions(-) create mode 100644 src/funcs.cpp create mode 100644 src/funcs.hpp diff --git a/src/funcs.cpp b/src/funcs.cpp new file mode 100644 index 0000000..bda3344 --- /dev/null +++ b/src/funcs.cpp @@ -0,0 +1,34 @@ +#include "funcs.hpp" + +double f(double x) { + return 100*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); +} + diff --git a/src/funcs.hpp b/src/funcs.hpp new file mode 100644 index 0000000..94053ed --- /dev/null +++ b/src/funcs.hpp @@ -0,0 +1,20 @@ +#ifndef __FUNCS__ +#define __FUNCS__ + +#include +#include + +#define PRECISION 8 + +double f(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 diff --git a/src/generalAlgorithm.cpp b/src/generalAlgorithm.cpp index bb52895..cd90384 100644 --- a/src/generalAlgorithm.cpp +++ b/src/generalAlgorithm.cpp @@ -1,37 +1,34 @@ +#include "funcs.hpp" #include "generalAlgorithm.hpp" #include -#include -#include -#include -#include -arma::vec* general_algorithm( - arma::vec* sub_diag, - arma::vec* main_diag, - arma::vec* sup_diag, - arma::vec* g_vec +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; + 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); + 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); + 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); + 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; + arma::vec sub_diag, main_diag, sup_diag, g_vec, v_vec; std::ofstream ofile; int steps; double step_size; @@ -39,23 +36,49 @@ void general_algorithm_main() for (int i = 0; i < 6; i++) { steps = std::pow(10, i+1); step_size = 1./(double) steps; - main_diag.resize(steps - 1); - g_vec.resize(steps - 1); - sub_diag.resize(steps - 2); - sup_diag.resize(steps - 2); - sub_diag.fill(-1.); - sup_diag.fill(-1.); - main_diag.fill(2.); + build_arrays(steps, sub_diag, main_diag, sup_diag, g_vec); - v_vec = general_algorithm(&sub_diag, &main_diag, &sup_diag, &g_vec); + v_vec = general_algorithm(sub_diag, main_diag, sup_diag, g_vec); - ofile.open("output/general_algorithm_" + std::to_string(steps) + ".txt"); + ofile.open("output/general/out_" + std::to_string(steps) + ".txt"); - for (int j=0; j < v_vec->n_elem; j++) { - ofile << std::setprecision(4) << std::scientific << step_size*(j+1) << "," - << std::setprecision(4) << std::scientific << (*v_vec)(j) << std::endl; + 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 < 6; 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 = f(step_size*(j+1)); + v_i = v_vec(j); + ofile << std::setprecision(PRECISION) << std::scientific + << step_size*(j+1) << "," + << std::setprecision(PRECISION) << std::scientific + << std::log10(std::abs(u_i - v_i)) << "," + << std::setprecision(PRECISION) << std::scientific + << std::log10(std::abs((u_i - v_i)/u_i)) << std::endl; + } + ofile.close(); + + } +} diff --git a/src/generalAlgorithm.hpp b/src/generalAlgorithm.hpp index cf9fd35..216e065 100644 --- a/src/generalAlgorithm.hpp +++ b/src/generalAlgorithm.hpp @@ -2,14 +2,17 @@ #define __GENERAL_ALG__ #include +#include -arma::vec* general_algorithm( - arma::vec* sub_diag, - arma::vec* main_diag, - arma::vec* sup_diag, - arma::vec* g_vec +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 diff --git a/src/specialAlgorithm.cpp b/src/specialAlgorithm.cpp index fdaf2a1..f831ff9 100644 --- a/src/specialAlgorithm.cpp +++ b/src/specialAlgorithm.cpp @@ -1,26 +1,26 @@ +#include "funcs.hpp" #include "specialAlgorithm.hpp" -#include -arma::vec* special_algorithm( +arma::vec& special_algorithm( double sub_sig, double main_sig, double sup_sig, - arma::vec* g_vec + arma::vec& g_vec ) { - int n = g_vec->n_elem; + 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); + 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); + 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); + g_vec(i) = (g_vec(i) + g_vec(i+1))/ diag(i); } return g_vec; @@ -28,7 +28,7 @@ arma::vec* special_algorithm( void special_algorithm_main() { - arma::vec g_vec, *v_vec; + arma::vec g_vec, v_vec; std::ofstream ofile; int steps; double sub_sig, main_sig, sup_sig, step_size; @@ -36,15 +36,15 @@ void special_algorithm_main() for (int i = 0; i < 6; i++) { steps = std::pow(10, i+1); step_size = 1./(double) steps; - g_vec.resize(steps - 1); + build_g_vec(steps, g_vec); - v_vec = special_algorithm(sub_sig, main_sig, sup_sig, &g_vec); + v_vec = special_algorithm(sub_sig, main_sig, sup_sig, g_vec); ofile.open("output/special_algorithm_" + std::to_string(steps) + ".txt"); - for (int j=0; j < v_vec->n_elem; j++) { - ofile << std::setprecision(4) << std::scientific << step_size*(j+1) << "," - << std::setprecision(4) << std::scientific << (*v_vec)(j) << std::endl; + 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(); } diff --git a/src/specialAlgorithm.hpp b/src/specialAlgorithm.hpp index 1578e64..68c442b 100644 --- a/src/specialAlgorithm.hpp +++ b/src/specialAlgorithm.hpp @@ -2,12 +2,13 @@ #define __SPECIAL_ALG__ #include +#include -arma::vec* special_algorithm( +arma::vec& special_algorithm( double sub_sig, double main_sig, double sup_sig, - arma::vec* g_vec + arma::vec& g_vec ); void special_algorithm_main();