From 641470e473d0e73a32dce90c1a880dd57fb96e46 Mon Sep 17 00:00:00 2001
From: crupest <crupest@outlook.com>
Date: Fri, 15 Apr 2022 22:19:23 +0800
Subject: import(life): Add Gaussian Elimination.

---
 .../numerical-analysis-lab/GaussianElimination.cpp | 120 +++++++++++++++++++++
 1 file changed, 120 insertions(+)
 create mode 100644 works/life/numerical-analysis-lab/GaussianElimination.cpp

(limited to 'works/life/numerical-analysis-lab')

diff --git a/works/life/numerical-analysis-lab/GaussianElimination.cpp b/works/life/numerical-analysis-lab/GaussianElimination.cpp
new file mode 100644
index 0000000..777eee7
--- /dev/null
+++ b/works/life/numerical-analysis-lab/GaussianElimination.cpp
@@ -0,0 +1,120 @@
+#include <cstdlib>
+#include <iomanip>
+#include <iostream>
+
+void PrintMatrix(int n, double **a) {
+  for (int i = 0; i < n; ++i) {
+    for (int j = 0; j < n; ++j) {
+      std::cout << a[i][j] << ' ';
+    }
+    std::cout << '\n';
+  }
+}
+
+void PrintEquation(int n, double **a, double *b) {
+  for (int i = 0; i < n; ++i) {
+    for (int j = 0; j < n; ++j) {
+      std::cout << a[i][j] << ' ';
+    }
+    std::cout << b[i] << '\n';
+  }
+}
+
+bool HasSolution(int n, double **a, double *b) {
+  for (int r = 0; r < n; r++) {
+    bool all_zero = true;
+    for (int c = 0; c < n; c++) {
+      if (a[r][c] != 0) {
+        all_zero = false;
+        break;
+      }
+    }
+
+    if (all_zero && b[r] != 0) {
+      return false;
+    }
+  }
+  return true;
+}
+
+void GaussianElimination(int n, double **a, double *b) {
+  for (int c = 0; c < n; c++) {
+    int pivot_r = -1;
+    for (int r = c; r < n; r++) {
+      if (a[r][c] != 0) {
+        pivot_r = r;
+        break;
+      }
+    }
+
+    if (pivot_r == -1) {
+      continue;
+    }
+
+    double pivot = a[pivot_r][c];
+
+    for (int r = 0; r < n; r++) {
+      if (r == pivot_r) {
+        continue;
+      }
+
+      double factor = a[r][c] / pivot;
+
+      for (int i = 0; i < n; i++) {
+        a[r][i] -= factor * a[pivot_r][i];
+        b[i] -= factor * b[pivot_r];
+      }
+    }
+
+    for (int i = 0; i < n; i++) {
+      a[pivot_r][i] /= pivot;
+      b[pivot_r] /= pivot;
+    }
+  }
+}
+
+double GetDiff(int n, double *x, double **a, double *b) {
+  double diff = 0;
+  for (int r = 0; r < n; r++) {
+    for (int c = 0; c < n; c++) {
+      diff += a[r][c] * x[c];
+    }
+    diff -= b[r];
+  }
+  return diff;
+}
+
+int main() {
+  std::cout << std::showpos << std::fixed << std::setprecision(3);
+
+  const int n = 6;
+  double **a = new double *[n];
+  double *b = new double[n];
+
+  for (int i = 0; i < n; i++) {
+    a[i] = new double[n];
+    for (int j = 0; j < n; j++) {
+      std::cin >> a[i][j];
+    }
+    std::cin >> b[i];
+  }
+
+  std::cout << "Equation:\n";
+  PrintEquation(n, a, b);
+
+  GaussianElimination(n, a, b);
+  std::cout << "After Gaussian Elimination:\n";
+  PrintEquation(n, a, b);
+
+  if (!HasSolution(n, a, b)) {
+    std::cout << "No solution!\n";
+  }
+
+  if (std::abs(GetDiff(n, b, a, b)) > 1e-6) {
+    std::cout << "Solution examination failed!\n";
+  } else {
+    std::cout << "Solution examination passed!\n";
+  }
+
+  return 0;
+}
-- 
cgit v1.2.3


