1import numpy as np
2import numpy.linalg as npla
3
4def newton(F, J, x0,
5 maxiter=100, tol=1.49e-8):
6 print(f'\n{0}: x = {x0}, ' + \
7 f'norm = {npla.norm(F(x0)):.1e}')
8 info = -1
9 for k in range(maxiter):
10 x = x0 - npla.inv(J(x0))@F(x0)
11 print(f'{k+1}: x = {x}, ' + \
12 f'norm = {npla.norm(F(x)):.1e}')
13 if (npla.norm(x - x0) < tol):
14 info = 0
15 break
16 x0 = x.copy()
17 return x, info
18
19def F(x):
20 n = x.size
21 y = np.empty(n)
22 y[0] = x[0]*x[1]**2 - x[0]**2*x[1] + 6
23 y[1] = x[0] + x[0]**2*x[1]**3 - 7
24 return y
25
26def J(x):
27 n = x.size
28 y = np.empty((n,n))
29 y[0,0] = x[1]**2 - 2*x[0]*x[1]
30 y[0,1] = 2*x[0]*x[1] - x[0]**2
31 y[1,0] = 1 + 2*x[0]*x[1]**3
32 y[1,1] = 3*x[0]**2*x[1]**2
33 return y
34
35x0 = np.array([-1.5, 1.5])
36x, info = newton(F, J, x0)