Breaking Ice with Air cushion Vehicle - Find minimum pressure with Newton-Raphson method
AIM: To solve for the minimum value of pressure using the newton raphson method using python.
Theory:
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) = 0. It uses the idea that a continuous and differentiable function can be approximated by a straight line tangent to it.
We draw a tangent line to the graph of f(x) at the point `x=x_n`. This line has slope `f'(x_n)` and goes through the point `(x_n,f(x_n))`. Therefore it has the equation `y=f'(x_n)(x-x_n)+f(x_n)`. Now, we find the root of this tangent line by setting y=0 and `x=x_n+1` for our new approximation. Solving this equation gives us our new approximation, which is `x_(n+1)=x_n-f(x_n)/(f'(x_n))`.
Implementation:
The given function be minimized is as follows which is considered as f(x):
, where
`p` - cushion pressure
`beta` - width of the ice wedge
`h` - thickness of ice field(height)
`sigma` - tensile strength of ice
`r` - size of the air cushion
The following is the code for the minimization of the function by Newton's Raphson Method.
import math
import matplotlib.pyplot as plt
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)*p/(3*pow(r,4)))-pow((sigma*pow(h,2)/(3*pow(r,2))),3)
#differentiate wrt p
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)))
''' performing newton iteration '''
#width of the ice wedge
beta=0.5
#size of the air cushion
r=40
#tensile strength
sigma=150
#height range
h_array=[0.6,1.2,1.8,2.4,3.0,3.6,4.2]
p_guess=0.2
pressure=[]
alpha =1.1
tol = 1e-4
iter=0
for h in h_array:
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)
print(iter)
print('height pressure')
for i in range(0,len(h_array)):
print(h_array[i],' ',pressure[i])
for i in range(0,len(h_array)):
if (h_array[i] == 0.6):
print('at ',h_array[i],'ft of height, pressure =',pressure[i])
plt.plot(h_array,pressure,label=alpha)
plt.xlabel('height')
plt.ylabel('pressure')
plt.legend()
plt.show()Here in the code, initial function and its derivative function is defined. Then the input parameters of r,`beta, sigma` are given. As per the requirement the h values are defined as an array with an initial guess value of pressure as p_guess. A tolerancce value was set so that the evaluation is done only till a lowest particular value after a number of iterations. The relaxation factor (`alpha`) is to be less than 1 so it is called as under relaxed factor, else called as over relaxed factor which is not so good to use as it over shoots the values.
Using a loop the Newton's Raphson method is initiated till the value is less than the tolerance value there by the minimum value for the function is found.
This is carried out for each value of the h array and the respective minimum pressure values found at each iteration are appended to an empty array of pressure in order to plot the h array and pressure array.
To find the optium relaxation factor value the following results from the above code would be helpful.
The graphs are plotted for height v/s pressure, and each curve is for different relaxation factor whose value is mentioned in graph. The following is the table showing the relaxation factor and the respective iterations to find the minimum value of pressure for the given function.
| Relaxation Factor (`alpha`) | Iterations |
| 0.2 | 116 |
| 0.5 | 40 |
| 0.8 | 22 |
| 1 | 13 |
| 1.1 | 11 |
| 1.2 | 16 |
From the above results it is observed that as the relaxation factor increases to 1 the number of iterations are reduced and the minimum value for the function is found immediately but beyond 1, the number of iterations shoots up and the minimum value obtained is comprised along with error. So it suitable for having the factor to be under relaxed i.e. between 0 to 1 with less number of iterations. So we can conclude that the optimum value of the relaxation factor is 1 for the above case.
So for the relaxation factor being 1, the following are the results for the above code and are tabulated within the code and also the minimum pressure for the height of 0.6 feet is about 0.3009 as shown below in the screenshot.
Conclusion:
The minimum value of pressure, at an optimum relaxation factor value of 1, is 0.03009 for the height of 0.6 ft respectively. Here the minimization is conculed to be optimum for the least number of iterations and the relaxation factor to be under relaxation.
References: