Upload files to "numethods"

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 *

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)

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)