Breaking ice with Air Cushion Vehicle, Finding Min pressure using Newton-Raphson Method

Introduction:

Named after the Issac Newton and Joseph Rapson . This method is one of the best method to solve the complex numerical equation, this method is based on the linear approximation. This method has the fast rate of convergence than other methods .

This method uses Tangent line of curve to locate the approximate range of roots of the equation.

Both tangents are always intersect at the axis except when are parallel ,in that case there is no root possible. The intersection points of tangents are converge to the true root

A closer approximation to the root is given by:

The above formula is known as the Newton-Raphson formula for the newton iteration.

`x_(n+1)=x_n-(f(x_n )/(f^' (x_n ) ))`

n determining the minimum cushion pressure needed to break a given thickness of ice
using an air cushion vehicle, Muller (“Ice Breaking with an Air Cushion Vehicle”) in
Mathematical Modeling: Classroom Notes in Applied Mathematics, SIAM 1987) derived
the equation

`p^3 (1-ß^2 )+(0.4hß^2-((σh^2)/r^2 )) p^2+(σ^2 h^4)/(3r^4 ) p-((σh^2)/(3r^2 ))^3=0`

We are going to perform Newton-Raphson method for the above function using python programming.

The programme for the above function with the following inputs:

where p denotes the cushion pressure, h the thickness of the ice field, R the size of the air
cushion, σ the tensile strength of the ice, and β is related to the width of the ice wedge.
Take β = 0.5, r = 40 feet and σ = 150 pounds per square inch (psi). Construct a table
with h in the first column, and p in the second column for h = 0.6, 1.2, 1.8, 2.4, 3.0, 3.6
and 4.2 feet.

And the relaxation factor of 1 is used .

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



#defining f(x)
def f(p):
	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))

	y = term1+term2+term3-term4
	return y

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

#co efficients are
r = 40 							#sze of the air cushion in feet
sigma = 150*144					#tensile strength of ice in psi
beta = 0.5 						#width of the ice wedge in feet
height = [0.6,1.2,1.8,2.4,3.0,3.6,4.2] #thickness of the ice field

# performing neeton iteration
p_guess = 80
alpha = 1
tol = 1e-4
iter = 1
P_array = []


for h in height:
	while (abs(f(p_guess))>tol):
		p_guess = p_guess - alpha*f(p_guess)/fprime(p_guess)
		iter = iter+1

	P_array.append(p_guess)


print(P_array)
print(iter)

#tabulating the pressure and thickness values

headers = ["Thickness (h)","Pressure (p)"]
tabel = tuple(zip(height,P_array))
print(tabulate(tabel,headers,tablefmt = "orgtbl",numalign="right"))

#plotting
plt.plot(height,P_array)
plt.xlabel('Thickness (in ft)')
plt.ylabel('Pressure (in psf)')
plt.show()

Output of the above program:

The plot of graph of the pressure vs thickness of the ice is :

Explanation of the code:

Importing matplotlib.pyploy to plot the graph ,

Importing tabulate to table the values of the pressure output and height values.

The above function is defined in the function as f(x)

And again an another function is defined as the derivative the function f(x) as fprime(x)

The co efficients that are in the equation are defined with the given values β, r, σ, h,

For loop is used to perform the newton iteration for all the values of the thickness of the ice. i.e. h=[0.6,1.2,1.8,2.4,3.0,3.6,4.2]

The values are tabulated using the tabulate command

First table values are defined by the table variable which are the values of pressure and thickness values.

Header command for the labels of the respective column

Tuple and zip command is used in the table variable to get the tabular values in the correct order , tuple aarranges the values of an array , as the tabulate command takes the values of rows only, to zip those values accordingly the zip command is used

Then later the tabulated values are plotted.

The programme for finding the optimum relaxation factor value is as below:

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



#defining f(x)
def f(p):
	term1 =  pow(p,3)*float((1-pow(beta,2)))
	term2 = ((0.4*h*pow(beta,2) - float((sigma*pow(h,2)/pow(r,2)))*pow(p,2)))
	term3 = float((pow(sigma,2)*pow(h,4)/(3*pow(r,4))))*p
	term4 = float(pow((sigma*pow(h,2)/(3*pow(r,2))),3))

	y = term1+term2+term3-term4
	return y

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

#co efficients are
r = 40 							#sze of the air cushion in feet
sigma = 150*144					#tensile strength of ice in psi
beta = 0.5 						#width of the ice wedge in feet
h = 0.6 #thickness of the ice field

# performing newton iteration
p_guess = 80
alpha = np.linspace(0.05,2,50)
tol = 1e-4
iter = 1
iterations = []


for i in alpha:
	p_guess = 80
	iter = 1
	while (abs(f(p_guess))>tol):
		p_guess = p_guess - i*f(p_guess)/fprime(p_guess)
		iter = iter+1

	iterations.append(iter)

x = min(iterations)

for i in range(len(iterations)):
	if iterations[i] == x:
		print(alpha[i],x)

print(p_guess)
print(iterations)
print(alpha)
print(x)

plt.plot(alpha,iterations)
plt.xlabel('relaxation factor')
plt.ylabel('no of iterations')
plt.show()

The output foer the programme obtained as follows:

The graph for the relaxation factor vs iterations is :

The explanation for the above code is:

The only difference in this programme is that we are trying find the optimum relaxation factor in this.

So we have created an empty array of iterations to fill with values obtained after the for loop is executed,

The for loop for the linearly spaced values of alpha array is executed for the newton-raphson method for the given function.

Here X is the optimum relaxation factor value

That relaxation factor is ran in a for loop for the length of the iterations optained by the newton-rapson method.