Lab 2: Solution of Systems of Non-linear Equations using Newton-Raphson Method

import numpy as np

def f(x):
    return np.array([x[0]**2 - x[1]**2 - 5,
                    x[0]**2 + x[1]**2 - 25])

def j(x):
    return np.array([[2*x[0], -2*x[1]],
                   [2*x[0], 2*x[1]])

x = np.array([4.0, 3.0])  # Initial guess
tol = 50e-5  # Tolerance
count = 0
max_iter = 100  # Maximum iterations

while True:
    count += 1
    if count > max_iter:
        raise Exception("Maximum iterations reached without convergence")
    
    # Calculate Jacobian and function value
    A = j(x)
    B = -f(x)
    
    # Solve the linear system
    try:
        H = np.linalg.solve(A, B)
    except np.linalg.LinAlgError:
        print("Singular Jacobian matrix encountered")
        break
        
    # Update solution
    x += H
    
    # Calculate error
    err = max(abs(H))
    if err < tol:
        break

print(f"Solution found in {count} iterations:")
print(f"x1 = {x[0]:.6f}, x2 = {x[1]:.6f}")