Numerical Methods

# ==========================================
# HARDCODED EQUATION: f(x) = x^3 - 4x - 9 = 0
# ==========================================
def f(x):
    return x**3 - 4*x - 9


def bisection(f, a, b, tol=1e-5, max_iter=100):
    if f(a) * f(b) >= 0:
        raise ValueError(f"Bisection fails: f({a}) and f({b}) must have opposite signs.")
    
    print(f"\n{'Iter':<6} {'a':<12} {'b':<12} {'c (mid)':<12} {'f(c)':<12} {'Error':<12}")
    print("-" * 66)
    
    for i in range(1, max_iter + 1):
        c = (a + b) / 2.0
        fc = f(c)
        error = abs(b - a) / 2.0
        
        print(f"{i:<6} {a:<12.6f} {b:<12.6f} {c:<12.6f} {fc:<12.6f} {error:<12.6f}")
        
        if abs(fc) < 1e-12 or error < tol:
            print("-" * 66)
            print(f"--> Root found at x = {c:.6f} in {i} iterations.\n")
            return c
        
        if f(a) * fc < 0:
            b = c
        else:
            a = c
            
    print(f"\n--> Reached maximum iterations ({max_iter}). Approximate Root = {c:.6f}\n")
    return c


if __name__ == "__main__":
    print("--- 01. Bisection Method ---")
    print("Hardcoded Equation: f(x) = x^3 - 4x - 9 = 0\n")
    
    # User inputs for interval, tolerance, and iterations
    a = float(input("Enter lower interval bound (a) [default 2.0]: ") or 2.0)
    b = float(input("Enter upper interval bound (b) [default 3.0]: ") or 3.0)
    tol_str = input("Enter tolerance (e.g. 0.0001, 1e-4, or number of decimal places like 4 or 8) [default 1e-4]: ").strip()
    if not tol_str:
        tol = 1e-4
    else:
        val = float(tol_str)
        tol = 10**(-val) if val >= 1 else val
    max_iter = int(input("Enter max iterations [default 50]: ") or 50)
    
    bisection(f, a, b, tol, max_iter)

# ==========================================
# HARDCODED EQUATION: f(x) = x^3 - 2x - 5 = 0
# ==========================================
def f(x):
    return x**3 - 2*x - 5


def regula_falsi(f, a, b, tol=1e-5, max_iter=100):
    if f(a) * f(b) >= 0:
        raise ValueError(f"Regula Falsi fails: f({a}) and f({b}) must have opposite signs.")
    
    print(f"\n{'Iter':<6} {'a':<12} {'b':<12} {'c (chord)':<12} {'f(c)':<12}")
    print("-" * 54)
    
    c_prev = a
    for i in range(1, max_iter + 1):
        fa, fb = f(a), f(b)
        # Linear interpolation false position formula
        c = (a * fb - b * fa) / (fb - fa)
        fc = f(c)
        error = abs(c - c_prev) if i > 1 else abs(b - a)
        
        print(f"{i:<6} {a:<12.6f} {b:<12.6f} {c:<12.6f} {fc:<12.6f}")
        
        if abs(fc) < 1e-12 or (i > 1 and error < tol):
            print("-" * 54)
            print(f"--> Root found at x = {c:.6f} in {i} iterations.\n")
            return c
        
        if fa * fc < 0:
            b = c
        else:
            a = c
        c_prev = c
            
    print(f"\n--> Reached maximum iterations ({max_iter}). Approximate Root = {c:.6f}\n")
    return c


if __name__ == "__main__":
    print("--- 02. Regula Falsi Method ---")
    print("Hardcoded Equation: f(x) = x^3 - 2x - 5 = 0\n")
    
    # User inputs for interval, tolerance, and iterations
    a = float(input("Enter lower interval bound (a) [default 2.0]: ") or 2.0)
    b = float(input("Enter upper interval bound (b) [default 3.0]: ") or 3.0)
    tol_str = input("Enter tolerance (e.g. 0.0001, 1e-4, or number of decimal places like 4 or 8) [default 1e-4]: ").strip()
    if not tol_str:
        tol = 1e-4
    else:
        val = float(tol_str)
        tol = 10**(-val) if val >= 1 else val
    max_iter = int(input("Enter max iterations [default 50]: ") or 50)
    
    regula_falsi(f, a, b, tol, max_iter)

