Breaking Ice with Air cushion Vehicle - Finding Minimum Pressure with Newton Raphson Method - Python

OBJECTIVE

1. Write a python code to find out the minimum value of pressure at h=0.6ft by using Newton Raphson method.
2. Write a python code to find the optimal relaxation factor for this problem with the help of a suitable plot.
3. Write a python code to tabulate the results of pressure for h = [0.6,1.2,1.8,2.4,3,3.6,4.2] assuming a suitable relaxation factor.

THEORY

• Newton Raphson Method

The Newton-Raphson method (also known as Newton's method) is a way to quickly find a good approximation for the root of a real-valued function f ( x ) = 0. It uses the idea that a continuous and differentiable function can be approximated by a straight line tangent to it.

Suppose you need to find the root of a continuous, differentiable function f(x), and you know the root you are looking for is near the point `x=x0`​. Then Newton's method tells us that a better approximation for the root is,

`x_1 = x _0 − f ( x_0 )/( f ' ( x_0 ))`

This process may be repeated as many times as necessary to get the desired accuracy. In general, for any x-value x_n​, the next value is given by

`x_ n + 1 = x_ n − f ( x_ n )/( f ' ( x_ n ))`

Here is a picture to demonstrate what Newton's method actually does:

PYTHON CODE

• Python code to find out the minimum value of pressure at h=0.6ft by using Newton Raphson method

import matplotlib.pyplot as plt

"""
Breaking Ice with Air cushion Vehicle - Finding minimum pressure with Newton-Raphson method
"""


# beta - width of the ice wedge
b=0.5

# Size of the air cushion
r=40

# sigma - tensile strength of the ice
s=150*144

# h - thickness of the icefield
h = 0.6



# defining the function
def f(p):
	return (pow(p,3)*(1-pow(b,2)))+(((0.4*h*pow(b,2))-(s*pow(h,2)/pow(r,2)))*pow(p,2))+(pow(s,2)*pow(h,4)*p)/(3*pow(r,4))-(pow((s*pow(h,2))/(3*pow(r,2)),3))

# defining the derivative
def fprime(p):
	return (3*pow(p,2)*(1-pow(b,2)))+(((0.4*h*pow(b,2))-(s*pow(h,2)/pow(r,2)))*2*p)+((pow(s,2)*pow(h,4))/(3*pow(r,4)))

# guess value for pressure
p_guess=10

# relaxation factor
alpha=0.5

# tolerance
tol=1e-4


iter=1

while ( abs(f(p_guess)) > tol ):
	p_guess=p_guess-(alpha*(f(p_guess)/fprime(p_guess)))
	iter=iter+1

# Minimum Pressure
print('Minimum pressure for h = 0.6 ft=',p_guess)

• Python code to find the optimal relaxation factor

import matplotlib.pyplot as plt
import numpy as np

"""
Finding optimum relaxation factor - Breaking Ice with Air cushion Vehicle
"""


# beta - width of the ice wedge
b=0.5

# Size of the air cushion
r=40

# sigma - tensile strength of the ice
s=150*144

# h - thickness of the icefield
h = 0.6



# defining the function
def f(p):
	return (pow(p,3)*(1-pow(b,2)))+(((0.4*h*pow(b,2))-(s*pow(h,2)/pow(r,2)))*pow(p,2))+((pow(s,2)*pow(h,4)*p)/(3*pow(r,4)))-(pow((s*pow(h,2))/(3*pow(r,2)),3))

# defining the derivative
def fprime(p):
	return (3*pow(p,2)*(1-pow(b,2)))+(((0.4*h*pow(b,2))-(s*pow(h,2)/pow(r,2)))*2*p)+((pow(s,2)*pow(h,4))/(3*pow(r,4)))


# relaxation factor
a=np.arange(0.02,2,0.02)

iteration=[]

for alpha in a:

	iter=1

	# tolerance
	tol=1e-4

	# guess value for pressure
	p_guess=10

	while ( abs(f(p_guess)) > tol ):
		p_guess=p_guess-(alpha*(f(p_guess)/fprime(p_guess)))
		iter=iter+1
	iteration.append(iter)

# Optimum relaxation factor
print('Optimum relaxation factor=',a[iteration.index(min(iteration))])

# Plotting
plt.plot(a,iteration)
plt.xlabel('alpha')
plt.ylabel('number of iterations')
plt.title('Number of iterations vs Relaxation factor')
plt.show()

• Python code to tabulate the results of pressure for h = [0.6,1.2,1.8,2.4,3,3.6,4.2]

import matplotlib.pyplot as plt


"""
Plotting Pressure vs Thickness - Breaking Ice with Air cushion Vehicle - Finding minimum pressure with Newton-Raphson method
"""


# beta - width of the ice wedge
b=0.5

# Size of the air cushion
r=40

# sigma - tensile strength of the ice
s=150*144

# h - thickness of the icefield
thickness = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]


