feat: simple test and non-numpy sequential Richardson

code/main.py

@@ -0,0 +1,87 @@
+import unittest
+import numpy as np # For testing ONLY
+from richardson import modified_richardson
+
+class TestModifiedRichardson(unittest.TestCase):
+    def setUp(self):
+        self.A_2x2 = np.random.rand(2, 2).tolist()
+        self.b_2x2 = np.random.rand(2).tolist()
+        self.x0_2x2 = np.random.rand(2).tolist()
+        self.alpha_2x2 = 0.1
+
+        self.A_3x3 = np.random.rand(3, 3).tolist()
+        self.b_3x3 = np.random.rand(3).tolist()
+        self.x0_3x3 = np.random.rand(3).tolist()
+        self.alpha_3x3 = 0.15
+
+    def test_convergence_2x2(self):
+        print("Testing 2x2 Convergence")
+        print(f"A: {self.A_2x2}")
+        print(f"b: {self.b_2x2}")
+        print(f"x0: {self.x0_2x2}")
+        result = modified_richardson(self.A_2x2, self.b_2x2, self.x0_2x2, self.alpha_2x2)
+        expected_solution = np.linalg.solve(np.array(self.A_2x2), np.array(self.b_2x2))
+        print(f"Result: {result}")
+        print(f"Expected: {expected_solution}")
+        for r, e in zip(result, expected_solution):
+            self.assertAlmostEqual(r, e, places=4)
+
+    def test_convergence_3x3(self):
+        print("Testing 3x3 Convergence")
+        print(f"A: {self.A_3x3}")
+        print(f"b: {self.b_3x3}")
+        print(f"x0: {self.x0_3x3}")
+        result = modified_richardson(self.A_3x3, self.b_3x3, self.x0_3x3, self.alpha_3x3)
+        expected_solution = np.linalg.solve(np.array(self.A_3x3), np.array(self.b_3x3))
+        print(f"Result: {result}")
+        print(f"Expected: {expected_solution}")
+        for r, e in zip(result, expected_solution):
+            self.assertAlmostEqual(r, e, places=2)
+
+    def test_invalid_alpha(self):
+        with self.assertRaises(ValueError):
+            modified_richardson(self.A_2x2, self.b_2x2, self.x0_2x2, alpha=-0.1)
+
+    def test_non_square_matrix(self):
+        A = [[1, 2, 3], [4, 5, 6]]  # Not a square matrix
+        b = [7, 8]
+        with self.assertRaises(ValueError):
+            modified_richardson(A, b, self.x0_2x2, self.alpha_2x2)
+
+    def test_dimension_mismatch(self):
+        b = [1, 2, 3]  # Length mismatch with A_2x2
+        with self.assertRaises(ValueError):
+            modified_richardson(self.A_2x2, b, self.x0_2x2, self.alpha_2x2)
+
+    def test_zero_matrix(self):
+        A = [[0, 0], [0, 0]]
+        b = [0, 0]
+        result = modified_richardson(A, b, self.x0_2x2, self.alpha_2x2)
+        # Solution should be [0, 0]
+        print("Testing Zero Matrix")
+        print(f"A: {A}")
+        print(f"b: {b}")
+        print(f"Result: {result}")
+        self.assertEqual(result, [0, 0])
+
+    def test_large_system(self):
+        # A large test case designed to take a long time to converge
+        size = 1000
+        A = np.random.rand(size, size) + size * np.eye(size)  # Large diagonally dominant matrix
+        b = np.random.rand(size)
+        x0 = np.random.rand(size)
+        alpha = 0.01 / size  # Small alpha to ensure convergence
+
+        print("Testing Large System")
+        #print(f"A: {A}")
+        #print(f"b: {b}")
+        #print(f"x0: {x0}")
+        result = modified_richardson(A.tolist(), b.tolist(), x0.tolist(), alpha, tol=1e-6, max_iter=500000)
+        expected_solution = np.linalg.solve(A, b)
+        print(f"Result: {result}")
+        print(f"Expected: {expected_solution}")
+        for r, e in zip(result, expected_solution):
+            self.assertAlmostEqual(r, e, places=2)
+
+if __name__ == '__main__':
+    unittest.main()

code/richardson.py

@@ -0,0 +1,59 @@
+def modified_richardson(A, b, x0, alpha, tol=1e-6, max_iter=1000):
+    """
+    Solves the system of linear equations Ax = b using the Modified Richardson iteration method.
+
+    Parameters:
+    A : list of lists
+        Coefficient matrix (n x n).
+    b : list
+        Right-hand side vector (n).
+    x0 : list
+        Initial guess for the solution (n).
+    alpha : float
+        Relaxation parameter (0 < alpha < 2 / max(eigenvalue(A))).
+    tol : float, optional
+        Tolerance for the stopping criterion (default is 1e-6).
+    max_iter : int, optional
+        Maximum number of iterations (default is 1000).
+
+    Returns:
+    x : list
+        Approximate solution to the system of equations.
+    """
+    n = len(A)
+    x = x0[:]
+
+    if len(A) != len(A[0]):
+        raise ValueError("Matrix A must be square.")
+    if len(b) != n:
+        raise ValueError("Dimension mismatch between A and b.")
+    if alpha <= 0:
+        raise ValueError("Alpha must be greater than 0.")
+
+    def vector_norm(v):
+        return sum(vi ** 2 for vi in v) ** 0.5
+
+    def mat_vec_mult(mat, vec):
+        return [sum(mat[i][j] * vec[j] for j in range(len(vec))) for i in range(len(mat))]
+
+    def vec_sub(v1, v2):
+        return [v1[i] - v2[i] for i in range(len(v1))]
+
+    def vec_add(v1, v2):
+        return [v1[i] + v2[i] for i in range(len(v1))]
+
+    def vec_scalar_mult(scalar, vec):
+        return [scalar * vi for vi in vec]
+
+    r = vec_sub(b, mat_vec_mult(A, x))
+    iteration = 0
+
+    while vector_norm(r) > tol and iteration < max_iter:
+        x = vec_add(x, vec_scalar_mult(alpha, r))
+        r = vec_sub(b, mat_vec_mult(A, x))
+        iteration += 1
+
+    if iteration == max_iter:
+        raise ValueError("Maximum number of iterations reached before convergence")
+
+    return x