# ==========================================
# HARDCODED EQUATION: f(x) = x^3 - 3x - 5 = 0
# HARDCODED DERIVATIVE: f'(x) = 3x^2 - 3
# ==========================================
def f(x):
    return x**3 - 3*x - 5

def df(x):
    return 3*x**2 - 3


def newton_raphson(f, df, a, b, tol=1e-5, max_iter=100):
    fa, fb = f(a), f(b)
    
    # 1. Verify interval [a, b] contains a root
    if fa * fb > 0:
        print(f"Warning: f({a}) = {fa:.4f} and f({b}) = {fb:.4f} have the same sign.")
        print("Root is not bracketed in this interval. Try intervals like [2, 3].\n")
    else:
        print(f"Interval [{a}, {b}] verified: f({a}) = {fa:.4f} and f({b}) = {fb:.4f} have opposite signs.")
    
    # 2. Initial guess x0 selected as midpoint of interval [a, b]
    x = (a + b) / 2.0
    print(f"Selected initial guess x0 = (a + b)/2 = {x:.4f}\n")
    
    print(f"{'Iter':<6} {'x_k':<12} {'f(x_k)':<14} {'f\'(x_k)':<14} {'x_{k+1}':<12} {'Error':<12}")
    print("-" * 72)
    
    for i in range(1, max_iter + 1):
        fx = f(x)
        dfx = df(x)
        
        # Guard against horizontal tangent
        if abs(dfx) < 1e-12:
            print("-" * 72)
            print(f"Error: Derivative f'({x}) ≈ 0 (Horizontal Tangent).")
            print("--> Please re-run and choose a different interval [a, b].\n")
            return None
            
        x_next = x - (fx / dfx)
        error = abs(x_next - x)
        
        print(f"{i:<6} {x:<12.6f} {fx:<14.6f} {dfx:<14.6f} {x_next:<12.6f} {error:<12.6f}")
        
        if abs(fx) < 1e-12 or error < tol:
            print("-" * 72)
            print(f"--> Root found at x = {x_next:.6f} in {i} iterations.\n")
            return x_next
        
        x = x_next
        
    print(f"\n--> Reached maximum iterations ({max_iter}). Approximate Root = {x:.6f}\n")
    return x


if __name__ == "__main__":
    print("--- 03. Newton-Raphson Method ---")
    print("Hardcoded Equation: f(x) = x^3 - 3x - 5 = 0\n")
    
    # Step 1: User enters the interval [a, b]
    a = float(input("Enter lower interval bound (a) [default 2.0]: ") or 2.0)
    b = float(input("Enter upper interval bound (b) [default 3.0]: ") or 3.0)
    
    # Step 2: User enters tolerance and max iterations
    tol_str = input("Enter tolerance (e.g. 0.0001, 1e-4, or number of decimal places like 4 or 8) [default 1e-4]: ").strip()
    if not tol_str:
        tol = 1e-4
    else:
        val = float(tol_str)
        tol = 10**(-val) if val >= 1 else val
        
    max_iter = int(input("Enter max iterations [default 50]: ") or 50)
    
    print(f"\nRunning with Tolerance = {tol:.1e}, Max Iterations = {max_iter}")
    newton_raphson(f, df, a, b, tol, max_iter)

# ==========================================================
# HARDCODED POLYNOMIAL: P(x) = x^4 - 3x^3 + 3x^2 - 3x + 2 = 0
# Coefficients: [a0, a1, a2, a3, a4] = [1, -3, 3, -3, 2]
# ==========================================================
POLYNOMIAL_COEFFS = [1.0, -3.0, 3.0, -3.0, 2.0]


def birge_vieta(coeffs, x0, tol=1e-5, max_iter=100):
    n = len(coeffs) - 1
    x = x0
    
    print(f"\n{'Iter':<6} {'x_k':<12} {'P(x_k)=b_n':<14} {'P\'(x_k)=c_{n-1}':<16} {'x_{k+1}':<12}")
    print("-" * 62)
    
    for iteration in range(1, max_iter + 1):
        # 1st synthetic division for P(x) -> b
        b = [0.0] * (n + 1)
        b[0] = coeffs[0]
        for i in range(1, n + 1):
            b[i] = coeffs[i] + x * b[i - 1]
            
        # 2nd synthetic division for P'(x) -> c
        c = [0.0] * n
        c[0] = b[0]
        for i in range(1, n):
            c[i] = b[i] + x * c[i - 1]
            
        P_x = b[n]
        P_prime_x = c[n - 1]
        
        if abs(P_prime_x) < 1e-12:
            raise ZeroDivisionError("Synthetic derivative P'(x) ≈ 0.")
            
        x_next = x - (P_x / P_prime_x)
        error = abs(x_next - x)
        
        print(f"{iteration:<6} {x:<12.6f} {P_x:<14.6f} {P_prime_x:<16.6f} {x_next:<12.6f}")
        
        if abs(P_x) < 1e-12 or error < tol:
            print("-" * 62)
            print(f"--> Polynomial Root found at x = {x_next:.6f} in {iteration} iterations.")
            print(f"--> Deflated Polynomial coefficients: {b[:-1]}\n")
            return x_next, b[:-1]
            
        x = x_next
        
    print(f"\n--> Reached maximum iterations ({max_iter}). Approximate Root = {x:.6f}\n")
    return x, b[:-1]


