# Newton's method with Gaussian elimination

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))
else:
• @Bernard thanks for the edit, would you be able to help me with my problem? Nov 6 '18 at 15:47
• @moo any help is appreciated I am stuck 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

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