new module differentiation, orthogonal small fix, roots pretty table (print trace method)

introduced differentiation module
built-in l2 norm implementation
root finding problems now print #of iteration, new value, error, etc.

numethods/differentiation.py

@@ -0,0 +1,51 @@
+from __future__ import annotations
+from typing import Callable
+
+
+# ----------------------------
+# First derivative approximations
+# ----------------------------
+
+
+def ForwardDiff(f: Callable[[float], float], x: float, h: float = 1e-5) -> float:
+    """Forward finite difference approximation of f'(x)."""
+    return (f(x + h) - f(x)) / h
+
+
+def BackwardDiff(f: Callable[[float], float], x: float, h: float = 1e-5) -> float:
+    """Backward finite difference approximation of f'(x)."""
+    return (f(x) - f(x - h)) / h
+
+
+def CentralDiff(f: Callable[[float], float], x: float, h: float = 1e-5) -> float:
+    """Central finite difference approximation of f'(x) (2nd-order accurate)."""
+    return (f(x + h) - f(x - h)) / (2 * h)
+
+
+def CentralDiff4th(f: Callable[[float], float], x: float, h: float = 1e-5) -> float:
+    """Fourth-order accurate central difference approximation of f'(x)."""
+    return (-f(x + 2 * h) + 8 * f(x + h) - 8 * f(x - h) + f(x - 2 * h)) / (12 * h)
+
+
+# ----------------------------
+# Second derivative
+# ----------------------------
+
+
+def SecondDerivative(f: Callable[[float], float], x: float, h: float = 1e-5) -> float:
+    """Central difference approximation of second derivative f''(x)."""
+    return (f(x + h) - 2 * f(x) + f(x - h)) / (h**2)
+
+
+# ----------------------------
+# Richardson Extrapolation
+# ----------------------------
+
+
+def RichardsonExtrap(f: Callable[[float], float], x: float, h: float = 1e-2) -> float:
+    """Richardson extrapolation to improve derivative accuracy.
+    Combines estimates with step h and h/2.
+    """
+    D_h = CentralDiff(f, x, h)
+    D_h2 = CentralDiff(f, x, h / 2)
+    return (4 * D_h2 - D_h) / 3

numethods/linalg.py

@@ -1,5 +1,5 @@
 from __future__ import annotations
-from typing import Iterable, Tuple, List
+from typing import Iterable, Tuple, List, Union
 from .exceptions import NonSquareMatrixError, SingularMatrixError
 
 Number = float  # We'll use float throughout
@@ -24,6 +24,9 @@ class Vector:
     def norm_inf(self) -> Number:
         return max(abs(x) for x in self.data) if self.data else 0.0
 
+    def norm2(self) -> Number:
+        return sum(x * x for x in self.data) ** 0.5
+
     def __add__(self, other: "Vector") -> "Vector":
         assert len(self) == len(other)
         return Vector(a + b for a, b in zip(self.data, other.data))
@@ -94,25 +97,35 @@ class Matrix:
 
     T = property(transpose)
 
-    def __matmul__(self, other: "Matrix") -> "Matrix":
+    def __matmul__(self, other: Union["Matrix", "Vector"]):
-        """Matrix multiplication with @ operator."""
+        if isinstance(other, Matrix):
-        if self.n != other.m:
+            if self.n != other.m:
-            raise ValueError("Matrix dimensions do not align for multiplication")
+                raise ValueError("dims")
-        return Matrix(
+            return Matrix(
-            [
                 [
-                    sum(self.data[i][k] * other.data[k][j] for k in range(self.n))
+                    [
-                    for j in range(other.n)
+                        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)
                 ]
-                for i in range(self.m)
+            )
-            ]
+        elif isinstance(other, Vector):
-        )
+            if self.n != len(other):
+                raise ValueError("dims")
+            return Vector(
+                [
+                    sum(self.data[i][k] * other[k] for k in range(self.n))
+                    for i in range(self.m)
+                ]
+            )
+        else:
+            raise TypeError("Unsupported @")
 
-    def __mul__(self, other):
+    def __mul__(self, s):
-        """Overload * for scalar multiplication."""
+        if isinstance(s, (int, float)):
-        if isinstance(other, (int, float)):
+            return Matrix([[v * s for v in row] for row in self.data])
-            return Matrix([[val * other for val in row] for row in self.data])
+        raise TypeError("Use @ for matrix multiply; * is scalar")
-        raise TypeError("Use @ for matrix multiplication, * only supports scalars")
 
     __rmul__ = __mul__
 
@@ -142,6 +155,7 @@ class Matrix:
 
 
 def forward_substitution(L: Matrix, b: Vector) -> Vector:
+    """Solve Lx = b for x using forward substitution"""
     if not L.is_square():
         raise NonSquareMatrixError("L must be square")
     n = L.n
@@ -155,6 +169,7 @@ def forward_substitution(L: Matrix, b: Vector) -> Vector:
 
 
 def backward_substitution(U: Matrix, b: Vector) -> Vector:
+    """Solve Ux = b for x using backward substitution"""
     if not U.is_square():
         raise NonSquareMatrixError("U must be square")
     n = U.n

numethods/orthogonal.py

@@ -5,7 +5,7 @@ from .exceptions import SingularMatrixError
 
 
 class QRGramSchmidt:
