Newton-Raphson Algorithm

newtonraphson.py

import numpy as np

def function(xy):
    x, y = xy
    return [2*x + y**2 - 8,
            x**2 - y**2 + x*y - 3]

def jacobian(xy):
    x, y = xy
    return [[2, 2*x],
            [2*x + y, x - 2*y]]


def iterative_newton(fun, x_init, jacobian, epsilon, max_iter):

    x_last = x_init

    for k in range(max_iter):
        # Solve J(xn)*( xn+1 - xn ) = -F(xn):
        J = np.array(jacobian(x_last))
        F = np.array(fun(x_last))

        diff = np.linalg.solve(J, -F)
        x_last = x_last + diff

        # Stop condition:
        if np.linalg.norm(diff) < epsilon:
            print('Convergence!, nre iter:', k )
            break

    else: # If loop ends without convergence
        print('No convergence')

    return x_last

xresult = iterative_newton(function, [0, 1], jacobian, 1e-4, 50)
print(xresult)