Breaking Ice with Air cushion Vehicle - Find minimum pressure with Newton-Raphson method and finding the optimum relaxation factor for the analysis

#By Abhishek Parmar
import math
import matplotlib.pyplot as plt
import numpy as np 
beta=0.5
r=40
sigma=150*144
#for i in len(h): 
def f_p(p,h):
	t1=((1-pow(beta,2))*pow(p,3))   
	t2=((0.4*h*pow(beta,2))-((sigma*pow(h,2))/(pow(r,2))))*(pow(p,2))
	t3=((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))*p
	t4=(sigma*pow(h,2))/(3*pow(r,2))
	t4=pow(t4,3)
	t5=t1+t2+t3-t4
	return(t5)
def fprime_p(p,h):
	t6=(3*(1-pow(beta,2))*(pow(p,2)))
	t7=(2*(0.4*h*(pow(beta,2)))-((sigma*pow(h,2))/(pow(r,2))))*p
	t8=((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))
	t9=t6+t7+t8
	return(t9)
h=[0.6,1.2,1.8,2.4,3.0,3.6,4.2]
p_guess=100
tol=1e-4
alpha=1
itr=1
p=[]
print(\'pressure\',\'psf |\',\'height\',\'\\n\')
for t in h:
	while(abs(f_p(p_guess,t)) > tol):
		p_guess=p_guess-alpha*(f_p(p_guess,t)/fprime_p(p_guess,t))
		itr= itr+1
	p.append(p_guess)
		
	print(p_guess,\'psf |\',t,\'\\n\')
print(p)
plt.plot(h,p)
plt.xlabel(\'height\')
plt.ylabel(\'pressure\')
plt.title(\'pressure vs height plot\')
plt.xlim([0,4.5])
plt.ylim([0,250])
plt.show()

Seperate code for choosing optimum value of relaxation factor

#to get optimium relaxation factor
import math
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
import numpy as np  
beta=0.5
r=40
sigma=150*144
def f_p(p,h):
	t1=((1-pow(beta,2))*pow(p,3))   
	t2=((0.4*h*pow(beta,2))-((sigma*pow(h,2))/(pow(r,2))))*(pow(p,2))
	t3=((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))*p
	t4=(sigma*pow(h,2))/(3*pow(r,2))
	t4=pow(t4,3)
	t5=t1+t2+t3-t4
	return(t5)
def fprime_p(p,h):
	t6=(3*(1-pow(beta,2))*(pow(p,2)))
	t7=(2*(0.4*h*(pow(beta,2)))-((sigma*pow(h,2))/(pow(r,2))))*p
	t8=((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))
	t9=t6+t7+t8
	return(t9)
h=[0.6,1.2,1.8,2.4,3.0,3.6,4.2]

p_guess=100
tol=1e-4
alpha1=0.15
itr=1
p=[]
i=[]
alpha=[]
while (alpha1<2):
	alpha.append(alpha1)
	while(abs(f_p(p_guess,0.6)) > tol):
		p_guess=p_guess-alpha1*(f_p(p_guess,0.6)/fprime_p(p_guess,0.6))
		itr= itr+1
	p.append(p_guess)
	i.append(itr)
	print(alpha1,\'|\',itr,\'\\n\')
	alpha1=alpha1+0.1
	p_guess=100 #this is again restore to 100 as new p_guess has been obtained
	itr=1 #this is again restore to 1 as it has been updated according to previous loop
	

#print(alpha)
#print(len(i))
ax = plt.axes(projection=\'3d\')
ax.scatter3D(alpha,p,i)
plt.show()
#print(p)
plt.figure(1)
plt.plot(alpha,i)
plt.show()
plt.figure(2)
plt.plot(alpha,p)
plt.show()
# choosing alpha =1.95 in the analysis
p1=[]
for t in h:
	while(abs(f_p(p_guess,t)) > tol):
		p_guess=p_guess-1.95*(f_p(p_guess,t)/fprime_p(p_guess,t))
		itr= itr+1
	p1.append(p_guess)
plt.plot(h,p1)
plt.show()

OR For obtaining the optimum relaxation factor using user defined function(other method)

#creating function to do this
import math
import matplotlib.pyplot as plt
import numpy as np 
from mpl_toolkits import mplot3d
beta=0.5
r=40
sigma=150*144
def f_p(p,h):
	t1=((1-pow(beta,2))*pow(p,3))   
	t2=((0.4*h*pow(beta,2))-((sigma*pow(h,2))/(pow(r,2))))*(pow(p,2))
	t3=((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))*p
	t4=(sigma*pow(h,2))/(3*pow(r,2))
	t4=pow(t4,3)
	t5=t1+t2+t3-t4
	return(t5)
def fprime_p(p,h):
	t6=(3*(1-pow(beta,2))*(pow(p,2)))
	t7=(2*(0.4*h*(pow(beta,2)))-((sigma*pow(h,2))/(pow(r,2))))*p
	t8=((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))
	t9=t6+t7+t8
	return(t9)
h=[0.6,1.2,1.8,2.4,3.0,3.6,4.2]


alpha=np.linspace(0.15,2,40)



def asd(t):
	p_guess=100
	tol=1e-4
	itr=1
	while(abs(f_p(p_guess,0.6)) > tol):
		p_guess=p_guess-t*(f_p(p_guess,0.6)/fprime_p(p_guess,0.6))
		itr=itr+1
	return(p_guess,itr)
p=[]
i=[]	
for t in alpha:
	[a,b]=asd(t)
	p.append(a)
	i.append(b)
#print(p)
#print(i)
plt.figure(1)
plt.plot(alpha,i)
plt.show()
plt.figure(2)
plt.plot(alpha,p)
plt.show()
ax = plt.axes(projection=\'3d\')
ax.scatter3D(alpha,p,i)
plt.show()


	

		

Above image show the variation of No. of iteration,Relaxation factor,Pressure

Above image is for variation for No. of iteration with Relaxation factor

This shows the variation of pressure with relaxation factor

From above plot we see that at relaxation factor of 1.95 the No. of iteration and presure obtained is mininmum so choosing relaxation factor of 1.95 we obtined result as shown