I am trying to solve this exercise for College. I have already submitted the code bellow. However, I am not completely satisfied with it.
The task is to build an implementation of Newton's method to solve the following non-linear system of equations:
In order to learn the Newton's method, besides the classes, I watched this YouTube video: https://www.youtube.com/watch?v=zPDp_ewoyhM
The guy on the video explained the math process behind Newton's method and did, manually, two iterations.
I did a Python implementation for that and the code went fine for the example on the video. Nonetheless, the example on the video deals with 2 variables and my homework deals with 3 variables. Hence, I adapted it.
That's the code:
import numpy as np
#### example from youtube https://www.youtube.com/watch?v=zPDp_ewoyhM
def jacobian_example(x,y):
return [[1,2],[2*x,8*y]]
def function_example(x,y):
return [(-1)*(x+(2*y)-2),(-1)*((x**2)+(4*(y**2))-4)]
####################################################################
### agora com os dados do exercício
def jacobian_exercise(x,y,z):
return [[1,1,1],[2*x,2*y,2*z],[np.exp(x),x,-x]]
#print (jacobian_exercise(1,2,3))
jotinha = (jacobian_exercise(1,2,3))
def function_exercise(x,y,z):
return [x+y+z-3, (x**2)+(y**2)+(z**2)-5,(np.exp(x))+(x*y)-(x*z)-1]
#print (function_exercise(1,2,3))
bezao = (function_exercise(1,2,3))
def x_delta_by_gauss(J,b):
return np.linalg.solve(J,b)
print (x_delta_by_gauss(jotinha, bezao))
x_delta_test = x_delta_by_gauss(jotinha,bezao)
def x_plus_1(x_delta,x_previous):
x_next = x_previous + x_delta
return x_next
print (x_plus_1(x_delta_test,[1,2,3]))
def newton_method(x_init):
first = x_init[0]
second = x_init[1]
third = x_init[2]
jacobian = jacobian_exercise(first, second, third)
vector_b_f_output = function_exercise(first, second, third)
x_delta = x_delta_by_gauss(jacobian, vector_b_f_output)
x_plus_1 = x_delta + x_init
return x_plus_1
def iterative_newton(x_init):
counter = 0
x_old = x_init
print ("x_old", x_old)
x_new = newton_method(x_old)
print ("x_new", x_new)
diff = np.linalg.norm(x_old-x_new)
print (diff)
while diff>0.000000000000000000000000000000000001:
counter += 1
print ("x_old", x_old)
x_new = newton_method(x_old)
print ("x_new", x_new)
diff = np.linalg.norm(x_old-x_new)
print (diff)
x_old = x_new
convergent_val = x_new
print (counter)
return convergent_val
#print (iterative_newton([1,2]))
print (iterative_newton([0,1,2]))
I am pretty sure this code is definitely not totally wrong. If I input the initial values as a vector [0,1,2], my code returns as an output [0,1,2]. This is a correct answer, it solves the three equations above.
Moreover, if a input [0,2,1], a slightly different input, the code also works and the answer it returns is also a correct one.
However, if I change my initial value to something like [1,2,3] I get a weird result: 527.7482, -1.63 and 2.14.
This result does not make any sense. Look at the first equation, if you input these values, you can easily see that (527)+(-1.63)+(2.14) does not equal to 3. This is false.
If I change the input value close to a correct solution, like [0.1,1.1,2.1] it also crashes.
OK, Newton's method does not guarantee the correct convergence. I know. It depends on the initial value, among other stuff.
Is my implementation wrong in any way? Or is the vector [1,2,3] just a "bad" initial value?
Thanks.
The guys that answered this question helped me. However, modifying one line of code made everything work in my implementation.
Since I am using the approach described on the YouTube video that I mentioned, I need to multiply the Vector-valued function by (-1), which modifies the value of each element of the vector.
I did this for the function_example
. However, when I coded function_exercise
, the one that I needed to solve for my homework without the negative sign. I missed it.
Now, it is fixed and it works fully, even with very diverse starting vectors.
import numpy as np
#### example from youtube https://www.youtube.com/watch?v=zPDp_ewoyhM
def jacobian_example(x,y):
return [[1,2],[2*x,8*y]]
def function_example(x,y):
return [(-1)*(x+(2*y)-2),(-1)*((x**2)+(4*(y**2))-4)]
####################################################################
### agora com os dados do exercício
def jacobian_exercise(x,y,z):
return [[1,1,1],[2*x,2*y,2*z],[np.exp(x),x,-x]]
#print (jacobian_exercise(1,2,3))
jotinha = (jacobian_exercise(1,2,3))
def function_exercise(x,y,z):
return [(-1)*(x+y+z-3),(-1)*((x**2)+(y**2)+(z**2)-5),(-1)*((np.exp(x))+(x*y)-(x*z)-1)]
#print (function_exercise(1,2,3))
bezao = (function_exercise(1,2,3))
def x_delta_by_gauss(J,b):
return np.linalg.solve(J,b)
print (x_delta_by_gauss(jotinha, bezao))
x_delta_test = x_delta_by_gauss(jotinha,bezao)
def x_plus_1(x_delta,x_previous):
x_next = x_previous + x_delta
return x_next
print (x_plus_1(x_delta_test,[1,2,3]))
def newton_method(x_init):
first = x_init[0]
second = x_init[1]
third = x_init[2]
jacobian = jacobian_exercise(first, second, third)
vector_b_f_output = function_exercise(first, second, third)
x_delta = x_delta_by_gauss(jacobian, vector_b_f_output)
x_plus_1 = x_delta + x_init
return x_plus_1
def iterative_newton(x_init):
counter = 0
x_old = x_init
#print ("x_old", x_old)
x_new = newton_method(x_old)
#print ("x_new", x_new)
diff = np.linalg.norm(x_old-x_new)
#print (diff)
while diff>0.0000000000001:
counter += 1
#print ("x_old", x_old)
x_new = newton_method(x_old)
#print ("x_new", x_new)
diff = np.linalg.norm(x_old-x_new)
#print (diff)
x_old = x_new
convergent_val = x_new
#print (counter)
return convergent_val
#print (iterative_newton([1,2]))
print (list(map(float,(iterative_newton([100,200,3])))))