Solution of Nonlinear Equations by Newton Raphson Approach in Python

In this tutorial, we are going to apply the widely popular, well-suited Newton-Raphson approach in Python to solve a pair of nonlinear equations. Let's say, we have the following two nonlinear equations:

y - cosh(x) = 0

x^2 + y^2 - 2 = 0

We will solve them in Python by coding the Newton-Raphson algorithm. Then, we will plot the functions together in a general X-Y plots. As we see that both x and y are variables in the above equations.

The Newton-Raphson method will be written in Python in detail. You will find several applications of this method with MATLAB and Mathematica in this blog that had been discussed earlier [please click the numbers to see previous discussions 1, 2, 3, 4].

Python Codes for the Newton-Raphson Method

import math # Importing the math library

import matplotlib.pyplot as plt # Importing 'matplotlib' library for plotting results

import numpy as np # Importing 'Numpy' Library for array operations

def fun_1(x): # Defining the given function

    return math.cosh(x) 

def fun_2(x): # Defining another given function

    return math.sqrt(2-x**2) 

step= 0.005 # Tolerance/step size

x = np.arange(-1.4,1.4,step) 

# Initializing arrays for the two given functions

f1_x = np.zeros(len(x)) 

f2_x = np.zeros(len(x)) 

# Evaluating f1_x and f2_x 

count= 0 

for i in x:

    f1_x[count] = fun_1(i) 

    f2_x[count] = fun_2(i) 

    count+= 1 

# Plotting the functions     

plt.plot(x, f1_x, label="$y - \cosh(x) = 0$") 

plt.plot(x, f2_x, label="$x^2 + y^2 - 2 = 0$") 

plt.xlabel("$x$"); plt.ylabel("$y$") 

plt.legend() 

plt.title("$X-Y$ Graph") 

# Formatting the plot

plt.rcParams["figure.figsize"] = (15,5) 

plt.rcParams['lines.linewidth'] = 4 

plt.rcParams.update({'font.size': 14}) 

plt.show() 

def func_sys(x_array): 

    x = x_array[0,0].copy() 

    y = x_array[1,0].copy() 

    return np.array([[y-math.cosh(x)],[x**2 + y**2 - 2]])

# Defining the jacobian of the equations 

def jacobian(x_arr):

    x = x_array[0,0].copy() 

    y = x_array[1,0].copy() 

    return np.array([[-math.sinh(x), 1],[2*x, 2*y]])

def fun_NR(F,J,x0,tol,iter_max):

    error = 1; iter = 0 

    while abs(error) > tol: 

        iter += 1 

        F_J, er_GS, iter_GS = fun(J(x0),F(x0),x0,tol,iter_max) 

        xnew = x0 - F_J 

        error = np.linalg.norm(xnew-x0)/np.linalg.norm(xnew) 

        x0 = xnew.copy() 

        if iter > iter_max:

            return 'Maximum number of iterations done'

        return xnew, error, iter

    # First initial guess 

    x01 = np.array([[-1.],[1.]]) 

    x_NRl, er_NRl, iter_NRl = fun_NR(func_sys,jacobian,x01,0.00001,100) 

    # Second initial guess 

    x02 = np.array([[1.],[1.]]) 

    x_NR2, er_NR2, iter_NR2 = fun_NR(func_sys,jacobian,x02,0.00001,100) 

    func_sys(x_NRl)

    func_sys(x_NR2)

Results:

Once you compile and run the above codes with a Python distribution package, it will generate the following plot:

The above Python codes are written in Spyder (Python 3.7). You may use any other distributions to execute the codes.

#Python #Anaconda #Spyder #NewtonRaphson #NonlinearFunction