From 2f0e9ebbbdaa4710846ee4606964c1bbf77698ed Mon Sep 17 00:00:00 2001
From: crupest <crupest@outlook.com>
Date: Fri, 15 Apr 2022 22:35:05 +0800
Subject: import(life): Add jacobian method

---
 .../numerical-analysis-lab/JacobianIteration.cpp   | 101 +++++++++++++++++++++
 1 file changed, 101 insertions(+)
 create mode 100644 works/life/numerical-analysis-lab/JacobianIteration.cpp

(limited to 'works/life/numerical-analysis-lab')

diff --git a/works/life/numerical-analysis-lab/JacobianIteration.cpp b/works/life/numerical-analysis-lab/JacobianIteration.cpp
new file mode 100644
index 0000000..20aa923
--- /dev/null
+++ b/works/life/numerical-analysis-lab/JacobianIteration.cpp
@@ -0,0 +1,101 @@
+#include <cstdlib>
+#include <ctime>
+#include <iomanip>
+#include <iostream>
+
+void PrintMatrix(int n, double **a) {
+  for (int i = 0; i < n; ++i) {
+    for (int j = 0; j < n; ++j) {
+      std::cout << a[i][j] << ' ';
+    }
+    std::cout << '\n';
+  }
+}
+
+void PrintEquation(int n, double **a, double *b) {
+  for (int i = 0; i < n; ++i) {
+    for (int j = 0; j < n; ++j) {
+      std::cout << a[i][j] << ' ';
+    }
+    std::cout << b[i] << '\n';
+  }
+}
+
+void JacobianIteration(int n, double *x, double **a, double *b) {
+  double *x_new = new double[n];
+
+  for (int r = 0; r < n; r++) {
+    x_new[r] = b[r];
+    for (int c = 0; c < n; c++) {
+      if (c != r) {
+        x_new[r] -= a[r][c] * x[c];
+      }
+    }
+    x_new[r] /= a[r][r];
+  }
+
+  for (int r = 0; r < n; r++) {
+    x[r] = x_new[r];
+  }
+
+  delete[] x_new;
+}
+
+double GetDiff(int n, double *x, double **a, double *b) {
+  double diff = 0;
+  for (int r = 0; r < n; r++) {
+    for (int c = 0; c < n; c++) {
+      diff += a[r][c] * x[c];
+    }
+    diff -= b[r];
+  }
+  return diff;
+}
+
+int main() {
+  std::cout << std::showpos << std::fixed << std::setprecision(3);
+
+  const int n = 6;
+  double **a = new double *[n];
+  double *b = new double[n];
+  double *x = new double[n];
+
+  for (int i = 0; i < n; i++) {
+    x[i] = 0;
+    a[i] = new double[n];
+    for (int j = 0; j < n; j++) {
+      std::cin >> a[i][j];
+    }
+    std::cin >> b[i];
+  }
+
+  std::cout << "Equation:\n";
+  PrintEquation(n, a, b);
+
+  int count = 0;
+
+  while (true) {
+    std::cout << "Solution:\n";
+    for (int i = 0; i < n; i++) {
+      std::cout << x[i] << '\n';
+    }
+
+    double diff = GetDiff(n, x, a, b);
+    std::cout << "diff: " << diff << '\n';
+    if (std::abs(diff) < 1e-6)
+      break;
+    count++;
+    std::cout << "Now doing iteration " << count << '\n';
+    JacobianIteration(n, x, a, b);
+
+    if (count > 100)
+      break;
+  }
+
+  std::cout << "Final Solution:\n";
+  for (int i = 0; i < n; i++) {
+    std::cout << x[i] << '\n';
+  }
+
+  return 0;
+}
-- 
cgit v1.2.3