From 48e54fde8d6071d4eef85c903d3dd1c36f4e1e69 Mon Sep 17 00:00:00 2001 From: Deniz Donmez Date: Thu, 11 Sep 2025 16:38:17 +0300 Subject: [PATCH] Upload files to "numethods" --- numethods/exceptions.py | 27 +++++++ numethods/interpolation.py | 57 ++++++++++++++ numethods/roots.py | 77 +++++++++++++++++++ numethods/solvers.py | 151 +++++++++++++++++++++++++++++++++++++ numethods/utils.py | 3 + 5 files changed, 315 insertions(+) create mode 100644 numethods/exceptions.py create mode 100644 numethods/interpolation.py create mode 100644 numethods/roots.py create mode 100644 numethods/solvers.py create mode 100644 numethods/utils.py diff --git a/numethods/exceptions.py b/numethods/exceptions.py new file mode 100644 index 0000000..9881ad1 --- /dev/null +++ b/numethods/exceptions.py @@ -0,0 +1,27 @@ +class NumericalError(Exception): + """Base class for numerical method errors.""" + pass + +class NonSquareMatrixError(NumericalError): + """Raised when a non-square matrix is provided where a square matrix is required.""" + pass + +class SingularMatrixError(NumericalError): + """Raised when a matrix is singular to working precision.""" + pass + +class NotSymmetricError(NumericalError): + """Raised when a matrix expected to be symmetric is not.""" + pass + +class NotPositiveDefiniteError(NumericalError): + """Raised when a matrix expected to be SPD is not.""" + pass + +class ConvergenceError(NumericalError): + """Raised when an iterative method fails to converge within limits.""" + pass + +class DomainError(NumericalError): + """Raised when inputs violate a method's domain assumptions.""" + pass diff --git a/numethods/interpolation.py b/numethods/interpolation.py new file mode 100644 index 0000000..370170b --- /dev/null +++ b/numethods/interpolation.py @@ -0,0 +1,57 @@ +from __future__ import annotations +from typing import List + +class NewtonInterpolation: + """Polynomial interpolation via divided differences.""" + def __init__(self, x: List[float], y: List[float]): + if len(x) != len(y): + raise ValueError("x and y must have same length") + if len(set(x)) != len(x): + raise ValueError("x values must be distinct") + self.x = [float(v) for v in x] + self.coeffs = self._divided_differences(self.x, [float(v) for v in y]) + + def _divided_differences(self, x: List[float], y: List[float]) -> List[float]: + n = len(x) + coef = y[:] + for j in range(1, n): + for i in range(n-1, j-1, -1): + denom = x[i] - x[i-j] + if abs(denom) < 1e-20: + raise ZeroDivisionError("Repeated x values in divided differences") + coef[i] = (coef[i] - coef[i-1]) / denom + return coef + + def evaluate(self, t: float) -> float: + n = len(self.x) + result = 0.0 + for i in reversed(range(n)): + result = result * (t - self.x[i]) + self.coeffs[i] + return result + + +class LagrangeInterpolation: + """Lagrange-form polynomial interpolation.""" + def __init__(self, x: List[float], y: List[float]): + if len(x) != len(y): + raise ValueError("x and y must have same length") + if len(set(x)) != len(x): + raise ValueError("x values must be distinct") + self.x = [float(v) for v in x] + self.y = [float(v) for v in y] + + def evaluate(self, t: float) -> float: + x, y = self.x, self.y + n = len(x) + total = 0.0 + for i in range(n): + L = 1.0 + for j in range(n): + if i == j: + continue + denom = x[i] - x[j] + if abs(denom) < 1e-20: + raise ZeroDivisionError("Repeated x values in Lagrange basis") + L *= (t - x[j]) / denom + total += y[i]*L + return total diff --git a/numethods/roots.py b/numethods/roots.py new file mode 100644 index 0000000..ea2428b --- /dev/null +++ b/numethods/roots.py @@ -0,0 +1,77 @@ +from __future__ import annotations +from typing import Callable +from .exceptions import ConvergenceError, DomainError + +class Bisection: + def __init__(self, f: Callable[[float], float], a: float, b: float, tol: float = 1e-10, max_iter: int = 10_000): + if a >= b: + raise ValueError("Require a < b") + fa, fb = f(a), f(b) + if fa * fb > 0: + raise DomainError("f(a) and f(b) must have opposite signs") + self.f, self.a, self.b = f, a, b + self.tol, self.max_iter = tol, max_iter + + def solve(self) -> float: + a, b, f = self.a, self.b, self.f + fa, fb = f(a), f(b) + for _ in range(self.max_iter): + c = 0.5*(a+b) + fc = f(c) + if abs(fc) <= self.tol or 0.5*(b-a) <= self.tol: + return c + if fa*fc < 0: + b, fb = c, fc + else: + a, fa = c, fc + raise ConvergenceError("Bisection did not converge") + + +class FixedPoint: + def __init__(self, g: Callable[[float], float], x0: float, tol: float = 1e-10, max_iter: int = 10_000): + self.g, self.x0, self.tol, self.max_iter = g, x0, tol, max_iter + + def solve(self) -> float: + x = self.x0 + for _ in range(self.max_iter): + x_new = self.g(x) + if abs(x_new - x) <= self.tol * (1.0 + abs(x_new)): + return x_new + x = x_new + raise ConvergenceError("Fixed-point iteration did not converge") + + +class Secant: + def __init__(self, f: Callable[[float], float], x0: float, x1: float, tol: float = 1e-10, max_iter: int = 10_000): + self.f, self.x0, self.x1, self.tol, self.max_iter = f, x0, x1, tol, max_iter + + def solve(self) -> float: + x0, x1, f = self.x0, self.x1, self.f + f0, f1 = f(x0), f(x1) + for _ in range(self.max_iter): + denom = (f1 - f0) + if abs(denom) < 1e-20: + raise ConvergenceError("Secant encountered nearly zero denominator") + x2 = x1 - f1*(x1 - x0)/denom + if abs(x2 - x1) <= self.tol * (1.0 + abs(x2)): + return x2 + x0, x1 = x1, x2 + f0, f1 = f1, f(x1) + raise ConvergenceError("Secant did not converge") + + +class NewtonRoot: + def __init__(self, f: Callable[[float], float], df: Callable[[float], float], x0: float, tol: float = 1e-10, max_iter: int = 10_000): + self.f, self.df, self.x0, self.tol, self.max_iter = f, df, x0, tol, max_iter + + def solve(self) -> float: + x = self.x0 + for _ in range(self.max_iter): + dfx = self.df(x) + if abs(dfx) < 1e-20: + raise ConvergenceError("Derivative near zero in Newton method") + x_new = x - self.f(x)/dfx + if abs(x_new - x) <= self.tol * (1.0 + abs(x_new)): + return x_new + x = x_new + raise ConvergenceError("Newton method did not converge") diff --git a/numethods/solvers.py b/numethods/solvers.py new file mode 100644 index 0000000..f9aa4db --- /dev/null +++ b/numethods/solvers.py @@ -0,0 +1,151 @@ +from __future__ import annotations +from .linalg import Matrix, Vector, forward_substitution, backward_substitution +from .exceptions import NonSquareMatrixError, SingularMatrixError, NotSymmetricError, NotPositiveDefiniteError, ConvergenceError + +class LUDecomposition: + """LU with partial pivoting: PA = LU""" + def __init__(self, A: Matrix): + if not A.is_square(): + raise NonSquareMatrixError("A must be square") + self.n = A.n + self.L = Matrix.identity(self.n) + self.U = A.copy() + self.P = Matrix.identity(self.n) + self._decompose() + + def _decompose(self) -> None: + n = self.n + for k in range(n): + pivot_row = self.U.max_abs_in_col(k, k) + if abs(self.U.data[pivot_row][k]) < 1e-15: + raise SingularMatrixError("Matrix is singular to working precision") + self.U.swap_rows(k, pivot_row) + self.P.swap_rows(k, pivot_row) + if k > 0: + self.L.data[k][:k], self.L.data[pivot_row][:k] = self.L.data[pivot_row][:k], self.L.data[k][:k] + for i in range(k+1, n): + m = self.U.data[i][k] / self.U.data[k][k] + self.L.data[i][k] = m + for j in range(k, n): + self.U.data[i][j] -= m * self.U.data[k][j] + + def solve(self, b: Vector) -> Vector: + Pb = Vector([sum(self.P.data[i][j]*b[j] for j in range(self.n)) for i in range(self.n)]) + y = forward_substitution(self.L, Pb) + x = backward_substitution(self.U, y) + return x + + +class GaussJordan: + def __init__(self, A: Matrix): + if not A.is_square(): + raise NonSquareMatrixError("A must be square") + self.n = A.n + self.A = A.copy() + + def solve(self, b: Vector) -> Vector: + n = self.n + Ab = self.A.augment(b) + for col in range(n): + pivot = Ab.max_abs_in_col(col, col) + if abs(Ab.data[pivot][col]) < 1e-15: + raise SingularMatrixError("Matrix is singular or nearly singular") + Ab.swap_rows(col, pivot) + pv = Ab.data[col][col] + Ab.data[col] = [v / pv for v in Ab.data[col]] + for r in range(n): + if r == col: + continue + factor = Ab.data[r][col] + Ab.data[r] = [rv - factor*cv for rv, cv in zip(Ab.data[r], Ab.data[col])] + return Vector(row[-1] for row in Ab.data) + + +class Jacobi: + def __init__(self, A: Matrix, b: Vector, tol: float = 1e-10, max_iter: int = 10_000): + if not A.is_square(): + raise NonSquareMatrixError("A must be square") + if A.n != len(b): + raise ValueError("Dimension mismatch") + self.A = A.copy() + self.b = b.copy() + self.tol = tol + self.max_iter = max_iter + + def solve(self, x0: Vector | None = None) -> Vector: + n = self.A.n + x = Vector([0.0]*n) if x0 is None else x0.copy() + for _ in range(self.max_iter): + x_new = [0.0]*n + for i in range(n): + diag = self.A.data[i][i] + if abs(diag) < 1e-15: + raise SingularMatrixError("Zero diagonal entry in Jacobi") + s = sum(self.A.data[i][j]*x[j] for j in range(n) if j != i) + x_new[i] = (self.b[i] - s) / diag + x_new = Vector(x_new) + if (x_new - x).norm_inf() <= self.tol * (1.0 + x_new.norm_inf()): + return x_new + x = x_new + raise ConvergenceError("Jacobi did not converge within max_iter") + + +class GaussSeidel: + def __init__(self, A: Matrix, b: Vector, tol: float = 1e-10, max_iter: int = 10_000): + if not A.is_square(): + raise NonSquareMatrixError("A must be square") + if A.n != len(b): + raise ValueError("Dimension mismatch") + self.A = A.copy() + self.b = b.copy() + self.tol = tol + self.max_iter = max_iter + + def solve(self, x0: Vector | None = None) -> Vector: + n = self.A.n + x = Vector([0.0]*n) if x0 is None else x0.copy() + for _ in range(self.max_iter): + x_old = x.copy() + for i in range(n): + diag = self.A.data[i][i] + if abs(diag) < 1e-15: + raise SingularMatrixError("Zero diagonal entry in Gauss-Seidel") + s1 = sum(self.A.data[i][j]*x[j] for j in range(i)) + s2 = sum(self.A.data[i][j]*x_old[j] for j in range(i+1, n)) + x[i] = (self.b[i] - s1 - s2) / diag + if (x - x_old).norm_inf() <= self.tol * (1.0 + x.norm_inf()): + return x + raise ConvergenceError("Gauss-Seidel did not converge within max_iter") + + +class Cholesky: + """A = L L^T for SPD matrices.""" + def __init__(self, A: Matrix): + if not A.is_square(): + raise NonSquareMatrixError("A must be square") + n = A.n + for i in range(n): + for j in range(i+1, n): + if abs(A.data[i][j] - A.data[j][i]) > 1e-12: + raise NotSymmetricError("Matrix is not symmetric") + self.n = n + self.L = Matrix.zeros(n, n) + self._decompose(A.copy()) + + def _decompose(self, A: Matrix) -> None: + n = self.n + for i in range(n): + for j in range(i+1): + s = sum(self.L.data[i][k]*self.L.data[j][k] for k in range(j)) + if i == j: + val = A.data[i][i] - s + if val <= 0.0: + raise NotPositiveDefiniteError("Matrix is not positive definite") + self.L.data[i][j] = val**0.5 + else: + self.L.data[i][j] = (A.data[i][j] - s) / self.L.data[j][j] + + def solve(self, b: Vector) -> Vector: + y = forward_substitution(self.L, b) + x = backward_substitution(self.L.transpose(), y) + return x diff --git a/numethods/utils.py b/numethods/utils.py new file mode 100644 index 0000000..2f4342f --- /dev/null +++ b/numethods/utils.py @@ -0,0 +1,3 @@ +from __future__ import annotations +def relative_tolerance_reached(delta: float, value: float, tol: float) -> bool: + return abs(delta) <= tol * (1.0 + abs(value))