Newton-Raphson method in Python

AIM:

To write a python script to solve for the minimum value of pressure using the Newton-Raphson method.

OBJECTIVE:

1. To use the Newton-Raphson method to find out the value of pressure for h = 0.6 .
2. To find the optimal relaxation factor for this problem with the help of a suitable plot.
3. To tabulate the results of p for h = [0.6,1.2,1.8,2.4,3,3.6,4.2] assuming a suitable relaxation factor.

THEORY:

Newton–Raphson method:

Newton–Raphson method, named after Isaac Newton and Joseph Raphson, is a root-finding algorithm which produces successively better approximations to the roots (or zeroes) of a real-valued function. The most basic version starts with a single-variable function f defined for a real variable x, the function's derivative f ′, and an initial guess x₀ for a root of f. If the function satisfies sufficient assumptions and the initial guess is close, then

`x_1= x_0-f(x_0 )/(f^' (x_0 ) )`

is a better approximation of the root than x₀. Geometrically, (x₁, 0) is the intersection of the x-axis and the tangent of the graph of f at (x₀, f (x₀)): that is, the improved guess is the unique root of the linear approximation at the initial point. The process is repeated as

`x_(i+1)= x_i-f(x_i )/(f^' (x_i ) )`

until a sufficiently precise value is reached.

How it works:

The idea is to start with an initial guess which is reasonably close to the true root, then to approximate the function by its tangent line using calculus, and finally to compute the x-intercept of this tangent line by elementary algebra. This x-intercept will typically be a better approximation to the original function's root than the first guess, and the method can be iterated.

More formally, suppose f : (a, b) → ℝ is a differentiable function defined on the interval (a, b) with values in the real numbers ℝ, and we have some current approximation x_i. Then we can derive the formula for a better approximation, x_i + 1 by referring to the diagram on the right. The equation of the tangent line to the curve y = f (x) at x = x_i is

`y=f'(x_i)(x-x_i)+f(x_i)`

where f′ denotes the derivative. The x-intercept of this line (the value of x which makes y = 0) is taken as the next approximation, x_i + 1, to the root, so that the equation of the tangent line is satisfied when (x,y)=(`x_(i+1)`,0)

`0=f'(x_{n})(x_{n+1}-x_{n})+f(x_{n})`

Solving for x_i + 1 gives

`x_{i+1}=x_{i}-{\frac {f(x_{i})}{f'(x_{i})}}`

We start the process with some arbitrary initial value x₀. (The closer to the zero, the better. But, in the absence of any intuition about where the zero might lie, a "guess and check" method might narrow the possibilities to a reasonably small interval by appealing to the intermediate value theorem.) The method will usually converge, provided this initial guess is close enough to the unknown zero, and that f ′(x₀) ≠ 0.

Ice breaking by Air Cushion Vehicle:

Ice sheet can be broken into small slabs with a heavy air cushion vehicle (ACV) traveling on oat ice sheet at the critical speed. Ice-breaking with ACV is a new ice-breaking method in inland river which has important application in prevention of ice jam flood in the Yellow River. Dependent on ice thickness, an ACV can readily break ice in flexure by making a depression and send cracks through the ice. Once a crack is made, the path of least resistance for the ACV's air curtain gets under the ice and reduces the supporting water under the ice. The ice then fails under it’s own weight. At speed an ACV can setup a wave in the ice and send cracks through it. This would act similar to the dynamics of a transport truck on an ice road.

GOVERNING EQUATION:

Where,

p denotes the cushion pressure,

h the thickness of the ice field,

r the size of the air cushion,

σ the tensile strength of the ice ,

β is related to the width of the ice wedge.

PROGRAM:

Programming language used is Python.

# Program to calculate pressure for given thickness of ice using Newton Raphson Method
import math
import numpy as np
import matplotlib.pyplot as plt

# Taking inputs

# Beta is related to the width of the ice wedge
beta = 0.5
# r is the size of the air cushion (feet)
r = 40
# sigma the tensile strength of the ice (pounds per square inch (psi))
sigma = 150
# sigma the tensile strength of the ice (pounds per square foot), 1foot = 12inches, 1 square feet = 12*12 square inches
sigma = 21600
# h the thickness of the ice field
h = [0.6, 1.2, 1.8, 2.4, 3.0, 3.6, 4.2]
# Relaxation factor assumed
alpha = 1.1

# Defining the governing function
def f(p,h,r,beta,sigma):

	term_1 = pow(p, 3)*(1-pow(beta, 2))
	term_2a = 0.4*h*pow(beta, 2)
	term_2b = (sigma*(pow(h,2))/pow(r,2))
	term_3 = pow(sigma,2)*(pow(h,4))/(pow(r,4))
	term_4 = (sigma/3)*(pow(h,2))/(pow(r,2))
	return term_1 + (term_2a - term_2b)*pow(p,2) + term_3*(p/3) - pow(term_4, 3)

# Defining the derivative of the governing function
def fprime(p,h,r,beta,sigma):

	term_5 = 3*pow(p,2)*(1-pow(beta, 2))
	term_6a = 0.4*h*pow(beta, 2)
	term_6b = (sigma*(pow(h,2))/pow(r,2))
	term_7 = pow(sigma,2)*(pow(h,4))/(pow(r,4))
	return term_5 + 2*(term_6a - term_6b)*p + term_7/3


