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} $$
from numpy import array,zeros,exp
from numpy.linalg import norm,solve
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
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
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
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