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

Introduction:

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) = 0f(x)=0. It uses the idea that a continuous and differentiable function can be approximated by a straight line tangent to it.

Suppose you need to find the root of a continuous, differentiable function f(x)f(x), and you know the root you are looking for is near the point x = x_0x=x0​. Then Newton's method tells us that a better approximation for the root is

`x 1 ​ =x 0 ​ − f (x 0 ​) /(f'(x 0 ​ ) ​ .`

This process may be repeated as many times as necessary to get the desired accuracy. In general, for any x-value `x_n`​, the next value is given by

xn+1​=xn​−f′(xn​)/f(xn​)​

The term "near" is used loosely because it does not need a precise definition in this context. However, x0​ should be closer to the root you need than to any other root (if the function has multiple roots).

Formula:

In 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

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

Program:

To find the cushion temperature for every thickness of ice field using relaxation factor 1

import tabulate as tbt
#p is cushion pressure
#h is thickness of ice field
def f(p,h):#defining the governing equation
	term1=pow(p,3)*(1-pow(b,2))
	term2 = ((0.4*h*b**2)-(s*h**2/r**2))*p**2
	term3 = s**2*pow(h,4)*p/(3*pow(r,4))
	term4 = pow((s*h**2/(3*r**2)),3)
	return term1+term2+term3-term4


def fprime(p,h):#differentiated equation with respect to pressure
	term1 = 3*p**2*(1-b**2)
	term2 = 2*p*((0.4*h*b**2)-(s*h**2/r**2))
	term3 =  s**2*pow(h,4)/(3*pow(r,4))
	return term1+term2+term3


h = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]
b =0.5
r =40
s = 150#sigma
x_guess =10#initial guess
tol =1e-5#tolerance of error
alpha=1#relaxation factor
headers = ['thickness','cushion pressure']#headers for table
itera =[]#counting iteration while solving
P=[]#used to contain p values
for a in h:
	i=1
	while(abs(f(x_guess,a))>tol):
	     x_guess = x_guess - (alpha*f(x_guess,a)/fprime(x_guess,a))
	     i=i+1
	     itera.append(i)
	     P.append(x_guess)
#to tabulate thickness and cushion pressure
data =zip(h,P)
print(tbt.tabulate(data,headers,tablefmt="grid"))

Output:

Explanation:

• The equation is splitted into 4 terms for simplicity and function f(p,h) is defined
• The differentiated function f'(p,h) is defined
• A tolerance for the calculated pressure value is given
• Headers for the table is given
• a table is formed for p & h values using tabulate module
• The table contains the pressure values for each thickness of ice field
• A relaxation factor of 1 is used here
• It is inferred that least pressure is required to break the ice of least thickness and vice versa

To find the optimum(minimum) relaxation factor for h=0.6

import tabulate as tbt
import matplotlib.pyplot as plt
import numpy as np
#p is cushion pressure
#h is thickness of ice field
def f(p,h):#defining the governing equation
	term1=pow(p,3)*(1-pow(b,2))
	term2 = ((0.4*h*b**2)-(s*h**2/r**2))*p**2
	term3 = s**2*pow(h,4)*p/(3*pow(r,4))
	term4 = pow((s*h**2/(3*r**2)),3)
	return term1+term2+term3-term4


def fprime(p,h):#differentiated equation with respect to pressure
	term1 = 3*p**2*(1-b**2)
	term2 = 2*p*((0.4*h*b**2)-(s*h**2/r**2))
	term3 =  s**2*pow(h,4)/(3*pow(r,4))
	return term1+term2+term3


h = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]
b =0.5
r =40
s = 150#sigma
x_guess =10#initial guess
tol =1e-5#tolerance of error
alpha=1#relaxation factor
headers = ['thickness','cushion pressure']
itera =[]#iteration for solving
P=[]
for a in h:
	i=1
	while(abs(f(x_guess,a))>tol):
		x_guess = x_guess - (alpha*f(x_guess,a)/fprime(x_guess,a))
		i=i+1
	itera.append(i)
	P.append(x_guess)
#to tabulate thickness and cushion pressure
data =zip(h,P)
print(tbt.tabulate(data,headers,tablefmt="grid"))
#to find the optimum relaxation factor using h=0.6
Alpha =np.linspace(0.1,5,50)
J=[]
CT =[]
L=[]
for j in Alpha:
	x_p=10
	ct=1
	while(abs(f(x_p,0.6))>tol):
		x_p = x_p - (j*f(x_p,0.6)/fprime(x_p,0.6))
		ct=ct+1
	J.append(j)
	CT.append(ct)
	L.append(ct)
W=zip(L,Alpha)
w=sorted(W)#finding the minimum no.of iteration
mydata1 =tuple(zip(J,CT))
headers1 =["relaxation factor","iteration"]
print(tbt.tabulate(mydata1,headers=headers1,tablefmt="grid"))
print(w[0])
plt.figure()
plt.title('cushion pressure vs thickness')
plt.xlabel('Thickness')
plt.ylabel('Pressure')
plt.plot(h,P)
plt.show()
plt.title('minimum relaxation factor is 3 having 2 iterations ')
plt.xlabel('Relaxation factor')
plt.ylabel('Iteration')
plt.plot(J,CT)
plt.show()

Output:

Explanation:

The same procedure as the first program is followed except h=0.6, taken

In order to find the optimum relaxation factor, an array is created of many values of factors betweem 0.1 and 5 in steps of 0.1

using each factor, the iterations are run to obtain different p values

The no.of iterations for each factor is found and the factor that is having least iterations is the optimum value. In this way optimum relaxatio factor is chosen.

Here Optimum relaxation factor is 3 showing 2 factors for h=0.6

The graph is plotted between no.of iterations and relaxation factor.

When using h=0.6 and changing alpha to 3, the simulation ran in 2 iterations as expected. The code is

import tabulate as tbt
#p is cushion pressure
#h is thickness of ice field
def f(p,h):#defining the governing equation
	term1=pow(p,3)*(1-pow(b,2))
	term2 = ((0.4*h*b**2)-(s*h**2/r**2))*p**2
	term3 = s**2*pow(h,4)*p/(3*pow(r,4))
	term4 = pow((s*h**2/(3*r**2)),3)
	return term1+term2+term3-term4


def fprime(p,h):#differentiated equation with respect to pressure
	term1 = 3*p**2*(1-b**2)
	term2 = 2*p*((0.4*h*b**2)-(s*h**2/r**2))
	term3 =  s**2*pow(h,4)/(3*pow(r,4))
	return term1+term2+term3


#h = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]
h=0.6
b =0.5
r =40
s = 150#sigma
x_guess =10#initial guess
tol =1e-5#tolerance of error
alpha=3.004#relaxation factor
headers = ['thickness','cushion pressure']#headers for table
itera =[]#counting iteration while solving
P=[]#used to contain p values
i=1
while(abs(f(x_guess,0.6))>tol):
	     x_guess = x_guess - (alpha*f(x_guess,0.6)/fprime(x_guess,0.6))
	     i=i+1
	     itera.append(i)
	     P.append(x_guess)
itera=str(itera)
print('no. of iteration is '+ itera)

Output:

no. of iteration is [2].