Upload files to "numethods"

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

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

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")

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

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))