'''Newton's method in two (or more) variables.

The parameter tol determines the desired accuracy, the vector x contains the starting values. 
Both have to be chosen suitably.
'''

from __future__ import print_function
import numpy as np

'''Example function in two variables.
Input:
    x ... 2D-array (current approximation to the zero of the function)
Output:
    f(x) ... column with two components, evaluated at x
    jac(x) ... the corresponding 2x2 Jacobian, evaluated at x
'''
def f(x):
   return np.array([x[0]**2 + x[1]**2 - 4, 
                    x[0]*x[1] - 1])

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

tol = 1e-10  # initialisation of tol
x = np.array([2, 1])  # starting value
print('Newton\'s method')
print('{0:.15f}'.format(x[0]), '{0:.15f}'.format(x[1]))

for k in range (0, 10):
    dx = np.linalg.solve(-jac(x),f(x))
    x = x + dx
    print('{0:.15f}'.format(x[0]), '{0:.15f}'.format(x[1]))
    error = np.linalg.norm(dx)
    print('norm of increment = ', '{0:.15f}'.format(error))
    if error < tol:
        break

raw_input()