Make some changes

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);
+}
+

src/funcs.hpp

@@ -0,0 +1,20 @@
+#ifndef __FUNCS__
+#define __FUNCS__
+
+#include <armadillo>
+#include <cmath>
+
+#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

src/generalAlgorithm.cpp

@@ -1,37 +1,34 @@
+#include "funcs.hpp"
 #include "generalAlgorithm.hpp"
 #include <cmath>
-#include <fstream>
-#include <iomanip>
-#include <ios>
-#include <string>
 
-arma::vec* general_algorithm(
+arma::vec& general_algorithm(
-    arma::vec* sub_diag,
+    arma::vec& sub_diag,
-    arma::vec* main_diag,
+    arma::vec& main_diag,
-    arma::vec* sup_diag,
+    arma::vec& sup_diag,
-    arma::vec* g_vec
+    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);
+        d = sub_diag(i-1) / main_diag(i-1);
-        (*main_diag)(i) -= d*(*sup_diag)(i-1);
+        main_diag(i) -= d*sup_diag(i-1);
-        (*g_vec)(i) -= d*(*g_vec)(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.);
+        build_arrays(steps, sub_diag, main_diag, sup_diag, g_vec);
-        sup_diag.fill(-1.);
-        main_diag.fill(2.);
 
-        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++) {
+        for (int j=0; j < v_vec.n_elem; j++) {
-            ofile << std::setprecision(4) << std::scientific << step_size*(j+1) << ","
+            ofile << std::setprecision(PRECISION) << std::scientific << step_size*(j+1) << ","
-                << std::setprecision(4) << std::scientific << (*v_vec)(j) << std::endl;
+                << 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();
+        
+    }
+}

src/generalAlgorithm.hpp

@@ -2,14 +2,17 @@
 #define __GENERAL_ALG__
 
 #include <armadillo>
+#include <iomanip>
 
-arma::vec* general_algorithm(
+arma::vec& general_algorithm(
-    arma::vec* sub_diag,
+    arma::vec& sub_diag,
-    arma::vec* main_diag,
+    arma::vec& main_diag,
-    arma::vec* sup_diag,
+    arma::vec& sup_diag,
-    arma::vec* g_vec
+    arma::vec& g_vec
 );
 
 void general_algorithm_main();
 
+void general_algorithm_error();
+
 #endif

src/specialAlgorithm.cpp

@@ -1,26 +1,26 @@
+#include "funcs.hpp"
 #include "specialAlgorithm.hpp"
-#include <iomanip>
 
-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++) {
+        for (int j=0; j < v_vec.n_elem; j++) {
-            ofile << std::setprecision(4) << std::scientific << step_size*(j+1) << ","
+            ofile << std::setprecision(PRECISION) << std::scientific << step_size*(j+1) << ","
-                << std::setprecision(4) << std::scientific << (*v_vec)(j) << std::endl;
+                << std::setprecision(PRECISION) << std::scientific << v_vec(j) << std::endl;
         }
         ofile.close();
     }

src/specialAlgorithm.hpp

@@ -2,12 +2,13 @@
 #define __SPECIAL_ALG__
 
 #include <armadillo>
+#include <iomanip>
 
-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();