Welcome to the World of Modelling and Simulation

What is Modelling?

This blog is all about system dynamics modelling and simulation applied in the engineering field, especially mechanical, electrical, and ...

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