Newton Raphson Method - Breaking Ice with Air Cushion Vehicle

Aim :

To find the minimum pressure required to break the ice using Air Cushion Vehicle using Newton-Raphson Method and find optimal relaxation factor using plots.

Objective :

1. Use the newton raphson method to find out the value of pressure for h = 0.6 ft.
2. Find the optimal relaxation factor for this problem with the help of a suitable plot.
3. 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.

Description :

Equation used for determining the pressure known as Muller's Equation is as follows:

where,

p - cushion pressure

h - thickness of the ice field

r - size of the air cushion

σ - tensile strength of the ice

β - related to the width of the ice wedge

The following values for the various parameters are considered,

β = 0.5, r = 40 feet and σ = 150 pounds per square inch (psi)

The above equation can be solved by Newton - Raphson method by initially guessing a value and following the below mentioned equation for every iteration.

`Xng = Xguess -alpha.(F(Xguess))/(Fderivative(Xguess))`

where,

`Xng` = Xnew_guess

`alpha`= relaxation factor

Program Codes :

Code for finding minimum pressure :

"""
Program to find the minimum pressure for various heights
	using Newton - Raphson method
"""

import math
import matplotlib.pyplot as plt
import texttable as tb 

s = 150*144  #conversion from inch to foot
r = 40       
b = 0.5

#function definition
def func(h,p):
	f1 = (pow(p,3)*(1-pow(b,2)))
	f2 = (((0.4*h*pow(b,2))-((s*pow(h,2))/pow(r,2)))*pow(p,2))
	f3 = ((pow(s,2)*pow(h,4)*p)/(3*pow(r,4)))
	f4 = pow(((s*pow(h,2))/(3*pow(r,2))),3)
	final = f1+f2+f3-f4
	return final

#function definition of derivative func
def func1(h,p):
	f1 = (3*pow(p,2)*(1-pow(b,2)))
	f2 = (((0.4*h*pow(b,2))-((s*pow(h,2))/pow(r,2)))*p*2)
	f3 = ((pow(s,2)*pow(h,4))/(3*pow(r,4)))
	final1 = f1+f2+f3
	return final1


height = [0.6,1.2,1.8,2.4,3,3.6,4.2]
pressure=[]
iteration=[]
a=0.5
tol = 0.0001

#iteration process 
for h in height:
	p=15
	it=1
	while(abs(func(h,p))>tol):
		p = p - a*(func(h,p)/func1(h,p))	
		it=it+1
	iteration.append(it)
	pressure.append(p)

#plotting
plt.plot(height,pressure)
plt.xlabel('Thickness of Ice')
plt.ylabel('Pressure')
plt.show()

print("Minimum pressure required for thickness =",height[0],' : ',round(pressure[0],3),'\n')
print('Pressure required for various other thickness are given below :\n')

#tabulation 
tab = tb.Texttable()
x = [[]]
ap=len(pressure)

for i in range(ap):
	x.append([(i+1),height[i],pressure[i],iteration[i]])

tab.add_rows(x)
tab.header(['Sl.No','Thickness', 'Pressure', 'Iterations'])
print(tab.draw())
print('\n\n')

Output for minimum Pressure :

Code for finding Optimal relaxation factor :

"""
Program to find the optimal relaxation factor
	with the help of a plot
"""

import math
import matplotlib.pyplot as plt
import numpy as np 
#import texttable as tb 

h = 0.6
s = 150*144  #conversion from inch to foot
r = 40       
b = 0.5

#function definition
def func(p):
	f1 = (pow(p,3)*(1-pow(b,2)))
	f2 = (((0.4*h*pow(b,2))-((s*pow(h,2))/pow(r,2)))*pow(p,2))
	f3 = ((pow(s,2)*pow(h,4)*p)/(3*pow(r,4)))
	f4 = pow(((s*pow(h,2))/(3*pow(r,2))),3)
	final = f1+f2+f3-f4
	return final

#function definition of derivative func
def func1(p):
	f1 = (3*pow(p,2)*(1-pow(b,2)))
	f2 = (((0.4*h*pow(b,2))-((s*pow(h,2))/pow(r,2)))*p*2)
	f3 = ((pow(s,2)*pow(h,4))/(3*pow(r,4)))
	final1 = f1+f2+f3
	return final1

#creating various no. of values for alpha
alpha = np.linspace(0.1,1.9,100)
iteration=[]
pressure=[]
tol = 0.0001

for a in alpha:
	p=10
	i=1
	while(abs(func(p))>tol):
		p = p - a*((func(p)/func1(p)))
		i = i + 1
	iteration.append(i)
	pressure.append(p)

m = min(iteration)
for i in range(0,len(iteration)):
	if (m == iteration[i]):
		j=i
		break
"""
	the value for alpha in which the no. of 
	iterations is minimum, practically the
	optimal relaxation factor
"""

#plotting
print("Optimal Relaxation Factor : ",round((alpha[j]),2),'\n\n')
plt.plot(alpha[i],m,'*')
plt.legend(['Optimum Relaxation Factor'],loc='best')
plt.plot(alpha,iteration)
plt.xlabel('Relaxation factors')
plt.ylabel('Iterations')
plt.savefig('Figure.png')
plt.show()

Output for optimal relaxation factor :

The optimal relaxation factor is calculated keeping in mind that pressure is approached with less number of iterations for resulting in faster approach of solution and reducing calculation time.