Merge pull request #16 from FYS3150-G2-2023/7-solve-problem-7

7 solve problem 7

src/analyticPlot.cpp

@@ -0,0 +1,52 @@
+#include <iostream>
+#include <cmath>
+#include <vector>
+#include <string>
+#include <numeric>
+#include <fstream>
+#include <iomanip>
+
+double u(double x);
+void generate_range(std::vector<double> &vec, double start, double stop, int n);
+
+int main() {
+    int n = 1000;
+
+    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);
+
+    // Parameters for formatting 
+    int width = 12;
+    int prec = 4;
+
+    // Calculate u(x) and write to file
+    for (int i = 0; i <= x.size(); i++) {
+        y[i] = u(x[i]);
+        outfile << std::setw(width) << std::setprecision(prec) << std::scientific << x[i] 
+            << std::setw(width) << std::setprecision(prec) << std::scientific << y[i] 
+            << 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;
+    }
+}

src/analyticPlot.py

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

src/simpleFile.cpp

@@ -15,7 +15,7 @@ arma::vec* general_algorithm(
     arma::vec* g_vec
 )
 {
-    int n = main_diag->n_elem;
+    int n = g_vec->n_elem;
     double d;
 
     for (int i = 1; i < n; i++) {
@@ -40,6 +40,21 @@ arma::vec* special_algorithm(
     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;
 }