mirror of
https://github.com/kuhyx/WUT_Computer_Science.git
synced 2026-07-04 14:43:08 +02:00
kuhyx
GitHub
commit
76815542f5
27
code/eigenvalue_methods.py
Normal file
27
code/eigenvalue_methods.py
Normal file
@ -0,0 +1,27 @@
|
||||
from linear_algebra_utils import LinearAlgebraUtils
|
||||
|
||||
class EigenvalueMethods:
|
||||
@staticmethod
|
||||
def power_method(A, max_iter, tol=1e-6):
|
||||
n = len(A)
|
||||
x = [1] * n
|
||||
lambda_old = 0
|
||||
|
||||
for _ in range(max_iter):
|
||||
x = LinearAlgebraUtils.matrix_vector_multiply(A, x)
|
||||
lambda_new = LinearAlgebraUtils.vector_norm(x)
|
||||
x = LinearAlgebraUtils.vector_scalar_divide(x, lambda_new)
|
||||
if abs(lambda_new - lambda_old) < tol:
|
||||
break
|
||||
lambda_old = lambda_new
|
||||
|
||||
return lambda_new
|
||||
|
||||
@staticmethod
|
||||
def inverse_power_method(A, max_iter, tol=1e-6):
|
||||
n = len(A)
|
||||
I = [[1 if i == j else 0 for j in range(n)] for i in range(n)]
|
||||
A_inv = [LinearAlgebraUtils.gaussian_elimination(A.tolist(), I_col) for I_col in I]
|
||||
A_inv = list(map(list, zip(*A_inv)))
|
||||
|
||||
return 1 / EigenvalueMethods.power_method(A_inv, max_iter, tol)
|
||||
71
code/linear_algebra_utils.py
Normal file
71
code/linear_algebra_utils.py
Normal file
@ -0,0 +1,71 @@
|
||||
import math
|
||||
|
||||
class LinearAlgebraUtils:
|
||||
@staticmethod
|
||||
def dot_product(v1, v2):
|
||||
return sum(x*y for x, y in zip(v1, v2))
|
||||
|
||||
@staticmethod
|
||||
def matrix_vector_multiply(A, x):
|
||||
return [LinearAlgebraUtils.dot_product(row, x) for row in A]
|
||||
|
||||
@staticmethod
|
||||
def vector_norm(v):
|
||||
return sum(x*x for x in v)**0.5
|
||||
|
||||
@staticmethod
|
||||
def vector_scalar_divide(x, scalar):
|
||||
return [xi / scalar for xi in x]
|
||||
|
||||
@staticmethod
|
||||
def matrix_scalar_multiply(A, w):
|
||||
return [[w * A[i][j] for j in range(len(A[0]))] for i in range(len(A))]
|
||||
|
||||
@staticmethod
|
||||
def vector_vector_subtraction(v1, v2):
|
||||
return [x-y for x, y in zip(v1, v2)]
|
||||
|
||||
@staticmethod
|
||||
def vector_vector_addition(v1, v2):
|
||||
return [x+y for x, y in zip(v1, v2)]
|
||||
|
||||
@staticmethod
|
||||
def scalar_matrix_multiply(omega, vector):
|
||||
return [omega * x for x in vector]
|
||||
|
||||
|
||||
@staticmethod
|
||||
def matrix_norm(A):
|
||||
return math.sqrt(sum(sum(element ** 2 for element in row) for row in A))
|
||||
|
||||
@staticmethod
|
||||
def matrix_matrix_subtraction(A, B):
|
||||
return [[A[i][j] - B[i][j] for j in range(len(A[0]))] for i in range(len(A))]
|
||||
|
||||
@staticmethod
|
||||
def gaussian_elimination(A, b):
|
||||
n = len(A)
|
||||
M = [row[:] for row in A]
|
||||
|
||||
for i in range(n):
|
||||
M[i].append(b[i])
|
||||
|
||||
for k in range(n):
|
||||
if M[k][k] == 0:
|
||||
for i in range(k + 1, n):
|
||||
if M[i][k] != 0:
|
||||
M[k], M[i] = M[i], M[k]
|
||||
break
|
||||
|
||||
for i in range(k + 1, n):
|
||||
factor = M[i][k] / M[k][k]
|
||||
for j in range(k, n + 1):
|
||||
M[i][j] -= factor * M[k][j]
|
||||
|
||||
x = [0] * n
|
||||
for i in range(n - 1, -1, -1):
|
||||
x[i] = M[i][-1] / M[i][i]
|
||||
for k in range(i - 1, -1, -1):
|
||||
M[k][-1] -= M[k][i] * x[i]
|
||||
|
||||
return x
|
||||
90
code/main.py
90
code/main.py
@ -1,87 +1,5 @@
|
||||
import unittest
|
||||
import numpy as np # For testing ONLY
|
||||
from richardson import modified_richardson
|
||||
import pytest
|
||||
|
||||
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()
|
||||
if __name__ == "__main__":
|
||||
# Run pytest and exit with the appropriate status code
|
||||
pytest.main(["-v", "tests.py"])
|
||||
|
||||
@ -1,101 +0,0 @@
|
||||
import unittest
|
||||
import numpy as np # For testing ONLY
|
||||
from richardson_abstract 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}")
|
||||
richardson = modified_richardson(self.A_2x2, self.b_2x2, self.x0_2x2, self.alpha_2x2)
|
||||
result = richardson()
|
||||
# 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}")
|
||||
richardson = modified_richardson(self.A_3x3, self.b_3x3, self.x0_3x3, self.alpha_3x3)
|
||||
result = richardson()
|
||||
# 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):
|
||||
richardson = modified_richardson(self.A_2x2, self.b_2x2, self.x0_2x2, alpha=-0.1)
|
||||
with self.assertRaises(ValueError):
|
||||
richardson()
|
||||
# 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]
|
||||
richardson = modified_richardson(A, b, self.x0_2x2, self.alpha_2x2)
|
||||
with self.assertRaises(ValueError):
|
||||
richardson()
|
||||
# modified_richardson(A, b, self.x0_2x2, self.alpha_2x2)
|
||||
|
||||
def test_dimension_mismatch(self):
|
||||
b = [1, 2, 3] # Length mismatch with A_2x2
|
||||
richardson = modified_richardson(self.A_2x2, b, self.x0_2x2, self.alpha_2x2)
|
||||
with self.assertRaises(ValueError):
|
||||
richardson()
|
||||
# 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]
|
||||
richardson = modified_richardson(A, b, self.x0_2x2, self.alpha_2x2)
|
||||
result = richardson()
|
||||
# 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 = 10 #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}")
|
||||
richardson = modified_richardson(A.tolist(), b.tolist(), x0.tolist(), alpha, tol=1e-6, max_iter=500000)
|
||||
result = richardson()
|
||||
# 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()
|
||||
26
code/matrix_generator.py
Normal file
26
code/matrix_generator.py
Normal file
@ -0,0 +1,26 @@
|
||||
import numpy as np
|
||||
|
||||
class MatrixGenerator:
|
||||
@staticmethod
|
||||
def generate_spd_matrix(n: int) -> np.ndarray:
|
||||
"""
|
||||
Generates a random symmetric positive definite matrix of size n x n.
|
||||
|
||||
Parameters:
|
||||
n (int): The size of the matrix to generate.
|
||||
|
||||
Returns:
|
||||
np.ndarray: A symmetric positive definite matrix of size n x n.
|
||||
"""
|
||||
A = np.random.rand(n, n)
|
||||
spd_matrix = np.dot(A, A.T) + n * np.eye(n) # Adding n*I ensures positive definiteness
|
||||
return spd_matrix
|
||||
|
||||
@staticmethod
|
||||
def generate_random_matrix_and_vector(size):
|
||||
A = MatrixGenerator.generate_spd_matrix(size)
|
||||
b = np.random.uniform(-1, 1, size)
|
||||
return A, b
|
||||
|
||||
def generate_identity_matrix(size):
|
||||
return np.eye(size)
|
||||
@ -1,59 +0,0 @@
|
||||
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
|
||||
@ -1,128 +0,0 @@
|
||||
from abc import ABC, abstractmethod
|
||||
import threading
|
||||
|
||||
class modified_richardson_base(ABC):
|
||||
"""
|
||||
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.
|
||||
"""
|
||||
def __init__(self, A, b, x0, alpha, tol=1e-6, max_iter=1000):
|
||||
self.A = A
|
||||
self.b = b
|
||||
self.x0 = x0
|
||||
self.alpha = alpha
|
||||
self.tol = tol
|
||||
self.max_iter = max_iter
|
||||
self.n = len(A)
|
||||
self.x = self.x0[:]
|
||||
|
||||
def check_input_data(self):
|
||||
if len(self.A) != len(self.A[0]):
|
||||
raise ValueError("Matrix A must be square.")
|
||||
if len(self.b) != self.n:
|
||||
raise ValueError("Dimension mismatch between A and b.")
|
||||
if self.alpha <= 0:
|
||||
raise ValueError("Alpha must be greater than 0.")
|
||||
|
||||
@abstractmethod
|
||||
def vector_norm(self, v):
|
||||
pass
|
||||
|
||||
@abstractmethod
|
||||
def mat_vec_mult(self, mat, vec):
|
||||
pass
|
||||
|
||||
@abstractmethod
|
||||
def vec_sub(self, v1, v2):
|
||||
pass
|
||||
|
||||
@abstractmethod
|
||||
def vec_add(self, v1, v2):
|
||||
pass
|
||||
|
||||
@abstractmethod
|
||||
def vec_scalar_mult(self, scalar, vec):
|
||||
pass
|
||||
|
||||
def __call__(self):
|
||||
self.check_input_data()
|
||||
x = self.x
|
||||
|
||||
r = self.vec_sub(self.b, self.mat_vec_mult(self.A, x))
|
||||
iteration = 0
|
||||
|
||||
while self.vector_norm(r) > self.tol and iteration < self.max_iter:
|
||||
x = self.vec_add(x, self.vec_scalar_mult(self.alpha, r))
|
||||
r = self.vec_sub(self.b, self.mat_vec_mult(self.A, x))
|
||||
iteration += 1
|
||||
|
||||
if iteration == self.max_iter:
|
||||
raise ValueError("Maximum number of iterations reached before convergence")
|
||||
|
||||
return x
|
||||
|
||||
|
||||
class modified_richardson(modified_richardson_base):
|
||||
def vector_norm(self, v):
|
||||
return sum(vi ** 2 for vi in v) ** 0.5
|
||||
|
||||
def mat_vec_mult(self, mat, vec):
|
||||
return [sum(mat[i][j] * vec[j] for j in range(len(vec))) for i in range(len(mat))]
|
||||
|
||||
def vec_sub(self, v1, v2):
|
||||
return [v1[i] - v2[i] for i in range(len(v1))]
|
||||
|
||||
def vec_add(self, v1, v2):
|
||||
return [v1[i] + v2[i] for i in range(len(v1))]
|
||||
|
||||
def vec_scalar_mult(self, scalar, vec):
|
||||
return [scalar * vi for vi in vec]
|
||||
|
||||
|
||||
class modified_richardson_with_threads(modified_richardson_base):
|
||||
def compute_with_threads(self, v1, v2, function: function):
|
||||
""" Executes the provided function on many threads """
|
||||
result = [0] * len(v1)
|
||||
threads = []
|
||||
|
||||
for i in range(len(v1)):
|
||||
t = threading.Thread(target=function, args=(v1, v2, result, i))
|
||||
threads.append(t)
|
||||
t.start()
|
||||
|
||||
for t in threads:
|
||||
t.join()
|
||||
|
||||
return result
|
||||
|
||||
def subtract_elements(self, v1, v2, result, index):
|
||||
""" This function is executed by single thread. It calculates one row in matrix """
|
||||
result[index] = v1[index] - v2[index]
|
||||
|
||||
def add_elements(self, v1, v2, result, index):
|
||||
""" As above """
|
||||
result[index] = v1[index] - v2[index]
|
||||
|
||||
def vec_sub(self, v1, v2):
|
||||
return self.compute_with_threads(v1, v2, self.subtract_elements)
|
||||
|
||||
def vec_add(self, v1, v2):
|
||||
return self.compute_with_threads(v1, v2, self.add_elements)
|
||||
|
||||
43
code/richardson_method.py
Normal file
43
code/richardson_method.py
Normal file
@ -0,0 +1,43 @@
|
||||
from linear_algebra_utils import LinearAlgebraUtils
|
||||
from eigenvalue_methods import EigenvalueMethods
|
||||
from matrix_generator import MatrixGenerator
|
||||
|
||||
class RichardsonMethod:
|
||||
def __init__(self, A, b, max_iterations, size: int, x0=None, tol=1e-5):
|
||||
self.A = A
|
||||
self.b = b
|
||||
self.x0 = x0 if x0 is not None else [0.0] * len(b)
|
||||
self.max_iterations = max_iterations
|
||||
self.tol = tol
|
||||
self.I = MatrixGenerator.generate_identity_matrix(size)
|
||||
self.lambda_min, self.lambda_max = RichardsonMethod.calculate_eigenvalues(self.A, max_iterations)
|
||||
if self.lambda_min < 0:
|
||||
raise ValueError("Matrix A is not positive semi-definite.")
|
||||
self.omega = RichardsonMethod.calculate_omega(self.lambda_min, self.lambda_max)
|
||||
|
||||
@staticmethod
|
||||
def calculate_eigenvalues(A, max_iterations):
|
||||
return EigenvalueMethods.inverse_power_method(A, max_iterations), EigenvalueMethods.power_method(A, max_iterations)
|
||||
|
||||
@staticmethod
|
||||
def calculate_omega(lambda_min, lambda_max):
|
||||
return 2 / (lambda_min + lambda_max)
|
||||
|
||||
@staticmethod
|
||||
def convergence_norm(A, omega, I) -> bool:
|
||||
wA = LinearAlgebraUtils.matrix_scalar_multiply(A, omega)
|
||||
IMinuswA = LinearAlgebraUtils.matrix_matrix_subtraction(I, wA)
|
||||
norm = LinearAlgebraUtils.matrix_norm(IMinuswA)
|
||||
return norm
|
||||
|
||||
def solve(self):
|
||||
x = self.x0[:]
|
||||
if RichardsonMethod.convergence_norm(self.A, self.omega, self.I) >= 1:
|
||||
return RichardsonMethod.convergence_norm(self.A, self.omega, self.I), "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)
|
||||
x = LinearAlgebraUtils.vector_vector_addition(x, LinearAlgebraUtils.scalar_matrix_multiply(self.omega, residual))
|
||||
|
||||
return x, 0
|
||||
80
code/tests.py
Normal file
80
code/tests.py
Normal file
@ -0,0 +1,80 @@
|
||||
import pytest
|
||||
import numpy as np
|
||||
from scipy.sparse.linalg import cg
|
||||
from matrix_generator import MatrixGenerator
|
||||
from richardson_method import RichardsonMethod
|
||||
|
||||
def calculate_norm_numpy(I, w, A):
|
||||
# Calculate the difference between I and w * A
|
||||
difference = I - w * A
|
||||
|
||||
# Calculate the Euclidean norm of the difference
|
||||
norm = np.linalg.norm(difference)
|
||||
|
||||
return norm
|
||||
|
||||
def calculate_eigenvalues(A):
|
||||
# Calculate the eigenvalues of matrix A
|
||||
eigenvalues = np.linalg.eigvals(A)
|
||||
|
||||
# Find the minimum and maximum eigenvalues
|
||||
lambda_min = np.min(eigenvalues)
|
||||
lambda_max = np.max(eigenvalues)
|
||||
|
||||
return lambda_min, lambda_max
|
||||
|
||||
def calcualte_norm_from_matrix_numpy(A, n, max_iterations):
|
||||
lambda_min, lambda_max = calculate_eigenvalues(A)
|
||||
omega = 2 / (lambda_min + lambda_max)
|
||||
I = np.eye(n)
|
||||
return calculate_norm_numpy(I, omega, A)
|
||||
|
||||
|
||||
|
||||
@pytest.mark.parametrize("n", [2, 3, 4, 5, 10, 20, 50, 100])
|
||||
def test_richardson_vs_cg(n: int):
|
||||
print("matrix size: ", n)
|
||||
tolerance = 1e-5
|
||||
max_iterations=1000
|
||||
A, b = MatrixGenerator.generate_random_matrix_and_vector(n)
|
||||
richardson_solver = RichardsonMethod(A, b, max_iterations, size=n, tol=1e-7)
|
||||
solution_richardson, info_richardson = richardson_solver.solve()
|
||||
|
||||
solution_cg, info = cg(A, b)
|
||||
|
||||
if info == 0: # SciPy CG converged
|
||||
assert_scipy_converged(solution_richardson, info_richardson, solution_cg, tolerance, A, b, max_iterations, n)
|
||||
else: # SciPy CG did not converge
|
||||
assert_scipy_not_converged(solution_richardson, info_richardson, A, b)
|
||||
|
||||
def assert_scipy_converged(solution_richardson, info_richardson, solution_cg, tolerance, A, b, max_iterations, n):
|
||||
if info_richardson == "Richardson method for those values will NOT converge":
|
||||
print("Richardson did not converge, while SciPy did")
|
||||
numpy_norm = calcualte_norm_from_matrix_numpy(A, n, max_iterations)
|
||||
print("Numpy norm: ", numpy_norm, " Richardson norm: ", solution_richardson)
|
||||
assert False, "Richardson did not converge, while SciPy did"
|
||||
else:
|
||||
difference = np.linalg.norm(solution_richardson - solution_cg)
|
||||
print(f"Difference between Richardson and CG solutions: {difference:.8f}")
|
||||
if difference < tolerance:
|
||||
print("Both Richardson and CG converged and calculated correct values.")
|
||||
print("Solution CG:\n", solution_cg)
|
||||
print("Solution Richardson:\n", solution_richardson)
|
||||
else:
|
||||
print("Matrix A:\n", A)
|
||||
print("Vector b:\n", b)
|
||||
assert difference < tolerance, f"The solutions are different! Difference: {difference:.8f}"
|
||||
|
||||
def assert_scipy_not_converged(solution_richardson, info_richardson, A, b):
|
||||
if info_richardson == "Richardson method for those values will NOT converge":
|
||||
print("Richardson and SciPy did not converge")
|
||||
else:
|
||||
print("Richardson converged while SciPy did not:", solution_richardson)
|
||||
print("Matrix A:\n", A)
|
||||
print("Vector b:\n", b)
|
||||
assert False, "Richardson converged while SciPy did not"
|
||||
|
||||
if __name__ == "__main__":
|
||||
# Run pytest and exit with the appropriate status code
|
||||
for n in [2, 3, 4, 5, 10, 20, 50, 100]:
|
||||
test_richardson_vs_cg(n)
|
||||
Loading…
Reference in New Issue
Block a user