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

Equation for Ice Breaking with an Air Cushion Vehicle:- 

`p^3*(1-beta^3)+(0.4*h*beta^2-(sigma*h^2)/r^2)*p^2+(sigma^2*h^4)/(3*r^4)*p-((sigma*h^2)/(3*r^2))^3=0`

p = pressure

h = thickness of the ice

`sigma= " tensile strength of the ice"`

r = size of the air cushion

`beta="Function of width of the ice wedge"`

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

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

# Equation

def eqn(beta, sigma, r, h, p):
	e = np.arange(0.0,2.0)
	term1 = 1-np.power(beta,2)
	term2 = 0.4*h*np.power(beta,2) - sigma*np.power(h/r,2)
	term3 = np.power(sigma,2)*np.power(h/r,4)/3
	term4 = np.power(sigma*np.power(h/r,2)/3,3)
	e[0] = np.power(p,3)*term1+np.power(p,2)*term2+p*term3-term4 # Equation
	e[1] = 3*np.power(p,2)*term1+2*p*term2+term3                 # Derivative of the Equation
	return e


# Newton Raphson	

def new_rap(ini_Val, iter, TolVal, alpha, beta, sigma, r, h):
	p = ini_Val
	for i in np.arange(1,iter+1):
		e = eqn(beta, sigma, r, h, p)
		if (abs(e[0])>TolVal): p = p - alpha*(e[0]/e[1])
		else: break	
	return p,e,i


# main	
ini_Val = 23040.0      # Initial estimate of pressure (psf)
iter = 100.0        # Maximum Number of iterations	
TolVal = 1e-5      # Allowed deviation of funtion value from 0
alpha = 1         # Relaxation
beta = 0.5
sigma = 150.0*144.0       # tensile strength of the ice (psf-pounds per square foot)
r = 40.0            # size of the air cushion (feet)
h = np.array([0.6, 1.2, 1.8, 2.4, 3.0, 3.6, 4.2]) # thickness of the ice field (feet)	
p = np.arange(0.0,len(h))	

print('\t\tTable')
print('\n\t h\t|\tp(psf)\t | Value of Eqn | Number of iterations')
for i in np.arange(0,len(h)):
	p[i],e,ni = new_rap(ini_Val, iter, TolVal, alpha, beta, sigma, r, h[i]) 
	# print('\nNumber of iterations = {}, Equation value = {}'.format(ni,e[0]))	
	print('\n\t{}\t|\t{:.6f} | {:.6f}\t\t\t | {}'.format(h[i],p[i],e[0],ni))


# Plot
pressure = np.arange(-200.0,200.1,0.1)
e = np.zeros((len(h), len(pressure), 2))

for j in np.arange(0,len(h)):
	for i in np.arange(0,len(pressure)):
		x = eqn(beta, sigma, r, h[j], pressure[i])
		e[j,i,0] = x[0]
		e[j,i,1] = x[1] 

plt.plot(h,p)
plt.xlabel('h')
plt.ylabel('p')
plt.show()

plt.figure()
for i in np.arange(0,len(h)):
	plt.plot(pressure,e[i,:,0])
	plt.plot([-100,100],[p[i],p[i]])
	plt.title('h = {} , p = {}'.format(h[i],p[i]))
	plt.xlabel('Pressure')
	plt.ylabel('Value of Equation')
	plt.show()

Output:-

From the table, the minimum value of pressure occurs for value of h = 0.6.

The code will give 7 such plots for 7 values of h.