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

Objective: To solve for the minimum value of pressure using the newton raphson method.

Problem Statement:

For this challenge, a python script is to be written to solve for the minimum value of pressure using the newton raphson method.

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

The result which is to be calculated is,

1. Use the newton raphson method to find out the value of pressure for h = 0.6 ft.

2. Find the optimal relaxation factor for this problem with the help of a suitable plot.

3. Tabulate the results of p for h = [0.6,1.2,1.8,2.4,3,3.6,4.2] assuming a suitable relaxation factor.

Theory:

Newton Raphson Method:

The Newton Raphson method is for solving equations of the form f(x) = 0. We make
an initial guess for the root we are trying to find, and we call this initial guess x0.
The sequence x0, x1, x2, x3, . . . generated in the manner described below should converge to the exact root.
To implement it analytically we need a formula for each approximation in terms of
the previous one, i.e. we need xn+1 in terms of xn.
The equation of the tangent line to the graph y = f(x) at the point (x0, f(x0)) is

The tangent line intersects the x-axis when y = 0 and x = x1, so

Solving this for x1 gives

and, more generally,

Its application to solving equations of the form f(x) = 0, as we now demonstrate, is called the Newton Raphson method. It is guaranteed to converge if the initial guess x0 is close enough, but it is hard to make a clear statement about what we mean by ‘close enough’ because this is highly problem-specific. A sketch of the graph of f(x) can help us decide on an appropriate initial guess x0 for a particular problem.

Procedure:

Case 1: When h = 0.6 ft the value of pressure is mentioned below in the result

import matplotlib.pyplot as plt 
# Using the given equation and defining it as an function
def f(p):
	term1 = pow(p,3)*(1-pow(b,2))
	term2 = ((0.4)*(h)*(pow(b,2))-(a*pow(h,2)/(pow(r,2))))*pow(p,2)
	term3 = ((pow(a,2)*pow(h,4))*p)/(3*pow(r,4))
	term4 = pow(((a*pow(h,2))/(3*pow(r,2))),3)
	all_terms = term1 + term2 + term3 - term4

	return all_terms
# Finding the dervative of the equation and defining in the function
def fprime(p):
	term1 = 3*pow(p,2)*(1-pow(b,2))
	term2 = 2*((0.4)*(h)*(pow(b,2))-(a*pow(h,2)/(pow(r,2))))*p
	term3 = ((pow(a,2)*pow(h,4)))/(3*pow(r,4))
	all_terms = term1 + term2 + term3

	return all_terms

# Inputs

b = 0.5 # It is related to the width of the ice wedge
r = 40	# The size of the air cushion
a = 150 # The tensile strength of the ice
h = 0.6 # Height

# Guess Value
x_guess = 20

# Relaxation Factor
alpha = 0.7

# Pressure
pressure = []

# Tolerance
tol = 1e-4
iter = 1
# Various height values are used in the for loop



while (abs(f(x_guess))>tol):
	x_guess = x_guess - alpha*(f(x_guess)/fprime(x_guess))
	iter = iter+1
pressure.append(x_guess)

print('Height    Pressure')
print(h,"    ",pressure)

# Printing the iterations
print("Total no. of Iterations are ",iter)

In this case, we create the two function which is used in the following steps to emulate the formula as shown in the theory. We pass on the inputs then these inputs are run in the function which is in turn executed in the while loop until the guess value is close to the actual value. A random guess value can be chosen. Then when we get the appropriate value which is stored in pressure variable and displayed at the end along with the no. of iterations.

Case 2: Finding the optimal relaxation factor

import matplotlib.pyplot as plt 
import numpy as np
# Using the given equation and defining it as an function
def f(p):
	term1 = pow(p,3)*(1-pow(b,2))
	term2 = ((0.4)*(h)*(pow(b,2))-(a*pow(h,2)/(pow(r,2))))*pow(p,2)
	term3 = ((pow(a,2)*pow(h,4))*p)/(3*pow(r,4))
	term4 = pow(((a*pow(h,2))/(3*pow(r,2))),3)
	all_terms = term1 + term2 + term3 - term4

	return all_terms
# Finding the dervative of the equation and defining in the function
def fprime(p):
	term1 = 3*pow(p,2)*(1-pow(b,2))
	term2 = 2*((0.4)*(h)*(pow(b,2))-(a*pow(h,2)/(pow(r,2))))*p
	term3 = ((pow(a,2)*pow(h,4)))/(3*pow(r,4))
	all_terms = term1 + term2 + term3

	return all_terms

