From 921d08e4e03268cf93da1656fe798e6e0ed4683a Mon Sep 17 00:00:00 2001 From: Deniz Donmez Date: Thu, 11 Sep 2025 16:38:39 +0300 Subject: [PATCH] Upload files to "numethods" --- numethods/__init__.py | 12 +++ numethods/linalg.py | 167 ++++++++++++++++++++++++++++++++++++++++ numethods/orthogonal.py | 128 ++++++++++++++++++++++++++++++ 3 files changed, 307 insertions(+) create mode 100644 numethods/__init__.py create mode 100644 numethods/linalg.py create mode 100644 numethods/orthogonal.py diff --git a/numethods/__init__.py b/numethods/__init__.py new file mode 100644 index 0000000..03532d6 --- /dev/null +++ b/numethods/__init__.py @@ -0,0 +1,12 @@ +from .linalg import Matrix, Vector +from .orthogonal import ( + QRGramSchmidt, + QRModifiedGramSchmidt, + QRHouseholder, + QRSolver, + LeastSquaresSolver, +) +from .solvers import LUDecomposition, GaussJordan, Jacobi, GaussSeidel, Cholesky +from .roots import Bisection, FixedPoint, Secant, NewtonRoot +from .interpolation import NewtonInterpolation, LagrangeInterpolation +from .exceptions import * diff --git a/numethods/linalg.py b/numethods/linalg.py new file mode 100644 index 0000000..02176c9 --- /dev/null +++ b/numethods/linalg.py @@ -0,0 +1,167 @@ +from __future__ import annotations +from typing import Iterable, Tuple, List +from .exceptions import NonSquareMatrixError, SingularMatrixError + +Number = float # We'll use float throughout + + +class Vector: + def __init__(self, data: Iterable[Number]): + self.data = [float(x) for x in data] + + def __len__(self) -> int: + return len(self.data) + + def __getitem__(self, i: int) -> Number: + return self.data[i] + + def __setitem__(self, i: int, value: Number) -> None: + self.data[i] = float(value) + + def copy(self) -> "Vector": + return Vector(self.data[:]) + + def norm_inf(self) -> Number: + return max(abs(x) for x in self.data) if self.data else 0.0 + + def __add__(self, other: "Vector") -> "Vector": + assert len(self) == len(other) + return Vector(a + b for a, b in zip(self.data, other.data)) + + def __sub__(self, other: "Vector") -> "Vector": + assert len(self) == len(other) + return Vector(a - b for a, b in zip(self.data, other.data)) + + def __mul__(self, scalar: Number) -> "Vector": + return Vector(scalar * x for x in self.data) + + __rmul__ = __mul__ + + def dot(self, other: "Vector") -> Number: + assert len(self) == len(other) + return sum(a * b for a, b in zip(self.data, other.data)) + + def __repr__(self): + return f"Vector({self.data})" + + +class Matrix: + def __init__(self, rows: List[Iterable[Number]]): + data = [list(map(float, row)) for row in rows] + if not data: + self.m, self.n = 0, 0 + else: + n = len(data[0]) + for r in data: + if len(r) != n: + raise ValueError("All rows must have the same length") + self.m, self.n = len(data), n + self.data = data + + @staticmethod + def zeros(m: int, n: int) -> "Matrix": + return Matrix([[0.0] * n for _ in range(m)]) + + @staticmethod + def identity(n: int) -> "Matrix": + A = Matrix.zeros(n, n) + for i in range(n): + A.data[i][i] = 1.0 + return A + + def copy(self) -> "Matrix": + return Matrix([row[:] for row in self.data]) + + def shape(self) -> Tuple[int, int]: + return self.m, self.n + + def __getitem__(self, idx): + i, j = idx + return self.data[i][j] + + def __setitem__(self, idx, value): + i, j = idx + self.data[i][j] = float(value) + + def row(self, i: int) -> Vector: + return Vector(self.data[i][:]) + + def col(self, j: int) -> Vector: + return Vector(self.data[i][j] for i in range(self.m)) + + def transpose(self) -> "Matrix": + return Matrix([[self.data[i][j] for i in range(self.m)] for j in range(self.n)]) + + T = property(transpose) + + def __matmul__(self, other: "Matrix") -> "Matrix": + """Matrix multiplication with @ operator.""" + if self.n != other.m: + raise ValueError("Matrix dimensions do not align for multiplication") + return Matrix( + [ + [ + sum(self.data[i][k] * other.data[k][j] for k in range(self.n)) + for j in range(other.n) + ] + for i in range(self.m) + ] + ) + + def __mul__(self, other): + """Overload * for scalar multiplication.""" + if isinstance(other, (int, float)): + return Matrix([[val * other for val in row] for row in self.data]) + raise TypeError("Use @ for matrix multiplication, * only supports scalars") + + __rmul__ = __mul__ + + def is_square(self) -> bool: + return self.m == self.n + + def augment(self, b: Vector) -> "Matrix": + if self.m != len(b): + raise ValueError("Dimension mismatch for augmentation") + return Matrix([self.data[i] + [b[i]] for i in range(self.m)]) + + def max_abs_in_col(self, col: int, start_row: int = 0) -> int: + max_i = start_row + max_val = abs(self.data[start_row][col]) + for i in range(start_row + 1, self.m): + v = abs(self.data[i][col]) + if v > max_val: + max_val, max_i = v, i + return max_i + + def swap_rows(self, i: int, j: int) -> None: + if i != j: + self.data[i], self.data[j] = self.data[j], self.data[i] + + def __repr__(self): + return f"Matrix({self.data})" + + +def forward_substitution(L: Matrix, b: Vector) -> Vector: + if not L.is_square(): + raise NonSquareMatrixError("L must be square") + n = L.n + x = [0.0] * n + for i in range(n): + s = sum(L.data[i][j] * x[j] for j in range(i)) + if abs(L.data[i][i]) < 1e-15: + raise SingularMatrixError("Zero pivot in forward substitution") + x[i] = (b[i] - s) / L.data[i][i] + return Vector(x) + + +def backward_substitution(U: Matrix, b: Vector) -> Vector: + if not U.is_square(): + raise NonSquareMatrixError("U must be square") + n = U.n + x = [0.0] * n + for i in reversed(range(n)): + s = sum(U.data[i][j] * x[j] for j in range(i + 1, n)) + if abs(U.data[i][i]) < 1e-15: + raise SingularMatrixError("Zero pivot in backward substitution") + x[i] = (b[i] - s) / U.data[i][i] + return Vector(x) diff --git a/numethods/orthogonal.py b/numethods/orthogonal.py new file mode 100644 index 0000000..24ece91 --- /dev/null +++ b/numethods/orthogonal.py @@ -0,0 +1,128 @@ +from __future__ import annotations +from typing import List +from .linalg import Matrix, Vector, backward_substitution +from .exceptions import SingularMatrixError + + +class QRGramSchmidt: + """Classical Gram–Schmidt orthogonalization.""" + + def __init__(self, A: Matrix): + self.m, self.n = A.shape() + self.Q = Matrix.zeros(self.m, self.n) + self.R = Matrix.zeros(self.n, self.n) + self._decompose(A) + + def _decompose(self, A: Matrix) -> None: + m, n = self.m, self.n + Qcols: List[Vector] = [] + for j in range(n): + v = A.col(j) + for k in range(j): + qk = Qcols[k] + r = qk.dot(v) + self.R.data[k][j] = r + v = Vector([vi - r * qi for vi, qi in zip(v.data, qk.data)]) + norm = sum(vi * vi for vi in v.data) ** 0.5 + if abs(norm) < 1e-15: + raise SingularMatrixError("Linearly dependent columns in Gram-Schmidt") + self.R.data[j][j] = norm + qj = Vector([vi / norm for vi in v.data]) + Qcols.append(qj) + for i in range(m): + self.Q.data[i][j] = qj[i] + + +class QRModifiedGramSchmidt: + """Modified Gram–Schmidt orthogonalization.""" + + def __init__(self, A: Matrix): + self.m, self.n = A.shape() + self.Q = Matrix.zeros(self.m, self.n) + self.R = Matrix.zeros(self.n, self.n) + self._decompose(A) + + def _decompose(self, A: Matrix) -> None: + m, n = self.m, self.n + V = [A.col(j).data for j in range(n)] + for i in range(n): + vi = Vector(V[i]) + norm = sum(v * v for v in vi.data) ** 0.5 + if abs(norm) < 1e-15: + raise SingularMatrixError("Linearly dependent columns in MGS") + self.R.data[i][i] = norm + qi = Vector([v / norm for v in vi.data]) + for r in range(m): + self.Q.data[r][i] = qi[r] + for j in range(i + 1, n): + r = qi.dot(Vector(V[j])) + self.R.data[i][j] = r + V[j] = [vj - r * qi_k for vj, qi_k in zip(V[j], qi.data)] + + +class QRHouseholder: + """Stable QR decomposition using Householder reflectors.""" + + def __init__(self, A: Matrix): + self.m, self.n = A.shape() + self.R = A.copy() + self.Q = Matrix.identity(self.m) + self._decompose() + + def _decompose(self) -> None: + m, n = self.m, self.n + for k in range(min(m, n)): + x = [self.R.data[i][k] for i in range(k, m)] + normx = sum(xi * xi for xi in x) ** 0.5 + if normx < 1e-15: + continue + sign = 1.0 if x[0] >= 0 else -1.0 + u1 = x[0] + sign * normx + v = [xi / u1 if i > 0 else 1.0 for i, xi in enumerate(x)] + normv = sum(vi * vi for vi in v) ** 0.5 + v = [vi / normv for vi in v] + for j in range(k, n): + s = sum(v[i] * self.R.data[k + i][j] for i in range(len(v))) + for i in range(len(v)): + self.R.data[k + i][j] -= 2 * s * v[i] + for j in range(m): + s = sum(v[i] * self.Q.data[j][k + i] for i in range(len(v))) + for i in range(len(v)): + self.Q.data[j][k + i] -= 2 * s * v[i] + self.Q = self.Q.transpose() + + +class QRSolver: + """Solve Ax=b given QR (square A).""" + + def __init__(self, qr: QRHouseholder | QRGramSchmidt | QRModifiedGramSchmidt): + self.Q, self.R = qr.Q, qr.R + + def solve(self, b: Vector) -> Vector: + Qtb = Vector( + [ + sum(self.Q.data[i][j] * b[i] for i in range(self.Q.m)) + for j in range(self.Q.n) + ] + ) + return backward_substitution(self.R, Qtb) + + +class LeastSquaresSolver: + """Solve overdetermined system Ax ≈ b in least squares sense using QR.""" + + def __init__(self, A: Matrix, b: Vector): + self.A, self.b = A, b + + def solve(self) -> Vector: + qr = QRHouseholder(self.A) + Q, R = qr.Q, qr.R + # Compute Q^T b (dimension m) + Qtb_full = [ + sum(Q.data[i][j] * self.b[i] for i in range(Q.m)) for j in range(Q.n) + ] + # Take only first n entries + Qtb = Vector(Qtb_full[: self.A.n]) + # Extract leading n×n block of R + Rtop = Matrix([R.data[i][: self.A.n] for i in range(self.A.n)]) + return backward_substitution(Rtop, Qtb)