Breaking ice with air cushion vehicle by Newton-Raphson method using PYTHON

OBJECTIVE:

To determine the minimum cushion pressure needed to break a given thickness of ice
using an air cushion vehicle through Newton-Raphson method using Python.

Inquest:

• Use relaxation factor of 1 and explain the results
• How will you choose the optimum relaxation factor for this case? and what is the optimum relaxation factor for this case? Show with the help of a graph, assume h = 0.6
• Explain your results with respect to the appropriate relaxation factor.

THEORY:

The Newton-Raphson method (also known as Newton's method) is a way to quickly find a good approximation for the root of a real-valued function`f(x)=0`. This method is based on linear approximation where it uses the idea that a continuous and differentiable function can be approximated by a straight line tangent to it.

Here is a picture to demosterate what Newton's method actudally does:

We draw a tangent line to the graph of f(x) at the point `x = x_n`. This line has slope `f'(x_n)`and goes through the point `(x_n, f(x_n)`.Therefore it has the equation`y=f′(x_n)(x−x _n)+f(x _n)`.Now, we find the root of this tangent line by setting `y = 0 and x=x_{n+1}`.for our new approximation. Solving this equation gives us our new approximation, which is

`x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}`.

Both tangents are always intersected at the axis, except when they are in parallel (in such cases, root is not possible).

Equation to be solved:

To determine the minimum pressure Muller (“Ice Breaking with an Air Cushion Vehicle”) in Mathematical Modeling: Classroom Notes in Applied Mathematics, SIAM 1987) derived the equation

Target Equation:

1st order derivative:

Newton-Raphson is performed using:

Here p = cushion pressure

h = thickness of the ice field

r = size of the air cushion

σ = tensile strength of the ice

β = related to the width of the ice wedge.

PROGRAM - 1:

To determine minimum pressure using plot

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

#Defining muller's equation f(x)
def f(p,h):
	term1 = pow(p,3)*((1-pow(beta,2)))
	term2 = ((0.4*h*pow(beta,2) - ((sigma*pow(h,2)/pow(r,2)))*pow(p,2)))
	term3 = ((pow(sigma,2)*pow(h,4)/(3*pow(r,4))))*p
	term4 = (pow((sigma*pow(h,2)/(3*pow(r,2))),3))

	equation = term1 + term2 + term3 - term4
	return equation


def fprime(p,h):
	term1 = 3*pow(p,2)*(1-pow(beta,2))
	term2 = 2*p*((0.4*h*pow(beta,2)) - (sigma*pow(h,2)/(pow(r,2))))
	term3 = ((pow(sigma,2)*pow(h,4)/(3*pow(r,4))))

	equation_prime = term1 + term2 + term3
	return equation_prime

#Input Parameters
beta = 0.5
r = 40
sigma = 150*144
height = [0.6,1.2,1.8,2.4,3.0,3.6,4.2] 

"""
Where beta   = width of ice (feet)
	  r      = size of air cushion (feet)
	  sigma  = tensile strength of ice (psi)
	  height = thickness of ice 
"""

#Performing Newton iteration
p_guess = 40
alpha = 1
tolerance = 1e-4
i = 1

#Appending pressure & iteration
pressure = []
iteration = []

for h in height:
	while(abs(f(p_guess,h))>tolerance):
		p_guess = p_guess - alpha*(f(p_guess,h)/fprime(p_guess,h))
		i = i + 1
	pressure.append(p_guess)

print(pressure)
print(i)

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

#Tabulating
headers = ["Thickness(h)","Pressure(p)"]
table = tuple(zip(height,pressure))
print(tabulate(table,headers,tablefmt = "orgtbl",numalign = "right"))

#Plotting
plt.plot(height,pressure)
plt.title('Pressure vs Thickness')
plt.xlabel('Height(feet)')
plt.ylabel('Pressure(psi)')
plt.show()

Explanation:

• import function is used to import all the required functions such as math, matplotlib.py, numpy & scipy.

Here

math library is used for importing trignometry functions.

matplotlib.py is used for plotting 2d graphs.

numpy is used for assigning same values of same common difference between two values.

tabulate is used for creating tables to represent values.

