Breaking Ice with Air Cushion Vehicle - Find minimum pressure with Newton Raphson method using Python
Aim:
To write programs in Python for Air Cushion vehicle using Newton Raphson method to:
- Find out the value of pressure for h = 0.6 ft.
- Find the optimal relaxation factor for this problem with the help of a suitable plot.
- 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.
Introduction:
(1) Air Cushion Vehicle:
An Air cushion vehicle or Hovercraft is used to break ice with the help of air pressure. To evaluate this pressure, the following equation is used.
`p^3(1- beta^2) + (0.4hbeta^2 - (sigmah^2)/(r^2))p^2+ (sigma^2h^4)/(3r^4)p - ((sigmah^2)/(3r^2))^3`
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.
The parameters to be used for the programs are given are:
Take β = 0.5, r = 40 feet and σ = 150 pounds per square inch (psi).
(2) Newton Raphson method:
The Newton-Raphson method (N-R) is a way to quickly find a good approximation for the root of a real-valued function f(x)=0. It uses the idea that a continuous and differentiable function can be approximated by a straight line tangent to it.
The N-R technique requires only one initial value `x_0` which we will refer to as the initial guess for the root. The iteration formula for getting a better estimate of the true root is given by:
`x_1 = x_0 - (f(x))/(f^'(x))`
where,
f'(x) denotes the first derivative of f(x) with respect to x.
Thus N-R method finds the tangent to the function f(x) at `x = x_0` and extrapolates it to intersect the x axis to get `x_1`
This point of intersection is taken as the new approximation to the root and the procedure is repeated until convergence is obtained whenever possible. Mathematically, given the value of `x = x_i` at the end of the `i`th iteration, we obtain `x_(i+1)` as,
`x_(i+1) = x_i - (f(x_i))/(f^'(x_i))`
To increase the accuracy of the obtained solution, a relaxation factor (`alpha`) is used:
`x_(i+1) = x_i - alpha((f(x_i))/(f^'(x_i)))`
Python Programs:
(1) To find the value of pressure for h = 0.6ft:
"""
Program to estimate the value of pressure at h = 0.6ft
using Newton Raphson method
By Sree Vignesh
"""
#Import all required modules
import numpy as np
import matplotlib.pyplot as plt
#Function definitions
def f(p,h):
return (pow(p,3)*(1-pow(beta,2)))+(((0.4*h*pow(beta,2))-((sigma*pow(h,2))/pow(r,2))*pow(p,2))+((pow(sigma,2)*pow(h,4))/(3*pow(r,4))*p)-pow((sigma*pow(h,2))/(3*pow(r,2)),3))
def fprime(p,h):
return (3*pow(p,2)*(1-pow(beta,2)))+((0.4*h*pow(beta,2))-((sigma*pow(h,2))/pow(r,2))*(2*p))+((pow(sigma,2)*pow(h,4)))/(3*pow(r,4))
#Relaxation factor
alpha = 1
#Tensile Strength of the ice (in pounds per square feet)
sigma = 150*144
#Size of the air cushion (in feet)
r = 40
#Width of the ice wedge
beta = 0.5
#Thickness of the ice field (in feet)
h = 0.6
#Tolerance
tol = 1e-4
#Iterations
itr = 1
#Guess value for pressure
p_guess = 100
i = 1
while(abs(f(p_guess,h))>tol):
p_guess = p_guess - (alpha*(f(p_guess,h)/fprime(p_guess,h)))
i = i+1
itr = itr + 1
print("The value of pressure for h = 0.6ft is",p_guess)
print("The number of iterations = ",itr)Explanation:
The title of the program is entered as comments and all the required modules are imported at the start.
The function used in the program such as f(p,h) is defined based on the formulae discussed in the Introduction part. fprime(p,h) is the derivative of the function f(p,h).
To find out the value of pressure, relaxation factor alpha is taken as 1 , tolerance tol as 1e-4. The rest of the parameters are taken as stated in the problem statement. The value of sigma is multiplied by 144 to convert it into pounds per square feet.
I have taken the guess value of pressure p_guess as 100 and to make note of the number of iterations, a variable itr is defined with initial value of 1.
Using while loop, the value of pressure is calculated by performing iterations based on the N-R formula and it is stored in p_guess
The value of pressure and the number of iterations are printed out using the print() statement.
Result:
In the output window, the value of pressure is printed out as 4.370779506857015 and has taken 14 iterations.
This value of pressure is for h = 0.6 ft
(2) To find the optimal relaxation factor for the above problem:
"""
Program to find out the optimum relaxation factor
using Newton Raphson method
By Sree Vignesh
"""
#Import all required modules
import numpy as np
import matplotlib.pyplot as plt
#Function definitions
def f(p,h):
return (pow(p,3)*(1-pow(beta,2)))+(((0.4*h*pow(beta,2))-((sigma*pow(h,2))/pow(r,2))*pow(p,2))+((pow(sigma,2)*pow(h,4))/(3*pow(r,4))*p)-pow((sigma*pow(h,2))/(3*pow(r,2)),3))
def fprime(p,h):
return (3*pow(p,2)*(1-pow(beta,2)))+((0.4*h*pow(beta,2))-((sigma*pow(h,2))/pow(r,2))*(2*p))+((pow(sigma,2)*pow(h,4)))/(3*pow(r,4))
#Tensile Strength of the ice (in pounds per square feet)
sigma = 150*144
#Size of the air cushion (in feet)
r = 40
#Width of the ice wedge
beta = 0.5
#Tolerance
tol = 1e-4
#Relaxation factor
alpha_array = np.linspace(0.1,1.95,500)
#Iterations array
iterations = []
for j in alpha_array:
i = 1
#Guess value for pressure
p_guess = 100
while(abs(f(p_guess,0.6))>tol):
p_guess = p_guess - (j*(f(p_guess,0.6)/fprime(p_guess,0.6)))
i = i+1
iterations.append(i)
print("The value of pressure for h = 0.6ft is",p_guess)
for i in range(len(iterations)):
if iterations[i]==min(iterations):
x = alpha_array[i]
y = iterations[i]
print("The optimal relaxation factor is",x)
print("Number of iterations for the optimal relaxation factor is",y)
#Plotting
plt.plot(alpha_array,iterations)
plt.xlabel('Relaxation factor (Alpha)')
plt.ylabel('Iterations')
plt.title('Iterations vs Alpha')
plt.savefig('Iterations vs Alpha.png')
plt.show()Explanation:
The title of the program is entered as comments and all the required modules are imported at the start.
The function used in the program such as f(p,h) is defined based on the formulae discussed in the Introduction part. fprime(p,h) is the derivative of the function f(p,h).
To find out the optimum value of relaxation factor, tolerance tol is taken as 1e-4. The rest of the parameters are taken as stated in the problem statement. The value of sigma is multiplied by 144 to convert it into pounds per square feet.
Here, an array of values is defined for alpha, which is the variable alpha_array . Using np.linspace() a range of values from 0.1 to 1.95 is defined. To take into account the number of iterations, a variable called iterations is defined as an empty array.
For every value of alpha, the pressure is evaluated and hence there is a while loop to evaluate pressure value and it is inside a for loop. Here also, the guess value for pressure p_guess is taken as 100. The iterations are added to the array using append() command.
After this loop, the value of pressure is printed out.
Using another for loop, each value of iterations array is compared and a variable called x is assigned the least value of alpha everytime when the loop is executed. Once the loop ends, the variable x will have the least value of the alpha array, which is the optimal relaxation factor. A variable y is used to hold the value of the least number of iterations.
After this loop, the value of x and y are printed out.
A plot of Iterations vs alpha is plotted with appropriate labels and title.
Result:
From this plot, we can see that, the value of the relaxation factor converges at a value close to 1.25 and further increase in this value, say 1.75 or 1.8 results in Overshoot.
This screenshot from my output window shows the value of pressure for h = 0.6ft and the exact value of optimal relaxation factor for the pressure at h = 0.6ft
On using this value, the number of iterations was found to be 11
(3) To find pressure for different values of h:
"""
To find the values of pressure at different values of h
using Newton Raphson method
By Sree Vignesh
"""
#Import all required modules
import matplotlib.pyplot as plt
#Function definitions
def f(p,h):
return (pow(p,3)*(1-pow(beta,2)))+(((0.4*h*pow(beta,2))-((sigma*pow(h,2))/pow(r,2))*pow(p,2))+((pow(sigma,2)*pow(h,4))/(3*pow(r,4))*p)-pow((sigma*pow(h,2))/(3*pow(r,2)),3))
def fprime(p,h):
return (3*pow(p,2)*(1-pow(beta,2)))+((0.4*h*pow(beta,2))-((sigma*pow(h,2))/pow(r,2))*(2*p))+((pow(sigma,2)*pow(h,4)))/(3*pow(r,4))
#Relaxation factor
alpha = 1.1
#Tensile Strength of the ice
sigma = 150*144
#Size of the air cushion
r = 40
#Width of the ice wedge
beta = 0.5
#Thickness of the ice field
h = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]
#Tolerance
tol = 1e-5
#Iterations
itr = 1
#Guess value for pressure
p_guess = 150
p = []
i = 0
for i in range(len(h)):
H = h[i]
while (abs(f(p_guess,H)) > tol):
p_guess = p_guess - alpha*((f(p_guess,H))/(fprime(p_guess,H)))
itr = itr + 1
i = i + 1
p.append(p_guess)
i = 0
for i in range(len(h)):
if i==0:
print('Height Pressure')
print(h[i]," ",p[i]," ")
#Plotting
plt.plot(h,p)
plt.xlabel('Height')
plt.ylabel('Pressure')
plt.title('Pressure vs Height')
plt.savefig('Pressure vs Height.png')
plt.show()Explanation:
The title of the program is entered as comments and all the required modules are imported at the start.
The function used in the program such as f(p,h) is defined based on the formulae discussed in the Introduction part. fprime(p,h) is the derivative of the function f(p,h).
To find out the value of pressure, relaxation factor alpha is taken as 1.1 , tolerance tol as 1e-5. The rest of the parameters are taken as stated in the problem statement. The value of sigma is multiplied by 144 to convert it into pounds per square feet.
I have taken the guess value of pressure p_guess as 150 and to make note of the number of iterations, a variable itr is defined with initial value of 1. An empty array p is defined to hold the values for pressure at different values of h
For every value of h, the pressure is evaluated and hence there is a while loop to evaluate pressure value and it is inside a for loop. The values of pressure are added to the p array using append() command.
To list down the values, a for loop is used and the values of pressure and height are printed accordingly.
A plot of pressure vs height is plotted with appropriate labels and title.
Result:
This is the plot of Pressure vs Height.
The above is the screenshot of my output window which shows the values of pressure for each respective value of height, identical to a table of values.
Note that, the relaxation factor is assumed a particular value which is not optimal and hence the number of iterations is not taken into account while printing.
Conclusion:
The values of pressure under each given condition were estimated using Newton Raphson method programmed in Python and the results were displayed as stated in the problem statements. A separate program to calculate the optimal relaxation factor was also coded. A detailed explanation was given to the programs and results obtained.