Implement Newton's method for a system of nonlinear equations $f(x) = 0$, where $f = (f^1,...,f^n)^T = 0$ and $x = (x^1,...,x^n)^T$. Both the function $f$ and the Jacobian $J$ are given as lambda functions. To solve the linear system requires at each iteration step, use the Gaussian elimination with partial pivoting. Test the method by finding a root of the nonlinear system. $$3x_1^2 + x_1x_2 -1 =0, x_1x_2+x_2^2 - 2 = 0$$

Looking for some help with the execution of putting this together. I was able to code Gaussian elimination and Newtons method separately but I am not sure how to put them together into one code. I have included all my workings below.

Both the function $f$ and the Jacobian $J$ are given as lambda functions.

For this part I coded the following:

def f(x):
    return np.array([3*x[1]**2+x[1]*x[2]-1, x[1]*x[2]+x[2]**2-2])
def j(x):
    return np.array([[6*x[1]+x[2], x[1]],[x[2], x[1]+2*x[2]]])

It is the part where at each iteration step use the gaussian elimination with partial pivoting is where I am getting confused. In a previous question I was able to code gaussian elimination with partial pivoting for a random $100\times 100$ matrix :

import numpy as np
def GAUSSPARTIALPIVOT(A, b, doTest = True):
    n = len(A)
    if b.size != n:
       raise ValueError("Invalid argument: incompatible    sizes between"+ "A & b.", b.size, n)
    for k in range(n-1):
       if doTest:
            maxindex = abs(A[k:,k]).argmax() + k
            if A[maxindex, k] == 0:
                raise ValueError("Matrix is singular.")
            if maxindex != k:
                A[[k,maxindex]] = A[[maxindex, k]]
                b[[k,maxindex]] = b[[maxindex, k]]
            if A[k, k] == 0:
                raise ValueError("Pivot element is zero. Try setting doPricing to True.")
        for row in range(k+1, n):
            multiply = A[row,k]/A[k,k]
            A[row, k:] = A[row, k:] - multiply*A[k, k:]
            b[row] = b[row] - multiply*b[k]
    x = np.zeros(n)
    for k in range(n-1, -1, -1):
        x[k] = (b[k] - np.dot(A[k,k+1:],x[k+1:]))/A[k,k]
    return x

if __name__ == "__main__":
    A = np.round(np.random.rand(100, 100)*10)
    b =  np.round(np.random.rand(100)*10)
    #Prints Ax=b, where A is a random 100x100 matrix and b is a random 100x1 vector
    print (GAUSSPARTIALPIVOT(np.copy(A), np.copy(b), doTest=False))

I was able to also code Newtons method for a function $f(x) = \cos x - \sin x $:

def Newton(f, dfdx, x, eps):
    f_value = f(x)
    iteration_counter = 0
    while abs(f_value) > eps and iteration_counter < 100:
            x = x - float(f_value)/dfdx(x)
        except ZeroDivisionError:
            print ("Error! - derivative zero for x = ", x)
        f_value = f(x)
        iteration_counter += 1
    if abs(f_value) > eps:
        iteration_counter = -1
    return x, iteration_counter

def f(x):
    return (math.cos(x)-math.sin(x))

def dfdx(x):
    return (-math.sin(x)-math.cos(x))

solution, no_iterations = Newton(f, dfdx, x=1, eps=1.0e-14)

if no_iterations > 0:    
    print ("Number of iterations: %d" % (no_iterations))
    print ("Answer: %f" % (solution))
    print ("Solution not found!")
  • $\begingroup$ @Bernard thanks for the edit, would you be able to help me with my problem? $\endgroup$
    – jh123
    Nov 6 '18 at 15:47
  • $\begingroup$ @moo any help is appreciated I am stuck $\endgroup$
    – jh123
    Nov 6 '18 at 15:58

Let's take a step back and look at the big picture. Newton's method says:

$x_{n+1}$ is the zero of the derivative at $x_n$

The familiar one-dimensional formula, which you implemented, is $$ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$$ and is gotten by solving the equation $0 = f'(x_n)(x_{n+1} - x_n) + f(x_n)$.

In the multivariable case, $x_n$ is a sequence of vectors, obtained by fixing an initial guess $x_0$ and then solving the system $$ df(x_n)(x_{n+1} - x_n) + f(x_n) = 0 $$ This is why you need an implementation of Gaussian elimination: instead of manually solving, as in the one-dimensional case, we're letting a computer solve for us.

So the pseudocode for the implementation is something like:

def newtonMethod(f, df, x0, num_iterations=1000):
    x_curr = x0
    x_next = gaussElimination(df(x_curr), -f(x_curr)) + x_curr
    for _ in range(num_iterations):
        x_curr = x_next
        x_next = gaussElimination(df(x_curr), -f(x_curr)) + x_curr
    return x_next

PS- if you want to check your code, you might find scipy's implementation useful: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.newton.html

  • $\begingroup$ thanks for the post, so how could I put this together with my codes above to get it working? $\endgroup$
    – jh123
    Nov 6 '18 at 16:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.