Week 4 -Breaking Ice with Air cushion Vehicle - Find minimum pressure with Newton-Raphson method

OBJECTIVE :

• Write a python code to find out the minimum value of pressure at h=0.6ft by using the Newton Raphson method.
• Write a python code to find the optimal relaxation factor for this problem with the help of a suitable plot.
• 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 :

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.

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

`x_("n+1")=x_("n")-alpha*(f(x_"n")/(f'(x_"n")))`

where,

`x_"n+1"` is the next estimate or new guess value

`x_"n"` is the current estimate or old guess value or starting value

`alpha` is the relaxation factor

`f(x_"n")` is the value of the function using the current estimate, and

`f'(x_"n")` is the value of the derivative of the given function using the current estimate

CODE 1:

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

import matplotlib.pyplot as plt

beta = 0.5         # width of the ice wedge (-)
r = 40		       # size of the air cushion (feet)
sigma = 150*144    # tensile strength of the ice (psf)
h = 0.6            # thickness of the icefield (feet)

def f(p,h):
	f = (1-beta**2)*p**3 + (0.4*h*beta**2 - sigma*h**2/r**2)*p**2 + (sigma**2*h**4)*p/(3*r**4) - ((sigma*h**2)/(3*r**2))**3
	return f

def fprime(p,h):
	fprime = (3*p**2*(1-beta**2) + 2*p*(0.4*h*beta**2 - sigma*h**2/r**2) + sigma**2*h**4/(3*r**4))
	return fprime

p_guess=10
alpha=1
tol=1e-4
iter=1
iterations=[]
P=[]

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

iterations.append(iter)
P.append(p_guess)

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

OUTPUT 1:

CODE 2:

#Plotting Iterations to converge vs Relaxation factor - Finding  the optimum relaxation factor

import matplotlib.pyplot as plt
import numpy as np

beta = 0.5       # width of the ice wedge (-)
r = 40		     # size of the air cushion (feet)
sigma = 150*144  # tensile strength of the ice (psf)
h = 0.6          # thickness of the icefield (feet)

def f(p,h):
	f = (1-beta**2)*p**3 + (0.4*h*beta**2 - sigma*h**2/r**2)*p**2 + (sigma**2*h**4)*p/(3*r**4) - ((sigma*h**2)/(3*r**2))**3
	return f

def fprime(p,h):
	fprime = (3*p**2*(1-beta**2) + 2*p*(0.4*h*beta**2 - sigma*h**2/r**2) + sigma**2*h**4/(3*r**4))
	return fprime


#alpha_range=np.linspace(0.02,2,50)  
alpha_range = np.arange(0.02,2,0.02)     



iterations=[]
P=[]

for alpha in alpha_range:
	iter=1
	p_guess=10
	tol=1e-4

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

	iterations.append(iter)
	#P.append(p_guess)

#determination of the optimum relaxation factor
for m in range(len(iterations)):

	if iterations[m]== min(iterations):
		alpha_opt=str(alpha_range[m])
		print("The optimum relaxtion factor is:",alpha_opt)
		

# Optimum relaxation factor
#print('Optimum relaxation factor=',alpha_range[iterations.index(min(iterations))])

plt.figure(1)
plt.plot(alpha_range,iterations)
#plt.plot(alpha_opt,iter_min,marker='o',color='blue',label='point corresponding to minimum number of iterations')
#plt.legend()
plt.xlabel(' Relaxation factor (alpha)')
plt.ylabel('Iterations to converge')
plt.title('Iterations vs Relaxation factor')
plt.grid ('on')
plt.savefig('Optimum alpha')
plt.show()

OUTPUT 2:

CODE 3:

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

import matplotlib.pyplot as plt

beta = 0.5                               # width of the ice wedge (-)
r = 40		                             # size of the air cushion (feet)
sigma = 150*144                          # tensile strength of the ice (psf)
h_range = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]  # thickness of the icefield (feet)

def f(p,h):
	f = (1-beta**2)*p**3 + (0.4*h*beta**2 - sigma*h**2/r**2)*p**2 + (sigma**2*h**4)*p/(3*r**4) - ((sigma*h**2)/(3*r**2))**3
	return f

def fprime(p,h):
	fprime = (3*p**2*(1-beta**2) + 2*p*(0.4*h*beta**2 - sigma*h**2/r**2) + sigma**2*h**4/(3*r**4))
	return fprime

p_guess=10
alpha=1.16   #optimum relaxation factor 
tol=1e-4
iter=1
iterations=[]
P=[]

for h in h_range:

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

	iterations.append(iter)
	P.append(p_guess)

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

for i in range(0,7):
	if i==0:
		print("Thickness of the ice field ","|","Minimum Pressure")
		print("--------------------------------------------------")

	print(h_range[i],"                        |",P[i])

plt.figure(1)
plt.plot(h_range,P,marker='o',markerfacecolor='red')
plt.xlabel('Thickness of ice (in ft)')
plt.ylabel(' Minimum cushion pressure  ' r'($ frac{pounds}{sqft} $)')
plt.title('Cushion pressure vs Ice thickness')
plt.grid ('on')
plt.savefig('Pressure vs Thickness')
plt.show()

OUTPUT 3: