Breaking Ice with Air cushion Vehicle - Find minimum pressure with Newton-Raphson method Python
#Plotting h vs pressure table for alpha = 1
import math
import matplotlib.pyplot as plt
def f(p,h,b,r,rho):
return pow(p,3)*(1-pow(b,2)) + pow(p,2)*(0.4*h*pow(b,2)-rho*pow(h,2)/pow(r,2)) + p*pow(rho,2)*pow(h,4)/3/pow(r,4) - pow(rho*pow(h,2)/3/pow(r,2),3)
def f_prime(p,h,b,r,rho):
return 3*pow(p,2)*(1-pow(b,2))+2*p*(0.4*h*pow(b,2)-rho*pow(h,2)/pow(r,2))+pow(rho,2)*pow(h,4)/3/pow(r,4)
#Given values
b = 0.5
r = 40 #feet
r = r*12 #feet to inch
rho = 150 #pounds per square inch
h_array = [0.6,1.2,1.8,2.4,3.0,3.6,4.2] #feet
tol = 1e-4 #tolerance
alpha = 1 #relaxation factor
pressure = []
print(\"h(feet) minimum pressure (psi)\")
#Finding the minimum pressure
for j in h_array:
iteration = 0
x_guess = 10
j = j*12 #feet to inch
while(abs(f(x_guess,j,b,r,rho))>tol):
x_guess = x_guess - alpha*f(x_guess,j,b,r,rho)/f_prime(x_guess,j,b,r,rho)
iteration = iteration + 1
print(\'%.1f %.20f\' %(j/12,x_guess))
pressure.append(x_guess)
plt.plot(h_array,pressure,marker=\"o\")
plt.title(\"Minimum pressure required, relaxation factor = 1\")
plt.xlabel(\"Thickness of ice (feet)\")
plt.ylabel(\"Pressure (psi)\")
plt.show()#Finding optimum relaxation factor for h = 0.6
import math
import matplotlib.pyplot as plt
def f(p,h,b,r,rho):
return pow(p,3)*(1-pow(b,2)) + pow(p,2)*(0.4*h*pow(b,2)-rho*pow(h,2)/pow(r,2)) + p*pow(rho,2)*pow(h,4)/3/pow(r,4) - pow(rho*pow(h,2)/3/pow(r,2),3)
def f_prime(p,h,b,r,rho):
return 3*pow(p,2)*(1-pow(b,2))+2*p*(0.4*h*pow(b,2)-rho*pow(h,2)/pow(r,2))+pow(rho,2)*pow(h,4)/3/pow(r,4)
#Given values
b = 0.5
r = 40 #feet
r = r*12 #feet to inch
rho = 150 #pounds per square inch
h = 0.6 #feet
tol = 1e-6
iteration_alpha = []
pressure = []
alpha_array = []
print(\"h x_ans f(x_ans) iteration alpha\")
h = h*12 #feet to inch
alpha = 0.1 #initial alpha
while(alpha<4):
iteration = 0
x_guess = 100
#Finding the minimum pressure
while(abs(f(x_guess,h,b,r,rho))>tol):
x_guess = x_guess - alpha*f(x_guess,h,b,r,rho)/f_prime(x_guess,h,b,r,rho)
iteration = iteration + 1
alpha_array.append(alpha)
iteration_alpha.append(iteration)
alpha = alpha + 0.1
print(h/12,\" \",x_guess,\" \",f(x_guess,h,b,r,rho),\" \",iteration,\" \",alpha-0.1)
plt.plot(alpha_array,iteration_alpha,marker=\"o\",markersize=3)
plt.plot(alpha_array[25],iteration_alpha[25],marker=\"o\",markersize=7,markeredgewidth=0.5,markeredgecolor=\"g\",markerfacecolor=\"None\")
plt.annotate(\'2.6\',xy=(alpha_array[25],iteration_alpha[25]+10))
plt.plot([0,4],[iteration_alpha[25],iteration_alpha[25]],\'--\',linewidth=0.5)
plt.title(\'Tolerance = 1e-6\')
plt.xlabel(\'Relaxation factor\')
plt.ylabel(\'Number of iterations\')
plt.show()Part 1 of the challenge: Plotting minimum pressure required for various thickness of ice for given relaxation factor (alpha = 1)
Define the function given as a definition (‘f’). Define the partial derivative of the function with respect to pressure (‘f_prime’)
Set the initial values of constants as given.
Define an array of various ice thickness value (‘h_array’)
Define the tolerance and the alpha (‘relaxation factor’)
For every thickness value in h_array, find the minimum pressure.
The value of the function is zero at minimum pressure. However, achieving zero is not feasible mathematically. We find the solution as close to zero as possible. We thus give a tolerance value.
Thus, while the absolute value of the function is greater than the tolerance, we find new solution for every iteration:
`x_\"new\" =x_\"old\"-α f(x_\"old\" )/(f^\' (x_\"old\") )`
Thus, find the minimum pressure for each thickness value and print a table as well as a figure.
Thus, we can see that the minimum pressure required to break the ice increases with the increase in thickness of the ice.
Part 2 of the challenge: Find optimum relaxation factor for ice thickness = 0.6 feet.
Optimum relaxation factor can be decided as the one which can find the solution within minimum number of iterations and the one which gives more accuracy in the answer.
Since, we are giving the tolerance value we need, the accuracy of the answer can’t be considered as a significant factor.
Thus, we find the number of iterations required for each factor.
The following graph was plotted by recording the number of iterations required.
It can be observed that the number of iterations is high for 0.1 and it keeps decreasing till factor = 2.6 and again starts increasing.
Testing the relaxation factor beyond 4 is impractical.
The result has been found consistent even with changing tolerance.
Results:
.................................................................................................................................
For tolerance = 1e-3, Minimum iterations found at ,`alpha =2.6`
.................................................................................................................................
.................................................................................................................................
For tolerance = 1e-4, Minimum iterations found at ,`alpha =2.6`
.................................................................................................................................
.................................................................................................................................
For tolerance = 1e-5, Minimum iterations found at ,`alpha =2.6`
.................................................................................................................................
.................................................................................................................................
For tolerance = 1e-6, Minimum iterations found at ,`alpha =2.4` which is very close to 2.6 and the number of iterations difference is also close.
15 iterations for ,`alpha =2.4`
5 iterations for,`alpha =2.6`
.................................................................................................................................
.................................................................................................................................
Thus, the optimum relaxation factor is 2.6