diff --git a/numethods/differentiation.py b/numethods/differentiation.py new file mode 100644 index 0000000..515e40f --- /dev/null +++ b/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 diff --git a/numethods/linalg.py b/numethods/linalg.py index 02176c9..af70b95 100644 --- a/numethods/linalg.py +++ b/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": - """Matrix multiplication with @ operator.""" - if self.n != other.m: - raise ValueError("Matrix dimensions do not align for multiplication") - return Matrix( - [ + def __matmul__(self, other: Union["Matrix", "Vector"]): + if isinstance(other, Matrix): + if self.n != other.m: + raise ValueError("dims") + 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): - """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") + def __mul__(self, s): + if isinstance(s, (int, float)): + return Matrix([[v * s for v in row] for row in self.data]) + raise TypeError("Use @ for matrix multiply; * is scalar") __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 diff --git a/numethods/orthogonal.py b/numethods/orthogonal.py index 27283df..a8992fa 100644 --- a/numethods/orthogonal.py +++ b/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) diff --git a/numethods/roots.py b/numethods/roots.py index ea2428b..734a4fc 100644 --- a/numethods/roots.py +++ b/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 + ) + )