Aim: Breaking ice with air cushion vehicle and find minimum pressure with netwon raphson method.

Introduction: the netwon raphson method simply called as newtons method, is fast way to find the good approximate for root of a real valued function.

in the following program the ice as to break by cushion vehicle, the program written to find the minimum presuure to break ice.

where,

p= cushion pressure

h= thickness of ice field

R=size of air cushion

program to create a table for different height and pressure

# finding mim pressure at difernt heights
import matplotlib.pyplot as plt


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))
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)))

beta=0.5
r=40
sigma=150
height=[0.6,1.2,1.8,2.4,3.0,3.6,4.2]
p_guess=0.1
alpha=1
tol= 1e-4
pressure=[]
iter=1


for h in height:
	while(abs(f(p_guess,h))>tol):
		p_guess=p_guess-alpha*((f(p_guess,h)/fprime(p_guess,h)))
		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\')
plt.ylabel(\'pressure\')
plt.show()

Explaination

1] function command is used for equation.

2] the value of h is varied with increment of 0.6 for 7 values.

3] for each value of h pressure is found as showm in below table afoe relaxiation factor alpha=1

4] the graph is obtained for pressure vs height.

program to minimum relaxation factor

# to find optimum relaxation factor for h=0.6
import math
import matplotlib.pyplot as plt
def f(p,h,beta,r,sigma):
	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))
def fprime(p,h,beta,r,sigma):
	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)))

beta=0.5
r=40
sigma=150
h=0.6
tol=1e-5
iter_alpha=[]
pressure=[]
alpha_array=[]
print(\"h x_ans f(x_ans) iteration alpha\")
alpha=0.1 #intial 
while(alpha<4):
	iter=0
	p_guess=100
	while(abs(f(p_guess,h,beta,r,sigma))>tol):
		p_guess=p_guess-alpha*f(p_guess,h,beta,r,sigma)/fprime(p_guess,h,beta,r,sigma)
		iter=iter+1
		alpha_array.append(alpha)
		iter_alpha.append(iter)
		alpha=alpha+0.1
		print(h/6,\"\",p_guess,\"\",f(p_guess,h,beta,r,sigma),\"\",iter,\"\",alpha-0.1)
plt.plot(alpha_array,iter_alpha,marker=\"o\",markersize=3)
plt.plot(alpha_array[25],iter_alpha[25],marker=\'o\',markersize=7,markeredgewidth=0.5,markeredgecolor=\"g\",markerfacecolor=\"none\")
plt.annotate(\'2.6\',xy=(alpha_array[25],iter_alpha[25]+10))
plt.plot([0,4],[iter_alpha[25],iter_alpha[25]],\"--\",linewidth=0.5)
plt.xlabel(\'relaxation factor\')
plt.ylabel(\'number of iter\')
plt.show()

expalination

1] optimum realxation factor can be decided as the one will find mimimum number of iterations

2] if number of iterations is hig for 0.1 and it keeps reducing till factor 2.6 and again increasing

Conclusion:

program wriiten and found solution for given condition.