# Null pressure array for storing calculated values after each iterations
pressure = [] 		# pound per square feet
pressure_0 = [] 	# pound per square inches(psi)

# Null array for storing the number of iterations in future
num = []

# Creating a for loop for interation and performing the newtion-raphson technique
for i in h:
	# Assuming initial pressure value
	p_guess = 100
	# Defining Tolerance
	tolerance = 1e-6
	# Loop counter for incrementing the number of iterations
	counter = 1
	while(abs(f(p_guess, i, r, beta, sigma)) > tolerance):
		p_guess = p_guess - alpha*(f(p_guess, i, r, beta, sigma)/fprime(p_guess, i, r, beta, sigma))
		counter = counter+1

	pressure.append(p_guess)
	pressure_0.append(p_guess/144)
	num.append(counter)

print('Thickness of the Ice wedge :',h)
print('Cushion Pressure (pound per square feet) :',pressure)
print('Cushion Pressure (psi) :',pressure_0)
print('Number of iterations :',num)

# Plotting graph
plt.figure(1)
plt.plot(h, pressure, marker = '.',color='blue',markerfacecolor='black',markersize=10)
plt.title('Cushion pressure versus Ice field thickness')
plt.xlabel('Thickness of ice field (feet)')
plt.ylabel('Cushion Pressure (pound per square feet)')
plt.grid()
plt.show()

plt.figure(2)
plt.plot(h, pressure_0, marker = '.',color='blue',markerfacecolor='red',markersize=10)
plt.title('Cushion pressure versus Ice field thickness')
plt.xlabel('Thickness of ice field (feet)')
plt.ylabel('Cushion Pressure (pound per square inches)')
plt.grid()
plt.show()

# Part 2 of the program to calculate optimal relaxation factor for given thickness of ice field
# h_1 the thickness of the ice field
h_1 = 0.6
# Relaxation factor array assumed considering underrelaxing and overrelaxing values
alpha_array = np.linspace(0.1,1.9,50)
# Null array for storing the number of iterations in future
iterations = []

# Creating a for loop for interation and performing the newtion-raphson technique
for j in alpha_array:
	# Assuming initial pressure value
	p_guess = 100
	# Tolerance
	tolerance = 1e-6
	# Loop counter for incrementing the number of iterations
	counter = 1
	while (abs(f(p_guess,h_1,r,beta,sigma)) > tolerance):

		p_guess = p_guess - (j*f(p_guess,h_1,r,beta,sigma))/fprime(p_guess,h_1,r,beta,sigma)

		counter = counter + 1
	iterations.append(counter)
	
print(":"" alpha"" : ""No. of iterations")
for j in range(0,len(alpha_array)):
	print(": ",'%.1f' %alpha_array[j]," : "'%0.1f' %iterations[j])
for k in range(len(iterations)):
	if iterations[k] == min(iterations):
		print('The optimum value of under Relaxation Factor: ', str(alpha_array[k]))

# Plotting graph
plt.figure(3)
plt.plot(alpha_array,iterations,marker = '.',color='black',markerfacecolor='green',markersize=7.5)
plt.xlabel('Relaxation factor (alpha)')
plt.ylabel('No. of iterations')
plt.title('Iterations versus Relaxation factor')
plt.grid()
plt.show()

EXPLANATION OF THE PROGRAM:

1. At first necessary modules such as math, matplotlib, numpy are imported and required inputs are given. For the first part of code an array of and an assumed value of is considered according to the challenge statement.
2. Two different functions named and are defined according to the formula given, and an empty array of number of interations and pressure(in psi and psft) are created.
3. For calculating minimum pressure for different values of ice thickness a for loop is executed taking an initial guessed value of pressure, tolerance, and a loop counter =1 which will increment after every iterations. Then a while loop is initiated inside the same for loop for checking the whether all inputs in the formula if is greater than tolerance provided, a new value of pressure will be calculated according to Newton Raphson method.
4. The calculated values of pressure and number of iterations are stored in their respective arrays using 'append' command and necessary statements are printed in command window.
5. Graphs for pressure versus ice thickness is plotted using calculated data set.
6. For calculating optimal relaxation factor( ) for given thickness of ice same steps as above is repeated for except that now an array of from 0.1 to 1.9 is taken with 50 points in between them and value is taken. Same for loop and while loop is again executed for checking the whether all inputs in the formula if is greater than tolerance provided, a new value of pressure will be calculated according to Newton Raphson method.
7. Number of iterations calculated is stored as an array and a graph is plotted between the same and relaxation factors.

RESULTS:

1. Output window :

1. In the output window we can see the cushion pressure for the respective heights of the ice wedge.

2. When ice thickness was 0.6, pressure corresponding to it was 4.24 pounds per square feet or 0.0294 psi.

2. Variation of cushion pressure with change in thickness of ice field.

3. No of iterations required for different relaxation factors.

It can be concluded that curve rapidly decreases as number of iterations decreases and relaxation factor increases and become constant around `alpha`=1 and again increases with number of iterations and hence optimal relaxation factor lies in the vicinity of 1. By looking at this graph we can conclude that optimal relaxation factor for this problem is 1.01836 .

CONCLUSION:

Newton-Raphson method is simplest method of calculating roots for Non-linear problems. Optimum choice under Relaxation factor helps reducing calculation time by reducing no. of iterations required to calculate the root. Also to restrict minimum value root a tolerance value must be chosen for simplicity in resolving problems.