# Inputs

b = 0.5 # It is related to the width of the ice wedge
r = 40	# The size of the air cushion
a = 150 # The tensile strength of the ice
h = 0.6 # Height

# Relaxation Factor
alpha = []
ctr =[]
addn = 0.03

# Tolerance
tol = 1e-4
iter = 1
# Various height values are used in the for loop
min_iter = 12000
min_alpha  = 2.5

while(addn<2):
	alpha.append(addn)
	iter = 1
	# Guess Value
	x_guess = 20
	while (abs(f(x_guess))>tol):
		x_guess = x_guess - addn*(f(x_guess)/fprime(x_guess))
		iter = iter+1
	ctr.append(iter)
	if (min_iter>iter):
		min_iter = iter
		min_alpha = addn
	addn = addn + 0.03


# Printing the iterations
print("Total no. of Iterations are ",min_iter)
# Plotting the data
print(" Relaxation Factor: ",min_alpha)
plt.plot(alpha,ctr)
plt.title("Optimal Relaxation Factor")
plt.xlabel('Relaxation Factor')
plt.ylabel('No. of Iterations')
plt.show()

In this case, we create the two functions similar to the above case which is used in the following steps to emulate the formula as shown in the theory. We pass on the inputs then these inputs are run in the function which is in turn executed in the while loop until the guess value is close to the actual value. Here, we use multiple relaxation factor values out of which the appropriate value is chosen.A random guess value can be chosen. Then when the if condition satisfies we get the relaxation factor and the no. of iterations which are plotted in the end. The plots are displayed below.

Case 3: Tabulate the values of height and pressure by assuming relaxation factor

import matplotlib.pyplot as plt 
# Using the given equation and defining it as an function
def f(p,h):
	term1 = pow(p,3)*(1-pow(b,2))
	term2 = ((0.4)*(h)*(pow(b,2))-(a*pow(h,2)/(pow(r,2))))*pow(p,2)
	term3 = ((pow(a,2)*pow(h,4))*p)/(3*pow(r,4))
	term4 = pow(((a*pow(h,2))/(3*pow(r,2))),3)
	all_terms = term1 + term2 + term3 - term4

	return all_terms
# Finding the dervative of the equation and defining in the function
def fprime(p,h):
	term1 = 3*pow(p,2)*(1-pow(b,2))
	term2 = 2*((0.4)*(h)*(pow(b,2))-(a*pow(h,2)/(pow(r,2))))*p
	term3 = ((pow(a,2)*pow(h,4)))/(3*pow(r,4))
	all_terms = term1 + term2 + term3

	return all_terms

# Inputs

b = 0.5 # It is related to the width of the ice wedge
r = 40	# The size of the air cushion
a = 150 # The tensile strength of the ice
heights = [0.6,1.2,1.8,2.4,3.0,3.6,4.2] # Height

# Guess Value
x_guess = 20

# Relaxation Factor
alpha = 0.8

# Pressure
pressure = []

# Tolerance
tol = 1e-4
iter = 1
# Various height values are used in the for loop


for height in heights:
	while (abs(f(x_guess,height))>tol):
		x_guess = x_guess - alpha*(f(x_guess,height)/fprime(x_guess,height))
		iter = iter+1
	pressure.append(x_guess)

j = 0
print('Height    Pressure')
for h in heights:
	print(h,"    ",pressure[j])
	j+=1

# Printing the iterations
print("Total no. of Iterations are ",iter)

# Plot Results
plt.plot(heights,pressure)
plt.title('Pressure vs Height')
plt.xlabel('Height')
plt.ylabel('Pressure')
plt.show()	

In this case, we create the two functions similar to the above case which is used in the following steps to emulate the formula as shown in the theory. We pass on the inputs then these inputs are run in the function which is in turn executed in the while loop until the guess value is close to the actual value. A random guess value can be chosen. We use a for loop here to obtain each value of height because there are multiple heights in this case. Then when we get the appropriate value which is stored in the pressure variable and displayed at the end along with the no. of iterations and height. The plots are displayed below.

Result:

The result for case 1.

The plot and result for case 2.

The plot and result for case 3.

Conclusion:

We can conclude from our observation that as the depth or height keeps increasing the amount of pressure required to break the ice also increases.

References:

http://personal.maths.surrey.ac.uk/st/S.Gourley/NewtonRaphson.pdf