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

Ojective: To find the minimum cushion pressure needed to break a given thickness of ice
using an air cushion vehicle

case1 : Use relaxation factor of 1

case2 : assume thickness of the ice field h, = 0.6

INTRODUCTION:

In numerical analysis, as the Newton–Raphson method, is a root-finding algorithm which produces successively better approximations to the roots (or zeroes) of a real-valued function. The most basic version starts with a single-variable function f defined for a real variable x, the function's derivative f ′, and an initial guess x₀ for a root of f. If the function satisfies sufficient assumptions and the initial guess is close, then

is a better approximation of the root than x₀. Geometrically, (x₁, 0) is the intersection of the x-axis and the tangent of the graph of f at (x₀, f (x₀)): that is, the improved guess is the unique root of the linear approximation at the initial point. The process is repeated as

CASE1:

code

import numpy as np
from scipy.integrate import odeint
import math
import matplotlib.pyplot as plt
from tabulate import tabulate

#defining muuler's equation

def fun(p):
	T1 = (1-pow(B,2))*pow(p,3)
	T2 = ((0.4*h*pow(B,2)) -(t *pow(h,2)/pow(r,2)))*pow(p,2)
	T3 = ((pow(t,2)*pow(h,4))/(3*pow(r,4)))*p
	T4 = t*pow(h,2)/(3*pow(r,2))
	T5 = pow(T4,3)

	x = T1+T2+T3-T5

	return x

def funprime(p):
	T1 = 3*pow(p,2)*(1- pow(B,2))
	T2 = 2*p*((0.4*h*pow(B,2)) -(t *pow(h,2)/pow(r,2)))
	T3 = ((pow(t,2)*pow(h,4))/(3*pow(r,4)))

	x = T1 + T2 + T3

	return x

#inputs

B = 0.5
r = 40
t = 150*144
height = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]
p = 10
tol = 1e-4
i = 1

iterations = []
pressure = []


for h in height :
	while(abs(fun(p))> tol):
		p = p - (fun(p)/funprime(p))
		i = i +1


	iterations.append(i)
	pressure.append(p)

print(pressure)
print(iterations)

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

headers = ["thickness","iterations","pressure"]

table = tuple(zip(height,iterations,pressure))

print(tabulate(table,headers=headers,tablefmt = "github"))

plt.plot(height,pressure)
plt.xlabel("thickness")
plt.ylabel("pressure")
plt.show()

OUTPUT:

EXPLANATION:

As the thickness of the ice feild increases, the minimum pressure acting on the cushion required to break the ice also increases.

CASE2 :

CODE:

import numpy as np
from scipy.integrate import odeint
import math
import matplotlib.pyplot as plt
from tabulate import tabulate


#defining muuler's equation

def fun(p):
	T1 = (1-pow(B,2))*pow(p,3)
	T2 = ((0.4*h*pow(B,2)) -(t *pow(h,2)/pow(r,2)))*pow(p,2)
	T3 = ((pow(t,2)*pow(h,4))/(3*pow(r,4)))*p
	T4 = t*pow(h,2)/(3*pow(r,2))
	T5 = pow(T4,3)

	x = T1+T2+T3-T5

	return x

def funprime(p):
	T1 = 3*pow(p,2)*(1- pow(B,2))
	T2 = 2*p*((0.4*h*pow(B,2)) -(t *pow(h,2)/pow(r,2)))
	T3 = ((pow(t,2)*pow(h,4))/(3*pow(r,4)))

	x = T1 + T2 + T3

	return x

#inputs
B = 0.5
r = 40
t = 150*144
tol = 1e-4
h = 0.6
alpha = np.linspace(0.01,1.99,200)

a =[]
iterations = []


for j in alpha :

	p =10
	i =1

	while(abs(fun(p))> tol):
		p = p - (j*(fun(p)/funprime(p)))
		i = i +1

	iterations.append(i)
	a.append(j)

print(iterations)
print("the minimum iterations is : ", min(iterations))

table = tuple(zip(iterations,a))
s = sorted(table)
print(s[0])

plt.plot(a,iterations,color = "blue",linewidth = 5)
plt.xlabel("alpha")
plt.ylabel("no of iterations")
plt.show()

OUTPUT:

EXPLANATION:

In order to find the minimum pressure to break a ice of thickness( 0.6) at minimum number of iterations,the relaxation factor should be 1.0547