From 76ace2b76c8fae65be8df4f91bc21d657eb5431c Mon Sep 17 00:00:00 2001
From: Krzysztof Rudnicki <krzysztofrudnicki0@gmail.com>
Date: Sun, 20 Oct 2024 18:09:06 +0200
Subject: [PATCH] feat: reworked convergence calculations

---
 code/Richardson_numpy.py  | 70 ---------------------------------------
 code/richardson_method.py | 28 ++++++++--------
 code/tests.py             |  5 +--
 3 files changed, 16 insertions(+), 87 deletions(-)
 delete mode 100644 code/Richardson_numpy.py

diff --git a/code/Richardson_numpy.py b/code/Richardson_numpy.py
deleted file mode 100644
index 9de5a394..00000000
--- a/code/Richardson_numpy.py
+++ /dev/null
@@ -1,70 +0,0 @@
-import numpy as np
-from scipy.sparse.linalg import cg
-
-def generate_random_matrix_and_vector(size):
-    A = np.random.uniform(-1, 1, (size, size))
-    A = np.dot(A.T, A) + np.eye(size) * size  # dodanie `size` do przekątnej zwiększa wartości własne -> obejście problemu z overflow
-    b = np.random.uniform(-1, 1, size)
-    return A, b
-
-def richardson(A, b, x0=None, max_iterations=1000, tol=1e-5):
-
-    if x0 is None:
-        x0 = np.zeros_like(b, dtype=float)
-    
-    x = np.array(x0, dtype=float)
-    n = len(b)
-
-    eigenvalues = np.linalg.eigvals(A)
-    lambda_min = np.min(eigenvalues)
-    lambda_max = np.max(eigenvalues)
-
-    if lambda_min < 0:
-        raise ValueError("Matrix A is not positive semi-definite.")
-    
-    omega = 2/(lambda_min + lambda_max)
-
-    for iteration in range(max_iterations):
-        Ax = np.dot(A, x)
-        residual = b - Ax
-
-        if np.linalg.norm(residual) < tol:
-            print('Convergence achieved.')
-            return x
-        
-        x = x + omega * residual
-
-    print('Maximum number of iterations reached without convergence.')
-    return x
-
-def run_tests():
-    test_sizes = [2, 3, 4, 5, 10, 20, 50, 100]
-    tolerance = 1e-5
-    
-    for n in test_sizes:
-        print(f"\nRunning test for n = {n}")
-        
-        A, b = generate_random_matrix_and_vector(n)
-        # print("\nMatrix A:")
-        # print(A)
-        # print("\nVector b:")
-        # print(b)
-
-        solution_richardson = richardson(A, b, max_iterations=10000000, tol=1e-7)
-        print("Richardson Method Solution:", solution_richardson)
-
-        solution_cg, info = cg(A, b)
-        if info == 0:
-            print("SciPy Conjugate Gradient solution:", solution_cg)
-        else:
-            print("SciPy Conjugate Gradient did not converge.")
-
-        difference = np.linalg.norm(solution_richardson - solution_cg)
-        print(f"Difference between Richardson and CG solutions: {difference:.8f}")
-        
-        if difference < tolerance:
-            print("The solutions are effectively the same.")
-        else:
-            print("The solutions are different!")
-
-run_tests()
diff --git a/code/richardson_method.py b/code/richardson_method.py
index 0dfc753b..59aa2964 100644
--- a/code/richardson_method.py
+++ b/code/richardson_method.py
@@ -10,27 +10,25 @@ class RichardsonMethod:
         self.max_iterations = max_iterations
         self.tol = tol
         self.I = MatrixGenerator.generate_identity_matrix(size)
+        self.lambda_min = EigenvalueMethods.inverse_power_method(self.A)
+        if self.lambda_min < 0:
+            raise ValueError("Matrix A is not positive semi-definite.")
+        self.lambda_max = EigenvalueMethods.power_method(self.A)
+        self.omega = 2 / (self.lambda_min + self.lambda_max)
+    
+    def will_converge(self) -> bool:
+        wA = LinearAlgebraUtils.matrix_scalar_multiply(self.A, self.omega)
+        IMinuswA = LinearAlgebraUtils.matrix_matrix_subtraction(self.I, wA)
+        return LinearAlgebraUtils.matrix_norm(IMinuswA) < 1
 
     def solve(self):
         x = self.x0[:]
-        lambda_min = EigenvalueMethods.inverse_power_method(self.A)
-        lambda_max = EigenvalueMethods.power_method(self.A)
-
-        if lambda_min < 0:
-            raise ValueError("Matrix A is not positive semi-definite.")
-        
-        omega = 2 / (lambda_min + lambda_max)
-
+        if not self.will_converge():
+            print("Richardson method for those values will NOT converge")
         for iteration in range(self.max_iterations):
             Ax = LinearAlgebraUtils.matrix_vector_multiply(self.A, x)
             residual = LinearAlgebraUtils.vector_vector_subtraction(self.b, Ax)
-            wA = LinearAlgebraUtils.matrix_scalar_multiply(self.A, omega)
-            IMinuswA = LinearAlgebraUtils.matrix_matrix_subtraction(self.I, wA)
-            if LinearAlgebraUtils.matrix_norm(IMinuswA) < 1:
-                print('Convergence achieved.')
-                return x
-            
-            x = LinearAlgebraUtils.vector_vector_addition(x, LinearAlgebraUtils.scalar_matrix_multiply(omega, residual))
+            x = LinearAlgebraUtils.vector_vector_addition(x, LinearAlgebraUtils.scalar_matrix_multiply(self.omega, residual))
 
         print('Maximum number of iterations reached without convergence.')
         return x
diff --git a/code/tests.py b/code/tests.py
index 0dc23ccc..80852a32 100644
--- a/code/tests.py
+++ b/code/tests.py
@@ -11,8 +11,9 @@ def run_tests():
         print(f"\nRunning test for n = {n}")
         
         A, b = MatrixGenerator.generate_random_matrix_and_vector(n)
-
-        richardson_solver = RichardsonMethod(A, b, size=n, max_iterations=10000000, tol=1e-7)
+        print("A: ", A)
+        print("b: ", b)
+        richardson_solver = RichardsonMethod(A, b, size=n, max_iterations=1000, tol=1e-7)
         solution_richardson = richardson_solver.solve()
         print("Richardson Method Solution:", solution_richardson)