import numpy as np
import numpy.linalg as la
def f(xvec):
x, y = xvec
return np.array([
x + 2*y -2,
x**2 + 4*y**2 - 4
])
def Jf(xvec):
x, y = xvec
return np.array([
[1, 2],
[2*x, 8*y]
])
Pick an initial guess.
x = np.array([1, 2])
Now implement Newton's method.
x = x - la.solve(Jf(x), f(x))
print(x)
[ 1.50295992e-16 1.00000000e+00]
Check if that's really a solution:
f(x)
array([ 0., 0.])
Let's keep an error history and check.
xtrue = np.array([0, 1])
errors = []
x = np.array([1, 2])
x = x - la.solve(Jf(x), f(x))
errors.append(la.norm(x-xtrue))
print(x)
[ 1.50295992e-16 1.00000000e+00]
for e in errors:
print(e)
0.931694990625 0.211748861506 0.0168589857887 0.000125221235922 7.01168369152e-09 1.50295991741e-16 1.50295991741e-16 1.50295991741e-16 1.50295991741e-16
for i in range(len(errors)-1):
print(errors[i+1]/errors[i]**2)
0.243934688455 0.376001239529 0.440570178174 0.447163497456 3.05705157878 6.65353738591e+15 6.65353738591e+15 6.65353738591e+15