125 lines
3.3 KiB
Python
125 lines
3.3 KiB
Python
import sys
|
|
|
|
# sys.path.append("../")
|
|
from numethods import *
|
|
|
|
|
|
def demo_linear_solvers():
|
|
A = Matrix([[4, -1, 0], [-1, 4, -1], [0, -1, 3]])
|
|
b = Vector([15, 10, 10])
|
|
|
|
print("LU:", LUDecomposition(A).solve(b))
|
|
print("Gauss-Jordan:", GaussJordan(A).solve(b))
|
|
print("Cholesky:", Cholesky(A).solve(b))
|
|
print("Jacobi:", Jacobi(A, b, tol=1e-12).solve())
|
|
print("Gauss-Seidel:", GaussSeidel(A, b, tol=1e-12).solve())
|
|
|
|
|
|
def demo_qr():
|
|
A = Matrix([[2, -1], [1, 2], [1, 1]])
|
|
|
|
b = Vector([1, 2, 3])
|
|
|
|
# Factorization
|
|
qr = QRHouseholder(A)
|
|
Q, R = qr.Q, qr.R
|
|
print("Q =", Q)
|
|
print("R =", R)
|
|
|
|
qrm = QRModifiedGramSchmidt(A)
|
|
Qm, Rm = qrm.Q, qrm.R
|
|
print("Qm =", Qm)
|
|
print("Rm =", Rm)
|
|
print("Q^T Q =", Q.T @ Q)
|
|
print("Qm^T Qm =", Qm.T @ Qm)
|
|
print("A=Qm Rm =", Qm @ Rm)
|
|
print("A=Q R =", Q @ R)
|
|
|
|
# Solve Ax = b (least squares, since A is tall)
|
|
x_ls = LeastSquaresSolver(A, b).solve()
|
|
print("Least squares solution:", x_ls)
|
|
|
|
|
|
def demo_roots():
|
|
f = lambda x: x**2 - 2
|
|
df = lambda x: 2 * x
|
|
|
|
# Newton
|
|
steps = NewtonRoot(f, df, x0=1.0).trace()
|
|
print("Newton Method Trace (x^2 - 2):")
|
|
print_trace(steps)
|
|
|
|
# Secant
|
|
steps = Secant(f, 0, 2).trace()
|
|
print("\nSecant Method Trace (x^2 - 2):")
|
|
print_trace(steps)
|
|
|
|
# Bisection
|
|
steps = Bisection(f, 0, 2).trace()
|
|
print("\nBisection Method Trace (x^2 - 2):")
|
|
print_trace(steps)
|
|
|
|
# Fixed-point: solve
|
|
g = lambda x: 0.5 * (x + 2 / x)
|
|
steps = FixedPoint(g, 1.0).trace()
|
|
print("\nFixed-Point Iteration Trace (x^2 - 2):")
|
|
print_trace(steps)
|
|
|
|
|
|
def demo_interpolation():
|
|
x = [0, 1, 2, 3]
|
|
y = [1, 2, 0, 5]
|
|
newt = NewtonInterpolation(x, y)
|
|
lagr = LagrangeInterpolation(x, y)
|
|
t = 1.5
|
|
print("Newton interpolation at", t, "=", newt.evaluate(t))
|
|
print("Lagrange interpolation at", t, "=", lagr.evaluate(t))
|
|
|
|
|
|
def demo_differentiation():
|
|
f = lambda x: x**3 # f'(x) = 3x^2, f''(x) = 6x
|
|
x0 = 2.0
|
|
|
|
print("Forward :", ForwardDiff(f, x0))
|
|
print("Backward :", BackwardDiff(f, x0))
|
|
print("Central :", CentralDiff(f, x0))
|
|
print("4th order:", CentralDiff4th(f, x0))
|
|
print("Richardson:", RichardsonExtrap(f, x0))
|
|
print("Second derivative:", SecondDerivative(f, x0))
|
|
|
|
|
|
def demo_eigen():
|
|
# Not sure about eigenvalue methods but let's include some demos
|
|
|
|
A = Matrix([[4, 1, 1], [1, 3, 0], [1, 0, 2]])
|
|
print("\n=== Power Iteration ===")
|
|
solver_pi = PowerIteration(A, tol=1e-12, max_iter=100)
|
|
lam, x = solver_pi.solve()
|
|
solver_pi.trace()
|
|
print(f"Dominant eigenvalue ≈ {lam:.6f}, eigenvector ≈ {x}\n")
|
|
|
|
print("\n=== Inverse Power Iteration (shift=0) ===")
|
|
solver_ip = InversePowerIteration(A, shift=0.0, tol=1e-12, max_iter=100)
|
|
mu, x = solver_ip.solve()
|
|
solver_ip.trace()
|
|
print(f"Smallest eigenvalue ≈ {mu:.6f}, eigenvector ≈ {x}\n")
|
|
|
|
print("\n=== Rayleigh Quotient Iteration ===")
|
|
solver_rqi = RayleighQuotientIteration(A, tol=1e-12, max_iter=20)
|
|
mu, x = solver_rqi.solve()
|
|
solver_rqi.trace()
|
|
print(f"Eigenvalue ≈ {mu:.6f}, eigenvector ≈ {x}\n")
|
|
|
|
M = Matrix([[3, 1, 1], [-1, 3, 1], [1, 1, 3], [0, 2, 1]])
|
|
U, S, V = SVD(M).solve()
|
|
print("Singular values:", S)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
demo_linear_solvers()
|
|
demo_roots()
|
|
demo_interpolation()
|
|
demo_qr()
|
|
demo_differentiation()
|
|
demo_eigen()
|