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**2+x*x-1, x*x+x**2-2]) def j(x): return np.array([[6*x+x, x],[x, x+2*x]])
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!")