From c865b8b1a053f1ae91a16a0e9997fd7beae865ff Mon Sep 17 00:00:00 2001 From: Krzysztof Rudnicki Date: Fri, 18 Oct 2024 12:33:00 +0200 Subject: [PATCH] feat: simple test and non-numpy sequential Richardson --- code/main.py | 87 ++++++++++++++++++++++++++++++++++++++++++++++ code/richardson.py | 59 +++++++++++++++++++++++++++++++ 2 files changed, 146 insertions(+) create mode 100644 code/main.py create mode 100644 code/richardson.py diff --git a/code/main.py b/code/main.py new file mode 100644 index 00000000..ee207c61 --- /dev/null +++ b/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() diff --git a/code/richardson.py b/code/richardson.py new file mode 100644 index 00000000..a94b48e2 --- /dev/null +++ b/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 \ No newline at end of file