-    """Classical Gram-Schmidt orthogonalization."""
+    """Classical Gram–Schmidt orthogonalization."""
 
     def __init__(self, A: Matrix):
         self.m, self.n = A.shape()
@@ -34,7 +34,7 @@ class QRGramSchmidt:
 
 
 class QRModifiedGramSchmidt:
-    """Modified Gram-Schmidt orthogonalization."""
+    """Modified Gram–Schmidt orthogonalization."""
 
     def __init__(self, A: Matrix):
         self.m, self.n = A.shape()
@@ -79,7 +79,7 @@ class QRHouseholder:
             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
+            normv = Vector(v).norm2()
             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)))
@@ -89,7 +89,7 @@ class QRHouseholder:
                 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()
+        # self.Q = self.Q.transpose()
 
 
 class QRSolver:
@@ -123,6 +123,6 @@ class LeastSquaresSolver:
         ]
         # Take only first n entries
         Qtb = Vector(Qtb_full[: self.A.n])
-        # Extract leading nxn block of R
+        # 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)

numethods/roots.py

@@ -1,9 +1,17 @@
 from __future__ import annotations
-from typing import Callable
+from typing import Callable, List, Dict, Any
 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):
+    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)
@@ -16,11 +24,38 @@ class Bisection:
         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)
+            c = 0.5 * (a + b)
             fc = f(c)
-            if abs(fc) <= self.tol or 0.5*(b-a) <= self.tol:
+            if abs(fc) <= self.tol or 0.5 * (b - a) <= self.tol:
                 return c
-            if fa*fc < 0:
+            if fa * fc < 0:
+                b, fb = c, fc
+            else:
+                a, fa = c, fc
+        raise ConvergenceError("Bisection did not converge")
+
+    def trace(self) -> List[Dict[str, Any]]:
+        steps = []
+        a, b, f = self.a, self.b, self.f
+        fa, fb = f(a), f(b)
+        for k in range(self.max_iter):
+            c = 0.5 * (a + b)
+            fc = f(c)
+            steps.append(
+                {
+                    "iter": k,
+                    "a": a,
+                    "b": b,
+                    "c": c,
+                    "f(a)": fa,
+                    "f(b)": fb,
+                    "f(c)": fc,
+                    "interval": b - a,
+                }
+            )
+            if abs(fc) <= self.tol or 0.5 * (b - a) <= self.tol:
+                return steps
+            if fa * fc < 0:
                 b, fb = c, fc
             else:
                 a, fa = c, fc
@@ -28,7 +63,13 @@ class Bisection:
 
 
 class FixedPoint:
-    def __init__(self, g: Callable[[float], float], x0: float, tol: float = 1e-10, max_iter: int = 10_000):
+    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:
@@ -40,28 +81,79 @@ class FixedPoint:
             x = x_new
         raise ConvergenceError("Fixed-point iteration did not converge")
 
+    def trace(self) -> List[Dict[str, Any]]:
+        steps = []
+        x = self.x0
+        for k in range(self.max_iter):
+            x_new = self.g(x)
+            steps.append({"iter": k, "x": x, "x_new": x_new, "error": abs(x_new - x)})
+            if abs(x_new - x) <= self.tol * (1.0 + abs(x_new)):
+                return steps
+            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):
+    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)
+            denom = f1 - f0
             if abs(denom) < 1e-20:
                 raise ConvergenceError("Secant encountered nearly zero denominator")
-            x2 = x1 - f1*(x1 - x0)/denom
+            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")
 
+    def trace(self) -> List[Dict[str, Any]]:
+        steps = []
+        x0, x1, f = self.x0, self.x1, self.f
+        f0, f1 = f(x0), f(x1)
+        for k 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
+            steps.append(
+                {
+                    "iter": k,
+                    "x0": x0,
+                    "x1": x1,
+                    "x2": x2,
+                    "f(x0)": f0,
+                    "f(x1)": f1,
+                    "error": abs(x2 - x1),
+                }
+            )
+            if abs(x2 - x1) <= self.tol * (1.0 + abs(x2)):
+                return steps
+            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):
+    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:
@@ -70,8 +162,50 @@ class NewtonRoot:
             dfx = self.df(x)
             if abs(dfx) < 1e-20:
                 raise ConvergenceError("Derivative near zero in Newton method")
-            x_new = x - self.f(x)/dfx
+            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")
+
+    def trace(self) -> List[Dict[str, Any]]:
+        steps = []
+        x = self.x0
+        for k 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
+            steps.append(
+                {
+                    "iter": k,
+                    "x": x,
+                    "f(x)": self.f(x),
+                    "df(x)": dfx,
+                    "x_new": x_new,
+                    "error": abs(x_new - x),
+                }
+            )
+            if abs(x_new - x) <= self.tol * (1.0 + abs(x_new)):
+                return steps
+            x = x_new
+        raise ConvergenceError("Newton method did not converge")
+
+
+def print_trace(steps: List[Dict[str, Any]]):
+    if not steps:
+        print("No steps recorded.")
+        return
+    # Get headers from dict keys
+    headers = list(steps[0].keys())
+    # Print header
+    print(" | ".join(f"{h:>10}" for h in headers))
+    print("-" * (13 * len(headers)))
+    # Print rows
+    for row in steps:
+        print(
+            " | ".join(
+                f"{row[h]:>10.6g}" if isinstance(row[h], (int, float)) else str(row[h])
+                for h in headers
+            )
+        )