if __name__ == "__main__":
    print("--- 04. Virge Beta (Birge-Vieta) Method ---")
    print("Hardcoded Polynomial: P(x) = x^4 - 3x^3 + 3x^2 - 3x + 2 = 0\n")
    
    # User inputs for initial guess, tolerance, and iterations
    x0 = float(input("Enter initial guess (x0) [default 0.5]: ") or 0.5)
    tol = float(input("Enter tolerance (e.g. 0.0001) [default 1e-4]: ") or 1e-4)
    max_iter = int(input("Enter max iterations [default 50]: ") or 50)
    
    birge_vieta(POLYNOMIAL_COEFFS, x0, tol, max_iter)

import numpy as np

# ==========================================
# HARDCODED LINEAR SYSTEM: Ax = b
# 10x1 +   x2 +   x3 = 12
#  2x1 + 10x2 +   x3 = 13
#  2x1 +  2x2 + 10x3 = 14
# ==========================================
HARDCODED_A = [[10, 1, 1],
               [2, 10, 1],
               [2, 2, 10]]
HARDCODED_B = [12, 13, 14]


def gauss_seidel(A, b, x0=None, tol=1e-5, max_iter=100):
    A = np.array(A, dtype=float)
    b = np.array(b, dtype=float)
    n = len(b)
    
    if x0 is None:
        x = np.zeros(n)
    else:
        x = np.array(x0, dtype=float)
        
    # Check Strict Diagonal Dominance condition
    is_dominant = all(abs(A[i, i]) > sum(abs(A[i, j]) for j in range(n) if j != i) for i in range(n))
    print(f"Strict Diagonal Dominance: {'Yes (Convergence Guaranteed)' if is_dominant else 'No (May diverge)'}\n")
    
    header = f"{'Iter':<6} " + " ".join([f"x_{j+1:<10}" for j in range(n)]) + " Max Error"
    print(header)
    print("-" * len(header))
    
    for iteration in range(1, max_iter + 1):
        x_new = np.copy(x)
        
        for i in range(n):
            sum_val = b[i]
            for j in range(n):
                if j != i:
                    sum_val -= A[i, j] * x_new[j]  # Uses newly computed components immediately
            x_new[i] = sum_val / A[i, i]
            
        error = np.max(np.abs(x_new - x))
        vals = " ".join([f"{val:<12.6f}" for val in x_new])
        print(f"{iteration:<6} {vals} {error:.6e}")
        
        if error < tol:
            print("-" * len(header))
            print(f"--> Converged after {iteration} iterations.\n")
            return x_new
            
        x = x_new
        
    print(f"\n--> Reached maximum iterations ({max_iter}). Approximate Vector = {x}\n")
    return x


if __name__ == "__main__":
    print("--- 07. Gauss Sedal (Gauss-Seidel Method) ---")
    print("Hardcoded System Ax = b:")
    print("A =", HARDCODED_A)
    print("b =", HARDCODED_B, "\n")
    
    # User inputs for initial guess, tolerance, and iterations
    raw_x0 = input("Enter initial guess vector separated by commas [default 0, 0, 0]: ") or "0, 0, 0"
    x0 = [float(x.strip()) for x in raw_x0.split(",")]
    tol = float(input("Enter tolerance (e.g. 0.0001) [default 1e-4]: ") or 1e-4)
    max_iter = int(input("Enter max iterations [default 50]: ") or 50)
    
    sol = gauss_seidel(HARDCODED_A, HARDCODED_B, x0=x0, tol=tol, max_iter=max_iter)
    print("--> Final Solution Vector x =", sol, "\n")