Newton's method with Gaussian elimination

Question

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]]
        else:
            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:
        try:
            x = x - float(f_value)/dfdx(x)
        except ZeroDivisionError:
            print ("Error! - derivative zero for x = ", x)
            sys.exit(1)     
        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))
else:
    print ("Solution not found!")

Answer 1

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