Solving a non-linear system of equations in Python using Newton's Method

Question

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.

Hard to check your code, but if you have time to play, you might try plotting the surface defined by the three curves, and then you could make an animation of the solution point (x,y,z) moving in space as Newton's method iterates — kevinkayaks, Aug 25 '18 at 19:57
@kevinkayaks, thanks. Can I do something in order to make the code more readble? — Pedro Delfino, Aug 25 '18 at 19:59
In the function `jacobian_exercise(x,y,z)`, the derivative of the third equation with respect to `x` is missing some terms. — Warren Weckesser, Aug 25 '18 at 21:02

kevinkayaks · Answer 1 · 2018-08-25T23:04:49.563

To make your code more readable, I would suggest reducing the number of function definitions. They obscure the relatively simple computations which are happening.

I rewrote my own version:

def iter_newton(X,function,jacobian,imax = 1e6,tol = 1e-5):
    for i in range(int(imax)):
        J = jacobian(X) # calculate jacobian J = df(X)/dY(X) 
        Y = function(X) # calculate function Y = f(X)
        dX = np.linalg.solve(J,Y) # solve for increment from JdX = Y 
        X -= dX # step X by dX 
        if np.linalg.norm(dX)<tol: # break if converged
            print('converged.')
            break
    return X

I don't find the same behavior:

>>>X_0 = np.array([1,2,3],dtype=float)
>>>iter_newton(X_0,function_exercise,jacobian_exercise)
converged.
array([9.26836542e-18, 2.00000000e+00, 1.00000000e+00])

even works for far worse guesses

>>>X_0 = np.array([13.4,-2,31],dtype=float)
>>>iter_newton(X_0,function_exercise,jacobian_exercise)
converged.
array([1.59654153e-18, 2.00000000e+00, 1.00000000e+00])

xdze2 · Answer 2 · 2018-08-26T19:28:54.563

I tried to rewrite your code in a more Pythonic way. I hope it helps. Maybe the error is the sign of vector_b_f_output in x_delta_by_gauss(jacobian, vector_b_f_output)? or some missing term in the Jacobian.

import numpy as np

# Example from the video:
# from youtube https://www.youtube.com/watch?v=zPDp_ewoyhM
def jacobian_example(xy):
    x, y = xy
    return [[1, 2],
            [2*x, 8*y]]

def function_example(xy):
    x, y = xy
    return [x + 2*y - 2, x**2 + 4*y**2 - 4]

# From the exercise:
def function_exercise(xyz):
    x, y, z = xyz
    return [x + y + z - 3,
            x**2 + y**2 + z**2 - 5,
            np.exp(x) + x*y - x*z - 1]

def jacobian_exercise(xyz):
    x, y, z = xyz
    return [[1, 1, 1],
            [2*x, 2*y, 2*z],
            [np.exp(x) + y - z, x, -x]]



def iterative_newton(fun, x_init, jacobian):
    max_iter = 50
    epsilon = 1e-8

    x_last = x_init

    for k in range(max_iter):
        # Solve J(xn)*( xn+1 - xn ) = -F(xn):
        J = np.array(jacobian(x_last))
        F = np.array(fun(x_last))

        diff = np.linalg.solve( J, -F )
        x_last = x_last + diff

        # Stop condition:
        if np.linalg.norm(diff) < epsilon:
            print('convergence!, nre iter:', k )
            break

    else: # only if the for loop end 'naturally'
        print('not converged')

    return x_last

# For the exercice:
x_sol = iterative_newton(function_exercise, [2.0,1.0,2.0], jacobian_exercise)
print('solution exercice:', x_sol )
print('F(sol)', function_exercise(x_sol) )

# For the example:
x_sol = iterative_newton(function_example, [1.0,2.0], jacobian_example)
print('solution example:', x_sol )
print( function_example(x_sol) )

If you want to verify using fsolve:

# Verification using fsvole from Scipy
from scipy.optimize import fsolve

x0 = [2, 2, 2]
sol = fsolve(function_exercise, x0, fprime=jacobian_exercise, full_output=1)
print('solution exercice fsolve:', sol)

score 1 · Accepted Answer · answered Aug 26 '18 at 00:11

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])))))

Solving a non-linear system of equations in Python using Newton's Method

3 Answers3