Finding the minimum pressure for breaking ice with air cushion vehicle using the Newton - Raphson method .
AIM:
To find the minimum pressure to break the ice with an air cushion vehicle using Newton - Raphson method in python programming.
THEORY:
Using the newton raphson method, we must determine the minimum pressure that the air cushion vehicle needs in order to break the ice. This has to be shown using the python programming. The equation for finding the minimum pressure need by the air cushion vehicle is as follows:
`p^3(1-beta^2)+(0.4cdotbeta^2-frac{sigmacdoth^2}{r^2})cdotp^2+(frac{sigma^2cdoth^4}{3cdotr^4})cdotp-(frac{sigmacdoth^2}{3cdotr^2})^3=0`_ _ _ (1)
In equation (1) the following parameters represent:
1. p = denotes the cushion pressure on the Hovercraft.
2. h = thickness of the ice field.
3. r = radius of the size of the air cushion.
4. `sigma` = tensile strength of the ice.
5. `beta` = width of the ice wedge.
The respective values assigned to these parameters are
1. `beta` = 0.5.
2. r = 40 feet.
3. `sigma` = 150 (psi)
The equation is defined and we must now apply the newton raphson method to it.
PROBLEM STATEMENT:
1. Use Newton Raphson method to find out the value of pressure for h = 0.6 ft.
2. Find the optimal relaxation factor for this problem using a suitable plot.
3. Tabulate the results for pressure p for h =[0.6,1.2,1.8,2.4,3.0,3.6,4.2] assuming a suitable relaxation factor.
METHODOLOGY:
1. Newton raphson method is an iterative technique. Its objective is to find out the roots.
2. The idea behind using newton raphson method is to linearize the problem and it also gives us them approximated answer in smaller number of iterations.
3. The formula to find out the roots of the given function is :`f(x_n +1) = x_n -frac{f(x_n)}{f'(x_n)}`
4. In order to optimize the solution we need a non- linear solver which is newton raphson method. Its objective is to find the roots for `x` which evaluates the function to zero.
5. To start with it, we need a guess value which is `f(xguess)` which can be any value for `x` which can be the solution.
6. `frac{delf}{delx} = frac{f(xguess)-f(x n e w)}{xguess-xn e w}`
7. `f'(xguess) = frac{f(xguess)-f(x n e w)}{xguess-xn e w}``
8. `xguess - xn e w = frac{f(xguess)}{f'(xguess)}`
9. On rearranging the above equation we get :
`xn ew = xguess -frac{f(xguess)}{f'(xguess)}`
10.Using this formula for finding the minimum pressure to break the ice is sufficient.
11. On using this formula, we might end up doing a lot of iterations to find the exact value at which the function `f(x)` becomes zero. In order to minimize it we are using an under relaxation factor called 'Alpha'. Alpha is going to be less than 1 since is has to be multiplied to `frac{f(xguess)}{f'(xguess)}`, we are using alpha here because it helps us to get a closer value for xnew to minimixe the final answer obtained.
12. In the case of a multidimensional equation like (1), we need to use an under relaxation factor to approximate our final value obtained, which makes us close to the value of x needed to make the function equal to 0.
13. There are chances when Alpha value can be sometimes large and hence it is called over relaxation factor which exceeds the output required.
14. The final formula for newton raphson method is as follows:
`(xn e w) = (xguess) - alphacdotfrac{f(xguess)}{f'(xguess)}`
ALGORITHM:
1. For programming the minimum pressure required by the air cushioned vehicle, we are using python programming. At first we are importing necessary modules for the program. We are using the math module for math calculations and the matplotlib library for plotting our final results obtained.
2. According to the 1st problem statement, we must find the value of pressure at h(Thickness) = 0.6. We are defining the functions as `f(x)` and `fprime(x)`.
#DEFINING FUNCTION DEFINITION
#GIVEN EQUATION
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)
#ON DIFFERENCIATING GIVEN EQUATION
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))In the given function above, we are defining the given equation `p^3(1-beta^2)+(0.4cdotbeta^2-frac{sigmacdoth^2}{r^2})cdotp^2+(frac{sigma^2cdoth^4}{3cdotr^4})cdotp-(frac{sigmacdoth^2}{3cdotr^2})^3=0` as f(p,h) since we need to find the value of pressure for h = 0.6 . We are returning the function to give a the value obtained after compilation. Next we are defining fprime(p,h) which is the differenciation of f(p,h).
3. The value of fprime is obtained by differenciating the function given. The required inputs are defined as follows p_guess = 10,000 (Value of x which is guessed) , alpha = 1 (Underrelaxation factor) , beta(Width of the ice wedge) = 0.5 given , r(Radius of the size of the cushion) = 40 , sigma(Tensile strength of the ice) = 150(psi) , tolerance = 1e-4 , Thickness(Given as an array) h =[0.6,1.2,1.8,2.4,3.0,3.6,4.2].
4. We are using an empty array for pressure = [] so that we can find the value of pressure for each value of thickness provided, it also includes h = 0.6. Now we are using a variable iter to be incremented on each iteration for newton raphson method.
5. Use of newton raphson method takes place over here,
for h in thickness:
while(abs(f(p_guess,h))>tol):
p_guess = p_guess - alpha*(f(p_guess,h)/fprime(p_guess,h))
iter = iter + 1
pressure.append(p_guess)
for i in range(0,7):
print("Thickness: Pressure")
print("Thickness",thickness[i],"pressure",pressure[i])We are using h as a variable in the for loop to run in the thickness array given. A while loop is used in order to find the value of f(p_guess,h) with the value of p_guess = 10,000 in this case and the value of h in taken from the thickness array. We are using the abs function here since there is a chance for the function to generate negative values on different iterations. Hence we are using the absolute function in order to give us only positive values as results. We are setting the tolerance value as 1e-4 since the function must not generate or exceed the value of tolerance. Hence for various values of of p_guess , thickness we are getting solutions. The iter variable is incremented on each loop to show us the number of iterations taking place and the number of times the newton raphson method is used.
6. We are appending the value of p_guess found to the empty pressure array. Then we are using a for loop which runs from (0,7) because the number of values in the thickness array is 7. Hence we are running the for loop and finding the value of pressure for each value of thickness.
7. For plotting the result we used the following code:
print("No of iterations",iter)
plt.plot(thickness,pressure)
plt.xlabel('Thickness')
plt.ylabel('Pressure')
plt.savefig('NEWTON RAPHSON GRAPH')
plt.show()
We are printing the value for the number of iterations. Then we are plotting thickness vs pressure graph using the plt.plot() functions. We have provided xlabel and ylabel as Thicknedd and Pressure . We saved the plot we obtained using the plt.savefig(). To print the result we used plt.show().
8. The final code is as follows:
#FINDING THE MINIMUM PRESSURE ON ICE - HOVERCRAFT USING NEWTON RAPHSON METHOD IN PYTHON PROGRAMMING
#IMPORTING NECESSARY MODULES
import math
import matplotlib.pyplot as plt
#DEFINING FUNCTION DEFINITION
#GIVEN EQUATION
def f(p,h):
return pow(p,3)*(1-pow(beta,2)) + (0.4*h*pow(beta,2) - sigma*pow(h,2)*pow(p,2)/(pow(r,2))) + pow(sigma,2)*pow(h,4)*p/(3*pow(r,4)) - pow((sigma*pow(h,2))/(3*pow(r,2)),3)
#ON DIFFERENCIATING GIVEN EQUATION
def fprime(p,h):
return 3*pow(p,2)*(1-pow(beta,2)) + (0.4*h*pow(beta,2) - sigma*pow(h,2)*2*p/(pow(r,2))) + (pow(sigma,2)*pow(h,4)/(3*pow(r,4)))
#INPUTS
p_guess = 10000 #Value of x which is guessed
alpha = 1 #Under relaxation factor
tol = 1e-3 #Tolerance provided since output must not exceed limit
beta = 0.5 #Width of the ice wedge
r = 40 #Radius of the size of air cushion
sigma = 150*144 #Tensile strength of the ice.
pressure = [] #Empty array for pressure.
thickness = [0.6,1.2,1.8,2.4,3.0,3.6,4.2] #Value of thickness given in question.
iter = 1
for h in thickness:
while(abs(f(p_guess,h))>tol):
p_guess = p_guess - alpha*(f(p_guess,h)/fprime(p_guess,h))
iter = iter + 1
pressure.append(p_guess)
for i in range(0,7):
print("Thickness: Pressure")
print("Thickness",thickness[i],"pressure",pressure[i])
print("No of iterations",iter)
plt.plot(thickness,pressure)
plt.xlabel('Thickness')
plt.ylabel('Pressure')
plt.savefig('NEWTON RAPHSON GRAPH')
plt.show()
On execution we get the following result as
From the result obtained we have found out the pressure = 5.72204589843692e-16 at a thickness of 0.6. Various values for pressure have also been found out.
OUTPUT GRAPH:
From the output obtained we can say at thickness = 1.8 , there is a sudden change in the value of pressure which means that the pressure has increased and the thickness of the ice field has decreased with the increase in pressure . At thickness = 1.8 the pressure was found out to be 5.72204589843692e-16 pascal. It is also found from the calculations that when thickness =1.6 we obtained pressure equal to 5.72204589843692e-16 pascal. This is the method used in newton raphson method. When we obtained the same output after 2 or more iterations. That result is considered to be the approximate answer. From this we can conclude that at h = 0.6 we found pressure = 5.72204589843692e-16 pascal.
PROGRAM TO FIND OPTIMAL RELAXATION FACTOR:
1. Function definition:
#Function definition
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))We have used the same function definition as the previous program. We have returned the function `p^3(1-beta^2)+(0.4cdotbeta^2-frac{sigmacdoth^2}{r^2})cdotp^2+(frac{sigma^2cdoth^4}{3cdotr^4})cdotp-(frac{sigmacdoth^2}{3cdotr^2})^3=0` over here and the function after being differenciated is fprime.
2. Inputs:
beta = 0.5
r = 40
sigma = 150 *144
h = 0.6
tol = 1e-4
these are the given inputs for the following program. Sigma is multiplies by 144 since we are converting the unit per square inch(psi) into per square foot(psf).
3.For defining the under relaxation factor , we are using the numpy module to generate values for alpha in the range of (0.1,1.8,100). We are generating values of alpha from 0.1 to 1.8 and generating values for 100 iterations in newton raphson method.
4. We are using an empty array for iterations and for pressure in this code.
for i in alpha:
counter = 1
p_guess = 10000
while(abs(f(p_guess,h))>tol):
p_guess = p_guess - i*(f(p_guess,h)/fprime(p_guess,h))
counter = counter + 1
pressure.append(p_guess)
iterations.append(counter)
#printing result
print("relaxation factor",':',"minimum pressure",":",'no. of iterations')
for j in range(0,100):
print(alpha[j],':',pressure[j],':',iterations[j])Here we are using a for loop to run in alpha . The value of alpha is taken from the linspace command to generate values for alpha , which is used in newton raphson method. We are then using a while loop to check the value generated by f(p_guess,h) is greater than or lesser than the tolerance value generated. We are also using the abs() function to only return positive values for f(p_guess,h) . Using the newton raphson formula we are generating values for p_guess. The counter variable gets appended on each iteration and it is being returned in the iterations array by appending the value in it. Hence , we are using a for loop to print the values of pressure obtained and no. of iterations generated.
5. For plotting the array we are using the matplotlib library to generate plots for our solutions.
plt.plot(alpha,iterations)
plt.xlabel('Relaxation factor')
plt.ylabel('No of iterations')
plt.savefig('Relaxation factor vs iterations')
plt.show()
6.Final program
#IMPORTING NECESSASRY MODULES
import math
import matplotlib.pyplot as plt
import numpy as np
#Function definition
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))
#INPUTS
beta = 0.5
r = 40
sigma = 150*144 #psf value = psi value x 144
h = 0.6
#Defining alpha
alpha = np.linspace(0.1,1.8,100)
tol = 1e-4
pressure = []
iterations =[]
for i in alpha:
counter = 1
p_guess = 10000
while(abs(f(p_guess,h))>tol):
p_guess = p_guess - i*(f(p_guess,h)/fprime(p_guess,h))
counter = counter + 1
pressure.append(p_guess)
iterations.append(counter)
#printing result
print("relaxation factor",':',"minimum pressure",":",'no. of iterations')
for j in range(0,100):
print(alpha[j],':',pressure[j],':',iterations[j])
plt.plot(alpha,iterations)
plt.xlabel('Relaxation factor')
plt.ylabel('No of iterations')
plt.savefig('Relaxation factor vs iterations')
plt.show() OUTPUT GRAPH:
From the graph obtained we can say that the value of relaxation factor was found to approximately equal to 1.
RESULTS:
From this report we conclude that value of pressure at thickness = 0.6 is 4.370806844 pascal which is sufficient for the air cushion vehicle to break the ice. We also plotted the graph for various values of thickness and pressure. In this report, we found the optimal value of relaxation factor using the newton raphson method.
CONCLUSION:
Newton raphson method is used for finding the roots of a particular function. In this case we are finding the roots equation to acquire the minimum pressure required by the air cushion vehicle to break the ice. This method helps us to approximate our final result by running various values for pressure and thickness . If we get the same value after 2 or more iterations. That result is considered to be the solution. Newton raphson method is a faster approach to finding the required roots to the exact decimal place in smaller number of iterations. Hence, Newton raphson method is an easy way to find the roots of a given equation. In the case of the relaxation factor, the value should be less than one since in that case the there wont be any over relaxation taking place eg. the number of iterations might exceed even after reaching the final answer.