new examples
This commit is contained in:
@@ -2,27 +2,6 @@ import sys
|
||||
sys.path.append('../')
|
||||
from numethods import *
|
||||
|
||||
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)
|
||||
|
||||
# Solve Ax = b (least squares, since A is tall)
|
||||
x_ls = LeastSquaresSolver(A, b).solve()
|
||||
print("Least squares solution:", x_ls)
|
||||
|
||||
|
||||
def demo_linear_solvers():
|
||||
A = Matrix([[4, -1, 0], [-1, 4, -1], [0, -1, 3]])
|
||||
@@ -34,16 +13,56 @@ def demo_linear_solvers():
|
||||
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**3 - x - 2
|
||||
df = lambda x: 3 * x**2 - 1
|
||||
print("Bisection:", Bisection(f, 1, 2).solve())
|
||||
print("Secant:", Secant(f, 1.0, 2.0).solve())
|
||||
print("Newton root:", NewtonRoot(f, df, 1.5).solve())
|
||||
# A simple contraction for demonstration; not general-purpose
|
||||
g = lambda x: (x + 2 / x**2) / 2
|
||||
print("Fixed point (demo):", FixedPoint(g, 1.5).solve())
|
||||
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():
|
||||
@@ -55,8 +74,22 @@ def demo_interpolation():
|
||||
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))
|
||||
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
demo_linear_solvers()
|
||||
demo_roots()
|
||||
demo_interpolation()
|
||||
demo_qr()
|
||||
demo_differentiation()
|
||||
|
Reference in New Issue
Block a user