# defining the function
def f(p,h):
	return (pow(p,3)*(1-pow(b,2)))+(((0.4*h*pow(b,2))-(s*pow(h,2)/pow(r,2)))*pow(p,2))+(pow(s,2)*pow(h,4)*p)/(3*pow(r,4))-(pow((s*pow(h,2))/(3*pow(r,2)),3))

# defining the derivative
def fprime(p,h):
	return (3*pow(p,2)*(1-pow(b,2)))+(((0.4*h*pow(b,2))-(s*pow(h,2)/pow(r,2)))*2*p)+((pow(s,2)*pow(h,4))/(3*pow(r,4)))

# guess value for pressure
p_guess=10

# optimum relaxation factor
alpha=1.16

# tolerance
tol=1e-4

# creating array for pressure
Pressure=[]

for h in thickness:

	while ( abs(f(p_guess,h)) > tol ):
		p_guess=p_guess-(alpha*(f(p_guess,h)/fprime(p_guess,h)))

	Pressure.append(p_guess)

for i in range(0,7):
	if i==0:
		print("Thickness","|","Pressure")
	print(thickness[i],"      |",Pressure[i])


#plotting

plt.plot(thickness,Pressure)
plt.title('Pressure - Thickness Graph')
plt.xlabel('Thickness')
plt.ylabel('Pressure')
plt.show()

EXPLANATION OF CODES

• import matplotlib.pyplot as plt - It is used to import matplotlib.pyplot, which is a collection of functions that make matplotlib work like MATLAB. Each pyplot function makes some change to a figure: e.g., creates a figure, creates a plotting area in a figure, plots some lines in a plotting area, decorates the plot with labels, etc.
• import numpy as np - It is used to import NumPy, which is a Python library used for working with arrays. It also has functions for working in domain of linear algebra, fourier transform, and matrices.
• def - It is used to define a function.
• abs - The abs() function returns the absolute value of the given number. If the number is a complex number, abs() returns its magnitude.
• while - With the while loop we can execute a set of statements as long as a condition is true.
• while ( abs(f(p_guess)) > tol ): - We use this while loop, instead of for loop to get the required number of iterations till the condition gets false.
• print - The print() function prints the specified message to the screen.
• arange() - The arange() function is used to get evenly spaced values within a given interval. Inside that bracket it will contain 3 values. The first one will be the starting value, the 2nd one will be the ending value and the last one will be the spacing between values(step).
• plt.plot() - The plot() function is used to plot graphs and diagrams.
• plt.xlabel() - The xlabel() function in pyplot module of matplotlib library is used to set the label for the x-axis.
• plt.ylabel() - The ylabel() function in pyplot module of matplotlib library is used to set the label for the y-axis.
• plt.title() - The title() function in pyplot module of matplotlib library is used to set tiltle for the figures.
• plt.show() - It is used to display the figures.
• return - A return statement is used to end the execution of the function call and “returns” the result to the caller.
• min() - The min() function returns the item with the lowest value.
• index() - The index() method returns the position at the first occurrence of the specified value.

RESULT

• Python code to find out the minimum value of pressure at h=0.6ft by using Newton Raphson method

We get the minimum pressure as, 4.2390 Pa

• Python code to find the optimal relaxation factor

Graph - Number of iterations vs Relaxation factor

We get the optimum relaxation factor as, 1.16

• Python code to tabulate the results of pressure for h = [0.6,1.2,1.8,2.4,3,3.6,4.2]

Graph - Pressure vs Thickness

Tabulated Result

COCLUSION

• From the Number of iterations vs Relaxation factor graph we can see that, the iteration getting decreased till it reaches 1.5. After that the number of iterations begin to increase. This happens because of overshooting.
• From the Pressure vs Thickness graph we can see that, the cushion pressure and thickness of ice field are directly proportional. As the thickness increases the cushion pressure also started to increse.