Newton-Raphson Method

In determining the minimum cushion pressure needed to break a given thickness of ice
using an air cushion vehicle, Muller (“Ice Breaking with an Air Cushion Vehicle”) in
Mathematical Modeling: Classroom Notes in Applied Mathematics, SIAM 1987) derived
the 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, and β is related to the width of the ice wedge.
Take β = 0.5, r = 40 feet and σ = 150 pounds per square inch (psi). Construct a table
with h in the first column, and p in the second column for h = 0.6, 1.2, 1.8, 2.4, 3.0, 3.6
and 4.2 feet

#Program to find minimum pressure with Newton-Raphson Mathod
# By Rajeeh Rasheed

import matplotlib.pyplot as plt
import numpy as np 
import math



# Function 

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

#Derivation of Function wrt to Pressure (p)

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


#Performing newton Raphson Method  iteration

# Width of the ice wedge
beta = 0.5   

# Size of the air cushion (in feet)
r = 40  
                      

# Tensile strength of the ice (psi)
sigma = 150    

# Thickness of the ice field
height = [0.6,1.2,1.8,2.4,3.0,3.6,4.2] 

# Cushion Pressure
p_guess = 0.1  

# Relaxation Factor
alpha = 0.8                    

#Tolerance  Error
tol = 1e-4                            

#Pressure
Pressure = []

iter = 1
for h in height :
	while(abs(f(p_guess,h))>tol):
		p_guess = p_guess - ((f(p_guess,h)/fprime(p_guess,h))*alpha)
		iter=iter +1
	Pressure.append(p_guess)
	


for i in range(0,7) :
	if i == 0:
		print ( \": height : Pressure\"  )
		
	print(\":\" , height[i] , \":\" ,  Pressure[i]  )

plt.plot(height,Pressure)
plt.xlabel(\'Height (feet)\')
plt.ylabel(\'Pressure (psi)\')
plt.title(\'Pressure Vs Height -  Newton-Raphson Method Solution\')
plt.xlim([0,5])
plt.ylim([0,0.3])
plt.show()