• Now the muller's equation (given in theory section) is defined as a function where the equation is split into 4 terms and derived as f.
• Muller's equation (given in theory section) is derivated for the 1st order and again created as a function as fprime.
• Input parameters are assigned as β = 0.5, r = 40 feet, σ = 150 pounds per square inch (psi) and h = [0.6, 1.2, 1.8, 2.4, 3.0, 3.6, 4.2] feet.
• We are now assuming p_guess = 40(first guess value of pressure, relaxation factor i.e alpha = 1, tolerance = 1e-4 and i = 1.
• pressure & iteration are created as empty arrays for storing the values.
• Since the value of h is an array, for loop is used to call each value and while loop to determine the greater value comparing with tolerance to do the N-R iteration for the output.
• N-R is performed using the formula stated in the theory section.
• Table is now created to with h in the first column, and p in the second column by using the command tuple(zip(height,pressure)) which helps in arranging the values in desired order.
• Finally Graph has been plotted for Pressure against thickness.

Output:

Plot of Pressure vs Thickness for assuming minimum pressure

Tabulation of Thickness (h) and Pressure (p)

PROGRAM - 2:

To find the min pressure by taking relaxation factor(alpha) as 1 & h =0.6

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

#Defining muller's equation f(x)
def f(p,h):
	term1 = pow(p,3)*((1-pow(beta,2)))
	term2 = ((0.4*h*pow(beta,2) - ((sigma*pow(h,2)/pow(r,2)))*pow(p,2)))
	term3 = ((pow(sigma,2)*pow(h,4)/(3*pow(r,4))))*p
	term4 = (pow((sigma*pow(h,2)/(3*pow(r,2))),3))

	equation = term1 + term2 + term3 - term4
	return equation


def fprime(p,h):
	term1 = 3*pow(p,2)*(1-pow(beta,2))
	term2 = 2*p*((0.4*h*pow(beta,2)) - (sigma*pow(h,2)/(pow(r,2))))
	term3 = ((pow(sigma,2)*pow(h,4)/(3*pow(r,4))))

	equation_prime = term1 + term2 + term3
	return equation_prime

#Input Parameters
beta = 0.5
r = 40
sigma = 150*144
height = [0.6] 

"""
Where beta   = width of ice (feet)
	  r      = size of air cushion (feet)
	  sigma  = tensile strength of ice (psi)
	  height = thickness of ice 
"""

#Performing Newton iteration
p_guess = 10
alpha = np.linspace(0.05,2,50)
tolerance = 1e-4
i = 1
h = 0.6

#Appending pressure & iteration
pressure = []
iteration = []

for a in alpha:
	p_guess = 10
	i = 1
	while(abs(f(p_guess,h))>tolerance):
		p_guess = p_guess - a*(f(p_guess,h)/fprime(p_guess,h))
		i = i + 1
	pressure.append(p_guess)
	iteration.append(i)

x = min(iteration)

for i in range(0,100):
	if i == 0:
		print(alpha[i],x)

#Tabulating
header = ['alpha','iteration']
table = tuple(zip(pressure,iteration))
s = sorted(table)
print(tabulate(table,headers = header, tablefmt = 'pipe'))

#Plotting
plt.plot(alpha,iteration)
plt.title('Optimization of Relaxation factor')
plt.xlabel('Relaxation factor')
plt.ylabel('No. of iteration')
plt.show()

Explanation:

• Importing module, defining functions such as f & fprime is identical as the previous program.
• Input is derived for β = 0.5, r = 40 feet and σ = 150 pounds per square inch (psi) which is same except height which is taken as standard value 0.6 .
• To perform Newton Iteration we assume p_guess = 10, tolerance = 1e-4, i = 1 and use linspace command for alpha for different relaxation factor values.
• pressure & iteration are created as empty arrays for storing the values.
• for loop is now executed for alpha and while for determining greater value for the linearly spaced values for given function.
• x is assigned as optimum relaxation factor.
• It is now ran in a for loop for the length of iterations obtained by the Newton-Raphson method.
• Table is now created using tuple command for alpha & iteration.
• Values of alpha and iteration are plotted in a graph.

Output:

Plot for Optimization of Relaxation factor (alpha)

Tabulation for alpha values for corresponding iterations

RESULTS:

Newton Raphson method is used in Python to find minimum cusion pressure needed to break a given thickness of ice using an air cushion vehicle.