Newton method for system

For $x = (x_0,x_1,x_2) \in R^3$, define $f : R^3 \to R^3$ by $$ f(x) = \begin{bmatrix} x_0 x_1 - x_2^2 - 1 \\ x_0 x_1 x_2 - x_0^2 + x_1^2 - 2 \\ \exp(x_0) - \exp(x_1) + x_2 - 3 \end{bmatrix} $$ The gradient of $f$ is given by $$ f'(x) = \begin{bmatrix} x_1 & x_0 & -2 x_2 \\ x_1 x_2 - 2 x_0 & x_0 x_2 + 2 x_1 & x_0 x_1 \\ \exp(x_0) & -\exp(x_1) & 1 \end{bmatrix} $$

In [17]:
from numpy import array,zeros,exp
from numpy.linalg import norm,solve
In [34]:
def f(x):
    y = zeros(3)
    y[0] = x[0]*x[1] - x[2]**2 - 1.0
    y[1] = x[0]*x[1]*x[2] - x[0]**2 + x[1]**2 - 2.0
    y[2] = exp(x[0]) - exp(x[1]) + x[2] - 3.0
    return y
In [35]:
def df(x):
    y = array([[x[1],               x[0],               -2.0*x[2]], \
               [x[1]*x[2]-2.0*x[0], x[0]*x[2]+2.0*x[1], x[0]*x[1]], \
               [exp(x[0]),          -exp(x[1]),        1.0]])
    return y
In [36]:
def newton(fun,dfun,x0,M=100,eps=1.0e-14,debug=False):
    x = x0
    for i in range(M):
        g = fun(x)
        J = dfun(x)
        h = solve(J,-g)
        x = x + h
        if debug:
            print i,x,norm(g)
        if norm(h) < eps * norm(x):
            return x
In [37]:
x0 = array([1.0,1.0,1.0])
x = newton(f,df,x0,debug=True)
0 [ 2.1893261   1.59847516  1.39390063] 2.44948974278
1 [ 1.85058965  1.44425142  1.278224  ] 2.52438353632
2 [ 1.7801612   1.42443598  1.23929244] 0.412336169632
3 [ 1.77767471  1.42396093  1.23747382] 0.0148129839519
4 [ 1.77767192  1.4239606   1.23747112] 1.83106597611e-05
5 [ 1.77767192  1.4239606   1.23747112] 2.43794280754e-11
6 [ 1.77767192  1.4239606   1.23747112] 6.66133814775e-16