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.title("$X-Y$ Graph"
# Formatting the plot
plt.rcParams["figure.figsize"] = (15,5
plt.rcParams['lines.linewidth'] = 4 
plt.rcParams.update({'font.size': 14}) 

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) 

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

Plot of the functions

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

No comments:

Post a Comment